MODULE initMUSTANG

#ifdef key_sedim_MUSTANG

   !&E==========================================================================
   !&E                   ***  MODULE  initMUSTANG  ***
   !&E
   !&E
   !&E ** Purpose : concerns all subroutines related to sediment initialization
   !&E              to be used in MARSiAM system
   !&E 
   !&E ** Description :
   !&E     subroutine MUSTANG_param          ! defines parameterization of sediment
   !&E     subroutine MUSTANG_compatibility  ! verification of compatibility between the various parameters
   !&E
   !&E     subroutine MUSTANG_init_hsed      ! initializes the sediment thickness (key_sedim_MUSTANG_mixsed)
   !&E     subroutine MUSTANG_init_ws_sand   ! initializes the settling velocity of sand variables and number of sub time step
   !&E
   !&E ** History :
   !&E     ! 2015-12  (B.Thouvenin ) : extract from old sedim.F90  for reorganization of module SEDIMARS
   !&E
   !&E==========================================================================

#include "toolcpp.h"
#include "coupleur_define_MUSTANG.h"
#include "coupleur_dimhydro_MUSTANG.h"

   !! * Modules used

   USE comMUSTANG
   USE coupleur_MUSTANG

#ifdef key_MARS
   USE ionc4,       ONLY : ionc4_openr,ionc4_read_subxyt,ionc4_read_subzxyt,ionc4_close,  &
                           ionc4_createfile,ionc4_read_xy 
   USE comvars2d,    ONLY : iscreenlog,ierrorlog,iwarnlog,dt,ig,id
   USE comsubstance, ONLY : nvpc,igrav1,igrav2,imud1,imud2,isand1,isand2,irk_fil,ros_r, &
                            diam_r,nb_ws_max,iv_ws_max,ndt_part,l_subs2d,name_var,tocd,  &
                            dhmin,unit_modif_mudbio_N2dw,irkm_var_assoc
   USE parameters,   ONLY : grav,valmanq
   USE_MPI toolmpi,  ONLY : MASTER
#endif

   IMPLICIT NONE

   !! * Accessibility

   ! functions & routines of this module, called outside :
   PUBLIC MUSTANG_param,MUSTANG_init_ws_sand,MUSTANG_compatibility
   PUBLIC MUSTANG_init_hsed,MUSTANG_sedinit

   PRIVATE

 CONTAINS
 
   !!===========================================================================

   SUBROUTINE MUSTANG_param(rw)
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE MUSTANG_param  ***
   !&E
   !&E ** Purpose : reads parameterization of sediment dynamics
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :  init ou main, after initialisation of physical variables
   !&E                  and before sed_init and MUSTANG_alloc
   !&E
   !&E ** External calls : tool_datosec(), tool_sectodat()
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (Pierre Le Hir) extracted from init.F90
   !&E       !  2008-10  (F. Cayocca)
   !&E       !  2009-02  (P. Le Hir)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used


   !! * Arguments
   CHARACTER(LEN=1), INTENT( in ) :: rw


   !! * Local declarations
   LOGICAL               :: l_varassoc
   CHARACTER(len=lchain) :: filepc
   CHARACTER(len=19)     :: tool_sectodat
   INTEGER               :: iv,ivpc,isubs,IERR_MPI
   REAL(KIND=rlg)        :: tool_datosec

   NAMELIST/namsedim_mixsed/ l_fricwave,l_morphocoupl,fricwav,l_repsed, &
                       filrepsed,l_consolid,l_diffused,l_dzsmaxuni,     &
                       dt_sed,date_start_sed,dtcal_sed,                 &
                       choice_flxdiss_diffsed,                          &
                       l_unised,z0seduni,hseduni,cseduni,aref_sand,     &
                       ksmiuni,ksmauni,sini_sed,tini_sed,l_init_hsed,   &
                       cvolmaxsort,cvolmaxmel,activlayer,frmudcr2,      &
                       coef_frmudcr1,x1toce_mud,x2toce_mud,l_bioturb,   &
                       E0_sand_para,n_eros_sand,E0_mud,n_eros_mud,      &
                       ero_option,xexp_ero,E0_sand_option,              &
                       corfluer1,corfluer2,cfreshmud,csed_mud_ini,      &
                       dzsmin,dzsmaxuni,csedmin ,cmudcr,fexcs,          &
                       xperm1,xperm2,xsigma1,xsigma2,fimpls,            &
                       xdifs1,xdifs2,xdifsi1,xdifsi2,epdifi,MF,         &
                       htncrit_eros,slopefac,coef_erolat,coef_tenfon_lat,&
                       l_erolat_wet_cell,l_bathy_smoothing,             &
                       csegreg,csandseg,l_initsed_vardiss,              &
                       xbioturbmax,xbioturbk,dbiotu0,dbiotum,           &
                       frmud_db_min,frmud_db_max
                       

   NAMELIST/namsedoutput/ name_out_eptot,name_out_nbniv,name_out_dzs,name_out_tenfon,     &
                          name_out_tenfonc,name_out_tenfonw,choice_nivsed_out,            &
                          riog_valid_min_eptot,riog_valid_max_eptot,                      &
                          riog_valid_min_nbniv,riog_valid_max_nbniv,riog_valid_min_dzs,   &
                          riog_valid_max_dzs,riog_valid_min_tenfon,riog_valid_max_tenfon, &
                          nk_nivsed_out,ep_nivsed_out,epmax_nivsed_out,l_outsed_saltemp
