The TELEMAC-MASCARET system  trunk
cflvf.f
Go to the documentation of this file.
1 ! ****************
2  SUBROUTINE cflvf
3 ! ****************
4 !
5  &(dtmax,hstart,fxmat,fxmatpar,mas,dt,fxbor,smh,yasmh,tab1,nseg,
6  & npoin,nptfr,gloseg,sizglo,mesh,msk,maskpt,rain,pluie,fc,
7  & nelem,ikle,limtra,kdir,kddl,fbor,fscexp,train,nbor,minfc,maxfc,
8  & secu,opt,withabs)
9 !
10 !***********************************************************************
11 ! BIEF V7P3
12 !***********************************************************************
13 !
14 !brief COMPUTES THE MAXIMUM TIMESTEP THAT ENABLES
15 !+ MONOTONY OF THE N ADVECTION SCHEME.
16 !
17 !history JMH
18 !+ 11/04/2008
19 !+
20 !+ ADDED YASMH
21 !
22 !history JMH
23 !+ 10/06/2008
24 !+
25 !+ ADDED SIZGLO
26 !
27 !history JMH
28 !+ 02/10/2008
29 !+
30 !+ PARALLEL MODE (ADDED FXMATPAR, ETC.)
31 !
32 !history C-T PHAM (LNHE)
33 !+ 30/11/2009
34 !+ V6P0
35 !+ REFINED COMPUTATION OF DTMAX (AS IN 3D)
36 !
37 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
38 !+ 13/07/2010
39 !+ V6P0
40 !+ Translation of French comments within the FORTRAN sources into
41 !+ English comments
42 !
43 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
44 !+ 21/08/2010
45 !+ V6P0
46 !+ Creation of DOXYGEN tags for automated documentation and
47 !+ cross-referencing of the FORTRAN sources
48 !
49 !history J-M HERVOUET (LNHE)
50 !+ 24/02/2012
51 !+ V6P2
52 !+ Rain and evaporation added (after initiative by O. Boutron, from
53 !+ Tour du Valat, and O. Bertrand, Artelia group)
54 !
55 !history S. PAVAN & J-M HERVOUET (EDF R&D, LNHE)
56 !+ 06/06/2013
57 !+ V6P3
58 !+ New stability criterion based on local min and max of function.
59 !+ See hardcoded option OPT.
60 !
61 !history S. PAVAN & J-M HERVOUET (EDF LAB, LNHE)
62 !+ 24/10/2013
63 !+ V7P0
64 !+ Old stability criterion simplified.
65 !
66 !history S. PAVAN & J-M HERVOUET (EDF LAB, LNHE)
67 !+ 29/04/2014
68 !+ V7P0
69 !+ Security coefficient added (for predictor-corrector scheme).
70 !
71 !history S. PAVAN & J-M HERVOUET (EDF LAB, LNHE)
72 !+ 13/10/2014
73 !+ V7P0
74 !+ More general formula with parameters COEMIN and COESOU to account
75 !+ for predictor-corrector scheme.
76 !
77 !history S. PAVAN & J-M HERVOUET (EDF LAB, LNHE)
78 !+ 16/06/2015
79 !+ V7P1
80 !+ Option OPT added in argument. Implementation of OPT=2 changed.
81 !
82 !history J-M HERVOUET (EDF LAB, LNHE)
83 !+ 27/09/2016
84 !+ V7P2
85 !+ Coefficients COEMIN and COESOU non longer useful. Stability is
86 !+ done here only for the N scheme.
87 !
88 !history J-M HERVOUET (EDF LAB, LNHE)
89 !+ 29/09/2017
90 !+ V7P3
91 !+ Argument WITHABS added to have a new computation of the stability
92 !+ condition for all predictor-corrector based schemes. This corrects
93 !+ a possible hack in monotonicity (which never occurred in tests).
94 !
95 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
96 !| DT |-->| TIME STEP
97 !| DTMAX |<--| MAXIMUM TIME STEP FOR STABILITY
98 !| FXBOR |-->| BOUNDARY FLUXES
99 !| FXMAT |-->| FLUXES
100 !| FXMATPAR |-->| FLUXES ASSEMBLED IN PARALLEL
101 !| GLOSEG |-->| GLOBAL NUMBER OF THE 2 POINTS OF A SEGMENT
102 !| HSTART |-->| H AT BEGINNING OF SUB TIME STEP
103 !| MAS |-->| INTEGRAL OF TEST FUNCTIONS (=AREA AROUND POINTS)
104 !| MASKPT |-->| ARRAY FOR MASKING POINTS
105 !| MESH |-->| MESH STRUCTURE
106 !| MSK |-->| IF YES, THERE IS MASKED ELEMENTS.
107 !| NPOIN |-->| NUMBER OF POINTS IN THE MESH
108 !| NPTFR |-->| NUMBER OF BOUNDARY POINTS
109 !| NSEG |-->| NUMBER OF SEGMENTS
110 !| OPT |-->| OPTION: 1=CLASSIC 2=BASED ON MIN AND MAX
111 !| | | 3=BASED ON GIVEN MIN AND MAX
112 !| PLUIE |-->| RAIN OR EVAPORATION IN M/S
113 !| RAIN |-->| IF YES: RAIN OR EVAPORATION
114 !| SIZGLO |-->| FIRST DIMENSION OF GLOSEG
115 !| SMH |-->| RIGHT HAND SIDE OF CONTINUITY EQUATION
116 !| TAB1 |-->| WORK ARRAY
117 !| WITHABS |-->| IF YES, STABILITY COMPUTED WITH ABSOLUTE VALUE
118 !| | | OF FLUXES.
119 !| YASMH |-->| IF YES, TAKE SHM INTO ACCOUNT
120 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
121 !
122  USE bief, ex_cflvf => cflvf
123 !
124  USE declarations_special
125  IMPLICIT NONE
126 !
127 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
128 !
129  INTEGER, INTENT(IN) :: NSEG,NPOIN,NPTFR,SIZGLO,OPT
130  INTEGER, INTENT(IN) :: GLOSEG(sizglo,2),NELEM
131  INTEGER, INTENT(IN) :: IKLE(nelem,3),NBOR(nptfr)
132  INTEGER, INTENT(IN) :: LIMTRA(nptfr),KDIR,KDDL
133  DOUBLE PRECISION, INTENT(INOUT) :: DTMAX
134  DOUBLE PRECISION, INTENT(IN) :: DT,HSTART(npoin)
135  DOUBLE PRECISION, INTENT(IN) :: MAS(npoin),SMH(*)
136  DOUBLE PRECISION, INTENT(IN) :: PLUIE(*)
137 ! NOT NPTFR, SEE TRACVF
138  DOUBLE PRECISION, INTENT(IN) :: FXBOR(npoin)
139  DOUBLE PRECISION, INTENT(IN) :: FXMAT(nseg),FXMATPAR(nseg)
140  DOUBLE PRECISION, INTENT(INOUT) :: FC(npoin)
141  DOUBLE PRECISION, INTENT(IN) :: FSCEXP(*)
142  DOUBLE PRECISION, INTENT(IN) :: FBOR(nptfr),TRAIN,SECU
143  LOGICAL, INTENT(IN) :: YASMH,MSK,RAIN,WITHABS
144  TYPE(bief_obj), INTENT(INOUT) :: TAB1,MINFC,MAXFC
145  TYPE(bief_mesh), INTENT(INOUT) :: MESH
146  TYPE(bief_obj), INTENT(IN) :: MASKPT
147 !
148 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
149 !
150  INTEGER I,J,IELEM,I1,I2,I3,N,ISEG
151  DOUBLE PRECISION DI,MINEL,MAXEL
152 !
153 ! HARDCODED OPTION (1: CLASSICAL 2: BASED ON LOCAL MIN AND MAX
154 ! 3: MIN AND MAX GIVEN)
155 ! OPT=1
156 !
157 !-----------------------------------------------------------------------
158 ! DETERMINES MAX AND MIN FOR THE NEW MONOTONICITY CRITERION
159 !-----------------------------------------------------------------------
160 !
161  IF(opt.EQ.2.OR.opt.EQ.3) THEN
162  CALL cpstvc(maskpt,minfc)
163  CALL cpstvc(maskpt,maxfc)
164  ENDIF
165 !
166  IF(opt.EQ.3) THEN
167 !
168 ! HARDCODED MIN AND MAX
169 !
170  DO i=1,npoin
171  minfc%R(i)=0.d0
172  maxfc%R(i)=1000.d0
173 ! CLIPPING IS THERE NECESSARY
174  fc(i)=max(fc(i),minfc%R(i))
175  fc(i)=min(fc(i),maxfc%R(i))
176  ENDDO
177 !
178  ELSEIF(opt.EQ.2) THEN
179 !
180  DO i=1,npoin
181  minfc%R(i)=fc(i)
182  maxfc%R(i)=fc(i)
183  ENDDO
184 !
185  DO i=1,nptfr
186  IF(limtra(i).EQ.kdir) THEN
187  n=nbor(i)
188  minfc%R(n)=min(minfc%R(n),fbor(i))
189  maxfc%R(n)=max(maxfc%R(n),fbor(i))
190  ENDIF
191  ENDDO
192 !
193  DO ielem=1,nelem
194  i1=ikle(ielem,1)
195  i2=ikle(ielem,2)
196  i3=ikle(ielem,3)
197  minel=min(fc(i1),fc(i2),fc(i3))
198  maxel=max(fc(i1),fc(i2),fc(i3))
199  minfc%R(i1)=min(minfc%R(i1),minel)
200  minfc%R(i2)=min(minfc%R(i2),minel)
201  minfc%R(i3)=min(minfc%R(i3),minel)
202  maxfc%R(i1)=max(maxfc%R(i1),maxel)
203  maxfc%R(i2)=max(maxfc%R(i2),maxel)
204  maxfc%R(i3)=max(maxfc%R(i3),maxel)
205  ENDDO
206 !
207  IF(yasmh.AND.rain) THEN
208  DO i=1,npoin
209  minfc%R(i)=min(minfc%R(i),fscexp(i),train)
210  maxfc%R(i)=max(maxfc%R(i),fscexp(i),train)
211  ENDDO
212  ELSEIF(yasmh) THEN
213  DO i=1,npoin
214  minfc%R(i)=min(minfc%R(i),fscexp(i))
215  maxfc%R(i)=max(maxfc%R(i),fscexp(i))
216  ENDDO
217  ELSEIF(rain) THEN
218  DO i=1,npoin
219  minfc%R(i)=min(minfc%R(i),train)
220  maxfc%R(i)=max(maxfc%R(i),train)
221  ENDDO
222  ENDIF
223 !
224 ! IN PARALLEL MODE: GLOBAL EXTREMA
225 !
226  IF(ncsize.GT.1) THEN
227 ! GLOBAL EXTREMA
228  CALL parcom(maxfc,3,mesh)
229  CALL parcom(minfc,4,mesh)
230  ENDIF
231 !
232  ENDIF
233 !
234 !---------------------------------------------------------------------
235 !
236 ! COMPUTES THE CRITERION FOR COURANT NUMBER
237 !
238 ! USES HERE FXMAT ASSEMBLED IN PARALLEL FOR UPWINDING
239 !
240  IF(opt.EQ.1) THEN
241 !
242 ! CLASSICAL CRITERION, WITH ALTERNATE OPTION IF WITHABS=.TRUE.
243 ! WITHABS=.FALSE. : MAX(FLUX,0.D0) TAKEN, OR -MIN(FLUX,0.D0)
244 ! WITHABS=.TRUE. : ABS(FLUX) TAKEN (FOR PREDICTOR-CORRECTOR AND LIPS)
245 !
246  DO i = 1,npoin
247  tab1%R(i) = 0.d0
248  ENDDO
249 !
250  IF(withabs) THEN
251  DO i = 1,nseg
252 ! HERE WILL BE WRONG IN // IF THE 2 FLUXES OF A SAME SEGMENT
253 ! ARE CROSSING. THIS IS FORBIDDEN ANYWAY.
254  tab1%R(gloseg(i,1)) = tab1%R(gloseg(i,1)) + abs(fxmat(i))
255  tab1%R(gloseg(i,2)) = tab1%R(gloseg(i,2)) + abs(fxmat(i))
256  ENDDO
257  ELSE
258  DO i = 1,nseg
259  IF(fxmatpar(i).GT.0.d0) THEN
260 ! FXMAT(I) POSITIVE = MAX(...,0.D0) FOR 1
261  tab1%R(gloseg(i,1)) = tab1%R(gloseg(i,1)) + fxmat(i)
262  ELSE
263 ! - FXMAT(I) POSITIVE = MAX(...,0.D0) FOR 2
264  tab1%R(gloseg(i,2)) = tab1%R(gloseg(i,2)) - fxmat(i)
265  ENDIF
266  ENDDO
267  ENDIF
268 !
269  IF(ncsize.GT.1) CALL parcom(tab1,2,mesh)
270 !
271 ! MASKS TAB1
272  IF(msk) THEN
273  CALL os('X=XY ',x=tab1,y=maskpt)
274  ENDIF
275 !
276 ! STABILITY (AND MONOTONICITY) CRITERION
277  dtmax = dt
278 !
279 ! SEE RELEASE NOTES 5.7, CRITERION AT THE END OF 4.4 PAGE 33
280 ! BUT HERE THE FINAL H IS NOT H(N+1) BUT A FUNCTION OF DTMAX ITSELF
281 ! H FINAL = HSTART + DTMAX/DT *(H-HSTART)
282 !
283  IF(yasmh.AND.rain) THEN
284  IF(withabs) THEN
285  DO i = 1,npoin
286  tab1%R(i)=tab1%R(i)+abs(fxbor(i))
287  & +abs(smh(i)+pluie(i))
288  ENDDO
289  ELSE
290  DO i = 1,npoin
291  tab1%R(i)=tab1%R(i)+max(fxbor(i),0.d0)
292  & -min(smh(i)+pluie(i),0.d0)
293  ENDDO
294  ENDIF
295  ELSEIF(yasmh) THEN
296  IF(withabs) THEN
297  DO i = 1,npoin
298  tab1%R(i)=tab1%R(i)+abs(fxbor(i))+abs(smh(i))
299  ENDDO
300  ELSE
301  DO i = 1,npoin
302  tab1%R(i)=tab1%R(i)+max(fxbor(i),0.d0)-min(smh(i),0.d0)
303  ENDDO
304  ENDIF
305  ELSEIF(rain) THEN
306  IF(withabs) THEN
307  DO i = 1,npoin
308  tab1%R(i)=tab1%R(i)+abs(fxbor(i))+abs(pluie(i))
309  ENDDO
310  ELSE
311  DO i = 1,npoin
312  tab1%R(i)=tab1%R(i)+max(fxbor(i),0.d0)-min(pluie(i),0.d0)
313  ENDDO
314  ENDIF
315  ELSE
316  IF(withabs) THEN
317  DO i = 1,npoin
318  tab1%R(i)=tab1%R(i)+abs(fxbor(i))
319  ENDDO
320  ELSE
321  DO i = 1,npoin
322  tab1%R(i)=tab1%R(i)+max(fxbor(i),0.d0)
323  ENDDO
324  ENDIF
325  ENDIF
326 !
327  DO i = 1,npoin
328  IF(tab1%R(i).GT.1.d-20) THEN
329  minfc%R(i)=secu*mas(i)*hstart(i)/tab1%R(i)
330  dtmax=min(dtmax,minfc%R(i))
331  ELSE
332  minfc%R(i)=dt
333  ENDIF
334  ENDDO
335 !
336  ELSEIF(opt.EQ.2.OR.opt.EQ.3) THEN
337 !
338  dtmax = dt
339 !
340 ! NEW CRITERION, COMPUTING BI (SEE RELEASE NOTES 7.1)
341 !
342 ! WITH MAXFC%R
343 !
344  DO i = 1,npoin
345  tab1%R(i) = 0.d0
346  ENDDO
347 !
348  DO iseg=1,nseg
349  i=gloseg(iseg,1)
350  j=gloseg(iseg,2)
351  IF(fxmatpar(iseg).LT.0.d0) THEN
352  tab1%R(i)=tab1%R(i)
353  & -fxmat(iseg)*(fc(j)-maxfc%R(i))
354  tab1%R(j)=tab1%R(j)
355  & +fxmat(iseg)*(fc(j)-maxfc%R(j))
356  ELSE
357  tab1%R(i)=tab1%R(i)
358  & -fxmat(iseg)*(fc(i)-maxfc%R(i))
359  tab1%R(j)=tab1%R(j)
360  & +fxmat(iseg)*(fc(i)-maxfc%R(j))
361  ENDIF
362  ENDDO
363 !
364  DO i=1,nptfr
365  IF(limtra(i).EQ.kdir) THEN
366  n=nbor(i)
367 ! FXBOR IS HERE IN GLOBAL NUMBERING
368  tab1%R(n)=tab1%R(n)
369  & -min(fxbor(n),0.d0)*(fbor(i)-maxfc%R(n))
370  ELSEIF(limtra(i).EQ.kddl) THEN
371  n=nbor(i)
372  tab1%R(n)=tab1%R(n)
373  & -max(fxbor(n),0.d0)*(fc(n) -maxfc%R(n))
374  ENDIF
375  ENDDO
376 !
377  IF(yasmh) THEN
378  DO i = 1,npoin
379  tab1%R(i)=tab1%R(i)
380  & +max(smh(i),0.d0)*(fscexp(i)-maxfc%R(i))
381  & +min(smh(i),0.d0)*(fc(i) -maxfc%R(i))
382  ENDDO
383  ENDIF
384 !
385  IF(rain) THEN
386  DO i = 1,npoin
387  tab1%R(i)=tab1%R(i)
388  & +max(pluie(i),0.d0)*(train-maxfc%R(i))
389  & +min(pluie(i),0.d0)*(fc(i)-maxfc%R(i))
390  ENDDO
391  ENDIF
392 !
393  IF(ncsize.GT.1) CALL parcom(tab1,2,mesh)
394 !
395 ! MASKS TAB1
396 !
397  IF(msk) CALL os('X=XY ',x=tab1,y=maskpt)
398 !
399  DO i= 1,npoin
400  di=tab1%R(i)/mas(i)
401  IF(di.GT.1.d-12) THEN
402  dtmax = min(dtmax,hstart(i)*(maxfc%R(i)-fc(i))/di)
403  ENDIF
404  ENDDO
405 !
406 ! WITH MINFC%R
407 !
408  DO i = 1,npoin
409  tab1%R(i) = 0.d0
410  ENDDO
411 !
412  DO iseg=1,nseg
413  i=gloseg(iseg,1)
414  j=gloseg(iseg,2)
415  IF(fxmatpar(iseg).LT.0.d0) THEN
416  tab1%R(i)=tab1%R(i)
417  & -fxmat(iseg)*(fc(j)-minfc%R(i))
418  tab1%R(j)=tab1%R(j)
419  & +fxmat(iseg)*(fc(j)-minfc%R(j))
420  ELSE
421  tab1%R(i)=tab1%R(i)
422  & -fxmat(iseg)*(fc(i)-minfc%R(i))
423  tab1%R(j)=tab1%R(j)
424  & +fxmat(iseg)*(fc(i)-minfc%R(j))
425  ENDIF
426  ENDDO
427 !
428  DO i=1,nptfr
429  IF(limtra(i).EQ.kdir) THEN
430  n=nbor(i)
431 ! FXBOR IS HERE IN GLOBAL NUMBERING
432  tab1%R(n)=tab1%R(n)
433  & -min(fxbor(n),0.d0)*(fbor(i)-minfc%R(n))
434  ELSEIF(limtra(i).EQ.kddl) THEN
435  n=nbor(i)
436  tab1%R(n)=tab1%R(n)
437  & -max(fxbor(n),0.d0)*(fc(n) -minfc%R(n))
438  ENDIF
439  ENDDO
440 !
441  IF(yasmh) THEN
442  DO i = 1,npoin
443  tab1%R(i)=tab1%R(i)
444  & +max(smh(i),0.d0)*(fscexp(i)-minfc%R(i))
445  & +min(smh(i),0.d0)*(fc(i) -minfc%R(i))
446  ENDDO
447  ENDIF
448 !
449  IF(rain) THEN
450  DO i = 1,npoin
451  tab1%R(i)=tab1%R(i)
452  & +max(pluie(i),0.d0)*(train-minfc%R(i))
453  & +min(pluie(i),0.d0)*(fc(i)-minfc%R(i))
454  ENDDO
455  ENDIF
456 !
457  IF(ncsize.GT.1) CALL parcom(tab1,2,mesh)
458 !
459 ! MASKS TAB1
460 !
461  IF(msk) CALL os('X=XY ',x=tab1,y=maskpt)
462 !
463  DO i= 1,npoin
464  di=tab1%R(i)/mas(i)
465  IF(di.LT.-1.d-12) THEN
466  dtmax = min(dtmax,hstart(i)*(minfc%R(i)-fc(i))/di)
467  ENDIF
468  ENDDO
469 !
470  ELSE
471 !
472  WRITE(lu,*) 'UNKNOWN OPTION IN CFLVF: ',opt
473  CALL plante(1)
474  stop
475 !
476  ENDIF
477 !
478 !-----------------------------------------------------------------------
479 !
480  RETURN
481  END
482 
declarations_special
Definition: declarations_special.F:3
declarations_special::lu
integer lu
Definition: declarations_special.F:45
cflvf
subroutine cflvf(DTMAX, HSTART, FXMAT, FXMATPAR, MAS, DT, FXBOR, SMH, YASMH, TAB1, NSEG, NPOIN, NPTFR, GLOSEG, SIZGLO, MESH, MSK, MASKPT, RAIN, PLUIE, FC, NELEM, IKLE, LIMTRA, KDIR, KDDL, FBOR, FSCEXP, TRAIN, NBOR, MINFC, MAXFC, SECU, OPT, WITHABS)
Definition: cflvf.f:10
cpstvc
subroutine cpstvc(X, Y)
Definition: cpstvc.f:7
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
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:30 for The TELEMAC-MASCARET system by doxygen 1.8.13