MODULE sed_MUSTANG

#if defined key_sedim_MUSTANG

   !&E==========================================================================
   !&E                   ***  MODULE  sed_MUSTANG (option mixsed) ***
   !&E
   !&E
   !&E ** Purpose : concerns subroutines related to sediment dynamics 
   !&E              if it is  used in MARS system (use key_MARS)
   !&E              [ONLY option key_mixsed in MUSTANG (implicitly)
   !&E 
   !&E ** Description :
   !&E     subroutine sed_MUSTANG_first      ! updates forcings and remoddeling sediment  
   !&E     subroutine sed_MUSTANG_deposition     ! sedimentation  sediment
   !&E 
   !&E     subroutine sed_MUSTANG_settlveloc     ! settling velocity in the water column
   !&E     subroutine sed_MUSTANG_depflx         ! computation of deposition tendancy
   !&E     subroutine sed_MUSTANG_sandconcextrap ! extrapolation of sand concentration
   !&E     subroutine sed_MUSTANG_erosion        ! computes erosion thickness
   !&E     subroutine sed_MUSTANG_comp_tocr_mixsed    ! critical shear stress for erosion (Pa)
   !&E     subroutine sed_MUSTANG_comp_eros_flx ! comp_eros_flx : erosion flux (kg/m2/s)
   !&E     subroutine sed_MUSTANG_effdep         ! effective deposition
   !&E     subroutine sed_MUSTANG_fusion         ! fusion of two first bottom layers
   !&E     subroutine sed_MUSTANG_slipdepo       ! slip lateral deposit
   !&E     subroutine sed_MUSTANG_consol_diff_bioturb   
   !&E     subroutine sed_MUSTANG_coefbioturb        ! calcul coefficient bioturbation coef in the sediment
   !&E     subroutine sed_constitutivrel
   !&E
   !&E ** History :
   !&E     ! 2016-03  (B.Thouvenin) : extraction of routines for reorganization of module SEDIMARS
   !&E
   !&E  **********************************************************************
   !&E  *** requires to have advected substances in a specific order       ***
   !&E  ***       in water : cw (iv,k) k=1,kmax bottom to surface          ***
   !&E  ***       in sediment : cs (iv,k) k=ksdmin,ksdmax                  ***
   !&E  ***                                                                ***
   !&E  ***       First nvp particulate substances                         ***
   !&E  ***    first GRAVEL (indices igrav1 to igrav2) : coarser to finer  ***
   !&E  ***    then  SAND (isand1=igrav1+1 to isand2) : coarser to finer   ***
   !&E  ***    then MUDS (imud1=isand2+1 to imud2                          ***
   !&E  ***                                                                ***
   !&E  ***    then nvd Dissolved substances                               *** 
   !&E  ***                                                                ***
   !&E  **********************************************************************
   !&E
   !&E     Used with specific date issued for hydrodynamic module : conversion in coupleur_define_MUSTANG.h
   !&E       kmax : number of layers in water column
   !&E       h0fond : residual water thickness (in m). This thickness is added to the computed water level (xe)
   !&E                to avoid possible mass loss during drying (due to the fact that the scheme for water continuity is not positive)
   !&E       rhoref : reference sea water density (in kg/m3)
   !&E       turb_nbeq : number of equations in the turbulence model(0, 1 ou 2)
   !&E       turb_2eq_option :  option number for the two-equation model : (1 : k-kl; 2 : k-epsilon ; 3 : k-omega ; 4 : modele generique)
   !&E       dt : time step 
   !&E
   !&E     Used with specific date issued for substances module : conversion in coupleur_define_MUSTANG.h
   !&E       nv_adv : number of advected (transported) state variables
   !&E       nv_tot : total number of  variables (state+intermediate+driving..)
   !&E       nv_mud : total number of  variables type MUDS
   !&E       nvp :  number of particulate advected variables
   !&E       nvpc :  number of constitutive particulate advected variables
   !&E       igrav1,igrav2 : first and last indices of gravel
   !&E       isand1,isand2 : first and last indices of sand (first : corsaer sand - last : finer sand)
   !&E       imud1,imud2 :first and last indices of mud
   !&E       typart(1:nv_state) : =0 for non-constitutive part. variable, typart=1 for constitutive particulate variable
   !&E       ws_free_opt,ws_free_para,ws_free_min,ws_free_max,ws_hind_opt,ws_hind_para : parameters defininf settling velocities
   !&E       water concentrations [cw(iv,k)] iv=1:nv_adv & k=1:kmax  (Unity : Masse/m3 of water)
   !&E       sediment concentrations [cs(iv,k)] iv=-1:nv_tot & k=ksdmin:ksdmax 
   !&E                                iv=-1 : temperature  - iv=0 : salinity 
   !&E                                iv=1,nvpc => Unity : kg/m3 of total sediment      
   !&E                                iv=nvpc+1,nvp => Unity : Masse/m3 of total sediment     
   !&E                                iv=nvp+1,nv_state => Unity : Masse/m3 of interstitial water     
   !&E       tocd : critical stress of deposition
   !&E       irkm_var_assoc : number of the constitutive particulate variable associated with a non constitutive particulate variable (both have same settling velocity)
   !&E                        If = 0, particulate variable with its own settling velocity (or dissolved variable)
   !&E
   !&E
   !&E===================================================================================================================
#include "coupleur_define_MUSTANG.h"
#ifdef key_MARS
#include "toolcpp.h"
#endif
#include "coupleur_dimhydro_MUSTANG.h"

#ifdef key_MARS
   USE comvars2d,      ONLY : ig,id,jb,jh,fwet,hm
   USE comvars3d,    ONLY : l_stateeq_lin
   USE parameters,     ONLY : epsilon,grav,valmanq
   !! variables du module SUBSTANCE utilisees et necessaires pour le module MUSTANG (pour MARS)
   USE comsubstance, ONLY : nv_adv,nv_tot,nvpc,nvp,igrav1,igrav2,isand1,isand2,imud1,imud2
   USE comsubstance, ONLY : ws_free_opt,ws_free_para,ws_free_min,ws_free_max,ws_hind_opt,ws_hind_para
   USE comsubstance, ONLY : tocd,irkm_var_assoc,l_subs2D,typart
   
!#else
   !! variables  SUBSTANCE known via 
#endif

   USE comMUSTANG 
   USE coupleur_MUSTANG 

   IMPLICIT NONE

   !! * Interface
   
   
   !! * Accessibility

   ! functions & routines of this module, called outside :
   ! PUBLIC functions                                ! none ?
   PUBLIC sed_MUSTANG_first,sed_MUSTANG_deposition

   PRIVATE

   !! * Shared or public module variables (SEDIM variables 1DV)

   
 CONTAINS
 
    !!==============================================================================
 
  SUBROUTINE sed_MUSTANG_first(ifirst,ilast,jfirst1,jfirst2,jlast,      &
               WATER_CONCENTRATION,SALINITY_MOD,TEMPERATURE_MOD ,       &
               THICK_LAYER_WAT,VERT_VISC,h0fond,                        &
#if ! defined key_MARS
               WATER_FLUX_INPUTS,SETTL_FLUX_w2s,EROS_FLUX_s2w,WAT_SETTL,&
               EROS_FLUX_TEMP_s2w,EROS_FLUX_SAL_s2w, &
               CORFLUX_SAND,CORFLUY_SAND, &
#endif
               h0,ssh,rhoref,roswat_bot,surf,saliref_lin,tetaref_lin,dt_true,t,u,v)

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_MUSTANG_first ***
   !&E
   !&E ** Purpose : update forcings (bottom stress, u* , settling velocities..)
   !&E               + remodelling (erosion, consolidation, diffusion, bioturbation
   !&E         
   !&E                 ATTENTION ON EST EN ZONE PARALLEL OMP
   !&E         
   !&E ** Description :
   !&E
   !&E  arguments IN : 
   !&E         loops  :ifirst,ilast,jfirst1,jfirst2,jlast
   !&E         bathy  :h0,surf,h0fond
   !&E         parametres ref  :rhoref
   !&E         time  :dt_true,t (DOUBLE PRECISION)
   !&E         hydro  :ssh,u,v,THICK_LAYER_WAT,VERT_VISC
   !&E         concentrations  : WATER_CONCENTRATION,SALINITY_MOD,TEMPERATURE_MOD
   !&E
   !&E  arguments OUT:
   !&E      pour modeles autres que MARS : (pour MARS : ils sont les memes dans MARS et dans MUSTANG)
   !&E      vaariables dont les dimensions et les noms peuvent etre differents dans le modele hydro/MUSTANG 
   !&E         flux d eau apporte dans chaque maille de la colonne d'eau : WATER_FLUX_INPUTS
   !&E         flux de depot pour chaque variable : SETTL_FLUX_w2s
   !&E         correction des flux horizontaux pour les sables traites en 2D :  CORFLUX_SAND,CORFLUY_SAND
   !&E
   !&E  en sortie en plus des arguments OUT ci dessus, on a remodelle le sediment
   !&E         nouvelles concentrations dans le sediment (cv_sed, poro)
   !&E         nombre de couches et epaisseurs des couches modifies (ksma, dzs)
   !&E         
   !&E        
   !&E ** Called by :  step before mouvement and transport equations
   !&E
   !&E ** Reference : 
   !&E
   !&E ** History :
    !&E       !  2015-12  (B.Thouvenin) reorganization of module SEDIMARS=MUSTANG
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used


#if defined key_MARS 
   USE sed_MUSTANG_MARS,    ONLY :  sed_skinstress,sed_gradvit,sed_obc_corflu,sed_exchange_s2w_MARS
#if defined key_parsub_newdt
   USE comsubstance,  ONLY   : ws3max
#endif
   USE comsubstance,  ONLY   : ws3 
   USE comfluxdo,    ONLY : phieau
#endif

   !! * Arguments
   !! jfirst1= ljminp1=ljmin for intern proc , or = jmin+1 for south limit proc
   !! jfirst2= ljminp2=ljmin for intern proc , or = jmin+2 for south limit proc
   INTEGER, INTENT(IN)             :: ifirst,ilast,jfirst1,jfirst2,jlast                           
   REAL(KIND=rsh),INTENT(IN)                                 :: h0fond,rhoref,saliref_lin,tetaref_lin 
   REAL(KIND=rlg),INTENT(IN)                                 :: dt_true,t  ! !  (dt_true=halfdt in MARS)
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_m2p2),INTENT(IN)   :: h0,ssh,surf,u,v                         
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_cw), INTENT(IN) :: WATER_CONCENTRATION         
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_TEMPSAL), INTENT(IN)   :: SALINITY_MOD,TEMPERATURE_MOD  
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_eps), INTENT(IN):: THICK_LAYER_WAT                
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_nz), INTENT(IN) :: VERT_VISC     
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY),INTENT(IN)        :: roswat_bot
#if ! defined key_MARS
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_ws), INTENT(INOUT) :: WAT_SETTL
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_waterflx), INTENT(OUT) :: WATER_FLUX_INPUTS 
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_flx_SW), INTENT(OUT) :: SETTL_FLUX_w2s,EROS_FLUX_s2w   
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_corflx), INTENT(OUT) :: CORFLUX_SAND, CORFLUY_SAND
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_flx_2D), INTENT(OUT) :: EROS_FLUX_TEMP_s2w,EROS_FLUX_SAL_s2w
#endif

   !! * Local declarations
   INTEGER        ::  i,j,k,iv,ksmin,ksmax,iappel
   REAL(KIND=rsh) ::  dtinv

 
   !!---------------------------------------------------------------------------
   !! * Executable part

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!! PRELIMINARY EVALUATIONS FOR COUPLING !!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! coupleur pour qq calcul propre a chaque modele ==> htot, alt_cw1, epn_bottom
    ! + pour modele autre que MARS : eps (THICK_LAYER_WAT), ws3 (WAT_SETTL)
    ! + calcul des concentrations au fond + epaisseur de la couche du fond
    iappel=1
    CALL coupl_conv2MUSTANG(ifirst,ilast,jfirst1,jlast,iappel,h0,ssh,   &
#if ! defined key_MARS
                 THICK_LAYER_WAT,VERT_VISC,WAT_SETTL,SETTL_FLUXSUM_w2s, &
#endif
                       WATER_CONCENTRATION,SALINITY_MOD,TEMPERATURE_MOD )
            

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!! FORCING :    tension de fond (bottom shear stress) (tenfon)!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef key_MARS
         ! depend du modele hydro et echanges ==> dans MUSTANG_MARS.F90
         ! ==> tenfon (comMUSTANG)
    CALL sed_skinstress
#else
         ! calcul a programmer pour autre modele ?
#endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!! FORCING :     u*=ustarbot            !!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!$OMP DO SCHEDULE(RUNTIME)
      DO j=jfirst1,jlast
#ifdef key_MARS
      DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
#else
      DO i=ifirst,ilast
#endif

          IF(htot(i,j).GT.h0fond)  THEN
            IF(tenfon(i,j) < 0.) THEN
              ustarbot(i,j)=0.0_rsh
            ELSE
              ustarbot(i,j)=(tenfon(i,j)/rhoref)**0.5_rsh
            ENDIF
          ENDIF
#ifdef key_MARS
!          DO iv=1,nv_adv
!             IF(cw_bottom_MUSTANG(iv,i,j) < -1.0e10_rsh)write(*,*)'cwb<0',t,i,j,iv,cw_bottom_MUSTANG(iv,i,j)
!          ENDDO
        ENDIF
#endif
      ENDDO
      ENDDO
!$OMP END DO


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!! FORCING : G dU/dz (taux de cisaillement : shear rate)   !!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     ! calcul de gradvit et ustarbot - peut dependre du code..
      ! en sortie : gradvit pas en argument car dans les common comMUSTANG
#ifdef key_MARS
    CALL sed_gradvit(ifirst,ilast,jfirst1,jlast,h0fond,rhoref)
#else
         ! calcul a programmer pour autre modele ?
#endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!! vitesse de chute : settling rate     !!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#if defined key_parsub_newdt && defined key_MARS
!$OMP SINGLE
   ! rajoute pour calcul de ndt_part en fonction du temps
      ws3max(:)=0.0_rsh
!$OMP END SINGLE
#endif
    CALL sed_MUSTANG_settlveloc(ifirst,ilast,jfirst1,jlast,      &
                           h0fond,rhoref,cv_wat,    &
                           THICK_LAYER_WAT_MUSTANG,VERT_VISC_MUSTANG, ws3)
                           
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!! tendance au depot : tendency to deposit  flx_w2s (m.s-1) !!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         
    CALL sed_MUSTANG_depflx(ifirst,ilast,jfirst2,jlast,ws3)
    
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!! Correction for sand horizontal transport  corflux,corfluy ,flx_w2s(sand)  !!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    CALL sed_MUSTANG_sandconcextrap(ifirst,ilast,jfirst2,jlast, &
                           h0fond,ws3)

   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!! traitement corflux, corfluy aux coins et echange entre processeurs MPI  !!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#ifdef key_MARS
    IF (isand2.GE.1) THEN
              ! routine MARS dependent dans MUSTANG_MARS.F90  
         CALL sed_obc_corflu
    ENDIF
#else
      ! a programmer pour autre modele ??
#endif     

        
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!! update ws3max (MARS)  !!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef key_MARS
#ifdef key_parsub_newdt
!$OMP DO REDUCTION (MAX:ws3max) private(i,j,iv,k) 
     DO j=jfirst1,jlast
        DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
          DO iv=1,nvp
            DO k=1,kmax
              IF (ws3(k,iv,i,j)>ws3max(iv))ws3max(iv)=ws3(k,iv,i,j)
            ENDDO
          ENDDO
        ENDDO
     ENDDO
!$OMP END DO
#endif
#endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! CONSOLIDATION & DIFFUSION  every half time step !!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


   dtinv=1.0_rsh/REAL(dt_true,rsh)
   
   IF(l_dyn_insed)CALL sed_MUSTANG_consol_diff_bioturb(ifirst,ilast,jfirst1,jlast,  &
               surf,h0fond,roswat_bot,rhoref,saliref_lin,tetaref_lin,t,dt_true)
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! EROSION every half time step !!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! initialisation EROSION LATERALE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!$OMP SINGLE
   IF(coef_erolat .NE. 0.0_rsh) THEN
      flx_s2w_corip1(:,:,:)=0.0_rsh
      flx_s2w_corim1(:,:,:)=0.0_rsh
      flx_s2w_corjp1(:,:,:)=0.0_rsh
      flx_s2w_corjm1(:,:,:)=0.0_rsh
   ENDIF  ! fin if coef_erolat
!$OMP END SINGLE
   
!!!!!!!!!!!!!!!!!!!!
!!! EROSION    !!!!!
!!!!!!!!!!!!!!!!!!!!   
   CALL sed_MUSTANG_erosion(ifirst,ilast,jfirst1,jlast,h0fond,dtinv,u,v, &
                             surf,dt_true)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Fluxes correction if EROSION LATERALE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   IF(coef_erolat .NE. 0.0_rsh) THEN
     ! lateral erosion of dry cell
     ! echange MPI
#if defined key_MARS && defined key_MPI_2D
      CALL sed_exchange_s2w_MARS
#else
      !! pour autre modele
#endif

      ! correction : neighboring cells of eroded laterally  dry cell receive one fraction of eroded sediment 
!$OMP DO SCHEDULE(RUNTIME) private(i,j,iv) 
      DO j=jfirst1,jlast
#ifdef key_MARS
        DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
           IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
#else
        DO i=ifirst,ilast
#endif

!   attention  if may be different for complex grid
             DO iv=-1,nv_adv
               flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+ dtinv*(                            &
                       +flx_s2w_corip1(iv,i-1,j)+flx_s2w_corim1(iv,i+1,j)   &
                       +flx_s2w_corjp1(iv,i,j-1)+flx_s2w_corjm1(iv,i,j+1))
             ENDDO
#ifdef key_MARS 
          ENDIF  ! fin boucle sur jb,jh 
#endif
        END DO
      END DO
!$OMP END DO
                          
   ENDIF  ! fin if coef_erolat

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! conversion des flux de depot, des vitesses de chute et des corflux pour le modele hydro        !!!!!
! + conversion des flux d erosion et des corflux pour le modele hydro                             !!!!
! + flux d eau a l interface sediment/eau resultant du depot au demi dt precedent et de l erosion  !!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#if ! defined key_MARS
    CALL coupl_MUSTANG2hydro(ifirst,ilast,jfirst1,jlast,SETTL_FLUX_w2s,  &
                            CORFLUX_SAND,CORFLUY_SAND,WAT_SETTL,
                            EROS_FLUX_s2w,EROS_FLUX_TEMP_s2w,EROS_FLUX_SAL_s2w)
#endif

#if ! defined key_nofluxwat_IWS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    prise en compte des flux d eau du sediment vers eau dans la couche du fond du a l erosion et a la consolidation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !! WATER_FLUX_INPUT_BOTCELL = WATER_FLUX_INPUTS (k,i,j) in bottom cell = phieau(1,:,:)
!$OMP SINGLE
    WATER_FLUX_INPUT_BOTCELL=WATER_FLUX_INPUT_BOTCELL+phieau_s2w(:,:)/dt_true
    phieau_s2w(:,:)=0.0_rlg
!$OMP END SINGLE
#endif

  END SUBROUTINE sed_MUSTANG_first
  

  !!==============================================================================
  
  SUBROUTINE sed_MUSTANG_deposition(ifirst,ilast,jfirst1,jfirts2,jlast,  &
                        xeloc,surf,h0,dx,dy,h0fond,WATER_CONCENTRATION, &
#if ! defined key_MARS
                        SETTL_FLUXSUM_w2s,                              &
#endif
                        SALINITY_MOD,TEMPERATURE_MOD)
              
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_MUSTANG_deposition  ***
   !&E
   !&E ** Purpose : settling from water to sediment, sliding processes 
   !&E                    and sediment sedimentation 
   !&E
   !&E                 ATTENTION ON EST EN ZONE PARALLEL OMP
   !&E
   !&E ** Description : 
   !&E  arguments IN : 
   !&E         loops  :ifirst,ilast,jfirst1,jfirst2,jlast
   !&E         bathy  :h0,surf,h0fond,dx,dy
   !&E         hydro  :xeloc (=ssh)
   !&E         concentrations  : WATER_CONCENTRATION,SALINITY_MOD,TEMPERATURE_MOD
   !&E      pour modeles autres que MARS : (pour MARS : ils sont les memes dans MARS et dans MUSTANG)
   !&E         flux de depot effectif apres transport pour chaque variable : SETTL_FLUXSUM_w2s
   !&E
   !&E  arguments OUT:
   !&E      pour modeles autres que MARS : (pour MARS : ils sont les memes dans MARS et dans MUSTANG)
   !&E      variables dont les dimensions et les noms peuvent etre differents dans le modele hydro/MUSTANG 
   !&E       SETTL_FLUXSUM_w2s : flux de depot pour chaque substance effective apres transport
   !&E
   !&E  en sortie en plus des arguments OUT ci dessus, on a remodelle le sediment
   !&E         nouvelles concentrations dans le sediment (cv_sed, poro)
   !&E         nombre de couches et epaisseurs des couches modifies (ksma, dzs)
   !&E         
   !&E
   !&E ** Called by :  step after transport equations
   !&E
   !&E ** Reference : 
   !&E
   !&E ** History :
    !&E       !  2015-12  (B.Thouvenin) reorganization of module SEDIMARS
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
#ifdef key_MARS
   USE sed_MUSTANG_MARS,    ONLY :  sed_exchange_w2s_MARS
