cvdftr.F

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE cvdftr
!                   *****************
!
     &(f,ftild,fn,fscexp,dift,iconvf,conv,
     & h,hn,hprop,uconv,vconv,dm1,zconv,solsys,
     & visc,visc_s,sm,smh,yasmh,smi,yasmi,am1,am2,
     & zf,fbor,afbor,bfbor,limtra,masktr,mesh,tb,
     & t1,t2,t4,t10,te1,te2,te3,kdir,kddl,dt,entet,
     & tetat,agglot,infogt,bilan,optadv,
     & isousi,opdtra,optban,msk,maskel,maskpt,mbor,
     & s,massou,optsou,slvtra,flbor,volu2d,v2dpar,unsv2d,optvf,flbortra,
     & flulim,yaflulim,flulimebe,yaflulimebe,dirflu,rain,pluie,train,
     & given_flux,flux_given,maxadv,tb2,nco_dist,nsp_dist)
!
!***********************************************************************
! BIEF   V7P3
!***********************************************************************
!
!brief    DIFFUSION, ADVECTION AND SOURCE TERMS FOR A TRACER.
!code
!+  THE EQUATION SOLVED IS :
!+
!+
!+          N+1                                            TILD
!+         F           1                                  F   + DT*SM
!+      ---------  -  ---  DIV ( H VISC * GRAD ( F   )) = ____________
!+         DT          H                                      DT
!+
!+                                                      N+1
!+                                  + SOURCES  + SMI * F
!+                                                     ___
!+                                                     H
!+
!+     WITH :    N+1  TILD   N
!+              F   ,F     ,F  =    DIFFUSED FUNCTION
!+              VISC           =    TURBULENT VISCOSITY
!+              SM             =    SECOND MEMBER (SOURCE TERMS)
!+              TETAT          =    IMPLICITATION COEFFICIENT
!+              DT             =    TIME STEP
!+                                         N+1              N
!+              F              =    TETAT F  + (1-TETAT) * F
!+              SMI            =    IMPLICIT SOURCE TERM
!+
!+
!+                    TILD       N
!+     DISTINGUISHES F     FROM F   IN CASE A FRACTIONAL STEPS METHOD
!+
!+     HAD BEEN PREVIOUSLY PERFORMED (ADVECTION FOR EXAMPLE) GIVING
!+
!+      TILD       N
!+     F     FROM F
!+
!+-----------------------------------------------------------------------
!+
!+      BOUNDARY CONDITIONS :
!+
!+      ==>   NEUMANN CONDITION
!+
!+      VISC DF/DN = AFBOR . F  +  BFBOR
!+
!+
!+      ==>   DIRICHLET CONDITION
!+
!+            TREATED BY MODIFICATION OF THE EQUATIONS IN THE
!+            SUBROUTINE DIRICH
!
!note     JMH : W IS NOT USED.
!
!warning  MATDIF DOES NOT GIVE THE DIFFUSION MATRIX. IT MISSES THE
!+            BOUNDARY TERMS AND THERE IS A MINUS SIGN WHICH IS TAKEN
!+            INTO ACCOUNT HERE
!warning  AFBOR AND BFBOR MUST BE 0 FOR THE BOUNDARY ELEMENTS
!+            WITH NO FRICTION
!warning  BEWARE DISCRETISATION OF VISC
!
!history  JMH
!+        27/02/2009
!+
!+   CALLS CVTFVF_POS, OPTION 14
!
!history  J-M HERVOUET (LNHE)     ; C MOULIN (LNH)
!+        29/12/2009
!+        V6P0
!+
!
!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  J-M HERVOUET (LNHE)
!+        17/07/2012
!+        V6P2
!+   Mass balance programmed in case of implicit source term.
!+   Arguments GIVEN_FLUX and FLUX_GIVEN added.
!
!history  J-M HERVOUET (LNHE)
!+        09/04/2013
!+        V6P3
!+   DIMGLO=MESH%GLOSEG%DIM1 used in call to CVTRVF_POS. Strangely
!+   avoids an "array temporary created" with Intel compiler.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        01/04/2014
!+        V7P0
!+   Now written to enable different numbering of boundary points and
!+   boundary segments.
!+   OPTSUP changed into OPTADV for a more general use.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        16/05/2014
!+        V7P0
!+   Argument VOLU2D added for calling CVTRVF.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        30/09/2014
!+        V7P0
!+   Checking mass in case of implicit source terms SMI depends on the
!+   advection scheme, a section was added for distributive schemes.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        20/08/2015
!+        V7P1
!+   Adaptation to new subroutine cvtrvf with locally implicit
!+   predictor-corrector distributive schemes.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        20/11/2015
!+        V7P1
!+   Implicitation on diffusion forced to 1. with parameter TETATD. It
!+   should be a keyword but there is no reason so far to do an explicit
!+   diffusion.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        03/06/2016
!+        V7P2
!+   Adapting to new scheme ERIA, ADV_PSI_TF = 15. Computation of HPROP
!+   removed, HPROP must now be given and is not rebuilt with TETA,
!+   which was a mistake. Adding OPTADV in cvtrvf_pos.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        06/09/2016
!+        V7P2
!+   Splitting cvtrvf_pos into cvtrvf_nerd and cvtrvf_eria + cleaning.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AFBOR,BFBOR    |-->| COEFFICIENTS OF NEUMANN CONDITION
!|                |   | VISC*DF/DN = AFBOR*F + BFBOR
!|                |   | GIVEN FOR EVERY BOUNDARY POINT
!| AGGLOT         |-->| MASS-LUMPING COEFFICIENT FOR T.
!| AM1            |<->| MATRIX.
!| AM2            |<->| MATRIX.
!| BILAN          |-->| LOGICAL, IF YES A BALANCE OF MASS EXCHANGES HAS
!|                |   | TO BE DONE.
!| CONV           |-->| IF YES ADVECTION OF F
!| DIFT           |-->| IF YES, DIFFUSION IS DONE
!| DIRFLU         |-->| 1: PRIORITY TO DIRICHLET VALUES (CALL DIRICH DONE)
!|                |   | 2: PRIORITY TO FLUXES (CALL DIRICH NOT DONE)
!| 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
!| FLUX_GIVEN     |-->| IF GIVEN_FLUX=YES, THE FLUX IS GIVEN IN
!|                |   | GIVEN_FLUX
!| 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
!| FTILD          |-->| F AFTER ADVECTION
!| GIVEN_FLUX     |-->| IF GIVEN_FLUX=YES, THE FLUX IS GIVEN IN
!|                |   | GIVEN_FLUX AND WILL NOT BE COMPUTED HERE
!| HPROP          |-->| WORK ARRAY
!| ICONVF         |-->| OPTION FOR ADVECTION TERMS
!|                |   | ICONVF = 1 : CHARACTERISTICS.
!|                |   | ICONVF = 2 : S.U.P.G.
!|                |   | ICONVF = 3 : CONSERVATIVE FINITE VOLUMES
!|                |   | ICONVF = 4 : IDEM
!|                |   | ICONVF = 6 : NON CONSERVATIVE PSI SCHEME.
!|                |   | ICONVF = 7 : NON CONSERVATIVE N SCHEME.
!|                |   | ICONVF =13 : EDGE BY EDGE FORM OF 3
!|                |   | ICONVF =14 : IDEM
!| INFOGT         |-->| LOGICAL, IF YES INFORMATION ON SOLVER WILL BE
!|                |   | PRINTED.
!| ISOUSI         |-->| SUB-ITERATION NUMBER
!| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
!| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
!| KENT           |-->| CONVENTION FOR INFLOW POINT
!| LIMTRA         |-->| BOUNDARY CONDITIONS ON BOUNDARY POINTS
!| MASKEL         |-->| MASKING OF ELEMENTS
!|                |   | =1. : NORMAL   =0. : MASKED ELEMENT
!| MASKPT         |-->| MASKING PER POINT.
!| MASSOU         |-->| MASS OF TRACER ADDED BY SOURCE TERM
!|                |   | SEE DIFSOU
!| MAXADV         |-->| MAXIMUM NUMBER OF ITERATIONS OF ADVECTION SCHEMES
!| MBOR           |-->| BOUNDARY MATRIX
!| 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
!| OPDTRA         |-->| OPTION FOR THE DIFFUSION OF TRACERS
!| OPTADV         |-->| ADVECTION SCHEME OPTION, THE MEANING DEPENDS ON
!|                |   | THE SCHEME. IF SCHEME IS SUPG:
!|                |   | 0: NO SUPG UPWIND
!|                |   | 1: CLASSIC SUPG
!|                |   | 2: MODIFIED SUPG
!|                |   | IF SCHEME IS PSI:
!|                |   | 1: EXPLICIT
!|                |   | 2: PREDICTOR-CORRECTOR
!|                |   | 3: SECOND ORDER PREDICTOR-CORRECTOR
!|                |   | 4: LOCALLY IMPLICIT
!| OPTBAN         |-->| OPTION FOR THE TREATMENT OF TIDAL FLATS
!|                |   | 1:NORMAL   2:WITH MASKING
!| OPTSOU         |-->| TYPE OF SOURCES
!|                |   | 1: NORMAL
!|                |   | 2: DIRAC
!| OPTVF          |-->| OPTIONS FOR FINITE VOLUMES (SEE CVTRVF)
!| PLUIE          |-->| RAIN OR EVAPORATION, IN M/S
!| RAIN           |-->| IF YES: RAIN OR EVAPORATION
!| S              |-->| VOID STRUCTURE
!| SLVTRA         |-->| SOLVER CONFIGURATION (SLVCFG) STRUCTURE
!|                |   | CONTAINING DATA FOR CALLING SOLVE
!| 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)
!| T1             |<->| WORK BIEF_OBJ STRUCTURE
!| T10            |<->| WORK BIEF_OBJ STRUCTURE
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| T4             |<->| WORK BIEF_OBJ STRUCTURE
!| TB             |<->| BLOCK OF WORK BIEF_OBJ STRUCTURES (CONTAINS T1,...)
!| TB2            |<->| SECOND BLOCK OF WORK BIEF_OBJ STRUCTURES
!|                |   | FOR IMPLICIT SCHEMES ONLY
!| TE1,TE2,TE3    |<->| WORK BIEF_OBJ STRUCTURE FOR ELEMENTS
!| TETAT          |-->| IMPLICITATION COEFFICIENT OF ADVECTION
!| TRAIN          |-->| VALUE OF TRACER IN THE RAIN
!| UCONV,VCONV    |-->| ADVECTION VELOCITY FIELD
!| UNSV2D         |-->| =1/V2DPAR
!| VOLU2D         |-->| INTEGRAL OF TEST FUNCTIONS (NOT ASSEMBLED IN PARALLEL)
!| V2DPAR         |-->| INTEGRAL OF TEST FUNCTIONS (ASSEMBLED IN PARALLEL)
!| VISC           |-->| VISCOSITY COEFFICIENTS ALONG X,Y AND Z .
!|                |   | IF P0 : PER ELEMENT
!|                |   | IF P1 : PER POINT
!| VISC_S         |<->| WORK ARRAY FOR SAVING VISC
!| YASMH          |-->| IF YES SMH TAKEN INTO ACCOUNT
!| YASMI          |-->| IF YES SMI TAKEN INTO ACCOUNT
!| ZCONV          |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
!| ZF             |-->| BOTTOM ELEVATION.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_cvdftr => cvdftr
      USE interface_parallel
      USE declarations_telemac, ONLY : adv_car,adv_sup,adv_nsc,adv_psi,
     &    adv_lpo,adv_nsc_tf,adv_psi_tf,adv_lpo_tf,adv_void
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)           :: ICONVF,ISOUSI,OPTADV,OPDTRA
      INTEGER, INTENT(IN)           :: OPTBAN,OPTSOU,KDIR,SOLSYS
      INTEGER, INTENT(IN)           :: KDDL,OPTVF,DIRFLU,MAXADV
      INTEGER, INTENT(IN)           :: NCO_DIST,NSP_DIST
      DOUBLE PRECISION, INTENT(IN)  :: TETAT,AGGLOT,DT,TRAIN
      DOUBLE PRECISION, INTENT(INOUT)  :: MASSOU
      LOGICAL, INTENT(IN)           :: INFOGT,BILAN,CONV,YASMH,RAIN
      LOGICAL, INTENT(IN)           :: DIFT,MSK,ENTET,YASMI,YAFLULIM
      LOGICAL, INTENT(IN)           :: FLUX_GIVEN,YAFLULIMEBE
      TYPE(slvcfg), INTENT(INOUT)   :: SLVTRA
      TYPE(bief_obj), INTENT(IN)    :: MASKEL,MASKPT,H,HN,AFBOR,BFBOR
      TYPE(bief_obj), INTENT(IN), TARGET :: HPROP
      TYPE(bief_obj), INTENT(INOUT) :: F,SM,FLBORTRA
      TYPE(bief_obj), INTENT(IN)    :: UCONV,VCONV,ZF
      TYPE(bief_obj), INTENT(IN)    :: FTILD,FN,SMI,FLULIM,PLUIE
      TYPE(bief_obj), INTENT(IN)    :: FLULIMEBE
      TYPE(bief_obj), INTENT(INOUT) :: SMH,FBOR
      TYPE(bief_obj), INTENT(INOUT) :: TE1,TE2,TE3
      TYPE(bief_obj), INTENT(INOUT) :: T1,T2,T4
      TYPE(bief_obj), INTENT(INOUT), TARGET :: T10
      TYPE(bief_obj), INTENT(IN)    :: FSCEXP,DM1,ZCONV
      TYPE(bief_obj), INTENT(IN)    :: S,V2DPAR,UNSV2D,VOLU2D
      TYPE(bief_obj), INTENT(INOUT) :: FLBOR,LIMTRA
      TYPE(bief_obj), INTENT(INOUT) :: VISC_S,VISC
      TYPE(bief_obj), INTENT(INOUT) :: AM1,AM2,MBOR
      TYPE(bief_obj), INTENT(INOUT) :: TB,TB2
      TYPE(bief_obj), INTENT(IN)    :: MASKTR,GIVEN_FLUX
      TYPE(bief_mesh)               :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION C,TETATD
