suspension_computation.f

Go to the documentation of this file.

!                   *********************************
                    SUBROUTINE suspension_computation
!                   *********************************
!
     &(slvtra, hn,hn_tel,uconv, vconv, mu,tob,fdm, fd90, ksp,ksr,ks,
     & elay, ava, afbor, bfbor, limdif, clt, maskel, masktr,
     & maskpt, ifamas, npoin, ielm, nptfr, itra, lt, nit, resol,
     & optban, kent,kddl,kdir,ksort,klog,kneu,
     & optadv, opdtra, debug,csf_sable,
     & teta_susp, dt, mased0, zero, xwc, karman, xmve, xmvs, vce,grav,
     & hmin, vitcd, partheniades, entets,
     & bilma,msk,charr,imp_inflow_c,mesh,zf,cs,
     & cst,ctild,cbor,disp,it1,it2,it3,it4,tb,t1,t2,t3,
     & t4, t8, t9, t10, t11, t12, t14, te1, te2, te3, s,
     & am1_s, am2_s, mbor,masten, mastou, masini, ac,
     & zfcl_s, fludpt, fludp, fluer, hprop, disp_c, cstaeq, csratio,
     & masfin, masdept, masdep, massou,qs_c,icq, zref,
     & corr_conv,u2d,v2d,sedco,dift,dm1,zconv,uconv_tel,vconv_tel,
     & solsys,flbor_tel,flbor_sis,flbortra,code,
     & volu2d,v2dpar,unsv2d,numliq,nfrliq,licbor,mixte,avail,nsicla,
     & es,es_sable,es_vase,nomblay,conc,toce_vase,toce_sable,
     & fluer_vase,toce_mixte,ms_sable,ms_vase,dirflu,qsclxs,qsclys,
     & maxadv)