#if defined key_siggen || defined key_gencoord
  USE toolgeom,     ONLY : f_hzu
#endif
#endif

   !! * Arguments
   INTEGER, INTENT(IN)  :: ifirst,ilast,jfirst1,jfirts2,jlast
   REAL(KIND=rsh),INTENT(IN) :: h0fond
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_m2p2),INTENT(IN) :: xeloc,surf,h0,dx,dy   ! PROC_IN_ARRAY_m2p2
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_TEMPSAL),INTENT(IN) :: SALINITY_MOD,TEMPERATURE_MOD  ! PROC_IN_ARRAY_TEMPSAL
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_cw), INTENT(IN)  :: WATER_CONCENTRATION  !  DIM_MODELE_HYDRO_cw
#if ! defined key_MARS
   REAL(KIND=rsh),DIMENSION(DIM_MODELE_HYDRO_flx_SW), INTENT(IN)  :: SETTL_FLUXSUM_w2s  !! DIM_MODELE_HYDRO_flx_SW
#endif

   !! * Local declarations
   INTEGER        :: iappel
   !! * Executable part 
  
   !deposit if strong bottom slope 
!$OMP SINGLE  
     flx_w2s_corim1(:,:,:)=0.0_rsh   
     flx_w2s_corip1(:,:,:)=0.0_rsh
     flx_w2s_corjm1(:,:,:)=0.0_rsh
     flx_w2s_corjp1(:,:,:)=0.0_rsh
     flx_w2s_corin(:,:,:)=0.0_rsh 
!$OMP END SINGLE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! coupleur pour conversion des nouvelles valeurs de concentrations et 
!!!! hauteurs d eau et flux de depot effectif
!!!! + calcul des concentrations au fond + epaisseur de la couche du fond
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    iappel=2
    CALL coupl_conv2MUSTANG(ifirst,ilast,jfirst1,jlast,iappel,h0,xeloc, &
#if ! defined key_MARS
                     THICK_LAYER_WAT,VERT_VISC,WAT_SETTL,SETTL_FLUXSUM_w2s, &
#endif
                       WATER_CONCENTRATION,SALINITY_MOD,TEMPERATURE_MOD )        
        
!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! initialisation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#if ! defined key_nodiss_insed
!$OMP SINGLE  
   flx_w2s(nvp+1:nv_adv,:,:)=0.0_rsh
   flx_w2s(-1:0,:,:)=0.0_rsh
!$OMP END SINGLE
#endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! glissement des depots si pente forte (slidepo)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   IF(slopefac .NE. 0.0_rsh) THEN
     CALL sed_MUSTANG_slipdepo(ifirst,ilast,jfirst1,jlast,       &
             dx,dy,surf,h0fond,h0)  

! echange MPI
#if defined key_MARS && defined key_MPI_2D 
         CALL sed_exchange_w2s_MARS
             
#else
 !! pour autre modele
#endif
   ENDIF 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! effective deposit after variables transport and settling => sediment sedimentation 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   call sed_MUSTANG_effdep(ifirst,ilast,jfirst1,jlast,surf,h0fond)
                   
  ! rajouter conversion pour les flux  deau et les flux a l interface pour les dissoutes
  

   PRINT_DBG*, 'FIN sed_MUSTANG_deposition'   


  END SUBROUTINE sed_MUSTANG_deposition
  
!=========================================================================== 


  SUBROUTINE sed_MUSTANG_settlveloc(ifirst,ilast,jfirst,jlast,    &
                     h0fond,rhoref,cw,eps,nz,ws3)

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_MUSTANG_settlveloc  ***
   !&E
   !&E ** Purpose : settling velocity computation
   !&E
   !&E ** Description : use arguments and common variable 
   !&E  arguments IN : 
   !&E         h0fond, : RESIDUAL_THICKNESS_WAT (m) 
   !&E         rhoref : RHOREF
   !&E         cw : WATER_CONCENTRATION_converti for MUSTANG
   !&E         eps : THICK_LAYER_WAT_MUSTANG
   !&E         nz: VERT_VISC_MUSTANG
   !&E         ws3: WAT_SETTL_MUSTANG
   !&E  arguments OUT:
   !&E         ws3: WAT_SETTL_MUSTANG
   !&E
   !&E  need to be know by hydrodynamic code:
   !&E         kmax  : connu via coupleur_dimhydro_MUSTANG.h
   !&E         grav :
   !&E  need to be know by code treated substance si pas MARS ==> coupleur_MUSTANG.F90
   !&E         imud1, nvpc, nvp, nv_adv, isand1,isand2
   !&E         f_ws(iv) (if key_flocmod)
   !&E         ws_free_opt,ws_free_para,ws_free_min,ws_free_max,ws_hind_opt,ws_hind_para :    
   !&E     
   !&E  use module MUSTANG variables  :
   !&E         ros(iv)
   !&E         ws_sand(iv)
   !&E     
   !&E     
   !&E     
   !&E ** Called by :  sed_MUSTANG_first
   !&E
   !&E ** Reference : Settling velocity computed after Van Leussen (1994)
   !&E                Ws=Ws*(1+aG)/(a+bG**2) with a=0.3 and b=0.09 (Ems estuary)
   !&E                G=sqrt(eps(k,i,j)/nz(k,i,j))
   !&E
   !&E ** History :
   !&E       !  2008-09  (Pierre Le Hir) extracted from init.F90
   !&E       !  2008-10  (F. Cayocca) inserted into module_sedimento.F90
   !&E       !  2014-05  (R. Verney) modifying available settling velocity formulations
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
#if defined key_MARS 
#if defined key_parsub_newdt
   USE comsubstance,  ONLY   : ws3max 
#endif
   USE parameters,   ONLY : l_testcase
#endif

   !! * Arguments
   INTEGER, INTENT(IN)                               :: ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),INTENT(IN)                         :: h0fond,rhoref
   REAL(KIND=rsh),DIMENSION(nv_adv,kmax,PROC_IN_ARRAY_m1p2), INTENT(IN)  :: cw     
   REAL(KIND=rsh),DIMENSION(0:kmax,PROC_IN_ARRAY_0p1), INTENT(IN)        :: nz      
   REAL(KIND=rsh),DIMENSION(kmax,PROC_IN_ARRAY), INTENT(IN)              :: eps     
   REAL(KIND=rsh),DIMENSION(kmax,nv_adv,PROC_IN_ARRAY), INTENT(INOUT)    :: ws3     

   
   !! * Local declarations
   INTEGER                            :: iv,k,ivpc,i,j
   REAL(KIND=rsh)                     :: cmes,phi,phiv,De
   REAL(KIND=rsh),PARAMETER           :: nuw=0.00000102_rsh
!MD 10/12/2018
   INTEGER, PARAMETER                 :: sal_lim=5
!END MD

   !!---------------------------------------------------------------------------
   !! * Executable part
#if defined key_MARS && defined key_parsub_newdt
!$OMP SINGLE
   ! rajoute pour calcul de ndt_part en fonction du temps
      ws3max(:)=0.0_rsh
!$OMP END SINGLE
#endif

!$OMP DO SCHEDULE(RUNTIME) private(j,i,k,ivpc,iv,phi &
!$OMP& ,cmes,De,phiv)  
      DO j=jfirst,jlast
#ifdef key_MARS
      DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
#else
      ! autres bornes pour autre modele ?
      DO i=ifirst,ilast
#endif
           
         IF(htot(i,j) > h0fond)  THEN
           DO k=1,kmax
              cmes=0.0_rsh
              DO ivpc=imud1,nvpc
                cmes=cmes+cw(ivpc,k,i,j)
              ENDDO
              cmes=MAX(0.0_rsh,cmes)

! first calculating sand settling velocity
              DO iv=isand1,isand2
                ws3(k,iv,i,j)=ws_sand(iv)
              ENDDO
              
! next mud settling velocity 
#ifdef key_flocmod
              DO iv=imud1,nvp         
                 ws3(k,iv,i,j)=f_ws(iv)
#ifdef key_MARS
                 if (l_testcase) ws3(k,iv,i,j)=0.0_rsh    
#endif  
              ENDDO
      
#else
      
              DO iv=imud1,nvp
! free settling velocity - flocculation
    
                IF(ws_free_opt(iv) == 0) THEN ! constant settling velocity
                  ws3(k,iv,i,j)=ws_free_max(iv)
                ELSEIF (ws_free_opt(iv) == 1) THEN ! Van Leussen 1994
                  ws3(k,iv,i,j)=ws_free_para(1,iv)*cmes**ws_free_para(2,iv) &
                              *(1._rsh+ws_free_para(3,iv)*gradvit(k,i,j))/  &
                               (1._rsh+ws_free_para(4,iv)*gradvit(k,i,j)**2._rsh)
                ELSEIF (ws_free_opt(iv) == 2) THEN ! Winterwerp 1999
                  De=ws_free_para(1,iv)+ws_free_para(2,iv)*cmes/(ws_free_para(3,iv)*sqrt(gradvit(k,i,j)))
                  IF (De.GT.sqrt(nuw/gradvit(k,i,j))) THEN 
                    De=sqrt(nuw/gradvit(k,i,j)) ! in case of large C/low G limit floc size to kolmogorov microscale
                  ENDIF
                  ws3(k,iv,i,j)=(ros(iv)-rhoref)*grav/(18._rsh*rhoref*nuw)  &
                                *ws_free_para(1,iv)**(3._rsh-ws_free_para(4,iv))  &
                                *De**(ws_free_para(4,iv)-1._rsh)
                
                ELSEIF (ws_free_opt(iv) == 3) THEN ! Wolanski et al., 1989
                  ws3(k,iv,i,j)=ws_free_para(1,iv)*cmes**ws_free_para(2,iv)

                ENDIF

! MD 10/12/2018
                IF(cv_sed(0,k,i,j).LT.sal_lim) THEN  !Si dans l'eau douce
                  ws3(k,iv,i,j)=max(ws_free_min(iv),cv_sed(0,k,i,j)*ws3(k,iv,i,j)/sal_lim)
                ENDIF
! END MD
! Hindered settling
! if ws_hind_opt.EQ.0 : no hindered settling...ws3 unchanged
                IF (ws_hind_opt(iv) == 1) THEN ! Scott, 1984
                  phi=MIN(1.0_rsh,cmes/ws_hind_para(1,iv))
                  ws3(k,iv,i,j)=ws3(k,iv,i,j)*(1._rsh-phi)**ws_hind_para(2,iv)
                ELSEIF (ws_hind_opt(iv) == 2) THEN ! Winterwerp, 2002 - ATTENTION : ros(iv) must be the same for all MUDS variables
                  phi=cmes/ros(iv)
                  IF (ws_free_opt(iv) == 2) THEN
                    phi=cmes/ros(iv)
                    phiv=phi*(De/ws_free_para(1,iv))**(3._rsh-ws_free_para(4,iv))
                  ELSE
                    phiv=cmes/ws_hind_para(1,iv)
                  ENDIF

                  ws3(k,iv,i,j)=ws3(k,iv,i,j)*(1._rsh-phiv)**ws_hind_para(2,iv)*(1._rsh-phi)  &
                                /(1._rsh+2.5_rsh*phiv)
                ELSEIF (ws_hind_opt(iv) == 3) THEN ! wolanski et al., 1989
                  IF (ws_free_opt(iv) == 3) THEN
                    ws3(k,iv,i,j)=ws3(k,iv,i,j)/(cmes**2._rsh+ws_hind_para(1,iv)**2._rsh)**ws_hind_para(2,iv)
                  ELSE
                    
                  ENDIF
                ENDIF

! limiting ws with min/max values...
! necessary if ndt_part not updated during the simulation, ndt_part calculated in subreaddat from a maximum settling velocity
                ws3(k,iv,i,j)=max(ws_free_min(iv),min(ws_free_max(iv),ws3(k,iv,i,j)))

              ENDDO
              

            ! rverney : attention update parsub_newdt with flocculation! 

#endif
              DO iv=nvp+1,nv_adv
                ws3(k,iv,i,j)=0.0_rsh
              ENDDO
              
           ENDDO
         ELSE
           ws3(:,:,i,j)=0.0_rsh
         ENDIF

#ifdef key_MARS
        ENDIF
#endif
       ENDDO
       ENDDO
!$OMP END DO


  END SUBROUTINE sed_MUSTANG_settlveloc      

   !!==============================================================================

  SUBROUTINE sed_MUSTANG_depflx(ifirst,ilast,jfirst,jlast,ws3)
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_MUSTANG_depflx  ***
   !&E
   !&E ** Purpose : computes deposition flux trends  before advection computations
   !&E
   !&E ** Description : use arguments and common variable 
   !&E  arguments IN : 
   !&E         ws3 : settling velocity converted for MUSTANG
   !&E          
   !&E  variables OUT:
   !&E         flx_w2s : tendance au depôt (m/s) 
   !&E         
   !&E  need to be know by hydrodynamic code:
   !&E         kmax  (known from coupleur_dimhydro_MUSTANG.h)
   !&E
   !&E  need to be know by code treated substance (si pas MARS known by coupleur_MUSTANG.F90
   !&E         igrav2,isand2, nvpc, nvp, nv_adv : 
   !&E         irkm_var_assoc (ivp)
   !&E         tocd(iv)
   !&E      
   !&E
   !&E ** Called by :  sed_MUSTANG_first
   !&E 
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-03  (F. Cayocca) extracted from marsedim.F90 (from F. Dufois)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used


   !! * Arguments
   INTEGER, INTENT(IN)  :: ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),DIMENSION(kmax,nv_adv,PROC_IN_ARRAY), INTENT(IN)     :: ws3 
   
   !! * Local declarations
   INTEGER                          :: ivp,i,j
   REAL(KIND=rlg)                   :: tocdpe,depo
   REAL(KIND=rlg),PARAMETER         :: epsdepot=1.e-14

   !!---------------------------------------------------------------------------
   !! * Executable part

!$OMP DO SCHEDULE(RUNTIME) PRIVATE(i,j,ivp,tocdpe,depo)
      DO j=jfirst,jlast
#ifdef key_MARS
      DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
#else
      ! autres bornes pour autre modele ?
      DO i=ifirst,ilast
#endif

       ! Separation des boucles car pas forcement la meme unite pour les non constitutives
       ! pas de calcul pour les graviers actuellement
       DO ivp=igrav2+1,nvpc
         IF(ivp.LE.isand2) THEN 
            ! Sanford et Halka (1993), Winterwerp (2004), Dufois (2008)
            flx_w2s(ivp,i,j)=ws3(1,ivp,i,j)*fwet(i,j)
         ELSE
            tocdpe=tocd(ivp)+0.00001_rsh
            flx_w2s(ivp,i,j)=ws3(1,ivp,i,j)*fwet(i,j)                &
                          *MAX(0.0_rsh,1.0_rsh-(tenfon(i,j)/tocdpe)) &  
                          *tocd(ivp)/tocdpe
         ENDIF
         depo=flx_w2s(ivp,i,j)*cw_bottom_MUSTANG(ivp,i,j)
         IF(depo.LE.epsdepot)flx_w2s(ivp,i,j)=0.0_rsh
       ENDDO
#ifdef key_sedim_partial
       DO ivp=nvpc+1,nvp
          flx_w2s(ivp,i,j)=0.0_rsh
       ENDDO
#else
       DO ivp=nvpc+1,nvp
         tocdpe=tocd(ivp)+0.00001_rsh
         flx_w2s(ivp,i,j)=ws3(1,ivp,i,j)*fwet(i,j)                   &
                          *MAX(0.0_rsh,1.0_rsh-(tenfon(i,j)/tocdpe)) &  
                          *tocd(ivp)/tocdpe
         depo=flx_w2s(ivp,i,j)*cw_bottom_MUSTANG(ivp,i,j)
         IF(irkm_var_assoc(ivp) > 0 ) THEN
          IF (flx_w2s(irkm_var_assoc(ivp),i,j)==0.0_rsh) flx_w2s(ivp,i,j)=0.0_rsh
         ELSE
          IF(depo.LE.epsdepot)flx_w2s(ivp,i,j)=0.0_rsh
         END IF
       ENDDO
#endif

#ifdef key_MARS
        ENDIF
#endif
       ENDDO
       ENDDO
!$OMP END DO
         
  END SUBROUTINE sed_MUSTANG_depflx

!!==============================================================================

  SUBROUTINE sed_MUSTANG_sandconcextrap(ifirst,ilast,jfirst,jlast,h0fond,ws3)

   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_MUSTANG_sandconcextrap  ***
   !&E
   !&E ** Purpose : interpolates sand concentration in the first layer
   !&E              assuming a Rouse profile 
   !&E
   !&E ** Description : use arguments and common variable 
   !&E  arguments IN : 
   !&E         sig (COORD_SIGMA)
   !&E         h0fond (RESIDUAL_THICKNESS_WAT)
   !&E         ws3 :
   !&E  variables OUT : 
   !&E     flx_w2s : tendances aux depots corriges pour sables 
   !&E     corflux, corfluy : corrections des flux horizontaux pour les sables
   !&E     
   !&E     
   !&E  need to be know by hydrodynamic code:
   !&E         kmax (known from coupleur_dimhydro_MUSTANG.h)
   !&E         alt_cw1 , htot: evaluated in coupleur_conv2MUSTANG
   !&E         
   !&E  need to be know by code treated substance:
   !&E         isand1,isand2,nv_adv 
   !&E     
   !&E  use module MUSTANG variables  :
   !&E         aref_sand
   !&E         ustarbot (evaluated in sed_gravit)
   !&E
   !&E ** Called by :  sed_MUSTANG_first
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (P. Le Hir) extracted from init.F90
   !&E       !  2008-10  (F. Cayocca) inserted into module_sedimento.F90
   !&E       !  2009-02 (P. Le Hir)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   INTEGER, INTENT(IN)                        :: ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),INTENT(IN)                  :: h0fond
   REAL(KIND=rsh),DIMENSION(kmax,nv_adv,PROC_IN_ARRAY), INTENT(IN)  :: ws3     

   !! * Local declarations
   INTEGER                              :: izz,ivp,i,j
   REAL(KIND=rsh)                       :: altc1,extrap,rouse,som,dzcche,hzed,     &
                                           einstein,alogaltc1sz0,ashmapwrouse,rouse1
   INTEGER,PARAMETER                    :: nbcche=20
   REAL(KIND=rsh),DIMENSION(1:nbcche+1) :: hzi
#ifdef key_sand2D
   REAL(KIND=rsh),DIMENSION(1:nbcche+1) :: hzisdbot
#endif


   !!---------------------------------------------------------------------------
   !! * Executable part


!$OMP SINGLE
    corflux(:,:,:)=1.0_rsh
    corfluy(:,:,:)=1.0_rsh
!$OMP END SINGLE

!$OMP DO SCHEDULE(RUNTIME) PRIVATE(i,j,altc1,extrap,ivp,alogaltc1sz0  &
!$OMP ,hzi,izz        &
#ifdef key_sand2D
!$OMP ,hzisdbot       &
#endif
!$OMP ,rouse,rouse1,som,einstein,hzed,dzcche,ashmapwrouse)
      DO j=jfirst,jlast
#ifdef key_MARS
      DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
        IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
#else
      ! autres bornes pour autre modele ?
      DO i=ifirst,ilast
#endif

          altc1=alt_cw1(i,j)
          IF(altc1.LE.aref_sand .OR. htot(i,j).LE.h0fond)THEN
            extrap=1.0_rsh
            ! in this case, flx_w2s(ivp,i,j) is not changed
            DO ivp=isand1,isand2
              corflux(ivp,i,j)=1.0_rsh
              corfluy(ivp,i,j)=1.0_rsh
            ENDDO
          ELSE
            alogaltc1sz0=LOG(altc1/z0sed(i,j))
            hzi(1)=epn_bottom_MUSTANG(i,j)
            DO izz=1,nbcche
              hzi(izz+1)=hzi(izz)-(hzi(1)-aref_sand)/2.0_rsh**izz
            ENDDO     
            hzi(nbcche+1)=aref_sand
#ifdef key_sand2D
            hzisdbot(1)=htot(i,j)
            DO izz=1,nbcche
              hzisdbot(izz+1)=hzisdbot(izz)-(hzisdbot(1)-aref_sand)/2.0_rsh**izz
            ENDDO     
            hzisdbot(nbcche+1)=aref_sand