!
      INTEGER IELMF,IELMH,IELMS,MSKNEU,I,N,IOPT,DIMGLO
!
      LOGICAL MSQ,FV_SCHEME
!
      CHARACTER(LEN=16) FORMUL
!
      TYPE(bief_obj), POINTER :: HPR
!
!-----------------------------------------------------------------------
!
!     HARDCODED IMPLICITATION COEFFICIENT FOR DIFFUSION
!     WAITING FOR A KEYWORD... BUT NO REAL INTEREST TO HAVE AN
!     EXPLICIT FORM, EXCEPT TO AVOID A LINEAR SYSTEM...
!
      tetatd=1.d0
!
!-----------------------------------------------------------------------
!
      ielmf = f%ELM
      ielmh = h%ELM
      ielms = sm%ELM
      dimglo=mesh%GLOSEG%DIM1
!
!-----------------------------------------------------------------------
!
!     IS IT A FINITE VOLUME SCHEME FOR ADVECTION ?
!
      fv_scheme=.false.
      IF(  iconvf.EQ.adv_lpo.OR.iconvf.EQ.adv_lpo_tf.OR.
     &     iconvf.EQ.adv_nsc.OR.iconvf.EQ.adv_nsc_tf.OR.
     &     iconvf.EQ.adv_psi.OR.iconvf.EQ.adv_psi_tf     ) THEN
        fv_scheme=.true.
      ENDIF
