suspension_main.f

Go to the documentation of this file.

!                       **************************
                        SUBROUTINE suspension_main
!                       **************************
!
     &(slvtra,hn,hn_tel,mu,tob,fdm,fd90,ksp,ksr,ks,volu2d,v2dpar,unsv2d,
     & afbor,bfbor,zf,licbor,ifamas,maskel,maskpt,u2d,v2d,nsicla,npoin,
     & nptfr,ielmt,optdif,resol,lt,nit,optban,optadv,opdtra,
     & kent,ksort,klog,kneu,kdir,kddl,debug,
     & dts,csf_sable,zero,grav,xkx,xky,karman,
     & xmve,xmvs,vce,hmin,xwc,vitcd,partheniades,bilma,msk,
     & charr,imp_inflow_c,mesh,zf_s,cs,cst,ctild,cbor,disp,
     & it1,it2,it3,it4,tb,t1,t2,t3,t4,t8,t9,t10,t11,t12,t14,
     & te1,clt,te2,te3,s,am1_s,am2_s,mbor,elay,limdif,
     & masktr, teta_susp, ac, mased0, masini, masten,
     & mastou, es,es_sable, es_vase,avail,  entets, pass, zfcl_s,
     & hprop, fludpt, fludp, fluer, disp_c, kx, ky,
     & kz, uconv, vconv,qsxs, qsys, qsclxs, qsclys, qscl_s,
     & qs_s,qs_c,cstaeq,csratio,icq,mastcp,masfin,masdept,masdep,massou,
     & corr_conv,zref,sedco,visc_tel,code,
     & dift,dm1,uconv_tel,vconv_tel,zconv,solsys,flbor_tel,flbor_sis,
     & flbortra,numliq,nfrliq,mixte,nomblay,conc,
     & toce_vase,toce_sable,fluer_vase,toce_mixte,ms_sable,ms_vase,
     & dirflu,maxadv)