#endif
            DO ivp=isand1,isand2
            ! Rouse number depends on settling velocity, then on the sand variable
            ! modif Julie 2012 to taking into account flat bottom
              IF(ws3(1,ivp,i,j)>1.34_rsh*ustarbot(i,j)) THEN
                 rouse=1.2_rsh
              ELSE
                rouse1=ws3(1,ivp,i,j)/ustarbot(i,j)
                IF(rouse1>=0.75_rsh) THEN
                   rouse=0.353_rsh*rouse1+0.727
                ELSEIF(rouse1>0.1_rsh .AND. rouse1<0.75_rsh) THEN
                   rouse=rouse1/(0.4_rsh+0.8_rsh*rouse1*rouse1)   ! rouse=rouse1/(0.4*beta) avec beta=1+2*rouse1*rouse1
                ELSE
                  rouse=rouse1/0.4_rsh
                ENDIF
              ENDIF
            ! integrale concentration et concentration*vitesse log
            ! dans premiere couche
              som=0.0_rsh
              einstein=0.0_rsh
              DO izz=1,nbcche
#ifdef key_sand2D
               IF(l_subs2D(ivp)) THEN
                hzed=0.5_rsh*(hzisdbot(izz)+hzisdbot(izz+1))
                dzcche=hzisdbot(izz)-hzisdbot(izz+1) 
               ELSE
#endif
                hzed=0.5_rsh*(hzi(izz)+hzi(izz+1))
                dzcche=hzi(izz)-hzi(izz+1) 
#ifdef key_sand2D
               ENDIF
#endif
                som=som+dzcche*((htot(i,j)-hzed)/hzed)**rouse
                einstein=einstein+dzcche*((htot(i,j)-hzed)/hzed)**rouse      &  
                                  *LOG(hzed/z0sed(i,j))
              ENDDO
#ifdef key_sand2D
              IF(l_subs2D(ivp)) THEN
                rouse2D(ivp,i,j)=rouse
                sum_tmp(ivp,i,j)=som
              ENDIF
#endif 
!!
              ashmapwrouse=(aref_sand/(htot(i,j)-aref_sand))**rouse
              som=som*ashmapwrouse
              !IF (som==0) PRINT*, 'hzi(1),ashmapwrouse rouse',hzi(1),ashmapwrouse,rouse

              einstein=einstein*ashmapwrouse/som
              ! the computation of einstein and USE of corflux should be checked
              !??? modif P.LeHir              if(hex(i,j).le.hm)then
              !??? modif P.LeHir              if(hey(i,j).le.hm)then

#ifdef key_sand2D
              IF(l_subs2D(ivp)) THEN
                extrap=(hzisdbot(1)-aref_sand)/som 
                corflux(ivp,i,j)=einstein/LOG(htot(i,j)/(2.718_rsh*z0sed(i,j)))
                corfluy(ivp,i,j)=einstein/LOG(htot(i,j)/(2.718_rsh*z0sed(i,j)))
              ELSE
#endif
                extrap=(hzi(1)-aref_sand)/som
#ifdef key_MARS
                IF(htot(i,j).LE.hm)THEN
                  corflux(ivp,i,j)=einstein/LOG(htot(i,j)/(2.718_rsh*z0sed(i,j)))
                  corfluy(ivp,i,j)=einstein/LOG(htot(i,j)/(2.718_rsh*z0sed(i,j)))
                ELSE
#endif
                  corflux(ivp,i,j)=einstein/alogaltc1sz0
                  corfluy(ivp,i,j)=einstein/alogaltc1sz0
#ifdef key_MARS
                ENDIF
#endif
                
#ifdef key_sand2D
              ENDIF
#endif 
              flx_w2s(ivp,i,j)=flx_w2s(ivp,i,j)*extrap              
            ENDDO              
          ENDIF

#ifdef key_MARS
        ENDIF
#endif
       ENDDO
       ENDDO
!$OMP END DO 
         
#ifdef key_MARS
    PRINT_DBG*, 'FIN SED_SANDCONCEXTRAP'
#endif

  END SUBROUTINE sed_MUSTANG_sandconcextrap
      
   !!==============================================================================
  SUBROUTINE sed_MUSTANG_erosion(ifirst,ilast,jfirst,jlast,h0fond,dtinv,u,v,surf,dt_true) 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_MUSTANG_erosion  ***
   !&E
   !&E ** Purpose : computes erosion fluxes
   !&E
   !&E ** Description : 
   !&E       arguments IN :
   !&E          h0fond : RESIDUAL_THICKNESS_WAT
   !&E          dtinv, dt_true : 1/dt  et dt
   !&E          surf : surface des mailles (SURF)
   !&E          u,v :   BAROTROP_ZONAL_VELOCITY_U, BAROTROP_MERIDIONAL_VELOCITY_V
   !&E      
   !&E     variables IN (used) : 
   !&E         
   !&E        
   !&E        
   !&E        
   !&E     variables OUT : 
   !&E          flx_s2w : flux sediment to water
   !&E          phieau_s2w : flux d eau sediment vers eau
   !&E          remaniement deu sediment (cv_sed, poro, dzs, ksma)
   !&E
   !&E
   !&E ** Called by :  stepu, stepv
   !&E 
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (P. Le Hir) extracted from subsedim.F90 (from F. Dufois)
   !&E       !  2008-10  (F. Cayocca) included into module_sedimento
   !&E       !  2008-02  (P. Le Hir) replaced by new "Marennes-Oleron" sand/mud sediment model
   !&E       !  2011-04  (F. Cayocca) introduction of the fraction of wet cell (fwetp)
   !&E       !  2015-07  (B.Thouvenin  P. Le Hir) : correction lateral erosion of dry cell  
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   USE comMUSTANG

   INTEGER, INTENT(IN)                        :: ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),INTENT(IN)                  :: h0fond,dtinv
   REAL(KIND=rlg),INTENT(IN)                  :: dt_true  ! =halfdt in MARS
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_m2p2),INTENT(IN)   :: surf,u,v 

   !! * Local declarations
   INTEGER        ::  i,j,k,iv,ksmax,ksup,ksmaxa
   REAL(KIND=rsh) ::  dt1,toce,cvolgrv,csanmud,                          &
                      erodab,ero,erosi,ddzs,dzsa,dzsisu,dzsam1,dflx_s2w,       &
                      cvolp,volerod,profsed,poroa,poroam1,dflusve,      &
                      xeros,sed_tocr_mixsed,sed_eros_flx,excespowr,dzpoi,      &
                      heauw,heaue,heaun,heaus,heau_milieu,eroe,erow,eros,eron

   !!---------------------------------------------------------------------------
   !! * Executable part

!$OMP DO SCHEDULE(RUNTIME) PRIVATE(i,j,k,iv,ksmax,ksup,ksmaxa &
!$OMP ,dt1,toce,cvolgrv,csanmud,erodab,ero,erosi,ddzs,dzsa,dzsisu,dzsam1,dflx_s2w       &
!$OMP ,cvolp,volerod,profsed,poroa,poroam1,dflusve, xeros,sed_tocr_mixsed   &
!$OMP ,sed_eros_flx,excespowr,dzpoi,heauw,heaue,heaun,heaus,heau_milieu,eroe,erow,eros,eron)
!!!$OMP SINGLE
      DO j=jfirst,jlast
#ifdef key_MARS
        DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
         IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
#else
        DO i=ifirst,ilast
#endif

          flx_s2w(-1,i,j)=0.0_rsh
          flx_s2w( 0,i,j)=0.0_rsh
          DO iv=1,nv_adv
            flx_s2w(iv,i,j)=0.0_rsh
          ENDDO
        
          ksmax=ksma(i,j)
          ero=0.0_rsh

          IF(ksmax.GE.ksmi(i,j))THEN

            dt1=REAL(dt_true,rsh)
            k=ksmax
                       
        2   CONTINUE

            CALL sed_MUSTANG_comp_tocr_mixsed(k,i,j,xeros,excespowr,sed_tocr_mixsed)
            toce=MAX(0.00005_rsh,sed_tocr_mixsed)
            cvolgrv=0.0_rsh
            DO iv=igrav1,igrav2
              cvolgrv=cvolgrv+cv_sed(iv,k,i,j)/ros(iv)
            ENDDO
            csanmud=0.0_rsh
            DO iv=igrav2+1,nvpc
              csanmud=csanmud+cv_sed(iv,k,i,j)
            ENDDO
            erodab=csanmud*dzs(k,i,j)
            ero=0.0_rsh

            IF(htot(i,j).GT.h0fond) THEN
              IF(tenfon(i,j).GT.toce)THEN
                CALL sed_MUSTANG_comp_eros_flx(tenfon(i,j),toce,c_sedtot(ksmax,i,j),excespowr,xeros,sed_eros_flx)
                ero=sed_eros_flx*fwet(i,j)    

                ! lateral erosion :  wet cell (cellule mouillee)
                IF (coef_erolat .NE. 0.0_rsh .AND. l_erolat_wet_cell) THEN
#if defined key_MARS && defined key_castest_estuary
                  !PLH : precaution for test case ESTUAR (secteur amont ou le courant fluvial peut etre fort)
                  IF(i.le.id(j)-4) THEN
#endif           
                                                           
                    heaue=HTOT_NEAR_E-htncrit_eros
                    heauw=HTOT_NEAR_W-htncrit_eros
                    heaus=HTOT_NEAR_S-htncrit_eros
                    heaun=HTOT_NEAR_N-htncrit_eros
                    heau_milieu=htot(i,j)-h0fond
                    eroe=max(0.0_rsh,coef_tenfon_lat/4.0_rsh*V_NEAR_E**2-toce)*max(0.0_rsh,heaue-heau_milieu)
                    erow=max(0.0_rsh,coef_tenfon_lat/4.0_rsh*V_NEAR_W**2-toce)*max(0.0_rsh,heauw-heau_milieu)
                    eros=max(0.0_rsh,coef_tenfon_lat/4.0_rsh*U_NEAR_S**2-toce)*max(0.0_rsh,heaus-heau_milieu)
                    eron=max(0.0_rsh,coef_tenfon_lat/4.0_rsh*U_NEAR_N**2-toce)*max(0.0_rsh,heaun-heau_milieu)
                    ero=ero+coef_erolat*(eroe+erow+eros+eron)
#ifdef key_castest_estuary
                  ENDIF
#endif                                                                         
                ENDIF
              ENDIF
            ELSE
              ! lateral erosion :  dry cell 
              IF (coef_erolat .NE. 0.0_rsh) THEN
#ifdef key_castest_estuary
               !PLH : precaution for test case ESTUAR (secteur amont ou le courant fluvial peut etre fort)
               IF(i.le.id(j)-2) THEN                                                                         
#endif
                 heaue=max(0.,HTOT_NEAR_E-htncrit_eros)
                 heauw=max(0.,HTOT_NEAR_W-htncrit_eros)
                 heaus=max(0.,HTOT_NEAR_S-htncrit_eros)
                 heaun=max(0.,HTOT_NEAR_N-htncrit_eros)
                 eroe=coef_erolat*max(0.0_rsh,coef_tenfon_lat/4.0_rsh*V_NEAR_E**2-toce)*heaue
                 erow=coef_erolat*max(0.0_rsh,coef_tenfon_lat/4.0_rsh*V_NEAR_W**2-toce)*heauw
                 eros=coef_erolat*max(0.0_rsh,coef_tenfon_lat/4.0_rsh*U_NEAR_S**2-toce)*heaus
                 eron=coef_erolat*max(0.0_rsh,coef_tenfon_lat/4.0_rsh*U_NEAR_N**2-toce)*heaun
                 ero=eroe+erow+eros+eron
#ifdef key_castest_estuary
               ENDIF
#endif
              ENDIF
            ENDIF

            erosi=MF*dt1*ero

            IF(erosi >  0.0_rsh)THEN
              IF(erosi.LE.erodab)THEN
                ddzs=erosi/csanmud
                dzsa=dzs(k,i,j)
                dzs(k,i,j)=MAX(dzsa-ddzs,dzsa*cvolgrv/cvolmaxsort)
                IF(dzs(k,i,j).LE.0.0_rsh)THEN
                  !WRITE(*,*)'en kij:',k,i,j,dzs(k,i,j),dzsa,ddzs,csanmud,cvolgrv
                  dzs(k,i,j)=epsilon
                ENDIF
                dzsisu=1.0_rsh/dzs(k,i,j)
                DO iv=igrav1,igrav2
                  cv_sed(iv,k,i,j)=cv_sed(iv,k,i,j)*dzsa*dzsisu
                ENDDO
                DO iv=igrav2+1,nv_use ! nv_use=nvp (or nvpc if key_sedim_partial)
                  dflx_s2w=ddzs*cv_sed(iv,k,i,j)
                  flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+dflx_s2w/MF 
                  cv_sed(iv,k,i,j)=(cv_sed(iv,k,i,j)*dzsa-dflx_s2w)*dzsisu
                ENDDO
                c_sedtot(k,i,j)=0.0_rsh
                cvolp=0.0_rsh
                DO iv=1,nvpc
                  c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
                  cvolp=cvolp+cv_sed(iv,k,i,j)/ros(iv)
                ENDDO
                poroa=poro(k,i,j)
                poro(k,i,j)=1.0_rsh-cvolp

#if ! defined key_nodiss_insed
                ! dissolved variable in pore waters
#if ! defined key_sedim_partial
                DO iv=nvp+1,nv_adv
                  dflusve=ddzs*poroa*cv_sed(iv,k,i,j)
                  flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+dflusve
                  ! Fev2017 BT : qd il y a erosion, pourquoi conc dissouste changerait elle dans les eaux intertertitielles ?
                  ! cv_sed(iv,k,i,j)=(cv_sed(iv,k,i,j)*dzsa*poroa-dflusve)*dzsisu/poro(k,i,j) 
                ENDDO
#endif

                ! Temperature and salinity
                DO iv=-1,0
                  dflusve=ddzs*poroa*cv_sed(iv,k,i,j)
                  flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+dflusve
                  ! Fev2017 BT : qd il y a erosion, pourquoi conc dissouste changerait elle dans les eaux intertertitielles ?
                  ! cv_sed(iv,k,i,j)=(cv_sed(iv,k,i,j)*dzsa*poroa-dflusve)*dzsisu/poro(k,i,j)
                ENDDO
#endif

#if ! defined key_nofluxwat_IWS
                ! flux d eau du a l erosion
                phieau_s2w(i,j)=phieau_s2w(i,j)+REAL(ddzs*poroa*surf(i,j),rlg)
#endif
              ELSE IF(csanmud /= 0.0_rsh) THEN  !erosi/erodab
                ! correction PLH Juillet 2015
                !volerod=erodab*fwet(i,j)/csanmud
                volerod=erodab/csanmud
                DO iv=igrav2+1,nv_use ! nv_use=nvp (or nvpc if key_sedim_partial)
                  flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+volerod*cv_sed(iv,k,i,j)/MF
                  cv_sed(iv,k,i,j)=0.0_rsh
                ENDDO
#if ! defined key_nodiss_insed
#if ! defined key_sedim_partial
                ! dissolved variable in pore waters
                DO iv=nvp+1,nv_adv
                  flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+volerod*poro(k,i,j)*cv_sed(iv,k,i,j)
                  cv_sed(iv,k,i,j)=0.0_rsh
                ENDDO
#endif
                ! Temperature and salinity
                DO iv=-1,0
                  flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+volerod*poro(k,i,j)*cv_sed(iv,k,i,j)
                  cv_sed(iv,k,i,j)=0.0_rsh
                ENDDO
#endif
#if ! defined key_nofluxwat_IWS
               ! flux d eau du a l erosion
                phieau_s2w(i,j)=phieau_s2w(i,j)+ REAL(volerod*poro(k,i,j)*surf(i,j),rlg)
#endif
                IF(cvolgrv.LE.0.0_rsh)THEN
                  ! mise a zero des dzs et des poro
                  dzs(ksmax,i,j)=0.0_rsh
                  poro(ksmax,i,j)=0.0_rsh
                  ksmax=ksmax-1
                ELSE
                  dzsa=dzs(k,i,j)
                  dzs(k,i,j)=dzsa*cvolgrv/cvolmaxsort
                  !??? dzs(k,i,j)=dzsa*cvolgrv/cvolmaxmel
                  c_sedtot(k,i,j)=0.0_rsh
                  DO iv=igrav1,igrav2
                      cv_sed(iv,k,i,j)=cv_sed(iv,k,i,j)*dzsa/dzs(k,i,j)
                      c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
                  ENDDO
                  poroa=poro(k,i,j)
                  poro(k,i,j)=1.0_rsh-cvolmaxsort
                  !??? poro(k,i,j)=1.0_rsh-cvolmaxmel
#if ! defined key_nodiss_insed
#if ! defined key_sedim_partial
                  ! dissolved variable in pore waters
                  DO iv=nvp+1,nv_adv
                    flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+cv_sed(iv,k,i,j)*(dzsa*poroa   &
                                   -dzs(k,i,j)*poro(k,i,j))
                  ENDDO
#endif

                  ! Temperature and salinity
                  DO iv=-1,0
                    flx_s2w(iv,i,j)=flx_s2w(iv,i,j)+cv_sed(iv,k,i,j)*(dzsa*poroa   &
                                   -dzs(k,i,j)*poro(k,i,j))
                  ENDDO
#endif
#if ! defined key_nofluxwat_IWS
                  ! flux d eau du a l erosion
                  phieau_s2w(i,j)=phieau_s2w(i,j)+                         &
                                 REAL((dzsa*poroa-dzs(k,i,j)*poro(k,i,j))*surf(i,j),rlg)
#endif                                 
                ENDIF  ! fin boucle sur cvolgrv
              
                ! memorisation of lateral erosion for dry cell
                IF(htot(i,j).LE.h0fond .AND. coef_erolat .NE. 0.0_rsh) THEN
                   DO iv=-1,nv_adv
                     flx_s2w_corip1(iv,i,j)=flx_s2w_corip1(iv,i,j)+flx_s2w(iv,i,j)*eroe/ero*surf(i,j)/SURF_NEAR_E
                     flx_s2w_corim1(iv,i,j)=flx_s2w_corim1(iv,i,j)+flx_s2w(iv,i,j)*erow/ero*surf(i,j)/SURF_NEAR_W
                     flx_s2w_corjm1(iv,i,j)=flx_s2w_corjm1(iv,i,j)+flx_s2w(iv,i,j)*eros/ero*surf(i,j)/SURF_NEAR_S
                     flx_s2w_corjp1(iv,i,j)=flx_s2w_corjp1(iv,i,j)+flx_s2w(iv,i,j)*eron/ero*surf(i,j)/SURF_NEAR_N
                   ENDDO
                ENDIF
                ! test on active layer...
                profsed=0.0_rsh
                DO ksup=k,ksmax
                  profsed=profsed+dzs(ksup,i,j)
                ENDDO
                k=k-1
                IF(profsed.LE.activlayer.AND.k.GE.ksmi(i,j))THEN
                  dt1=dt1-erodab/ero
                  GOTO 2
                ENDIF ! fin test activlayer

              ENDIF   ! fin test erosi > erodab
            ENDIF    ! fin test erosi=0
          ENDIF       ! fin test ksmax > ksmi(i,j) 
 
          ! fusion des 2 couches superf., tant que dzs(ksmax) < dzsmin
          !     et si l_consolid=TRUE et couche diminuee par erosion
          !     dans ce cas, il faudrait rajouter un test 
          !    sur la difference entre les concentrations des 2 couches superficielles
          ! pas encore fait
          ! ----------------------------------------------------------
          IF(ero > 0.0_rsh) THEN
 1000       IF(ksmax.GT.ksmi(i,j))THEN
              IF(dzs(ksmax,i,j).LE.dzsmin .AND. (dzs(ksmax-1,i,j)< dzsmax(i,j)  .OR. l_consolid))THEN
                ksmaxa=ksmax
                dzsa=dzs(ksmax,i,j)
                dzsam1=dzs(ksmax-1,i,j)
                poroa=poro(ksmax,i,j)
                poroam1=poro(ksmax-1,i,j)
                ksmax=ksmax-1
                dzs(ksmax,i,j)=dzsa+dzsam1
                dzsisu=1.0_rsh/dzs(ksmax,i,j)
                c_sedtot(ksmax,i,j)=0.0_rsh
                cvolp=0.0_rsh
                DO iv=1,nv_use ! nv_use=nvp (or nvpc if key_sedim_partial)
                  cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmaxa,i,j)*dzsa             &
                                     +cv_sed(iv,ksmax,i,j)*dzsam1)*dzsisu
                  c_sedtot(ksmax,i,j)=c_sedtot(ksmax,i,j)                      &
                                    +cv_sed(iv,ksmax,i,j)*typart(iv)
                  cvolp=cvolp+cv_sed(iv,ksmax,i,j)*typart(iv)/ros(iv)
                ENDDO
                poro(ksmax,i,j)=1.0_rsh-cvolp
                dzpoi=dzsisu/poro(ksmax,i,j)
#if ! defined key_nodiss_insed
#if ! defined key_sedim_partial
                ! dissolved variable in pore waters
                DO iv=nvp+1,nv_adv
                    ! in .../m3 of interstial waters
                   cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmaxa,i,j)*poroa*dzsa+        &
                                    cv_sed(iv,ksmax,i,j)*poroam1*dzsam1)*dzpoi
                ENDDO
