!------------------------------------------------------------------------------
MODULE initMUSTANG
!------------------------------------------------------------------------------

#include "cppdefs.h"

#ifdef MUSTANG

!&E============================================================================
!&E                   ***  MODULE  initMUSTANG  ***
!&E
!&E ** Purpose : concerns all subroutines related to module MUSTANG
!&E              initialization
!&E 
!&E ** Description :
!&E     subroutine MUSTANG_init           ! initialize the 
!&E                                         parameters of sediment module
!&E                                       ! and the sediment layers
!&E     subroutine MUSTANG_readnml        ! read namelists
!&E     subroutine MUSTANG_param_log      ! write parameters for information
!&E     subroutine MUSTANG_compatibility  ! check compatibility 
!&E                                         between the various parameters
!&E     subroutine MUSTANG_sedinit        ! initialize the sediment if not 
!&E                                         from file
!&E     subroutine MUSTANG_init_hsed      ! initialize the sediment thickness 
!&E                                         function sediment parameters 
!&E                                         (MUSTANG_mixsed)
!&E     subroutine MUSTANG_init_param     ! initialize the settling velocity 
!&E                                         of sand variables and erosion 
!&E                                         parameters
!&E     subroutine MUSTANG_init_output    ! initialize dimensions of output 
!&E                                         tables
!&E     subroutine MUSTANG_morphoinit      ! initialize 
!&E
!&E============================================================================
    !! * Modules used
#include "coupler_define_MUSTANG.h"

    USE comMUSTANG
    USE sed_MUSTANG_HOST,  ONLY : sedinit_fromfile
    USE sed_MUSTANG,  ONLY : MUSTANG_E0sand
    !USE sed_MUSTANG_HOST,  ONLY : sed_exchange_hxe_HOST
    USE comsubstance
    USE module_substance
#ifdef key_MUSTANG_flocmod
    USE flocmod, ONLY : flocmod_alloc, flocmod_init
    USE flocmod, ONLY : f_ws, f_diam, f_vol, f_rho, f_mass
#endif

    IMPLICIT NONE

    !! * Accessibility
    PUBLIC MUSTANG_init

    PRIVATE

    !! * Local declarations
    ! declaration of namelists, variables described in comMUSTANG
    namelist /namsedim_init/ l_repsed, filrepsed, l_unised, fileinised,       &
                             date_start_morpho, date_start_dyninsed,          &
                             hseduni, cseduni, ksmiuni, ksmauni,              &
                             sini_sed, tini_sed,                              &
                             l_init_hsed, csed_mud_ini,                       &
                             l_initsed_vardiss, poro_mud_ini

    namelist /namsedim_layer/ l_dzsminuni, dzsminuni,                         &
                              l_dzsmaxuni, dzsmaxuni,                         &
                              dzsmax_bottom, dzsmin,                          &
                              nlayer_surf_sed,                                &
                              k1HW97, k2HW97,                                 &
                              fusion_para_activlayer

    namelist /namsedim_erosion/ activlayer, frmudcr2, coef_frmudcr1,          &
                                x1toce_mud, x2toce_mud,                       &
                                E0_sand_option, E0_sand_para, n_eros_sand,    &
                                E0_mud, n_eros_mud,                           &
                                ero_option, xexp_ero,                         &
                                E0_sand_Cst,                    &
                                tau_cri_option,                               &
                                tau_cri_mud_option_eroindep,                  &
                                l_peph_suspension, l_xexp_ero_cst,            &
                                l_eroindep_mud, l_eroindep_noncoh,            &
                                E0_mud_para_indep

    namelist /namsedim_bottomstress/ l_z0seduni,                              &
                                     z0seduni, z0sedmud, z0sedbedrock,        &
                                     l_fricwave, fricwav,                     &
                                     l_z0hydro_coupl_init,                    & 
                                     l_z0hydro_coupl,                         &
                                     coef_z0_coupl,                           &
                                     z0_hydro_mud,                            &
                                     z0_hydro_bed

    namelist /namsedim_deposition/ cfreshmud, csedmin, cmudcr, aref_sand,     &
                                   cvolmaxsort, cvolmaxmel, slopefac

    namelist /namsedim_lateral_erosion/ coef_erolat, coef_tauskin_lat,        &
                                        l_erolat_wet_cell, htncrit_eros 

    namelist /namsedim_consolidation/ l_consolid, xperm1, xperm2, xsigma1,    &
                                      xsigma2, csegreg, csandseg,             &
                                      dt_consolid, subdt_consol

    namelist /namsedim_diffusion/ l_diffused, choice_flxdiss_diffsed,         &
                                  xdifs1, xdifsi1,           &
                                  epdifi, fexcs, dt_diffused

    namelist /namsedim_bioturb/ l_bioturb, l_biodiffs,                        &
                                xbioturbmax_part, xbioturbk_part,             &
                                dbiotu0_part, dbiotum_part,                   &
                                xbioturbmax_diss, xbioturbk_diss,             &
                                dbiotu0_diss, dbiotum_diss,                   &
                                frmud_db_min, frmud_db_max,                   &
                                dt_bioturb, subdt_bioturb

    namelist /namsedim_morpho/ l_morphocoupl, MF,               &
                               l_MF_dhsed, l_bathy_actu,                      &
                               dt_morpho                                  

    namelist /namsedoutput/ choice_nivsed_out,                                &
                            nk_nivsed_out, ep_nivsed_out, epmax_nivsed_out,   &
                            l_outsed_flx_WS_int, l_outsed_flx_WS_all,         &
                            l_outsed_saltemp

#ifdef key_MUSTANG_V2
    namelist /namsedim_poro/ poro_option, poro_min,                           &
                             Awooster, Bwooster, Bmax_wu 

    namelist /namsedim_bedload/ l_peph_bedload, l_slope_effect_bedload,       &
                                alphabs, alphabn, hmin_bedload, l_fsusp

#ifdef key_MUSTANG_debug
    namelist /namsedim_debug/ lon_debug, lat_debug,                           &
                              i_MUSTANG_debug, j_MUSTANG_debug,               &
                              l_debug_erosion, l_debug_effdep,                &
                              date_start_debug
#endif
#endif

#ifdef key_MUSTANG_flocmod
    namelist /namflocmod/ l_ADS, l_ASH, l_COLLFRAG,                           &
                          f_dp0, f_nf, f_nb_frag, f_alpha, f_beta, f_ater,    &
                          f_ero_frac, f_ero_nbfrag, f_ero_iv, f_mneg_param,   &
                          f_collfragparam, f_dmin_frag, f_cfcst, f_fp, f_fy,  &
                          f_clim
#endif
#if !defined key_noTSdiss_insed
    namelist /namtempsed/ mu_tempsed1, mu_tempsed2, mu_tempsed3,              &
                          epsedmin_tempsed,                                   &
                          epsedmax_tempsed
#endif

CONTAINS
  
!!=============================================================================
    SUBROUTINE MUSTANG_init(ifirst, ilast, jfirst, jlast,                     &
            WATER_ELEVATION,                                                  &
#if defined MORPHODYN  
            dhsed,                                                            &
#endif
            h0fondin, z0hydro, WATER_CONCENTRATION )
    !&E------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_init  ***
    !&E
    !&E ** Purpose : initialize the sediment parameters and variables 
    !&E
    !&E ** Description : called at the beginning of the simulation
    !&E
    !&E ** Called by : mustang_init_main
    !&E
    !&E------------------------------------------------------------------------
   !! * Modules used
    USE coupler_MUSTANG,  ONLY : coupl_conv2MUSTANG
    USE sed_MUSTANG,  ONLY : sed_MUSTANG_comp_z0hydro
#ifdef key_MUSTANG_splitlayersurf
    USE sed_MUSTANG,  ONLY : sed_MUSTANG_split_surflayer
#endif
#ifdef key_MUSTANG_V2
    USE sed_MUSTANG,  ONLY : MUSTANGV2_comp_poro_mixsed
#ifdef key_MUSTANG_bedload
    USE sed_MUSTANG_HOST,  ONLY : sed_bottom_slope
#endif
#endif

    !! * Arguments
    INTEGER, INTENT(IN)                    :: ifirst, ilast, jfirst, jlast
    REAL(KIND=rsh),INTENT(IN)              :: h0fondin
    REAL(KIND=rsh),DIMENSION(ARRAY_Z0HYDRO),INTENT(INOUT)        :: z0hydro                         
    REAL(KIND=rsh),DIMENSION(ARRAY_WATER_ELEVATION),INTENT(INOUT):: WATER_ELEVATION                         
    REAL(KIND=rsh),DIMENSION(ARRAY_WATER_CONC), INTENT(IN)       :: WATER_CONCENTRATION  
#if defined MORPHODYN_MUSTANG_byHYDRO  
    REAL(KIND=rsh),DIMENSION(ARRAY_DHSED),INTENT(INOUT)          :: dhsed                       
#endif

   !! * Local declarations
    INTEGER   :: i, j, k, iv, isplit
    CHARACTER(len=lchain) :: filepc
    CHARACTER(len=lchain) :: filepc_user
    INTEGER               :: lstr
    INTEGER               :: lenstr
#ifdef key_MUSTANG_V2
    REAL(KIND=rsh)                  :: mass_tot, poro_kij, crel_mud_ini
    REAL(KIND=rsh),DIMENSION(1:nvpc):: frac_sed
    REAL(KIND=rsh),DIMENSION(1:nvp) :: mass_sed
#else
    REAL(KIND=rsh)                  :: somalp
#endif
#ifdef key_MUSTANG_flocmod
    LOGICAL :: l_0Dcase
#endif

    !! * Executable part
    h0fond = h0fondin

    ! Read default and then user namelist
    lstr = lenstr(sedname_must)
    filepc_user = sedname_must(1:lstr)
    filepc = REPFICNAMELIST//'/paraMUSTANG_default.txt'
    CALL MUSTANG_readnml(filepc) ! read default namelist
    CALL MUSTANG_readnml(filepc_user) ! read user namelist (replace default param)

    ! Allocate tables 
    CALL MUSTANG_alloc()

    ! Initialize MUSTANG parameters from namelist
    CALL MUSTANG_init_param()

    ! Floculation module
#ifdef key_MUSTANG_flocmod
    CALL flocmod_alloc(nv_mud)
#ifdef SED_TOY_FLOC_0D
    ! for 0D test case, we need to suppress all settling process
    l_0Dcase = .true.
#else
    l_0Dcase = .false.
#endif
    CALL flocmod_init(l_ADS, l_ASH, l_COLLFRAG,             &
        f_dp0, f_nf, f_nb_frag, f_alpha, f_beta, f_ater,    &
        f_ero_frac, f_ero_nbfrag, f_ero_iv, f_mneg_param,   &
        f_collfragparam, f_dmin_frag, f_cfcst, f_fp, f_fy,  &
        f_clim, diam_sed(imud1:nvpc), ros(imud1:nvpc),      &
        RHOREF, l_0Dcase, ierrorlog)