!
!-----------------------------------------------------------------------
!
!     CASE WHERE H AND T DON'T HAVE THE SAME DISCRETISATION
!
      IF(ielmf.NE.ielms) THEN
        CALL chgdis(sm ,ielms,ielmf,mesh)
      ENDIF
      IF(ielmf.NE.ielmh.AND.yasmh) THEN
        CALL chgdis(smh ,ielmh,ielmf,mesh)
      ENDIF
!
!
!-----------------------------------------------------------------------
!
!     HPROP IS THE DEPTH IN THE DIVERGENCE TERM IN THE CONTINUITY
!     EQUATION. IT IS REPLACED HERE BY POINTER HPR, IN CASE IT HAS
!     TO BE EXTENDED.
!
      IF(ielmf.NE.ielmh) THEN
        CALL os( 'X=Y     ' , x=t10   , y=hprop )
        CALL chgdis(t10,ielmh,ielmf,mesh)
        hpr=>t10
      ELSE
        hpr=>hprop
      ENDIF
!
!-----------------------------------------------------------------------
!
!     INITIALISES THE VARIABLES
!
!     SOLUTION INITIALISED TO F AT TIME N
      IF(isousi.EQ.1) CALL os( 'X=Y     ' , x=f , y=fn )
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!  IF SUPG, BUILDS THE SEMI-IMPLICIT MATRIX + SUPG IN AM2
!
      IF(iconvf.EQ.adv_sup.AND.conv) THEN