!
!***********************************************************************
! SISYPHE   V7P1
!***********************************************************************
!
!brief    MAIN SUBROUTINE FOR THE COMPUTATION OF THE
!+                CONCENTRATION AND THE ELEVATION SOLVING THE EQUATION :
!code
!+      D(ZF)
!+      ----  + DIV(QS) + (E-D)ZA = 0
!+      DT
!
!note     IF COUPLING, DIV(QS) ALREADY COMPUTED
!+         ELSE,        DIV(QS) = 0
!
!history  F. HUVELIN
!+        22/12/2004
!+
!+
!
!history  JMH:
!+        10/11/2010
!+
!+   ENTET CHANGED INTO ENTETS IN THE CALL TO CVDFTR
!
!history
!+        05/05/2008
!+
!+   ADAPTED FOR FINITE VOLUME ADVECTION
!
!history
!+        09/05/2008
!+
!+   FLUDP REMOVED FROM SUSPENSION_FLUX, SUSPENSION_NERBED DELETED
!
!history
!+        28/05/2008
!+
!+   NEW SUSPENSION_BILAN WITH FLUXES THROUGH BOUNDARIES
!
!history
!+        09/06/2008
!+
!+   NEW SUSPENSION_BILAN WITH FLBORTRA GIVEN BY CVDFTR
!
!history
!+        12/06/2008
!+
!+   SECTIONS "TREATING SMALL DEPTHS" AND
!
!history
!+        25/06/2008
!+
!+   CALLS DIFFIN (USED TO BE IN SUSPENSION_MAIN)
!
!history
!+        31/07/2008
!+
!+   CALLS SUSPENSION_FLUX SPLIT IN 2 : DEPOSITION + EROSION
!
!history
!+        16/09/2009
!+
!+   AVAIL(NPOIN,10,NSICLA)
!
!history
!+        05/04/2010
!+        V6P0
!+   CSTAEQ TAKES INTO ACCOUNT THE % OF LAYER QQ OR THE SELECTED FORMULATION
!
!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)
!+        19/04/2011
!+        V6P1
!+   Adaptation to the new call of Sisyphe in Telemac-2D, and various
!+   modifications for mass-conservation.
!
!history  C.VILLARET (EDF-LNHE), P.TASSI (EDF-LNHE)
!+        19/07/2011
!+        V6P1
!+   Name of variables
!
!history  MAK (HRW)
!+        31/05/2012
!+        V6P2
!+  Include CSRATIO
!
!history  PAT (P. TASSI) (EDF & LNHE)
!+        18/06/2012
!+        V6P2
!+   updated version with HRW's development Soulsby-van Rijn's concentration
!
!history  C. VILLARET (EDF & LNHE)
!+        18/06/2012
!+        V6P2
!+   modification call to suspension_erosion_coh (simplified)
!+   arguments for M_VASE in double precision
!+   modification for limitation of FLUER with rigid lid
!+   calling to suspension_bilan_coh
!
!history  J.-M. HERVOUET (LNHE)
!+        15/04/2013
!+        V6P3
!+   YAFLULIM was not initialised in one case.
!
!history  J.-M. HERVOUET (LNHE)
!+        18/11/2013
!+        V6P3
!+   Pointers GLOSEG1 and GLOSEG2 added to avoid an unwanted copy of
!+   arrays in call to suspension_conv (after message by Intel compiler)
!
!history  J.-M. HERVOUET (LNHE)
!+        28/04/2014
!+        V7P0
!+   Call to diffin modified.
!+   OPTSUP replaced by OPTADV in the call to cvdftr.
!+   (see keyword SCHME OPTION FOR ADVECTION)
!
!history  J.-M. HERVOUET (LNHE)
!+        28/05/2015
!+        V7P1
!+   Call to cvdftr modified, 3 extra arguments
!+   (for implicit distributive schemes)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AC             |<->| CRITICAL SHIELDS PARAMETER
!| ACLADM         |-->| MEAN DIAMETER OF SEDIMENT
!| AFBOR          |-->| BOUNDARY CONDITION ON F: NU*DF/DN=AFBOR*F+BFBOR
!| AM1_S          |<->| MATRIX OBJECT
!| AM2_S          |<->| MATRIX OBJECT
!| AVAIL          |<->| VOLUME PERCENT OF EACH CLASS
!| BFBOR          |-->| BOUNDARY CONDITION ON F: NU*DF/DN=AFBOR*F+BFBOR
!| BILMA          |-->| MASS BALANCE
!| CBOR           |<->| IMPOSED SUSPENDED SAND CONCENTRATION AT THE BOUNDARY
!| CF             |-->| QUADRATIC FRICTION COEFFICIENT
!| CHARR          |-->| LOGICAL, BEDLOAD OR NOT
!| CLT            |<->| BOUNDARY CONDITIONS FOR TRACER (MODIFIED LITBOR)
!| CODE           |-->| HYDRODYNAMIC CODE IN CASE OF COUPLING
!| CONC_VASE      |<->| MUD CONCENTRATION FOR EACH LAYER
!| CORR_CONV      |-->| CORRECTION ON CONVECTION VELOCITY
!| CS             |<->| CONCENTRATION AT TIME N
!| CSF_SABLE      |-->| VOLUME CONCENTRATION OF THE SAND BED
!| CSRATIO        |<->| EQUILIBRIUM CONCENTRATION FOR SOULSBY-VAN RIJN EQ.
!| CST            |<->| CONCENTRATION AT TIME T(N+1)
!| CSTAEQ         |<->| EQUILIBRIUM CONCENTRATION
!| CTILD          |<->| CONCENTRATION AFTER ADVECTION
!| DEBUG          |-->| FLAG FOR DEBUGGING
!| DIFT           |-->| DIFFUSION OF SUSPENDED SEDIMENT CONCENTRATION
!| DISP           |-->| VISCOSITY COEFFICIENTS ALONG X,Y AND Z .
!|                |   | IF P0 : PER ELEMENT
!|                |   | IF P1 : PERR POINT
!| DISP_S         |<->| WORK ARRAY FOR SAVING DISPC
!| DM1            |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV)
!| DTS            |-->| TIME STEP FOR SUSPENSION
!| ELAY           |<->| THICKNESS OF EACH LAYER
!| ENTET          |<->| LOGICAL, IF YES INFORMATION IS GIVEN ON MASS CONSERVATION
!| ENTETS         |-->| LOGICAL, IF YES INFORMATION IS GIVEN ON MASS CONSERVATION FOR SUSPENSION
!| ES             |<->| LAYER THICKNESSES AS DOUBLE PRECISION
!| FLBORTRA       |<->| FLUXES AT BOUNDARIES TRACER
!| FLBOR_SIS      |<->| FLUXES AT BOUNDARIES SISYPHE
!| FLBOR_TEL      |-->| FLUXES AT BOUNDARIES TELEMAC
!| FLUDP          |<->| DEPOSITION FLUX
!| FLUDPT         |<->| DEPOSITION FLUX (IMPLICIT)
!| FLUER          |<->| EROSION FLUX
!| FLUER_VASE     |<->| FOR MIXED SEDIMENTS
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| HMIN           |-->| MINIMUM VALUE OF WATER DEPTH
!| HN             |-->| WATER DEPTH
!| HN_TEL         |-->| WATER DEPTH AS SENT BY TELEMAC OR CALLING CODE
!| HPROP          |<->| PROPAGATION DEPTH (DONE IN CVDFTR)
!| ICQ            |-->| REFERENCE CONCENTRATION FORMULA
!| IELM           |-->| TYPE OF ELEMENT
!| IFAMAS         |-->| A MODIFIED IFABOR WHEN ELEMENTS ARE MASKED
!| IMP_INFLOW_C   |-->| IMPOSED CONCENTRATION IN INFLOW
!| IT1            |<->| INTEGER WORK ARRAY IN A BIEF_OBJ STRUCTURE
!| IT2            |<->| INTEGER WORK ARRAY IN A BIEF_OBJ STRUCTURE
!| IT3            |<->| INTEGER WORK ARRAY IN A BIEF_OBJ STRUCTURE
!| IT4            |<->| INTEGER WORK ARRAY IN A BIEF_OBJ STRUCTURE
!| KARMAN         |-->| VON KARMAN CONSTANT
!| KDDL           |-->| CONVENTION FOR DEGREE OF FREEDOM
!| KDIR           |-->| CONVENTION FOR DIRICHLET POINT
!| KENT           |-->| CONVENTION FOR LIQUID INPUT WITH PRESCRIBED VALUE
!| KINC           |-->| CONVENTION FOR INCIDENT WAVE BOUNDARY CONDITION
!| KLOG           |-->| CONVENTION FOR SOLID BOUNDARY
!| KNEU           |-->| CONVENTION FOR NEUMANN CONDITION
!| KS             |-->| TOTAL BED ROUGHNESS
!| KSORT          |-->| CONVENTION FOR FREE OUTPUT
!| KSP            |-->| SKIN BED ROUGHNESS
!| KSR            |-->| RIPPLE BED ROUGHNESS
!| KX             |<->| COEFFICIENTS OF THE DISPERSION TENSOR (DIM. NPOIN)
!| KY             |<->| COEFFICIENTS OF THE DISPERSION TENSOR (DIM. NPOIN)
!| KZ             |<->| COEFFICIENTS OF THE DISPERSION TENSOR (DIM. NPOIN)
!| LICBOR         |-->| BOUNDARY CONDITIONS FOR SEDIMENT
!| LIMDIF         |<->| BOUNDARY CONDITIONS FOR DIFFUSION
!| LT             |-->| ITERATION
!| MASDEP         |<--| TOTAL DEPOSITED MASS
!| MASDEPT        |<--| DEPOSITED MASS DURING THE TIME STEP
!| MASED0         |<->| SUSPENDED MASS BALANCE
!| MASFIN         |<--| MASS AT THE END
!| MASINI         |<->| INITIAL MASS
!| MASKEL         |-->| MASKING OF ELEMENTS
!| MASKPT         |-->| MASKING PER POINT
!| MASKTR         |<->| MASKING FOR TRACERS, PER POINT
!| MASSOU         |<--| MASS OF TRACER ADDED BY SOURCE TERM
!|                |   | SEE DIFSOU
!| MASTCP         |<--| ??? NE SERT A RIEN, A SUPPRIMER
!| MASTEN         |<->| MASS ENTERED THROUGH LIQUID BOUNDARY
!| MASTOU         |<->| MASS CREATED BY SOURCE TERM
!| MAXADV         |<->| MAXIMUM NUMBER OF ITERATIONS OF ADVECTION SCHEMES
!| MBOR           |<->| MATRIX OBJECT
!| MESH           |<->| MESH STRUCTURE
!| MIXTE          |-->| MIXTURE OF COHESIVE AND NON COHESIVE SEDIMENT
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS
!| MS_SABLE       |<->| MASS OF SAND PER LAYER (KG/M2)
!| MS_VASE        |<->| MASS OF MUD PER LAYER (KG/M2)
!| MU             |-->| CORRECTION FACTOR FOR BED ROUGHNESS
!| NOMBLAY        |-->| NUMBER OF LAYERS
!| NFRLIQ         |-->| NUMBER OF LIQUID BOUNDARIES
!| NIT            |-->| TOTAL NUMBER OF ITERATIONS
!| NPOIN          |-->| NUMBER OF POINTS
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| NSICLA         |-->| NUMBER OF SIZE CLASSES FOR BED MATERIALS
!| NUMLIQ         |-->| LIQUID BOUNDARY NUMBER OF BOUNDARY POINTS
!| OPDTRA         |-->| OPTION FOR THE DIFFUSION OF TRACERS
!| OPTBAN         |-->| OPTION FOR THE TREATMENT OF TIDAL FLATS
!| OPTDIF         |-->| OPTION FOR THE DISPERSION
!| OPTADV         |-->| SCHEME OPTION FOR ADVECTION
!| PARTHENIADES   |-->| CONSTANT OF THE KRONE AND PARTHENIADES EROSION LAW (M/S)
!| PASS           |<->| IN FACT PASS_SUSP IN SISYPHE.F, ARRIVES AS .TRUE.
!|                |   | AT FIRST CALL AND IS CHANGED INTO .FALSE. BELOW
!| QSCLXS         |<->| TRANSPORT RATE FOR EACH CLASS X-DIRECTION
!| QSCLYS         |<->| TRANSPORT RATE FOR EACH CLASS Y-DIRECTION
!| QSCL_S         |<->| SUSPENDED LOAD TRANSPORT RATE
!| QSXS           |<->| SOLID DISCHARGE X (SUSPENSION)
!| QSYS           |<->| SOLID DISCHARGE Y (SUSPENSION)
!| QS_C           |-->| BEDLOAD TRANSPORT RATE
!| QS_S           |<->| SUSPENDED TRANSPORT RATE
!| RESOL          |-->| CHOICE OF ADVECTION SCHEME
!| S              |<->| VOID STRUCTURE
!| SEDCO          |-->| LOGICAL, SEDIMENT COHESIVE OR NOT
!| SLVTRA         |<->| SLVCFG STRUCTURE
!| SOLSYS         |-->| SLVCFG STRUCTURE
!| T1             |<->| WORK BIEF_OBJ STRUCTURE
!| T10            |<->| WORK BIEF_OBJ STRUCTURE
!| T11            |<->| WORK BIEF_OBJ STRUCTURE
!| T12            |<->| WORK BIEF_OBJ STRUCTURE
!| T2             |<->| WORK BIEF_OBJ STRUCTURE
!| T3             |<->| WORK BIEF_OBJ STRUCTURE
!| T4             |<->| WORK BIEF_OBJ STRUCTURE
!| T5             |<->| WORK BIEF_OBJ STRUCTURE
!| T6             |<->| WORK BIEF_OBJ STRUCTURE
!| T7             |<->| WORK BIEF_OBJ STRUCTURE
!| T8             |<->| WORK BIEF_OBJ STRUCTURE
!| T9             |<->| WORK BIEF_OBJ STRUCTURE
!| TASS           |-->| CONSOLIDATION TAKEN INTO ACCOUNT
!| TB             |-->| BLOCK OF WORKING ARRAYS
!| TE1            |<->| WORKING ARRAY FOR ELEMENTS
!| TE2            |<->| WORKING ARRAY FOR ELEMENTS
!| TE3            |<->| WORKING ARRAY FOR ELEMENTS
!| TETA_SUSP      |<->| IMPLICITATION FACTOR FOR THE DEPOSITION FLUX AND DIFFUSION
!| TOB            |-->| BED SHEAR STRESS (TOTAL FRICTION)
!| TOCE_SABLE     |<->| CRITICAL BED SHEAR STRESS OF SAND
!| TOCE_MIXTE     |<->| CRITICAL BED SHEAR STRESS OF THE MIXED SEDUIMENT PER LAYER
!| TOCE_VASE      |<->| CRITICAL EROSION SHEAR STRESS OF THE MUD PER LAYER (N/M2)
!| U2D            |-->| MEAN FLOW VELOCITY X-DIRECTION
!| UCONV          |<->| X-COMPONENT ADVECTION FIELD (SISYPHE)
!| UCONV_TEL      |-->| X-COMPONENT ADVECTION FIELD (TELEMAC)
!| UNSV2D         |-->| INVERSE OF INTEGRALS OF TEST FUNCTIONS
!| V2D            |-->| MEAN FLOW VELOCITY Y-DIRECTION
!| V2DPAR         |-->| INTEGRAL OF TEST FUNCTIONS, ASSEMBLED IN PARALLEL
!| VCE            |-->| FLOW VISCOSITY
!| VCONV          |<->| Y-COMPONENT ADVECTION FIELD (SISYPHE)
!| VCONV_TEL      |-->| Y-COMPONENT ADVECTION FIELD (TELEMAC)
!| VISC_TEL       |-->| VELOCITY DIFFUSIVITY (TELEMAC)
!| VITCD          |-->| CRITICAL SHEAR VELOCITY FOR MUD DEPOSITION
!| VOLU2D         |-->| INTEGRAL OF BASES
!| W1             |<->| WORKING ARRAY
!| XKX            |-->| COEFFICIENT USED FOR COMPUTING THE DISPERSION
!|                |   | DEPENDS OF OPTIONS
!| XKY            |-->| COEFFICIENT USED FOR COMPUTING THE DISPERSION
!|                |   | DEPENDS OF OPTIONS
!| XMVE           |-->| FLUID DENSITY
!| XMVS           |-->| SEDIMENT DENSITY
!| XWC            |-->| SETTLING VELOCITIES
!| ZCONV          |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV)
!| ZERO           |-->| ZERO
!| ZF             |-->| ELEVATION OF BOTTOM
!| ZFCL_S         |<->| BED EVOLUTION PER CLASS, DUE TO SUSPENDED SEDIMENT
!| ZF_S           |<->| ACCUMULATED BED EVOLUTION DUE TO SUSPENDED SEDIMENT
!| ZREF           |<->| REFERENCE ELEVATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe,
     &    ex_suspension_computation => suspension_computation
      USE bief
      USE declarations_sisyphe, ONLY : flulim,tb2,nco_dist,nsp_dist
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      type(slvcfg),    INTENT(INOUT) :: slvtra
      type(bief_obj),  INTENT(IN)    :: zf,volu2d,v2dpar,unsv2d
      type(bief_obj),  INTENT(IN), TARGET    :: hn,hn_tel
      type(bief_obj),  INTENT(INOUT) :: uconv,vconv
      type(bief_obj),  INTENT(IN)    :: mu,ksp,ksr,ks
      type(bief_obj),  INTENT(IN)    :: tob,licbor
      type(bief_obj),  INTENT(INOUT) :: elay
      type(bief_obj),  INTENT(IN)    :: afbor,bfbor
      type(bief_obj),  INTENT(IN)    :: maskel,maskpt,ifamas
      type(bief_obj),  INTENT(INOUT) :: masktr,limdif,clt
      INTEGER,          INTENT(IN)    :: NPOIN,IELM,NPTFR,ITRA,LT
      INTEGER,          INTENT(IN)    :: NIT,RESOL,OPTBAN,KENT,KDDL
      INTEGER,          INTENT(IN)    :: KDIR,OPTADV,OPDTRA,SOLSYS
      INTEGER,          INTENT(IN)    :: KSORT,KLOG,KNEU
      INTEGER,          INTENT(IN)    :: NFRLIQ,NSICLA,NOMBLAY
      INTEGER,          INTENT(IN)    :: DEBUG,DIRFLU,MAXADV
      INTEGER,          INTENT(IN)    :: NUMLIQ(*)
      DOUBLE PRECISION, INTENT(IN)    :: TETA_SUSP, DT, MASED0
      DOUBLE PRECISION, INTENT(IN)    :: XWC,FDM,FD90
      DOUBLE PRECISION, INTENT(IN)    :: CSF_SABLE,AVA(npoin)
      DOUBLE PRECISION, INTENT(IN)    :: KARMAN, XMVE, XMVS,VCE, GRAV
      DOUBLE PRECISION, INTENT(IN)    :: VITCD,PARTHENIADES,HMIN
      LOGICAL,          INTENT(IN)    :: ENTETS,BILMA,MSK,SEDCO
      LOGICAL,          INTENT(IN)    :: CHARR, IMP_INFLOW_C,CORR_CONV
      LOGICAL,          INTENT(IN)    :: DIFT,MIXTE
      type(bief_mesh), INTENT(INOUT) :: mesh
      type(bief_obj),  INTENT(INOUT) :: cs,cst,ctild,cbor,flbor_sis
      type(bief_obj),  INTENT(INOUT) :: disp,it1,it2,it3,it4,tb
      type(bief_obj),  INTENT(INOUT) :: t2, t3, t4, t8
      type(bief_obj),  INTENT(INOUT), TARGET :: t1
      type(bief_obj),  INTENT(INOUT) :: t9, t10, t11, t12, t14, te1
      type(bief_obj),  INTENT(INOUT) :: te2, te3, s, am1_s, am2_s
      type(bief_obj),  INTENT(INOUT) :: mbor,zref
      DOUBLE PRECISION, INTENT(INOUT) :: MASTEN, MASTOU, MASINI, AC
      type(bief_obj),  INTENT(INOUT) :: zfcl_s
      type(bief_obj),  INTENT(IN)    :: uconv_tel,vconv_tel
      type(bief_obj),  INTENT(INOUT) :: fludpt,fludp,fluer,flbortra
      type(bief_obj),  INTENT(INOUT) :: hprop, disp_c, cstaeq,csratio
      type(bief_obj),  INTENT(INOUT) :: fluer_vase,toce_mixte
      type(bief_obj),  INTENT(INOUT) :: qsclxs,qsclys
      DOUBLE PRECISION, INTENT(INOUT) :: MS_SABLE(*)
      DOUBLE PRECISION, INTENT(INOUT) :: MS_VASE(*)
      DOUBLE PRECISION, INTENT(INOUT) :: ES_SABLE(*)
      DOUBLE PRECISION, INTENT(INOUT) :: ES_VASE(*)
      DOUBLE PRECISION, INTENT(INOUT) :: MASFIN,MASDEPT,MASDEP
      DOUBLE PRECISION, INTENT(IN)    :: ZERO
      DOUBLE PRECISION, INTENT(INOUT) :: MASSOU
      DOUBLE PRECISION, INTENT(INOUT) :: AVAIL(npoin,nomblay,nsicla)
      DOUBLE PRECISION, INTENT(INOUT) :: ES(npoin,nomblay),TOCE_SABLE
      DOUBLE PRECISION, INTENT(INOUT) :: CONC(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: TOCE_VASE(nomblay)
      type(bief_obj),  INTENT(IN)    :: qs_c,u2d,v2d,dm1,zconv
      type(bief_obj),  INTENT(IN)    :: flbor_tel
      INTEGER,          INTENT(IN)    :: ICQ
      CHARACTER(LEN=24), INTENT(IN)   :: CODE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER          :: I,K,SOLSYS_SIS,OPTVF,BID(1),RESOL_MOD,IELMT
      DOUBLE PRECISION :: TETAH,AGGLOT
      LOGICAL          :: YASMI2,YAFLULIM
      type(bief_obj),  POINTER :: hold
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: SAVE_UCONV,SAVE_VCONV
      DOUBLE PRECISION :: MSTOT
      DOUBLE PRECISION :: CONC_SABLE(npoin, nomblay)
      INTEGER          :: J
!
      INTEGER, POINTER, DIMENSION(:) :: GLOSEG1,GLOSEG2
!
!-----------------------------------------------------------------------
!
!     IN CHARAC IELMT IS INTENT(INOUT)
      ielmt=ielm
!
!     UCONV POINTER SAVED BEFORE PLAYING WITH IT
!
      save_uconv=>uconv%R
      save_vconv=>vconv%R
      gloseg1=>mesh%GLOSEG%I(1:mesh%GLOSEG%DIM1)
      gloseg2=>mesh%GLOSEG%I(mesh%GLOSEG%DIM1+1:2*mesh%GLOSEG%DIM1)
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
!     TAKES DETAILS OF THE CONTINUITY EQUATION INTO ACCOUNT
!     IN TELEMAC-2D OR 3D, WITH SOLSYS=2, DM1 AND ZCONV ARE USED.
!
      IF(code(1:9).EQ.'TELEMAC2D') THEN
        solsys_sis=solsys
      ELSEIF(lt.GT.1.AND.code(1:9).EQ.'TELEMAC3D') THEN
!       CALL TO SISYPHE TO BE MOVED IN THE TELEMAC-3D TIME LOOP
        solsys_sis=solsys
      ELSE
        solsys_sis=1
      ENDIF
!
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! 1.  COMPUTES THE REFERENCE ELEVATION  -->  ZREF
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!     THREE OPTIONS : ICQ=1: FREDSOE REFERENCE CONC. ZREF = 2.D50
!                     ICQ=2: BIJKER METHOD ZREF = MAX(KSP,KS)
!                     ICQ=3: VAN RIJN ZREF= 0.5 KS
!
      IF(icq.EQ.1) THEN
        CALL os('X=Y     ', x=zref, y=ksp)
      ELSEIF(icq.EQ.2) THEN
        CALL os('X=Y     ', x=zref, y=ksr)
      ELSEIF(icq.EQ.3) THEN
        CALL os('X=CY    ', x=zref, y=ks,c=0.5d0)
      ELSEIF(icq.EQ.4) THEN
        CALL os('X=CY    ', x=zref, y=ks,c=0.5d0)
      ELSE
        WRITE(lu,201) icq
201     FORMAT(1x,'SUSPENSION_COMPUTATION:',/,1x,
     &            'REFERENCE CONCENTRATION FORMULA',/,1x,
     &            'UNEXPECTED VALUE:',1i6)
        CALL plante(1)
        stop
      ENDIF
!
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! 2.  ADVECTION VELOCITY -->  UCONV, VCONV
!     TAKING INTO ACCOUNT THE VERTICAL PROFILE
!     OF CONCENTRATIONS AND VELOCITIES
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!     OPTVF : TENS                  0 : NORMAL
!                                   1 : ADVECTION FIELD DOES NOT SATISFY
!                                       CONTINUITY
!
!     OPTVF : UNITS                 0 : CONSTANT = 0
!                                   1 : CHI-TUAN CONSTANT
!                                   2 : LEO POSTMA CONSTANT
!                                   SEE CVTRVF IN BIEF AND
!                                   V5.7 RELEASE NOTES
!
      yaflulim=.false.
!
      IF(corr_conv.AND.(.NOT.sedco)) THEN
!
        CALL cpstvc(u2d,t12)
        CALL suspension_conv(tob,xmve, ksr,npoin,zref,u2d,v2d,hn,
     &                       uconv,vconv,karman,zero,xwc,t12,resol,
     &                       gloseg1,gloseg2,mesh%NSEG,flulim,
     &                       yaflulim,solsys_sis,solsys,
     &                       uconv_tel,vconv_tel)
!
!       ADVECTION FORM WHICH ACCEPTS AN ADVECTION FIELD
!       THAT DOES NOT SATISFY CONTINUITY + LEO-POSTMA CONSTANT
!
!       WITH 12: MASS CONSERVATION BUT NO MONOTONICITY
!                THE CORRECT THEORY
        optvf=12
!
!       WITH 2: MONOTONICITY BUT NO MASS CONSERVATION
!               WRONG THEORY
!       OPTVF=2
!
!       OPTVF=2 IS POSSIBLE BUT WITH MASS CONSERVATION SPOILED
!       THE UNIT (HERE 2) IS REDONE IN CVDFTR ACCORDING TO THE
!       VALUE OF RESOL, SO IT IS NOT IMPORTANT HERE.
!
      ELSE
!
!       POINTERS ARE USED TO AVOID COPY
!
        IF(solsys_sis.EQ.1) THEN
          uconv%R=>u2d%R
          vconv%R=>v2d%R
        ELSE
!         HERE UCONV_TEL IS PASSED ON
          uconv%R=>uconv_tel%R
          vconv%R=>vconv_tel%R
        ENDIF
!       ADVECTION FORM THAT REQUIRES AN ADVECTION FIELD
!       THAT SATISFIES CONTINUITY + LEO-POSTMA CONSTANT
        optvf=2
!
      ENDIF
!
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! 3.  EROSION FLUX   : FLUER
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!
!     SKIN FRICTION TAUP  --> T4
!
      CALL os('X=CYZ   ', x= t4, y= tob, z= mu, c=1.d0)
      CALL os('X=+(Y,C)', x=t4, y=t4, c=zero)
!
!     SAND ONLY
!
      IF(.NOT.mixte) THEN
        IF(.NOT.sedco) THEN
          IF (debug > 0) WRITE(lu,*) 'SUSPENSION_EROSION'
          CALL suspension_erosion(t4,hn,fdm,fd90,ava,npoin,charr,xmve,
     &                            xmvs,vce,grav,xwc,zero,
     &                            zref,ac,fluer,cstaeq,qs_c,icq,u2d,v2d,
     &                            csratio,debug)
          IF (debug > 0) WRITE(lu,*) 'END_SUSPENSION_EROSION'
!
!         TODO: NOTE JMH : THIS SHOULD BE INCLUDED IN SUSPENSION_EROSION
!
          DO i=1,npoin
            fluer%R(i)=min(fluer%R(i),elay%R(i)*ava(i)/dt*csf_sable)
          ENDDO
!
!       MUD ONLY
!
        ELSE
          CALL suspension_erosion_coh(t4,npoin,xmvs,
     &                      partheniades,fluer,
     &                      toce_vase,nomblay,dt,ms_vase)
!
          IF(nomblay.EQ.1) THEN
            DO i=1,npoin
              fluer%R(i)=min(fluer%R(i),ms_vase(i)/dt/xmvs)
            ENDDO
          ELSE
            DO i=1,npoin
              mstot=0.d0
              DO j=1,nomblay
                mstot=mstot+ms_vase(i+(j-1)*npoin)
              ENDDO
              fluer%R(i)=min(fluer%R(i),mstot/dt/xmvs)
            ENDDO
          ENDIF
        ENDIF
!
!       MIXED SEDIMENT
!       FIRST CLASS= SAND, SECOND CLASS = MUD
!
      ELSE
        IF(.NOT.sedco) THEN
          IF(debug > 0) WRITE(lu,*) 'SUSPENSION_FLUX_MIXTE'
          CALL suspension_flux_mixte(t4,fdm,npoin,charr,xmve,xmvs,
     &                               vce,grav,xwc,zero,
     &                               partheniades,fluer,fluer_vase,
     &                               zref,ac,cstaeq,qs_c,icq,debug,
     &                               avail,nsicla,es,toce_vase,
     &                               toce_sable,nomblay,
     &                               dt,toce_mixte%R,ms_sable,
     &                               ms_vase)
          IF (debug > 0) WRITE(lu,*) 'END_SUSPENSION_FLUX_MOY'
        ENDIF
        IF(sedco) CALL os('X=Y     ',x=fluer, y=fluer_vase)
      ENDIF
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!  4. DEPOSITION FLUX   : FLUDPT =WC*T2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!---> FLUDPT: IMPLICIT TERM
! TO ADD ? mak change ---> T2    : RATIO BETWEEN BOTTOM CONCENTRATION AND AVERAGE
!--->  CSRATIO   : RATIO BETWEEN BOTTOM CONCENTRATION AND AVERAGE
!             CONCENTRATION
!
      IF (debug > 0) WRITE(lu,*) 'SUSPENSION_DEPOT'
      CALL suspension_depot(tob,hn,npoin,hmin,xwc,vitcd,zero,karman,
     &    fdm,fd90,xmve,t1,t2,zref,fludpt,debug,sedco)
! TO ADD? mak   &    FDM,FD90,XMVE,T1,CSRATIO,T14,ZREF,FLUDPT,DEBUG,SEDCO,U2D,V2D,
! TO ADD? mak   &    CSTAEQ,DT)
!     &                      XMVE,T1,T2,ZREF,FLUDPT,DEBUG,SEDCO)
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!  5. DIFFIN A SPECIFIC TREATMENT IS DONE IF THE ADVECTION METHOD
!     IS THE CHARACTERISTICS: FREE OUTPUTS ARE TREATED LIKE DIRICHLET.
!     THIS SPECIFIC TREATMENT IS CANCELLED HERE BY SENDING A MODIFIED
!     VALUE FOR RESOL : RESOL_MOD (IN DIFFIN THE ONLY TEST IS:
!     IF(RESOL.EQ.1) THEN .... ELSE ....  ENDIF)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      resol_mod=resol
      IF(resol_mod.EQ.1) resol_mod=2
      IF (debug > 0) WRITE(lu,*) 'DIFFIN'
      CALL diffin(masktr,limdif%I,licbor%I,clt%I,u2d%R,v2d%R,
     &            mesh%XNEBOR%R,mesh%YNEBOR%R,
     &            mesh%NBOR%I,nptfr,
     &            kent,ksort,klog,kneu,kdir,kddl,resol_mod,
     &            mesh%NELBOR%I,npoin,
!                              NFRLIQ
     &            msk,maskel%R,0,
!                  THOMFR FRTYPE
     &            .false.,bid,    cs,cbor,numliq,
     &            mesh%IKLBOR%I,mesh%NELEB,mesh%NELEBX)
      IF (debug > 0) WRITE(lu,*) 'END DIFFIN'