#endif

                ! Temperature and salinity
                DO iv=-1,0
                  cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmaxa,i,j)*poroa*dzsa+        &
                                    cv_sed(iv,ksmax,i,j)*poroam1*dzsam1)*dzpoi
                ENDDO
#endif
              
                GOTO 1000
              ENDIF
            ENDIF
          ENDIF

          ! pour retrouver un flux d erosion en .../m2/s (pour les particulaires seulement):
          DO iv=1,nv_use ! nv_use=nvp (or nvpc if key_sedim_partial)
            flx_s2w(iv,i,j)=flx_s2w(iv,i,j)*dtinv
          ENDDO
          !==============================================
#if ! defined key_nodiss_insed
          ! for dissolved variable in pore waters
          ! add water-to-sediment flux explicitly treated with the sediment-to-water flux also explicitly treated
#if ! defined key_sedim_partial
          !    pour les dissoutes, on ajoute le flux diffusif a l interface
          !    et celui qui resulte de l expulsion d eau interst. par tasst:
          DO iv=nvp+1,nv_adv
            flx_s2w(iv,i,j)=(flx_s2w(iv,i,j)-flx_w2s(iv,i,j)+fluconsol(iv,i,j)-fludif(iv,i,j))*dtinv
          ENDDO
#endif

          DO iv=-1,0
            flx_s2w(iv,i,j)=(flx_s2w(iv,i,j)-flx_w2s(iv,i,j)+fluconsol(iv,i,j)-fludif(iv,i,j))*dtinv
          ENDDO
#endif
          ! updating ksma (as ksmax can be modified in the routine)
          ksma(i,j)=ksmax

          ! no correction of water column height (could be !)

#ifdef key_MARS
         ENDIF  ! fin boucle sur jb,jh 
#endif

      END DO
   END DO
!$OMP END DO

#ifdef key_MARS
   PRINT_DBG*, 'FIN SED_EROSION'
#endif

  END SUBROUTINE sed_MUSTANG_erosion
   !!==============================================================================

  SUBROUTINE sed_MUSTANG_comp_tocr_mixsed(k,i,j,xeros,excespowr,sed_tocr_mixsed)
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  SUBROUTINE sed_comp_tocr_mixsed  ***
   !&E
   !&E ** Purpose : the critical shear stress for erosion in sediment layer k
   !&E
   !&E ** Description : 
   !&E
   !&E ** Called by :  sed_erosion
   !&E 
   !&E ** External calls : none
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-02  (P. Le Hir) extracted from (SiAM) new "Marennes-Oleron"
   !&E                               sand/mud sediment model
   !&E       !  2016-02  (B. Mengual, B. Thouvenin) : modification of erosion formulations - several choices
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used


   !! * Arguments
   INTEGER,INTENT(in) :: k,i,j
   REAL(KIND=rsh),INTENT(out) :: xeros,excespowr,sed_tocr_mixsed

   !! * Local declarations
   INTEGER            :: iv
   REAL(KIND=rsh)     :: sommud,somsan,diamsan,critstressan,frvolsan,frvolgrv,somgrav,   &
                         somsedsusp,frvolsangrv,frmudsup,frmudcr1,pinterp,cmudr,frgrav
   ! utilisation de la relation de van Rijn pour flux d erosion (valable seulement si sable seul)
   REAL(KIND=rsh)     :: rossan,xeromud,coef_tmp,E0_sand
   !!---------------------------------------------------------------------------
   !! * Executable part

   sommud=0.0_rsh     
   DO iv=imud1,imud2
     sommud=sommud+cv_sed(iv,k,i,j)
   ENDDO
   somsan=0.0_rsh
   somgrav=0.0_rsh
   diamsan=0.0_rsh
   critstressan=0.0_rsh
   frvolsan=0.0_rsh   
   frvolgrv=0.0_rsh         
   !rossan=0.0_rsh
   !!!!!!!!!!!!!!!!!
!!! ATTENTION : MEME SI PLUSIEURS SABLES, ON SUPPOSE QU ILS ONT TOUS LA MEME DENSITE
   rossan=ros_sand_homogen
!!!!!!!!!!!!!!!!!
!MD 12/09/2017
   DO iv=isand1+1,isand2
!END MD
     somsan=somsan+cv_sed(iv,k,i,j)
     diamsan=diamsan+diam_sed(iv)*cv_sed(iv,k,i,j)
     critstressan=critstressan+stresscri0(iv)*cv_sed(iv,k,i,j)
     frvolsan=frvolsan+cv_sed(iv,k,i,j)/ros(iv)
     !rossan=rossan+ros(iv)*cv_sed(iv,k,i,j)
   ENDDO
   somsedsusp=sommud+somsan
   IF(somsedsusp.LE.0.0_rsh)THEN
     ! surficial sediment non transportable in suspension
     sed_tocr_mixsed=1000.0_rsh
   ELSE
     IF (isand2>0 .AND. somsan>0.0_rsh) THEN
       diamsan=MAX(diamsan/(somsan+epsilon),diam_sed(isand2))
       critstressan=MAX(critstressan/(somsan+epsilon),stresscri0(isand2))
     ENDIF
    
     DO iv=igrav1,igrav2
       somgrav=somgrav+cv_sed(iv,k,i,j)
       frvolgrv=frvolgrv+cv_sed(iv,k,i,j)/ros(iv)
     ENDDO
     frvolsangrv=frvolgrv+frvolsan
     frgrav=somgrav/(somgrav+sommud+somsan)
     IF(frgrav == 1._rsh) THEN
!****************************************!
! Il n y a que du gravier ! 
!****************************************!
       xeros=0.0_rsh
       excespowr=n_eros_sand
       sed_tocr_mixsed=critstressan
     ELSE      
!**************************************************************!
! Il y a  potentiellement du sable et/ou de la vase            !
!                                                              ! 
!   ATTENTION : les formules proposees pour l erosion          !
!               des melanges sablo-vaseux sont discutables     !
!               et doivent etre testees, ameliorees, changees  !
!               chaque utilisateur peut ecrire sa loi          !  
!**************************************************************!
       frmudsup=sommud/(sommud+somsan)
       frmudcr1=MIN(coef_frmudcr1*diamsan,frmudcr2)
       
       IF(frmudsup.LE.frmudcr1)THEN
       !************************************!
       !************************************!
       ! I) Comportement sableux du melange !
       !************************************!
       !************************************!
       !  Si E0_sand_option==0 :  utilisation de E0_sand lu dans parasedim
          IF (E0_sand_option == 2) THEN
            ! E0_sand calcule par formulation erodimetrie en fonction du diametre moyen des sables
              E0_sand = E0_sand_para * critstressan**n_eros_sand * MIN(0.27,1000*diamsan-0.01)
          ELSE IF(E0_sand_option == 0) THEN    !utilisation de E0_sand_para lu dans parasedim
              E0_sand = E0_sand_para
          ELSE IF(E0_sand_option == 1 .AND. rossan.GT.1000_rsh) THEN
            ! E0_sand calcule a partir de la formule de Van Rijn (1984)
              E0_sand=0.000005_rsh*rossan*((rossan/1000.0_rsh-1.0_rsh)*9.81*diamsan)**0.5_rsh &
                     *(diamsan*((rossan/1000.0_rsh-1)*9.81_rsh/(0.000001_rsh)**2)**(1.0_rsh/3.0_rsh))**(0.3_rsh)
          ENDIF
          IF (ero_option == 0) THEN
              cmudr=sommud/(1.0_rsh-frvolsangrv)
              xeros=E0_mud
              sed_tocr_mixsed=x1toce_mud*cmudr**x2toce_mud
              excespowr=n_eros_mud
          ELSE
              xeros=E0_sand
              sed_tocr_mixsed=critstressan
              excespowr=n_eros_sand
          ENDIF
          

       ELSE IF(frmudsup.LE.frmudcr2)THEN
       !****************************************!
       !****************************************!
       ! II) Intermediaire sable / vase         !
       !****************************************!
       !****************************************!
          cmudr=sommud/(1.0_rsh-frvolsangrv)
       !  Si E0_sand_option==0 :  utilisation de E0_sand lu dans parasedim
          IF (E0_sand_option == 2) THEN
            ! E0_sand calcule par formulation erodimetrie en fonction du diametre moyen des sables
              E0_sand = E0_sand_para * critstressan**n_eros_sand * MIN(0.27,1000*diamsan-0.01)
          ELSE IF(E0_sand_option == 0) THEN    !utilisation de E0_sand_para lu dans parasedim
              E0_sand = E0_sand_para
          ELSE IF(E0_sand_option == 1 .AND. rossan.GT.1000_rsh) THEN
            ! E0_sand calcule a partir de la formule de Van Rijn (1984)
              E0_sand=0.000005_rsh*rossan*((rossan/1000.0_rsh-1.0_rsh)*9.81*diamsan)**0.5_rsh &
                     *(diamsan*((rossan/1000.0_rsh-1)*9.81_rsh/(0.000001_rsh)**2)**(1.0_rsh/3.0_rsh))**(0.3_rsh)
          ENDIF
          pinterp=(frmudcr2-frmudsup)/(frmudcr2-frmudcr1)

          IF (ero_option == 0) THEN

            !!! On garde des parametres constants (par defaut ceux definis pour la vase)
            xeros=E0_mud
            sed_tocr_mixsed=x1toce_mud*cmudr**x2toce_mud
            excespowr=n_eros_mud

          ELSE IF (ero_option == 1) THEN

            !!! Formulation initiale avec transition lineaire dans le cas dun melange intermediaire
            !!! Le Hir et al (2011)
            
            xeros=pinterp*E0_sand+(1.0_rsh-pinterp)*E0_mud
            excespowr=pinterp*n_eros_sand+(1.0_rsh-pinterp)*n_eros_mud
            sed_tocr_mixsed=critstressan*pinterp+(1.0_rsh-pinterp)*x1toce_mud*cmudr**x2toce_mud

          ELSE IF (ero_option .EQ. 2) THEN

            !!! Formulation de Julie Vareilles (2013) 
            !!! En accord avec Van Ledden (2001) et Carniello et al (2012)
            xeromud=x1toce_mud*cmudr**x2toce_mud
            xeros = ( (E0_sand/E0_mud)**pinterp )*E0_mud
            sed_tocr_mixsed = ( (critstressan/xeromud)**pinterp )*xeromud
            excespowr = ( (n_eros_sand/n_eros_mud)**pinterp )*n_eros_mud

          ELSE IF (ero_option .EQ. 3) THEN

            !!! Formulation de Baptiste Mengual (2015)
            !!! Formulations exponentielles des parametres pour un melange sablovaseux            
            !xeromud=x1toce_mud*cmudr**x2toce_mud
            !coef_tmp=((-5.0_rsh*frmudsup+(frmudcr1/0.2_rsh)) / (frmudcr2-0.1_rsh))*xexp_ero 
            ! nouvelle version Baptiste Mengual du 04/05/2016 introduit en decembre 2016
            !coef_tmp=(frmudcr1-frmudsup)/(frmudcr2-frmudcr1)
            !!! Parametre d erodabilite
            !xeros = (E0_sand-E0_mud)*EXP(coef_tmp) + E0_mud       
            !!! Contrainte critique pour l erosion
            !sed_tocr_mixsed = ABS((-critstressan +xeromud) * EXP( coef_tmp ) - xeromud)
            !!! Puissance
            !excespowr = (n_eros_sand-n_eros_mud)*EXP( coef_tmp )+ n_eros_mud
            
           !!! Formulation de Baptiste Mengual (2017)
            !!!Mengual, B.; Le Hir, P.; Cayocca, F.; Garlan, T. Modelling Fine Sediment Dynamics: Towards a Common Erosion Law for Fine Sand, Mud and Mixtures. Water 2017, 9, 564.
            !!! Formulations exponentielles des parametres pour un melange sablovaseux            
            xeromud=x1toce_mud*cmudr**x2toce_mud
            coef_tmp=(frmudcr1-frmudsup)/(frmudcr2-frmudcr1) 
            !!! Parametre d erodabilite
            xeros = (E0_sand-E0_mud)*EXP(coef_tmp*xexp_ero) + E0_mud       
            !!! Contrainte critique pour l erosion
            sed_tocr_mixsed = (critstressan -xeromud) * EXP( coef_tmp*xexp_ero ) + xeromud
            !!! Puissance
            excespowr = (n_eros_sand-n_eros_mud)*EXP( coef_tmp*xexp_ero)+ n_eros_mud


          END IF

       ELSE
       !****************************************!
       !****************************************!
       ! III) Comportement vaseux               !
       !****************************************!
       !****************************************!

          cmudr=sommud/(1.0_rsh-frvolsangrv)
          sed_tocr_mixsed=x1toce_mud*cmudr**x2toce_mud
          excespowr=n_eros_mud
          xeros=E0_mud

       ENDIF

     ENDIF
   ENDIF
   END SUBROUTINE sed_MUSTANG_comp_tocr_mixsed
   !!==============================================================================

  SUBROUTINE sed_MUSTANG_comp_eros_flx(tenfo,toce,csurface,excespowr,xeros,sed_eros_flx)
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  SUBROUTINE sed_erosion_flux  ***
   !&E
   !&E ** Purpose : computes erosion flux for given stress, critical stress and
   !&E              surficial sediment concentration
   !&E              uses toce, excespowr and xeros computed by sed_tocr
   !&E
   !&E ** Description : Partheniades or adapted Partheniades
   !&E
   !&E ** Called by :  sed_erosion
   !&E 
   !&E ** External calls : nonee
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-02  (P. Le Hir) extracted from (SiAM) new "Marennes-Oleron"
   !&E                               sand/mud sediment model
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
   !! * Arguments
   REAL(KIND=rsh),INTENT(in) :: csurface,tenfo,toce
   REAL(KIND=rsh),INTENT(in) :: excespowr,xeros
   REAL(KIND=rsh),INTENT(out) :: sed_eros_flx

   !! * Local declarations   
   REAL(KIND=rsh) ::  corfluero

   !!---------------------------------------------------------------------------
   !! * Executable part
   
   corfluero=1.0_rsh+corfluer1*MAX(0.0_rsh,corfluer2-csurface)     
   ! Partheniades :
   sed_eros_flx=xeros*corfluero*(tenfo/toce-1.0_rsh)**excespowr
!choisir une des formules
!  sed_eros_flx=xeros*corfluero*(tenfo-toce)**excespowr
    
  END SUBROUTINE sed_MUSTANG_comp_eros_flx

   !!==============================================================================

  SUBROUTINE sed_MUSTANG_effdep(ifirst,ilast,jfirst,jlast,surf,h0fond)
    
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_effdep  ***
   !&E
   !&E ** Purpose : computes effective deposition
   !&E
   !&E ** Description : 
   !&E       arguments IN :
   !&E          surf,h0fond
   !&E       variables IN :
   !&E           flx_w2s   
   !&E           flx_w2s_corim1, _cip1, _corjp1, _corjm1, corin  
   !&E           flx_w2s_sum
   !&E          
   !&E       variables OUT :
   !&E           phieau_s2w
   !&E           flx_w2s
   !&E
   !&E
   !&E ** Called by :  sed_MUSTANG_deposition
   !&E 
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (P. Le Hir) extracted from subsedim.F90 (from F. Dufois)
   !&E       !  2008-10  (F. Cayocca) included into module_sedimento
   !&E       !  2008-02  (P. Le Hir)replaced by new "Marennes-Oleron" sand/mud sediment model
   !&E       !  2010-07  (B. Thouvenin) taking into account of non constitutive particulate variables
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2015-10  (F. Grasso) Modification from B. Mengual on mud deposition
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
#ifdef key_MARS
   USE comvars3d,    ONLY : dsigu
#if defined key_siggen || defined key_gencoord
   USE toolgeom,     ONLY : f_hzu
#endif
#endif

   !! * Arguments
   INTEGER, INTENT(IN)  :: ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),INTENT(IN) :: h0fond
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_m2p2),INTENT(IN) :: surf

   !! * Local declarations
   REAL(KIND=rsh),DIMENSION(nvp)  :: flx_w2s_loc
   REAL(KIND=rsh),DIMENSION(nvpc) :: frdep          
   INTEGER                        :: i,j,k,iv,ksmin,ksmax,kl,ivp_assoc
   REAL(KIND=rsh) :: fludep,ddzs,ddzs1,ddzs2,dzsa,dzsisu,dzsi,dzsaici,      &
                     dzsgrv,dzssan,dzsmud,dzsnew,dzsa2,dzsa3,ddzs3,ddzs4,ddzs5,ddzs6,    &
                     poroa,somalp,somala,cdepiv,porosi,dzpoi,porodep,fluxdissous,        &
                     ddzsici,voldepgrv,voldepsan,masdepmud,dzskmanew,porewatera,         &
                     sommud,cvolinigrv,cvolinisan,cvolp,cmudr,porewater,                 &
                     dvolgrv,dvolsan,dmasmud,dmasmudgrav,dmasmudsand,cordepflu,          &
                     cordepfluw,cordepflue,cordepflus,cordepflun,                        &
                     dmasmud1,dmasmud2,dmasmud3,dmasmud4
   REAL(KIND=rlg),PARAMETER       :: epsdepef=1.e-14,epsd=1.e-16

   !!----------------------------------------------------------------------
   !! * Executable part

   ddzsici=-1000.0_rsh
   flx_w2s_loc(:)=0.0_rsh
   frdep(:)=0.0_rsh

!$OMP DO SCHEDULE(RUNTIME)   
      DO j=jfirst,jlast
#ifdef key_MARS
        DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
#else
        DO i=ifirst,ilast
#endif
          IF(htot(i,j) > h0fond) THEN

            ksmax=ksma(i,j)
            ksmin=ksmi(i,j)

            ! updating effective deposition :  (flx_w2s) is implicit in the vertical advection scheme
            ! hence multiplied here by the newly computed concentration at each sub time step and cumulated
            DO iv=igrav2+1,nv_use   ! nv_use=nvp (or nvpc if key_sedim_partial)
              flx_w2s_loc(iv)=MAX(0.0_rsh,flx_w2s_sum(iv,i,j)+flx_w2s_corin(iv,i,j)      &
                           +flx_w2s_corip1(iv,i-1,j)+flx_w2s_corim1(iv,i+1,j)         &
                           +flx_w2s_corjp1(iv,i,j-1)+flx_w2s_corjm1(iv,i,j+1))

              flx_w2s_sum(iv,i,j)=0.0_rsh

              ! MF /= 1 only if l_morphocoupl
              flx_w2s_loc(iv)=MF*flx_w2s_loc(iv)
            ENDDO

            fludep=0.0_rsh
            DO iv=1,nvpc
             fludep=fludep+flx_w2s_loc(iv)
            ENDDO

            voldepgrv=0.0_rsh
            voldepsan=0.0_rsh
            masdepmud=0.0_rsh
            dmasmud1 = 0.0_rsh
            dmasmud2 = 0.0_rsh
            dmasmud3 = 0.0_rsh
            dmasmud4 = 0.0_rsh
              
            IF(fludep.GT.0.0_rsh)THEN

             ! case 1:  there is deposition
             ! ****************************

             ! characterising deposition:
             DO iv=igrav1,igrav2
              voldepgrv=voldepgrv+flx_w2s_loc(iv)/ros(iv)
             ENDDO
             DO iv=igrav1,igrav2
              frdep(iv)=flx_w2s_loc(iv)/(voldepgrv+epsd)
             ENDDO
             DO iv=isand1,isand2
              voldepsan=voldepsan+flx_w2s_loc(iv)/ros(iv)
             ENDDO
             DO iv=isand1,isand2
              frdep(iv)=flx_w2s_loc(iv)/(voldepsan+epsd)
             ENDDO
             DO iv=imud1,imud2
              masdepmud=masdepmud+flx_w2s_loc(iv)
             ENDDO
             DO iv=imud1,imud2
              frdep(iv)=flx_w2s_loc(iv)/(masdepmud+epsd)
             ENDDO

             IF(ksmax.GE.ksmi(i,j))THEN
          
              ! case 1.1: sediment already exists
              ! ---------------------------------
              ! characterising previous surficial sediment:

              sommud=0.0_rsh      
              DO iv=imud1,imud2
                sommud=sommud+cv_sed(iv,ksmax,i,j)
              ENDDO
              cvolinigrv=0.0_rsh      
              DO iv=igrav1,igrav2
                cvolinigrv=cvolinigrv+cv_sed(iv,ksmax,i,j)/ros(iv)
              ENDDO
              cvolinisan=0.0_rsh
              DO iv=isand1,isand2
                cvolinisan=cvolinisan+cv_sed(iv,ksmax,i,j)/ros(iv)
              ENDDO
              cmudr=sommud/(1.0_rsh-cvolinisan-cvolinigrv)
            
              IF(cmudr.LE.cmudcr)THEN
            
                ! case 1.1.1: surficial sediment is not consolidated
                ! --------------------------------------------------
                ! mixing deposits with upper layer, until completion:

                dzsa=dzs(ksmax,i,j)

                ! first, gravels deposition :
                ! ---------------------------
                ! dvolgrv represente le volume de graviers qui s insere dans
                ! le sediment existant, en le densifiant
                ! il restera ddzs a deposer a la concentration volumique
                ! cvolmaxsort, dont ddzs2 dans la maille existante, ddzs1 en
                ! une nouvelle couche
                ! attention: a confirmer:
                dvolgrv=MIN(voldepgrv,dzsa*(cvolmaxsort-cvolinigrv),dzsa*(cvolmaxmel-cvolinigrv-cvolinisan))
                dvolgrv=MAX(dvolgrv,0.0_rsh)
                ddzs=(voldepgrv-dvolgrv)/cvolmaxsort
                ddzsici=ddzs
                dzsaici=dzsa
                IF((dzsa+ddzs).LE.dzsmax(i,j))THEN
                  ddzs1=0.0_rsh
                ELSEIF(dzsa.GE.dzsmax(i,j))THEN
                  ddzs1=ddzs
                ELSE
                  ddzs1=dzsa+ddzs-dzsmax(i,j)
                ENDIF              
                !ddzs1=max(0.0_rsh,dzsa+ddzs-dzsmax(i,j))t
                ddzs2=ddzs-ddzs1
                dzsa2=dzsa+ddzs2
                voldepgrv=ddzs1*cvolmaxsort

                ! sand deposition :
                ! -----------------
                ! attention: a confirmer:
                ! modification du calcul de dvolsan (9 janvier 2012)
                !dvolsan=MIN(voldepsan,dzsa2*cvolmaxsort-dzsa*cvolinisan, &
                !                      dzsa2*cvolmaxmel-dzsa*(cvolinigrv+cvolinisan)-ddzs2*cvolmaxsort)
                ! ajout dvolgrv
                dvolsan=MIN(voldepsan,dzsa2*cvolmaxsort-dzsa*cvolinisan, &
                                      dzsa2*cvolmaxmel-dzsa*(cvolinigrv+cvolinisan)-ddzs2*cvolmaxsort)
                dvolsan=MAX(dvolsan,0.0_rsh)
                ! une partie des sandes dvolsab (avec repartition frdep) 
                ! vient densifier l epaisseur dzsa2,l autre partie vient creer 
                ! une surepaisseur ddzs4 (repartition frdep) et
                ! contribuer a une nouvelle couche pour ddzs3 (repartition frdep)
                ddzs=(voldepsan-dvolsan)/cvolmaxsort
                ! suppression de la ligne suivante (9 janvier 2012)
                !voldepsan=0.0_rsh
                IF((dzsa2+ddzs).LE.dzsmax(i,j))THEN
                  ddzs3=0.0_rsh
                ELSEIF(dzsa2.GE.dzsmax(i,j))THEN
                  ddzs3=ddzs
                ELSE
                  ddzs3=dzsa2+ddzs-dzsmax(i,j)
                ENDIF              
                !ddzs3=max(0.0_rsh,dzsa2+ddzs-dzsmax(i,j))
                ddzs4=ddzs-ddzs3
                dzsa3=dzsa2+ddzs4
                voldepsan=ddzs3*cvolmaxsort

                ! mud deposition :
                ! ----------------
                ! une partie des vases dmasmud (avec repartition frdep) 
                ! vient densifier l epaisseur dzsa3, l autre partie vient creer 
                ! une surepaisseur ddzs6 (repartition frdep) et
                ! contribuer a une nouvelle couche pour ddzs5 (repartition frdep)