#endif

    ! Information on screen
    CALL MUSTANG_param_log()  
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Recovery of concentrations at the bottom layer
    ! which could be used to initiate concentrations in the sediment
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! here iappel = 0 
    ! only WATER_CONCENTRATION, SALINITY_MOD, TEMPERATURE_MOD are used     
    CALL coupl_conv2MUSTANG(ifirst,ilast,jfirst,jlast,0,BATHY_H0,             &
                            WATER_ELEVATION, WATER_CONCENTRATION )

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Reading new bathy issued from a previous run with morphocoupl 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !IF(l_bathy_actu) THEN
    ! **TODO** has not been implemented in CROCO yet
    ! this routine depends on hydrodynamic host model because reading save file 
    ! from a revious run
    ! CALL bathy_actu_fromfile(BATHY_H0) 
    !ENDIF

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Evaluation of slope for bedload
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#if defined key_MUSTANG_V2 && defined key_MUSTANG_bedload
#if defined MORPHODYN_MUSTANG_byHYDRO
    it_morphoYes=0
#endif
    CALL sed_bottom_slope(ifirst, ilast, jfirst, jlast, BATHY_H0)
#endif


    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Definition of initial conditions in sediment
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    IF (l_repsed) THEN
        CALL sedinit_fromfile
    ELSE
        CALL MUSTANG_sedinit(ifirst, ilast, jfirst, jlast, BATHY_H0)
    END IF
      
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Porosity estimation
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      DO j = jfirst, jlast
        DO i = ifirst, ilast
          IF (ksma(i,j) .NE. 0) THEN

#ifdef key_MUSTANG_splitlayersurf
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            !! Splitting surface layers if too thick
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            isplit = 0
            DO k = ksma(i,j), ksma(i,j)-nlayer_surf_sed+1, -1
                IF (k > ksmi(i,j)) THEN
                    IF (dzs(k,i,j) > dzsmax(i,j) + 5.0_rsh* dzsmin) isplit=1
                ENDIF
            ENDDO
            IF(isplit == 1) THEN            
                CALL sed_MUSTANG_split_surflayer(i,j,ksma(i,j))
            ENDIF
#endif

            DO k=ksmi(i,j),ksma(i,j)
#ifdef key_MUSTANG_V2
             mass_tot = 0.0_rsh
             DO iv=igrav1,imud2
               mass_sed(iv)=cv_sed(iv,k,i,j)*dzs(k,i,j) 
               mass_tot=mass_tot+mass_sed(iv)
             END DO
             DO iv=igrav1,imud2
               frac_sed(iv)=mass_sed(iv)/mass_tot
             END DO
             crel_mud_ini=(1.0_rsh-poro_mud_ini)*ros(1)
             CALL MUSTANGV2_comp_poro_mixsed(frac_sed,poro_mud_ini,  &
                                crel_mud_ini,poro_kij)
             poro(k,i,j)=poro_kij
             crel_mud(k,i,j)=crel_mud_ini
             poro_mud(k,i,j)=poro_mud_ini
             c_sedtot(k,i,j)=(1.0_rsh-poro(k,i,j))*ros(1)
             DO iv=igrav1,imud2
               cv_sed(iv,k,i,j)=c_sedtot(k,i,j)*frac_sed(iv)
             END DO
#else
             somalp=0.0_rsh
             DO iv=1,nvpc
                somalp=somalp+cv_sed(iv,k,i,j)/ros(iv)
             ENDDO
             poro(k,i,j)=1.0_rsh-somalp
#endif

            ENDDO
          ELSE
            poro(:,i,j)=0.0_rsh
#ifdef key_MUSTANG_V2
            crel_mud(:,i,j)=0.0_rsh
            poro_mud(:,i,j)=0.0_rsh
#endif 
          ENDIF
        ENDDO
      ENDDO
      
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Estimation of sediment heights and bathy bedrock for morphodynamic
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
    ! call even if no morpho, initialization of hsed in all cases, 
    ! even if initfromfile
    CALL MUSTANG_morphoinit(ifirst, ilast, jfirst, jlast, BATHY_H0, WATER_ELEVATION   &
#if defined MORPHODYN_MUSTANG_byHYDRO  
                                                ,dhsed               &
#endif
            )

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Evaluation of Z0_hydro if l_z0_coupl_init
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
    IF(l_z0hydro_coupl_init) THEN
        CALL sed_MUSTANG_comp_z0hydro(ifirst, ilast, jfirst, jlast, z0hydro)
    ENDIF

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Initialization of output tables
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
    CALL MUSTANG_init_output

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Check CPP keys and parameters for compatibility
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
    CALL MUSTANG_compatibility(ifirst, ilast, jfirst, jlast)

    END SUBROUTINE MUSTANG_init
!!===========================================================================
 
    SUBROUTINE MUSTANG_readnml(filein)
    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_readnml  ***
    !&E
    !&E ** Purpose : reads namelist file from filein
    !&E
    !&E ** Description : read namelist file paraMUSTANG
    !&E
    !&E ** Called by :  MUSTANG_init
    !&E
    !&E ** External calls : tool_datosec()
    !&E
    !&E--------------------------------------------------------------------------
    !! * Arguments
    CHARACTER(len=lchain), INTENT(IN) :: filein

    !! * Local declarations
    REAL(KIND=rlg)        :: tool_datosec

    !! * Executable part
    OPEN(unit = 50, file = filein, status = 'old', action = 'read')
 
    MPI_master_only write(*,*) '*****************************************************'
    MPI_master_only write(*,*) 'READING MUSTANG input file'
    MPI_master_only write(*,*) TRIM(filein)
    MPI_master_only write(*,*) '*****************************************************'

    READ(50, namsedim_init); rewind(50)
    READ(50, namsedim_layer); rewind(50)
    READ(50, namsedim_bottomstress); rewind(50)
    READ(50, namsedim_deposition); rewind(50)
    READ(50, namsedim_erosion); rewind(50)
#ifdef key_MUSTANG_V2
    READ(50, namsedim_poro); rewind(50)
    READ(50, namsedim_bedload); rewind(50)
#endif    
    READ(50, namsedim_lateral_erosion); rewind(50)
    READ(50, namsedim_consolidation); rewind(50)
    READ(50, namsedim_diffusion); rewind(50)
    READ(50, namsedim_bioturb); rewind(50)
    READ(50, namsedim_morpho); rewind(50)
#if !defined key_noTSdiss_insed
    READ(50, namtempsed); rewind(50)
#endif
    READ(50, namsedoutput); rewind(50)
#if defined key_MUSTANG_V2 && defined key_MUSTANG_debug
    READ(50, namsedim_debug); rewind(50)
    t_start_debug = tool_datosec(date_start_debug)
#endif
#ifdef key_MUSTANG_flocmod
    ! module FLOCULATION
    READ(50,namflocmod); rewind(50)
#endif
    CLOSE(50) 
   
    END SUBROUTINE MUSTANG_readnml
!!===========================================================================

    SUBROUTINE MUSTANG_param_log()
    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_param_log  ***
    !&E
    !&E ** Purpose : write parameters in log file
    !&E
    !&E ** Description : namelists writing for information
    !&E
    !&E ** Called by :  MUSTANG_init
    !&E
    !&E--------------------------------------------------------------------------
    !! * Local declarations
    INTEGER   :: iv

    MPI_master_only WRITE(iscreenlog, namsedim_init)
    MPI_master_only WRITE(iscreenlog, namsedim_layer)
    MPI_master_only WRITE(iscreenlog, namsedim_bottomstress)
    MPI_master_only WRITE(iscreenlog, namsedim_deposition)
#ifdef key_MUSTANG_V2
    MPI_master_only WRITE(iscreenlog, namsedim_poro)
    MPI_master_only WRITE(iscreenlog, namsedim_bedload)
#else
    MPI_master_only WRITE(iscreenlog, namsedim_erosion)
#endif
    MPI_master_only WRITE(iscreenlog, namsedim_lateral_erosion)
    MPI_master_only WRITE(iscreenlog, namsedim_consolidation)
    MPI_master_only write(iscreenlog, *) 'MUSTANG_param_log dt_consolid=', dt_consolid
    MPI_master_only WRITE(iscreenlog, namsedim_diffusion)
    MPI_master_only WRITE(iscreenlog, namsedim_bioturb)
    MPI_master_only WRITE(iscreenlog, namsedim_morpho)
#if defined key_MUSTANG_V2 && defined key_MUSTANG_debug
    MPI_master_only WRITE(iscreenlog, namsedim_debug)
#endif
#if !defined key_noTSdiss_insed
    MPI_master_only WRITE(iscreenlog, namtempsed)
#endif

#ifdef key_MUSTANG_flocmod
    !! module floculation
    MPI_master_only WRITE(iscreenlog, *) ' '
    MPI_master_only WRITE(iscreenlog, *) '    FLOCMOD'
    MPI_master_only WRITE(iscreenlog, *) '***********************'
    MPI_master_only WRITE(iscreenlog, *) 'class  diameter  volume  density  mass Ws'
    DO iv = 1, nv_mud
        MPI_master_only WRITE(iscreenlog, *) iv, f_diam(iv), f_vol(iv), f_rho(iv), f_mass(iv), f_ws(iv)
    ENDDO
    MPI_master_only WRITE(iscreenlog, *) ' '
    MPI_master_only WRITE(iscreenlog, *) ' *** PARAMETERS ***'
    MPI_master_only WRITE(iscreenlog, *) 'Primary particle size (f_dp0)                                : ', f_dp0
    MPI_master_only WRITE(iscreenlog, *) 'Fractal dimension (f_nf)                                     : ', f_nf
    MPI_master_only WRITE(iscreenlog, *) 'Flocculation efficiency (f_alpha)                            : ', f_alpha
    MPI_master_only WRITE(iscreenlog, *) 'Floc break up parameter (f_beta)                             : ', f_beta
    MPI_master_only WRITE(iscreenlog, *) 'Nb of fragments (f_nb_frag)                                  : ', f_nb_frag
    MPI_master_only WRITE(iscreenlog, *) 'Ternary fragmentation (f_ater)                               : ', f_ater
    MPI_master_only WRITE(iscreenlog, *) 'Floc erosion (% of mass) (f_ero_frac)                        : ', f_ero_frac
    MPI_master_only WRITE(iscreenlog, *) 'Nb of fragments by erosion (f_ero_nbfrag)                    : ', f_ero_nbfrag
    MPI_master_only WRITE(iscreenlog, *) 'fragment class (f_ero_iv)                                    : ', f_ero_iv
    MPI_master_only WRITE(iscreenlog, *) 'negative mass tolerated before redistribution (f_mneg_param) : ', f_mneg_param
    MPI_master_only WRITE(iscreenlog, *) 'Boolean for differential settling aggregation (L_ADS)        : ', l_ADS
    MPI_master_only WRITE(iscreenlog, *) 'Boolean for shear aggregation (L_ASH)                        : ', l_ASH
    MPI_master_only WRITE(iscreenlog, *) 'Boolean for collision fragmenation (L_COLLFRAG)              : ', l_COLLFRAG
    MPI_master_only WRITE(iscreenlog, *) 'Collision fragmentation parameter (f_collfragparam)          : ', f_collfragparam
    MPI_master_only WRITE(iscreenlog, *) 'Min concentration below which flocculation is not calculated : ', f_clim
    MPI_master_only WRITE(iscreenlog, *) ' '
    MPI_master_only WRITE(iscreenlog, *) '*** END FLOCMOD INIT *** '    