!
!+++++++++++++++++++++++++++++++++++ELEVATION OF BOTTOM++++++++++++++++++++++++++++++++++++
!  6. BOUNDARY CONDITIONS : CBORCONC_V
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!     IMPOSES THE EQUILIBRIUM CONCENTRATION FOR THE INFLOW NODES  !
!     HERE CBOR FROM BOUNDARY CONDITIONS FILE OR SUBROUTINE CONLIT
!     OVERWRITTEN
!
!     T2 = RATIO BETWEEN BOTTOM CONC.
!     AND AVERAGE CONC. MUST BE KEPT UNTIL THIS STAGE
!
      IF (debug > 0) WRITE(lu,*) 'IMP_INFLOW_C'
      IF(imp_inflow_c) THEN
!
        DO k = 1, nptfr
          IF(clt%I(k).EQ.kent) THEN
            i = mesh%NBOR%I(k)
            IF(.NOT.sedco) THEN
              cbor%R(k) = cstaeq%R(i)/t2%R(i)
              IF(mixte) cbor%R(k) = fluer%R(i)/t2%R(i)/xwc
            ELSE
              cbor%R(k) = fluer%R(i)/xwc
            ENDIF
!           THIS IS THE CONDITION TO HAVE NO EVOLUTION
!           CS%R(I) MAY BE DIFFERENT FROM CBOR%R(K) IF UNSTEADY FLOW
!           OR IF DIRFLU.EQ.2 (CASE OF PRIORITY TO FLUXES)
            fluer%R(i)=fludpt%R(i)*cs%R(i)
          ENDIF
        ENDDO