!cval variable restee locale                          riog_valid_min_cvsed,riog_valid_max_cvsed


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

IF(rw == 'r')THEN

! Lecture des namelists
!==========================

    filepc='./parasedim.txt'
    OPEN(50,file=filepc,status='old',form='formatted',access='sequential')
    READ(50,namsedim_mixsed)
    xsigma1sg=xsigma1/grav
    IF(csed_mud_ini==0.0_rsh)csed_mud_ini=cfreshmud
    tstart_sed=tool_datosec(date_start_sed)
    IF(.NOT.l_morphocoupl) MF=1.0_rsh
    cimpls=1.0_rsh-fimpls
    cexcs=1.0_rsh-fexcs
    fws2=fricwav*0.5_rsh
    l_dyn_insed=.FALSE.
    IF(l_consolid .OR. l_bioturb .OR. l_diffused)l_dyn_insed=.TRUE.

    READ(50,namsedoutput)
    CLOSE(50) 

    ! conversion mm lus in 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

    ! fin de lecture de la namelist

ELSE
! ecriture des namelists
!==========================
     IF_MPI (MASTER) THEN
    
        WRITE(iscreenlog,namsedim_mixsed)

    ENDIF_MPI
    
ENDIF
     PRINT_DBG*, 'FIN MUSTANG_param'   
   
  END SUBROUTINE MUSTANG_param

  !!===========================================================================
 
   SUBROUTINE MUSTANG_compatibility
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE MUSTANG_compatibility  ***
   !&E
   !&E ** Purpose : verifies compatibility between the various parameters (sedim)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : cas1DV apres cas_init_sedim
   !&E                init_param
   !&E
   !&E ** External calls : tool_datosec(), tool_sectodat()
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2015 (B. Thouvenin)
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used


   !! * Arguments

   !! * Local declarations
   INTEGER               :: iv,ivpc,isubs,IERR_MPI,iv2
   REAL(rsh)             :: rohomog
 
    ! test MUSTANG + key_MARS
#if ! defined key_MARS
            WRITE(ierrorlog,*)
            WRITE(ierrorlog,*)'****************************'
            WRITE(ierrorlog,*)' You are using key_MUSTANG with MARS3D model '
            WRITE(ierrorlog,*)' you must set CPP key : key_MARS in Makefile_datarmor'
            WRITE(ierrorlog,*)' It is not the case : the model is stopped'

            CALL_MPI MPI_FINALIZE(IERR_MPI)
            STOP
#endif
    ! test dissolved variables in sediment must be used without  key_nofluxwat_IWS
#if defined key_nofluxwat_IWS && ! defined key_nodiss_insed
       WRITE(iwarnlog,*)
       WRITE(iwarnlog,*)'****************************'
       WRITE(iwarnlog,*)'        ATTENTION           ' 
       WRITE(iwarnlog,*)'****************************'
       WRITE(iwarnlog,*)' You are using CPP key : key_nofluxwat_IWS ' 
       WRITE(iwarnlog,*)' So you are not taking into account water fluxes threw water/sediment interface ' 
       WRITE(iwarnlog,*)' which are generated by erosion/ settling/ consolidation..'
       WRITE(iwarnlog,*)' AND you are NOT using CPP key : key_nodiss_insed ' 
       WRITE(iwarnlog,*)' so you want to simulate dissolved variables in sediment'
       WRITE(iwarnlog,*)' You will have probably conservativity problems '
       WRITE(iwarnlog,*)' If you simulate dissolved variables in water AND in sediment, don t use key_nofluxwat_IWS '