#endif
   
    END SUBROUTINE MUSTANG_param_log
!!===========================================================================
 
    SUBROUTINE MUSTANG_compatibility(ifirst, ilast, jfirst, jlast)
    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_compatibility  ***
    !&E
    !&E ** Purpose : verifies compatibility between the various parameters (sedim)
    !&E
    !&E ** Description :
    !&E
    !&E ** Called by : MUSTANG_init
    !&E
    !&E ** External calls : tool_datosec()
    !&E
    !&E--------------------------------------------------------------------------
    !! * Arguments
    INTEGER, INTENT(IN) :: ifirst, ilast, jfirst, jlast
   
    !! * Local declarations
    INTEGER   :: iv, iv2, isubs, i, j 
    REAL(rsh) :: rohomog
 
    ! test dissolved variables in sediment must be used without  key_nofluxwat_IWS
#if defined key_nofluxwat_IWS && ! defined key_noTSdiss_insed
    MPI_master_only WRITE(iwarnlog, *)
    MPI_master_only WRITE(iwarnlog, *) '****************************'
    MPI_master_only WRITE(iwarnlog, *) '        WARNING           ' 
    MPI_master_only WRITE(iwarnlog, *) '****************************'
    MPI_master_only WRITE(iwarnlog, *) ' You are using CPP key : key_nofluxwat_IWS ' 
    MPI_master_only WRITE(iwarnlog, *) ' So you are not taking into account water fluxes threw water/sediment interface ' 
    MPI_master_only WRITE(iwarnlog, *) ' which are generated by erosion/ settling/ consolidation..'
    MPI_master_only WRITE(iwarnlog, *) ' AND you are NOT using CPP key : key_noTSdiss_insed ' 
    MPI_master_only WRITE(iwarnlog, *) ' so you want to simulate dissolved variables in sediment'
    MPI_master_only WRITE(iwarnlog, *) ' You will have probably conservativity problems '
    MPI_master_only WRITE(iwarnlog, *) ' If you simulate dissolved variables in water AND in sediment, don t use key_nofluxwat_IWS '
#endif

#ifdef key_MUSTANG_flocmod
    IF (.not.l_ADS .and. .not.l_ASH) THEN
        MPI_master_only write(ierrorlog, *) 'CAUTION : incompatible flocculation kernel options : '
        MPI_master_only write(ierrorlog, *) '*****************************************************'
        MPI_master_only write(ierrorlog, *) 'l_ADS=', l_ADS
        MPI_master_only write(ierrorlog, *) 'l_ASH=', l_ASH
        MPI_master_only write(ierrorlog, *) 'simulation stopped'
        STOP
    ENDIF
#endif

    ! test compatibility tocd /consolidation
    IF (l_consolid) THEN
       MPI_master_only WRITE(iwarnlog, *)
       MPI_master_only WRITE(iwarnlog, *) '****************************'
       MPI_master_only WRITE(iwarnlog, *) ' You are taking into account CONSOLIDATION PROCESS'
       MPI_master_only WRITE(iwarnlog, *) ' and fresh deposit concentration (cfreshmud) = ', cfreshmud
       MPI_master_only WRITE(iwarnlog, *) ' and critical stress of deposition (tocd) of each sand and mud variable are equal to :'
       DO isubs = imud1, imud2
         MPI_master_only WRITE(iwarnlog, *) ' variable', NAME_SUBS,'  tocd=', tocd(isubs)
       ENDDO
       MPI_master_only WRITE(iwarnlog, *) ' If cfreshmud is big (> 100), you have to choose a tocd between 1 and 10'
       MPI_master_only WRITE(iwarnlog, *) ' otherwise you have to choose a tocd > 10 or 20'
       
    ELSE
       DO isubs = isand2+1, nvp
         IF (tocd(isubs) > 1.5 ) THEN
           MPI_master_only WRITE(iwarnlog, *) '****************************'
           MPI_master_only WRITE(iwarnlog, *) '  WARNING    '
           MPI_master_only WRITE(iwarnlog, *) ' You are not taking into account CONSOLIDATION PROCESS'
           MPI_master_only WRITE(iwarnlog, *) ' fresh deposit concentration (cfreshmud) = ', cfreshmud
           MPI_master_only WRITE(iwarnlog, *) ' and critical stress of deposition (tocd) of isusb=', isubs, &
                                              ' exceed 1.5 Pa :', tocd(isubs)
           MPI_master_only WRITE(iwarnlog, *) ' it should be smaller (0.5 to 1 ?) '
         ENDIF
       ENDDO
    ENDIF
    
    ! test ros homogen for sand in all cases
    IF (l_consolid) THEN
       rohomog = ros(irk_fil(1))
       DO iv = 2, nvpc
        IF(ros(irk_fil(iv)) .NE. rohomog) THEN
            MPI_master_only WRITE(ierrorlog, *)
            MPI_master_only WRITE(ierrorlog, *) '****************************'
            MPI_master_only WRITE(ierrorlog, *) ' You are taking into account CONSOLIDATION PROCESS'
            MPI_master_only WRITE(ierrorlog, *) ' and you must have same density for all constitutive particles'
            MPI_master_only WRITE(ierrorlog, *) ' It is not the case : (see in substance file) '
            MPI_master_only WRITE(ierrorlog, *) ' ros(1) = ', ros(irk_fil(1))
            MPI_master_only WRITE(ierrorlog, *) ' ros(', irk_fil(iv), ') = ', ros(irk_fil(iv))
            STOP
         ENDIF
       ENDDO
     ELSE
       IF(isand2 > isand1) THEN
        rohomog = ros(irk_fil(isand1))
        DO iv = isand1+1, isand2
         IF(ros(irk_fil(iv)) .NE. rohomog) THEN
            MPI_master_only WRITE(ierrorlog, *)
            MPI_master_only WRITE(ierrorlog, *) '****************************'
            MPI_master_only WRITE(ierrorlog, *) ' you must have same density for all sand particles'
            MPI_master_only WRITE(ierrorlog, *) ' It is not the case : (see in substance file) '
            MPI_master_only WRITE(ierrorlog, *) ' ros(1) = ', ros(irk_fil(isand1))
            MPI_master_only WRITE(ierrorlog, *) ' ros(', irk_fil(iv), ') = ', ros(irk_fil(iv))
            STOP
          ENDIF
          DO iv2 = isand1, iv
           IF(diam_sed(irk_fil(iv)) > diam_sed(irk_fil(iv2))) THEN
             MPI_master_only WRITE(ierrorlog, *)
             MPI_master_only WRITE(ierrorlog, *) '****************************'
             MPI_master_only WRITE(ierrorlog, *) ' The "SAND" variables should be stored in substance file '
             MPI_master_only WRITE(ierrorlog, *) ' in order of decreasing diameters (the coarsest sand to the finest sand)'
             MPI_master_only WRITE(ierrorlog, *) ' It is not the case : (see in substance file) '
             STOP
           ENDIF
          ENDDO
        ENDDO
       ENDIF
     ENDIF

! test compatibility keys sedimento+biolo
#if defined key_biolo && defined key_benthos && ! defined key_Pconstitonly_insed
    MPI_master_only WRITE(ierrorlog, *)
    MPI_master_only WRITE(ierrorlog, *) '****************************'
    MPI_master_only WRITE(ierrorlog, *) ' You are using key_MUSTANG + key_biolo + key_benthos'
    MPI_master_only WRITE(ierrorlog, *) ' It is not compatible'
    MPI_master_only WRITE(ierrorlog, *) ' You must add the key key_Pconstitonly_insed and biological variables '
    MPI_master_only WRITE(ierrorlog, *) ' will not be included in sedim module but in benthos module -in biolo)'
    MPI_master_only WRITE(ierrorlog, *)
    MPI_master_only WRITE(ierrorlog, *) ' OR you want simulate biological variables in sediment module'
    MPI_master_only WRITE(ierrorlog, *) ' and you have to use : key_MUSTANG + key_biolo without key_benthos'
    STOP
#endif

! test compatibility key_BLOOM_insed
#if defined key_BLOOM_insed && ! defined key_oxygen
    MPI_master_only WRITE(ierrorlog, *)
    MPI_master_only WRITE(ierrorlog, *) '****************************'
    MPI_master_only WRITE(ierrorlog, *) ' It is not possible to use key_BLOOM_insed  without the key_oxygen'
    MPI_master_only WRITE(ierrorlog, *) ' WARNING : use CPP key : key_oxygen in Makefile'
    STOP
#endif

   ! test compatibility bedload with MUSTANG Version 2 
#if defined key_MUSTANG_bedload && ! defined key_MUSTANG_V2
    MPI_master_only WRITE(ierrorlog, *)
    MPI_master_only WRITE(ierrorlog, *) '****************************'
    MPI_master_only WRITE(ierrorlog, *) ' You are using key_MUSTANG_bedload but '
    MPI_master_only WRITE(ierrorlog, *) ' You must add the key key_MUSTANG_V2 '
    MPI_master_only WRITE(ierrorlog, *) ' in order to use the new version V2 of MUSTANG which include bedload behavior'
    MPI_master_only WRITE(ierrorlog, *)
    STOP
#endif

#ifdef key_MUSTANG_V2
! test version V2 consolidation/poro_option
    IF (l_consolid .AND. poro_option .NE. 2) THEN
        MPI_master_only WRITE(iwarnlog, *) '****************************'
        MPI_master_only WRITE(iwarnlog, *) ' WARNING    '
        MPI_master_only WRITE(iwarnlog, *) ' You are  taking into account CONSOLIDATION PROCESS'
        MPI_master_only WRITE(iwarnlog, *) ' and have choose poro_optio n= 1'
        MPI_master_only WRITE(iwarnlog, *) ' poro_option must be = 2. it has been changed automatically'
        MPI_master_only WRITE(iwarnlog, *)
        poro_option = 2
    ENDIF