!
      ENDIF
      IF (debug > 0) WRITE(lu,*) 'FIN IMP_INFLOW_C'
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!  7. SOLVING TRANSPORT EQUATION IF METHOD OF CHARACTERISTICS
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      IF(resol == 1) THEN
        IF (debug > 0) WRITE(lu,*) 'CHARAC'
        CALL charac(cs,ctild,1,uconv,vconv,s,s,s,s,dt,ifamas,
     &              ielmt,npoin,1,1,1,
     &              msk,am1_s%X,am1_s%D,am1_s%D,
     &              tb,it1%I,it2%I,it2%I,it3%I,it4%I,it2%I,
     &              mesh,mesh%NELEM,mesh%NELMAX,mesh%IKLE,mesh%SURDET,
     &              am2_s,t14,slvtra,1.d0,entets,3,unsv2d,1)
        IF (debug > 0) WRITE(lu,*) 'END_CHARAC'
      ENDIF
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!  8. SOURCE AND SINKS
!     IMPLICIT SOURCE TERM FOR THE DEPOSITION       : T9
!     EXPLICIT SOURCE TERM WITHOUT PUNCTUAL SOURCES : T11
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      IF(optban.EQ.2) THEN
        CALL os('X=XY    ',x=fluer ,y=maskpt)
      ENDIF