!
!       TERM IN U.GRAD(T) CENTERED:
!
        CALL matrix(am2,'M=N     ','MATVGR          ',ielmf,ielmf,
     &              1.d0,s,s,s,uconv,vconv,s,
     &              mesh,msk,maskel)
!
!       ADDS SUPG CONTRIBUTION TO AM2
!
        IF(optadv.EQ.1) THEN
!         CLASSICAL SUPG
          CALL ksupg(te1,te2,1.d0,uconv,vconv,mesh)
          CALL matrix(am2,'M=M+N   ','MASUPG          ',ielmf,ielmf,
     &                1.d0,te1,te2,s,uconv,vconv,s,mesh,msk,maskel)
!
        ELSEIF(optadv.EQ.2) THEN
!         MODIFIED SUPG
          CALL matrix(am2,'M=M+N   ','MAUGUG          ',ielmf,ielmf,
     &                0.5d0*dt,s,s,s,uconv,vconv,s,mesh,msk,maskel)
        ELSEIF(optadv.NE.0) THEN
          WRITE(lu,*) 'CVDFTR: UNKNOWN OPTION FOR SUPG =',optadv
          CALL plante(1)
          stop
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!   COMPUTES AM1: MASS MATRIX MULTIPLIED BY 1/DT
!
      IF(dift.OR..NOT.fv_scheme.OR..NOT.conv.OR.bilan) THEN
!
      IF(optban.EQ.1.OR.optban.EQ.3) THEN
        CALL os('X=Y+Z   ',x=t2,y=zf,z=hn)
        CALL decvrt(te3,t2,zf,mesh)
!       MASS MATRIX LOCALLY LUMPED ON THE TIDAL FLATS
        formul='MSLUMP          '
!       WILL USE MSQ TO MASK THE DIFFUSION
        msq=.true.
        IF(msk) CALL os('X=XY    ',x=te3,y=maskel)
      ELSE
!       NORMAL MASS MATRIX
        formul='MATMAS          '
!       MASK FOR THE DIFFUSION = MASKEL
        IF(msk) CALL os('X=Y     ',x=te3,y=maskel)
        msq=msk
      ENDIF
      CALL matrix(am1,'M=N     ',formul,ielmf,ielmf,
     &            1.d0/dt,te3,s,s,s,s,s,mesh,msk,maskel)
!
!   POSSIBLE MASS-LUMPING
!
      IF(agglot.GT.0.001d0) THEN
        CALL lump(t1,am1,mesh,agglot)
        CALL om('M=CM    ', m=am1, c=1.d0-agglot, mesh=mesh)
        CALL om('M=M+D   ', m=am1, d=t1, mesh=mesh)
      ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!   COMPUTES THE SECOND MEMBERS
!
!     COMPUTES DT * SM IN T2
!
!     CVTRVF AND CVTRVF_POS WILL TREAT SM IN A DIFFERENT WAY
!
      IF(.NOT.fv_scheme) THEN
        CALL os( 'X=CY    ' , x=t2 , y=sm , c=dt )
      ENDIF
!
!=======================================================================
! TREATS THE VARIOUS TYPES OF ADVECTION:
!-----------------------------------------------------------------------
!
      IF(iconvf.EQ.adv_car.OR..NOT.conv) THEN
!
        CALL os( 'X=X+Y   ' , x=t2 , y=ftild )
        CALL matvec( 'X=AY    ',sm,am1,t2,c,mesh)
!
!-----------------------------------------------------------------------
!
      ELSEIF(iconvf.EQ.adv_sup.AND.conv) THEN
!
!       AM1 MADE NONSYMMETRICAL IF IT WAS NOT ALREADY
!
        IF(am1%TYPEXT.NE.'Q') THEN
          CALL om('M=X(M)  ', m=am1, mesh=mesh)
        ENDIF
