cflvf.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\cflvf.f
00002 !
00105                      SUBROUTINE CFLVF
00106 !                    ****************
00107 !
00108      &(DTMAX,HSTART,H,FXMAT,FXMATPAR,MAS,DT,FXBOR,SMH,YASMH,TAB1,NSEG,
00109      & NPOIN,NPTFR,GLOSEG,SIZGLO,MESH,MSK,MASKPT,RAIN,PLUIE,FC,
00110      & NELEM,IKLE,LIMTRA,KDIR,FBOR,FSCEXP,TRAIN,NBOR,MINFC,MAXFC,
00111      & SECU,COEMIN,COESOU)
00112 !
00113 !***********************************************************************
00114 ! BIEF   V7P0
00115 !***********************************************************************
00116 !
00117 !
00118 !
00119 !
00120 !
00121 !
00122 !
00123 !
00124 !
00125 !
00126 !
00127 !
00128 !
00129 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00130 !| COEMIN         |-->| COEFFICIENT IN THE STABILITY CONDITION
00131 !| COESOU         |-->| COEFFICIENT IN THE STABILITY CONDITION
00132 !| DT             |-->| TIME STEP
00133 !| DTMAX          |<--| MAXIMUM TIME STEP FOR STABILITY
00134 !| FXBOR          |-->| BOUNDARY FLUXES
00135 !| FXMAT          |-->| FLUXES
00136 !| FXMATPAR       |-->| FLUXES ASSEMBLED IN PARALLEL
00137 !| GLOSEG         |-->| GLOBAL NUMBER OF THE 2 POINTS OF A SEGMENT
00138 !| H              |-->| H AT THE END OF FULL TIME STEP
00139 !| HSTART         |-->| H AT BEGINNING OF SUB TIME STEP
00140 !| MAS            |-->| INTEGRAL OF TEST FUNCTIONS (=AREA AROUND POINTS)
00141 !| MASKPT         |-->| ARRAY FOR MASKING POINTS
00142 !| MESH           |-->| MESH STRUCTURE
00143 !| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
00144 !| NPOIN          |-->| NUMBER OF POINTS IN THE MESH
00145 !| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
00146 !| NSEG           |-->| NUMBER OF SEGMENTS
00147 !| PLUIE          |-->| RAIN OR EVAPORATION IN M/S
00148 !| RAIN           |-->| IF YES: RAIN OR EVAPORATION
00149 !| SIZGLO         |-->| FIRST DIMENSION OF GLOSEG
00150 !| SMH            |-->| RIGHT HAND SIDE OF CONTINUITY EQUATION
00151 !| TAB1           |-->| WORK ARRAY
00152 !| YASMH          |-->| IF YES, TAKE SHM INTO ACCOUNT
00153 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00154 !
00155       USE BIEF, EX_CFLVF => CFLVF
00156 !
00157       IMPLICIT NONE
00158       INTEGER LNG,LU
00159       COMMON/INFO/LNG,LU
00160 !
00161 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00162 !
00163       INTEGER, INTENT(IN)             :: NSEG,NPOIN,NPTFR,SIZGLO
00164       INTEGER, INTENT(IN)             :: GLOSEG(SIZGLO,2),NELEM
00165       INTEGER, INTENT(IN)             :: IKLE(NELEM,3),NBOR(NPTFR)
00166       INTEGER, INTENT(IN)             :: LIMTRA(NPTFR),KDIR
00167       DOUBLE PRECISION, INTENT(INOUT) :: DTMAX
00168       DOUBLE PRECISION, INTENT(IN)    :: DT,HSTART(NPOIN)
00169       DOUBLE PRECISION, INTENT(IN)    :: H(NPOIN),MAS(NPOIN),SMH(NPOIN)
00170       DOUBLE PRECISION, INTENT(IN)    :: PLUIE(NPOIN)
00171 !                                              NOT NPTFR, SEE TRACVF
00172       DOUBLE PRECISION, INTENT(IN)    :: FXBOR(NPOIN)
00173       DOUBLE PRECISION, INTENT(IN)    :: FXMAT(NSEG),FXMATPAR(NSEG)
00174       DOUBLE PRECISION, INTENT(INOUT) :: FC(NPOIN)
00175       DOUBLE PRECISION, INTENT(IN)    :: FSCEXP(NPOIN)
00176       DOUBLE PRECISION, INTENT(IN)    :: FBOR(NPTFR),TRAIN,SECU
00177       DOUBLE PRECISION, INTENT(IN)    :: COEMIN,COESOU
00178       LOGICAL, INTENT(IN)             :: YASMH,MSK,RAIN
00179       TYPE(BIEF_OBJ), INTENT(INOUT)   :: TAB1,MINFC,MAXFC
00180       TYPE(BIEF_MESH), INTENT(INOUT)  :: MESH
00181       TYPE(BIEF_OBJ), INTENT(IN)      :: MASKPT
00182 !
00183 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00184 !
00185       INTEGER I,IELEM,I1,I2,I3,OPT,N
00186       DOUBLE PRECISION DENOM,A,B,DIFF,DI,MINEL,MAXEL
00187 !
00188 !     HARDCODED OPTION (1: CLASSICAL 2: BASED ON LOCAL MIN AND MAX
00189 !                                    3: MIN AND MAX GIVEN)
00190       OPT=1
00191 !
00192 !-----------------------------------------------------------------------
00193 ! DETERMINES MAX AND MIN FOR THE NEW MONOTONICITY CRITERION
00194 !-----------------------------------------------------------------------
00195 !
00196       IF(OPT.EQ.3) THEN
00197 !
00198 !       HARDCODED MIN AND MAX
00199 !
00200         CALL CPSTVC(MASKPT,MINFC)
00201         CALL CPSTVC(MASKPT,MAXFC)
00202         DO I=1,NPOIN
00203           MINFC%R(I)=0.D0
00204           MAXFC%R(I)=1000.D0
00205 !         CLIPPING IS THERE NECESSARY
00206           FC(I)=MAX(FC(I),MINFC%R(I))
00207           FC(I)=MIN(FC(I),MAXFC%R(I))
00208         ENDDO
00209 !
00210       ELSEIF(OPT.EQ.2) THEN
00211 !
00212         DO I=1,NPOIN
00213           MINFC%R(I)=FC(I)
00214           MAXFC%R(I)=FC(I)
00215         ENDDO
00216 !
00217         DO I=1,NPTFR
00218           IF(LIMTRA(I).EQ.KDIR) THEN
00219             N=NBOR(I)
00220             MINFC%R(N)=MIN(MINFC%R(N),FBOR(I))
00221             MAXFC%R(N)=MAX(MAXFC%R(N),FBOR(I))
00222           ENDIF
00223         ENDDO
00224 !
00225         DO IELEM=1,NELEM
00226           I1=IKLE(IELEM,1)
00227           I2=IKLE(IELEM,2)
00228           I3=IKLE(IELEM,3)
00229           MINEL=MIN(FC(I1),FC(I2),FC(I3))
00230           MAXEL=MAX(FC(I1),FC(I2),FC(I3))
00231           MINFC%R(I1)=MIN(MINFC%R(I1),MINEL)
00232           MINFC%R(I2)=MIN(MINFC%R(I2),MINEL)
00233           MINFC%R(I3)=MIN(MINFC%R(I3),MINEL)
00234           MAXFC%R(I1)=MAX(MAXFC%R(I1),MAXEL)
00235           MAXFC%R(I2)=MAX(MAXFC%R(I2),MAXEL)
00236           MAXFC%R(I3)=MAX(MAXFC%R(I3),MAXEL)
00237         ENDDO
00238 !
00239         IF(YASMH.AND.RAIN) THEN
00240           DO I=1,NPOIN
00241             MINFC%R(I)=MIN(MINFC%R(I),FSCEXP(I),TRAIN)
00242             MAXFC%R(I)=MAX(MAXFC%R(I),FSCEXP(I),TRAIN)
00243           ENDDO
00244         ELSEIF(YASMH) THEN
00245           DO I=1,NPOIN
00246             MINFC%R(I)=MIN(MINFC%R(I),FSCEXP(I))
00247             MAXFC%R(I)=MAX(MAXFC%R(I),FSCEXP(I))
00248           ENDDO
00249         ELSEIF(RAIN) THEN
00250           DO I=1,NPOIN
00251             MINFC%R(I)=MIN(MINFC%R(I),TRAIN)
00252             MAXFC%R(I)=MAX(MAXFC%R(I),TRAIN)
00253           ENDDO
00254         ENDIF
00255 !
00256 !       IN PARALLEL MODE: GLOBAL EXTREMA
00257 !
00258         IF(NCSIZE.GT.1) THEN
00259 !         ENSURING THAT THE WORK BIEF_OBJ HAVE THE RIGHT DISCRETISATION
00260           CALL CPSTVC(MASKPT,MINFC)
00261           CALL CPSTVC(MASKPT,MAXFC)
00262 !         GLOBAL EXTREMA
00263           CALL PARCOM(MAXFC,3,MESH)
00264           CALL PARCOM(MINFC,4,MESH)
00265         ENDIF
00266 !
00267       ENDIF
00268 !
00269 !---------------------------------------------------------------------
00270 !
00271 ! COMPUTES THE CRITERION FOR COURANT NUMBER
00272 !
00273       DO I = 1,NPOIN
00274         TAB1%R(I) = 0.D0
00275       ENDDO
00276 !
00277 !     USES HERE FXMAT ASSEMBLED IN PARALLEL FOR UPWINDING
00278 !
00279       IF(OPT.EQ.1) THEN
00280 !
00281 !       CLASSICAL CRITERION
00282 !
00283         DO I = 1,NSEG
00284           IF(FXMATPAR(I).GT.0.D0) THEN
00285 !           FXMAT(I) POSITIVE
00286 !           MAX(...,0.D0) FOR 1
00287             TAB1%R(GLOSEG(I,1)) = TAB1%R(GLOSEG(I,1)) + FXMAT(I)
00288 !           MIN(...,0.D0) FOR 2
00289             TAB1%R(GLOSEG(I,2)) =
00290      &      TAB1%R(GLOSEG(I,2)) - COEMIN * FXMAT(I)
00291           ELSE
00292 !           - FXMAT(I) POSITIVE
00293             TAB1%R(GLOSEG(I,2)) = TAB1%R(GLOSEG(I,2)) - FXMAT(I)
00294 !           MIN(...,0.D0) FOR 1
00295             TAB1%R(GLOSEG(I,1)) =
00296      &      TAB1%R(GLOSEG(I,1)) + COEMIN * FXMAT(I)
00297           ENDIF
00298         ENDDO
00299 !
00300       ELSEIF(OPT.EQ.2.OR.OPT.EQ.3) THEN
00301 !
00302 !       NEW CRITERION, COMPUTES MIN(FI_ij,0)*(Ci-Cj)
00303 !
00304         DO I=1,NSEG
00305           IF(FXMATPAR(I).LT.0.D0) THEN
00306             DIFF=FC(GLOSEG(I,1))-FC(GLOSEG(I,2))
00307             TAB1%R(GLOSEG(I,1))=TAB1%R(GLOSEG(I,1))-FXMAT(I)*DIFF
00308           ELSEIF(FXMATPAR(I).GT.0.D0) THEN
00309             DIFF=FC(GLOSEG(I,2))-FC(GLOSEG(I,1))
00310             TAB1%R(GLOSEG(I,2))=TAB1%R(GLOSEG(I,2))+FXMAT(I)*DIFF
00311           ENDIF
00312         ENDDO
00313 !
00314         DO I=1,NPTFR
00315           IF(LIMTRA(I).EQ.KDIR) THEN
00316             N=NBOR(I)
00317             DIFF=FC(N)-FBOR(I)
00318 !           FXBOR IS HERE IN GLOBAL NUMBERING
00319             TAB1%R(N)=TAB1%R(N)-MIN(FXBOR(N),0.D0)*DIFF
00320           ENDIF
00321         ENDDO
00322 !
00323         IF(YASMH) THEN
00324           DO I = 1,NPOIN
00325             DIFF=FC(I)-FSCEXP(I)
00326             TAB1%R(I)=TAB1%R(I)+MAX(SMH(I),0.D0)*DIFF
00327           ENDDO
00328         ENDIF
00329 !
00330         IF(RAIN) THEN
00331           DO I = 1,NPOIN
00332             DIFF=FC(I)-TRAIN
00333             TAB1%R(I)=TAB1%R(I)+MAX(PLUIE(I),0.D0)*DIFF
00334           ENDDO
00335         ENDIF
00336 !
00337       ENDIF
00338 !
00339       IF(NCSIZE.GT.1) CALL PARCOM(TAB1,2,MESH)
00340 !
00341 !     MASKS TAB1
00342 !
00343       IF(MSK) THEN
00344         CALL OS('X=XY    ',X=TAB1,Y=MASKPT)
00345       ENDIF
00346 !
00347 !     STABILITY (AND MONOTONICITY) CRITERION
00348 !
00349       DTMAX = DT
00350 !
00351       IF(OPT.EQ.1) THEN
00352 !
00353 !       SEE RELEASE NOTES 5.7, CRITERION AT THE END OF 4.4 PAGE 33
00354 !       BUT HERE THE FINAL H IS NOT H(N+1) BUT A FUNCTION OF DTMAX ITSELF
00355 !       H FINAL = HSTART + DTMAX/DT *(H-HSTART)
00356 !
00357         IF(YASMH.AND.RAIN) THEN
00358           DO I = 1,NPOIN
00359             DENOM=TAB1%R(I)+MAX(FXBOR(I),0.D0)
00360      &                     -MIN(SMH(I)+PLUIE(I),0.D0)
00361      &                     +COESOU*(MIN(FXBOR(I),0.D0)
00362      &                             -MAX(SMH(I)+PLUIE(I),0.D0))
00363             IF(DENOM.GT.1.D-20) THEN
00364               DTMAX = MIN(DTMAX,SECU*MAS(I)*HSTART(I)/DENOM)
00365             ENDIF
00366           ENDDO
00367         ELSEIF(YASMH) THEN
00368           DO I = 1,NPOIN
00369             DENOM=TAB1%R(I)+        MAX(FXBOR(I),0.D0)-MIN(SMH(I),0.D0)
00370      &                     +COESOU*(MIN(FXBOR(I),0.D0)-MAX(SMH(I),0.D0))
00371             IF(DENOM.GT.1.D-20) THEN
00372               DTMAX = MIN(DTMAX,SECU*MAS(I)*HSTART(I)/DENOM)
00373             ENDIF
00374           ENDDO
00375         ELSEIF(RAIN) THEN
00376           DO I = 1,NPOIN
00377             DENOM=TAB1%R(I)
00378      &                  +MAX(FXBOR(I),0.D0)-MIN(PLUIE(I),0.D0)
00379      &          +COESOU*(MIN(FXBOR(I),0.D0)-MAX(PLUIE(I),0.D0))
00380             IF(DENOM.GT.1.D-20) THEN
00381               DTMAX = MIN(DTMAX,SECU*MAS(I)*HSTART(I)/DENOM)
00382             ENDIF
00383           ENDDO
00384         ELSE
00385           DO I = 1,NPOIN
00386             DENOM=TAB1%R(I)+       MAX(FXBOR(I),0.D0)
00387      &                     +COESOU*MIN(FXBOR(I),0.D0)
00388             IF(DENOM.GT.1.D-20) THEN
00389               DTMAX = MIN(DTMAX,SECU*MAS(I)*HSTART(I)/DENOM)
00390             ENDIF
00391           ENDDO
00392         ENDIF
00393 !
00394       ELSEIF(OPT.EQ.2.OR.OPT.EQ.3) THEN
00395 !
00396 !       NEW CRITERION, SEE RELEASE NOTES 6.3
00397 !
00398         DO I= 1,NPOIN
00399           DI=TAB1%R(I)/MAS(I)
00400           IF(DI.GT.1.D-12) THEN
00401             A=(FC(I)-MINFC%R(I))/DI
00402             B=DT+A*(HSTART(I)-H(I))
00403             IF(B.GT.0.D0) THEN
00404               DTMAX = MIN(DTMAX,A*HSTART(I)*DT/B)
00405             ENDIF
00406           ELSEIF(DI.LT.-1.D-12) THEN
00407             A=(FC(I)-MAXFC%R(I))/DI
00408             B=DT+A*(HSTART(I)-H(I))
00409             IF(B.GT.0.D0) THEN
00410               DTMAX = MIN(DTMAX,A*HSTART(I)*DT/B)
00411             ENDIF
00412           ENDIF
00413         ENDDO
00414 !
00415       ENDIF
00416 !
00417 !-----------------------------------------------------------------------
00418 !
00419       RETURN
00420       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0