!
      CALL os('X=-Y    ',x=t9,y=fludpt)
      CALL os('X=Y     ',x=t11,y=fluer)
!
      DO i=1,npoin
        IF(hn%R(i).GT.hmin) THEN
          t11%R(i)=t11%R(i)/hn%R(i)
        ELSE
          t11%R(i)=0.d0
!         FLUER WILL BE USED AS T11*HN, SO IT MUST BE
!         CANCELLED ACCORDINGLY, OTHERWISE MASS BALANCE WRONG
          fluer%R(i)=0.d0
        ENDIF
      ENDDO
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!  9. ADVECTION-DISPERSION STEP
!     CONFIGURATION OF ADVECTION
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      tetah  = 1.d0 - teta_susp
      massou = 0.d0
      agglot=1.d0
      yasmi2 = .true.
!
!     BOUNDARY FLUXES MUST BE SPECIFIED TO CVDFTR (FINITE VOLUMES CASE)
!     AND TO SUSPENSION_BILAN
!     SISYPHE ALONE     : THEY MUST BE COMPUTED
!     WHEN COUPLING     : THEY ARE GIVEN BY THE CALLING SUBROUTINE
!                         EXCEPT AT THE 1ST ITERATION
!
!     IF(CODE(1:7).NE.'TELEMAC'.OR.LT.EQ.1) THEN
      IF(code(1:7).NE.'TELEMAC'.OR.
     &  (code(1:9).EQ.'TELEMAC3D'.AND.lt.EQ.1)) THEN
        IF (debug > 0) WRITE(lu,*) 'VECTOR'
        CALL vector(flbor_sis,'=','FLUBDF          ',ielbor(ielmt,1),
!                        HPROP (HERE HPROP=HN, INVESTIGATE)
     &              1.d0,hn   ,hn,hn,uconv,vconv,vconv,
     &              mesh,.true.,masktr%ADR(5)%P)
