cvtrvf.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE cvtrvf
!                   *****************
!
     &(f,fn,fscexp,h,hn,hprop,uconv,vconv,dm1,zconv,solsys,
     & sm,smh,yasmh,smi,yasmi,fbor,masktr,mesh,
     & aggloh,dt,entet,msk,maskel,s,massou,optsou,
     & limtra,kdir,kddl,nptfr,flbor,yaflbor,volu2d,v2dpar,unsv2d,iopt,
     & flbortra,maskpt,rain,pluie,train,optadv,tb,free,am2,tb2,
     & nco_dist,nsp_dist,yaflulim,flulim,yaflulimebe,flulimebe,
     & slvtra)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    DISTRIBUTIVE EXPLICIT OR IMPLICIT ADVECTOR.
!
!warning  AFBOR AND BFBOR MUST BE 0 FOR THE BOUNDARY ELEMENTS
!+            WITH NO FRICTION
!
!history  CHI-TUAN PHAM  (LNHE)
!+        09/10/09
!+        V6P0
!+   First version.
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET (LNHE)
!+        24/02/2012
!+        V6P2
!+   Rain and evaporation added (after initiative by O. Boutron, from
!+   Tour du Valat, and O. Bertrand, Artelia group)
!
!history  SARA PAVAN & J-M HERVOUET (LNHE)
!+        18/06/2013
!+        V6P3
!+   New call to CFLVF, for new monotonicity criterion.
!
!history  SARA PAVAN & J-M HERVOUET (EDF LAB, LNHE)
!+        05/05/2014
!+        V7P0
!+   New predictor-corrector PSI scheme (OPTADV=2). Security coefficient
!+   on maximum time step to avoid truncation errors that would give
!+   negative values.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        16/05/2014
!+        V7P0
!+   Boundary conditions LIMTRA redone (they may have been done with
!+   U.N in diffin.f, which is different from flbor here. Argument
!+   VOLU2D added.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        01/09/2015
!+        V7P1
!+   Now with the locally implicit predictor-corrector.
!+   Call of CFLVF changed, with option OPTCFL added in the arguments.
!+   For locally implicit schemes, fluxes may be reduced with array
!+   FLULIM. The solver configuration of tracers is copied but the
!+   choice of the solver, if not GMRES or direct, is set to Jacobi
!+   in TVF_IMP.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        30/06/2016
!+        V7P2
!+   Stability criterion for the LIPS scheme was too restrictive
!+   compared to theory. The real theoretical value is now taken.
!
!history  J-M HERVOUET (jubilado)
!+        27/09/2016
!+        V7P2
!+   Predictor and corrector of LIPS grouped in the same loop.
!+   Limitation of the predictor even for the first correction.
!+   Extrema FMIN and FMAX depend on the scheme, so their computation
!+   is slightly changed. Coefficients COEMIN and COESOU removed in the
!+   call to cflvf.
!
!history  J-M HERVOUET (jubilado)
!+        20/10/2016
!+        V7P2
!+   Reduction of time-step for predictor-corrector is no longer fixed
!+   to be 2, any value > 1 is possible, but here hardcoded at 2.D0.
!+   So no change except for further research (values less than 2 seem
!+   to give better results).
!
!history  J-M HERVOUET (jubilado)
!+        29/09/2017
!+        V7P3
!+   Theory of stability condition slightly changed for all predictor-
!+   corrector and LIPS schemes. Steady state flows unchanged.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AGGLOH         |-->| MASS-LUMPING IN CONTINUITY EQUATION
!| DM1            |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
!| DT             |-->| TIME STEP
!| ENTET          |-->| LOGICAL, IF YES INFORMATION IS GIVEN ON MASS
!|                |   | CONSERVATION.
!| F              |<--| F AT TIME T(N+1)
!| FBOR           |-->| DIRICHLET CONDITIONS ON F.
!| FLBOR          |-->| FLUXES AT BOUNDARIES
!| FLBORTRA       |<->| TRACER FLUXES AT BOUNDARIES
!| FLULIM         |-->| A LIMITATION OF FLUXES IF YAFLULIM=.TRUE.
!| FN             |-->| F AT TIME T(N)
!| FSCEXP         |-->| EXPLICIT PART OF THE SOURCE TERM
!|                |   | EQUAL TO ZERO EVERYWHERE BUT ON SOURCES
!|                |   | WHERE THERE IS FSCE - (1-TETAT) FN
!|                |   | SEE DIFSOU
!| HPROP          |-->| PROPAGATION DEPTH (DONE IN CVDFTR).
!| IOPT           |-->| OPTIONS FOR COMPUTATION (NUMBER BETWEEN 0 AND 13)
!|                |   | THE TENS (IOPT2, I.E. 0 OR 1):
!|                |   | 0: UCONV OBEYS THE CONTINUITY EQUATION
!|                |   | 1: UCONV DOES NOT OBEY THE CONTINUITY EQUATION
!|                |   | THE UNITS (IOPT1, I.E. 0 TO 3): VARIANT FOR FLUXES
!|                |   | 0: CONSTANT PER ELEMENT = 0
!|                |   | 1: CHI-TUAN PHAM'S CONSTANT
!|                |   | 2: N SCHEME
!|                |   | 3: PSI SCHEME
!| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
!| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
!| LIMTRA         |-->| BOUNDARY CONDITIONS ON BOOUNDARY POINTS
!| MASKEL         |-->| MASKING OF ELEMENTS
!|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
!| MASKPT         |-->| MASKING PER POINT.
!| MASSOU         |-->| MASS OF TRACER ADDED BY SOURCE TERM
!|                |   | SEE DIFSOU
!| MESH           |-->| MESH STRUCTURE
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS.
!| NCO_DIST       |-->| NUMBER OF CORRECTIONS OF DISTRIBUTIVE SCHEMES
!| NSP_DIST       |-->| NUMBER OF SUB-STEPS OF DISTRIBUTIVE SCHEMES
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| OPDADV         |-->| SCHEME OPTION FOR THE ADVECTION OF TRACERS
!|                |   | WITH N OR PSI SCHEME:
!|                |   |  1: EXPLICIT
!|                |   |  2: PREDICTOR-CORRECTOR 1ST ORDER IN TIME
!|                |   |  3: PREDICTOR-CORRECTOR 2ND ORDER IN TIME
!|                |   |  4: LOCALLY IMPLICIT PREDICTOR-CORRECTOR
!| OPTSOU         |-->| TYPE OF SOURCES
!|                |   | 1: NORMAL
!|                |   | 2: DIRAC
!| S              |-->| VOID STRUCTURE
!| SLVTRA         |-->| SOLVER CONFIGURATION FOR TRACER DIFFUSION.
!| SM             |-->| SOURCE TERMS.
!| SMH            |-->| SOURCE TERM IN CONTINUITY EQUATION
!| SMI            |-->| IMPLICIT SOURCE TERM
!| SOLSYS         |-->| 1 OR 2. IF 2 ADVECTION FIELD IS UCONV + DM1*GRAD(ZCONV)
!| TB             |<->| BLOCK OF WORK STRUCTURES
!| TB2            |<->| SECOND BLOCK OF AT LEAST 7 WORK STRUCTURES
!| TRAIN          |-->| VALUE OF TRACER IN THE RAIN
!| UCONV,VCONV    |-->| ADVECTION VELOCITY FIELD
!| UNSV2D         |-->| INVERSE OF INTEGRALS OF TEST FUNCTIONS
!| VOLU2D         |-->| INTEGRAL OF TEST FUNCTIONS, NOT ASSEMBLED IN PARALLEL
!| V2DPAR         |-->| INTEGRAL OF TEST FUNCTIONS, ASSEMBLED IN PARALLEL
!| YAFLBOR        |-->| IF YES FLBOR IS GIVEN
!| YAFLULIM       |-->| IF YES FLULIM IS GIVEN
!| YASMH          |-->| IF YES, SMH MUST BE TAKEN INTO ACCOUNT
!| YASMI          |-->| IF YES, SMI MUST BE TAKEN INTO ACCOUNT
!| ZCONV          |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_cvtrvf => cvtrvf
      USE interface_parallel
      USE declarations_special
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: OPTSOU,KDIR,NPTFR,SOLSYS
      INTEGER, INTENT(IN)             :: KDDL,IOPT,OPTADV,FREE
      INTEGER, INTENT(IN)             :: NCO_DIST,NSP_DIST
      INTEGER, INTENT(INOUT)          :: LIMTRA(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: DT,AGGLOH,TRAIN
      DOUBLE PRECISION, INTENT(IN)    :: FLULIM(*)
      DOUBLE PRECISION, INTENT(INOUT) :: MASSOU
      LOGICAL, INTENT(IN)             :: YASMH,YAFLBOR
      LOGICAL, INTENT(IN)             :: MSK,ENTET,YASMI,RAIN
      LOGICAL, INTENT(IN)             :: YAFLULIM,YAFLULIMEBE
      TYPE(bief_obj), INTENT(IN)      :: MASKEL,DM1,ZCONV,MASKPT
      TYPE(bief_obj), INTENT(IN), TARGET :: H,HN
      TYPE(bief_obj), INTENT(IN)      :: VOLU2D,V2DPAR,UNSV2D,HPROP
      TYPE(bief_obj), INTENT(INOUT)   :: F,SM,AM2
      TYPE(bief_obj), INTENT(IN)      :: UCONV,VCONV,FN,SMI,SMH
      TYPE(bief_obj), INTENT(INOUT)   :: FBOR,TB,TB2
      TYPE(bief_obj), INTENT(INOUT)   :: FLBORTRA
      TYPE(bief_obj), INTENT(IN)      :: FSCEXP,S,MASKTR
      TYPE(bief_obj), INTENT(IN)      :: PLUIE
      TYPE(bief_mesh)                 :: MESH
      DOUBLE PRECISION, INTENT(IN)    :: FLULIMEBE(*)
      TYPE(bief_obj), INTENT(IN)      :: FLBOR
      TYPE(slvcfg), INTENT(IN)        :: SLVTRA
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELMF,I,IOPT1,IOPT2,N,IELEM,I1,I2,I3,ICOR,NIT,OPTCFL
!
!-----------------------------------------------------------------------
!
      DOUBLE PRECISION DT_REMAIN,DDT,TDT,SECU,TETAFCOR
      DOUBLE PRECISION ADMASS,LOCALMIN,LOCALMAX,TWOTHIRDS
!
      CHARACTER(LEN=16) FORMUL
!
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: FXMAT,FXMATPAR
!
      DOUBLE PRECISION K,SURK,KM1,KM1SURK
      LOGICAL MASS_BAL,PREDICOR,LIPS,WITHABS
!
      TYPE(slvcfg)::SLVPSI
!
!     HARDCODED PARAMETER !!!!!!!!!!
!
!     MAXIMUM NUMBER OF ITERATIONS
      INTEGER, PARAMETER :: NITMAX = 200
!
!-----------------------------------------------------------------------
!
      TYPE(bief_obj), POINTER :: T1,T2,FLBOUND,T4,T6,FXBORPAR,T8,HNT,HT
      TYPE(bief_obj), POINTER :: HNP1MT,TETAF_VAR,FMIN,FMAX
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: DFDT
!
!-----------------------------------------------------------------------
!
!     HARDCODED PARAMETER FOR STABILITY OF PREDICTOR CORRECTOR
!     AND NUMBERS DEPENDING OF K. K MUST BE > 1.D0
!
      k=2.d0
!
      surk=1.d0/k
      km1=k-1.d0
      km1surk=km1/k
      twothirds=2.d0*km1/(2.d0*k-1.d0)
!
!-----------------------------------------------------------------------
!
!     SOLVER OPTIONS FOR IMPLICIT SCHEMES : SLVTRA TAKEN BUT EITHER
!                                           GMRES OR DIRECT ASKED
!                                           AND HARDCODED HERE
!                                           KRYLOV DIMENSION SET TO 3
!
!     SEE ALSO TVF_IMP WITH A JACOBI METHOD
!
      slvpsi%SLV   =slvtra%SLV
      slvpsi%NITMAX=slvtra%NITMAX
      slvpsi%PRECON=slvtra%PRECON
      slvpsi%KRYLOV=slvtra%KRYLOV
      slvpsi%EPS   =slvtra%EPS
      slvpsi%ZERO  =slvtra%ZERO
!
!     CHANGING THIS WILL TRIGGER CHANGING THE SIZE OF TB2 IN CALLING
!     PROGRAMMES
      slvpsi%KRYLOV=3
!
!-----------------------------------------------------------------------
!
!     BEWARE : OPTIMISATION OF MEMORY USE...
!
      hnp1mt    =>tb%ADR(free)%P
      t1        =>tb%ADR(free+1)%P
      t2        =>tb%ADR(free+2)%P
      flbound   =>tb%ADR(free+3)%P
      t4        =>tb%ADR(free+4)%P
      t6        =>tb%ADR(free+5)%P
      fxborpar  =>tb%ADR(free+6)%P
      t8        =>tb%ADR(free+7)%P
      hnt       =>tb%ADR(free+8)%P
      ht        =>tb%ADR(free+9)%P
      tetaf_var =>tb%ADR(free+10)%P
!
!     BEWARE : FMIN AND FMAX MIXED IN MEMORY WITH T6%R AND T8%R
!              LOOK HOW DFDT AND T8 ARE ALWAYS RECOMPUTED AFTER FMIN AND FMAX
!              ARE USED  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     BEWARE : DFDT=T6%R !!!!!!!!!!!
      dfdt      => tb%ADR(free+5)%P%R
      fmin      => t6
      fmax      => t8
!
!-----------------------------------------------------------------------
!
      fxmat=>mesh%MSEG%X%R(1:mesh%NSEG)
!     IN PARALLEL MODE, ASSEMBLED AND NON ASSEMBLED VERSIONS ARE DIFFERENT
      IF(ncsize.GT.1) THEN
        fxmatpar=>mesh%MSEG%X%R(mesh%NSEG+1:2*mesh%NSEG)
      ELSE
        fxmatpar=>mesh%MSEG%X%R(1:mesh%NSEG)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     CHECKING OPTIONS
!
      IF(optadv.LT.1.OR.optadv.GT.4) THEN
        WRITE(lu,*) 'CVTRVF: OPTION OF ADVECTION SCHEME'
        WRITE(lu,*) '        UNKNOWN: ',optadv
        CALL plante(1)
        stop
      ENDIF
!
      IF(am2%STO.NE.3.AND.optadv.EQ.4) THEN
        WRITE(lu,*) 'CVTRVF: EDGE-BASED STORAGE REQUESTED WITH'
        WRITE(lu,*) '        IMPLICIT N OR PSI SCHEME'
        CALL plante(1)
        stop
      ENDIF
!
!     EXTRACTS OTHER OPTIONS
!
      iopt2=iopt/10
      iopt1=iopt-10*iopt2
!
      IF(iopt1.LT.1.OR.iopt1.GT.3) THEN
        WRITE(lu,*) 'CVTRVF : OPTION IOPT1=',iopt1,' UNKNOWN'
        CALL plante(1)
        stop
      ENDIF
!
      IF(iopt2.LT.0.OR.iopt2.GT.1) THEN
        WRITE(lu,*) 'CVTRVF : OPTION IOPT2=',iopt2,' UNKNOWN'
        CALL plante(1)
        stop
      ENDIF
!
!     SIMPLIFYING TESTS: OPTIONS WITH AN EXPLICIT PREDICTOR-CORRECTOR
!
      IF((iopt1.EQ.2.OR.iopt1.EQ.3).AND.
     &  (optadv.EQ.2.OR.optadv.EQ.3)) THEN
        predicor=.true.
      ELSE
        predicor=.false.
      ENDIF
!
!     SIMPLIFYING TESTS: OPTION WITH LOCALLY IMPLICIT PSI (OR N) SCHEME
!
      IF((iopt1.EQ.2.OR.iopt1.EQ.3).AND.optadv.EQ.4) THEN
        lips=.true.
      ELSE
        lips=.false.
      ENDIF
!
!     SELECTING THE IMPLICITATION COEFFICIENT
!
      IF(lips) THEN
!       ADAPTIVE IMPLICIT, N OR PSI
      ELSEIF((iopt1.EQ.2.OR.iopt1.EQ.3).AND.optadv.EQ.3) THEN
!       EXPLICIT PREDICTOR, PSEUDO-IMPLICIT CORRECTOR
!       FROM 0 TO 1
        tetafcor=0.5d0
      ELSE
!       EXPLICIT N OR PSI, EXPLICIT PREDICTOR-CORRECTOR PSI
        tetafcor=0.d0
      ENDIF
!
!-----------------------------------------------------------------------
!
!     IELMF = F%ELM
!     FORCED TO LINEAR
      ielmf=11
!
!     TAKES MASS-LUMPING INTO ACCOUNT IN THE CONTINUITY EQUATION
!
      IF(abs(1.d0-aggloh).GT.1.d-8) THEN
        CALL vector(ht ,'=','MASVEC          ',ielmf,
     &              1.d0-aggloh,h ,s,s,s,s,s,mesh,msk,maskel)
        CALL vector(hnt,'=','MASVEC          ',ielmf,
     &              1.d0-aggloh,hn,s,s,s,s,s,mesh,msk,maskel)
        IF(ncsize.GT.1) THEN
          CALL parcom(ht ,2,mesh)
          CALL parcom(hnt,2,mesh)
        ENDIF
        CALL os('X=YZ    ',x=ht ,y=ht ,z=unsv2d)
        CALL os('X=YZ    ',x=hnt,y=hnt,z=unsv2d)
        CALL os('X=X+CY  ',x=ht ,y=h  ,c=aggloh)
        CALL os('X=X+CY  ',x=hnt,y=hn ,c=aggloh)
      ELSE
!       CALL OS('X=Y     ',X=HT ,Y=H )
!       CALL OS('X=Y     ',X=HNT,Y=HN)
        ht=>h
        hnt=>hn
      ENDIF
!
!     IF NO FLBOR IS GIVEN, IT IS COMPUTED
!
      IF(yaflbor) THEN
        CALL os('X=Y     ',x=flbound,y=flbor)
      ELSE
!       MASK=5 FOR NON NEUMANN BOUNDARIES IN DIFFIN
        CALL vector(flbound,'=','FLUBDF          ',1,1.d0,
     &              hprop,hprop,hprop,
     &              uconv,vconv,vconv,mesh,.true.,masktr%ADR(5)%P)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     CORRECTION OF THE BOUNDARY CONDITIONS
!
!-----------------------------------------------------------------------
!
!     A SIMILAR CORRECTION IS DONE IN DIFFIN, BUT IT MAY BE INCOMPATIBLE
!     AS U*N MAY NOT BE OF THE SAME SIGN AS FLBOR HERE (RARE BUT ALREADY
!     SEEN).
!
      CALL cpstvc(hn,fxborpar)
!     JUST IN CASE INTERNAL POINTS HAVE NON INITIALISED VALUES
      CALL os('X=0     ',x=fxborpar)
      DO i=1,mesh%NPTFR
        n=mesh%NBOR%I(i)
        fxborpar%R(n)=flbound%R(i)
      ENDDO
      IF(ncsize.GT.1) CALL parcom(fxborpar,2,mesh)
      DO i=1,mesh%NPTFR
        n=mesh%NBOR%I(i)
        IF(limtra(i).EQ.kdir.AND.fxborpar%R(n).GT.0.d0) THEN
          limtra(i)=kddl
        ELSEIF(limtra(i).EQ.kddl.AND.fxborpar%R(n).LT.0.d0) THEN
          limtra(i)=kdir
!         WHEN VELOCITIES RE-ENTER THROUGH AN EXIT, WE ARBITRARILY CHOOSE
!         THAT THE PRESCRIBED VALUE IS THE LAST KNOWN
          fbor%R(i)=fn%R(n)
        ENDIF
      ENDDO
!
!     MASKS FLBOR IF(MSK)
!
      IF(msk) CALL os('X=XY    ',x=fxborpar,y=maskpt)
!
!-----------------------------------------------------------------------
!
!     INITIALISES THE ADDED MASS AND THE TRACER FLUX AT THE BOUNDARY
!
      admass=0.d0
      DO i=1,mesh%NPTFR
        flbortra%R(i)=0.d0
      ENDDO
!
!     COMPUTES THE FLUXES PHIIJ = FXMAT
!
      formul='HUGRADP         '
      IF(solsys.EQ.2) formul(8:8)='2'
      CALL vector(t2,'=',formul,ielmf,-1.d0,
     &            hprop,dm1,zconv,uconv,vconv,vconv,mesh,msk,maskel)
!                 T2 AS HUGRADP IS NOT USED AS AN ASSEMBLED VECTOR
!                 BUT TO GET THE NON ASSEMBLED FORM MESH%W
      nit=0
      dt_remain=dt
      tdt=0.d0
      CALL cpstvc(h,hnp1mt)
      CALL cpstvc(h,t4)
      CALL cpstvc(f,t8)
!
!     T4 WILL BE F PROGRESSIVELY UPDATED
!     HNP1MT WILL BE THE DEPTH AT THE END OF THE PREVIOUS SUB-TIMESTEP
!     HENCE THE DEPTH At THE BEGINNING OF THE NEXT SUB TIME-STEP
!     (INITIALISED HERE TO CALL CFLVF)
!
      DO i=1,hn%DIM1
        t4%R(i)=fn%R(i)
        hnp1mt%R(i)=hnt%R(i)
      ENDDO
!
!     T1 WILL BE THE DEPTH ACCORDING TO THE CONTINUITY EQUATION
!
      IF(iopt2.EQ.1) THEN
        DO i=1,hn%DIM1
          t1%R(i)=hnt%R(i)
        ENDDO
      ENDIF
!
100   CONTINUE
      nit=nit+1
!
!---------------------------------------
! VARIOUS OPTIONS TO COMPUTE THE FLUXES
!---------------------------------------
!
      IF(nit.EQ.1.OR.(iopt1.EQ.3.AND..NOT.lips)) THEN
        CALL flux_ef_vf(fxmat,mesh%W%R,mesh%NSEG,mesh%NELEM,
     &                  mesh%NELMAX,
     &                  mesh%ELTSEG%I,mesh%ORISEG%I,
!                                          N SCHEME FORCED (HERE FOR CFL ONLY)
     &                  mesh%IKLE%I,.true.,2,
!                       OPTIONAL LIMITATION OF FLUXES FOR LIPS
!                       (OTHER SCHEMES DO NOT WORK WITH TIDAL FLATS)
     &                  yaflulim=yaflulim.AND.lips,flulim=flulim,
     &                  yaflulimebe=yaflulimebe.AND.lips,
     &                  flulimebe=flulimebe)
!    &                  MESH%IKLE%I,.TRUE.,2    ,FN=T4)
!                                          IOPT1 HERE FORCED TO N SCHEME
!                                          NO NEED OF FN
!
!       CANCELS FLUXES TO AND FROM MASKED POINTS
!
        IF(msk) THEN
          CALL flux_mask(fxmat,mesh%NSEG,
     &                   mesh%GLOSEG%I,mesh%GLOSEG%DIM1,maskpt%R)
        ENDIF
!       ASSEMBLES THE FLUXES AT INTERFACES IN PARALLEL MODE, THIS
!       IS FOR UPWINDING (STORED IN SECOND DIMENSION OF MESH%MSEG)
        IF(ncsize.GT.1) THEN
          CALL ov('X=Y     ', x=fxmatpar, y=fxmat, dim1=mesh%NSEG)
          CALL parcom2_seg(fxmatpar,fxmatpar,fxmatpar,
     &                     mesh%NSEG,1,2,1,mesh,1,11)
        ENDIF
      ENDIF
!
!--------------------------------------------
! DETERMINES THE LARGEST ADMISSIBLE TIMESTEP
!--------------------------------------------
!
!     COMPUTES THE MAXIMUM TIMESTEP ENSURING MONOTONICITY
!     ACCORDING TO THEORY OF N SCHEME, BUT POSSIBLY
!     DIVIDED BY 2 FOR PREDICTOR-CORRECTORS
!
      IF(predicor.AND.nco_dist.GT.0) THEN
!       SECU=SURK
        secu=1.d0
        withabs=.true.
      ELSEIF(lips) THEN
        secu=1.d0
        withabs=.true.
      ELSE
        secu=0.99d0
        withabs=.false.
      ENDIF
!
!     HARDCODED OPTION FOR THE CFL (AND OPTION 2 ONLY FOR THE N SCHEME)
!
      optcfl=1
!     OPTCFL=2
!     IF(OPTADV.NE.1.OR.IOPT1.NE.2) OPTCFL=1
!
      CALL cflvf(ddt,hnp1mt%R,fxmat,fxmatpar,
!                                   FLBOR%R(NPOIN)
     &           v2dpar%R,dt_remain,fxborpar%R   ,smh%R,
     &           yasmh,t8,mesh%NSEG,mesh%NPOIN,mesh%NPTFR,
     &           mesh%GLOSEG%I,mesh%GLOSEG%DIM1,mesh,msk,maskpt,
     &           rain,pluie%R,t4%R,mesh%NELEM,mesh%IKLE%I,
     &           limtra,kdir,kddl,fbor%R,fscexp%R,train,mesh%NBOR%I,
     &           t2,t6,secu,optcfl,withabs)
!
!     NOW RECOMPUTING THE PSI FLUXES (THE N FLUXES HAVE BEEN
!                                     USED FOR THE STABILITY CRITERION).
!
!     WITH THE SEMI-IMPLICIT SCHEME, THE PSI FLUXES CANNOT BE USED
!     AT THIS LEVEL AND WILL NOT BE COMPUTED.
!
      IF(iopt1.EQ.3.AND..NOT.lips) THEN
        CALL flux_ef_vf(fxmat,mesh%W%R,mesh%NSEG,mesh%NELEM,
     &                  mesh%NELMAX,
     &                  mesh%ELTSEG%I,mesh%ORISEG%I,
     &                  mesh%IKLE%I,.true.,iopt1,t4)
!       CANCELS FLUXES TO AND FROM MASKED POINTS
        IF(msk) THEN
          CALL flux_mask(fxmat,mesh%NSEG,
     &                   mesh%GLOSEG%I,mesh%GLOSEG%DIM1,maskpt%R)
        ENDIF
!       ASSEMBLES THE FLUXES AT INTERFACES IN PARALLEL MODE, THIS
!       IS FOR UPWINDING (STORED IN SECOND DIMENSION OF MESH%MSEG)
        IF(ncsize.GT.1) THEN
          CALL ov('X=Y     ', x=fxmatpar, y=fxmat, dim1=mesh%NSEG)
          CALL parcom2_seg(fxmatpar,fxmatpar,fxmatpar,
     &                     mesh%NSEG,1,2,1,mesh,1,11)
        ENDIF
      ENDIF
!
!
      IF(ncsize.GT.1) ddt=p_min(ddt)
      ddt=min(ddt,dt_remain)
!
!     CASE OF ADAPTIVE IMPLICIT N OR PSI SCHEME
!
      IF(lips) THEN
        DO i=1,hn%DIM1
! V7P2
!         TETAF_VAR%R(I)=MAX(0.D0,1.D0-SURK*NSP_DIST*T2%R(I)/DT)
! V7P3
          tetaf_var%R(i)=max(0.d0,1.d0-nsp_dist*t2%R(i)/dt)
        ENDDO
        ddt=dt/nsp_dist/0.999999d0
      ENDIF
!
      ddt=min(ddt,dt_remain)
!
!     T2 WILL TAKE THE SUCCESSIVE VALUES OF H
!     AT THE BEGINNING OF THE SUB-TIMESTEP
!     WARNING: T2 ALSO USED WITH IOPT2=1, BUT SO FAR PREDICTOR-CORRECTOR
!              NOT USED WITH THIS OPTION
!
      IF(predicor) THEN
        DO i=1,hn%DIM1
          t2%R(i)=hnt%R(i)+tdt*(ht%R(i)-hnt%R(i))/dt
        ENDDO
      ENDIF
!
      tdt=tdt+ddt
!
!     HNP1MT WILL TAKE THE SUCCESSIVE VALUES OF H
!     AT THE END OF THE SUB-TIMESTEP (EXPLICIT) OR IN BETWEEN (IMPLICIT)
!
      IF(lips) THEN
        DO i=1,hn%DIM1
          hnp1mt%R(i)=hnt%R(i)+
     &               (tdt-tetaf_var%R(i)*ddt)*(ht%R(i)-hnt%R(i))/dt
        ENDDO
      ELSE
        DO i=1,hn%DIM1
          hnp1mt%R(i)=hnt%R(i)+tdt*(ht%R(i)-hnt%R(i))/dt
        ENDDO
      ENDIF
!
!     IN TVF FACTOR HT/HLIN MAY TRIGGER DIVERGENCE FOR DRY POINTS
!
      IF(msk) THEN
        DO i=1,hn%DIM1
          IF(maskpt%R(i).LT.0.5d0) hnp1mt%R(i)=ht%R(i)
        ENDDO
      ENDIF
!
!------------------------------------
!  FINAL RESOLUTION OR PREDICTOR STEP
!------------------------------------
!
      IF(.NOT.lips) THEN
!
        CALL tracvf(f,fscexp,fxmat,fxmatpar,
     &              volu2d,unsv2d,
     &              ddt,flbound,fbor,smh,yasmh,t1,t2,t4,hnp1mt,
     &              fxborpar,t8,
     &              mesh,limtra,kdir,kddl,optsou,iopt2,flbortra,msk,
     &              dt,rain,pluie,train,admass,
!                   CASES WHERE THE MASS BALANCE MUST BE PREPARED
     &              optadv.EQ.1.OR.nco_dist.EQ.0)
!
      ELSE
!
!       LOOP WITH PREDICTOR AND ALL CORRECTIONS
!
        DO icor=0,nco_dist
!
!         UNLIKE WITH EXPLICIT SCHEMES, THE PREDICTOR IS NOT
!         GUARANTEED WITHOUT STABILITY PROBLEM IF NOT LIMITED
          IF(icor.GT.0) THEN
!           LIMITING THE PREDICTOR
            DO i=1,hn%DIM1
              fmin%R(i)=t4%R(i)
              fmax%R(i)=t4%R(i)
            ENDDO
            DO i=1,nptfr
              IF(limtra(i).EQ.kdir) THEN
                n=mesh%NBOR%I(i)
                fmin%R(n)=min(fmin%R(n),fbor%R(i))
                fmax%R(n)=max(fmax%R(n),fbor%R(i))
              ENDIF
            ENDDO
            DO ielem=1,mesh%NELEM
              i1=mesh%IKLE%I(ielem)
              i2=mesh%IKLE%I(ielem+  mesh%NELMAX)
              i3=mesh%IKLE%I(ielem+2*mesh%NELMAX)
              localmin=min( f%R(i1), f%R(i2), f%R(i3),
     &                     t4%R(i1),t4%R(i2),t4%R(i3))
              localmax=max( f%R(i1), f%R(i2), f%R(i3),
     &                     t4%R(i1),t4%R(i2),t4%R(i3))
              fmin%R(i1)=min(fmin%R(i1),localmin)
              fmax%R(i1)=max(fmax%R(i1),localmax)
              fmin%R(i2)=min(fmin%R(i2),localmin)
              fmax%R(i2)=max(fmax%R(i2),localmax)
              fmin%R(i3)=min(fmin%R(i3),localmin)
              fmax%R(i3)=max(fmax%R(i3),localmax)
            ENDDO
            IF(ncsize.GT.1) THEN
              CALL parcom(fmin,4,mesh)
              CALL parcom(fmax,3,mesh)
            ENDIF
            DO i=1,hn%DIM1
              f%R(i)=min(f%R(i),t4%R(i)+km1surk*(fmax%R(i)-t4%R(i)))
              f%R(i)=max(f%R(i),t4%R(i)+km1surk*(fmin%R(i)-t4%R(i)))
            ENDDO
          ENDIF
!
!         DERIVATIVE IN TIME
!
          IF(icor.EQ.0) THEN
!           FOR GETTING THE EXPLICIT PSI CONTRIBUTION
!           WITHOUT THE DERIVATIVE IN TIME
            DO i=1,hn%DIM1
              dfdt(i)=0.d0
            ENDDO
          ELSE
!           FOR GETTING THE EXPLICIT PSI CONTRIBUTION
!           WITH THE DERIVATIVE IN TIME
            DO i=1,hn%DIM1
!             HERE DFDT=H*DFDT WITH SEMI IMPLICIT H
              dfdt(i)=hnp1mt%R(i)*(f%R(i)-t4%R(i))/ddt
            ENDDO
          ENDIF
!
!         EXPLICIT CONTRIBUTION
!
          IF(icor.EQ.0) THEN
!           FN GIVEN FOR FSTAR
            CALL flux_ef_vf_3(mesh%W%R,mesh%NELEM,mesh%NELMAX,
     &                        mesh%ELTSEG%I,mesh%ORISEG%I,
     &                        fxmatpar,mesh%NSEG,
!                                                          FN
     &                        mesh%IKLE%I,mesh%NPOIN,t4,
!    &                        FI_I               HDFDT
     &                        t8%R,mesh%SURFAC%R,dfdt,tetaf_var%R,
!                                             FN FOR FSTAR
     &                        yaflulim,flulim,yaflulimebe,
     &                        flulimebe)
          ELSE
!           FSTAR REALLY GIVEN
            CALL flux_ef_vf_3(mesh%W%R,mesh%NELEM,mesh%NELMAX,
     &                        mesh%ELTSEG%I,mesh%ORISEG%I,
     &                        fxmatpar,mesh%NSEG,
!                                                          FN
     &                        mesh%IKLE%I,mesh%NPOIN,t4,
!    &                        FI_I               HDFDT
     &                        t8%R,mesh%SURFAC%R,dfdt,tetaf_var%R,
!                                             FSTAR
     &                        yaflulim,flulim,yaflulimebe,flulimebe)
          ENDIF
!         NO, WILL GO INTO SM THAT IS NOT ASSEMBLED
!         IF(NCSIZE.GT.1) CALL PARCOM(T8,2,MESH)
!
!         BUILDING AND SOLVING THE LINEAR SYSTEM (NORMAL OR PREDICTOR)
!
          CALL tvf_imp(f%R,t4%R,fxmat,
     &                 fxmatpar,ddt,
     &                 flbound%R,hnp1mt%R,
     &                 fbor%R,smh%R,yasmh,fscexp%R,
     &                 mesh%NSEG,mesh%NPOIN,mesh%NPTFR,
     &                 mesh%GLOSEG%I,mesh%GLOSEG%DIM1,
     &                 mesh%NBOR%I,limtra,kdir,kddl,
     &                 optsou,iopt2,flbortra%R,ddt/dt,mesh,f,
     &                 rain,pluie%R,train,tetaf_var%R,
     &                 entet,volu2d%R,
     &                 t6,t8%R,am2,tb2,slvpsi,
!                      PREDICTOR CORRECTOR
     &                 icor.EQ.0,icor.NE.0,icor,nco_dist,admass)
!
        ENDDO
!
      ENDIF
!
!-------------------------------------
!  CORRECTOR STEP FOR N AND PSI SCHEME
!-------------------------------------
!
      IF(predicor.AND.nco_dist.GT.0) THEN
!
        DO icor=1,nco_dist
!
!       CASES WITH A LIMITATION OF THE FIRST CORRECTOR
!
        IF(optadv.EQ.3.OR.icor.GT.1.OR.k.LT.2.d0) THEN
!         COMPUTING THE MINIMUM AND MAXIMUM
          DO i=1,hn%DIM1
            fmin%R(i)=t4%R(i)
            fmax%R(i)=t4%R(i)
          ENDDO
          DO i=1,nptfr
            IF(limtra(i).EQ.kdir) THEN
              n=mesh%NBOR%I(i)
              fmin%R(n)=min(fmin%R(n),fbor%R(i))
              fmax%R(n)=max(fmax%R(n),fbor%R(i))
            ENDIF
          ENDDO
          IF(optadv.EQ.2) THEN
!           ALL FSTAR(J) NOT INCLUDED
            DO ielem=1,mesh%NELEM
              i1=mesh%IKLE%I(ielem)
              i2=mesh%IKLE%I(ielem+  mesh%NELMAX)
              i3=mesh%IKLE%I(ielem+2*mesh%NELMAX)
              localmin=min(t4%R(i1),t4%R(i2),t4%R(i3))
              localmax=max(t4%R(i1),t4%R(i2),t4%R(i3))
              fmin%R(i1)=min(f%R(i1),fmin%R(i1),localmin)
              fmax%R(i1)=max(f%R(i1),fmax%R(i1),localmax)
              fmin%R(i2)=min(f%R(i2),fmin%R(i2),localmin)
              fmax%R(i2)=max(f%R(i2),fmax%R(i2),localmax)
              fmin%R(i3)=min(f%R(i3),fmin%R(i3),localmin)
              fmax%R(i3)=max(f%R(i3),fmax%R(i3),localmax)
            ENDDO
          ELSEIF(optadv.EQ.3) THEN
!           ALL FSTAR(J) INCLUDED
            DO ielem=1,mesh%NELEM
              i1=mesh%IKLE%I(ielem)
              i2=mesh%IKLE%I(ielem+  mesh%NELMAX)
              i3=mesh%IKLE%I(ielem+2*mesh%NELMAX)
              localmin=min( f%R(i1), f%R(i2), f%R(i3),
     &                     t4%R(i1),t4%R(i2),t4%R(i3))
              localmax=max( f%R(i1), f%R(i2), f%R(i3),
     &                     t4%R(i1),t4%R(i2),t4%R(i3))
              fmin%R(i1)=min(fmin%R(i1),localmin)
              fmax%R(i1)=max(fmax%R(i1),localmax)
              fmin%R(i2)=min(fmin%R(i2),localmin)
              fmax%R(i2)=max(fmax%R(i2),localmax)
              fmin%R(i3)=min(fmin%R(i3),localmin)
              fmax%R(i3)=max(fmax%R(i3),localmax)
            ENDDO
          ENDIF
          IF(ncsize.GT.1) THEN
            CALL parcom(fmin,4,mesh)
            CALL parcom(fmax,3,mesh)
          ENDIF
        ENDIF
!
!       LIMITATIONS OF THE PREDICTOR
!
        IF(optadv.EQ.3) THEN
          DO i=1,hn%DIM1
            f%R(i)=min(f%R(i),t4%R(i)+km1*(t4%R(i)-fmin%R(i)))
            f%R(i)=max(f%R(i),t4%R(i)+km1*(t4%R(i)-fmax%R(i)))
          ENDDO
          IF(icor.GT.1.OR.k.LT.2.d0) THEN
            DO i=1,hn%DIM1
              f%R(i)=min(f%R(i),t4%R(i)+twothirds*(fmax%R(i)-t4%R(i)))
              f%R(i)=max(f%R(i),t4%R(i)+twothirds*(fmin%R(i)-t4%R(i)))
            ENDDO
          ENDIF
        ELSEIF(optadv.EQ.2.AND.(icor.GT.1.OR.k.LT.2.d0)) THEN
          DO i=1,hn%DIM1
            f%R(i)=min(f%R(i),t4%R(i)+km1surk*(fmax%R(i)-t4%R(i)))
            f%R(i)=max(f%R(i),t4%R(i)+km1surk*(fmin%R(i)-t4%R(i)))
          ENDDO
        ENDIF
!
        DO i=1,hn%DIM1
          dfdt(i)=(f%R(i)-t4%R(i))/ddt
        ENDDO
!
        CALL flux_ef_vf_2(mesh%W%R,mesh%NELEM,mesh%NELMAX,
!                                                      FN
     &                    mesh%IKLE%I,iopt1,mesh%NPOIN,t4,
!    &                    FI_I,FSTAR, HN   H
     &                    t8%R,f%R   ,t2%R,hnp1mt%R,mesh%SURFAC%R,
     &                    tetafcor,dfdt)
        IF(ncsize.GT.1) CALL parcom(t8,2,mesh)
!
        IF(icor.EQ.nco_dist) THEN
          mass_bal=.true.
        ELSE
          mass_bal=.false.
        ENDIF
        DO i=1,hn%DIM1
!         WILL BE FSTAR IN TVF_2
          t6%R(i)=f%R(i)
        ENDDO
        CALL tvf_2(f%R,t6%R,t4%R,volu2d%R,unsv2d%R,ddt,
     &             flbound%R,fxborpar%R,fbor%R,smh%R,yasmh,fscexp%R,
     &             mesh%NPOIN,mesh%NPTFR,mesh%NBOR%I,limtra,kdir,kddl,
     &             optsou,hnp1mt%R,iopt2,flbortra%R,ddt/dt,rain,
     &             pluie%R,train,t8%R,tetafcor,mass_bal,admass)
!                                FI_I
!
        ENDDO
!
      ENDIF  ! IF(PREDICOR.AND.NCO_DIST.GT.1)
!
!-----------------
! END CORRECTOR STEP
!-----------------
!
!     PREPARING FOR FURTHER ITERATIONS
!
      IF(.NOT.(lips.AND.nsp_dist.EQ.1)) THEN
!
        DO i=1,hn%DIM1
!         T4 IS F(N+1)
          t4%R(i)=f%R(i)
        ENDDO
        IF(iopt2.EQ.1) THEN
          DO i=1,hn%DIM1
            t1%R(i)=t2%R(i)
          ENDDO
        ENDIF
!
      ENDIF
!
      dt_remain=dt_remain-ddt
!
      IF(dt_remain.NE.0.d0.AND.nit.LT.nitmax) GO TO 100
!
      IF(nit.GE.nitmax.AND.optadv.NE.4) THEN
        WRITE(lu,901) nit
901     FORMAT(1x,'CVTRVF: ',1i6,' SUB-ITERATIONS REQUIRED FOR THE'
     &   ,/,1x,   '        DISTRIBUTIVE SCHEME. DECREASE THE TIME-STEP')
        CALL plante(1)
        stop
      ELSEIF(entet) THEN
        WRITE(lu,903) nit
903     FORMAT(1x,'CVTRVF (BIEF): ',1i6,' SUB-ITERATIONS')
      ENDIF
!
!-----------------------------------------------------------------------
!
!     EXPLICIT SOURCE TERM
!
      DO i = 1,mesh%NPOIN
        f%R(i) = f%R(i)+dt*sm%R(i)
      ENDDO
!
!     IMPLICIT SOURCE TERM
!
      IF(yasmi) THEN
        DO i = 1,mesh%NPOIN
          f%R(i) = f%R(i)/(1.d0-dt*smi%R(i)/max(h%R(i),1.d-15))
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
!     UPDATING MASSOU
!
      IF(ncsize.GT.1) THEN
        massou=massou+p_sum(admass)
      ELSE
        massou=massou+admass
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END