#endif


    ! test compatibilite tocd /consolidation
    IF_MPI (MASTER) THEN
    IF (l_consolid) THEN
       WRITE(iwarnlog,*)
       WRITE(iwarnlog,*)'****************************'
       WRITE(iwarnlog,*)' You are taking into account CONSOLIDATION PROCESS'
       WRITE(iwarnlog,*)' and fresh deposit concentration (cfreshmud) = ',cfreshmud
       WRITE(iwarnlog,*)' and critical stress of deposition (tocd) of each sand and mud variable are equal to :'
       DO isubs=imud1,imud2
         WRITE(iwarnlog,*)' variable',name_var(irk_fil(isubs)),'  tocd=',tocd(isubs)
       ENDDO
       WRITE(iwarnlog,*)' If cfreshmud is big (> 100), you have to choose a tocd between 1 and 10'
       WRITE(iwarnlog,*)' otherwise you have to choose a tocd > 10 or 20'
       
    ELSE
       DO isubs=isand2+1,nvp
         IF (tocd(isubs) > 1.5 ) THEN
           WRITE(iwarnlog,*)'****************************'
           WRITE(iwarnlog,*)'  WARNING    '
           WRITE(iwarnlog,*)' You are not taking into account CONSOLIDATION PROCESS'
           WRITE(iwarnlog,*)' fresh deposit concentration (cfreshmud) = ',cfreshmud
           WRITE(iwarnlog,*)' and critical stress of deposition (tocd) of isusb=',isubs,' exceed 1 Pa :',tocd(isubs)
           WRITE(iwarnlog,*)' it should be smaller (0.5 to 1 ?) '
         !  CALL_MPI MPI_FINALIZE(IERR_MPI)
         !  STOP
         ENDIF
       ENDDO
    ENDIF
    ENDIF_MPI
    
    ! test ros homogene pour les sables dans tous les cas
    ! test ros homogene pour toutes les particules si consolidation

    IF_MPI (MASTER) THEN
    IF (l_consolid) THEN
       rohomog=ros_r(irk_fil(1))
       DO iv=2,nvpc
        IF(ros_r(irk_fil(iv)) .NE. rohomog) THEN
            WRITE(ierrorlog,*)
            WRITE(ierrorlog,*)'****************************'
            WRITE(ierrorlog,*)' You are taking into account CONSOLIDATION PROCESS'
            WRITE(ierrorlog,*)' and you must have same density for all constitutive particles'
            WRITE(ierrorlog,*)' It is not the case : (see in variable.dat) '
            WRITE(ierrorlog,*)' ros(1)= ',ros_r(irk_fil(1))
            WRITE(ierrorlog,*)' ros(',irk_fil(iv),')= ',ros_r(irk_fil(iv))

            CALL_MPI MPI_FINALIZE(IERR_MPI)
            STOP
         ENDIF
       ENDDO
     ELSE
       IF(isand2 > isand1) THEN
        rohomog=ros_r(irk_fil(isand1))
        DO iv=isand1+1,isand2
         IF(ros_r(irk_fil(iv)) .NE. rohomog) THEN
            WRITE(ierrorlog,*)
            WRITE(ierrorlog,*)'****************************'
            WRITE(ierrorlog,*)' you must have same density for all sand particles'
            WRITE(ierrorlog,*)' It is not the case : (see in variable.dat) '
            WRITE(ierrorlog,*)' ros(1)= ',ros_r(irk_fil(isand1))
            WRITE(ierrorlog,*)' ros(',irk_fil(iv),')= ',ros_r(irk_fil(iv))

            CALL_MPI MPI_FINALIZE(IERR_MPI)
            STOP
          ENDIF
          DO iv2=isand1,iv
           IF(diam_r(irk_fil(iv)) > diam_r(irk_fil(iv2))) THEN
             WRITE(ierrorlog,*)
             WRITE(ierrorlog,*)'****************************'
             WRITE(ierrorlog,*)' The "SAND" variables should be stored in variable.dat '
             WRITE(ierrorlog,*)'in order of decreasing diameters (the coarsest sand to the finest sand)'
             WRITE(ierrorlog,*)' It is not the case : (see in variable.dat) '
             WRITE(ierrorlog,*)' diam_r(',irk_fil(iv),')=',diam_r(irk_fil(iv)),' > diam_r(',irk_fil(iv2),')=',diam_r(irk_fil(iv2))

             CALL_MPI MPI_FINALIZE(IERR_MPI)
             STOP
           ENDIF
          ENDDO
        ENDDO
       ENDIF
     ENDIF

    ENDIF_MPI


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

    PRINT_DBG*, 'FIN MUSTANG_compatibility'   
   
  END SUBROUTINE MUSTANG_compatibility

  !!===========================================================================
 
   !!===========================================================================
 
  SUBROUTINE MUSTANG_init_ws_sand(rhoref)
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE MUSTANG_init_ws_sand  ***
   !&E
   !&E ** Purpose : initialize sand settlings velocities and number of sub time step
   !&E
   !&E ** Description : 
   !&E
   !&E ** Called by :  sed_init, MUSTANG_init_fromfile and cas_init_sedim (if testcase)
   !&E 
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2015-04  (B. Thouvenin) extracted from sed_init et sed_initfromfile in sedim.F90
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(KIND=rsh),INTENT(IN)          ::   rhoref

   !! * Local declarations
   INTEGER        :: iv,ivg                

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