!                                          5: MASK OF LIQUID BOUNDARIES
!                                             SEE DIFFIN IN BIEF 6.1
        IF (debug > 0) WRITE(lu,*) 'FIN VECTOR'
      ELSE
        CALL os('X=Y     ',x=flbor_sis,y=flbor_tel)
!       MUST ALSO CHANGE BOUNDARY FLUXES IF THE ADVECTION
!       FIELD IS CORRECTED (T12 MUST HAVE BEEN KEPT SINCE
!                           CALL TO SUSPENSION_CONV)
        IF(corr_conv.AND..NOT.sedco) THEN
          CALL osbd('X=CXY   ',flbor_sis,t12,t12,1.d0,mesh)
        ENDIF
      ENDIF
!
!     FINITE VOLUMES ADVECTION USES THE TRUE H FROM THE PREVIOUS STEP
      IF(code(1:7).EQ.'TELEMAC') THEN
        IF(optban.NE.0) THEN
          CALL cpstvc(cst,t1)
!         HN_TEL IS NOT CLIPPED
          DO i=1,npoin
            t1%R(i)=max(hn_tel%R(i),hmin)
          ENDDO
          hold=>t1
        ELSE
          hold=>hn_tel
        ENDIF
      ELSE
!       IN THIS CASE H AND HN ARE CONFUNDED
        hold=>hn
      ENDIF
