The TELEMAC-MASCARET system  trunk
positive_depths_eria.f
Go to the documentation of this file.
1 ! *******************************
2  SUBROUTINE positive_depths_eria
3 ! *******************************
4 !
5  &(t1,t2,t3,t4,h,hn,mesh,flodel,compute_flodel,flbor,dt,
6  & unsv2d,npoin,gloseg1,gloseg2,nbor,nptfr,
7  & smh,yasmh,pluie,rain,optsou,limpro,hbor,kdir,info,
8  & flopoint,namecode,nitmax,makeflulimebe,flulimebe)
9 !
10 !***********************************************************************
11 ! BIEF V7P3
12 !***********************************************************************
13 !
14 !brief Suppresses negative depths by a limitation of fluxes, with the
15 !+ ERIA technique.
16 !
17 !history J-M HERVOUET (EDF LAB, LNHE)
18 !+ 19/11/2016
19 !+ V7P3
20 !+ First version, taken out of previous positive_depths
21 !
22 !history J-M HERVOUET (EDF LAB, LNHE)
23 !+ 23/11/2016
24 !+ V7P3
25 !+ Adding an option (hardcoded so far) for sharing volumes of points
26 !+ between elements (see variable OPTPRE, OPTPRE=1 old strategy,
27 !+ OPTPRE=2, new and simpler strategy). Tests do not show much
28 !+ difference, kept here to be tested in other cases.
29 !
30 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
31 !| COMPUTE_FLODEL |-->| IF YES, COMPUTE FLODEL HERE
32 !| MAKEFLULIMEBE |-->| OPTIONAL, IF YES DOES ARRAY FLULIMEBE
33 !| DT |-->| TIME STEP
34 !| FLBOR |<->| BOUNDARY FLUXES
35 !| FLODEL |<->| FLUXES GIVEN BY SEGMENT
36 !| | | MAY BE COMPUTED HERE (SEE COMPUTE-FLODEL)
37 !| | | OR SIMPLY GIVEN. AT THE EXIT, THE REAL FLUX
38 !| | | TRANSMITTED IS GIVEN BACK.
39 !| FLOPOINT |-->| FLUXES GIVEN BY POINTS (ELEMENT BY ELEMENT)
40 !| FLULIMEBE |<--| AS FLULIM BUT GIVEN PER ELEMENT
41 !| | | HENCE FLULIMEBE(NELMAX,3)
42 !| GLOSEG1 |-->| FIRST POINT OF SEGMENTS
43 !| GLOSEG2 |-->| SECOND POINT OF SEGMENTS
44 !| H |<->| NEW DEPTH
45 !| HBOR |-->| PRESCRIBED DEPTHS AT BOUNDARIES
46 !| HN |-->| OLD DEPTH
47 !| INFO |-->| IF YES, PRINTING INFORMATION ON LISTING
48 !| KDIR |-->| CONVENTION FOR DIRICHLET BOUNDARY CONDITION
49 !| LIMPRO |-->| TYPE OF BOUNDARY CONDITIONS
50 !| | | IF EQUAL TO KDIR: PRESCRIBED DEPTH.
51 !| MESH |<->| MESH STRUCTURE
52 !| NAMECODE |-->| NAME OF CALLING CODE (SISYPHE, TELEMEC2D, ETC.)
53 !| NBOR |-->| GLOBAL NUMBERS OF BOUNDARY POINTS
54 !| NITMAX |-->| MAXIMUM NUMBER OF ITERATIONS
55 !| NPOIN |-->| NUMBER OF POINTS IN THE MESH
56 !| NPTFR |-->| NUMBER OF BOUNDARY POINTS
57 !| OPTSOU |-->| OPTION FOR SOURCES 1: NORMAL 2: DIRAC
58 !| PLUIE |-->| RAIN IN A BIEF_OBJ, IN M/S.
59 !| RAIN |-->| IF YES, THERE IS RAIN OR EVAPORATION
60 !| SMH |-->| SOURCE TERMS
61 !| T1 |-->| WORK ARRAY
62 !| T2 |-->| WORK ARRAY
63 !| T3 |-->| WORK ARRAY
64 !| T4 |-->| WORK ARRAY
65 !| UNSV2D |-->| INVERSE OF INTEGRAL OF BASIS FUNCTIONS
66 !| YASMH |-->| IF(YES) SMH MUST BE TAKEN INTO ACCOUNT
67 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
68 !
69  USE bief, ex_positive_depths_eria => positive_depths_eria
70  USE declarations_telemac, ONLY : deja_pdept_eria,indic_pdept_eria
71  USE interface_parallel
72 !
73  USE declarations_special
74  IMPLICIT NONE
75 !
76 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
77 !
78  INTEGER, INTENT(IN) :: NPOIN,NPTFR,OPTSOU,KDIR
79  INTEGER, INTENT(IN) :: NITMAX
80  INTEGER, INTENT(IN) :: GLOSEG1(*),GLOSEG2(*)
81  INTEGER, INTENT(IN) :: NBOR(nptfr)
82  INTEGER, INTENT(IN) :: LIMPRO(nptfr)
83  DOUBLE PRECISION, INTENT(IN) :: DT,HBOR(nptfr)
84  TYPE(bief_mesh),INTENT(INOUT) :: MESH
85  DOUBLE PRECISION, INTENT(INOUT) :: FLOPOINT(mesh%nelmax,3)
86  TYPE(bief_obj), INTENT(INOUT) :: T1,T2,T3,T4,FLODEL,H,FLBOR
87  TYPE(bief_obj), INTENT(INOUT) :: PLUIE
88  TYPE(bief_obj), INTENT(IN) :: UNSV2D,HN,SMH
89  LOGICAL, INTENT(IN) :: YASMH,INFO,RAIN
90  LOGICAL, INTENT(IN) :: COMPUTE_FLODEL
91  CHARACTER(LEN=24) :: NAMECODE
92  LOGICAL, INTENT(IN) :: MAKEFLULIMEBE
93  DOUBLE PRECISION, INTENT(INOUT) :: FLULIMEBE(mesh%nelmax,3)
94 !
95 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
96 !
97  INTEGER I,I1,I2,I3,IPTFR,REMAIN,NEWREMAIN,IR,NITER
98  INTEGER NELEM,NELMAX,NSEG,IELEM,OPTPRE
99  DOUBLE PRECISION C,CPREV,CINIT,VOL1,VOL2,VOL3,HFL1,F1,F2,F3
100  DOUBLE PRECISION DTLIM1,DTLIM2,DTLIM3,FP1,FP2,FP3,DT1,DT2,DT3
101  DOUBLE PRECISION SURDT,A1,A2,A3
102  DOUBLE PRECISION, PARAMETER :: TIERS=1.d0/3.d0
103 !
104  DOUBLE PRECISION, PARAMETER :: EPS_FLUX = 1.d-15
105  LOGICAL, PARAMETER :: TESTING = .false.
106 !
107 !-----------------------------------------------------------------------
108 !
109 ! VARIANTS
110 !
111 ! OPTPRE MUST BE EQUAL TO ITS VALUE IN CVTRVF_ERIA !!!!!!!!!!
112  optpre=1
113 !
114 !-----------------------------------------------------------------------
115 !
116 ! INDIC_PDEPT_ERIA WILL BE A LIST OF SEGMENTS WITH NON ZERO FLUXES
117 ! HSEG IS THE DEPTH SHARED BETWEEN SEGMENTS
118 !
119 !-----------------------------------------------------------------------
120 !
121  surdt=1.d0/dt
122  nelem=mesh%NELEM
123  nelmax=mesh%NELMAX
124  nseg=mesh%NSEG
125 !
126  IF(.NOT.deja_pdept_eria) THEN
127  ALLOCATE(indic_pdept_eria(nelem))
128  deja_pdept_eria=.true.
129  ENDIF
130 !
131 !-----------------------------------------------------------------------
132 !
133  IF(testing) THEN
134  c=1.d99
135  cinit=1.d99
136  DO i=1,npoin
137  c =min(c ,h%R(i))
138  cinit=min(cinit,hn%R(i))
139  ENDDO
140  IF(ncsize.GT.1) THEN
141  c=p_min(c)
142  cinit=p_min(cinit)
143  ENDIF
144  WRITE(lu,*) 'AVANT TRAITEMENT HAUTEURS NEGATIVES, H MIN=',c
145  WRITE(lu,*) 'AVANT TRAITEMENT HAUTEURS NEGATIVES, HN MIN=',cinit
146  c=0.d0
147  cinit=0.d0
148  IF(ncsize.GT.1) THEN
149  DO i=1,npoin
150  c =c +h%R(i) *mesh%IFAC%I(i)/unsv2d%R(i)
151  cinit=cinit+hn%R(i)*mesh%IFAC%I(i)/unsv2d%R(i)
152  ENDDO
153  c =p_sum(c )
154  cinit=p_sum(cinit)
155  ELSE
156  DO i=1,npoin
157  c =c +h%R(i) /unsv2d%R(i)
158  cinit=cinit+hn%R(i)/unsv2d%R(i)
159  ENDDO
160  ENDIF
161  WRITE(lu,*) 'AVANT TRAITEMENT MASSE INITIALE=',cinit
162  WRITE(lu,*) 'AVANT TRAITEMENT MASSE FINALE =',c
163  ENDIF
164 !
165 ! CALCUL DES FLUX PAR SEGMENT (T1 SUIVI DE FALSE NON UTILISE)
166 ! WHEN RECEIVED FLODEL IS NOT ASSEMBLED IN //
167 !
168  IF(compute_flodel) THEN
169 ! SO FAR HARDCODED OPTION 2
170  CALL flux_ef_vf(flodel%R,flopoint,nseg,nelem,nelmax,
171  & mesh%ELTSEG%I,mesh%ORISEG%I,
172  & mesh%IKLE%I,.true.,2)
173  ENDIF
174 !
175 ! ASSEMBLING FLODEL IN PARALLEL
176 !
177  IF(ncsize.GT.1) THEN
178  CALL parcom2_seg(flodel%R,flodel%R,flodel%R,
179  & nseg,1,2,1,mesh,1,11)
180  ENDIF
181 !
182 ! CHANGING FLUXES FROM POINTS INTO N FLUXES BETWEEN POINTS
183  DO ielem = 1,nelem
184  a1 = abs(flopoint(ielem,1))
185  a2 = abs(flopoint(ielem,2))
186  a3 = abs(flopoint(ielem,3))
187  IF(a1.GE.a2.AND.a1.GE.a3) THEN
188 ! ALL FLOW TO AND FROM NODE 1
189  flopoint(ielem,1)=-flopoint(ielem,2)
190  flopoint(ielem,2)=0.d0
191 ! FLOPOINT(IELEM,3)= UNCHANGED!
192  ELSEIF(a2.GE.a1.AND.a2.GE.a3) THEN
193 ! ALL FLOW TO AND FROM NODE 2
194 ! FLOPOINT(IELEM,1)= UNCHANGED!
195  flopoint(ielem,2)=-flopoint(ielem,3)
196  flopoint(ielem,3)=0.d0
197  ELSE
198 ! ALL FLOW TO AND FROM NODE 3
199  flopoint(ielem,3)=-flopoint(ielem,1)
200  flopoint(ielem,1)=0.d0
201 ! FLOPOINT(IELEM,2)= UNCHANGED!
202  ENDIF
203  ENDDO
204 !
205 ! SAVING THE ORIGINAL FLOPOINT
206 !
207  IF(makeflulimebe) THEN
208  DO ielem=1,nelem
209  flulimebe(ielem,1)=flopoint(ielem,1)
210  flulimebe(ielem,2)=flopoint(ielem,2)
211  flulimebe(ielem,3)=flopoint(ielem,3)
212  ENDDO
213  ENDIF
214 !
215  CALL cpstvc(h,t2)
216  CALL cpstvc(h,t4)
217  CALL cpstvc(h,t1)
218 !
219 ! ADDING THE SOURCES (SMH IS NATURALLY ASSEMBLED IN //)
220 ! FIRST THE POSITIVE SOURCES
221 !
222  IF(yasmh) THEN
223  IF(optsou.EQ.1) THEN
224  DO i=1,npoin
225  h%R(i)=hn%R(i)+dt*max(smh%R(i),0.d0)
226  ENDDO
227  ELSEIF(optsou.EQ.2) THEN
228  DO i=1,npoin
229  h%R(i)=hn%R(i)+dt*max(smh%R(i),0.d0)*unsv2d%R(i)
230  ENDDO
231  ENDIF
232  ELSE
233  DO i=1,npoin
234  h%R(i)=hn%R(i)
235  ENDDO
236  ENDIF
237 !
238 ! RAIN (POSITIVE PLUIE)
239 !
240  IF(rain) THEN
241  DO i=1,npoin
242  h%R(i)=h%R(i)+dt*max(pluie%R(i),0.d0)
243  ENDDO
244  ENDIF
245 !
246 ! BOUNDARY FLUXES : ADDING THE ENTERING (NEGATIVE) FLUXES
247 ! FIRST PUTTING FLBOR (BOUNDARY) IN T2 (DOMAIN)
248 !
249 ! T2 WILL BE THE ASSEMBLED FLBOR, INITIALISATION HERE
250 ! IS USELESS EXCEPT THAT PARCOM MAY ADD UNDEFINED
251 ! NUMBERS (THAT WILL NOT BE USED BUT THAT WILL STOP
252 ! A COMPILER... TOO BAD!)
253  IF(ncsize.GT.1) CALL os('X=0 ',x=t2)
254  CALL osdb( 'X=Y ' ,t2,flbor,flbor,0.d0,mesh)
255 ! ASSEMBLING T2 (FLBOR IS NOT ASSEMBLED)
256  IF(ncsize.GT.1) CALL parcom(t2,2,mesh)
257  DO iptfr=1,nptfr
258  i=nbor(iptfr)
259  h%R(i)=h%R(i)-dt*unsv2d%R(i)*min(t2%R(i),0.d0)
260  ENDDO
261 !
262 ! FOR OPTIMIZING THE LOOP ON ELEMENTS, ONLY ELEMENTS
263 ! WITH NON ZERO FLUXES WILL BE CONSIDERED, THIS LIST
264 ! WILL BE UPDATED. TO START WITH, ALL ELEMENTS IN THE LIST
265 !
266  remain=nelem
267 !
268  DO i=1,remain
269  indic_pdept_eria(i)=i
270  ENDDO
271 !
272  cprev=0.d0
273 ! MAXIMUM INITIAL FLUX
274  DO ir=1,nelem
275  cprev=cprev+abs(flopoint(ir,1))
276  & +abs(flopoint(ir,2))
277  & +abs(flopoint(ir,3))
278  ENDDO
279  IF(ncsize.GT.1) cprev=p_sum(cprev)
280  IF(testing) WRITE(lu,*) 'INITIAL SUM OF FLUXES=',cprev
281  cinit=cprev
282 !
283 ! LOOP OVER THE LOOP OVER THE ELEMENTS
284 !
285  niter = 0
286 777 CONTINUE
287  niter = niter + 1
288 !
289 ! T4 IS THE EVOLUTION OF VOLUME, HERE INITIALISED TO 0
290  CALL os('X=0 ',x=t4)
291 !
292 ! LOOP OVER THE ELEMENTS
293 !
294  newremain=0
295  c=0.d0
296 !
297 ! COMPUTING DEMAND (T1) AND OFFER (T3)
298 !
299  CALL os('X=0 ',x=t1)
300  CALL os('X=0 ',x=t3)
301  IF(optpre.EQ.1) THEN
302  DO ir=1,remain
303  i=indic_pdept_eria(ir)
304  i1=mesh%IKLE%I(i )
305  i2=mesh%IKLE%I(i +nelmax)
306  i3=mesh%IKLE%I(i+2*nelmax)
307 ! A PRIORI AVAILABLE VOLUMES
308  vol1=mesh%SURFAC%R(i)*h%R(i1)*tiers
309  vol2=mesh%SURFAC%R(i)*h%R(i2)*tiers
310  vol3=mesh%SURFAC%R(i)*h%R(i3)*tiers
311 ! FLUXES FROM POINTS
312  f1= flopoint(i,1)-flopoint(i,3)
313  f2=-flopoint(i,1)+flopoint(i,2)
314  f3=-flopoint(i,2)+flopoint(i,3)
315 ! DEMAND OR OFFER, COMPARED TO A PRIORI AVAILABLE VOLUMES
316  IF(f1*dt.GT.vol1) THEN
317  t1%R(i1)=t1%R(i1)+f1*dt-vol1
318  ELSE
319  t3%R(i1)=t3%R(i1)+min(vol1,vol1-f1*dt)
320  ENDIF
321  IF(f2*dt.GT.vol2) THEN
322  t1%R(i2)=t1%R(i2)+f2*dt-vol2
323  ELSE
324  t3%R(i2)=t3%R(i2)+min(vol2,vol2-f2*dt)
325  ENDIF
326  IF(f3*dt.GT.vol3) THEN
327  t1%R(i3)=t1%R(i3)+f3*dt-vol3
328  ELSE
329  t3%R(i3)=t3%R(i3)+min(vol3,vol3-f3*dt)
330  ENDIF
331  ENDDO
332  ELSEIF(optpre.EQ.2) THEN
333  DO ir=1,remain
334  i=indic_pdept_eria(ir)
335  i1=mesh%IKLE%I(i )
336  i2=mesh%IKLE%I(i +nelmax)
337  i3=mesh%IKLE%I(i+2*nelmax)
338 ! A PRIORI AVAILABLE VOLUMES
339  vol1=mesh%SURFAC%R(i)*h%R(i1)*tiers
340  vol2=mesh%SURFAC%R(i)*h%R(i2)*tiers
341  vol3=mesh%SURFAC%R(i)*h%R(i3)*tiers
342 ! FLUXES FROM POINTS
343  f1= flopoint(i,1)-flopoint(i,3)
344  f2=-flopoint(i,1)+flopoint(i,2)
345  f3=-flopoint(i,2)+flopoint(i,3)
346 ! DEMAND AND OFFER
347  IF(f1.GT.0.d0) t1%R(i1)=t1%R(i1)+f1
348  IF(f2.GT.0.d0) t1%R(i2)=t1%R(i2)+f2
349  IF(f3.GT.0.d0) t1%R(i3)=t1%R(i3)+f3
350  t3%R(i1)=t3%R(i1)+vol1
351  t3%R(i2)=t3%R(i2)+vol2
352  t3%R(i3)=t3%R(i3)+vol3
353  ENDDO
354  ENDIF
355  IF(ncsize.GT.1) THEN
356  CALL parcom(t1,2,mesh)
357  CALL parcom(t3,2,mesh)
358  ENDIF
359 !
360  DO ir=1,remain
361 !
362  i=indic_pdept_eria(ir)
363 !
364  i1=mesh%IKLE%I(i )
365  i2=mesh%IKLE%I(i +nelmax)
366  i3=mesh%IKLE%I(i+2*nelmax)
367 !
368 ! A PRIORI AVAILABLE VOLUMES
369 !
370  vol1=mesh%SURFAC%R(i)*h%R(i1)*tiers
371  vol2=mesh%SURFAC%R(i)*h%R(i2)*tiers
372  vol3=mesh%SURFAC%R(i)*h%R(i3)*tiers
373 !
374 ! FLUXES FROM POINTS
375 !
376  f1= flopoint(i,1)-flopoint(i,3)
377  f2=-flopoint(i,1)+flopoint(i,2)
378  f3=-flopoint(i,2)+flopoint(i,3)
379 !
380 ! DISTRIBUTION OF VOLUMES ACCORDING TO DEMAND AND OFFER
381 !
382  IF(optpre.EQ.1) THEN
383  IF(f1*dt.GT.vol1) THEN
384  IF(t1%R(i1).GT.t3%R(i1)) THEN
385  vol1=vol1+(f1*dt-vol1)*(t3%R(i1)/t1%R(i1))
386  ELSE
387  vol1=f1*dt
388  ENDIF
389  ELSE
390  IF(t3%R(i1).GT.t1%R(i1)) THEN
391  vol1=vol1-min(vol1,vol1-f1*dt)*(t1%R(i1)/t3%R(i1))
392  ELSE
393  vol1=max(f1,0.d0)*dt
394  ENDIF
395  ENDIF
396  IF(f2*dt.GT.vol2) THEN
397  IF(t1%R(i2).GT.t3%R(i2)) THEN
398  vol2=vol2+(f2*dt-vol2)*(t3%R(i2)/t1%R(i2))
399  ELSE
400  vol2=f2*dt
401  ENDIF
402  ELSE
403  IF(t3%R(i2).GT.t1%R(i2)) THEN
404  vol2=vol2-min(vol2,vol2-f2*dt)*(t1%R(i2)/t3%R(i2))
405  ELSE
406  vol2=max(f2,0.d0)*dt
407  ENDIF
408  ENDIF
409  IF(f3*dt.GT.vol3) THEN
410  IF(t1%R(i3).GT.t3%R(i3)) THEN
411  vol3=vol3+(f3*dt-vol3)*(t3%R(i3)/t1%R(i3))
412  ELSE
413  vol3=f3*dt
414  ENDIF
415  ELSE
416  IF(t3%R(i3).GT.t1%R(i3)) THEN
417  vol3=vol3-min(vol3,vol3-f3*dt)*(t1%R(i3)/t3%R(i3))
418  ELSE
419  vol3=max(f3,0.d0)*dt
420  ENDIF
421  ENDIF
422  ELSEIF(optpre.EQ.2) THEN
423  IF(t1%R(i1).GT.1.d-30) THEN
424 ! ( THIS IS IMPORTANT
425  vol1=t3%R(i1)*(max(f1,0.d0)/t1%R(i1))
426  ELSE
427  vol1=mesh%SURFAC%R(i)*h%R(i1)*tiers
428  ENDIF
429  IF(t1%R(i2).GT.1.d-30) THEN
430  vol2=t3%R(i2)*(max(f2,0.d0)/t1%R(i2))
431  ELSE
432  vol2=mesh%SURFAC%R(i)*h%R(i2)*tiers
433  ENDIF
434  IF(t1%R(i3).GT.1.d-30) THEN
435  vol3=t3%R(i3)*(max(f3,0.d0)/t1%R(i3))
436  ELSE
437  vol3=mesh%SURFAC%R(i)*h%R(i3)*tiers
438  ENDIF
439  ENDIF
440 !
441 ! NOW LIMITATION OF FLUXES
442 !
443  IF(f1*dt.GT.vol1) THEN
444  dt1=dt*(vol1/(f1*dt))
445  ELSE
446  dt1=dt
447  ENDIF
448  IF(f2*dt.GT.vol2) THEN
449  dt2=dt*(vol2/(f2*dt))
450  ELSE
451  dt2=dt
452  ENDIF
453  IF(f3*dt.GT.vol3) THEN
454  dt3=dt*(vol3/(f3*dt))
455  ELSE
456  dt3=dt
457  ENDIF
458 !
459 ! LIMITED VOLUMES TRANSITING BETWEEN POINTS (1/DT MISSING)
460 !
461  dtlim1=min(dt1,dt2)
462  dtlim2=min(dt2,dt3)
463  dtlim3=min(dt3,dt1)
464 !
465  fp1=flopoint(i,1)*dtlim1
466  fp2=flopoint(i,2)*dtlim2
467  fp3=flopoint(i,3)*dtlim3
468 !
469 ! CORRESPONDING VARIATIONS OF VOLUMES OF POINTS (DT MISSING SO OK)
470 !
471  t4%R(i1)=t4%R(i1)-( fp1-fp3)
472  t4%R(i2)=t4%R(i2)-(-fp1+fp2)
473  t4%R(i3)=t4%R(i3)-(-fp2+fp3)
474 !
475 ! IF REMAINING FLUXES, THE ELEMENT IS KEPT IN THE LIST
476 !
477  IF(dtlim1.EQ.dt.AND.dtlim2.EQ.dt.AND.dtlim3.EQ.dt) THEN
478  flopoint(i,1)=0.d0
479  flopoint(i,2)=0.d0
480  flopoint(i,3)=0.d0
481  ELSE
482  newremain=newremain+1
483 ! BEFORE NEWREMAIN: FOR NEXT ITERATION
484 ! AFTER NEWREMAIN: STILL VALID FOR NEXT ITERATION
485  indic_pdept_eria(newremain)=i
486  flopoint(i,1)=flopoint(i,1)*(1.d0-dtlim1*surdt)
487  flopoint(i,2)=flopoint(i,2)*(1.d0-dtlim2*surdt)
488  flopoint(i,3)=flopoint(i,3)*(1.d0-dtlim3*surdt)
489  c=c+abs(flopoint(i,1))+abs(flopoint(i,2))
490  & +abs(flopoint(i,3))
491  ENDIF
492 !
493  ENDDO
494 !
495  remain=newremain
496 !
497 ! ADDING THE EVOLUTIONS TO THE DEPTHS
498 ! AFTER ASSEMBLY AT INTERFACES AND AFTER
499 ! CHANGING VOLUMES INTO DEPTHS.
500 !
501  IF(ncsize.GT.1) CALL parcom(t4,2,mesh)
502  CALL os('X=XY ',x=t4,y=unsv2d)
503  CALL os('X=X+Y ',x=h,y=t4)
504  DO i=1,h%DIM1
505  h%R(i)=max(h%R(i),0.d0)
506  ENDDO
507  IF(ncsize.GT.1) c=p_sum(c)
508  IF(testing) WRITE(lu,*) 'FLUX NON PRIS EN COMPTE=',c
509 !
510 ! STOP CRITERION
511 !
512  IF(c.NE.cprev.AND.abs(c-cprev).GT.cinit*1.d-9
513  & .AND.c.NE.0.d0) THEN
514  cprev=c
515  IF(niter.LT.nitmax) GO TO 777
516  ENDIF
517 !
518 ! ADDING THE SOURCES (SMH IS NATURALLY ASSEMBLED IN //)
519 ! NOW THE NEGATIVE SOURCES
520 !
521  IF(yasmh) THEN
522  IF(optsou.EQ.1) THEN
523  DO i=1,npoin
524  h%R(i)=h%R(i)+dt*min(smh%R(i),0.d0)
525  ENDDO
526  ELSEIF(optsou.EQ.2) THEN
527  DO i=1,npoin
528  h%R(i)=h%R(i)+dt*min(smh%R(i),0.d0)*unsv2d%R(i)
529  ENDDO
530  ENDIF
531  ENDIF
532 !
533 ! EVAPORATION (NEGATIVE PLUIE), WITH A POSSIBLE LIMITATION
534 !
535  IF(rain) THEN
536  DO i=1,npoin
537  IF(-dt*pluie%R(i).LT.h%R(i)) THEN
538  h%R(i)=h%R(i)+dt*min(pluie%R(i),0.d0)
539  ELSE
540  pluie%R(i)=-h%R(i)/dt
541  h%R(i)=0.d0
542  ENDIF
543  ENDDO
544  ENDIF
545 !
546 ! BOUNDARY FLUXES : ADDING THE EXITING (POSITIVE) FLUXES
547 ! WITH A POSSIBLE LIMITATION
548 ! IMPORTANT: MUST BE DONE AFTER ALL OTHER SOURCES
549 !
550  DO iptfr=1,nptfr
551  i=nbor(iptfr)
552 ! T2 = // ASSEMBLED FLBOR
553  hfl1=dt*unsv2d%R(i)*max(t2%R(i),0.d0)
554  IF(hfl1.GT.h%R(i)) THEN
555 ! FLBOR ACTUALLY TAKEN INTO ACCOUNT
556  flbor%R(iptfr)=flbor%R(iptfr)*h%R(i)/hfl1
557  h%R(i)=0.d0
558  ELSE
559  h%R(i)=h%R(i)-hfl1
560  ENDIF
561  IF(limpro(iptfr).EQ.kdir) THEN
562  IF(hbor(iptfr).LT.0.d0) THEN
563  WRITE(lu,*) 'NEGATIVE DEPTH PRESCRIBED ON BOUNDARY'
564  WRITE(lu,*) 'CHECK YOUR SPECIFIC SUBROUTINE:'
565  IF(namecode(1:7).EQ.'SISYPHE') THEN
566  WRITE(lu,*) 'BEDLOAD_SOLVS_FE'
567  ELSEIF(namecode(1:4).EQ.'GAIA') THEN
568  WRITE(lu,*) 'BEDLOAD_SOLVS_FE'
569  ELSEIF(namecode(1:9).EQ.'TELEMAC2D') THEN
570  WRITE(lu,*) 'BORD'
571  ELSEIF(namecode(1:9).EQ.'TELEMAC3D') THEN
572  WRITE(lu,*) 'BORD3D'
573  ENDIF
574  CALL plante(1)
575  stop
576  ENDIF
577 ! HERE V2DPAR WOULD BE MORE CONVENIENT THAN UNSV2D...
578  flbor%R(iptfr)=flbor%R(iptfr)
579  & +(h%R(i)-hbor(iptfr))/(dt*unsv2d%R(i))
580  h%R(i)= hbor(iptfr)
581  ENDIF
582  ENDDO
583 !
584  IF(testing) THEN
585  c=1.d99
586  DO i=1,npoin
587  c=min(c,h%R(i))
588  ENDDO
589  IF(ncsize.GT.1) c=p_min(c)
590  WRITE(lu,*) 'APRES TRAITEMENT HAUTEURS NEGATIVES, HMIN=',c
591  c=0.d0
592  IF(ncsize.GT.1) THEN
593  DO i=1,npoin
594  c=c+h%R(i)*mesh%IFAC%I(i)/unsv2d%R(i)
595  ENDDO
596  c=p_sum(c)
597  ELSE
598  DO i=1,npoin
599  c=c+h%R(i)/unsv2d%R(i)
600  ENDDO
601  ENDIF
602  WRITE(lu,*) 'APRES TRAITEMENT MASSE FINALE =',c
603  ENDIF
604 !
605 !-----------------------------------------------------------------------
606 !
607 ! NOW WE WANT:
608 ! FLULIM TO BE THE PERCENTAGE OF FLUX THAT HAS BEEN TRANSMITTED
609 ! FLODEL TO BE THE ACTUAL FLUX THAT HAS BEEN TRANSMITTED
610 !
611 ! WE START WITH FLULIM=THE ORIGINAL FLODEL
612 !
613 ! WORKING ON ASSEMBLED VALUES TO FIND AN AVERAGE FLULIM
614 ! THAT WILL BE THE SAME AT INTERFACES
615 !
616 ! ASSEMBLING THE REMAINING FLUXES IN FLODEL
617 !
618  DO i=1,nelem
619  IF(abs(flulimebe(i,1)).GT.1.d-30) THEN
620  flulimebe(i,1)=(flulimebe(i,1)-flopoint(i,1))
621  & /flulimebe(i,1)
622  ELSE
623  flulimebe(i,1)=0.d0
624  ENDIF
625  IF(abs(flulimebe(i,2)).GT.1.d-30) THEN
626  flulimebe(i,2)=(flulimebe(i,2)-flopoint(i,2))
627  & /flulimebe(i,2)
628  ELSE
629  flulimebe(i,2)=0.d0
630  ENDIF
631  IF(abs(flulimebe(i,3)).GT.1.d-30) THEN
632  flulimebe(i,3)=(flulimebe(i,3)-flopoint(i,3))
633  & /flulimebe(i,3)
634  ELSE
635  flulimebe(i,3)=0.d0
636  ENDIF
637  ENDDO
638 !
639 !-----------------------------------------------------------------------
640 !
641 ! CHECKING: COMPARING H(N+1) WITH H RECONSTRUCTED WITH THE FLUXES
642 ! SOURCES LACKING...
643 !
644  IF(testing) THEN
645  DO i=1,npoin
646  t1%R(i)=0.d0
647  ENDDO
648  DO i=1,nseg
649  i1=gloseg1(i)
650  i2=gloseg2(i)
651  t1%R(i1)=t1%R(i1)-dt*unsv2d%R(i1)*flodel%R(i)
652  t1%R(i2)=t1%R(i2)+dt*unsv2d%R(i2)*flodel%R(i)
653  ENDDO
654  DO iptfr=1,nptfr
655  i=nbor(iptfr)
656  t1%R(i)=t1%R(i)-dt*unsv2d%R(i)*flbor%R(iptfr)
657  ENDDO
658  IF(ncsize.GT.1) CALL parcom(t1,2,mesh)
659  DO i=1,npoin
660  t1%R(i)=t1%R(i)+hn%R(i)-h%R(i)
661  ENDDO
662  WRITE(lu,*) 'ERREUR POSITIVE_DEPTHS=',p_dots(t1,t1,mesh)
663  ENDIF
664 !
665 !-----------------------------------------------------------------------
666 !
667  IF(info) THEN
668  IF(namecode(1:7).EQ.'SISYPHE'.OR.namecode(1:4).EQ.'GAIA') THEN
669  WRITE(lu,102) niter
670  ELSE
671  WRITE(lu,202) niter
672  ENDIF
673  ENDIF
674 !
675  IF(niter.EQ.nitmax) THEN
676  IF(namecode(1:7).EQ.'SISYPHE'.OR.namecode(1:4).EQ.'GAIA') THEN
677  WRITE(lu,103) niter
678  ELSE
679  WRITE(lu,203) niter
680  ENDIF
681  ENDIF
682 !
683 102 FORMAT(' BEDLOAD EQUATION SOLVED IN ',1i5,' ITERATIONS')
684 202 FORMAT(' POSITIVE DEPTHS OBTAINED IN ',1i5,' ITERATIONS')
685 103 FORMAT(' BEDLOAD EQUATION SOLVED IN ',1i5,' ITERATIONS = MAXIMUM',
686  & /,'INCREASE MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES')
687 203 FORMAT(' POSITIVE DEPTHS SOLVED IN ',1i5,' ITERATIONS = MAXIMUM',
688  & /,'INCREASE MAXIMUM NUMBER OF ITERATIONS FOR ADVECTION SCHEMES')
689 !
690 !-----------------------------------------------------------------------
691 !
692  RETURN
693  END
interface_parallel
Definition: interface_parallel.f:3
declarations_special
Definition: declarations_special.F:3
positive_depths_eria
subroutine positive_depths_eria(T1, T2, T3, T4, H, HN, MESH, FLODEL, COMPUTE_FLODEL, FLBOR, DT, UNSV2D, NPOIN, GLOSEG1, GLOSEG2, NBOR, NPTFR, SMH, YASMH, PLUIE, RAIN, OPTSOU, LIMPRO, HBOR, KDIR, INFO, FLOPOINT, NAMECODE, NITMAX, MAKEFLULIMEBE, FLULIMEBE)
Definition: positive_depths_eria.f:10
parcom2_seg
subroutine parcom2_seg(X1, X2, X3, NSEG, NPLAN, ICOM, IAN, MESH, OPT, IELM)
Definition: parcom2_seg.f:7
declarations_special::lu
integer lu
Definition: declarations_special.F:45
declarations_telemac
Definition: declarations_telemac.f:3
interface_parallel::p_sum
Definition: interface_parallel.f:247
declarations_telemac::deja_pdept_eria
logical deja_pdept_eria
Definition: declarations_telemac.f:286
interface_parallel::p_min
Definition: interface_parallel.f:235
cpstvc
subroutine cpstvc(X, Y)
Definition: cpstvc.f:7
osdb
subroutine osdb(OP, X, Y, Z, C, MESH)
Definition: osdb.f:7
flux_ef_vf
subroutine flux_ef_vf(FLOW, PHIEL, NSEG, NELEM, NELMAX, ELTSEG, ORISEG, IKLE, INIFLO, IOPT, FN, YAFLULIM, FLULIM, YAFLULIMEBE, FLULIMEBE)
Definition: flux_ef_vf.f:8
os
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
parcom
subroutine parcom(X, ICOM, MESH)
Definition: parcom.f:7
declarations_telemac::indic_pdept_eria
integer, dimension(:), allocatable indic_pdept_eria
Definition: declarations_telemac.f:287
p_dots
double precision function p_dots(X, Y, MESH)
Definition: p_dots.f:7
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:37:24 for The TELEMAC-MASCARET system by doxygen 1.8.13