! valable aussi pour les graviers ???
    DO iv=1,isand2

      ! here, for simplification, water density = rhoref 
      ! here, for simplifcation, 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*(grav*(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)
      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)*grav*(ros(iv)-rhoref)*diam_sed(iv)
      xnielsen(iv)=12.0_rsh*1000.0_rsh/grav/(1.0_rsh-rhoref/ros(iv))
    ENDDO
    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

#ifdef key_MARS
! reinitialization of ndt_part for sand variables
   DO ivg=1,nb_ws_max
     iv=iv_ws_max(ivg,1)
     IF ( iv<=isand2) THEN
         IF(l_subs2D(iv)) THEN
           ndt_part(ivg)=1
         ELSE
           IF(ws_sand(iv) .GT. 0.0_rsh)THEN
             ndt_part(ivg)=NINT(dt/MIN(dhmin/ws_sand(iv),dt))
           ELSE
             ndt_part(ivg)=1
             IF_MPI (MASTER) THEN
                write(iscreenlog,*)'VARIABLE SAND :',iv,'  : settling velocity =0 ==> It is not normal'
             ENDIF_MPI
           ENDIF
         ENDIF
         IF_MPI (MASTER) THEN
           write(iscreenlog,*)'VARIABLE : ',iv,'  : number of sub-time step  =',ndt_part(ivg)
         ENDIF_MPI

     ENDIF
   ENDDO
#endif

    PRINT_DBG*, 'FIN MUSTANG_INIT_WS_SAND'

  END SUBROUTINE MUSTANG_init_ws_sand

   !!===========================================================================
!
  SUBROUTINE MUSTANG_init_hsed(cv_sedini)