!
      IF(debug > 0) WRITE(lu,*) 'APPEL DE CVDFTR'
!
      CALL cvdftr
     & (cst, ctild, cs, t2,
!                            H         HTILD
     &  dift, resol, .true., hn, hold, hprop,
     &  uconv,vconv,dm1,zconv,solsys_sis,
!                     TEXP SMH  YASMH   TIMP
     &  disp, disp_c, t11, t2, .false., t9,  yasmi2,am1_s,am2_s,
     &  zf, cbor, afbor, bfbor, limdif, masktr, mesh,
     &  tb, t8, t12,  t4, t10, te1, te2, te3,
     &  kdir,kddl,dt,entets,teta_susp,
!                      BILAN
     &  agglot,entets,.false.,optadv,
     &  1, opdtra, optban, msk, maskel, maskpt, mbor, s,
!               OPTSOU
     &  massou, 1,     slvtra,flbor_sis,volu2d,v2dpar,unsv2d,
     &  optvf,flbortra,
     &  flulim,yaflulim,flulim,.false.,dirflu,.false.,t8    ,0.d0,
!                               RAIN  ,PLUIE ,TRAIN
     &  flulim     ,.false.,maxadv,tb2,nco_dist,nsp_dist)
!       GIVEN_FLUX   FLUX_GIVEN (NOW THE FLUX CAN BE GIVEN, THIS COULD
!       BE AN OPTIMISATION, AS HERE IT IS RECOMPUTED FOR EVERY CLASS...)
!
      IF(debug > 0) WRITE(lu,*) 'END_CVDFTR'
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! 10. BED EVOLUTION DUE TO NET EROSION/DEPOSITUON FLUX
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
      DO i=1,npoin
        fludp%R(i)=fludpt%R(i)*cst%R(i)
      ENDDO