!
!       SUPG CONTRIBUTION TO THE MASS MATRIX
!
        IF(optadv.EQ.1) THEN
!         CLASSICAL SUPG
!         TE1 AND TE2 ALREADY COMPUTED
          CALL matrix(am1,'M=M+TN    ','MATVGR          ',ielmf,ielmf,
     &                1.d0/dt,s,s,s,te1,te2,s,mesh,msk,maskel)
!
        ELSEIF(optadv.EQ.2) THEN
!         MODIFIED SUPG
          CALL matrix(am1,'M=M+TN    ','MATVGR          ',ielmf,ielmf,
     &                0.5d0,s,s,s,uconv,vconv,s,mesh,msk,maskel)
        ELSEIF(optadv.NE.0) THEN
          WRITE(lu,*) 'CVDFTR: UNKNOWN OPTION FOR SUPG =',optadv
          CALL plante(1)
          stop
        ENDIF
!
!       END OF THE SUPG CONTRIBUTION TO THE MASS MATRIX
!
        CALL os( 'X=X+Y   ' , x=t2 , y=fn )
        CALL matvec( 'X=AY    ',sm,am1,t2,c,mesh)
!
! EXPLICIT ADVECTION TERM:
!
        CALL matvec( 'X=X+CAY ',sm,am2,fn,tetat-1.d0,mesh)
!
! ADDS THE IMPLICIT ADVECTION PART IN AM2 TO AM1
!
        CALL om('M=M+CN  ', m=am1, n=am2, c=tetat, mesh=mesh)
!
!-----------------------------------------------------------------------
!
      ELSEIF( (iconvf.EQ.adv_lpo.OR.
     &         iconvf.EQ.adv_nsc.OR.
     &         iconvf.EQ.adv_psi    ).AND.conv ) THEN
!
! CONSERVATIVE EQUATION, DISTRIBUTIVE SCHEMES (LEO POSTMA, N AND PSI)
!                        LEO POSTMA AND N-SCHEME ARE THE SAME IN 2D
!
!       TO BE REMOVED WHEN ALL CALLS TO CVDFTR ARE CHECKED
!       OPTVF SHOULD BE 0 (VELOCITY FIELD OBEYS THE CONTINUITY EQUATION)
!       OR 10 (VELOCITY FIELD DOES NOT OBEY THE CONTINUITY EQUATION)
        iopt=10*(optvf/10)
!       OPTION TO DISTRIBUTE THE FLUXES (HERE 2 OR 3)
        IF(iconvf.EQ.adv_lpo) iopt=iopt+2
        IF(iconvf.EQ.adv_nsc) iopt=iopt+2
        IF(iconvf.EQ.adv_psi) iopt=iopt+3
!
        IF(tb%N.LT.22) THEN
          WRITE(lu,*) 'SIZE OF TB TOO SMALL IN CVDFTR'
          CALL plante(1)
          stop
        ENDIF
        CALL cvtrvf(f,fn,fscexp,h,hn,hprop,uconv,vconv,
     &              dm1,zconv,solsys,sm,smh,yasmh,smi,yasmi,
     &              fbor,masktr,mesh,agglot,dt,entet,
     &              msk,maskel,s,massou,optsou,
!                                                       YAFLBOR
     &              limtra%I,kdir,kddl,mesh%NPTFR,flbor,.true.,
     &              volu2d,v2dpar,unsv2d,iopt,flbortra,maskpt,
     &              rain,pluie,train,optadv,tb,12,am2,tb2,
     &              nco_dist,nsp_dist,yaflulim,flulim%R,
     &              yaflulimebe,flulimebe%R,slvtra)
!       IF EXITS AT THIS POINT, THE DIRICHLET ARE NOT DONE, ALSO WORKS
!       CAN THEN CHECK THE MASS CONSERVATION EXACTLY
        IF(.NOT.dift) RETURN
        CALL matvec( 'X=AY    ',sm,am1,f,c,mesh)
!
!-----------------------------------------------------------------------
!
      ELSEIF( (iconvf.EQ.adv_lpo_tf.OR.
     &         iconvf.EQ.adv_nsc_tf.OR.
     &         iconvf.EQ.adv_psi_tf    ).AND.conv ) THEN
!
!       SCHEMES NERD OR ERIA, FOR TIDAL FLATS
!       LEO POSTMA AND N-SCHEME ARE THE SAME IN 2D
!
!       OPTION TO DISTRIBUTE THE FLUXES (HERE 2 OR 3 AND 12 OR 13)
        iopt=10*(optvf/10)
        IF(tb%N.LT.20) THEN
          WRITE(lu,*) 'SIZE OF TB TOO SMALL IN CVDFTR'
          WRITE(lu,*) 'FOR CALLING CVTRVF_NERD AND CVTRVF_ERIA'
          CALL plante(1)
          stop
        ENDIF
        IF(iconvf.EQ.adv_lpo_tf.OR.
     &     iconvf.EQ.adv_nsc_tf) THEN
