The TELEMAC-MASCARET system  trunk
bedload_solvs_vf_gaia.f
Go to the documentation of this file.
1 ! ********************************
2  SUBROUTINE bedload_solvs_vf_gaia
3 ! ********************************
4 !
5  &(mesh,qsx,qsy,limtec,unsv2d,ebor,breach,nseg,nptfr,npoin,
6  & kent,kdir,kddl,dt,flux,flbcla,
7  & liqbor,qbor,nubo,vnoin,evcl_mb,ratio_sand,xmvs)
8 !
9 !***********************************************************************
10 ! GAIA
11 !***********************************************************************
12 !
14 ! with the finite volume method.
15 !
16 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
40 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
41 !
42  USE interface_gaia, ex_bedload_solvs_vf => bedload_solvs_vf_gaia
43  USE bief
46  IMPLICIT NONE
47 !
48 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
49 !
50  TYPE(bief_mesh), INTENT(INOUT) :: MESH
51  TYPE(bief_obj), INTENT(IN) :: QSX,QSY,LIMTEC,UNSV2D,EBOR
52  TYPE(bief_obj), INTENT(IN) :: BREACH,LIQBOR,QBOR
53  INTEGER, INTENT(IN) :: NSEG,NPTFR,NPOIN,KENT,KDIR,KDDL
54  DOUBLE PRECISION, INTENT(IN) :: DT
55  TYPE(bief_obj), INTENT(INOUT) :: FLBCLA,FLUX
56  TYPE(bief_obj), INTENT(INOUT) :: EVCL_MB
57  INTEGER, INTENT(IN) :: NUBO(2,nseg)
58  DOUBLE PRECISION, INTENT(IN) :: VNOIN(3,nseg)
59  DOUBLE PRECISION, INTENT(IN) :: RATIO_SAND(npoin),XMVS
60 !
61 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
62 !
63  INTEGER :: ISEGIN,K,N,IEL,IEL1,IEL2,IELEM,NELEM,ERR
64  DOUBLE PRECISION :: QSMOY1,QSMOY2,QSP,VNOIN1,VNOIN2,RNORM,XN,YN
65  DOUBLE PRECISION :: EVCL_MDIR,PROD_SCAL
66  LOGICAL, ALLOCATABLE :: YESNO(:)
67 !
68  nelem = mesh%NELEM
69  ALLOCATE(yesno(nseg),stat=err)
70 ! INITIALIZATION OF YESNO
71  DO k=1,nseg
72  yesno(k)=.false.
73  ENDDO
74 !
75 ! INITIALISES THE DIVERGENCE
76 !
77  CALL os('X=0 ', x=flux)
78 !
79 ! DETERMINES THE OUTGOING FLUX FOR EACH CELL
80 ! CONVENTION: POSITIVE FLUX IS OUTGOING FLUX
81 !
82  DO ielem=1,nelem
83  DO k=1,3
84  isegin = mesh%ELTSEG%I(ielem+(k-1)*nelem)
85  IF(.NOT.yesno(isegin)) THEN
86  iel1 = nubo(1,isegin)
87  iel2 = nubo(2,isegin)
88 !
89 ! 1 - SEGMENT LENGTH (RNORM)
90 !
91  vnoin1 = vnoin(1,isegin)
92  vnoin2 = vnoin(2,isegin)
93  rnorm = vnoin(3,isegin)
94  prod_scal= (mesh%X%R(iel2)-mesh%X%R(iel1))*vnoin1+
95  & (mesh%Y%R(iel2)-mesh%Y%R(iel1))*vnoin2
96  IF(prod_scal.LT.0.d0)THEN
97  iel1 = nubo(2,isegin)
98  iel2 = nubo(1,isegin)
99  ENDIF
100 !
101 ! 2 - QS FOR THE SEGMENT, BROKEN UP ACCORDING TO X AND Y
102 ! CENTRED FLUX: ONLY USED TO FIND THE SOLID TRANSPORT DIRECTION
103 !
104  qsmoy1 = 0.5d0*(qsx%R(iel1) + qsx%R(iel2))
105  qsmoy2 = 0.5d0*(qsy%R(iel1) + qsy%R(iel2))
106 !
107 ! 3 - PROJECTS QS FOR THE SEGMENT ONTO THE SEGMENT NORMAL
108 !
109  qsp = vnoin1*qsmoy1 + vnoin2*qsmoy2
110 !
111 ! 4 - UPWIND SCHEME OR CENTRED SCHEME
112 !
113  IF(qsp.GT.0.d0) THEN
114 ! NODE WITH A "PROBLEM": UPWIND
115  IF(breach%I(iel1).EQ.1) THEN
116  qsmoy1 = qsx%R(iel1)
117  qsmoy2 = qsy%R(iel1)
118  ELSE
119 ! UPWIND OR CENTRED ACCORDING TO DVF
120  qsmoy1 = dvf*qsx%R(iel1) + (1.d0-dvf)*qsx%R(iel2)
121  qsmoy2 = dvf*qsy%R(iel1) + (1.d0-dvf)*qsy%R(iel2)
122  ENDIF
123  ELSE !QSP.LE.0.D0
124 ! NODES WITH A "PROBLEM": UPWIND
125  IF(breach%I(iel2).EQ.1) THEN
126  qsmoy1 = qsx%R(iel2)
127  qsmoy2 = qsy%R(iel2)
128  ELSE
129 ! UPWIND OR CENTRED ACCORDING TO DVF
130  qsmoy1 = (1.d0-dvf)*qsx%R(iel1) + dvf*qsx%R(iel2)
131  qsmoy2 = (1.d0-dvf)*qsy%R(iel1) + dvf*qsy%R(iel2)
132  ENDIF
133  ENDIF
134  qsp = vnoin1*qsmoy1 + vnoin2*qsmoy2
135 !
136 ! 5 - INTEGRATES BY THE SEGMENT LENGTH:
137 ! ONLY IF NODE IS WET
138 !
139  IF(hn%R(iel1).GE.hmin_bedload.AND.
140  & hn%R(iel2).GE.hmin_bedload) THEN
141  IF(mesh%IFABOR%I(ielem+(k-1)*nelem).NE.-2) THEN
142  flux%R(iel1) = flux%R(iel1) + rnorm*qsp
143  flux%R(iel2) = flux%R(iel2) - rnorm*qsp
144  ELSE
145 ! IF IT IS AN INTERFACE EDGE IN PARALLEL, ONLY HALF
146  flux%R(iel1) = flux%R(iel1) + 0.5d0*rnorm*qsp
147  flux%R(iel2) = flux%R(iel2) - 0.5d0*rnorm*qsp
148  ENDIF
149  ENDIF
150  yesno(isegin)=.true.
151  ENDIF
152  ENDDO
153  ENDDO
154 !
155 ! BOUNDARIES
156 !
157  IF(ncsize.LE.1) THEN
158  DO k = 1 , nptfr
159  iel = mesh%NBOR%I(k)
160 ! VECTOR NORMAL TO A BOUNDARY NODE
161 ! VERSION WHICH IS NOT NORMED
162  xn = mesh%XNEBOR%R(k+nptfr)
163  yn = mesh%YNEBOR%R(k+nptfr)
164 !
165 ! ADDING BOUNDARY FLUX AT OPEN BOUNDARIES
166 ! QBOR HAS PRIORITY HERE, SO IN CASE OF LIQBOR=KENT
167 ! LIMTEC=KDIR WILL NOT BE CONSIDERED, SEE BEDLOAD_SOLVS_FE_GAIA
168 ! FOR MORE EXPLANATIONS
169 !
170 ! QBOR AND FLUX MUST HAVE THE SAME DIMENSION !
171 !
172  IF(liqbor%I(k).EQ.kent) THEN
173  flbcla%R(k) = -qbor%R(k)
174  flux%R(iel) = flux%R(iel) + flbcla%R(k)
175  ELSEIF(limtec%I(k).EQ.kddl.OR.limtec%I(k).EQ.kdir) THEN
176 ! IF KDIR WILL BE UPDATED LATER
177  flbcla%R(k)= qsx%R(iel)*xn + qsy%R(iel)*yn
178 ! ADDS THE CONTRIBUTION OF THE FLUX ON THE BOUNDARY SEGMENT
179  flux%R(iel) = flux%R(iel) + flbcla%R(k)
180  ELSE
181 ! NO SEDIMENT FLUX ACCROSS SOLID BOUNDARIES
182  flbcla%R(k)=0.d0
183  ENDIF
184  ENDDO
185  ELSE
186 ! COMPUTED ONLY ONCE WHEN THE LIQUID BOUNDARY NODE BELONGS TO
187 ! 2 SUBDDOMAINS
188  DO k = 1 , nptfr
189  iel = mesh%NBOR%I(k)
190 ! VECTOR NORMAL TO A BOUNDARY NODE
191 ! VERSION WHICH IS NOT NORMED
192  xn = mesh%XNEBOR%R(k+nptfr)
193  yn = mesh%YNEBOR%R(k+nptfr)
194 !
195 ! ADDING BOUNDARY FLUX AT OPEN BOUNDARIES
196 ! QBOR HAS PRIORITY HERE, SO IN CASE OF LIQBOR=KENT
197 ! LIMTEC=KDIR WILL NOT BE CONSIDERED, SEE BEDLOAD_SOLVS_FE_GAIA
198 ! FOR MORE EXPLANATIONS
199 !
200 ! QBOR AND FLUX MUST HAVE THE SAME DIMENSION !
201 !
202  IF(liqbor%I(k).EQ.kent) THEN
203  flbcla%R(k) = -qbor%R(k)
204  flux%R(iel) = flux%R(iel) + flbcla%R(k)*mesh%IFAC%I(iel)
205  ELSEIF(limtec%I(k).EQ.kddl.OR.limtec%I(k).EQ.kdir) THEN
206 ! IF KDIR WILL BE UPDATED LATER
207  flbcla%R(k)= qsx%R(iel)*xn + qsy%R(iel)*yn
208 ! ADDS THE CONTRIBUTION OF THE FLUX ON THE BOUNDARY SEGMENT
209  flux%R(iel) = flux%R(iel) + flbcla%R(k)*mesh%IFAC%I(iel)
210  ELSE
211 ! NO SEDIMENT FLUX ACCROSS SOLID BOUNDARIES
212  flbcla%R(k)=0.d0
213  ENDIF
214  ENDDO
215  ENDIF
216 !
217 ! ASSEMBLING IN PARALLEL
218 !
219  IF(ncsize.GT.1) CALL parcom(flux, 2, mesh)
220 !
221 ! SOLVING, NEGATIVE SIGN BECAUSE OUTGOING FLUX IS POSITIVE
222 ! NOTE JMH: FLUX MUST BE HERE DIV(QS)
223 !
224  CALL os('X=CYZ ', x=evcl_mb, y=flux, z=unsv2d, c=-dt)
225 !
226 ! ebor*xmvs replaced by a unique value given by the user? Not done
227 !
228  DO k=1,nptfr
229 ! PRIORITY OF LIQBOR, SEE ABOVE
230  IF(limtec%I(k).EQ.kdir.AND.liqbor%I(k).NE.kent) THEN
231  n = mesh%NBOR%I(k)
232  evcl_mdir = ratio_sand(n)*ebor%R(k)*xmvs
233 ! CORRECTION OF BOUNDARY FLUX TO GET EVCL_MDIR
234  flbcla%R(k)=flbcla%R(k)-(evcl_mdir-evcl_mb%R(n))/
235  & (dt*unsv2d%R(n))
236 ! EVCL_MDIR FINALLY IMPOSED
237  evcl_mb%R(n) = evcl_mdir
238  ENDIF
239 !
240  ENDDO
241 !
242  DEALLOCATE(yesno)
243 !
244 !======================================================================!
245 !======================================================================!
246 !
247  RETURN
248  END
subroutine breach
Definition: breach.f:4
double precision dvf
Upwinding for Exner FV.
type(bief_obj), target hn
Water depth.
subroutine bedload_solvs_vf_gaia(MESH, QSX, QSY, LIMTEC, UNSV2D, EBOR, BREACH, NSEG, NPTFR, NPOIN, KENT, KDIR, KDDL, DT, FLUX, FLBCLA, LIQBOR, QBOR, NUBO, VNOIN, EVCL_MB, RATIO_SAND, XMVS)
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
subroutine parcom(X, ICOM, MESH)
Definition: parcom.f:7
double precision hmin_bedload
Minimum depth for bedload.
Definition: bief.f:3