#endif

#ifdef key_MUSTANG_debug
#if defined SPHERICAL
    ! initialization of position of point where we want informations for debug  
    IF(lat_debug .NE. 0._rlg .OR. lon_debug .NE. 0._rlg) THEN
        ! locating i,j from latitude, longitude 
        i_MUSTANG_debug = 0
        j_MUSTANG_debug = 0            
        DO j = jfirst, jlast
            DO i = ifirst, ilast
                IF(LATITUDE(i,j) == lat_debug .AND. LONGITUDE(i,j) == lon_debug) THEN
                    i_MUSTANG_debug = i
                    j_MUSTANG_debug = j             
                ENDIF
            ENDDO
        ENDDO
        IF(i_MUSTANG_debug == 0) THEN
            MPI_master_only write(ierrorlog, *) 'CAUTION : the point for debug is not well known  '
            MPI_master_only write(ierrorlog, *) '*****************************************************'
            MPI_master_only write(ierrorlog, *) 'lat_debug and lon_debug do not correspond to one point in the grid'
            MPI_master_only write(ierrorlog, *) 'choose i_MUSTANG_debug et j_MUSTANG_debug directly in paraMUSTANGV2.txt'
            MPI_master_only write(ierrorlog, *) 'or change lat_debug and lon_debug'
            MPI_master_only write(ierrorlog, *) 'simulation stopped'
            stop    
        ENDIF
    ENDIF
#endif
#endif
   
    END SUBROUTINE MUSTANG_compatibility
!!===========================================================================
 
    SUBROUTINE MUSTANG_init_param()
    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_init_param  ***
    !&E
    !&E ** Purpose : initialize MUSTANG parameters  
    !&E
    !&E ** Description : 
    !&E    - compute for each sand and gravel (for 1 to isand2)  
    !&E       * diamstar
    !&E       * ws_sand
    !&E       * tetacri0
    !&E       * stresscri0
    !&E       * psi_sed (in V2)
    !&E    - compute for each sand (for isand1 to isand2) 
    !&E       * E0_sand (in V2)
    !&E    - compute for each gravel, dans and mud (for 1 to imud2) 
    !&E       * rosmrowsros
    !&E    - ros_sand_homogen = ros(isand1) if there is sand (isand1 > 0)
    !&E
    !&E ** Note : GRAVITY and RHOREF must be known 
    !&E
    !&E ** Called by :  MUSTANG_init
    !&E
    !&E--------------------------------------------------------------------------

    !! * Local declarations
    INTEGER        :: iv
    REAL(KIND=rlg) :: dtsedc, dtsedd, dtsedb  ! to compute consolidation subdt_consol
    REAL(KIND=rlg) :: tool_datosec
    REAL(KIND=rsh),PARAMETER :: shield_cri_wu = 0.03_rsh          

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

    IF (csed_mud_ini == 0.0_rsh ) csed_mud_ini = cfreshmud
    tstart_dyninsed = tool_datosec(date_start_dyninsed)
    t_dyninsed = tstart_dyninsed
    tstart_morpho = tool_datosec(date_start_morpho)
    dtsedc = 365._rlg * 86400._rlg
    MPI_master_only write(*,*) 'init dtsedc = ', dtsedc
    dtsedd = dtsedc
    dtsedb = dtsedc
    IF (l_consolid) THEN
        dtsedc = dt_consolid
        MPI_master_only write(*,*) 'init dt_consolid = ',dt_consolid
    ELSE
        dt_consolid = dtsedc
        MPI_master_only write(*,*) 'init dt_consolid = ',dt_consolid
    ENDIF
    IF (l_diffused) THEN
        dtsedd = dt_diffused
    ELSE
        dt_diffused = dtsedd
    ENDIF
    IF (l_bioturb .OR. l_biodiffs) THEN
        dtsedb = dt_bioturb
    ELSE
        dt_bioturb = dtsedb
    ENDIF
    dt_dyninsed = MIN(dtsedc,dtsedd,dtsedb)
    dtsedc = 365._rlg*86400._rlg
    dtsedb = dtsedc
    IF (l_consolid) THEN
        dtsedc = subdt_consol
    ELSE
        subdt_consol = dtsedc
    ENDIF
    IF (l_bioturb) THEN
        dtsedb = subdt_bioturb
    ELSE
        subdt_bioturb = dtsedb
    ENDIF    
    subdt_consol = MIN(dtsedb, dtsedc, dt_dyninsed)
#ifdef key_MUSTANG_V2    
    IF (l_dzsminuni) THEN 
        coeff_dzsmin = 0._rsh        
    ELSE           
        coeff_dzsmin = 1._rsh  
    ENDIF   
#endif  
    IF (l_morphocoupl) THEN
    t_morpho = tstart_morpho
        IF (l_MF_dhsed) THEN
            MF_dhsed = MF
            MF = 1._rsh
        ELSE
            MF_dhsed = 0._rsh  ! unused
            ! MF used, read from namelist
        ENDIF
    ELSE
        MF = 1.0_rsh
    ENDIF
    cexcs = 1.0_rsh - fexcs
    fws2 = fricwav * 0.5_rsh
    l_dyn_insed = .FALSE.
    IF (l_consolid .OR. l_bioturb .OR. l_diffused .OR. l_biodiffs) l_dyn_insed = .TRUE.



    ! conversion of thickness reading in mm to m for computation
    epmax_nivsed_out = epmax_nivsed_out / 1000.0_rsh
    ep_nivsed_out(:) = ep_nivsed_out(:) / 1000.0_rsh
    IF (choice_nivsed_out == 1) nk_nivsed_out = ksdmax
                               
    MPI_master_only write(*,*) 'MUSTANG_init_param dt_consolid=', dt_consolid

    DO iv = 1, isand2

      ! here, for simplification, water density = RHOREF 
      ! here, for simplification, viscosity = 1.e-6 constant
      ! 10000 = (1/viscosite**2)**(1/3)
      rosmrowsros(iv) = (ros(iv) - RHOREF) / ros(iv)
      diamstar(iv) = diam_sed(iv) * 10000.0_rsh * (GRAVITY * (ros(iv) / RHOREF - 1.0_rsh))**0.33_rsh
      ! according to Soulsby, 1997, and if viscosity = 10-6 m/s :
      ws_sand(iv)=.000001_rsh*((107.33_rsh+1.049_rsh*diamstar(iv)**3)**0.5_rsh-10.36_rsh)/diam_sed(iv)


       ! Critical shear stress in erosion law (N/m2)
      IF (tau_cri_option == 0) THEN
        tetacri0(iv)=0.3_rsh/(1.0_rsh+1.2_rsh*diamstar(iv))+0.055_rsh*(1.0_rsh-EXP(-0.02_rsh*diamstar(iv)))
        stresscri0(iv)=tetacri0(iv)*GRAVITY*(ros(iv)-RHOREF)*diam_sed(iv)
      ELSE IF (tau_cri_option == 1) THEN ! Wu and Lin (2014)
        stresscri0(iv)=GRAVITY*(ros(iv)-RHOREF)*diam_sed(iv)*shield_cri_wu
      END IF

#ifdef key_MUSTANG_V2
      psi_sed(iv)=LOG(diam_sed(iv)*1000.0_rsh)/LOG(2.0_rsh)
      IF (iv .ge. isand1) THEN
            E0_sand(iv) = MUSTANG_E0sand(diam_sed(iv), stresscri0(iv), ros(iv), ws_sand(iv))
      ENDIF
#endif
    ENDDO ! iv = 1, isand2

    DO iv=imud1,imud2
      rosmrowsros(iv)=(ros(iv)-RHOREF)/ros(iv)
    ENDDO

    ! homogeneous sand density
    IF (isand1 > 0) THEN
       iv=isand1
       ros_sand_homogen=ros(iv)
    ENDIF

    END SUBROUTINE MUSTANG_init_param
!!===========================================================================

    SUBROUTINE MUSTANG_init_hsed(cv_sedini)
    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_init_hsed  ***
    !&E
    !&E ** Purpose : initialize sediment thickness, in order to be coherent with sediment parameters
    !&E     function of cseduni, cvolmax, csed_ini of each sediment,
    !&E      valid only if cv_sed  and dzs not variable
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by :  MUSTANG_init
    !&E
    !&E--------------------------------------------------------------------------
    !! * Arguments
    REAL(KIND=rsh), INTENT(IN),DIMENSION(nvpc) :: cv_sedini

    !! * Local declarations
    INTEGER        :: iv                
    REAL(KIND=rsh) :: dzsgrv,dzssan,dzsmud,voldepgrv,voldepsan,masdepmud,     &
                        dvolsan,dmasmud,dmasmudgrav,dmasmudsand

    !! * Executable part

    ! **TODO** translate all comments in english
      !!! Volume a deposer
    voldepgrv=0.0_rsh
    voldepsan=0.0_rsh        
    DO iv=igrav1,igrav2
        voldepgrv=voldepgrv+cv_sedini(iv)/ros(iv)*hseduni  ! kg/m3 * m /(kg/m3)  (vol_g)
    ENDDO
    DO iv=isand1,isand2
        voldepsan=voldepsan+cv_sedini(iv)/ros(iv)*hseduni  ! kg/m3 * m /(kg/m3)  (vol_s)
    ENDDO

    !!! Masse a deposer
    masdepmud=0.0_rsh
    DO iv=imud1,imud2
        masdepmud=masdepmud+cv_sedini(iv)*hseduni  ! kg/m2  (m_v)
    ENDDO

    !!! nouvelle Couche de gravier
    IF(voldepgrv.GE.0.0_rsh)THEN
        dzsgrv=voldepgrv/cvolmaxsort
    ENDIF

    !!! Couche de sable
    IF(voldepsan.GE.0.0_rsh)THEN
        dvolsan=MIN(voldepsan,dzsgrv*(cvolmaxmel-cvolmaxsort)) ! volume de sable qui s insere dans le gravier (cvolmaxmel-cvolmaxsort)
        dzssan=MAX(0.,(voldepsan-dvolsan)/cvolmaxsort)
    ENDIF
    
    !!! Couche de vase
    IF(masdepmud.GE.0.0_rsh) THEN
        IF(dzsgrv >0.0_rsh .and. dzssan > 0.0_rsh) THEN
            ! melange sable+gravier  
            dmasmudgrav=dzsgrv*csed_mud_ini*(1.0_rsh-cvolmaxmel)  ! masse de vase qui s insere dans le gravier (1-cvolmaxmel)
            IF (isand2>isand1 ) THEN
                dmasmudsand=dzssan*csed_mud_ini*(1.0_rsh-cvolmaxmel) ! masse de vase qui s insere dans le sable (1-cvolmaxsort)
            ELSE
                dmasmudsand=dzssan*csed_mud_ini*(1.0_rsh-cvolmaxsort) ! masse de vase qui s insere dans le sable (1-cvolmaxsort)
            ENDIF
        ELSE IF (dzssan > 0.0_rsh .AND. isand2>isand1 ) THEN
            ! melange sables  
            dmasmudgrav=0.
            dmasmudsand=dzssan*csed_mud_ini*(1.0_rsh-cvolmaxmel) ! masse de vase qui s insere dans le sable (1-cvolmaxsort)
        ELSE IF(dzsgrv >0.0_rsh .AND. igrav2>igrav1 ) THEN
            ! melange graviers  
            dmasmudgrav=dzsgrv*csed_mud_ini*(1.0_rsh-cvolmaxmel)  ! masse de vase qui s insere dans le gravier (1-cvolmaxmel)
            dmasmudsand=0.
        ELSE
            ! un seul type de sediment 
            dmasmudgrav=dzsgrv*csed_mud_ini*(1.0_rsh-cvolmaxsort)  ! masse de vase qui s insere dans le gravier ou sable (1-cvolmaxsort)
            dmasmudsand=dzssan*csed_mud_ini*(1.0_rsh-cvolmaxsort) ! masse de vase qui s insere dans le sable (1-cvolmaxsort)
        ENDIF
        ! nouvelle couche de vase
        dzsmud=(masdepmud-dmasmudgrav-dmasmudsand)/csed_mud_ini
        IF(dzsmud <0.0_rsh) THEN
            dzsmud=0.0_rsh
            IF (masdepmud-dmasmudgrav < 0.0_rsh) THEN
                dmasmudgrav = masdepmud
                dmasmudsand = 0.0_rsh
            ELSE
                dmasmudsand = masdepmud-dmasmudgrav
            ENDIF
        ENDIF     
    dmasmud = dmasmudgrav +dmasmudsand
    ENDIF

    !!! Nouvelle Hauteur totale
    hsed_new = dzsgrv +dzssan +dzsmud

    END SUBROUTINE MUSTANG_init_hsed