!!! B.Mengual (17/09/2015): 
!!! Le depot de la vase se fait en commencant a melanger depuis la surface
!!! dans le cas ou un excedant de sables ou de graviers conduira forcement
!!! a la creation dune nouvelle couche
!!! BUT : ne pas pieger de la vase par melange dans la couche ksmax-1 suite
!!!       a la creation dune nouvelle couche


                IF ((voldepgrv+voldepsan) .GT. 0.0_rsh) THEN

                  !!! Part de vase qui va etre melangee avec la future couche
                  ! On enleve les cas a cvolmaxmel pour etre coherent avec les ddzs1 et ddzs3
                  ! calcules avec cvolmaxsort
                  dmasmud1 = MIN(masdepmud,(ddzs1+ddzs3)*cfreshmud*(1-cvolmaxsort))
                  !!! Part de vase melangee avec les depots en sables / graviers effectues 
                  !!! au cours du pas de temps dans la couche ksmax
                  ! Meme commentaire
                  dmasmud2 = MIN(masdepmud-dmasmud1,(dzsa3-dzsa)*cfreshmud*(1-cvolmaxsort))
                  !!! Part de vase melangee avec le sediment present initialement dans la couche ksmax
                  IF((cvolinigrv+cvolinisan) .GT. 0.0_rsh) THEN
                    dmasmud3 = MIN(masdepmud-dmasmud1-dmasmud2,dzsa*(cfreshmud-cmudr)*(1-cvolinigrv-cvolinisan))
                  END IF
                  dmasmud3 = MAX(0.0_rsh,dmasmud3)
                  dmasmud = dmasmud2 + dmasmud3
                  !!! Calcul de l excedant de vase a ajouter dans la futur couche une fois que tout a ete melange
                  dmasmud4 = masdepmud - dmasmud1 - dmasmud
                  !!! Calcul de l epaisseur totale de vase qui sera attribuee a la future couche
                  !!! Attention a ne pas confondre ddzs5 et dzsmud car une partie de ddzs5 sera melangee a lexcedant 
                  !!! de sable gravier avant daccroitre lepaisseur de la couche
                  ddzs5 = (dmasmud4/cfreshmud) + (dmasmud1/cfreshmud)
                  ddzs6 = 0.0_rsh

                  !IF (dmasmud1+dmasmud2+dmasmud3+dmasmud4 .LT. masdepmud) THEN
                  !  print*,'*** BUG ***'
                  !  print*,'La somme des dmasmud est inferieure a masdepmud au point (i,j) : ',i,' ',j
                  !  print*,'dmasmud1,dmasmud2,dmasmud3,dmasmud4,masdepmud : ',dmasmud1,dmasmud2,dmasmud3,dmasmud4,masdepmud
                  !END IF

                  !IF (dmasmud1+dmasmud2+dmasmud3+dmasmud4 .GT. masdepmud) THEN
                  !  print*,'*** BUG ***'
                  !  print*,'La somme des dmasmud est superieure a masdepmud au point (i,j) : ',i,' ',j
                  !  print*,'dmasmud1,dmasmud2,dmasmud3,dmasmud4,masdepmud : ',dmasmud1,dmasmud2,dmasmud3,dmasmud4,masdepmud
                  !END IF


                ELSE
                  dmasmud2 = MIN(masdepmud,(dzsa3-dzsa)*cfreshmud*(1-cvolmaxsort))
                  IF((cvolinigrv+cvolinisan) .GT. 0.0_rsh) THEN
                    dmasmud3 = MIN(masdepmud-dmasmud2,dzsa*(cfreshmud-cmudr)*(1-cvolinigrv-cvolinisan))
                  END IF
                  dmasmud3 = MAX(0.0_rsh,dmasmud3)
                  dmasmud = dmasmud2 + dmasmud3

!                 ! une partie des vases dmasmud (avec repartition frdep) 
!                 ! vient densifier l epaisseur dzsa3, l autre partie vient creer 
!                 ! une surepaisseur ddzs6 (repartition frdep) et
!                 ! contribuer a une nouvelle couche pour ddzs5 (repartition frdep)
!                 IF(cvolinigrv >0.0_rsh .and. cvolinisan > 0.0_rsh) THEN
!                   ! melange sable+gravier  
!                   dmasmud=MIN(masdepmud,(dzsa3-dzsa)*cfreshmud*(1.0_rsh-cvolmaxmel))  ! masse de vase qui s insere dans un melange
!                 ELSE IF ((cvolinisan > 0.0_rsh .AND. isand2>isand1) .OR. (cvolinigrv > 0.0_rsh .AND. igrav2>igrav1) ) THEN
!                   ! melange sables  ou graviers
!                  dmasmud=MIN(masdepmud,(dzsa3-dzsa)*cfreshmud*(1.0_rsh-cvolmaxmel))  ! masse de vase qui s insere dans un melange
!                 ELSE
!                   ! un seul type de sediment 
!                   dmasmud=MIN(masdepmud,(dzsa3-dzsa)*cfreshmud*(1.0_rsh-cvolmaxsort))  ! masse de vase qui s insere dans le gravier ou sable (1-cvolmaxsort)
!                 ENDIF

                  ddzs=(masdepmud-dmasmud)/cfreshmud
                  IF((dzsa3+ddzs).LE.dzsmax(i,j))THEN
                    ddzs5=0.0_rsh
                  ELSEIF(dzsa3.GE.dzsmax(i,j))THEN
                    ddzs5=ddzs
                  ELSE
                    ddzs5=dzsa3+ddzs-dzsmax(i,j)
                  ENDIF
                  !ddzs5=max(0.0_rsh,dzsa3+ddzs-dzsmax(i,j))
                  ddzs6=ddzs-ddzs5

                END IF

                masdepmud=ddzs5*cfreshmud

                dzs(ksmax,i,j)=dzsa3+ddzs6
                dzsi=1.0_rsh/dzs(ksmax,i,j)
                DO iv=igrav1,igrav2
                  cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+(dvolgrv    &
                                     +ddzs2*cvolmaxsort)*frdep(iv))*dzsi
                ENDDO
                DO iv=isand1,isand2
                  cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+(dvolsan    &
                                     +ddzs4*cvolmaxsort)*frdep(iv))*dzsi
                ENDDO
                DO iv=imud1,imud2
                  cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+(dmasmud    &
                                     +ddzs6*cfreshmud)*frdep(iv))*dzsi
                ENDDO

                DO iv=nvpc+1,nvp
#ifdef key_sedim_partial
                  cv_sed(iv,ksmax,i,j)=0.0_rsh
#else
                  ivp_assoc=irkm_var_assoc(iv)
                  IF (ivp_assoc == 0) THEN
                    ! its own settling velocity , totaly integrated in this layer
                    cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+flx_w2s_loc(iv))*dzsi
                  ELSE IF (ivp_assoc < isand1) THEN
                    cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+        &
                                     flx_w2s_loc(iv)/(flx_w2s_loc(ivp_assoc)+epsd)*    &
                                     (dvolgrv+ddzs2*cvolmaxsort)*frdep(ivp_assoc))*dzsi
                  ELSE IF (ivp_assoc < imud1) THEN
                    cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+         &
                                      flx_w2s_loc(iv)/(flx_w2s_loc(ivp_assoc)+epsd)*   &
                                    (dvolsan+ddzs4*cvolmaxsort)*frdep(ivp_assoc))*dzsi
                  ELSE 
                    cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa             &
                                     +flx_w2s_loc(iv)/(flx_w2s_loc(ivp_assoc)+epsd)*    &
                                     (dmasmud+ddzs6*cfreshmud)*frdep(ivp_assoc))*dzsi
                  END IF
#endif
                ENDDO
                k=ksmax
                c_sedtot(k,i,j)=0.0_rsh
                cvolp=0.0_rsh
                DO iv=1,nvpc
                  c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
                  cvolp=cvolp+cv_sed(iv,k,i,j)/ros(iv)
                ENDDO    
                poroa=poro(k,i,j)
                poro(k,i,j)=1.0_rsh-cvolp
          
                !  dissolved variable in pore waters
                !------------------------------------------
                porewatera=poroa*dzsa
                porewater=poro(k,i,j)*dzs(ksmax,i,j)
                IF(porewater > porewatera) THEN
                   !input of pore water from water column
                   dzpoi=dzsi/poro(ksmax,i,j)
#if ! defined key_nodiss_insed
#if ! defined key_sedim_partial
                   DO iv=nvp+1,nv_adv
                     fluxdissous=(porewater-porewatera)*cw_bottom_MUSTANG(iv,i,j)  
                     ! flux dissous ou flx_w2s  sera divise par dt qd pris en compte dans erosion, pour etre en quantite/m2/s
                     cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*porewatera +fluxdissous)*dzpoi
                     flx_w2s(iv,i,j)=flx_w2s(iv,i,j)+fluxdissous
                   ENDDO
#endif
                   ! temperature et salinite
                   fluxdissous=(porewater-porewatera)*temp_bottom_MUSTANG(i,j)  
                   cv_sed(-1,ksmax,i,j)=(cv_sed(-1,ksmax,i,j)*porewatera +fluxdissous)*dzpoi
                   flx_w2s(-1,i,j)=flx_w2s(-1,i,j)+fluxdissous
                   fluxdissous=(porewater-porewatera)*sal_bottom_MUSTANG(i,j)  
                   cv_sed(0,ksmax,i,j)=(cv_sed(0,ksmax,i,j)*porewatera +fluxdissous)*dzpoi
                   flx_w2s(0,i,j)=flx_w2s(0,i,j)+fluxdissous
#endif
#if ! defined key_nofluxwat_IWS
                   ! flux d eau 
                   phieau_s2w(i,j)=phieau_s2w(i,j)+REAL((porewatera-porewater)*surf(i,j),rlg)
#endif
                ELSE
                   !output of pore-water to water column ; cvsed are inchanged
#if ! defined key_nodiss_insed
#if ! defined key_sedim_partial
                   DO iv=nvp+1,nv_adv
                     flx_w2s(iv,i,j)=flx_w2s(iv,i,j)+(porewater-porewatera)*cv_sed(iv,k,i,j)  
                     !flx_w2s sera divise par dt qd pris en compte dans erosion, pour etre en quantite/m2/s
                   ENDDO
#endif
                   ! temperature et salinite
                   flx_w2s(-1,i,j)=flx_w2s(-1,i,j)+(porewater-porewatera)*cv_sed(-1,k,i,j)  
                   flx_w2s(0,i,j)=flx_w2s(0,i,j)+(porewater-porewatera)*cv_sed(0,k,i,j)
#endif
#if ! defined key_nofluxwat_IWS
                   ! flux d eau 
                   phieau_s2w(i,j)=phieau_s2w(i,j)+REAL((porewatera-porewater)*surf(i,j),rlg)
#endif
                ENDIF              
              ENDIF
             ENDIF
  
             ! case 1.2 or following 1.1: computation of remaining deposits
             ! -----------------------------------------------------------

             dzsgrv=0.0_rsh
             dzssan=0.0_rsh
             dzsmud=0.0_rsh
             ! gravel deposition:
             ! ------------------
             IF(voldepgrv.GE.0.0_rsh)THEN
              dzsgrv=voldepgrv/cvolmaxsort
             ENDIF

             ! sand deposition:
             ! ----------------
             IF(voldepsan.GE.0.0_rsh)THEN
              dvolsan=MIN(voldepsan,dzsgrv*(cvolmaxmel-cvolmaxsort))
              dzssan=(voldepsan-dvolsan)/cvolmaxsort
             ENDIF

             ! mud deposition:
             ! ---------------
             IF(masdepmud.GE.0.0_rsh)THEN
               ! depot de la vase
               dmasmud=MIN(masdepmud,(dzsgrv+dzssan)*cfreshmud*(1.0_rsh-cvolmaxmel))
               dzsmud=(masdepmud-dmasmud)/cfreshmud
             ENDIF
          
             dzsnew=dzsgrv+dzssan+dzsmud

             IF(dzsnew.GT.0.0_rsh)THEN
          
              ! case 1.2 or following 1.1: possible creation of a new layer
              ! ----------------------------------------------------------

              !we choose to do the same test with or whithout consolidation. 
              !We will manage very small dzs in consolidation module
              IF(ksmax == 0 .OR. l_consolid) THEN
                dzskmanew=0.0_rsh
              ELSE
                dzskmanew=dzs(ksmax,i,j)+dzsnew
              ENDIF
              ! creation only if dzsnew>dzsmin or if overtake dzsmax 
              IF(((dzsnew >= dzsmin .OR. dzskmanew > dzsmax(i,j)) .AND. ksdmax > 1) .OR. ksmax == 0)THEN

                ! actual constitution of the new layer:
                ! -------------------------------------
                IF(ksmax.EQ.ksdmax)CALL sed_MUSTANG_fusion(i,j,ksmax)
                dzsa=0.0_rsh
                ksmax=ksmax+1
                dzs(ksmax,i,j)=dzsnew
                dzsi=1.0_rsh/dzs(ksmax,i,j)
                DO iv=igrav1,igrav2
                  cv_sed(iv,ksmax,i,j)=voldepgrv*frdep(iv)*dzsi
                ENDDO
                DO iv=isand1,isand2
                  cv_sed(iv,ksmax,i,j)=voldepsan*frdep(iv)*dzsi
                ENDDO
                DO iv=imud1,imud2
                  cv_sed(iv,ksmax,i,j)=masdepmud*frdep(iv)*dzsi
                ENDDO

#if ! defined key_sedim_partial
                DO iv=nvpc+1,nvp
                  ivp_assoc=irkm_var_assoc(iv)
                  IF (ivp_assoc .NE. 0) THEN
                    ! pas de depot si pas associe car tout depose dans la couche juste avant
                    cv_sed(iv,ksmax,i,j)= flx_w2s_loc(iv)/(flx_w2s_loc(ivp_assoc)+epsd)*    &
                                        cv_sed(ivp_assoc,ksmax,i,j)
                  ELSE
                   cv_sed(iv,ksmax,i,j)=0.0_rsh
                  END IF 
                ENDDO
#endif
                ! computing total concentration and porosity of the surf.layer:
                ! -------------------------------------------------------------
                k=ksmax
                c_sedtot(k,i,j)=0.0_rsh
                cvolp=0.0_rsh
                DO iv=1,nvpc
                  c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
                  cvolp=cvolp+cv_sed(iv,k,i,j)/ros(iv)
                ENDDO

                poro(k,i,j)=1.0_rsh-cvolp
              
                !  dissolved variable in pore waters
                !------------------------------------------
                porewater=poro(k,i,j)*dzs(ksmax,i,j)
                !input of pore water from water column
                dzpoi=dzsi/poro(ksmax,i,j)
#if ! defined key_nodiss_insed
#if ! defined key_sedim_partial
                DO iv=nvp+1,nv_adv
                  cv_sed(iv,ksmax,i,j)=cw_bottom_MUSTANG(iv,i,j)
                  flx_w2s(iv,i,j)=flx_w2s(iv,i,j)+porewater*cw_bottom_MUSTANG(iv,i,j)   
                  !sera divise par dt qd pris en compte dans erosion, pour etre en quantite/m2/s
                ENDDO
#endif
                cv_sed(-1,ksmax,i,j)=temp_bottom_MUSTANG(i,j)
                flx_w2s(-1,i,j)=flx_w2s(-1,i,j)+porewater*temp_bottom_MUSTANG(i,j) 
                cv_sed(0,ksmax,i,j)=sal_bottom_MUSTANG(i,j)
                flx_w2s(0,i,j)=flx_w2s(0,i,j)+porewater*sal_bottom_MUSTANG(i,j) 
#endif
#if ! defined key_nofluxwat_IWS
                ! flux d eau 
                phieau_s2w(i,j)=phieau_s2w(i,j)-REAL(porewater*surf(i,j),rlg)
#endif                
              ELSE
            
                ! as the new layer is too small, back to an increase of the surf. layer:
                ! ----------------------------------------------------------------------
                dzsa=dzs(ksmax,i,j)
                dzs(ksmax,i,j)=dzs(ksmax,i,j)+dzsnew
                dzsi=1.0_rsh/dzs(ksmax,i,j)
                DO iv=igrav1,igrav2
                  cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa                 &
                                  +voldepgrv*frdep(iv))*dzsi
                ENDDO
                DO iv=isand1,isand2
                  cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa                 &
                                 +voldepsan*frdep(iv))*dzsi
                ENDDO
                DO iv=imud1,imud2
                  cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa                 &
                                  +masdepmud*frdep(iv))*dzsi
                ENDDO

#if ! defined key_sedim_partial
                DO iv=nvpc+1,nvp
                  ivp_assoc=irkm_var_assoc(iv)
                  IF (ivp_assoc == 0) THEN
                   ! its own settling velocity 

                  ELSE IF (ivp_assoc < isand1) THEN
                   cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+                    &
                                       flx_w2s_loc(iv)/(flx_w2s_loc(ivp_assoc)+epsd)*  &
                                       voldepgrv*frdep(ivp_assoc))*dzsi
                  ELSE IF (ivp_assoc < imud1) THEN
                   cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+                    &
                                     flx_w2s_loc(iv)/(flx_w2s_loc(ivp_assoc)+epsd)*   &
                                     voldepsan*frdep(ivp_assoc))*dzsi
                  ELSE 
                   cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*dzsa+                    &
                                      flx_w2s_loc(iv)/(flx_w2s_loc(ivp_assoc)+epsd)*   &
                                      masdepmud*frdep(ivp_assoc))*dzsi
                  END IF 
                ENDDO