!
!     COMPUTES EVOLUTION AND UPDATES DATA
!     TASS TO BE PASSED IN ARGUMENT
!
!
      IF(.NOT.sedco) THEN
        IF(.NOT.mixte) THEN
          CALL os('X=Y-Z   ', x=zfcl_s, y=fludp, z=fluer)
          CALL os('X=CX    ', x=zfcl_s, c=dt/csf_sable)
        ELSE
          DO i=1, npoin
            DO j= 1, nomblay
              conc_sable(i,j)=xmvs
            ENDDO
          ENDDO
          CALL suspension_evol(zfcl_s,fludp,fluer,dt,
     &                npoin,xmvs,t3,ms_sable,es_sable,
     &                   conc_sable,nomblay)
        ENDIF
      ELSE
        CALL suspension_evol(zfcl_s,fludp,fluer,dt,
     &      npoin,xmvs,t3,ms_vase,es_vase,
     &      conc,nomblay)
      ENDIF
!
!     WRITES OUT THE MIN/MAX VALUES TO THE LISTING
!
      IF(entets) THEN
        IF (debug > 0) WRITE(lu,*) 'SUSPENSION_LISTING'
        CALL suspension_listing(mesh,cst,zfcl_s,uconv,vconv,
     &                          maskel,ielmt,dt,msk,t1)
        IF(debug > 0) WRITE(lu,*) 'END_SUSPENSION_LISTING'
      ENDIF
!
!     MASS-BALANCE FOR THE SUSPENSION
!
      IF(bilma) THEN
        IF(sedco) THEN
          IF (debug > 0) WRITE(lu,*) 'SUSPENSION_BILAN_COH'
          CALL suspension_bilan_coh
     &         (mesh,cst,hn,maskel,ielmt,itra,lt,nit,dt,xmvs,
     &         ms_vase,nomblay,npoin,
     &         massou,mased0,msk,entets,masten,mastou,
     &         masini,t1,t2,t3,masfin,masdept,masdep,agglot,volu2d,
     &         numliq,nfrliq,nptfr,flbortra,sedco)
          IF(debug > 0) WRITE(lu,*) 'END_SUSPENSION_BILAN_COH'
        ELSE
!Modifs CVL
          IF(mixte) THEN
            IF (debug > 0) WRITE(lu,*) 'SUSPENSION_BILAN_COH'
            CALL suspension_bilan_coh
     &           (mesh,cst,hn,maskel,ielmt,itra,lt,nit,dt,xmvs,
     &           ms_sable,nomblay,npoin,
     &           massou,mased0,msk,entets,masten,mastou,
     &           masini,t1,t2,t3,masfin,masdept,masdep,agglot,volu2d,
     &           numliq,nfrliq,nptfr,flbortra,sedco)
            IF(debug > 0) WRITE(lu,*) 'END_SUSPENSION_BILAN_COH'
          ELSE
! fin modifs CVL
            IF (debug > 0) WRITE(lu,*) 'SUSPENSION_BILAN'
            CALL suspension_bilan
     &          (mesh,cst,hn,zfcl_s,maskel,ielmt,itra,lt,nit,
     &           dt,csf_sable,massou,mased0,msk,entets,masten,mastou,
     &           masini,t2,t3,masfin,masdept,masdep,agglot,volu2d,
     &           numliq,nfrliq,nptfr,flbortra)
            IF(debug > 0) WRITE(lu,*) 'END_SUSPENSION_BILAN'
          ENDIF
        ENDIF
      ENDIF
!
!
      CALL os('X=Y     ', x=cs, y=cst)
!
!     NOTE: ARE QSCLXS AND QSCLYS USEFUL ????????
!
      CALL os('X=YZ    ', x=t1, y=uconv, z=hn)
      CALL os('X=YZ    ', x=t2, y=vconv, z=hn)
      CALL os('X=YZ    ', x=qsclxs, y=cs, z=t1)
      CALL os('X=YZ    ', x=qsclys, y=cs, z=t2)
!
!     RESTORING UCONV POINTERS
!
      uconv%R=>save_uconv
      vconv%R=>save_vconv
!
!======================================================================!
!======================================================================!
!
      RETURN
      END