!
!     rinitialize sediment thickness, in order to be coherent with sediment parameters
!     function of cseduni, cvolmax, csed_ini of each sediment, 
! **********************************************************************
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used

   !! * 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

    
      !!! 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,rhoref,h0)
 
   !&E--------------------------------------------------------------------------
   !&E                 ***  ROUTINE sed_init  ***
   !&E
   !&E ** Purpose : initialize the sediment geometry & content
   !&E
   !&E ** Description : 
   !&E
   !&E ** Called by :  main, after reading MUSTANG_param
   !&E                 only if not.l_repsed
   !&E 
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (Pierre Le Hir) extracted from subsedim.F90 (from F. Dufois)
   !&E       !  2009-02  (P. Le Hir, F. Cayocca) modified for new sand/mud sediment 
   !&E       !  2015-04  (F. Grasso) modified initialization at -valmanq instead of valmanq
   !&E
   !&E--------------------------------------------------------------------------
   !! * Modules used


   !! * Arguments
   INTEGER,INTENT(IN)                 ::   ifirst,ilast,jfirst,jlast
   REAL(KIND=rsh),INTENT(IN)          ::   rhoref
   REAL(KIND=rsh),DIMENSION(PROC_IN_ARRAY_m2p2),INTENT(IN)       :: 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
   REAL(KIND=rsh),DIMENSION(5):: cini_sed_modif
   REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax)                        :: frac_grav_ini,frac_sand_ini,frac_mud_ini
   REAL(KIND=rsh),DIMENSION(1:5,imin:imax,jmin:jmax)                    :: fracini
   REAL(KIND=rsh),DIMENSION(5)                                          :: c_sed_ini
   CHARACTER(LEN=lchain)                                                :: file_fracini


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

    CALL MUSTANG_init_ws_sand(rhoref)

    dzsmax(PROC_IN_ARRAY)=dzsmaxuni
    z0sed(PROC_IN_ARRAY)=z0seduni
    ksmi(PROC_IN_ARRAY)=ksmiuni
    ksma(PROC_IN_ARRAY)=0
    hsed(PROC_IN_ARRAY_m1p2)=-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 (h0(PROC_IN_ARRAY) /= -valmanq) ksmi(PROC_IN_ARRAY)=ksmiuni
      WHERE (h0(PROC_IN_ARRAY) /= -valmanq) ksma(PROC_IN_ARRAY)=ksmauni
      DO k=ksdmin,ksdmax
        WHERE (h0(PROC_IN_ARRAY) /= -valmanq) dzs(k,PROC_IN_ARRAY)=0.0_rsh
        WHERE (h0(PROC_IN_ARRAY) /= -valmanq) c_sedtot(k,PROC_IN_ARRAY)=0.0_rsh
        DO iv=-1,nv_tot
          WHERE (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(tot_csed_init .NE. 1.0_rsh) THEN
        WRITE(ierrorlog,*)' '
        WRITE(ierrorlog,*)' inital concentrations of constitutive particulate matter'
        WRITE(ierrorlog,*)' have not been well defined'
        WRITE(ierrorlog,*)' See variable.dat to give fraction of total concentration '
        WRITE(ierrorlog,*)' for each constitutive variable (GRAV, SAND, MUD) (but not MUDB) '
        WRITE(ierrorlog,*)' : the sum of all must be =1 '
        CALL_MPI MPI_FINALIZE(IERR_MPI)
        STOP         
      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)
         IF_MPI (MASTER) THEN
           WRITE(iscreenlog,*)'NEW TOTAL SEDIMENT THICKNESS AFTER ADJUSTEMENT = ',hsed_new
         ENDIF_MPI
         cseduni=cseduni*hseduni/hsed_new         
         hseduni=hsed_new
         WHERE (h0(PROC_IN_ARRAY_0p1) /= -valmanq) hsed(PROC_IN_ARRAY_0p1)=hseduni
         DEALLOCATE(cv_sedini)
      ELSE
         WHERE (h0(PROC_IN_ARRAY_0p1) /= -valmanq) hsed(PROC_IN_ARRAY_0p1)=hseduni
      ENDIF

      DO j=jfirst,jlast
#ifdef key_MARS
        DO i=MAX0(ifirst,ig(j)),MIN0(ilast,id(j))
#else
        DO i=ifirst,ilast
#endif
          IF(h0(i,j) /= -valmanq) THEN
            DO k=ksmi(i,j),ksma(i,j)
              somalp=0.
              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)
                somalp=somalp+cv_sed(iv,k,i,j)/ros(iv)
              ENDDO
              poro(k,i,j)=1.0_rsh-somalp
#ifdef key_sedim_partial
              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 : modif FG(15/12/2020)   
      ! -----------------------

      PRINT*, 'NON UNIFORM BED COVERAGE'


      ! uniform initialisation
      WHERE (h0(PROC_IN_ARRAY) /= -valmanq) ksmi(PROC_IN_ARRAY)=ksmiuni
      WHERE (h0(PROC_IN_ARRAY) /= -valmanq) ksma(PROC_IN_ARRAY)=ksmauni
      DO k=ksdmin,ksdmax
        WHERE (h0(PROC_IN_ARRAY) /= -valmanq) dzs(k,PROC_IN_ARRAY)=0.0_rsh
        WHERE (h0(PROC_IN_ARRAY) /= -valmanq) c_sedtot(k,PROC_IN_ARRAY)=0.0_rsh
        DO iv=-1,nv_tot
          WHERE (h0(PROC_IN_ARRAY) /= -valmanq) cv_sed(iv,k,PROC_IN_ARRAY)=0.0_rsh
        ENDDO
      ENDDO

! fichier netcdf ou sont sont lus les pourcentages de chaque classe pour chaque maille du domaine (base facies MANGA du SHOM)
      file_fracini='../inputs/SED/fr_sed_ini.nc'

      CALL ionc4_openr(file_fracini)
      CALL ionc4_read_xy(file_fracini,"fr_grav_ini",frac_grav_ini,imin,imax,jmin,jmax)
      CALL ionc4_read_xy(file_fracini,"fr_sand_ini",frac_sand_ini,imin,imax,jmin,jmax)
      CALL ionc4_read_xy(file_fracini,"fr_mud_ini",frac_mud_ini,imin,imax,jmin,jmax)
      CALL ionc4_close(file_fracini)
      print*,'File fr_sed_ini read'


      DO j=jfirst,jlast