!
!         NERD
          iopt=iopt+2
          CALL cvtrvf_nerd(f,fn,fscexp,h,hn,hprop,uconv,vconv,
     &              dm1,zconv,solsys,sm,smh,yasmh,smi,yasmi,
     &              fbor,masktr,mesh,
     &              tb%ADR(13)%P,tb%ADR(14)%P,tb%ADR(15)%P,
     &              tb%ADR(16)%P,tb%ADR(17)%P,tb%ADR(18)%P,
     &              tb%ADR(19)%P,tb%ADR(20)%P,
     &              dt,entet,msk,maskel,optsou,
!                                                       YAFLBOR
     &              limtra%I,kdir,kddl,mesh%NPTFR,flbor,.true.,
     &              unsv2d,iopt,flbortra,
     &              mesh%GLOSEG%I(       1:  dimglo),
     &              mesh%GLOSEG%I(dimglo+1:2*dimglo),
     &              mesh%NBOR%I,flulim%R,yaflulim,rain,pluie,
     &              train,given_flux,flux_given,maxadv)
!
        ELSEIF(iconvf.EQ.adv_psi_tf) THEN
!
!         ERIA
          iopt=iopt+3
          CALL cvtrvf_eria(f,fn,fscexp,h,hn,hprop,uconv,vconv,
     &              dm1,zconv,solsys,sm,smh,yasmh,smi,yasmi,
     &              fbor,masktr,mesh,
     &              tb%ADR(13)%P,tb%ADR(14)%P,tb%ADR(15)%P,
     &              tb%ADR(16)%P,tb%ADR(17)%P,tb%ADR(18)%P,
     &              tb%ADR(19)%P,tb%ADR(20)%P,
     &              dt,entet,msk,maskel,optsou,
!                                                       YAFLBOR
     &              limtra%I,kdir,kddl,mesh%NPTFR,flbor,.true.,
     &              unsv2d,iopt,flbortra,
     &              mesh%NBOR%I,rain,pluie,
     &              train,maxadv,nco_dist,optadv)
!
        ENDIF
!

!       IF EXITS AT THIS POINT, THE DIRICHLET ARE NOT DONE, ALSO WORKS
!       CAN THEN CHECK THE MASS CONSERVATION EXACTLY
        IF(.NOT.dift) RETURN
        CALL matvec( 'X=AY    ',sm,am1,f,c,mesh)
!
!-----------------------------------------------------------------------
!
      ELSEIF( iconvf.EQ.adv_void.AND.conv ) THEN
!
        CALL os( 'X=Y     ' , x=f, y=fn)
        CALL os( 'X=X+Y   ' , x=t2 , y=fn )
        CALL matvec( 'X=AY    ',sm,am1,t2,c,mesh)
!
!-----------------------------------------------------------------------
!
      ELSE
        WRITE(lu,*) 'CVDFTR: UNKNOWN ADVECTION OPTION : ',iconvf
        CALL plante(1)
        stop
!
!-----------------------------------------------------------------------
!
      ENDIF
!
! END OF TREATMENT OF THE VARIOUS TYPES OF ADVECTION:
!=======================================================================
!
!                   COMPUTES THE MATRICES
!
!-----------------------------------------------------------------------
!   COMPUTES AM2 : - DIFFUSION MATRIX, AND BOUNDARY TERMS
!
      IF(dift) THEN
!
        IF(opdtra.EQ.2) THEN
!         SAVES THE DIFFUSION
          CALL os('X=Y     ',x=visc_s,y=visc)
!         MULTIPLIES THE DIFFUSION BY HPROP
          CALL ov_2('X=XY    ',visc%R,1,hpr%R,1,hpr%R,1,c,
     &                         visc%MAXDIM1,visc%DIM1)
          IF(visc%DIM2.EQ.3) THEN
          CALL ov_2('X=XY    ',visc%R,2,hpr%R,1,hpr%R,1,c,
     &                         visc%MAXDIM1,visc%DIM1)
          CALL ov_2('X=XY    ',visc%R,3,hpr%R,1,hpr%R,1,c,
     &                         visc%MAXDIM1,visc%DIM1)
          ENDIF
        ENDIF
!
!       COMPUTES THE DIFFUSION MATRIX (OPTION WITH MONOTONICITY)
!
        CALL matrix(am2,'M=N     ','MATDIF       MON',ielmf,ielmf,
     &              1.d0,s,s,s,visc,s,s,mesh,msq,te3)
!
        IF(opdtra.EQ.2) THEN
!         MULTIPLIES THE MATRIX BY 1/HPROP
          CALL os( 'X=1/Y   ',x=t4,y=hpr,
     &             iopt=2,infini=0.d0,zero=1.d-2)
          CALL om('M=X(M)  ', m=am2, mesh=mesh )
          CALL om('M=DM    ', m=am2, d=t4, mesh=mesh)
!         RETURNS THE DIFFUSION
          CALL os('X=Y     ',x=visc,y=visc_s)
        ENDIF
!
!   TAKES THE BOUNDARY TERMS INTO ACCOUNT IN THE DIFFUSION MATRIX
!
        mskneu=3
!
        IF(mesh%NELEB.GT.0) THEN
          CALL matrix(mbor,'M=N     ','FMATMA          ',
     &                ielbor(ielmf,1),ielbor(ielmf,1),
     &                -1.d0,afbor,s,s,s,s,s,
     &                mesh,.true.,masktr%ADR(mskneu)%P)
          CALL om('M=M+N   ', m=am2, n=mbor, mesh=mesh)
        ENDIF
!
!       DIFFUSION TERMS
!
        IF(am1%TYPEXT.NE.'Q'.AND.am2%TYPEXT.EQ.'Q') THEN
          CALL om('M=X(M)  ', m=am1, mesh=mesh)
        ENDIF