!
!***********************************************************************
! SISYPHE   V6P2                                   18/06/2012
!***********************************************************************
!
!brief    MAIN SUBROUTINE FOR THE SUSPENDED-LOAD TRANSPORT.
!
!history  F. HUVELIN
!+        22/12/2004
!+
!+
!history  JMH
!+        25/06/2008
!+        V5P9
!+   CALL TO DIFFIN MOVED IN SUSPENSION_COMPUTATION
!
!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  C. VILLARET
!+        20/03/2011
!+        V6P1
!+        BUG CORRECTION : SEND FDM(I) instead of ACLADM
!+        suppression of NSOUS, ISOUS
!+        adding VCE, TOCE_SABLE
!+
!history  J.-M. HERVOUET
!+        19/04/2011
!+        V6P1
!+        COMPUTATION OF INITIAL MASS CHANGED
!+
!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 (LNHE)
!+        18/06/2012
!+        V6P2
!+   updated version with HRW's development for Soulsby-van Rijn's concentration
!+       V6P2
!
!history  C. VILLARET
!+        21/08/2012
!+        V6P2
!+    Added call variable to suspension_computation
!
!history  C. VILLARET
!+        28/08/2012
!+        V6P2
!+    Added ES_SABLE and ES_VASE
!+    Replaced CONC_VASE by CONC
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        28/04/2014
!+        V7P0
!+   OPTSUP replaced by OPTADV in the call to suspension_computation.
!+   (see keyword SCHME OPTION FOR ADVECTION)
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        28/03/2017
!+        V7P3
!+   HPROP is not equal to HN in coupling, but to HN_TEL, in Sisyphe
!+   what is called HN is in fact H in Telemac-2D or 3D.
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| 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 AND PER LAYER
!| 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 (KG/M3)
!| CORR_CONV      |-->| LOGICAL, CORRECTION ON CONVECTION VELOCITY OR NOT
!| CS             |<->| CONCENTRATION AT TIME N
!| CSF_SABLE      |-->| VOLUME CONCENTRATION OF SAND (1-POROSITY)
!| 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 TOP ACTIVE LAYER (SANG GRADING ALGORITHM)
!|                |<->| THICKNESS OF THE WHOLE COHESIVE SEDIMENT BED (CONSOLIDATION)
!| ENTET          |<->| LOGICAL, IF YES INFORMATION IS GIVEN ON MASS CONSERVATION
!| ENTETS         |-->| LOGICAL, IF YES INFORMATION IS GIVEN ON MASS CONSERVATION FOR SUSPENSION
!| ES             |<->| THICKNESS OF EACH LAYER (M)
!| ES_VASE        |<->| THICKNESS OF THE MUD LAYER (M)
!| ES_SABLE       |<->| THICKNESS OF THE SAND LAYER (M)
!| FDM            |-->| GRAIN SIZE PER SEDIMENT CLASS
!| FLBORTRA       |<->| FLUXES AT BOUNDARIES TRACER
!| FLBOR_SIS      |<->| FLUXES AT BOUNDARIES SISYPHE
!| FLBOR_TEL      |-->| FLUXES AT BOUNDARIES TELEMAC
!| FLUDP          |<->| DEPOSITION FLUX (M/S)
!| FLUDPT         |<->| DEPOSITION FLUX (IMPLICIT)
!| FLUER          |<->| EROSION FLUX (M/S)
!| FLUER_VASE     |<->| EROSION FLUX (M/S)FOR MIXED SEDIMENTS
!| GRAV           |-->| ACCELERATION OF GRAVITY
!| HMIN           |-->| MINIMUM VALUE OF WATER DEPTH (M)
!| HN             |-->| WATER DEPTH (M)
!| HN_TEL         |-->| WATER DEPTH SENT BY TELEMAC OR CALLING CODE
!| HPROP          |<->| PROPAGATION DEPTH (DONE IN CVDFTR)
!| ICQ            |-->| FLAG FOR REFERENCE CONCENTRATION FORMULA
!| IELMT          |-->| NUMBER OF ELEMENTS
!| IFAMAS         |-->| A MODIFIED IFABOR WHEN ELEMENTS ARE MASKED
!| IMP_INFLOW_C   |-->| LOGICAL, IMPOSED EQUILIBRIUM CONCENTRATION AT THE INFLOW OR NOT
!| 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
!| 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 FOR ADVECTION SCHEMES
!| MBOR           |<->| MATRIX OBJECT
!| MESH           |<->| MESH STRUCTURE
!| MIXTE          |-->| LOGICAL, MIXTE SEDIMENT OR NOT
!| MSK            |-->| IF YES, THERE IS MASKED ELEMENTS
!| MS_SABLE       |<->| MASS OF SAND PER LAYER (KG/M2)
!| MS_VASE        |<->| MASS OF MUD PER LAYERv (KG/M2)
!| MU             |-->| CORRECTION FACTOR FOR BED ROUGHNESS
!| NOMBLAY        |-->| NUMBER OF LAYERS FOR CONSOLIDATION
!| 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         |<->| SUSPENDED LOAD TRANSPORT RATE FOR EACH CLASS X-DIRECTION
!| QSCLYS         |<->| SUSPENDED LOAD 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 LOAD 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
!| 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 SHEAR STRESS FOR SAND (N/M2)
!| TOCE_MIXTE     |<->| CRITICAL SHEAR STRESS FOR MIXED SEDIMENTS (N/M2)
!| 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
!| 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 PER CLASS OF SEDIMENT
!| 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_main => suspension_main
      USE bief
      USE declarations_special
      USE interface_parallel, ONLY : p_dsum
      IMPLICIT NONE