!!===========================================================================================
    
    SUBROUTINE MUSTANG_sedinit(ifirst, ilast, jfirst, jlast, BATHY_H0)
    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_sedinit  ***
    !&E
    !&E ** Purpose : initialize the sediment geometry & content
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by :  MUSTANG_init_sediment only if not.l_repsed
    !&E
    !&E--------------------------------------------------------------------------

    !! * Arguments
    INTEGER,INTENT(IN)                                        :: ifirst, ilast, jfirst, jlast
    REAL(KIND=rsh),DIMENSION(ARRAY_BATHY_H0),INTENT(IN)       :: BATHY_H0      

    !! * Local declarations
    CHARACTER*50   :: filesand, filemud
    INTEGER        :: i, j, k, iv, ivg, imud_nobio, IERR_MPI
    REAL(KIND=rsh) :: somalp,dzsf,dzsgrv,dzssan,dzsmud,voldepgrv,voldepsan,masdepmud,h0max_def,     &
                        dvolsan,dmasmud,dmasmudgrav,dmasmudsand,cini_modif,tot_csed_init,dzsmaxmin

    !! * Executable part

    dzsmax(PROC_IN_ARRAY)=dzsmaxuni
    z0sed(PROC_IN_ARRAY)=z0seduni
    ksmi(PROC_IN_ARRAY)=ksmiuni
    ksma(PROC_IN_ARRAY)=0
    hsed(PROC_IN_ARRAY)=-valmanq
    dzs(ksdmin:ksdmax,PROC_IN_ARRAY)=-valmanq
    cv_sed(-1:nv_tot,ksdmin:ksdmax,PROC_IN_ARRAY)=-valmanq
    c_sedtot(ksdmin:ksdmax,PROC_IN_ARRAY)=-valmanq

    IF (l_unised)THEN
   
        ! uniform bed :
        WHERE (BATHY_H0(PROC_IN_ARRAY) /= -valmanq) ksmi(PROC_IN_ARRAY)=ksmiuni
        WHERE (BATHY_H0(PROC_IN_ARRAY) /= -valmanq) ksma(PROC_IN_ARRAY)=ksmauni
        DO k=ksdmin,ksdmax
            WHERE (BATHY_H0(PROC_IN_ARRAY) /= -valmanq) dzs(k,PROC_IN_ARRAY)=0.0_rsh
            WHERE (BATHY_H0(PROC_IN_ARRAY) /= -valmanq) c_sedtot(k,PROC_IN_ARRAY)=0.0_rsh
            DO iv=-1,nv_tot
                WHERE (BATHY_H0(PROC_IN_ARRAY) /= -valmanq) cv_sed(iv,k,PROC_IN_ARRAY)=0.0_rsh
            ENDDO
        ENDDO
 
 ! if biological variable (detrital organic matter) are constitutive particulaite variable
 !     readjustment of cini_sed (fraction of csedtot : somme=1)
        cini_modif=0.0_rsh
        imud_nobio=1
        DO iv=1,isand2
            cini_modif=cini_modif+cini_sed(iv)
        ENDDO
        DO iv=imud1,imud2
            IF(unit_modif_mudbio_N2dw(irk_fil(iv)) .NE. 1.0_rsh) THEN
                ! variable bio (Ndet) mis en constitutive mais exprime dans la donnee initiale en mmole/kg de sediment total
                !cini_sed(iv)=cini_sed(iv)/cseduni
                cini_modif=cini_modif+cini_sed(iv)
            ELSE 
                ! vase non detritique
                IF(cini_sed(iv) > (1.0_rsh-cini_modif)) imud_nobio=iv
                cini_modif=cini_modif+cini_sed(iv)
            ENDIF    
        ENDDO
        IF(cini_modif > 1.0_rsh) cini_sed(imud_nobio)=cini_sed(imud_nobio)-cini_modif+1.0_rsh
   
        tot_csed_init=0.0_rsh
        DO iv=1,nvpc
            ! constitutive variables : cini_sed in fraction of cseduni
            tot_csed_init=tot_csed_init+cini_sed(iv)              
        ENDDO

        IF( abs(tot_csed_init - 1.0_rsh) > epsilon_MUSTANG ) THEN
            MPI_master_only WRITE(ierrorlog,*)' '
            MPI_master_only WRITE(ierrorlog,*)' inital concentrations of constitutive particulate matter'
            MPI_master_only WRITE(ierrorlog,*)' have not been well defined'
            MPI_master_only WRITE(ierrorlog,*)' See variable.dat to give fraction of total concentration '
            MPI_master_only WRITE(ierrorlog,*)' for each constitutive variable (GRAV, SAND, MUD) (but not MUDB) '
            MPI_master_only WRITE(ierrorlog,*)' : the sum of all must be =1 and it is :'
            MPI_master_only WRITE(ierrorlog,*)  tot_csed_init
            MPI_master_only WRITE(ierrorlog,*)' difference with 1 is :'
            MPI_master_only WRITE(ierrorlog,*)  abs(tot_csed_init - 1.0_rsh) 
            STOP   
        ELSE
            DO iv=1,nvpc
                ! force fraction sum to equal 1.
                cini_sed(iv) = cini_sed(iv)/tot_csed_init            
            ENDDO    
        ENDIF
      
        ! ajustement du hsed et cv_sed
        ! ==============================
        IF(l_init_hsed .AND. hseduni > 0.0_rsh) THEN
            ALLOCATE(cv_sedini(1:nvpc))
            !!! Concentration massique
            DO iv=1,nvpc
                ! constitutive variables : cini_sed in % of cseduni
                cv_sedini(iv)=cseduni*cini_sed(iv)
            ENDDO

            CALL MUSTANG_init_hsed(cv_sedini)
            MPI_master_only WRITE(iscreenlog,*)'NEW TOTAL SEDIMENT THICKNESS AFTER ADJUSTEMENT = ',hsed_new
            cseduni=cseduni*hseduni/hsed_new         
            hseduni=hsed_new
            WHERE (BATHY_H0(PROC_IN_ARRAY) /= -valmanq) hsed(PROC_IN_ARRAY)=hseduni
            DEALLOCATE(cv_sedini)
        ELSE
            WHERE (BATHY_H0(PROC_IN_ARRAY) /= -valmanq) hsed(PROC_IN_ARRAY)=hseduni
        ENDIF

        DO j=jfirst,jlast
        DO i=ifirst,ilast
            IF(BATHY_H0(i,j) /= -valmanq) THEN
                DO k=ksmi(i,j),ksma(i,j)
                    c_sedtot(k,i,j)=0.0_rsh
                    DO iv=1,nvpc
                        ! constitutive variables : cini_sed in % of cseduni
                        cv_sed(iv,k,i,j)=cseduni*cini_sed(iv)        
                        c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
                    ENDDO
                    ! porosity estimation moved after dzs evaluation below
#ifdef key_Pconstitonly_insed
                    cv_sed(nvpc+1:nv_state,k,i,j)=0.0_rsh 
#else
                    DO iv=nvpc+1,nvp
                        ! non constitutive variables : cini_sed in conc/kg of sediment
                        IF (irkm_var_assoc(iv) >0) THEN
                            cv_sed(iv,k,i,j)=cini_sed(iv)*cv_sed(irkm_var_assoc(iv),k,i,j)
                        ELSE
                            cv_sed(iv,k,i,j)=cini_sed(iv)*c_sedtot(k,i,j)
                        END IF
                    ENDDO
                    IF(l_initsed_vardiss) THEN
                    ! initalisation avec cv_wat (CW) du fond
                        DO iv=nvp+1,nv_adv
                            ! dissolved variables : cini_sed in conc/m3 of pore water
                            cv_sed(iv,k,i,j)=cw_bottom_MUSTANG(iv,i,j)
                        ENDDO
                    ELSE
                        ! initalisation avec cini (variable.dat)
                        DO iv=nvp+1,nv_adv
                            ! dissolved variables : cini_sed in conc/m3 of pore water
                            cv_sed(iv,k,i,j)=cini_sed(iv)
                        ENDDO
                    ENDIF
                    DO iv=nv_adv+1,nv_state
                        ! + fixed variables : cini_sed in conc/m3 of sediment
                        cv_sed(iv,k,i,j)=cini_sed(iv)
                    ENDDO
#endif
                    IF(l_initsed_vardiss) THEN
                    ! initalisation avec valeurs du fond
                        cv_sed(-1,k,i,j)=temp_bottom_MUSTANG(i,j)
                        cv_sed( 0,k,i,j)=sal_bottom_MUSTANG(i,j)
                    ELSE
                        cv_sed(-1,k,i,j)=tini_sed
                        cv_sed( 0,k,i,j)=sini_sed
                    ENDIF
                ENDDO
            ENDIF
        ENDDO
        ENDDO

    ELSE

        ! non uniform bed
        MPI_master_only write(ierrorlog,*) 'NON UNIFORM BED COVERAGE'
        MPI_master_only write(ierrorlog,*) 'you have to read a netcdf file describing for each grid cell (file_inised)'
        MPI_master_only write(ierrorlog,*) 'the values for hsed(i,j),z0sed(i,j),ksmi(i,j),ksma(i,j)'
        MPI_master_only write(ierrorlog,*) 'c_sedtot(i,j),cini_sed(iv,iv=1:nv_state)'
        STOP
      
    ENDIF