#ifdef key_MARS
        DO i=MAX0(ifirst,ig(j)),MIN0(ilast,id(j))
#else
        DO i=ifirst,ilast
#endif

!!! initialisation des fractions initiales
! 1 gravel
fracini(1,i,j)=frac_grav_ini(i,j)
c_sed_ini(1)=1500.0_rsh
! 3 sands
fracini(2:4,i,j)=frac_sand_ini(i,j)/3
c_sed_ini(2:4)=1500.0_rsh
! 1 mud
fracini(5,i,j)=frac_mud_ini(i,j)
c_sed_ini(5)=500.0_rsh

          DO k=ksmi(i,j),ksma(i,j)
            somalp=0.
            DO iv=1,nvpc
            ! - constitutive variables : cini_sed as percentage of sediment mass
              !cv_sed(iv,k,i,j)=cini_sed(iv)*c_sedtot(k,i,j)
              IF (fracini(iv,i,j) /= -valmanq) THEN
                cv_sed(iv,k,i,j)=fracini(iv,i,j)*c_sed_ini(iv)
              ELSE
                cv_sed(iv,k,i,j)=0.0
              END IF
              c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
              somalp=somalp+cv_sed(iv,k,i,j)/ros(iv)
            ENDDO
            poro(k,i,j)=1.-somalp
            cv_sed(-1,k,i,j)=tini_sed
            cv_sed( 0,k,i,j)=sini_sed
          ENDDO

!  reinitializes sediment thickness
    !!! Concentration massique
       IF(l_init_hsed) THEN
          ALLOCATE(cv_sedini(1:nvpc))
          DO iv=1,nvpc
         ! constitutive variables : cini_sed in % of cseduni
           cv_sedini(iv)=0.0_rsh
           DO k=ksmi(i,j),ksma(i,j)
              cv_sedini(iv)=cv_sedini(iv)+cv_sed(iv,k,i,j)
           ENDDO
          ENDDO
          CALL MUSTANG_init_hsed(cv_sedini)
          IF_MPI (MASTER) THEN
            WRITE(iscreenlog,*)'NEW TOTAL SEDIMENT THICKNESS AFTER ADJUSTEMENT = ',hsed_new
          ENDIF_MPI
          hsed(i,j)=hsed_new
          IF (hsed(i,j) == 0.0_rsh) ksma(i,j)=0
          DEALLOCATE(cv_sedini)
       ENDIF

        ENDDO
      ENDDO   !end of i,j loops

    ENDIF
!
!   FIN COUVERTURE SEDIMENTAIRE NON UNIFORME 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_MPI (MASTER) THEN
            IF(dzsmaxuni == 0.0_rsh) THEN
               write(ierrorlog,*)' ERROR in parasedim : dzsmaxuni =0. and l_repsed=.FALSE.'
               write(ierrorlog,*)' if dzsmaxuni =0. dzsmax must be read in an initial file (l_repsed=T)'
               write(ierrorlog,*)' if l_repsed=.FALSE. dzsmax must be evaluated from dzsmaxuni not null'

               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
            ENDIF
       ENDIF_MPI
       ! definition des dzsmax
       ! estime en fonction de la profondeur et de dzsmaxuni considere comme la valeur max de dzsmax
       ! ici valeur minimum de dzsmax en profondeur=dzsmaxuni/100
       ! profondeur max donnant la plus faible valeur de dzsma a informer, mise par defaut=50m
       dzsmaxmin=dzsmaxuni/100._rsh
       h0max_def=50.0_rsh
       write(iscreenlog,*)
       write(iscreenlog,*)'DZSMAX NON UNIFORM'
       write(iscreenlog,*)'by default linearly computed in MUSTANG_sedinit (sedim.F90)'
       write(iscreenlog,*)'from dzsmaxuni (',dzsmaxuni,'m) to dzsmaxuni/100 (',dzsmaxmin,'m) depending on water depth'
       write(iscreenlog,*)'dzsmax minimum for depth > ',h0max_def
       write(iscreenlog,*)'PROGRAM your own initialization of dzsmax (or READ a file) in MUSTANG_sedinit'
       dzsmax(PROC_IN_ARRAY)=(h0max_def-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

       ! ou lus dans un fichier (a programmer)
        
    ENDIF
    
    
    PRINT_DBG*, 'FIN sed_MUSTANG_sedINIT'

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

#endif

END MODULE initMUSTANG