#endif
            
                ! computing total concentration and porosity of the surf.layer:
                ! -------------------------------------------------------------
                k=ksmax
                c_sedtot(k,i,j)=0.0_rsh
                cvolp=0.0_rsh
                DO iv=1,nvpc
                  c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
                  cvolp=cvolp+cv_sed(iv,k,i,j)/ros(iv)
                ENDDO
                poroa=poro(k,i,j)
                poro(k,i,j)=1.0_rsh-cvolp
 
                !  dissolved variable in pore waters
                !------------------------------------------
                porewatera=poroa*dzsa
                porewater=poro(k,i,j)*dzs(ksmax,i,j)
                IF(porewater > porewatera) THEN
                  !input of pore water from water column
                  dzpoi=dzsi/poro(ksmax,i,j)
#if ! defined key_nodiss_insed
#if ! defined key_sedim_partial
                  DO iv=nvp+1,nv_adv
                    fluxdissous=(porewater-porewatera)*cw_bottom_MUSTANG(iv,i,j)  
                    !sera divise par dt qd pris en compte dans erosion, pour etre en quantite/m2/s
                    cv_sed(iv,ksmax,i,j)=(cv_sed(iv,ksmax,i,j)*porewatera +fluxdissous)*dzpoi
                    flx_w2s(iv,i,j)=flx_w2s(iv,i,j)+fluxdissous
                  ENDDO
#endif
                  ! salinite et temperature
                  fluxdissous=(porewater-porewatera)*temp_bottom_MUSTANG(i,j)  
                  cv_sed(-1,ksmax,i,j)=(cv_sed(-1,ksmax,i,j)*porewatera +fluxdissous)*dzpoi
                  flx_w2s(-1,i,j)=flx_w2s(-1,i,j)+fluxdissous
                  fluxdissous=(porewater-porewatera)*sal_bottom_MUSTANG(i,j)  
                  cv_sed(0,ksmax,i,j)=(cv_sed(0,ksmax,i,j)*porewatera +fluxdissous)*dzpoi
                  flx_w2s(0,i,j)=flx_w2s(0,i,j)+fluxdissous                          
#endif
#if ! defined key_nofluxwat_IWS
                  ! flux d eau 
                  phieau_s2w(i,j)=phieau_s2w(i,j)+REAL((porewatera-porewater)*surf(i,j),rlg)
#endif
                ELSE
                !output of pore-water to water column ; cvsed are inchanged
#if ! defined key_nodiss_insed
#if ! defined key_sedim_partial
                  DO iv=nvp+1,nv_adv
                    flx_w2s(iv,i,j)=flx_w2s(iv,i,j)+(porewater-porewatera)*cv_sed(iv,k,i,j)  
                    !sera divise par dt qd pris en compte dans erosion, pour etre en quantite/m2/s
                  ENDDO
#endif
                ! salinite et temperature
                  DO iv=-1,0
                    flx_w2s(iv,i,j)=flx_w2s(iv,i,j)+(porewater-porewatera)*cv_sed(iv,k,i,j)  
                  ENDDO
#endif
#if ! defined key_nofluxwat_IWS
                ! flux d eau 
                  phieau_s2w(i,j)=phieau_s2w(i,j)+REAL((porewatera-porewater)*surf(i,j),rlg)
#endif
                ENDIF  !  porewater > porewatera
              ENDIF  ! dwsnew > dzsmin
 
              ! control writing :            
              !if(poro(k,i,j).ge.0.9999_rsh)then
              !  write(*,*)'dans effdep, a',t,'  en kij:',k,i,j
              !  write(*,555)k,dzs(k,i,j),c_sedtot(k,i,j),poro(k,i,j)
              !  write(*,*)dzsa,dzsgrv,dzssan,dzsmud,voldepgrv,voldepsan, &
              !            masdepmud,frdep,(cv_sed(iv,k,i,j),iv=1,4)
              !  write(*,*)igrav1,igrav2,isand1,isand2,imud1,imud2
              !  write(*,*)(fluevs(iv,i,j),iv=1,5)
              !  write(*,*)cvasr,cvascr,ddzs,ddzs1,ddzs2,ddzs3,ddzs4,ddzs5
              !  write(*,*)ddzsici,dvolgrv,dvolsan,cvolinigrv,cvolinisan,dzsaici
              !endif

             ENDIF  ! dzsnew > 0

            ELSE   !

             ! case 2. : no deposition
             ! ***********************  
               !      write(*,*)'no deposition ',i,j,fludep

               ! computing flx_w2s for dissolved var., to be used in the next time step:
               ! -----------------------------------------------------------------------
#if ! defined key_nodiss_insed
               DO iv=nvp+1,nv_adv
                 flx_w2s(iv,i,j)=0.0_rsh
               ENDDO
               flx_w2s(-1,i,j)=0.0_rsh
               flx_w2s(0,i,j)=0.0_rsh
#endif

               ! accounting for possible deposition of non constitutive particulate variable:
               ! ----------------------------------------------------------------------------
               ! prise en compte des depots eventuels 
               ! de variables particulaires non constitutives:
               IF(ksmax.GE.ksmi(i,j))THEN
                 DO iv=nvpc+1,nvp
                   cv_sed(iv,ksmax,i,j)=cv_sed(iv,ksmax,i,j)+flx_w2s_loc(iv)/dzs(ksmax,i,j)
                 ENDDO
               ELSE
                 ! no sediment
                 DO iv=nvpc+1,nvp
#ifdef key_MARS
#if defined key_siggen || defined key_gencoord
                  cw_bottom_MUSTANG(iv,i,j)=cw_bottom_MUSTANG(iv,i,j)                                        &
                               +flx_w2s_loc(iv)/(epn_bottom_MUSTANG(i,j))
#else
                  cw_bottom_MUSTANG(iv,i,j)=cw_bottom_MUSTANG(iv,i,j)                                        &
                               +flx_w2s_loc(iv)/(epn_bottom_MUSTANG(i,j))
#endif
#else
                ! autre modele hydro?
#endif
 !   il manque ici echange des cw_bottom_MUSTANG pour MPI ???????
                 ENDDO
               ENDIF


            ENDIF

            !GOTO 51
            ! control writing :            
            !DO k=ksmi(i,j),ksmax
            !  IF(poro(k,i,j).GE.0.9999_rsh)THEN
            !    !WRITE(*,*)'a la fin de effdep, a',t,'  en kij:',k,i,j
            !    !DO kl=ksmax,ksmi(i,j),-1
            !      !WRITE(*,555)kl,dzs(kl,i,j),c_sedtot(kl,i,j),poro(kl,i,j)
            !      !555 FORMAT('k=',i3,'  dzs ',g8.2,'  csed ',g8.2,'  poro ',g8.2)
            !    !ENDDO
            !    GOTO 51
            !  ENDIF
            !ENDDO
        51  CONTINUE

            ! updating ksma ( ksmax can be modified in the routine)
            ksma(i,j)=ksmax
        
          END IF ! test on htot
        END DO  ! loop on i
    END DO    ! loop on j
!$OMP END DO

#ifdef key_MARS
   PRINT_DBG*, 'FIN SED_EFFDEP'   
#endif

  END SUBROUTINE sed_MUSTANG_effdep

  !!==============================================================================
  SUBROUTINE sed_MUSTANG_fusion(i,j,ksmaxhere)
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_fusion  ***
   !&E
   !&E ** Purpose : makes fusion of the 2 first layers of the sediment (at the bottom)
   !&E              when the max. number of layers is reached
   !&E              and updating of the sediment column
   !&E
   !&E ** Description : use arguments and common variable 
   !&E  arguments IN : 
   !&E         ksmin :

   !&E  arguments OUT:
   !&E         ksmax : 
   !&E         cs,c_sedtotij:
   !&E          
   !&E  need to be know by hydrodynamic code:
   !&E          
   !&E  need to be know by code treated substance:
   !&E         typart , 
   !&E         nv_adv,nvp                  
   !&E     
   !&E  use module SEDIM1DV variables  :
   !&E         ros
   !&E              
   !&E
   !&E
   !&E ** Called by :  sed_effdep
   !&E 
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-03 (P. Le Hir) reproduces s_fusion in new "Marennes-Oleron" sand/mud sediment model
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   INTEGER,INTENT(IN)     ::  i,j
   INTEGER,INTENT(INOUT)  ::  ksmaxhere

   !! * Local declarations
   INTEGER        ::  k,iv,ksmin
   REAL(KIND=rsh) ::  dzsa,poroa,cvolp,dzspwi
              
   
   !!----------------------------------------------------------------------
   !! * Executable part
      
   ksmin=ksmi(i,j)
   dzsa=dzs(ksmin,i,j)
   poroa=poro(ksmin,i,j)
   dzs(ksmin,i,j)=dzs(ksmin,i,j)+dzs(ksmin+1,i,j)
   cvolp=0.0_rsh
   c_sedtot(ksmin,i,j)=0.0_rsh
   DO iv=1,nvp
     cv_sed(iv,ksmin,i,j)=(cv_sed(iv,ksmin,i,j)*dzsa+cv_sed(iv,ksmin+1,i,j)*dzs(ksmin+1,i,j))/dzs(ksmin,i,j)
     c_sedtot(ksmin,i,j)=c_sedtot(ksmin,i,j)+cv_sed(iv,ksmin,i,j)*typart(iv)
     cvolp=cvolp+cv_sed(iv,ksmin,i,j)*typart(iv)/ros(iv)
   ENDDO
   poro(ksmin,i,j)=1.0_rsh-cvolp
#if ! defined key_nodiss_insed
   dzspwi=1.0_rsh/(dzs(ksmin,i,j)*poro(ksmin,i,j))
   DO iv=-1,0
      cv_sed(iv,ksmin,i,j)=dzspwi*(cv_sed(iv,ksmin,i,j)*dzsa*poroa &
                                +cv_sed(iv,ksmin+1,i,j)*dzs(ksmin+1,i,j)*poro(ksmin+1,i,j))
   ENDDO
#if ! defined key_sedim_partial
   DO iv=nvp+1,nv_adv
      cv_sed(iv,ksmin,i,j)=dzspwi*(cv_sed(iv,ksmin,i,j)*dzsa*poroa &
                                +cv_sed(iv,ksmin+1,i,j)*dzs(ksmin+1,i,j)*poro(ksmin+1,i,j))
   ENDDO
#endif
#endif
   DO k=ksmin+1,ksmaxhere-1
     dzs(k,i,j)=dzs(k+1,i,j)
     poro(k,i,j)=poro(k+1,i,j)
     c_sedtot(k,i,j)=c_sedtot(k+1,i,j)
     DO iv=-1,nv_adv
       cv_sed(iv,k,i,j)=cv_sed(iv,k+1,i,j)
     ENDDO
   ENDDO
   ksmaxhere=ksmaxhere-1
  

  END SUBROUTINE sed_MUSTANG_fusion
      

 
!!==============================================================================

  SUBROUTINE sed_MUSTANG_slipdepo(ifirst,ilast,jfirst,jlast, &
                         dx,dy,surf,h0fond,h0)
   
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_MUSTANG_slipdepo  ***
   !&E
   !&E ** Purpose : computes slipping settling if slope bathy
   !&E
   !&E ** Description :  
   !&E  arguments IN : 
   !&E         h0, dx,dy, h0fond, surf 
   !&E          
   !&E   variables in  (comMUSTANG)      flx_w2s_sum 
   !&E          
   !&E   variables OUT:
   !&E         flx_w2s_corip1,flx_w2s_corim1,flx_w2s_corjp1,flx_w2s_corjm1
   !&E         flx_w2s_corin
   !&E          
   !&E
   !&E ** Called by :  sed_MUSTANG_deposition
   !&E 
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (P. Le Hir) extracted from subsedim.F90 (from F. Dufois)
   !&E       !  2008-10  (F. Cayocca) included into module_sedimento
   !&E       !  2008-02  (P. Le Hir)replaced by new "Marennes-Oleron" sand/mud sediment model
   !&E       !  2010-07  (B. Thouvenin) taking into account of non constitutive particulate variables
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2015-10  (F. Grasso) Modification from B. Mengual on mud deposition
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   INTEGER, INTENT(IN)  :: ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),INTENT(IN)  :: h0fond
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_m2p2),INTENT(IN)  :: surf,dx,dy,h0


   !! * Local declarations
   INTEGER                :: iv,i,j
   REAL(KIND=rsh)         :: cordepfluw,cordepflue,cordepflus,cordepflun,cordepflu


   !!----------------------------------------------------------------------
   !! * Executable part 
!$OMP DO SCHEDULE(RUNTIME)   
   DO j=jfirst,jlast 
   
#ifdef key_MARS
     DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1)
       !!ARR test j.GE. jb+1 & j.LE.jh-1 pour coherence avec erosion
       IF(j.GE.jb(i)+1 .AND. j .LE. jh(i)-1) THEN
#else
     DO i=ifirst,ilast
#endif
        IF(htot(i,j) > h0fond) THEN

         cordepfluw=max(0.0_rsh,slopefac*SLOPE_W)
         cordepflue=max(0.0_rsh,slopefac*SLOPE_E)
         cordepflus=max(0.0_rsh,slopefac*SLOPE_S)
         cordepflun=max(0.0_rsh,slopefac*SLOPE_N)
         cordepflu=cordepfluw+cordepflue+cordepflus+cordepflun
         IF (cordepflu.gt.1.0_rsh)THEN
           cordepfluw=cordepfluw/cordepflu
           cordepflue=cordepflue/cordepflu
           cordepflus=cordepflus/cordepflu
           cordepflun=cordepflun/cordepflu
           cordepflu=1.0_rsh
         ENDIF
         DO iv=isand2+1,nvp
           flx_w2s_corim1(iv,i,j)=cordepfluw*flx_w2s_sum(iv,i,j)*surf(i,j)/SURF_NEAR_W
           flx_w2s_corip1(iv,i,j)=cordepflue*flx_w2s_sum(iv,i,j)*surf(i,j)/SURF_NEAR_E
           flx_w2s_corjm1(iv,i,j)=cordepflus*flx_w2s_sum(iv,i,j)*surf(i,j)/SURF_NEAR_S
           flx_w2s_corjp1(iv,i,j)=cordepflun*flx_w2s_sum(iv,i,j)*surf(i,j)/SURF_NEAR_N
           flx_w2s_corin(iv,i,j)=-cordepflu*flx_w2s_sum(iv,i,j)
         ENDDO
        ENDIF

#ifdef key_MARS
        ENDIF
#endif
      ENDDO
   ENDDO
!$OMP END DO

  END SUBROUTINE sed_MUSTANG_slipdepo
 
    !!==============================================================================

   SUBROUTINE sed_MUSTANG_consol_diff_bioturb(ifirst,ilast,jfirst,jlast,  &
                surf,h0fond,roswat_bot,rhoref,saliref_lin,tetaref_lin,t,dt_true)
! 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_consol_diff_bioturb ***
   !&E
   !&E ** Purpose : dynamic in sediment : processes of consolidation 
   !&E                   and/or diffusion and/or bioturbation
   !&E
   !&E ** Description :
   !&E        arguments IN :
   !&E            surf : surface des mailles (SURF)
   !&E            h0fond : RESIDUAL_THICKNESS_WAT
   !&E            t : temps courrant
   !&E            dt_true : pas de temps utilise par le modele hydro entre deux appels (DOUBLE PRECISION)
   !&E        variales OUT :
   !&E            phieau_s2w : flux d eau a l interface eau/sediment du a la consolidation
   !&E                        (en m3 eau echange pendant le pas de temps effectif de la consolidation)
   !&E            fluconsol : flux de substance a l interface eau/sediment du a la consolidation
   !&E            fludiff : flux de substance a l interface eau/sediment du a la diffusion
   !&E            + remaniement du sediment (cv_sed, poro, dzs, ksma.. modifies)
   !&E
   !&E ** Called by :  sed_MUSTANG_first
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (P. Le Hir) extracted from init.F90
   !&E       !  2008-10  (F. Cayocca) inserted into module_sedimento.F90
   !&E       !  2009-02  (P. Le Hir) new version from SiAM "sable-vase"
   !&E       !  2014-05  (F. Grasso) modifications based on "Grasso et al. (Ocean Dynynamics, 2015)"
   !&E       !  2015-07  (F. Grasso, P. Le Hir)  rewriting             
   !&E       !  2017-03  (B.Thouvenin) fusion of consolidation and diffusion for module mixsed and corrections for dissolved subst
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used


   !! * Arguments
   INTEGER, INTENT(IN)            :: ifirst,ilast,jfirst,jlast                           
   REAL(KIND=rsh),INTENT(IN)      :: h0fond,rhoref,saliref_lin,tetaref_lin
   REAL(KIND=rlg),INTENT(IN)      :: t,dt_true
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_m2p2),INTENT(IN)   :: surf
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY),INTENT(IN)        :: roswat_bot
   
   !! * Local declarations
   INTEGER                                   :: i,j,k,iv,ksmin,ksmax,celldry,ksmaxa,kp,kl,nv_use2
   REAL(KIND=rsh)                            :: volpwa_tot,volpwnv_tot,volpwk,roro,dt_sed_inv,dtsdzs,   &
                                                sigmadjge,somdeltaro,halfmaslayer,somcsan,somcgrav,   &
                                                somcmud,cvolp,rapcsed,csednew,poroam1,poroa,     &
                                                dzsa,dzsam1,dzsisu,rowinv,dtsdzsmin,  &
                                                volak,volakm1,vola,vola_inv
   !REAL(KIND=rsh)                            :: dvol,dvol_fusion 
   REAL(KIND=rsh),DIMENSION(ksdmin  :ksdmax) :: volpwak,stateconsol,permeab,sigmapsg,  &
                                                poroin,dzsiin,loadograv,sed_rate,hinder,cseda
   REAL(KIND=rsh),DIMENSION(ksdmin-1:ksdmax) :: winter
   REAL(KIND=rlg)                            :: dt_sed_cor,dt_sed_eff,dtrest,dtiter                                           
   REAL(KIND=rsh),DIMENSION(nvp,ksdmin:ksdmax)            :: cvsednew
   REAL(KIND=rsh),DIMENSION(ksdmin:ksdmax,-1:nv_adv)      :: difbio                   
#if ! defined key_nodiss_insed
 !  REAL(KIND=rsh),DIMENSION(-1:nv_adv)                    :: flu_fusion              
   REAL(KIND=rsh),DIMENSION(ksdmin  :ksdmax)              :: disvi,aa,bb,cc,dd
   REAL(KIND=rsh)                                         :: hcrit,dzssdt
#endif



   !!----------------------------------------------------------------------
   !! * Executable part
#if ! defined key_sedim_partial
   nv_use2=nv_adv
#else
   nv_use2=0
#endif

   IF(t.GE.tstart_sed .AND. t.GE.t_sed) THEN
   
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!!!! DYNAMIQUE INTERNE DANS LE SEDIMENT    !!!!!!
      !!!!  COSOLIDATION/BIOTURBATION/DIFFUSION   !!!!!!
      !!!!  MISE A JOUR TOUS LES dt_sed     !!!!!!
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     dt_sed_cor=MAX(dt_sed,dt_true)
     dt_sed_eff=dt_sed_cor+(t-t_sed)
     dt_sed_inv=1.0_rsh/REAL(dt_sed_eff,rsh)
!$OMP SINGLE
     phieau_s2w_consol(:,:)=0.0_rsh
!$OMP END SINGLE

     PRINT_DBG*, 'DEBUT DE DIFFUSION et CONSOL et Bioturb '

!$OMP DO SCHEDULE(RUNTIME) PRIVATE(i,j,k,kp,iv,volpwak,celldry,volpwa_tot,  &
!$OMP poroin,dzsiin,ksmax,ksmin,kl,winter,roro,rowinv,sed_rate,hinder,   &
!$OMP volak,volakm1,dtrest,volpwk,sigmapsg,    &
!$OMP cseda,dtiter,sigmadjge,somdeltaro,halfmaslayer,loadograv,permeab,csednew, &
!$OMP stateconsol,somcgrav,somcsan,somcmud,dtsdzs,dtsdzsmin,difbio,cvsednew, &
!$OMP cvolp,rapcsed,dzsa,ksmaxa,dzsam1,poroam1,poroa,dzsisu,vola,vola_inv,  &
#if ! defined key_nodiss_insed
!$OMP disvi,dzssdt,hcrit,aa,bb,cc,dd,                                         &
#endif
!$OMP volpwnv_tot)
!!!$OMP SINGLE
     DO j=jfirst,jlast
#ifdef key_MARS
      DO i=MAX0(ifirst,ig(j)+1),MIN0(ilast,id(j)-1) 
#else
      DO i=ifirst,ilast