!
!   END NON UNIFORM BED COVERAGE L_UNISED
!   -------------------------------------------------
!
!   regular vertical mesh
!   ---------------------

    DO j=jfirst,jlast
        DO i=ifirst,ilast
            IF(ksma(i,j).NE.0)THEN
                dzsf=hsed(i,j)/float(ksma(i,j)-ksmi(i,j)+1)
                DO k=ksmi(i,j),ksma(i,j)
                    dzs(k,i,j)=dzsf
                ENDDO
            ELSE
                dzs(:,i,j)=0.0_rsh
            ENDIF
        ENDDO
    ENDDO

! DZSmax non uniform
    IF(.NOT. l_dzsmaxuni) THEN
        IF(dzsmaxuni == 0.0_rsh) THEN
            MPI_master_only write(ierrorlog,*)' ERROR in paraMUSTANG : dzsmaxuni =0. and l_repsed=.FALSE.'
            MPI_master_only write(ierrorlog,*)' if dzsmaxuni =0. dzsmax must be read in an initial file (l_repsed=T)'
            MPI_master_only write(ierrorlog,*)' if l_repsed=.FALSE. dzsmax must be evaluated from dzsmaxuni not null'
            STOP
        ENDIF

        ! definition of dzsmax
        ! estimated as a function of depth and dzsmaxuni considered as the max value of dzsmax
        ! here minimum value of dzsmax in depth = dzsmaxuni / 100
        ! max depth giving the lowest value of dzsma to inform, setting by default = 50m
       dzsmaxmin=dzsmaxuni/100._rsh
       h0max_def=50.0_rsh
       MPI_master_only write(iscreenlog,*)
       MPI_master_only write(iscreenlog,*)'DZSMAX NON UNIFORM'
       MPI_master_only write(iscreenlog,*)'by default linearly computed in MUSTANG_sedinit (sedim.F90)'
       MPI_master_only write(iscreenlog,*)'from dzsmaxuni (',dzsmaxuni,'m) to dzsmaxuni/100 (',dzsmaxmin,'m) depending on water depth'
       MPI_master_only write(iscreenlog,*)'dzsmax minimum for depth > ',h0max_def
       MPI_master_only write(iscreenlog,*)'PROGRAM your own initialization of dzsmax (or READ a file) in MUSTANG_sedinit'
       dzsmax(PROC_IN_ARRAY)=(h0max_def-BATHY_H0(PROC_IN_ARRAY))/h0max_def*dzsmaxuni+dzsmaxmin
       WHERE (dzsmax(PROC_IN_ARRAY) > dzsmaxuni ) dzsmax(PROC_IN_ARRAY)=dzsmaxuni
       WHERE (dzsmax(PROC_IN_ARRAY) < dzsmaxmin ) dzsmax(PROC_IN_ARRAY)=dzsmaxmin

       ! or to be read in a data file (to be programmed) 
        
    ENDIF 

    END SUBROUTINE MUSTANG_sedinit
  
!!===========================================================================

    SUBROUTINE MUSTANG_init_output
    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_init_output  ***
    !&E
    !&E ** Purpose : define output arrays in sediment
    !&E
    !&E ** Description : 
    !&E
    !&E ** Called by :  MUSTANG_init
    !&E 
    !&E--------------------------------------------------------------------------
    !! * Modules used
#if defined key_BLOOM_insed
    USE bioloinit,  ONLY : ndiag_tot, ndiag_3d_sed, ndiag_2d_sed, ndiag_1d, ndiag_2d
#endif

    !! * Local declarations
    INTEGER        :: k, nk_nivsed_outlu, nv_out
    REAL(KIND=rsh) :: dzs_estim
                        
    !!--------------------------------------------------------------------------
    !! * Executable part

    MPI_master_only WRITE(iscreenlog, *)
    MPI_master_only WRITE(iscreenlog, *) '***************************************************************'
    MPI_master_only WRITE(iscreenlog, *) '*****************     SEDIMENT OUTPUT  ************************'
    MPI_master_only WRITE(iscreenlog, *) '***************************************************************'
    MPI_master_only WRITE(iscreenlog, *)
    nk_nivsed_outlu = nk_nivsed_out
    IF(choice_nivsed_out == 1 ) THEN
        ALLOCATE(nivsed_out(0:ksdmax+1))
        nivsed_out(0)=0
        DO k=ksdmin,ksdmax
            nivsed_out(k)=ksdmax+1-k
        END DO
        nivsed_out(ksdmax+1)=0
        MPI_master_only WRITE(iscreenlog,*)'results in sediment will be save on all the sediment layers (ksdmax)'
    ELSE IF(choice_nivsed_out == 2 ) THEN
        IF(nk_nivsed_outlu > ksdmax) THEN
            nk_nivsed_out=ksdmax
        ENDIF
        ALLOCATE(nivsed_out(0:nk_nivsed_out+1))
        nivsed_out(0)=0
        DO k=1,nk_nivsed_out
            nivsed_out(k)=nk_nivsed_out+1-k
        END DO
        nivsed_out(nk_nivsed_out+1)=0       
        MPI_master_only WRITE(iscreenlog,*)'results in sediment will be save only for the first ', nk_nivsed_out, &
                                ' layers from sediment surface, one average for the last'
    ELSE IF(choice_nivsed_out == 3 ) THEN
        dzs_estim=MIN(dzsmaxuni,hseduni/ksmauni)
        nk_nivsed_out = MIN(ksdmax, INT(epmax_nivsed_out / dzs_estim) +3) 
        ALLOCATE(nivsed_out(0:nk_nivsed_out+1))
        nivsed_out(0)=0
        DO k=1,nk_nivsed_out
            nivsed_out(k)=nk_nivsed_out+1-k
        END DO
        nivsed_out(nk_nivsed_out+1)=0

        MPI_master_only WRITE(iscreenlog,*)'results in sediment will be save from sediment surface till ',  &
                            epmax_nivsed_out*1000,'mmeters (integration below)'
        MPI_master_only WRITE(iscreenlog,*)'the last layer saved in the sediment output file will be that  ',  &
                            'the bottom of which exceeds the desired maximum thickness '          
!           MPI_master_only WRITE(iscreenlog,*)'number of sediment layers in output file :',nk_nivsed_out,'INT(',epmax_nivsed_out,'/',dzs_estim,'+3)'
        IF (.NOT. l_dzsmaxuni) THEN
            MPI_master_only WRITE(iscreenlog,*)