!
      ! 2/ GLOBAL VARIABLES
      ! -------------------
      type(slvcfg),    INTENT(INOUT) :: slvtra
      type(bief_obj),  INTENT(IN)    :: hn,hn_tel,mu,tob
      type(bief_obj),  INTENT(IN)    :: ksp,ksr,ks
      type(bief_obj),  INTENT(IN)    :: volu2d,afbor,bfbor,zf
      type(bief_obj),  INTENT(IN)    :: v2dpar,unsv2d
      type(bief_obj),  INTENT(IN)    :: licbor, ifamas, maskel, maskpt
      type(bief_obj),  INTENT(IN)    :: u2d, v2d,dm1,zconv,flbor_tel
      INTEGER,          INTENT(IN)    :: NSICLA, NPOIN, NPTFR, IELMT
      INTEGER,          INTENT(IN)    :: OPTDIF, RESOL,LT, NIT
      INTEGER,          INTENT(IN)    :: OPTBAN,OPTADV,OPDTRA,NFRLIQ
      INTEGER,          INTENT(IN)    :: KENT, KSORT, KLOG, KNEU
      INTEGER,          INTENT(IN)    :: KDIR,KDDL
      INTEGER,          INTENT(IN)    :: DEBUG,SOLSYS,NOMBLAY,MAXADV
      INTEGER,          INTENT(IN)    :: NUMLIQ(nptfr)
      DOUBLE PRECISION, INTENT(IN)    :: DTS,CSF_SABLE
      DOUBLE PRECISION, INTENT(IN)    :: ZERO,GRAV
      DOUBLE PRECISION, INTENT(IN)    :: FDM(nsicla),FD90(nsicla)
      DOUBLE PRECISION, INTENT(IN)    :: XKX,XKY,KARMAN,VCE
      DOUBLE PRECISION, INTENT(IN)    :: XMVE, XMVS, HMIN, XWC(nsicla)
      DOUBLE PRECISION, INTENT(IN)    :: VITCD
      DOUBLE PRECISION, INTENT(IN)    :: PARTHENIADES
      LOGICAL,          INTENT(IN)    :: BILMA, MSK, CHARR
      LOGICAL,          INTENT(IN)    :: IMP_INFLOW_C
      LOGICAL,          INTENT(IN)    :: SEDCO(nsicla),MIXTE
      TYPE(bief_mesh),  INTENT(INOUT) :: MESH
      type(bief_obj),  INTENT(INOUT) :: zf_s,cs,cst,ctild,cbor
      type(bief_obj),  INTENT(INOUT) :: disp,it1,it2,it3,it4,tb
      type(bief_obj),  INTENT(INOUT) :: t1,t2,t3,t4,t8
      type(bief_obj),  INTENT(INOUT) :: t9,t10,t11,t12,t14,te1,clt
      type(bief_obj),  INTENT(INOUT) :: te2,te3,s,am1_s,am2_s,mbor
      type(bief_obj),  INTENT(INOUT) :: elay, limdif,flbortra
      type(bief_obj),  INTENT(INOUT) :: masktr
      DOUBLE PRECISION, INTENT(INOUT) :: TETA_SUSP, AC(nsicla)
      DOUBLE PRECISION, INTENT(INOUT) :: MASED0(nsicla), MASINI(nsicla)
      DOUBLE PRECISION, INTENT(INOUT) :: MASTEN(nsicla), MASTOU(nsicla)
      DOUBLE PRECISION, INTENT(INOUT) :: ES(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: ES_SABLE(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: ES_VASE(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: CONC(npoin,nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: TOCE_VASE(nomblay)
      DOUBLE PRECISION, INTENT(INOUT) :: TOCE_SABLE
      DOUBLE PRECISION, INTENT(INOUT) :: AVAIL(npoin,nomblay,nsicla)
      LOGICAL,          INTENT(INOUT) :: ENTETS, PASS
      type(bief_obj),  INTENT(INOUT) :: zfcl_s,hprop,zref
      type(bief_obj),  INTENT(INOUT) :: fludpt,fludp,fluer
      type(bief_obj),  INTENT(INOUT) :: disp_c,kx,ky,kz,uconv
      type(bief_obj),  INTENT(INOUT) :: vconv,flbor_sis
      type(bief_obj),  INTENT(INOUT) :: qsxs,qsys,qsclxs,qsclys
      type(bief_obj),  INTENT(INOUT) :: qscl_s,qs_s,cstaeq,csratio
      type(bief_obj),  INTENT(INOUT) :: fluer_vase,toce_mixte
      type(bief_obj),  INTENT(INOUT) :: ms_sable,ms_vase
      type(bief_obj),  INTENT(IN)    :: qs_c,visc_tel
      type(bief_obj),  INTENT(IN)    :: uconv_tel,vconv_tel
      DOUBLE PRECISION, INTENT(OUT)   :: MASTCP(nsicla),MASFIN(nsicla)
      DOUBLE PRECISION, INTENT(OUT)   :: MASDEPT(nsicla),MASDEP(nsicla)
      DOUBLE PRECISION, INTENT(OUT)   :: MASSOU
      INTEGER, INTENT(IN)             :: ICQ,DIRFLU
      LOGICAL, INTENT (IN)            :: CORR_CONV,DIFT
      CHARACTER(LEN=24), INTENT(IN)   :: CODE
!
      ! 3/ LOCAL VARIABLES
      ! ------------------
      INTEGER :: I,J
!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
      IF(pass) THEN
!
        ! *************************  !
        ! III - INITIAL MASS-BALANCE !
        ! *************************  !
!
        IF(bilma) THEN
          DO i = 1,nsicla
!           MUST BE DONE LIKE IN SUSPENSION_BILAN
!           I.E. WITH MASS-LUMPING AGGLOT=1. WHICH IS
!           SET LATER IN SUSPENSION_COMPUTATION...
            CALL os('X=YZ    ',x=t1,y=volu2d,z=cs%ADR(i)%P)
!
            IF(code(1:7).EQ.'TELEMAC') THEN
!             WITH COUPLING, HN-TEL IS THE OLD DEPTH
!                            HN     IS THE NEW DEPTH
              mased0(i) = dots(t1,hn_tel)
            ELSE
!             SISYPHE WITHOUT COUPLING, MASS CONSERVATION
!             DIFFICULT TO CHECK...
              mased0(i) = dots(t1,hn)
            ENDIF
            IF(ncsize.GT.1) mased0(i)=p_dsum(mased0(i))
            masini(i) = mased0(i)
            masten(i) = 0.d0
            mastou(i) = 0.d0
            mastcp(i) = 0.d0
            WRITE(lu,2) i, mased0(i)
          ENDDO
        ENDIF
!
        !----------------------------------------------------------------!
002     FORMAT(1x,'INITIAL QUANTITY IN SUSPENSION FOR CLASS ',
     &         i2,' : ', g16.7, ' M3')
        !----------------------------------------------------------------!
!       END OF IF(PASS)
      ENDIF
!
      pass = .false.
!
      ! ********************************* !
      ! V - COMPUTES THE DISPERSION       !
      ! ********************************* !
      IF (debug > 0) WRITE(lu,*) 'SUSPENSION_DISPERSION'
      CALL suspension_dispersion
     &     (tob,xmve,hn,optdif,npoin,xkx,xky,t1,t2,t3,kx,ky,kz,disp,
     &      u2d,v2d,visc_tel,code)
      IF (debug > 0) WRITE(lu,*) 'END_SUSPENSION_DISPERSION'

      ! ************************************************ !
      ! VI  - COMPUTES THE CONCENTRATION AND EVOLUTION   !
      ! ************************************************ !

      IF(code(1:7).EQ.'TELEMAC') THEN
        CALL os('X=Y     ', x=hprop, y=hn_tel)
      ELSE
        CALL os('X=Y     ', x=hprop, y=hn)
      ENDIF
      DO i = 1, nsicla
        CALL os('X=0     ', x=zfcl_s%ADR(i)%P)
        IF(debug > 0) WRITE(lu,*)
     &               'SUSPENSION_COMPUTATION : ',i,'/',nsicla
        CALL suspension_computation(slvtra,hn,hn_tel,uconv,
     &vconv,mu,tob,fdm(i),fd90(i),ksp,ksr,ks,elay,avail(1:npoin,1,i),
     &afbor,bfbor,limdif,clt,maskel,masktr,maskpt,ifamas,npoin,ielmt,
     &nptfr,i,lt,nit,resol,optban,kent,kddl,kdir,ksort,klog,kneu,
     &optadv,opdtra,debug, csf_sable, teta_susp,dts,
     &mased0(i),zero,xwc(i),karman,xmve,xmvs,vce,grav,hmin,vitcd,
     &partheniades,entets,bilma,
     &msk,charr,imp_inflow_c,mesh,zf,cs%ADR(i)%P,
     &cst%ADR(i)%P,ctild%ADR(i)%P,cbor%ADR(i)%P,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(i),mastou(i),
     &masini(i),ac(i),zfcl_s%ADR(i)%P,fludpt%ADR(i)%P,fludp%ADR(i)%P,
     &fluer%ADR(i)%P, hprop,disp_c,cstaeq, csratio,
     &masfin(i),masdept(i),masdep(i),massou,qs_c,icq,zref,
     &corr_conv,u2d,v2d,sedco(i),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%R,ms_vase%R,dirflu,
     &qsclxs%ADR(i)%P,qsclys%ADR(i)%P,maxadv)
        IF (debug > 0) WRITE(lu,*) 'END_SUSPENSION_COMPUTATION'
!
      ENDDO
!
! FOR MIXTE OR COHESIVE SEDIMENTS ELAY UPDATED
!
! REACTUALISATION DU ELAY ET DES AVAI fait dans suspension_MAIN
      IF(sedco(1).OR.mixte) THEN
        DO i = 1, npoin
          elay%R(i)= 0.d0
          DO j= 1, nomblay
            es(i,j) = es_vase(i,j)
            IF(mixte) THEN
              es(i,j)= es_vase(i,j) + es_sable(i,j)
              IF(es(i,j).GT.1.d-04) THEN
                avail(i,j,1)= es_sable(i,j)/es(i,j)
                avail(i,j,2)= es_vase(i,j)/es(i,j)
!CVL          ELSE
!CVL            AVAIL(I,J,1)=0.5 D0
!CVL            AVAIL(I,J,2)=0.5 D0
              ENDIF
            ENDIF
            elay%R(i)=elay%R(i)+es(i,j)
          ENDDO
        ENDDO
      ENDIF
!
      ! *********************************************************** !
      ! VII  - UPDATES EVOLUTION, CONCENTRATION AND TRANSPORT RATE  !
      ! *********************************************************** !
!
      IF (debug > 0) WRITE(lu,*) 'UPDATING_DATA'
!
!     COULD BE OPTIMISED: FIRST CLASS ON ZF_S, THEN ADDING OTHERS...
!
      CALL os('X=0     ', x=qsxs)
      CALL os('X=0     ', x=qsys)
      CALL os('X=0     ', x=zf_s)
!
      DO i = 1, nsicla
        CALL os('X=X+Y   ', x=zf_s, y=zfcl_s%ADR(i)%P)
        CALL os('X=X+Y   ', x=qsxs, y=qsclxs%ADR(i)%P)
        CALL os('X=X+Y   ', x=qsys, y=qsclys%ADR(i)%P)
      ENDDO
      CALL os('X=N(Y,Z)', x=qscl_s, y=qsclxs, z=qsclys)
      CALL os('X=N(Y,Z)', x=qs_s, y=qsxs, z=qsys)
      IF (debug > 0) WRITE(lu,*) 'END_UPDATING_DATA'
!
!======================================================================!
!======================================================================!
!
      RETURN
      END