#endif
       difbio(ksdmin:ksdmax,-1:nv_adv)=0.0_rsh
       poroin(ksdmin:ksdmax)=0.0_rsh
       volpwak(:)=0.0_rsh
       IF (ksma(i,j).GE.ksmi(i,j))THEN
         celldry=0
         IF(htot(i,j) .LE. h0fond)celldry=1
           
         ! memorisation des porosites, du volume d eau interstitielle
         !  des masses de substances dissoutes
         ! et de l epaisseur de la couche de surface
         ! ----------------------------------------------------------
         volpwak(ksma(i,j))=poro(ksma(i,j),i,j)*dzs(ksma(i,j),i,j)
         volpwa_tot=volpwak(ksma(i,j))
         poroin(ksma(i,j))=poro(ksma(i,j),i,j)
         DO k=ksmi(i,j),ksma(i,j)-1
            dzsiin(k)=2.0_rsh/(dzs(k,i,j)+dzs(k+1,i,j))
            poroin(k)=0.5_rsh*(poro(k+1,i,j)+poro(k,i,j))
            volpwak(k)=poro(k,i,j)*dzs(k,i,j)
            volpwa_tot=volpwa_tot+volpwak(k)
         ENDDO

         winter(:)=0.0_rsh
                      
         ksmax=ksma(i,j)
         ksmin=ksmi(i,j)

      !!! calcul de la densite de l eau 
         ! densite de l eau interstitielle
         ! on la suppose egale a la densite de l eau au fond de la colonne d eau
         ! dans MARS celle ci n est calcule que si htot>hm (3D) par tsstateeqblumberg si l_stateeq_lin=F
         !                        tout le temps par tsstateeqlinear si l_stateeq_lin=T
         ! si il n y a pas d eau la salinite et temperature ne sont pas connus
         ! si defined key_nodiss_insed, temp et sal ne sont pas calcules dans le sediment 
         !  et sont faux (=valmanq) dans les nouvelles couches
         IF(celldry ==1) THEN
         ! pas d eau
         !  densite a partir de la maille du fond ou de surface de sediment
#ifdef key_nodiss_insed
           !roro=(rhoref*(1.0_rsh+0.0008_rsh*(cv_sed(0,1,i,j)-saliref_lin)  &
           !                            -0.00016_rsh*(cv_sed(-1,1,i,j)-tetaref_lin)))
           roro=rhoref
#else
           roro=(rhoref*(1.0_rsh+0.0008_rsh*(cv_sed(0,ksmax,i,j)-saliref_lin)  &
                                       -0.00016_rsh*(cv_sed(-1,ksmax,i,j)-tetaref_lin)))
#endif
         ELSE
         roro=roswat_bot(i,j)
         ENDIF
         rowinv=1.0_rsh/roro

         sed_rate(ksdmin:ksdmax)=0.0_rsh
         hinder(ksdmin:ksdmax)=0.0_rsh
         !dvol=0.0_rsh
         !dvol_fusion=0.0_rsh
         volak=0.0_rsh
         volakm1=0.0_rsh
!#if ! defined key_nodiss_insed
!         flu_fusion(:)=0.0_rsh
!#endif
         
         dtrest=dt_sed_eff
         
  100    CONTINUE
                   
         cseda(:)=0.0_rsh     
         DO k=ksmin,ksmax
              cseda(k)=c_sedtot(k,i,j)
         ENDDO
            
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         !!!!!!!!!!!  CONSOLIDATION   !!!!!!!!!!!!!!!!!!!!!!!!!!
         !!!!!  seulement si plus d une couche   !!!!!!!!!!!!!!!
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         dtiter=MIN(dtrest,dtcal_sed)

         IF(l_consolid .AND. ksma(i,j) > ksmi(i,j)) THEN
             
          ! computation of permeability, effective stress and hydraulic load
          ! ----------------------------------------------------------------
           CALL sed_constitutivrel(i,j,stateconsol,permeab,sigmapsg)        

           sigmadjge=0.0_rsh
           DO k=ksmax,ksmin,-1
              somdeltaro=0.0_rsh
              DO iv=1,nvpc
              ! rosmrowsros=(ros-row)/ros
                somdeltaro=somdeltaro+cv_sed(iv,k,i,j)*rosmrowsros(iv)
              ENDDO
              halfmaslayer=.5*dzs(k,i,j)*somdeltaro
              sigmadjge=sigmadjge+halfmaslayer
              ! loadograv : exces de pression d eau interstitielle au milieu de la couche
              ! sigmadjge : sigma dejauge (sans la part de l eau)
              ! sigmapsg : contrainte effective (transmise de grains a grains)
              loadograv(k)=MAX(0.0_rsh,sigmadjge-sigmapsg(k))
              ! poids total de toutes les couches au dessus de la couche k
              sigmadjge=sigmadjge+halfmaslayer
           ENDDO
          
           ! computing new concentrations from mass conservation
           !  only particulate variables are modified during consolidation
           ! ---------------------------------------------------       
           k=ksmin
           !sed_rate : vitesse d advection des particules de vase
           ! loi de Darcy
           ! stateconsol : indicateur d etat de la consolidation
           sed_rate(k)=-rowinv*(permeab(k+1)+permeab(k))                                 &
                *MIN(0.0_rsh,loadograv(k+1)-loadograv(k))*stateconsol(k)         &
                /(dzs(k+1,i,j)+dzs(k,i,j))
           somcgrav=0.0_rsh
           DO iv=igrav1,igrav2
              somcgrav=somcgrav+cv_sed(iv,k,i,j)
           ENDDO
           somcsan=0.0_rsh
           DO iv=isand1,isand2
              somcsan=somcsan+cv_sed(iv,k,i,j)
           ENDDO
           somcmud=0.0_rsh
           DO iv=imud1,imud2
             somcmud=somcmud+cv_sed(iv,k,i,j)
           ENDDO
           !hinder : entravement du sable/gravier (sans dimension) entre 0 et 1
           ! ATTENTION : hypothese de ros homogene quelque soit le type de sediment
           ! entrave par le trop de vase qui empeche la segregation sable/vase (1er terme)
           ! entrave par le trop de sablequi ne permet pas le passage entre les grains (2ieme terme)
           hinder(k)=(MAX(0.0_rsh,MIN(1.0_rsh-somcmud/(1-(somcsan+somcgrav)/ros(1))/csegreg,   &
                   1.0_rsh-(somcsan+somcgrav)/csandseg)))**4.65_rsh
                   
           DO k=ksmin+1,ksmax-1

              sed_rate(k)=-rowinv*(permeab(k+1)+permeab(k))                               &
                  *MIN(0.0_rsh,loadograv(k+1)-loadograv(k))*stateconsol(k)       &
                  /(dzs(k+1,i,j)+dzs(k,i,j))
              somcgrav=0.0_rsh
              DO iv=igrav1,igrav2
                somcgrav=somcgrav+cv_sed(iv,k,i,j)
              ENDDO
              somcsan=0.0_rsh
              DO iv=isand1,isand2
                somcsan=somcsan+cv_sed(iv,k,i,j)
              ENDDO
              somcmud=0.0_rsh
              DO iv=imud1,imud2
                somcmud=somcmud+cv_sed(iv,k,i,j)
              ENDDO
              hinder(k)=(MAX(0.0_rsh,MIN(1.0_rsh-somcmud/(1-(somcsan+somcgrav)/ros(1))/csegreg,   &
                   1.0_rsh-(somcsan+somcgrav)/csandseg)))**4.65_rsh
           ENDDO
        
           k=ksmax
           dtsdzs=REAL(dtiter,rsh)/dzs(k,i,j)
           !schema implicite decentre amont -cv_sed(k+1) connu juste avant donc implicite
           DO iv=igrav1,isand2
              cv_sed(iv,ksmax,i,j)=cv_sed(iv,ksmax,i,j)                               &
                                 /(1.0_rsh+dtsdzs*MAX(sed_rate(ksmax-1),hinder(ksmax-1)*ws_sand(iv)))
           ENDDO
           DO iv=imud1,nvp
              cv_sed(iv,ksmax,i,j)=cv_sed(iv,ksmax,i,j)/(1.0_rsh+dtsdzs*sed_rate(ksmax-1))
           ENDDO
           DO k=ksmax-1,ksmin+1,-1
              dtsdzs=REAL(dtiter,rsh)/dzs(k,i,j)
              DO iv=igrav1,isand2
                cv_sed(iv,k,i,j)=(cv_sed(iv,k,i,j)+dtsdzs*MAX(sed_rate(k),ws_sand(iv)*               &
                               hinder(k))*cv_sed(iv,k+1,i,j))                        &
                               /(1.0_rsh+dtsdzs*MAX(sed_rate(k-1),hinder(k-1)*ws_sand(iv)))
              ENDDO
              DO iv=imud1,nvp
                cv_sed(iv,k,i,j)=(cv_sed(iv,k,i,j)+dtsdzs*sed_rate(k)*cv_sed(iv,k+1,i,j))/(1.0_rsh+dtsdzs*sed_rate(k-1))
              ENDDO
           ENDDO
           k=ksmin
           dtsdzs=REAL(dtiter,rsh)/dzs(k,i,j)
           DO iv=igrav1,isand2
              cv_sed(iv,k,i,j)=(cv_sed(iv,k,i,j)+dtsdzs*MAX(sed_rate(k),ws_sand(iv)*hinder(k))*cv_sed(iv,k+1,i,j))
           ENDDO
           DO iv=imud1,nvp
              cv_sed(iv,k,i,j)=(cv_sed(iv,k,i,j)+dtsdzs*sed_rate(k)*cv_sed(iv,k+1,i,j))
           ENDDO
          
         ELSE
            dtiter=dtrest
         ENDIF 
         ! end of l_consolid

         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         !!!!!!!!!!!  BIOTURBATION PARTICULAIRE    !!!!!!!!!!!!!
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                                      
!  bioturbation particulaire A REVOIR
#ifdef key_sed_bioturbpart
          ! schema explicite par pas fractionnaire pour la diffusion : bioturbation par melange diffusif
         IF(l_bioturb) THEN
         
            ! calcul du coef de biodiffusion a reactualiser si nombre de couches a change pendant diter
            CALL sed_MUSTANG_coefbioturb(i,j,difbio) ! A REprogrammer pour differencier particulaires et dissoutes

            k=ksmax
            dtsdzs=REAL(dtiter,rsh)/dzs(k,i,j)
            dtsdzsmin=REAL(dtiter,rsh)/MIN(dzs(k,i,j),dzs(k-1,i,j))
            IF(dtsdzsmin>0.5_rsh*dzs(k,i,j)/(difbio(k-1,imud1)+1.e-10_rsh) .AND. difbio(k-1,imud1) > 1.e-13) THEN
            ! melange instantane
              DO iv=isand1,nvp
                cv_sed(iv,k,i,j)= (cv_sed(iv,k,i,j)*dzs(k,i,j)+cv_sed(iv,k-1,i,j)*dzs(k-1,i,j)) &
                                   /(dzs(k,i,j)+dzs(k-1,i,j))
                cv_sed(iv,k-1,i,j)= cv_sed(iv,k,i,j) 
                cvsednew(iv,k)= cv_sed(iv,k,i,j)                                                     
              ENDDO
            ELSE
              DO iv=isand1,nvp
                cvsednew(iv,k)= cv_sed(iv,k,i,j) + &
                                dtsdzs*difbio(k-1,iv)*dzsiin(k-1)*(cv_sed(iv,k-1,i,j)-cv_sed(iv,k,i,j))                                     
              ENDDO
            ENDIF
            DO k=ksmax-1,ksmin+1,-1
              dtsdzs=REAL(dtiter,rsh)/dzs(k,i,j)
              dtsdzsmin=REAL(dtiter,rsh)/MIN(dzs(k,i,j),dzs(k-1,i,j),dzs(k+1,i,j))
              IF(dtsdzsmin>0.5_rsh*dzs(k,i,j)/(difbio(k,imud1)+1.e-10_rsh).AND. difbio(k,imud1) > 1.e-13) THEN
               ! melange instantane
               DO iv=isand1,nvp
                  cv_sed(iv,k,i,j)= (cv_sed(iv,k,i,j)*dzs(k,i,j)+cv_sed(iv,k-1,i,j)*dzs(k-1,i,j)+ &
                          cv_sed(iv,k+1,i,j)*dzs(k+1,i,j))/(dzs(k,i,j)+dzs(k-1,i,j)+dzs(k+1,i,j))
                  cv_sed(iv,k-1,i,j)= cv_sed(iv,k,i,j) 
                  cv_sed(iv,k+1,i,j)= cv_sed(iv,k,i,j) 
                  cvsednew(iv,k)= cv_sed(iv,k,i,j)                                                     
                ENDDO
              ELSE
                DO iv=isand1,nvp
                  cvsednew(iv,k)= cv_sed(iv,k,i,j) + dtsdzs*(   &
                                difbio(k,iv)*dzsiin(k)*cv_sed(iv,k+1,i,j)  &
                              -(difbio(k,iv)*dzsiin(k)+difbio(k-1,iv)*dzsiin(k-1))*cv_sed(iv,k,i,j) &
                              + difbio(k-1,iv)*dzsiin(k-1)*cv_sed(iv,k-1,i,j))                                     
                ENDDO              
              ENDIF
            ENDDO              
            k=ksmin
            dtsdzs=REAL(dtiter,rsh)/dzs(k,i,j)    
            dtsdzsmin=REAL(dtiter,rsh)/MIN(dzs(k,i,j),dzs(k+1,i,j))
            IF(dtsdzsmin>0.5_rsh*dzs(k,i,j)/(difbio(k,imud1)+1.e-10_rsh).AND. difbio(k,imud1) > 1.e-13) THEN
            ! melange instantane
              DO iv=isand1,nvp
                cv_sed(iv,k,i,j)= (cv_sed(iv,k,i,j)*dzs(k,i,j)+cv_sed(iv,k+1,i,j)*dzs(k+1,i,j))&
                                   /(dzs(k,i,j)+dzs(k+1,i,j))
                cv_sed(iv,k+1,i,j)= cv_sed(iv,k,i,j) 
                cvsednew(iv,k)= cv_sed(iv,k,i,j)                                                     
              ENDDO
            ELSE
              DO iv=isand1,nvp
               cvsednew(iv,k)= cv_sed(iv,k,i,j) + &
                                dtsdzs*difbio(k,iv)*dzsiin(k)*(cv_sed(iv,k+1,i,j)-cv_sed(iv,k,i,j))                                     
              ENDDO              
            ENDIF
            DO k=ksmin,ksmax
              cv_sed(isand1:nvp,k,i,j)=cvsednew(isand1:nvp,k)
            ENDDO              
         ENDIF   
          ! end of l_bioturb
#endif

          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          !!!!!!!!!!!  MISE A JOUR DES CONCENTRATIONS PARTICULAIRES   !!!!!!!!!!!
          !!!!!!!!!!!       APRES CONSOLIDATION ET/OU BIOTURBATION       !!!!!!!!!!!
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef key_sed_bioturbpart
         IF(l_consolid .OR. l_bioturb) THEN
#else
         IF(l_consolid) THEN