!             MPI_master_only WRITE(iscreenlog,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
!             MPI_master_only WRITE(iscreenlog,*)'WARNING : the maximum number of saved sed. layers has been evaluated from dzsmaxuni=',dzsmaxuni
            MPI_master_only WRITE(iscreenlog,*)'            if dzsmax is not uniform, it could be smaller than dzsmaxuni and  &
                        then there may be points where '
            MPI_master_only WRITE(iscreenlog,*)'             this number is not sufficient to describe the maximum thickness'
            MPI_master_only WRITE(iscreenlog,*)
            MPI_master_only WRITE(iscreenlog,*)'This option is not recommended with l_dzsmaxuni=.False.'
!             MPI_master_only WRITE(iscreenlog,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        ENDIF

    ELSE IF(choice_nivsed_out == 4 ) THEN
        IF(nk_nivsed_outlu > 5 ) THEN
            MPI_master_only WRITE(iscreenlog,*)'ERROR in namelist namsedoutput (paraMUSTANG.txt) '
            MPI_master_only WRITE(iscreenlog,*)'nk_nivsed_out must be <=5'
            MPI_master_only WRITE(iscreenlog,*)'nk_nivsed_out is automatically set to 5 '
            MPI_master_only WRITE(iscreenlog,*)'and an latter layer (6th) is added to integrate till the bottom sediment  '
            nk_nivsed_outlu = 5
        ENDIF        
        nk_nivsed_out =  nk_nivsed_outlu+1   
        ALLOCATE(ep_nivsed_outp1(nk_nivsed_out))
        ep_nivsed_outp1(1:nk_nivsed_outlu)=ep_nivsed_out(1:nk_nivsed_outlu)
        ep_nivsed_outp1(nk_nivsed_out)=10.0_rsh  

        MPI_master_only WRITE(iscreenlog,*)'results in sediment will be save on ',nk_nivsed_out-1, &
            'integrated layers whom thickness are constant and given by user - first is  sediment surface'
        MPI_master_only WRITE(iscreenlog,*)'the ',nk_nivsed_out,'the layer will be an integrated layer till the bottom' 

        ALLOCATE(nivsed_out(0:nk_nivsed_out+1))
        nivsed_out(0)=0
        DO k=1,nk_nivsed_out
            nivsed_out(k)=nk_nivsed_out+1-k
        END DO
        nivsed_out(nk_nivsed_out+1)=0

    ENDIF !choice_nivsed_out

    ALLOCATE(var3D_dzs(nk_nivsed_out,PROC_IN_ARRAY))     
    ALLOCATE(var3D_TEMP(nk_nivsed_out,PROC_IN_ARRAY))
    ALLOCATE(var3D_SAL(nk_nivsed_out,PROC_IN_ARRAY))
#ifdef key_Pconstitonly_insed
    nv_out=nvpc
#else
    nv_out=nv_adv
#endif   
    ALLOCATE(var3D_cvsed(nk_nivsed_out,PROC_IN_ARRAY,nv_out))
#ifdef key_BLOOM_insed
    ALLOCATE(var2D_diagsed(PROC_IN_ARRAY,ndiag_1d+ndiag_2d-ndiag_2d_sed+1:ndiag_1d+ndiag_2d))
    ALLOCATE(var3D_diagsed(nk_nivsed_out,PROC_IN_ARRAY,ndiag_tot-ndiag_3d_sed+1:ndiag_tot))
#endif
#if defined key_MUSTANG_specif_outputs               
    ALLOCATE(var3D_specifout(nk_nivsed_out,PROC_IN_ARRAY,nv_out3Dk_specif))
#endif

    END SUBROUTINE MUSTANG_init_output
!!===========================================================================

    SUBROUTINE MUSTANG_morphoinit(ifirst, ilast, jfirst, jlast, BATHY_H0, WATER_ELEVATION  &
#if defined MORPHODYN_MUSTANG_byHYDRO  
                  , dhsed                                                &
#endif
                  )

    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE morphoinit  ***
    !&E
    !&E ** Purpose : initialization of hsed, hsed0, hsed_previous,
    !&E                 morpho0 
    !&E
    !&E ** Description :
    !&E
    !&E ** Called by :  MUSTANG_init_sediment 
    !&E
    !&E--------------------------------------------------------------------------

    !! * Arguments 
    INTEGER, INTENT(IN)                    :: ifirst, ilast, jfirst, jlast
    REAL(KIND=rsh),DIMENSION(ARRAY_BATHY_H0),INTENT(INOUT)        :: BATHY_H0                         
    REAL(KIND=rsh),DIMENSION(ARRAY_WATER_ELEVATION),INTENT(INOUT) :: WATER_ELEVATION
#if defined MORPHODYN_MUSTANG_byHYDRO  
    REAL(KIND=rsh),DIMENSION(ARRAY_DHSED),INTENT(INOUT)           :: dhsed                       
#endif
    !! * Local declarations
    INTEGER :: i, j, k

    !! * Executable part

    IF(l_morphocoupl) THEN
        ! morpho = 1 if morphodynamic effective
        ! morpho = 0 if depth not vary
        morpho0(ARRAY_morpho)=1.0_rsh
        !  **TODO** To Program
        !   morpho0(boundaries)=0.0_rsh
    ENDIF
 
 
    !!  initialisation of hsed  inside the domain  !!!
    !!  eliminate the meshes at open boundaries
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    DO j=jfirst,jlast
    DO i=ifirst,ilast
        !  WARNING : not in first and last mesh on i and j axis (boundaries)
        hsed(i,j)=0.0_rsh
        DO k=ksmi(i,j),ksma(i,j)
            hsed(i,j)=hsed(i,j)+dzs(k,i,j)
        ENDDO
    ENDDO
    ENDDO       

    IF(l_morphocoupl)THEN
        DO j=jfirst,jlast
        DO i=ifirst,ilast
            !  WARNING: not in first and last mesh on i and j axis (boundaries)
            hsed_previous(i,j)=hsed(i,j)
        ENDDO
        ENDDO       
        
        IF(.NOT.l_repsed)THEN   
            hsed0(:,:)=0.0_rsh
            DO j=jfirst,jlast
            DO i=ifirst,ilast
                !  WARNING : not in first and last mesh on i and j axis (boundaries)
                hsed0(i,j)=hsed(i,j)
            ENDDO
            ENDDO
        ENDIF                         

    ENDIF   ! endif l_morphocoupl

    IF(l_morphocoupl)THEN    

        IF(l_repsed)THEN          

        !!  update bathy and water elevation if initalisation from file !!!
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        DO j=jfirst,jlast
        DO i=ifirst,ilast
            !ssh(i,j)=MAX(ssh(i,j),-BATHY_H0(i,j))         
            SURF_ELEVATION_ij=MAX(SURF_ELEVATION_ij,-BATHY_H0(i,j))         
        ENDDO
        ENDDO
         
        !!  echange MPI of BATHY_H0 and WATER_ELEVATION for neighboring cells
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      
        !! To Program, has not been program for CROCO **TODO** CALL sed_exchange_hxe_MARS(1,xh0=BATHY_H0,xssh=WATER_ELEVATION)

        ENDIF   ! endif l_repsed



#if defined MORPHODYN_MUSTANG_byHYDRO
        DO j=jfirst,jlast
        DO i=ifirst,ilast
            dhsed(i,j)=hsed0(i,j)-hsed(i,j)
#ifdef key_MUSTANG_debug
            IF (i==i_MUSTANG_debug .AND. j==j_MUSTANG_debug ) THEN
                MPI_master_only write(*,*) 'dhsed(',i,',',j,') initial:',dhsed(i,j),hsed0(i,j),hsed(i,j)
            ENDIF
#endif
        ENDDO
        ENDDO
#endif

    ENDIF        ! endif l_morphocoupl

    END SUBROUTINE MUSTANG_morphoinit
!!===========================================================================

    SUBROUTINE MUSTANG_alloc()
    !&E--------------------------------------------------------------------------
    !&E                 ***  ROUTINE MUSTANG_alloc  ***
    !&E
    !&E ** Purpose : allocation of arrays relative to sediment
    !&E
    !&E ** Called by : MUSTANG_init
    !&E
    !&E--------------------------------------------------------------------------
    
    !! * Local declarations
    INTEGER :: iv

    !! * Executable part
    ALLOCATE(ws_sand(nvp))        
    ALLOCATE(diamstar(nvp))
    ALLOCATE(rosmrowsros(nvp))
    ALLOCATE(stresscri0(nvp))
    ALLOCATE(tetacri0(nvp))
    ALLOCATE (typart(-1:nv_adv))
    typart(-1:0)=0.0_rsh
    typart(1:nvpc)=1.0_rsh
    typart(nvpc+1:nv_adv)=0.0_rsh
#ifdef key_MUSTANG_V2
    ALLOCATE(E0_sand(nvp))
    E0_sand(1:nvp)=0.0_rsh
#endif   
#if  ! defined key_noTSdiss_insed
    ALLOCATE(ivdiss(-1:nv_adv-nvp))
#endif
 
 
    !  allocation of spatial variables  
    !  dimensions defined dans coupler_define_MUSTANG.h
    ALLOCATE(ksmi(PROC_IN_ARRAY))
    ALLOCATE(ksma(PROC_IN_ARRAY))
    ALLOCATE(hsed(PROC_IN_ARRAY))  
    ALLOCATE(z0sed(PROC_IN_ARRAY))
    ALLOCATE(tauskin(PROC_IN_ARRAY))
    ALLOCATE(tauskin_c(PROC_IN_ARRAY))
    ALLOCATE(tauskin_w(PROC_IN_ARRAY))
    ALLOCATE(ustarbot(PROC_IN_ARRAY))
    ALLOCATE(dzsmax(PROC_IN_ARRAY))
    ALLOCATE(htot(PROC_IN_ARRAY_m2p2))
    ALLOCATE(alt_cw1(PROC_IN_ARRAY))
    ALLOCATE(epn_bottom_MUSTANG(PROC_IN_ARRAY_m1p2))  
    ALLOCATE(sal_bottom_MUSTANG(PROC_IN_ARRAY_m1p2))  
    ALLOCATE(temp_bottom_MUSTANG(PROC_IN_ARRAY_m1p2))
    ALLOCATE(cw_bottom_MUSTANG(nv_tot,PROC_IN_ARRAY_m1p2))
    ALLOCATE(ws3_bottom_MUSTANG(nvp,PROC_IN_ARRAY_m1p2))  
    ALLOCATE(roswat_bot(PROC_IN_ARRAY))  
    ksmi(PROC_IN_ARRAY) = 0
    ksma(PROC_IN_ARRAY) = 0
    hsed(PROC_IN_ARRAY) = 0.0_rsh
    z0sed(PROC_IN_ARRAY) = 0.0_rsh
    tauskin(PROC_IN_ARRAY) = 0.0_rsh
    tauskin_c(PROC_IN_ARRAY) = 0.0_rsh
    tauskin_w(PROC_IN_ARRAY) = 0.0_rsh
    ustarbot(PROC_IN_ARRAY) = 0.0_rsh
    dzsmax(PROC_IN_ARRAY) = 0.0_rsh
    htot(PROC_IN_ARRAY_m2p2) = 0.0_rsh
    alt_cw1(PROC_IN_ARRAY) = 0.0_rsh
    epn_bottom_MUSTANG(PROC_IN_ARRAY_m1p2) = 0.0_rsh
    sal_bottom_MUSTANG(PROC_IN_ARRAY_m1p2) = 0.0_rsh
    temp_bottom_MUSTANG(PROC_IN_ARRAY_m1p2) = 0.0_rsh
    cw_bottom_MUSTANG(nv_tot,PROC_IN_ARRAY_m1p2) = 0.0_rsh
    ws3_bottom_MUSTANG(nvp,PROC_IN_ARRAY_m1p2) = 0.0_rsh
    roswat_bot(PROC_IN_ARRAY) = 0.0_rsh
#if ! defined key_noTSdiss_insed
    ALLOCATE(phitemp_s(PROC_IN_ARRAY))
    phitemp_s(PROC_IN_ARRAY) = 0.0_rsh
#endif
#ifdef key_MUSTANG_V2
    ALLOCATE(sigmapsg(ksdmin:ksdmax))
    ALLOCATE(stateconsol(ksdmin:ksdmax))
    ALLOCATE(permeab(ksdmin:ksdmax))
    ALLOCATE(psi_sed(nvp))
    psi_sed(1:nvp) = 0.0_rsh
#endif
#ifdef key_sand2D
    ALLOCATE(rouse2D(nv_adv,PROC_IN_ARRAY))
    ALLOCATE(sum_tmp(nv_adv,PROC_IN_ARRAY))
    rouse2D(1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
    sum_tmp(1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
#endif
    
    ALLOCATE(cv_sed(-1:nv_tot,ksdmin:ksdmax,PROC_IN_ARRAY))
    ALLOCATE(c_sedtot(ksdmin:ksdmax,PROC_IN_ARRAY))
    ALLOCATE(poro(ksdmin:ksdmax,PROC_IN_ARRAY))
    ALLOCATE(dzs(ksdmin:ksdmax,PROC_IN_ARRAY))
    ALLOCATE(flx_s2w(-1:nv_adv,PROC_IN_ARRAY))
    ALLOCATE(flx_w2s(-1:nv_adv,PROC_IN_ARRAY))
    ALLOCATE(flx_w2s_sum(-1:nv_adv,PROC_IN_ARRAY))
    ALLOCATE(corflux(nv_adv, PROC_IN_ARRAY_m1p1))
    ALLOCATE(corfluy(nv_adv, PROC_IN_ARRAY_m1p1))
    ALLOCATE(fludif(-1:nv_adv,PROC_IN_ARRAY))
    ALLOCATE(fluconsol(-1:nv_adv,PROC_IN_ARRAY))
    ALLOCATE(fluconsol_drycell(-1:nv_adv,PROC_IN_ARRAY))
    ALLOCATE(flu_dyninsed(-1:nv_adv,PROC_IN_ARRAY))
    ALLOCATE(gradvit(NB_LAYER_WAT,PROC_IN_ARRAY))

    cv_sed(-1:nv_tot,ksdmin:ksdmax,PROC_IN_ARRAY) = 0.0_rsh
    c_sedtot(ksdmin:ksdmax,PROC_IN_ARRAY) = 0.0_rsh
    poro(ksdmin:ksdmax,PROC_IN_ARRAY) = 0.0_rsh
    dzs(ksdmin:ksdmax,PROC_IN_ARRAY) = 0.0_rsh
    flx_s2w(-1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
    flx_w2s(-1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
    flx_w2s_sum(-1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
    corflux(1:nv_adv, PROC_IN_ARRAY_m1p1) = 1.0_rsh
    corfluy(1:nv_adv, PROC_IN_ARRAY_m1p1) = 1.0_rsh
    fludif(-1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
    fluconsol(-1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
    fluconsol_drycell(-1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
    flu_dyninsed(-1:nv_adv,PROC_IN_ARRAY) = 0.0_rsh
    gradvit(NB_LAYER_WAT,PROC_IN_ARRAY) = 0.0_rsh

#ifdef key_MUSTANG_specif_outputs
    ! outputs
    ! variables 3D /k
    nv_out3Dk_specif = 1
     ! 1 : poro_save  
#if defined key_MUSTANG_add_consol_outputs && defined key_MUSTANG_V2
    nv_out3Dk_specif = nv_out3Dk_specif + 8
     ! 2 : loadograv_save
     ! 3 : permeab_save
     ! 4 : sigmapsg_save
     ! 5 : dtsdzs_save
     ! 6 : hinder_save
     ! 7 : sed_rate_save
     ! 8 : sigmadjge_save
     ! 9 : stateconsol_save
#endif
    ALLOCATE(varspecif3Dk_save(nv_out3Dk_specif,ksdmin:ksdmax,PROC_IN_ARRAY))
    varspecif3Dk_save(1:nv_out3Dk_specif,ksdmin:ksdmax,PROC_IN_ARRAY) = 0.0_rsh

    nv_out3Dnv_specif = 3
     ! 1 : toce_save
     ! 2 : flx_s2w_save
     ! 3 : flx_w2s_save
#ifdef key_MUSTANG_V2
    nv_out3Dnv_specif = nv_out3Dnv_specif + 1
     ! 4 : pephm_fcor_save  
#ifdef key_MUSTANG_bedload
    nv_out3Dnv_specif = nv_out3Dnv_specif + 4
     ! 5 : flx_bx
     ! 6 : flx_by
     ! 7 : bil_bedload
     ! 8 : fsusp
#endif
#endif
    ALLOCATE(varspecif3Dnv_save(nv_out3Dnv_specif,nvpc,PROC_IN_ARRAY))
    ALLOCATE(varspecif3Dnv_out(nv_out3Dnv_specif,nvpc,PROC_IN_ARRAY))
    varspecif3Dnv_save(1:nv_out3Dnv_specif,1:nvpc,PROC_IN_ARRAY)= 0.0_rsh
  
    nv_out2D_specif = 2
     ! 1 : frmudsup 
     ! 2 : dzs_ksmax 
#ifdef key_MUSTANG_V2
    nv_out2D_specif = nv_out2D_specif + 13
     ! 3 : dzs_aclay_comp_save
     ! 4 : dzs_aclay_kept_save
     ! 5 : tero_noncoh (cumulated time (in hours) elapsed in non cohesive regime)
     ! 6 : tero_coh (cumulated time (in hours) elapsed in cohesive regime)
     ! 7 : pct_iter_noncoh
     ! 8 : pct_iter_coh
     ! 9 : niter_ero
     ! 10: z0sed
     ! 11 : flx_s2w_noncoh
     ! 12 : flx_w2s_noncoh
     ! 13 : flx_s2w_coh
     ! 14 : flx_w2s_coh
     ! 15: z0hydro (if l_z0hydro_coupl)
#ifdef key_MUSTANG_bedload
    nv_out2D_specif = nv_out2D_specif + 3
     ! 16 : flx_bx_int
     ! 17 : flx_by_int
     ! 18 : bil_bedload_int
#endif
!   end version MUSTANG V2
#endif

    ALLOCATE(varspecif2D_save(nv_out2D_specif,PROC_IN_ARRAY))
    ALLOCATE(varspecif2D_out(nv_out2D_specif,PROC_IN_ARRAY))
    varspecif2D_save(1:nv_out2D_specif,PROC_IN_ARRAY) = 0.0_rsh
   
!   end specifs output
#endif

#ifdef key_MUSTANG_V2
    ALLOCATE(poro_mud(ksdmin:ksdmax,PROC_IN_ARRAY))
    ALLOCATE(crel_mud(ksdmin:ksdmax,PROC_IN_ARRAY))
    ALLOCATE(l_isitcohesive(PROC_IN_ARRAY))
    poro_mud(ksdmin:ksdmax,PROC_IN_ARRAY) = 0.0_rsh
    crel_mud(ksdmin:ksdmax,PROC_IN_ARRAY) = 0.0_rsh
    l_isitcohesive(PROC_IN_ARRAY) = .FALSE.

#ifdef key_MUSTANG_bedload
    ALLOCATE( flx_bx(1:nvp,PROC_IN_ARRAY_m1p1)) ! Warning /Baptiste : m1p1 sur les 2 indices au lieu d1
    ALLOCATE( flx_by(1:nvp,PROC_IN_ARRAY_m1p1)) ! Warning /Baptiste : m1p1 sur les 2 indices au lieu d1
    ALLOCATE( slope_dhdx(PROC_IN_ARRAY))
    ALLOCATE( slope_dhdy(PROC_IN_ARRAY))
    ALLOCATE( sedimask_h0plusxe(PROC_IN_ARRAY_m1p1)) ! Warning /Baptiste : m1p1 sur les 2 indices au lieu d1
    flx_bx(1:nvp,PROC_IN_ARRAY_m1p1)=0.0_rsh 
    flx_by(1:nvp,PROC_IN_ARRAY_m1p1)=0.0_rsh
    slope_dhdx(PROC_IN_ARRAY)=0.0_rsh
    slope_dhdy(PROC_IN_ARRAY)=0.0_rsh
    sedimask_h0plusxe(PROC_IN_ARRAY_m1p1)=0.0_rsh
#endif
#endif

    ALLOCATE(phieau_s2w(PROC_IN_ARRAY))
    ALLOCATE(phieau_s2w_drycell(PROC_IN_ARRAY))
    ALLOCATE(phieau_s2w_consol(PROC_IN_ARRAY))
    phieau_s2w(PROC_IN_ARRAY)=0.0_rlg
    phieau_s2w_drycell(PROC_IN_ARRAY)=0.0_rlg
    phieau_s2w_consol(PROC_IN_ARRAY)=0.0_rlg

    ALLOCATE( raphbx(PROC_IN_ARRAY_m1p1), raphby(PROC_IN_ARRAY_m1p1) )
    ALLOCATE( tauskin_x(PROC_IN_ARRAY), tauskin_y(PROC_IN_ARRAY) )
    raphbx(PROC_IN_ARRAY_m1p1)=0.0_rsh
    raphby(PROC_IN_ARRAY_m1p1)=0.0_rsh

    ALLOCATE(flx_s2w_corim1(-1:nv_adv,PROC_IN_ARRAY_m1p1))
    ALLOCATE(flx_s2w_corip1(-1:nv_adv,PROC_IN_ARRAY_m1p1))
    ALLOCATE(flx_s2w_corjm1(-1:nv_adv,PROC_IN_ARRAY_m1p1))
    ALLOCATE(flx_s2w_corjp1(-1:nv_adv,PROC_IN_ARRAY_m1p1))
    ALLOCATE(flx_w2s_corin(nvp,PROC_IN_ARRAY))
    ALLOCATE(flx_w2s_corim1(nvp,PROC_IN_ARRAY_m1p1))
    ALLOCATE(flx_w2s_corip1(nvp,PROC_IN_ARRAY_m1p1))
    ALLOCATE(flx_w2s_corjm1(nvp,PROC_IN_ARRAY_m1p1))
    ALLOCATE(flx_w2s_corjp1(nvp,PROC_IN_ARRAY_m1p1))
#if ! defined key_nofluxwat_IWS
    ALLOCATE(phieau_s2w_corim1(PROC_IN_ARRAY_m1p1))
    ALLOCATE(phieau_s2w_corip1(PROC_IN_ARRAY_m1p1))
    ALLOCATE(phieau_s2w_corjm1(PROC_IN_ARRAY_m1p1))
    ALLOCATE(phieau_s2w_corjp1(PROC_IN_ARRAY_m1p1))
#endif
  
    ! Initialization
    ALLOCATE(cini_sed(nv_state))
    DO iv = 1, nv_adv
        cini_sed(iv) = cini_sed_r(irk_fil(iv)) * unit_modif_mudbio_N2dw(irk_fil(iv))
    ENDDO   

#ifdef key_Pconstitonly_insed 
    nv_use = nvpc
#else
    nv_use = nvp
#endif  

#if ! defined key_noTSdiss_insed
    ! counting of dissolved variables for diffusion in the sediment
    ivdiss(:) = 0
    ivdiss(-1) = -1
    ivdiss(0) = 0
#if ! defined key_Pconstitonly_insed
    DO iv = 1, nv_adv-nvp
        ivdiss(iv) = iv + nvp
    ENDDO  
#endif
#endif

    !  option morpho
    !!!!!!!!!!!!!!!!!!!!
    IF(l_morphocoupl) THEN
    !! Warning : no hx, hy in other model than MARS 
        ALLOCATE(morpho0(ARRAY_morpho))
        ALLOCATE(hsed0(PROC_IN_ARRAY))
        ALLOCATE(hsed_previous(PROC_IN_ARRAY))
        hsed0(PROC_IN_ARRAY) = 0.0_rsh
        hsed_previous(PROC_IN_ARRAY) = 0.0_rsh

    ENDIF

    !! declaration of the MUSTANG variables needed in the hydro model 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ALLOCATE(EROS_FLUX_s2w(ARRAY_EROS_FLUX_s2w))
    ALLOCATE(SETTL_FLUX_w2s(ARRAY_SETTL_FLUX_w2s))
    ALLOCATE(SETTL_FLUXSUM_w2s(ARRAY_SETTL_FLUXSUM_w2s))
    EROS_FLUX_s2w(ARRAY_EROS_FLUX_s2w) = 0.0_rsh
    SETTL_FLUX_w2s(ARRAY_SETTL_FLUX_w2s) = 0.0_rsh
    SETTL_FLUXSUM_w2s(ARRAY_SETTL_FLUXSUM_w2s) = 0.0_rsh
#if ! defined key_nofluxwat_IWS && ! defined key_noTSdiss_insed
    ALLOCATE(WATER_FLUX_INPUTS(ARRAY_WATER_FLUX_INPUTS)) ! not operationnal, stil to code **TODO**
    WATER_FLUX_INPUTS(ARRAY_WATER_FLUX_INPUTS) = 0.0_rsh
#endif    
    
    ALLOCATE(fwet(PROC_IN_ARRAY))
    fwet(:,:) = 1.0_rsh
  

    END SUBROUTINE MUSTANG_alloc



#endif /* MUSTANG */

END MODULE initMUSTANG