!
        IF(tetatd.NE.1.d0) THEN
!         EXPLICIT DIFFUSION TERM
          CALL matvec( 'X=AY    ',t1,am2,fn,c,mesh)
          CALL os( 'X=X+CY  ' , x=sm , y=t1 , c=tetatd-1.d0 )
!         IMPLICIT DIFFUSION TERM ( AM1 + TETAT * AM2 )
          CALL om('M=M+CN  ', m=am1, n=am2, c=tetatd, mesh=mesh)
        ELSE
!         IMPLICIT DIFFUSION TERM ( AM1 + AM2 )
          CALL om('M=M+N   ', m=am1, n=am2, mesh=mesh)
        ENDIF
!
!       BOUNDARY STRESS TERMS
!
        CALL vector(t2,'=','MASVEC          ',ielbor(ielmf,1),
     &              1.d0,bfbor,s,s,s,s,s,mesh,
     &              .true.,masktr%ADR(mskneu)%P)
        CALL osdb( 'X=X+Y   ' , sm , t2 , t2 , c , mesh )
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!     IMPLICIT AND EXPLICIT TERM DUE TO THE SOURCES:
!
      IF(yasmh.AND..NOT.(fv_scheme.AND.conv)) THEN
!
        IF(optsou.EQ.1) THEN
          CALL vector(t2,'=','MASVEC          ',ielmf,
     &                1.d0,smh,s,s,s,s,s,mesh,msk,maskel)
          CALL os('X=Y/Z   ',x=t1,y=t2,z=hpr,
     &              iopt=2,infini=0.d0,zero=1.d-3)