#endif
         ! computing new concentrations/porosities :
          DO k=ksmin,ksmax
            c_sedtot(k,i,j)=0.0_rsh
            cvolp=0.0_rsh
            DO iv=1,nvpc
              c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
              cvolp=cvolp+cv_sed(iv,k,i,j)/ros(iv)
            ENDDO
            poro(k,i,j)=1.0_rsh-cvolp
            
            ! if(porosity too small (<32%), it is increased until (1.0_rsh-cvolmaxmel):
            !  dissolved variables concentrations are inchanged
            !------------------------------------------------------------------------
            IF(poro(k,i,j).LT.1.0_rsh-cvolmaxmel-0.01_rsh)THEN
              poro(k,i,j)=1.0_rsh-cvolmaxmel
            !!! ATTENTION : on suppose ros identique quelque soit les particules
              csednew=ros(1)*cvolmaxmel
              rapcsed=csednew/c_sedtot(k,i,j)
              DO iv=1,nvp
                cv_sed(iv,k,i,j)=cv_sed(iv,k,i,j)*rapcsed
              ENDDO
              c_sedtot(k,i,j)=csednew
              dzs(k,i,j)=dzs(k,i,j)/rapcsed
            ENDIF 
            
          ENDDO
        
          ! again calcul pour ksmax car on fait varier dzs dans la couche de surface
          !  on a perdu de la matiere mais c_sedtot reste constant => dzs change
          k=ksmax
          dzsa=dzs(k,i,j)
          dzs(k,i,j)=c_sedtot(k,i,j)*dzsa/cseda(k)
          cvolp=0.0_rsh
          c_sedtot(k,i,j)=0.0_rsh
          DO iv=1,nvp
            cv_sed(iv,k,i,j)=cv_sed(iv,k,i,j)*dzsa/dzs(k,i,j)
            c_sedtot(k,i,j)=c_sedtot(k,i,j)+typart(iv)*cv_sed(iv,k,i,j)
            cvolp=cvolp+typart(iv)*cv_sed(iv,k,i,j)/ros(iv)
          ENDDO
          ! dans l hypothese ou ros est homogene quelque soit le type de sediment, 
          !  poro devrait rester inchange aussi
          poro(k,i,j)=1.0_rsh-cvolp
          
          ! mis en commentaire car inutile
          ! if(porosity too small (<32%), it is increased until (1.0_rsh-cvolmaxmel):
          !IF(poro(k,i,j) < 1.0_rsh-cvolmaxmel-0.01_rsh)THEN
          !   poro(k,i,j)=1.0_rsh-cvolmaxmel
          !   !!! ATTENTION : on suppose ros identique quelque soit les particules
          !   csednew=ros(1)*cvolmaxmel
          !   rapcsed=csednew/c_sedtot(k,i,j)
          !   DO iv=1,nvp
          !     cv_sed(iv,k,i,j)=cv_sed(iv,k,i,j)*rapcsed
          !   ENDDO
          !   c_sedtot(k,i,j)=csednew
          !   dzs(k,i,j)=dzs(k,i,j)/rapcsed
          !ENDIF 
 
          ksmaxa=ksmax
          DO k=ksmaxa,ksmin+1,-1
                   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                   !!!  FUSION DE 2 COUCHES SI TROP FINE OU TROP DE PEU DE MATIERE   !!!!!
                   !!! (HYP : csed_tot de la couche de dessous inchange)             !!!!!
                   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            IF(c_sedtot(k,i,j) .LE. csedmin .OR. dzs(k,i,j) .LE. dzsmin)THEN

              ! etat de la couche avant fusion
              dzsa=dzs(k,i,j)
              dzsam1=dzs(k-1,i,j)
              poroam1=poro(k-1,i,j)
              poroa=poro(k,i,j)
              
              ! csed_tot(k-1 Cste : l epaisseur de la couche k-1 n est donc pas la somme des deux epaisseurs.
              !dzs(k-1,i,j)=(dzsam1*c_sedtot(k-1,i,j)+c_sedtot(k,i,j)*dzsa)   &
              !             /c_sedtot(k-1,i,j)

              ! ou autre hypothese pour conserver eau et variables dissoutes
              ! csed_tot non cste mais volume d eau et quantite de particules restent constantes
              dzs(k-1,i,j)=dzsam1+dzsa
              dzsisu=1.0_rsh/dzs(k-1,i,j)
              
              c_sedtot(k-1,i,j)=0.0_rsh
              cvolp=0.0_rsh
              DO iv=1,nvp
                cv_sed(iv,k-1,i,j)=(cv_sed(iv,k,i,j)*dzsa+                           &
                                    cv_sed(iv,k-1,i,j)*dzsam1)*dzsisu
                c_sedtot(k-1,i,j)=c_sedtot(k-1,i,j)+cv_sed(iv,k-1,i,j)*typart(iv)
                cvolp=cvolp+cv_sed(iv,k-1,i,j)*typart(iv)/ros(iv)
              ENDDO
              poro(k-1,i,j)=1.0_rsh-cvolp
 
              volak=poroa*dzsa
              volakm1=poroam1*dzsam1
              vola=volak+volakm1
              
              ! mis en commentaire car inutile
              ! if(porosity too small (<32%), it is increased until (1.0_rsh-cvolmaxmel):
              !IF(poro(k-1,i,j) < 1.0_rsh-cvolmaxmel-0.01_rsh)THEN
              !  poro(k-1,i,j)=1.0_rsh-cvolmaxmel
            !!!! ATTENTION : on suppose ros identique quelque soit les particules
              !  csednew=ros(1)*cvolmaxmel
              !  rapcsed=csednew/c_sedtot(k-1,i,j)
              !  DO iv=1,nvp
              !    cv_sed(iv,k-1,i,j)=cv_sed(iv,k-1,i,j)*rapcsed
              !  ENDDO
              !  c_sedtot(k-1,i,j)=csednew
              !  dzs(k-1,i,j)=dzs(k-1,i,j)/rapcsed
              !ENDIF   
             
             ! calcul des flux d eau et de matiere dissoute dus a la fusion des deux couches
             ! dvol_fusion  en m3/m2; 
              !dvol=vola-poro(k-1,i,j)*dzs(k-1,i,j)
              !dvol_fusion=dvol +dvol_fusion            

#if ! defined key_nodiss_insed
              ! calcul des nouvelles concentrations de variables dissoutes dans l eau interstitielle
              ! dans ancien volume d eau melange dans la couche k-1
              vola_inv=1._rsh/vola

              DO iv=-1,0
                cv_sed(iv,k-1,i,j)=vola_inv*(cv_sed(iv,k-1,i,j)*volakm1 + &
                                   cv_sed(iv,k,i,j)*volak)
              !  flu_fusion(iv)=flu_fusion(iv)+cv_sed(iv,k-1,i,j)*dvol
              ENDDO
#if ! defined key_sedim_partial
              DO iv=nvp+1,nv_use2
                cv_sed(iv,k-1,i,j)=vola_inv*(cv_sed(iv,k-1,i,j)*volakm1 + &
                                   cv_sed(iv,k,i,j)*volak)
               ! flu_fusion(iv)=flu_fusion(iv)+cv_sed(iv,k-1,i,j)*dvol
              ENDDO
#endif
#endif

              volpwak(k-1)=volpwak(k-1)+volpwak(k)
              
              DO kp=k,ksmax-1
                dzs(kp,i,j)=dzs(kp+1,i,j)
                volpwak(kp)=volpwak(kp+1)
                DO iv=-1,nv_use2
                  cv_sed(iv,kp,i,j)=cv_sed(iv,kp+1,i,j)
                ENDDO
                c_sedtot(kp,i,j)=c_sedtot(kp+1,i,j)
                poro(kp,i,j)=poro(kp+1,i,j)
              ENDDO
              dzs(ksmax:ksdmax,i,j)=0.0_rsh
              c_sedtot(ksmax:ksdmax,i,j)=-valmanq
              cv_sed(-1:nv_use2,ksmax:ksdmax,i,j)=-valmanq
              poro(ksmax:ksdmax,i,j)=0.0_rsh
              ksmax=ksmax-1
            ENDIF  
            ! fin fusion

          ENDDO        
             
          ! updating ksma (ksmax can be modified in the routine)
          ksma(i,j)=ksmax
 
          ! plusieurs pas possibles de consolidation 
          dtrest=dtrest-dtiter
          IF(dtrest.GE.1.0_rlg)GOTO 100

          ! Exchange flux at the water/sediment interface
          ! and intertitial water vertical rates
          ! ---------------------------------------------
          !  fluconsol is flux resulting from consolidation, >0 upwards, centree a l interface superieure
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
          !!!! interstitial water vertical rates  : winter (m3/m2/s) k :between k and k+1
          !!!!flux total de substance dissoute due à la fusion : M/m2/s
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
          winter(ksmin-1)=0.0_rsh
          volpwnv_tot=0.0_rsh
          DO k=ksmin,ksmax
           volpwk=poro(k,i,j)*dzs(k,i,j)
           volpwnv_tot=volpwnv_tot+volpwk
           winter(k)=winter(k-1)+(volpwak(k)-volpwk)*dt_sed_inv
          ENDDO
         
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
          !!!! water fluxes between water and sediment due to consolidation and fusion (m3/s)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
#if ! defined key_nofluxwat_IWS
          phieau_s2w_consol(i,j)=(volpwa_tot-volpwnv_tot)*surf(i,j)
#endif
          !  flu_dyninsed est le flux par tassement a l interface, >0 vers le haut (M/m2)
          flu_dyninsed(-1,i,j)=winter(ksmax)*cv_sed(-1,ksmax,i,j)*REAL(dt_sed_eff,rsh)
          flu_dyninsed(0,i,j)=winter(ksmax)*cv_sed(0,ksmax,i,j)*REAL(dt_sed_eff,rsh) 
          flu_dyninsed(nvp+1:nv_use2,i,j)= &
                        winter(ksmax)*cv_sed(nvp+1:nv_use2,ksmax,i,j)*REAL(dt_sed_eff,rsh) 
                
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                      
          !!!!!!!  FIN consolidation -                               !!!
          !! actualisation epaisseurs et porosites intermediaires    !!!
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                      
           
          poroin(ksma(i,j))=poro(ksma(i,j),i,j)
          DO k=ksmi(i,j),ksma(i,j)-1
              dzsiin(k)=2.0_rsh/(dzs(k,i,j)+dzs(k+1,i,j))
              poroin(k)=0.5_rsh*(poro(k+1,i,j)+poro(k,i,j))          
          ENDDO
        ENDIF  ! fin test l_consolid or l_bioturb
        
#if ! defined key_nodiss_insed
        
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !!!!!!!!!!!  DIFFUSION -BIOTURBATION DES VARIABLES DISSOUTES       !!!!
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#if ! defined key_toseelater_bioturb
         IF(l_bioturb) THEN
            ! calcul du coef de biodiffusion a reactualiser si nombre de couches a change 
            CALL sed_MUSTANG_coefbioturb(i,j,difbio) ! A Reprogrammer pour differencier particulaires et dissoutes
         ENDIF
#endif

         hcrit=MAX(0.0_rsh,htot(i,j)-h0fond)
         DO iv=-1,nv_use2
          !only temperature and Salinity  (-1 and 0) and if NOT key_sedim_partial :  dissolved variables also
           IF(iv<1 .OR. iv>nvp) THEN
         
              ! calcul des vitesses de diffusion/biodiffusion dans le sediment
              !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            disvi(:)=0.0_rsh
            IF(l_diffused) THEN

              k=ksmax
              IF(choice_flxdiss_diffsed == 1) THEN
                disvi(k)=2.0_rsh*xdifsi1/(dzs(k,i,j)+epn_bottom_MUSTANG(i,j))
              ELSE IF(choice_flxdiss_diffsed == 2) THEN
                disvi(k)=xdifsi1/(dzs(k,i,j)*0.5_rsh+epdifi)
              ELSE
               ! dans la couche limite : formulation de Boudreau (1997)
               ! vitesse de diffusion= D/deltae=beta (page 182) eq 5.55 a 5.56)
               ! et beta =0.0889 u* Sc^(-0.704)
               ! formule de Shaw and Hanratty (1977) eq 5.58 in Boudreau 
               ! Sc nombre de Schmidt = viscosite cinematique/diffusivite massique
               ! de l ordre de 1000 mais peut varier ??   A REVOIR
                disvi(k)=0.0889*ustarbot(i,j)*(1000._rsh)**(-0.704)
              ENDIF
              disvi(k)=disvi(k)*hcrit/(hcrit+epsilon)
              DO k=ksmin,ksmax-1
               ! formulation en fonction de la tortuosite Dsed=Dpure/Tortuosite^2
               ! tortuosite fontion de la porosite (eq 4.120)
               ! in Boudreau 1997 p 132 ! totuosity^2=1-log(poro^2)
               !disvi(k)=xdifs1*dzsiin(k,i,j)/(1-LOG(poro(k,i,j)**2))
                disvi(k)=(xdifs1/(1.0_rsh-2.0_rsh*LOG(poro(k,i,j)))+difbio(k,iv))*dzsiin(k)
              ENDDO
            ENDIF  !(diffused)
!
              !resolution equation d advection/diffusion (advection par consolidation)
!
            IF(ksmax > ksmin)THEN
!
                !cas ou il y a plusieurs couches:
                !--------------------------------
                !coefficients de la matrice
!
                k=ksmin
                dzssdt=dzs(k,i,j)*dt_sed_inv
                aa(k)=0.0_rsh
                bb(k)=poro(k,i,j)*dzssdt+disvi(k)*poroin(k)+winter(k)*fexcs
                cc(k)=-disvi(k)*poroin(k)+winter(k)*cexcs
                dd(k)=cv_sed(iv,k,i,j)*volpwak(k)*dt_sed_inv
                DO k=ksmin+1,ksmax-1
                 dzssdt=dzs(k,i,j)*dt_sed_inv
                 aa(k)=-disvi(k-1)*poroin(k-1)-winter(k-1)*fexcs
                 bb(k)=poro(k,i,j)*dzssdt+disvi(k-1)*poroin(k-1)+disvi(k)*poroin(k) &
                                +winter(k)*fexcs-winter(k-1)*cexcs
                 cc(k)=-disvi(k)*poroin(k)+winter(k)*cexcs
                 dd(k)=cv_sed(iv,k,i,j)*volpwak(k)*dt_sed_inv
                ENDDO
                k=ksmax
                aa(k)=-disvi(k-1)*poroin(k-1)-winter(k-1)*fexcs
                bb(k)=dzs(k,i,j)*poro(k,i,j)*dt_sed_inv+disvi(k-1)*poroin(k-1)+disvi(k)+winter(k)-winter(k-1)*cexcs
                cc(k)=0.0_rsh
                IF (iv==-1) THEN
                  fludif(iv,i,j)=disvi(k)*temp_bottom_MUSTANG(i,j)
                ELSE IF (iv== 0) THEN
                  fludif(iv,i,j)=disvi(k)*sal_bottom_MUSTANG(i,j)
                ELSE
                  fludif(iv,i,j)=disvi(k)*cw_bottom_MUSTANG(iv,i,j)
                ENDIF
                dd(k)=cv_sed(iv,k,i,j)*volpwak(k)*dt_sed_inv+fludif(iv,i,j)
!
                !resolution du systeme tridiagonal
!
                dd(1)=dd(1)/bb(1)
                DO k=ksmin+1,ksmax
                  bb(k)=bb(k)-aa(k)*cc(k-1)/bb(k-1)
                  dd(k)=(dd(k)-aa(k)*dd(k-1))/bb(k)
                ENDDO
                cv_sed(iv,ksmax,i,j)=dd(ksmax)
                DO k=ksmax-1,ksmin,-1
                  cv_sed(iv,k,i,j)=dd(k)-cc(k)*cv_sed(iv,k+1,i,j)/bb(k)
                ENDDO
!
            ELSE
!
               !cas ou il n y a qu une couche (ksmax=ksmin):
               !--------------------------------------------       
                k=ksmax        
                IF (iv==-1) THEN
                  fludif(iv,i,j)=disvi(k)*temp_bottom_MUSTANG(i,j)
                ELSE IF (iv== 0) THEN
                  fludif(iv,i,j)=disvi(k)*sal_bottom_MUSTANG(i,j)
                ELSE
                  fludif(iv,i,j)=disvi(k)*cw_bottom_MUSTANG(iv,i,j)
                ENDIF        
                cv_sed(iv,k,i,j)=  &
                    (cv_sed(iv,k,i,j)*volpwak(k)+fludif(iv,i,j)*REAL(dt_sed_eff,rsh))           &
                   /(dzs(k,i,j)*poro(k,i,j)+(disvi(k)+winter(k))*REAL(dt_sed_eff,rsh))
            ENDIF
!
              !  flux d echange a l interface eau/sediment:
              !  ------------------------------------------
              !  fludif est le flux par diffusion a l interface, >0 vers le bas  ==> M/m2
              IF(l_diffused)fludif(iv,i,j)=  &
                           (fludif(iv,i,j)-disvi(ksmax)*cv_sed(iv,ksmax,i,j))*REAL(dt_sed_eff,rsh)
           

           ENDIF                     ! test sur iv ( var dissoutes seules)
         ENDDO                         ! fin boucle sur iv

#endif
        ENDIF    ! test si existence sediment
    
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       !!!!!!!!!!!! FIN DU PROCESSUS DE DYNAMIQUE DANS LE SEDIMENT      !!
       !!!!!!!!!!!!  CALCUL DES FLUX D EAU AVEC TERME CONSTANT          !!
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! si la cellule est mise a sec  ==> on cumule les apports d  eau et les flux de substances
       !  dissoutes pour ne pas perdre de la matiere
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         IF(htot(i,j) .LE. h0fond) THEN  
            ! il n y a pas d eau
#if ! defined key_nofluxwat_IWS
            ! phieau_s2w_drycell en m3  
            phieau_s2w_drycell(i,j)=phieau_s2w_drycell(i,j)+REAL(phieau_s2w_consol(i,j),rlg)
#endif
#if ! defined key_nodiss_insed
            !fluconsol_drycell : cumul => M/m2
            fluconsol_drycell(-1:0,i,j)=fluconsol_drycell(-1:0,i,j) +flu_dyninsed(-1:0,i,j)
            fluconsol_drycell(nvp+1:nv_use2,i,j)=fluconsol_drycell(nvp+1:nv_use2,i,j)  &
                                               +flu_dyninsed(nvp+1:nv_use2,i,j)
            fluconsol(-1:0,i,j)=0.0_rsh
            fluconsol(nvp+1:nv_use2,i,j)=0.0_rsh                                    
#endif 
        ELSE
#if ! defined key_nofluxwat_IWS
           ! phieau_s2w en m3  
            phieau_s2w(i,j)=phieau_s2w(i,j)+phieau_s2w_drycell(i,j) + &
                                          REAL(phieau_s2w_consol(i,j),rlg)
            phieau_s2w_drycell(i,j)=0.0_rlg
#endif
#if ! defined key_nodiss_insed
            ! fluconsol  => M/m2 sera distribue dans flx_s2w (M/m2/s)  
            fluconsol(-1:0,i,j)=flu_dyninsed(-1:0,i,j) +fluconsol_drycell(-1:0,i,j)
            fluconsol(nvp+1:nv_use2,i,j)=flu_dyninsed(nvp+1:nv_use2,i,j)  &
                                       +fluconsol_drycell(nvp+1:nv_use2,i,j)
            fluconsol_drycell(-1:0,i,j)=0.0_rsh
            fluconsol_drycell(nvp+1:nv_use2,i,j)=0.0_rsh
#endif 
        ENDIF
  
      ENDDO
    ENDDO
!$OMP END DO
!$OMP SINGLE
     t_sed=t+dt_sed_cor
!$OMP END SINGLE
    ELSE
!$OMP SINGLE
            fluconsol(:,:,:)=0.0_rsh
            fludif(:,:,:)=0.0_rsh
!$OMP END SINGLE
    ENDIF
   
  END SUBROUTINE sed_MUSTANG_consol_diff_bioturb


!!===========================================================================================
!
      SUBROUTINE sed_MUSTANG_coefbioturb(i,j,difbio)
!
!  !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_coefbioturb  ***
   !&E
   !&E ** Purpose :     Bioturbation in sediment
   !&E
   !&E ** Description : use arguments and common variable 
   !&E  arguments IN : 
   !&E         ksmin,ksmax
   !&E         cs
   !&E          
   !&E  arguments OUT:
   !&E            
   !&E  need to be know by hydrodynamic code:
   !&E           ksdmin,ksdmax
   !&E          
   !&E  need to be know by code treated substance:
   !&E         nv_adv ,nv_tot                  
   !&E         igrav1,isand2
   !&E       
   !&E  use module SEDIM1DV variables  :
   !&E         dbiotu0, dbiotum
   !&E         frmud_db_min
   !&E         frmud_db_max
   !&E         xbioturbk xbioturbmax
   !&E          
   !&E        
   !&E
   !&E ** Called by :  consoldif
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-02  new version from SiAM"sande-vase" by P. Le Hir
   !&E       !  2014-05  (F. Grasso) modifications based on "Grasso et al. (Ocean Dynynamics, 2014)"
   !&E
! **********************************************************************
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   INTEGER, INTENT(IN)            :: i,j   
   REAL(KIND=rsh),DIMENSION(ksdmin:ksdmax,-1:nv_adv),INTENT(OUT)  ::difbio                   
   
   !! * Local declarations
   INTEGER        :: iv,k,ksmib,ksmin,ksmax               
   REAL(KIND=rsh) :: dsurf,somcgravsand,somcmud,coefdifbio,frmud,hsedloc

   !!---------------------------------------------------------------------------
   !! * Executable part

        ksmax=ksma(i,j)
        ksmin=ksmi(i,j)
        difbio(ksdmin:ksdmax,-1:nv_adv)=0.0_rsh

        IF (ksmax.GE.ksmin)THEN
          hsedloc=0.0_rsh
          DO k=ksmin,ksmax
            hsedloc=hsedloc+dzs(k,i,j)
            difbio(k,:)=0.0_rsh
          ENDDO   
          IF(hsedloc < dbiotu0) THEN
            ksmib=ksmax
          ELSE
            dsurf=-dzs(ksmax,i,j) ! evaluation de Db(k) a la limite superieure
            DO k=ksmax,ksmin,-1
              dsurf=dsurf+dzs(k,i,j)
              IF(dsurf .GE. dbiotu0) THEN
                ksmib=k
                GOTO 11
              ENDIF
              !difbio dependent on mud fraction
              somcgravsand=0.0_rsh
              DO iv=igrav1,isand2
                somcgravsand=somcgravsand+cv_sed(iv,k,i,j)
              ENDDO
              somcmud=0.0_rsh
              DO iv=imud1,imud2
                somcmud=somcmud+cv_sed(iv,k,i,j)
              ENDDO
              frmud=somcmud/(somcgravsand+somcmud)
              IF(frmud < frmud_db_min) THEN
                coefdifbio=0.0_rsh
              ELSE IF(frmud > frmud_db_max) THEN
                coefdifbio=1.0_rsh
              ELSE
                coefdifbio=(frmud-frmud_db_min)/(frmud_db_max-frmud_db_min)
              ENDIF
              IF(dsurf .GE. dbiotum) THEN
               ! difbio=0 for gravels
               ! difbio same for sand, muds and dissolved variables here - to be changed furthermore ?
               coefdifbio=coefdifbio*EXP(-xbioturbk*(dsurf-dbiotum)/(dbiotu0-dbiotum))
              ENDIF
              difbio(k,isand1:nv_adv)=coefdifbio*xbioturbmax
              difbio(k,-1:0)=coefdifbio*xbioturbmax
            ENDDO    
            ksmib=1
          ENDIF
 11       DO k=ksmin,ksmib
            difbio(k,:)=0.0_rsh
          ENDDO
        ENDIF
 
    
  END SUBROUTINE sed_MUSTANG_coefbioturb
   
!!==============================================================================
  SUBROUTINE sed_constitutivrel(i,j,stateconsol,permeab,sigmapsg)
! 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_constitutivrel  ***
   !&E
   !&E ** Purpose : applies constitutive relationships for permeability, 
   !&E              effective stess
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  consoldif
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-02  new version from SiAM"sande-vase" by P. Le Hir
   !&E       !  2014-05  (F. Grasso) modifications based on "Grasso et al. (Ocean Dynynamics, 2014)"
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used
#ifdef key_MARS
   USE parameters,   ONLY : grav
#endif

   !! * Arguments
   INTEGER,INTENT(IN) :: i,j
   REAL(KIND=rsh),DIMENSION(ksdmin:ksdmax),INTENT(out) :: stateconsol,permeab,sigmapsg
   
  !! * Local declarations
   INTEGER                        :: k,iv,icond
   REAL(KIND=rsh),DIMENSION(nvpc) :: cvol
   REAL(KIND=rsh)                 :: cvolsed,cvolsangrv,fracvolcum,diarepres,       &
                                     voidratio,diamstarrep,wsoulsby,somcmud,crelmud,cvolrelmud, &
                                     permeab_e


   !! * executable statement
   !!   --------------------
   sigmapsg(:)=0.0_rsh
   stateconsol(:)=0.0_rsh

   DO k=ksmi(i,j),ksma(i,j)
     cvol(:)=0.0_rsh
     cvolsed=0.0_rsh
     DO iv=igrav1,isand2
       cvol(iv)=cv_sed(iv,k,i,j)/ros(iv)
       cvolsed=cvolsed+cvol(iv)
     ENDDO
     cvolsangrv=cvolsed
     somcmud=0.0_rsh
     DO iv=isand2+1,nvpc
       cvol(iv)=cv_sed(iv,k,i,j)/ros(iv)
       cvolsed=cvolsed+cvol(iv)
       somcmud=somcmud+cv_sed(iv,k,i,j)
     ENDDO
       
     !  permeability with voidratio:
     voidratio=poro(k,i,j)/(1.0000001_rsh-poro(k,i,j))
     voidratio=MAX(voidratio,0.0_rsh)
     permeab_e=2.0e-9_rsh*voidratio**3.7_rsh   ! best param for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)

     !  permeability with cvolrelmud (relative volumic concentration of mud, cvolrelmud=cvolmud/(1-cvolsand), cf. Merckelbach & Kranenburg, 2004)
     crelmud=somcmud/(1.0_rsh-cvolsangrv)
     cvolrelmud=crelmud/ros(1)
     !!!FG(06/03/2015) pour corriger le cas ou cvolrelmud est trop faible (Cvolrelmud = 0.01 -> Crelmud = 26 g/l -> permeab = 4 m/s)
     IF(cvolrelmud .LE. 0.01) THEN
       permeab(k)=permeab_e
     ELSE
       permeab(k)=xperm1*cvolrelmud**xperm2
     ENDIF

!     !  minimalisation of permeability
     permeab(k)=MIN(permeab(k),permeab_e,(1.0_rsh+voidratio)*0.0001_rsh/1.65_rsh)   ! equivalent permeability for ws=0.1 mm/s

     !  effective stress :
     icond=1
     DO iv=igrav1,isand2
       icond=icond*MAX(0.0_rsh,100.0_rsh*(cvolmaxsort-cvol(iv)-.02_rsh))
     ENDDO
     IF(cvolsed.LE.cvolmaxmel.AND.icond.NE.0)THEN
       sigmapsg(k)=xsigma1/grav*cvolrelmud**xsigma2 !FG(04/09/2013) to use Merckelback & Kranenburg s (2004) formulation
       stateconsol(k)=1.0_rsh
     ELSE
       !    in this case, the effective stress is maximum, 
       !    but porewater overpressure is zero :
       stateconsol(k)=0.0_rsh
     ENDIF
   ENDDO

  END SUBROUTINE sed_constitutivrel
                                
 !!==============================================================================


 
#endif

END MODULE sed_MUSTANG