!         IMPLICIT PART OF THE POINT SOURCE TERM
!         - TETAT T 1/HPROP SUM ( SCE PSI D(OMEGA)
          CALL os('X=CX    ', x=t1 , c=tetat )
          CALL om('M=M+D   ', m=am1, d=t1, mesh=mesh)
!         PREPARES THE EXPLICIT PART
          CALL os('X=YZ    ', x=t1 , y=smh , z=fscexp)
          CALL vector(t2,'=','MASVEC          ',ielmf,
     &                1.d0,t1,s,s,s,s,s,mesh,msk,maskel)
          CALL os('X=Y/Z   ',x=t1,y=t2,z=hpr,
     &             iopt=2,infini=0.d0,zero=1.d-3)
          CALL os('X=X+Y   ', x=sm , y=t1)
        ELSEIF(optsou.EQ.2) THEN
          CALL os('X=Y/Z   ',x=t1,y=smh,z=hpr,
     &              iopt=2,infini=0.d0,zero=1.d-3)
!         EXPLICIT PART OF THE POINT SOURCE TERM
!         1/HPROP (FSCE-(1-TETAT)FN) SMH
          CALL os('X=X+YZ  ', x=sm , y=t1 , z=fscexp)
!         IMPLICIT PART OF THE POINT SOURCE TERM
!         - TETAT T 1/HPROP SUM ( SCE PSI D(OMEGA)
          CALL os('X=CX    ', x=t1, c=tetat)
          CALL om('M=M+D   ', m=am1, d=t1, mesh=mesh)
        ENDIF
!
      ENDIF
!
!     IMPLICIT AND EXPLICIT TERM DUE TO RAIN AND EVAPORATION
!
      IF(rain.AND..NOT.(fv_scheme.AND.conv)) THEN
!
        CALL vector(t2,'=','MASVEC          ',ielmf,
     &              1.d0,pluie,s,s,s,s,s,mesh,msk,maskel)
        CALL os('X=Y/Z   ' ,x=t1,y=t2,z=hpr,
     &          iopt=2,infini=0.d0,zero=1.d-3)
!       IMPLICIT PART OF THE POINT SOURCE TERM
!       - TETAT T 1/HPROP SUM ( SCE PSI D(OMEGA)
        CALL os('X=CX    ',x=t1,c=tetat)
        CALL om('M=M+D   ', m=am1, d=t1, mesh=mesh)
!       PREPARES THE EXPLICIT PART
        CALL os('X=CYZ   ',x=t1,y=pluie,z=fn,c=tetat-1.d0)
!       VALUE OF TRACER IN RAIN CONSIDERED (ONLY EXPLICIT HENCE LIKE TETAT=0)
        IF(abs(train).GT.1.d-15) THEN
          DO i=1,fn%DIM1
!                           ONLY REAL RAIN CONSIDERED, NOT EVAPORATION
            t1%R(i)=t1%R(i)+max(pluie%R(i),0.d0)*train
          ENDDO
        ENDIF
        CALL vector(t2,'=','MASVEC          ',ielmf,
     &              1.d0,t1,s,s,s,s,s,mesh,msk,maskel)
        CALL os('X=Y/Z   ',x=t1,y=t2,z=hpr,
     &          iopt=2,infini=0.d0,zero=1.d-3)
        CALL os('X=X+Y   ',x=sm,y=t1)
!
      ENDIF
!
!   IMPLICIT TERM IF THERE IS ONE :
!
!   THE TREATMENT BELOW ENSURES THAT IF THE EXPLICIT SOURCE TERM
!   IS IN THE FORM  K*FN/H AND SMI EQUALS -K THEN THE TWO TERMS
!   WILL BE BALANCED (CASE OF EROSION AND DEPOSITION)
!
!                  FV_SCHEME : IMPLICIT SOURCE TERM HAS BEEN TREATED
!                  AND WILL NOT BE DONE TWICE
      IF(yasmi.AND..NOT.(fv_scheme.AND.conv)) THEN
        CALL matrix(am2,'M=N     ','MATMAS          ',ielmf,ielmf,
     &              -1.d0,s,s,s,s,s,s,mesh,msk,maskel)
!       POSSIBLE MASS-LUMPING
        IF(agglot.GT.0.001d0) THEN
          CALL lump(t1,am2,mesh,agglot)
          CALL om('M=CM    ', m=am2, c=1.d0-agglot, mesh=mesh)
          CALL om('M=M+D   ', m=am2, d=t1, mesh=mesh)
        ENDIF
!       COMPUTES SMI/H (DOES NOT CHECK IF H SIZE IS SUFFICIENT!!)
        IF(optban.GT.0) THEN
!         DIVIDES BY H WITH HARD-CODED CLIPPING AT 1.D-4
          CALL cpstvc(smi,t4)
          DO i=1,smi%DIM1
            t4%R(i)=smi%R(i)/max(h%R(i),1.d-4)
          ENDDO
        ELSE
!         DIVIDES WITHOUT CHECKING
          CALL os( 'X=Y/Z   ',x=t4,y=smi,z=h)
        ENDIF
        CALL om('M=X(M)  ', m=am2, mesh=mesh)
        CALL om('M=MD    ', m=am2, d=t4, mesh=mesh)
!       ADDS TO MATRIX AM1
        IF(am1%TYPEXT.NE.'Q') THEN
          CALL om('M=X(M)  ', m=am1, mesh=mesh)
        ENDIF
        CALL om('M=M+N   ', m=am1, n=am2, mesh=mesh)
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(iconvf.EQ.adv_car.AND..NOT.dift) THEN
        CALL os('X=Y     ', x=f, y=ftild)
      ENDIF
!
!-----------------------------------------------------------------------
!
!   BOUNDARY CONDITIONS (POINTS OF THE TYPE DIRICHLET)
!
!     DIRICHLET BOUNDARY CONDITIONS ARE TREATED IN A WEAK FORM BY
!     FINITE VOLUME SCHEMES, TO ENSURE THE CORRECT FLUX
!     THEN CALL DIRICH HERE WOULD SPOIL MASS CONSERVATION. HOWEVER
!     FORCING THE CALL IS ALLOWED IN THIS CASE WITH PARAMETER DIRFLU
!
!
      IF(.NOT.fv_scheme.OR.dirflu.EQ.1) THEN
        CALL dirich(f, am1, sm,fbor,limtra%I,tb,mesh,kdir,msk,maskpt)
      ENDIF
!
!-----------------------------------------------------------------------
!
!   SOLVES THE LINEAR SYSTEM:
!
#if defined COMPAD
      CALL ad_solve(f,am1,sm,tb,slvtra,infogt,mesh,am2)
!!                                               NOT USED HERE
#else
      CALL solve(f,am1,sm,tb,slvtra,infogt,mesh,am2)
!!                                               NOT USED HERE
#endif
!
!-----------------------------------------------------------------------
!
      IF(yasmi.AND.bilan) THEN
        IF(yasmi.AND..NOT.(fv_scheme.AND.conv)) THEN
!         WE ASSUME HERE THAT AM2 HAS NOT BEEN DESTROYED IN
!         SOLVE (THIS WOULD BE THE CASE WITH SOME PRECONDITIONING
!         BUT IN THIS CASE AM2 COULD NOT BE USED BEFORE)
          CALL matvec( 'X=AY    ',t2,am2,f,c,mesh)
          IF(ncsize.GT.1) THEN
            CALL parcom(t2,2,mesh)
            massou = massou - dt*p_dots(t2,hprop,mesh)
          ELSE
            massou = massou - dt*dots(t2,hprop)
          ENDIF
        ELSEIF(yasmi.AND.(fv_scheme.AND.conv)) THEN
!         THIS IS CONSISTENT WITH WHAT IS DONE
!         AT THE END OF CVTRVF_POS, FOR SCHEME 13 OR 14.
          c=0.d0
          DO i=1,smi%DIM1
            c = c + dt*smi%R(i)*f%R(i)*volu2d%R(i)
          ENDDO
          IF(ncsize.GT.1) c=p_sum(c)
          massou=massou+c
        ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
!
!     COMPUTES THE TRACER FLUX AT THE BOUNDARY
!
      IF(.NOT.fv_scheme) THEN
        DO i=1,mesh%NPTFR
          n=mesh%NBOR%I(i)
          flbortra%R(i)=flbor%R(i)*(tetat*f%R(n)+(1.d0-tetat)*fn%R(n))
        ENDDO
!     ELSE
!       FLBORTRA ALREADY COMPUTED
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(ielmf.NE.ielms) CALL chgdis(sm  ,ielmf,ielms,mesh)
      IF(ielmf.NE.ielmh.AND.yasmh) CALL chgdis(smh ,ielmf,ielmh,mesh)
!
!-----------------------------------------------------------------------
!
      RETURN
      END