  SUBROUTINE cas_estuar
#ifdef key_castest_estuary
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_estuar  ***
   !&E
   !&E ** Purpose : definition of a test case for a westward open channel
   !&E              containing sandy and muddy sediments and a tidal flat 
   !&E              on the eastward boundary
   !&E
   !&E ** Description : several possibilities according to the key selected, 
   !&E                  determining the number and quality of variables  :
   !&E       only key_sedim_MUSTANG (always_mixsed ) : 2 sand and 1 mud variable   !&E
   !&E       key_sedim_udsand no more operationnal with MUSTANG (since r2276)
   !&E                      + key_castestriver                           : with a river upstream and outflow 
   !&E
   !&E       several possiblities for BATHYMETRY
   !&E                         - key_castestestuary_bathy1 : printemps 2014 (Stage Clement Bouvier)
   !&E                         - key_castestestuary_bathy2 : ete  2014 (C3E2) 
   !&E                         - key_castestestuary_bathy3 : hiver 2015 (DCE)
   !&E
   !&E
   !&E ** Called by : cas_init
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2011-03 (F. Cayocca) new test case issued from test caschansed
   !&E       !  2011-11 (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2013-11 (S. Petton) Update of the code for new version
   !&E       !  2016-02 (B. Thouvenin, F. Cayocca) new test cas (estuar) from castidalflat2DV
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : rsh,rlg,riolg,valmanq,kmax,imin,imax,lchain,&
                          ljminm1,ljmaxp2,liminm1,limaxp2,liminp1,epsilon, &
                          liminm2,ljminm2,pi,limin,limax,ljmin,ljmax,limaxp1,ljmaxp1, &
                          obc_width,ishift_obc_west
   USE comvars3d,  ONLY : nzinit,kzinit,qmaxz,l_stateeq_lin,saliref_lin,tetaref_lin
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : theta_sig,b_sig,hc
#endif
   USE casinit,   ONLY : cas_init_param,cas_init_grid,cas_init_bathy, &
                         cas_init_io,cas_init_obc,cas_init_meteo,     &
                         cas_init_visc,cas_init_turb,cas_init_tspar,  &
                         cas_init_progvar,cas_init_misc,cas_init_river,cas_init_river
   USE output,    ONLY : pasor,l_out_sal,l_out_temp,fileresu,tdebsor,tfinsor, &
                         iextmin,iextmax,jextmin,jextmax,l_out_ssh,l_out_ubt,l_out_u3d,l_out_wind
   USE toolgeom,  ONLY : f_zu
   USE comvarp2d, ONLY : temp,sal,ssh
   USE comvars3d, ONLY : tetaref,saliref
   USE comvars2d,  ONLY : fisac,gwac,dtmin,dtmax,dt,suffix,fvisc,kx,ky,qmax,iscreenlog,dhx_g,dhy_g, &
                          dhx,dhy,h0_g,hx_g,hy_g,h0fond,hm,vis2d,vish_xe,vish_phi,dhxy_g,dhxy,halfdt, &
                          hrdef,dhj,hx,hy,dfiac,dgac,ig,id,jb,jh,l_initfromfile,l_saverestart_bydate, &
                          saverestart_step,l_restart_subs,file_init,l_init_rtemp,l_init_rsal,fwetmin, &
                          ierrorlog,vismin,vismax,cosmag,l_smagor,l_spongw,sponge_nbcell,icon,fu,fv,h0, &
                          tobs,out_dir,hminfrot,h0_g,hx_g,hy_g,h0fond,hm,dhj,hrdef,vis2d,vish_xe,vish_phi,z0b

   USE obccombine,ONLY : l_obc_mars,l_obc_west,l_obc_tide,l_obc_ogcm
   USE obcogcm,   ONLY : obc_sw,obc_tw,obc_coefrel
   USE comaree
   USE comvars3d, ONLY : turb_0eq_option,turb_nbeq,turb_2eq_option,bgdiff,OptionSF
   USE comdiag, ONLY : l_diag,file_diag,tdebdiag,tfindiag
#if defined key_substance
   USE comsubstance
   USE casinit,   ONLY : cas_init_subs,cas_init_subs_disspart,cas_init_outflow,cas_init_subs_budget
#endif
#ifdef key_sedim_MUSTANG
   USE parameters, ONLY : ksdmax
   USE comMUSTANG
   USE casinit,   ONLY : cas_init_sedim
   USE initMUSTANG,  ONLY : MUSTANG_compatibility
   USE sedimorpho, ONLY : sed_morphoinit,sed_alloc_morpho
   USE sed_MUSTANG_MARS,  ONLY :    sedinit_fromfile
#endif
#if defined key_castestriver
   USE comfluxdo,  ONLY : nbfl,iriv,jriv,kriv,debmoy,tmfle,amtfle,l_tempriv,l_flestat,river_thickness,deb_ampli_var                                     
#endif

   IMPLICIT NONE

   !! * Local declarations
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   CHARACTER(LEN=19) :: date_startout,date_endout
   CHARACTER(LEN=lchain) :: fileresutmp 
   INTEGER           :: i,j,k,isubs,iv,lnam
   INTEGER           :: ni_wall_e,ni_wall_w,nj_wall_n,nj_wall_s
   INTEGER           :: IERR_MPI
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst,zr,sinit_sed,tinit_sed,z0sedinit,dzsmaxinit
   REAL(KIND=rsh)    :: hshallow,hdeep,concssmriver,nu0
   REAL(KIND=rlg)    :: tool_datosec
   LOGICAL           :: l_filesubs,l_two_layers,l_nc4par
   LOGICAL           :: l_csedinit_varz,l_hsedinit_cste,l_casinit_hsed,l_dzsinit_varz
   INTEGER           :: jmilieu,i_mouth,jtempo,iterloc
   REAL(KIND=rsh),DIMENSION(imin:imax)    ::xfond
   REAL(KIND=rsh),DIMENSION(imin:imax)    ::yfond,yfondx,l1,l2,l3,l1y,l2y,l3y,               &
                                               Smm,hchenal,hpm,hbm,hpmshbm2
   REAL(KIND=rsh)    :: incr,dxval,dyval,x,m,largchenal_embouch,varexp_section,hlocchenal,   &
                        varexp_largchenal,largchenal_min,cotchenal_min,cotchenal_embouch,    &
                        pentechenal,section_embouch,gamma,gamma2,cotchenal_amont,l2i,l2n,    &
                        varexp_sect_off,l2ml1m,l3ml1m,cotfond_off  
   REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: hx_t,hy_t

#ifdef key_castestestuary_bathy3
   REAL(KIND=rsh)    :: lgcvg_section,lgcvg_largeur,lgcvg_cotchenal,larbm_embouch,           &
                        deltacorsect,corsectmax,corsectmin,corsect,l4sl2,heaubm,denom
   REAL(KIND=rsh),DIMENSION(imin:imax)    :: l4,l4y,sectmm
#elif key_castestestuary_bathy1
   REAL(KIND=rsh),DIMENSION(imin:imax)    :: l4
#endif


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

   PRINT_DBG*, 'enter CAS_ESTUAR routine'

   ! general parametrization (each parameter can be reassigned)
   ! -----------------------

    ! 1- general parametrization (each parameter can be reassigned)
   ! ==========================
   CALL cas_init_param
        ! over-assignment :
        h0fond=0.01_rsh   ! the cell bottom is h0fond lower than any of the edges (ref.: 0.)
        hm=1.00_rsh       ! transition 2D/3D (ref.: 0.0 ?)
        fwetmin=1.0_rsh   ! uncovered cell fraction : if =1, the cell cannot be partially uncovered
        dhj=0.2           ! additional water height for exchanges between cells   !PLH 0.1 avant
        hrdef=20.0        ! bottom altitude unnderwhich dhj ignored


   ! 2- Water depth  (Bathymetry)
   ! ===========================
   ! assignement : input variables for cas_init_bathy routine
      nj_wall_s = 2  ! land over nj_wall_s cells (from jmin row) at the southern part of the domain
      nj_wall_n = 1  ! land over nj_wall_n cells at the northern part of the domain
      ni_wall_e = 1
   ! provisory call for  cas_init_bathy :
      ! for definition of a 2DV slope ??? :
      hshallow=-1.0_rsh
      hdeep=16.0_rsh
      hx_g(imin,2)=hdeep
      h0_g(imin,2)=hdeep+h0fond
      DO i=imin+1,imax
        hx_g(i,2)=hdeep-REAL(i-imin-1,rsh)/REAL(imax-imin-1,rsh)*(hdeep-hshallow)
      END DO
      DO i=imin+1,imax
        h0_g(i,2)=MAX(hx_g(i,2),hx_g(i-1,2))+h0fond
      END DO
      depthcst=0.0_rsh  ! should be 0 to be ignored in cas_init_bathy

   CALL cas_init_bathy(depthcst,nj_wall_n=nj_wall_n,nj_wall_s=nj_wall_s,ni_wall_e=ni_wall_e)      

    ! for thix test case : specific bathymetry definition after grid parametrization

   ! 3- Grid parametrization in the water column
   ! ===========================================
   ! assignement : input variables for cas_init_grid routine
#ifdef key_castestestuary_bathy2
   dxval=600.0_rsh
#else
   dxval=300.0_rsh
#endif
   dyval=100.0_rsh
   fisac=0.0_riolg        ! lattitude (for Coriolis estimation ?)
   gwac=0.0_riolg         ! longitude (for Coriolis estimation ?)
   CALL cas_init_grid(dxval,dyval)

   ! over-assignement
        ! over-assignment
        fu(:,:)=0.0_rsh   ! no Coriolis effect
        fv(:,:)=0.0_rsh   ! no Coriolis effect 
#if defined key_siggen || defined key_gencoord
        theta_sig= 4.0_rsh ! entre 0 et 20
        b_sig= 0.750_rsh   ! entre 0 et 1    
        hc = 5.0_rsh
#endif

   ! 2b- Water depth  (Bathymetry)
   ! ===========================
   ! Schematic Estuarine bathymetry (Clement Bouvier / Pierre Le Hir)
   ! ----------------------------------------------------------------

      i_mouth=imin+40
      xfond(i_mouth)=0_rsh
      jmilieu=int((jmax-jmin)/2)
      m=2.0                           !demi marnage

      ! Definition de xfond
      ! -------------------
     DO i=i_mouth+1,imax
       xfond(i)=xfond(i-1)+dxval
     ENDDO
     DO i=i_mouth-1,imin,-1
       xfond(i)=xfond(i+1)-dxval
     ENDDO
  
     ! Definition de yfond
     ! -------------------
     yfond(jmilieu)=dyval/2.0_rsh
     yfondx(jmilieu)=0._rsh
     DO j=jmilieu+1,jmax
       yfond(j)=yfond(j-1)+dyval
       yfondx(j)=yfondx(j-1)+dyval
     ENDDO

   ! Parametres de la geometrie schematique
   ! --------------------------------------

#ifdef key_castestestuary_bathy1

!   1) calcul printemps 2014 (Stage Clement Bouvier) 
!   -----------------------------------------------

     largchenal_embouch=200_rsh        ! tests 1 à 4 : 205_rsh
     varexp_section=0.0000427_rsh
     varexp_largchenal=0.00002135_rsh
     largchenal_min=50_rsh             ! tests 1 à 4 : 75_rsh
     cotchenal_min=0.1_rsh+m           ! tests 1 à 4 : 5.1_rsh+m            !1.0_rsh+m
     cotchenal_embouch=6.0_rsh+m       ! tests 1 à 4 : 9.31_rsh+m       !10.0_rsh+m
     pentechenal=0.0001_rsh            ! tests 1 à 4 : 0.0002_rsh
     section_embouch=14740.0_rsh
     varexp_sect_off=1.0_rsh/3000.0_rsh
     cotfond_off=20.0_rsh
     
     DO i=i_mouth,imax
        ! Definition de l1
        l1(i)=amax1(largchenal_min,largchenal_embouch*exp(-varexp_largchenal*xfond(i)))
        ! Definition de la hauteur d'eau au fond (hl1)
        hx_g(i,jmilieu)=amax1(cotchenal_min,cotchenal_embouch-pentechenal*xfond(i))           
        ! Definition de la section a mi maree
        Smm(i)=section_embouch*exp(-varexp_section*xfond(i))
        ! Definition de l2
        hbm(i)=hx_g(i,jmilieu)-m               !hbm=m
        hpm(i)=hbm(i)+2.0_rsh*m                !hpm=-m
        hpmshbm2(i)=(hpm(i)/hbm(i))**2              
          ! l2(i)=(4.0_rsh*Smm(i)-l1(i)*(hx_g(i,jmilieu)-m)*(3.0_rsh+hpmshbm2))/ &
          !       (4.0_rsh*hx_g(i,jmilieu)+(hx_g(i,jmilieu)-m)*(hpmshbm2-1.0_rsh))
        l2(i)=(4.0_rsh*Smm(i)-l1(i)*hx_g(i,jmilieu)*(3.0_rsh+hpmshbm2(i)))/    &
                (4.0_rsh*(hx_g(i,jmilieu)+m)+hx_g(i,jmilieu)*(hpmshbm2(i)-1.0_rsh))
        ! si Yates  l2(i)=(Smm(i)-l1(i)*hx_g(i,jmilieu))/(hx_g(i,jmilieu)+4_rsh*m)    
        l2(i)=0.011+amax1(l1(i),l2(i))
        ! Definition de l3
        l3(i)=(2.0_rsh*Smm(i)-3.0_rsh*m*l2(i)-2.0_rsh*hx_g(i,jmilieu)*(l1(i)+l2(i)))/m
        ! si Yates  l3(i)=5*l3
        l3(i)=0.01+amax1(l2(i),l3(i))
        !creation de l4 pour adoucir la pente
        l4(i)=400.0_rsh+l3(i)
     ENDDO 
#endif

#ifdef key_castestestuary_bathy2

!   2) calcul ete  2014 (C3E2)  
!   -------------------------
    
     cotchenal_amont=2.5_rsh
     cotchenal_embouch=10.0_rsh
     section_embouch=10000.0_rsh
     varexp_section=1.0_rsh/24000.0_rsh      !1.0_rsh/12000.0_rsh
     largchenal_embouch=200_rsh              !100_rsh
     varexp_largchenal=1.0_rsh/36000.0_rsh   !1.0_rsh/18000.0_rsh
     gamma=1.1_rsh
     varexp_sect_off=1.0_rsh/3000.0_rsh
     cotfond_off=20.0_rsh
    
     pentechenal=(cotchenal_amont-cotchenal_embouch)/xfond(imax)
     gamma2=gamma**2
     DO i=i_mouth,imax
          ! Definition de l1
          l1(i)=largchenal_embouch*exp(-varexp_largchenal*xfond(i))
          ! Definition de la hauteur d'eau au fond (hl1)
          hx_g(i,jmilieu)=cotchenal_embouch+pentechenal*xfond(i)           
          ! Definition de la section a mi maree
          Smm(i)=section_embouch*exp(-varexp_section*xfond(i))
          ! Definition de l2
          l2i=100.
          hlocchenal=hx_g(i,jmilieu)-m
          DO iterloc = 1,1000
             l2n=(Smm(i)-l1(i)*hlocchenal)/(hlocchenal+1.5_rsh*m+0.5_rsh*m*(1.0_rsh   &
                  +4.0_rsh*(m/hlocchenal)*l2i/(l1(i)+l2i))**2/gamma2)
                      !!  IF (abs(l2n-l2i).le.0.5) THEN
                      !!    PRINT*,'pour i = ',i,' , iterloc = ',iterloc
                      !!    exit
                      !!  ENDIF
             l2i=l2n
          ENDDO
          l2(i)=l2n
          IF (l2(i).ge.l1(i)) THEN
            l3(i)=l2(i)*(1.0_rsh+4.0_rsh*(m/hlocchenal)*l2n/(l1(i)+l2n))**2/gamma2
          ELSE
            l2(i)=l1(i)
            l3(i)=(Smm(i)-hlocchenal*(l1(i)+l2(i))-2.0_rsh*m*l2(i))*2.0_rsh/m-l2(i)
          ENDIF
     ENDDO
#endif


#ifdef key_castestestuary_bathy3   

!   3) calcul hiver 2015 (DCE)   
!   -------------------------  

     cotchenal_embouch=10.0_rsh
     cotchenal_amont=-3.0_rsh
     section_embouch=10000.0_rsh
     lgcvg_section=8000.0_rsh
     larbm_embouch=600_rsh
     l4sl2=2.5_rsh
     lgcvg_largeur=10000.0_rsh
     deltacorsect=0.0
     corsectmin=0.7
     corsectmax=1.3
     gamma=1.1  
     gamma2=gamma**2
     !lgcvg_cotchenal=lgcvg_largbm*lgcvg_section/(0.1+lgcvg_largbm-lgcvg_section)
     varexp_sect_off=1.0_rsh/3000.0_rsh
     cotfond_off=20.0_rsh

     pentechenal=(cotchenal_amont-cotchenal_embouch)/xfond(imax)
     corsect=1.0_rsh
     DO i=i_mouth,imax
           ! definition de largeur BM (l2) :
           l2(i)=larbm_embouch*exp(-xfond(i)/lgcvg_largeur)
           l4(i)=l2(i)*l4sl2
           ! Definition de la hauteur d'eau au fond (hl1)
           !hx_g(i,jmilieu)=cotchenal_embouch*exp(-xfond(i)/lgcvg_cotchenal)           
           hx_g(i,jmilieu)=cotchenal_embouch+pentechenal*xfond(i)   
           ! Definition de la section a mi maree
           sectmm(i)=section_embouch*exp(-xfond(i)/lgcvg_section)
           IF((corsect.GE.corsectmax).OR.(corsect.LE.corsectmin))THEN
             deltacorsect=-deltacorsect
           ENDIF
           corsect=corsect+deltacorsect
           sectmm(i)=sectmm(i)*corsect
           ! Definition de l4
           !l4(i)=amax1(l2(i),(sectmm(i)-m*l2(i)-(hx_g(i,jmilieu)-m)*(l1(i)+l2(i)))/m)
           IF((hx_g(i,jmilieu)-m).LE.0.0_rsh)THEN
             l1(i)=l2(i)
           ELSE
             l1(i)=(sectmm(i)-m*(l2(i)+l4(i))-(hx_g(i,jmilieu)-m)*l2(i))/(hx_g(i,jmilieu)-m)
             l1(i)=amin1(l1(i),l2(i))
           ENDIF
           !l2(i)=l2(i)+0.05_rsh
           heaubm=(hx_g(i,jmilieu)-m)*0.5_rsh*(l1(i)+l2(i))/l2(i)
           IF(heaubm.le.0.1_rsh)THEN
             l3(i)=5.0_rsh*l2(i)
           ELSE
             l3(i)=l2(i)*((heaubm+2.0_rsh*m)/heaubm)**2/gamma2
             l3(i)=amin1(5.0_rsh*l2(i),l3(i))
           ENDIF
           l3(i)=amax1(l4(i)+1.0_rsh,l3(i))
     ENDDO    
#endif
 
     ! extension vers le large
     ! ----------------------- 
     ! nouvelle progression exponentielle des sections, maintien de la largeur d'estran et de l2-l1

     l2ml1m=l2(i_mouth)-l1(i_mouth)
     l3ml1m=l3(i_mouth)-l1(i_mouth)
     pentechenal=(cotchenal_embouch-cotfond_off)/(xfond(i_mouth)-xfond(imin))
#if defined key_castestestuary_bathy1 || defined key_castestestuary_bathy2
     DO i=i_mouth-1,imin,-1
        hx_g(i,jmilieu)=cotchenal_embouch+pentechenal*xfond(i)
        hlocchenal=hx_g(i,jmilieu)-m  
        ! Definition de la section a mi maree
        Smm(i)=section_embouch*exp(-varexp_sect_off*xfond(i))
        ! Definition de l1
        l1(i)=(Smm(i)-l2ml1m*(hlocchenal+1.5_rsh*m)-0.5_rsh*l3ml1m*m)/            & 
                            (2.0_rsh*hlocchenal+2.0_rsh*m)
        l2(i)=l1(i)+l2ml1m
        l3(i)=l1(i)+l3ml1m
     ENDDO
#elif defined key_castestestuary_bathy3
     DO i=i_mouth-1,imin,-1
        hx_g(i,jmilieu)=cotchenal_embouch+pentechenal*xfond(i)
        hlocchenal=hx_g(i,jmilieu)-m  
        ! Definition de la section a mi maree
        sectmm(i)=section_embouch*exp(-varexp_sect_off*xfond(i))
        ! Definition de l1
        l1(i)=(sectmm(i)-l2ml1m*(hlocchenal+1.5_rsh*m)-0.5_rsh*l3ml1m*m)/         &
                           (2.0_rsh*hlocchenal+2.0_rsh*m)
        l2(i)=l1(i)+l2ml1m
        l3(i)=l1(i)+l3ml1m
        l4(i)=0.5*(l2(i)+l3(i))
     ENDDO

#endif


    ! Definition des largeurs en y
     DO i=imin+1,imax
       l1y(i)=0.5_rsh*(l1(i)+l1(i-1))
       l2y(i)=0.5_rsh*(l2(i)+l2(i-1))
       l3y(i)=0.5_rsh*(l3(i)+l3(i-1))
#if defined key_castestestuary_bathy3
       l4y(i)=0.5*(l4(i)+l4(i-1))       ! largeur de mi-maree pour bathy 3 
#endif
       hchenal(i)=0.5_rsh*(hx_g(i,jmilieu)+hx_g(i-1,jmilieu))

      ! Definition des hauteurs     
       DO j=jmilieu,jmax
         !en y
         IF(yfond(j).le.l1y(i))THEN
            hy_g(i,j)=hchenal(i)
         ELSE IF (yfond(j).le.l2y(i))THEN
            hy_g(i,j)=hchenal(i)+(m-hchenal(i))/(l2y(i)-l1y(i))*(yfond(j)-l1y(i))
#if defined key_castestestuary_bathy3
         ELSE IF (yfond(j).le.l4y(i))THEN
            hy_g(i,j)=amax1(-5.0_rsh,m-m/(l4y(i)-l2y(i))*(yfond(j)-l2y(i)))
         ELSE 
            denom=(l3y(i)-l4y(i))*(yfond(j)-l4y(i))
            IF (denom.EQ.0.0_rsh) PRINT*,'denom hy',i,l3y(i),l4y(i),yfond(j)
            hy_g(i,j)=amax1(-5.0_rsh,-m/(l3y(i)-l4y(i))*(yfond(j)-l4y(i)))
#else
         ELSE 
            hy_g(i,j)=amax1(-7.0_rsh,m-2.0_rsh*m/(l3y(i)-l2y(i))*(yfond(j)-l2y(i)))
#endif
         ENDIF
         !en x
         IF(yfondx(j).le.l1(i))THEN
            hx_g(i,j)=hx_g(i,jmilieu)
         ELSE IF (yfondx(j).le.l2(i))THEN
            hx_g(i,j)=hx_g(i,jmilieu)+(m-hx_g(i,jmilieu))/(l2(i)-l1(i))*(yfondx(j)-l1(i))
#if defined key_castestestuary_bathy3
         ELSE IF (yfondx(j).le.l4(i))THEN
            hx_g(i,j)=amax1(-5.0_rsh,m-m/(l4(i)-l2(i))*(yfondx(j)-l2(i)))
         ELSE
            denom=(l3(i)-l4(i))*(yfondx(j)-l4(i)) 
            IF (denom.EQ.0.0_rsh) PRINT*,'denom hx', i,j,l3(i),l4(i),yfond(j)
            hx_g(i,j)=amax1(-5.0_rsh,-m/(l3(i)-l4(i))*(yfondx(j)-l4(i)))
#else
         ELSE 
            hx_g(i,j)=amax1(-7.0_rsh,m-2.0_rsh*m/(l3(i)-l2(i))*(yfondx(j)-l2(i)))
#endif
         ENDIF
       ENDDO
     ENDDO

#ifdef key_printdbg          
    ! Visualisation des variables
    DO i=imin,imax
      write(*,*) 'i',i,' l123= ',l1(i),l2(i),l3(i),' chenal=',hx_g(i,jmilieu)
      write(*,*) 'i',i,'l3(i)',l3(i)
      write(*,*) 'i',i,'l2(i)',l2(i)
      write(*,*) 'i',i,'l1(i)',l1(i)
#if defined key_castestestuary_bathy3
      write(*,*) 'i',i,'l4(i)',l4(i)
#endif
    ENDDO
    DO i=imin,imax
      write(*,*)'i',i,'hx_g(i,jmilieu)',hx_g(i,jmilieu)
    ENDDO
    DO i=imin,imax
      write(*,*)'i',i,'Smm(i)',Smm(i)
    ENDDO
    DO i=imin,imax
      write(*,*)'i',i,'hpmshbm2(i)',hpmshbm2(i)
    ENDDO
#endif

    i=imin 
    DO j=jmilieu,jmax
      hy_g(i,j)=hy_g(i+1,j)
      hx_g(i,j)=hx_g(i+1,j)
    ENDDO
    
    !Symetrie
    DO i=imin,imax
      incr=0
      DO j=jmilieu-1,jmin,-1
        incr=incr+1
        hx_g(i,j)=hx_g(i,jmilieu+incr)
        hy_g(i,j)=hy_g(i,jmilieu+incr-1)
      ENDDO
    ENDDO
    
    i=imax-1
    DO j=jmilieu+1,jmax
      hx_g(i,j)=-7.0_rsh 
    ENDDO    
    DO j=jmin,jmilieu-1
       hx_g(i,j)=-7.0_rsh 
    ENDDO    
    DO j=jmin,jmax
      hx_g(imax,j)=-7.0_rsh
      hx_g(imax-1,j)=-7.0_rsh
      hy_g(imax-1,j)=-7.0_rsh
      hy_g(imax,j)=-7.0_rsh
    ENDDO
    ! murs nord et sud :
    do i=imin,imax
      hx_g(i,jmin)=-7.0_rsh
      hy_g(i,jmin)=-7.0_rsh
      hx_g(i,jmin+1)=-7.0_rsh
      hy_g(i,jmin+1)=-7.0_rsh
      hx_g(i,jmax)=-7.0_rsh
      hy_g(i,jmax)=-7.0_rsh
      hx_g(i,jmax-1)=-7.0_rsh
      hy_g(i,jmax-1)=-7.0_rsh
      hy_g(i,jmax-2)=-7.0_rsh
    ENDDO
  
   ! Definition de h0
   ! ----------------
   DO i=imin+1,imax
   DO j=jmin+1,jmax 
     h0_g(i,j)=MAX(hx_g(i-1,j),hx_g(i,j),hy_g(i,j-1),hy_g(i,j))+h0fond   
   END DO
   END DO
   i=imin
   DO j=jmin+1,jmax 
     h0_g(i,j)=MAX(hx_g(i,j),hy_g(i,j-1),hy_g(i,j))+h0fond   
   END DO
   j=jmin
   DO i=imin+1,imax
     h0_g(i,j)=MAX(hx_g(i-1,j),hx_g(i,j),hy_g(i,j))+h0fond   
   END DO
   h0_g(imin,jmin)=MAX(hx_g(imin,jmin),hy_g(imin,jmin))+h0fond 
   

   ! repetition of the end of cas_init_bathy
   ! ---------------------------------------
   hx_t(:,:)=100000.0_rsh
   hy_t(:,:)=100000.0_rsh
   WHERE (hx_g(:,:) /= -valmanq -h0fond) hx_t(:,:)=hx_g(:,:)
   WHERE (hy_g(:,:) /= -valmanq -h0fond) hy_t(:,:)=hy_g(:,:)
   dhxy_g(:,:)=1.0_rsh
   DO j=jmin+1,jmax
      DO i=imin+1,imax
         dhxy_g(i,j)=MAX(hx_g(i,j),hx_g(i-1,j),hy_g(i,j),hy_g(i,j-1))- &
                     MIN(hx_t(i,j),hx_t(i-1,j),hy_t(i,j),hy_t(i,j-1))+epsilon
      END DO
   END DO

   dhx_g(:,:)=0.0_rsh
   dhy_g(:,:)=0.0_rsh
   WHERE (hx_g(:,:) < hrdef) dhx_g(:,:)=dhj
   WHERE (hy_g(:,:) < hrdef) dhy_g(:,:)=dhj

#if defined key_siggen || defined key_gencoord
   hc_sig_g(:,:)=0.0_rsh ; hcx_sig_g(:,:)=0.0_rsh ; hcy_sig_g(:,:)=0.0_rsh
   WHERE (h0_g(:,:)/=-valmanq) hc_sig_g(:,:)=MIN(h0_g(:,:),hc)
   WHERE (hx_g(:,:)/=-valmanq-h0fond) hcx_sig_g(:,:)=MIN(hx_g(:,:),hc)
   WHERE (hy_g(:,:)/=-valmanq-h0fond) hcy_sig_g(:,:)=MIN(hy_g(:,:),hc)
   hc_sig(:,:)  = hc_sig_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hcx_sig(liminm1:limaxp1,ljminm1:ljmaxp1)= hcx_sig_g(liminm1:limaxp1,ljminm1:ljmaxp1)
   hcy_sig(liminm1:limaxp1,ljminm1:ljmaxp1)= hcy_sig_g(liminm1:limaxp1,ljminm1:ljmaxp1)
#endif


   ! local variables
   ! ---------------
   h0(:,:) = h0_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hx(:,:) = hx_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hy(:,:) = hy_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   dhx(:,:) = dhx_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   dhy(:,:) = dhy_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   dhxy(:,:) = dhxy_g(limin:limax,ljmin:ljmax)
   

   ! specificity for morphosedimentary computation 
   ! (l_morphocupl need to be initialize before init_subs because of sed_alloc but after init_bathy) 
   ! --------------------------------------------- 
#ifdef key_sedim_MUSTANG
#if defined key_castestmorpho
   l_morphocoupl=.TRUE.
#else
   l_morphocoupl=.FALSE.
#endif
#endif


   ! 4- Open boundaries 
   ! ==================
   CALL cas_init_obc 
   
      ! over-assignment
      ! geography (where are the open boundaries ?) :
      l_obc_west = .TRUE. 
      ! used scheme to manage obc :
      l_obc_mars = .TRUE.
      ! physics introduced into the simulated area (tide or ogcm solution)
      l_obc_tide = .TRUE.
      !salinity & temperature relaxed to saliref & tetaref during input flow at the W boundary
      l_obc_ogcm=.FALSE.    !.TRUE.
      obc_coefrel=1.00_rlg/(2.0_rlg*3600.0_rlg)
   
   ALLOCATE( obc_tw(kmax,imin  :imin+obc_width+ishift_obc_west,ljmin:ljmax) )
   ALLOCATE( obc_sw(kmax,imin  :imin+obc_width+ishift_obc_west,ljmin:ljmax) )
   obc_sw(:,1,:)=saliref
   obc_tw(:,1,:)=tetaref
   
   tide_ampli_moy=2.0_rsh
   tide_ampli_var=30.0_rsh/70.0_rsh   ! If ME/VE 30.0_rsh/70.0_rsh ELSE 0.0_rsh
   surelev=2.0_rlg                    ! Sea level rise (m/year) 

 
   ! 5- Horizontal viscosity  (none by default)
   ! ===========================================
    CALL cas_init_visc
    
      ! over-assignment
      l_smagor=.TRUE.
      cosmag=1.0_rsh
      vismin=0.001_rsh
      vismax=4.0_rsh
      sponge_nbcell=15
      l_spongw=.TRUE.
      fvisc=5.0_rsh   !0.5_rsh  !1.0    ! ds la zone éponge, nu0=fvisc*0.01_rsh*dy(i,j)**1.15_rsh
      DO j=ljmin,ljmax
        IF (h0_g(ig(j),j) > -valmanq) THEN
          DO i=ig(j),ig(j)+sponge_nbcell
#ifdef key_MPI_2D
            IF ((i.GE.limin).AND.(i.LE.limax)) THEN
#endif
             nu0=fvisc*0.01_rsh*dyval**1.15_rsh*real(sponge_nbcell+2-i)
             vish_phi(:,i,j)=nu0
             vish_xe (:,i,j)=nu0
#ifdef key_MPI_2D
            END IF
#endif
          END DO
        END IF
      ENDDO
 
   ! 6- Vertical Turbulence
   ! ======================
   z0bottom=0.0005_rsh
   z0surface=0.0001_rsh
   nzinit=0.001_rsh
   kzinit=0.00001_rsh
   
   CALL cas_init_turb(kzinit,nzinit,z0surface,z0bottom)
   
   ! over-assignment
      bgdiff=0.000001_rsh
      turb_nbeq=0
      turb_0eq_option=3
      turb_2eq_option=2
      hminfrot=0.1_rsh
      DO j=ljmin,ljmax
        IF (h0_g(ig(j),j) > -valmanq) THEN
          DO i=ig(j),ig(j)+sponge_nbcell
#ifdef key_MPI_2D
           IF ((i.GE.limin).AND.(i.LE.limax)) THEN
#endif
            z0b(i,j)=z0b(i,j)*real(sponge_nbcell+2-i)*2.0_rsh
#ifdef key_MPI_2D
           END IF
#endif
          END DO
        END IF
      ENDDO


   ! 7- Forcing Dynamics : Meteorology (none by default)
   ! ====================================================================
 
   ! no forcing
     CALL cas_init_meteo


   ! 8- Temporal parameters
   ! ======================
   
   !  **************** dates - time steps - start/end of simulaton**********  
     date_start=   '01/01/1900 00:00:00'
     !date_end  =   '01/07/1900 06:00:00'
! test case validation
     date_end  =   '03/01/1900 06:00:00'
     date_start3d=date_start
     date_startagrif=date_start
     dtmin=5.0_rlg
     dtmax=40.0_rlg
     dt   =30.0_rlg
     halfdt=dt*0.5_rlg
#ifdef key_substance
     date_startsub=date_start
     tdebsubs=tool_datosec(date_startsub)
#ifdef key_sedim_MUSTANG
   ! consolidation only if  l_consolid=.true.  or l_diffused or l_bioturb
   date_start_sed=date_start
   tstart_sed=tool_datosec(date_start_sed)
   tstart_sed=MAX(tstart_sed,tdebsubs)
   t_sed=tstart_sed
   dt_sed=500.0_rlg
#endif
#endif
     ! suffixe run
     !------------ 
#if defined key_castestestuary_bathy3
     suffix='run_bathy3'           ! used in result, error_xxx.log, warning_xxx.log and simu_xxx.log files
#elif defined key_castestestuary_bathy2
     suffix='run_bathy2'           ! used in result, error_xxx.log, warning_xxx.log and simu_xxx.log files
#else
     suffix='run_bathy1'           ! used in result, error_xxx.log, warning_xxx.log and simu_xxx.log files
#endif
     l_nc4par = .TRUE.

!  **********************************************************************   
     date_startagrif=date_start3d
     tobs=0.1_rlg*3600.0_rlg                   ! pas de temps de réactualisation de dt  1/2h jusqu'au 29 juillet

   ! IO parametrization
   ! default values of casinit.F90 not changed 
   ! only one file out with all dates, nb_out = 1
   ! if to be changed , set nb_file with a choosed value and call init_io with nb_file in argument
   ! 1DV : all layers in water column are saved
     CALL cas_init_io(date_start,date_end,date_start3d,date_startagrif,suffix,l_nc4par)
     
     ! over-assignment:
             
         iextmin(1)=limin
         iextmax(1)=limax
         jextmin(1)=ljmin
         jextmax(1)=ljmax
         ig(:)=imin
         id(:)=imax
         jb(:)=jmin
         jh(:)=jmax

   ! 9- Saving parameters (results, and possible restart)
   ! ====================================================
 
         ! file name and time step for results :
         pasor(1)=3600.0_rlg                      
         !pasor(1)=200.0_rlg                      
         out_dir='./'
         fileresu(1)='result'//TRIM(suffix)//'.nc'
         fileresu(1)= TRIM(out_dir)//TRIM(fileresu(1))
         
         ! date_startout IF different from date-start:
         ! date_startout='01/01/1900 00:00:00'
         ! tdebsor(1)=tool_datosec(date_startout)
         
         ! date_endout IF different from date-end:
         ! date_endout=date_end
         ! tfinsor(1)=tool_datosec(date_endout)

         l_out_ssh(1)=.TRUE.
         l_out_u3d(1)=.FALSE.
         l_out_ubt(1)=.TRUE.
         l_out_temp(1)=.FALSE.
         l_out_sal(1)=.FALSE.
         l_out_wind(1)=.FALSE.
                  
         l_saverestart_bydate=.TRUE.
         saverestart_step=1344.0_rlg*3600.0_rlg   ! in hours, i.e. IF 1344h: 4 times 14d

        ! restart for this run 
        ! -------------------- 
        l_initfromfile=.FALSE.
        l_init_rtemp=l_initfromfile
        l_init_rsal=l_initfromfile
#if defined key_substance
        l_restart_subs=l_initfromfile
        file_init='../estur-conf36/save_18519000618050035.nc'       !save_12419021201050510.nc'     
#if defined key_sedim_MUSTANG
        l_repsed=l_initfromfile
        filrepsed='../estur-conf36/save_sedim_18519000618050035.nc'   !save_sedim_12419021201050510.nc'   
#endif
#endif         

   ! 10- Initialiation of temperature and salinity  
   ! ============================================

    CALL cas_init_tspar
        ! over-assignment
        kx=10.0_rsh        ! constant coef (m2/s) for tracer diffusion (ref. 0 ?)
        ky=10.0_rsh        ! constant coef (m2/s) for tracer diffusion (ref. 0 ?)
        qmax=8.0_rsh      ! used in quick upwind scheme   (ref. 120000 ?)
        qmaxz=10.0_rsh    ! used in quick upwind scheme   (ref. 1000 ?)
        saliref=35.00_rsh
        tetaref=10.00_rsh
        l_stateeq_lin=.TRUE.
        saliref_lin=saliref
        tetaref_lin=tetaref
   
   ! 11- Substances
   ! ============== 

   !   ATTENTION : set the good nb_var in parameters.F90_cas1DV
   !               INTEGER,PARAMETER :: nb_var_base=????
   !   --------------------------------------------------------
#if defined key_substance
   ! Initialization of the number of variables depending on their type
   ! (gravel,sand,mud,constitutive,non-constitutive,dissolved,fixed,driving or intermediate)
   !----------------------------------------------------------------------------------------
   ! definition of nv.. only if l_filesubs).=.FALSE. 
   ! but it is better here to read a variable.dat :l_filesubs=.TRUE.
     l_filesubs=.FALSE.
     name_filesubs = 'bidon'

   ! ******************* possibility to read a file********************
   IF (l_filesubs) THEN  ! we use "sub_read_var" instead of "cas_init_subs"
     name_filesubs = ''

!*****************************************************************   
   ELSE
     ! definition of nv_grav,nv_sand,nv_mud, nv_cp, nv_ncp, nv_sorb, nv_dis, nb_fix
     ! si different de 0, donner aussi nv_dri et nv_int
     !    nv_adv=nb_var
     !    nv_state=nb_var   nv_dri=0   nv_int=0    nv_tot=nv_state+nv_dri+nv_int
#if defined key_sedim_MUSTANG
     nv_grav=0
     nv_mud=1
     nv_sand=2
#if defined key_nodiss_insed
     nv_dis=0
#else
     nv_dis=1
#endif
     nvpc=3
     nv_ncp=0
     nv_sorb=0
     unit_var(1)='kg.m-3'         
     name_var(1)='sable_grossier'
     name_var(2)='sable_fin'
     unit_var(2)='kg.m-3'         
     name_var(3)='vase'
     unit_var(3)='kg.m-3'         
#if ! defined key_nodiss_insed
     name_var(4)='DISS'
     unit_var(4)='kg.m-3'         
#endif


     nvpc=nv_cp+nv_mud+nv_sand+nv_grav
     nvp=nvpc+nv_ncp+nv_sorb
       
#else
     nv_grav=0
     nv_mud=0
     nv_sand=0
     nv_dis=1
     nvpc=0
     nv_ncp=0
     nv_sorb=0
     unit_var(1)='kg.m-3'         
     name_var(1)='diss'


     nvpc=nv_cp+nv_mud+nv_sand+nv_grav
     nvp=nvpc+nv_ncp+nv_sorb
       
#endif
     ! allocation and initialisation provisoire of cv_wat,f90,t90,  et de tous les tableaux necessaires
     CALL cas_init_subs(date_start)  ! nvp must be defined before
 
!=========================== VARIABLES DEFINITION ======================================
! each variable has to be characterised :
! reminding : 	itypv = 1  for 'GRAV'   (see variables.dat for details...)
!		itypv = 2  for 'SAND'
!		itypv = 3  for 'MUDS'
!		itypv = 3  for 'PART'
!		itypv = 4  for 'NoCP'
!		itypv = 5  for 'SORB'
!		itypv = 6  for 'DISS'
!		itypv = 7  for 'FIXE'
!		itypv = 8  for 'DRIV'
!		itypv = 9  for 'INTE'
#if defined key_sedim_MUSTANG
     isubs=1
     name_var(isubs)='sable_grossier'
     itypv(isubs)=2
     typart(isubs)=1.
     unit_var(isubs)='kg.m-3'         
     ws_free_opt(isubs)=0
     ws_free_min(isubs)=.01_rsh
     ws_free_max(isubs)=.01_rsh
     ws_free_para(1,isubs)=0.0_rsh
     ws_free_para(2,isubs)=0.0_rsh
     ws_free_para(3,isubs)=0.0_rsh
     ws_free_para(4,isubs)=0.0_rsh
     ws_hind_opt(isubs)=0
     ws_hind_para(1,isubs)=.0_rsh
     ws_hind_para(2,isubs)=.0_rsh   
     tocd(isubs)=100.0_rsh
     cini_wat(isubs)=0.0_rsh
     t90(isubs)=0.0_rsh
     sub_flx_atm(isubs)=0.0_rsh
     cv_rain(isubs)=0.0_rsh
     cv_wat(isubs,:,:,:) = cini_wat(isubs)

     isubs=2
     unit_var(isubs)='kg.m-3'
     typart(isubs)=1.
     itypv(isubs)=2
     ws_free_opt(isubs)=0
     ws_free_min(isubs)=.0015_rsh
     ws_free_max(isubs)=.0015_rsh
     ws_free_para(1,isubs)=0.0_rsh
     ws_free_para(2,isubs)=0.0_rsh
     ws_free_para(3,isubs)=0.0_rsh
     ws_free_para(4,isubs)=0.0_rsh
     ws_hind_opt(isubs)=0
     ws_hind_para(1,isubs)=.0_rsh
     ws_hind_para(2,isubs)=.0_rsh
     tocd(isubs)=100.0_rsh
     cini_wat(isubs)=0.0_rsh
     t90(isubs)=0.0_rsh
     sub_flx_atm(isubs)=0.0_rsh
     cv_rain(isubs)=0.0_rsh
     cv_wat(isubs,:,:,:) = cini_wat(isubs)

     isubs=3
     unit_var(isubs)='kg.m-3'
     typart(isubs)=1.
     itypv(isubs)=3
     ws_free_opt(isubs)=0
     ws_free_min(isubs)=.0005_rsh
     ws_free_max(isubs)=.0005_rsh
     ws_free_para(1,isubs)=0.0_rsh
     ws_free_para(2,isubs)=0.0_rsh
     ws_free_para(3,isubs)=0.0_rsh
     ws_free_para(4,isubs)=0.0_rsh
     ws_hind_opt(isubs)=0
     ws_hind_para(1,isubs)=.0_rsh
     ws_hind_para(2,isubs)=.0_rsh
     tocd(isubs)=1.0_rsh
     cini_wat(isubs)=0.0_rsh
     t90(isubs)=0.0_rsh
     sub_flx_atm(isubs)=0.0_rsh
     cv_rain(isubs)=0.0_rsh
     cv_wat(isubs,:,:,:) = cini_wat(isubs)

#if ! defined key_nodiss_insed
     isubs=4
     unit_var(isubs)='kg.m-3'
     typart(isubs)=0.
     itypv(isubs)=6
     cini_wat(isubs)=10.0_rsh
     t90(isubs)=0.0_rsh
     sub_flx_atm(isubs)=0.0_rsh
     cv_rain(isubs)=0.0_rsh
     cv_wat(isubs,:,:,:) = cini_wat(isubs)
     obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
     l_subs2D(isubs)=.FALSE.
#endif

     ! condition a la limite ouverte W (avec relaxation eventuelle) :
     !****************************************************************************
     obc_cv_wat_w(1,:,:,:)=0.0_rsh
     obc_cv_wat_w(2,:,:,:)=0.0_rsh
     obc_cv_wat_w(3,:,:,:)=0.0_rsh
#else
     isubs=1
     name_var(isubs)='Diss'
     itypv(isubs)=6
     typart(isubs)=0.
     unit_var(isubs)='kg.m-3'         
     cini_wat(isubs)=10.0_rsh
     t90(isubs)=0.0_rsh
     sub_flx_atm(isubs)=0.0_rsh
     cv_rain(isubs)=0.0_rsh
     cv_wat(isubs,:,:,:) = cini_wat(isubs)
     obc_cv_wat_w(1,:,:,:)=cini_wat(isubs)
     l_subs2D(isubs)=.FALSE.

#endif

!================================ END of VARIABLES DEFINITION =====================

   ENDIF
     ! if there are particulate variables  or if l_filesubs=true
   IF (nvp>0 .OR. l_filesubs) THEN
       ! if l_filesubs=TRUE: cas_init_subs_disspart call sub-read-var, for direct reading of variable.dat 
       ! and initialise all variables   init_tssub, also initalisation of obc_cv_wat_w/_s/_n/_e
       ! if key_biolo : call biolo_param
       ! if key_sedim_MUSTANG : call sed_alloc
       ! if key_flocmod : call flocmod_alloc
       ! if l_filesubs=FALSE : cas_init_subs_disspart replaces sub_read_var
       CALL cas_init_subs_disspart(l_filesubs)
     
          ! over-assignment
          l_subs2D(1)=.true.
          l_subs2D(2)=.true.
          l_subs2D(3)=.FALSE.
#ifdef key_sand2D
          l_outsandrouse(1)=.true.
          l_outsandrouse(2)=.true.
          l_outsandrouse(3)=.FALSE.
#endif
   
   ENDIF

   ! ******************************** budget computation *************
     l_bilan=.TRUE.
     call cas_init_subs_budget(date_start)
   ! you can modify values of :
     dtout_budget=2.0_rlg  ! in hours
   ! 
   ! *****************************************************************   

   
#endif
! end substances

   ! 12- Initialization of the pronostic variables (ssh,u,v,uz,vz,wz,temp,sal and subsequent parameters)
   ! ===================================================================================================
   ! default values : ssh=u=v=uz=vz=wz=0., temp=tetaref, sal=saliref
   CALL cas_init_progvar
       ! over-assignment:
       DO j=ljminm1,ljmaxp2
       DO i=liminm1,limaxp2
          ssh(i,j)=MAX(0.0_rsh,-h0(i,j))
          IF (h0(i,j) > -valmanq) THEN
            DO k=1,kmax
              temp(k,i,j)= tetaref 
            END DO
          ELSE
            temp(:,i,j) = valmanq
          END IF
       END DO
       END DO
       sal(:,:,:)=saliref
       
   ! Miscellaneous  (ig,id...)
   ! Domain of calculation
   ! Percent of wet surface on each cell
   CALL cas_init_misc

#ifdef key_sedim_MUSTANG

   ! 13- Initialization of sedimentary parameters
   ! ============================================
   
   !******************* before calling cas_init_sedim ********************************
        l_consolid = .FALSE.
        hseduni = 4.0_rsh                  ! initial sediment height  
        cseduni = 1508.0_rsh               ! initial sediment concentration (ignored if l_csedinit_varz=.TRUE.)
        cini_sed(1)=0.1_rsh
        cini_sed(2)=0.9_rsh
        cini_sed(3)=0.0_rsh
#if ! defined key_nodiss_insed
        cini_sed(4)=10.0_rsh
#endif
        l_two_layers=.TRUE.
        ksmiuni = 1                        !  default value of ksmi
        ! if l_dzsinit_varz=TRUE, give  dzs of each initial sediment layer k (m) for i=j=2 [dzs(k,2,2)] 
        ! if l_hsedinit_cste=FALSE give jsed and ksma variable with i,j
        ! if l_csedinit_varz=TRUE, give cv_sed of each variable iv in each layer k for i=j=2 [cv_sed(iv,k,2,2)]     
        ! (including temperature (iv=-1) and salinity (iv=0)
        ! ATTENTION cv_sed in mass/m3 of sediment 
        !           cv_sed of dissolved variables in mass/m3 of pore water!!!       couche de fond sableuse

        IF (l_two_layers) THEN
           ksmauni = 2                     ! initial number of sediment layers
           l_csedinit_varz=.TRUE.           ! true if csed varies within the sediment 
           l_hsedinit_cste=.FALSE.           ! true if hsed does not vary i  spacee (and ksma also)
           l_casinit_hsed=.FALSE.            ! .true. if we want to adjust the sediment thickness in order to be coherent with sediment parameters
                                              ! possible only if l_csedinit_varz=.FALSE. and l_hsedinit_cste=.TRUE.
           l_dzsinit_varz=.TRUE.            ! true if dzs (layer thickness) varies within the sediment

           DO j=ljmin,ljmax
            DO i=limin,limax
            
    !!!     k=ksmi=1  couche de fond sableuse
    !!!     ---------------------------------
              k=ksmiuni
              DO iv=1,isand2
                ! constitutive variables : cini_sed in % of cseduni
                cv_sed(iv,k,i,j)=cseduni*cini_sed(iv)
              ENDDO
              cv_sed(3,k,i,j)=0.0_rsh
              cv_sed(-1,k,i,j)=tini_sed
              cv_sed( 0,k,i,j)=sini_sed
#if ! defined key_nodiss_insed
              cv_sed(4,k,i,j)=cini_sed(4)
#endif
              dzs(k,i,j)=0.5_rsh
              IF(h0_g(i,j).ge.h0_g(i,45) .OR. h0_g(i,j) == -7.0_rsh+h0fond)THEN
                hsed(i,j)=0.5_rsh
                ksma(i,j)=1   !une seule couche de sable au fond du chenal
              ELSE
                hsed(i,j)=h0_g(i,45)-h0_g(i,j)+0.5_rsh       !PLH avant 25 aout : min(hseduni,h0_g(i,45)-h0_g(i,j))
                ksma(i,j)=ksmauni
              ENDIF
           
    !!!     k=2 couche de surface vaseuse (sauf dans chenal)
    !!!     ------------------------------------------------
              k=ksma(i,j)
              IF(k.GE.2)THEN
                dzs(ksma(i,j),i,j)=hsed(i,j)-dzs(1,i,j)
                cv_sed(1,k,i,j)=0.0_rsh
                cv_sed(2,k,i,j)=0.0_rsh
                cv_sed(3,k,i,j)=400.0_rsh
                cv_sed(-1,k,i,j)=tini_sed
                cv_sed( 0,k,i,j)=sini_sed
#if ! defined key_nodiss_insed
                cv_sed(4,k,i,j)=cini_sed(4)
#endif
              ENDIF
              
            ENDDO
           ENDDO
           
        ELSE 
           !  une seule couche melangee mais hsed variable
          !!! ----- ---------------------------------------
           ksmauni = 1                     ! initial number of sediment layers
           cseduni=305.8*2.-0_rsh
           cini_sed(1)=0.5_rsh
           cini_sed(2)=cini_sed(1)
           cini_sed(3)=0.0_rsh
           l_csedinit_varz=.FALSE.           ! true if csed varies within the sediment 
           l_hsedinit_cste=.FALSE.           ! true if hsed does not vary in horizontal space (and ksma)
           l_casinit_hsed=.FALSE.            ! .true. if we want to adjust the sediment thickness in order to be coherent with sediment parameters
                                              ! possible only if l_csedinit_varz=.FALSE. and l_hsedinit_cste=.TRUE.
           l_dzsinit_varz=.FALSE.            ! true if dzs (layer thickness) varies within the sediment
           ! if l_casinit_hsed=FALSE, give hsed 
           DO i=limin,limax
           DO j=ljmin,ljmax
              ksma(i,j)=ksmauni
              IF (h0_g(i,j).ge.h0_g(i,45))then
                 hsed(i,j)=1.0_rsh
              ELSE
                 hsed(i,j)=h0_g(i,45)-h0_g(i,j)+1.0_rsh
              ENDIF
           ENDDO
           ENDDO
        ENDIF
        dzsmaxinit = 100.0_rsh             ! "theoretical" maximum layer thickness
        sinit_sed=0.0_rsh                  ! initial uniform salinity in sediment
        tinit_sed=10.0_rsh                 ! initial uniform temperature in sediment
        z0sedinit=0.0001_rsh              ! bottom roughness for skin shear stress computation

         ! In case l_filesubs=.FALSE. define ros and diam of each particle 
        IF(.NOT.l_filesubs) THEN
          ! first variable : sable_grossier
          ros(1)=2600.0_rsh
          diam_sed(1)=0.002_rsh
          ros(2:nvp)=2600.0_rsh
          diam_sed(2:nvp)=0.0002_rsh
          ros(3)=2600.0_rsh
          diam_sed(3)=0.00005_rsh          
        ENDIF


       CALL cas_init_sedim(l_filesubs,l_csedinit_varz=l_csedinit_varz,l_dzsinit_varz=l_dzsinit_varz, &
                       sinit_sed=sinit_sed,tinit_sed=tinit_sed,z0sedinit=z0sedinit,dzsmaxinit=dzsmaxinit, &
                       l_casinit_hsed=l_casinit_hsed,l_hsedinit_cste=l_hsedinit_cste)

       ! *************parameters to be checked (see parasedim.txt) ********
         dzsmin = 1.0e-5                  ! minimum layer thickness
         cfreshmud = 400.0_rsh            ! concentration of fresh muddy deposits
#if ! defined key_nodiss_insed
         l_diffused=.TRUE.               ! option for diffusion in sediment
#else
         l_diffused=.FALSE.               ! option for diffusion in sediment
#endif    
         l_fricwave=.FALSE.               ! option wave effect in shear stress computation (skin only ???)
         fricwav=0.06_rsh                 ! wave friction coef. (Jonson)
         ! required if l_diffused=.TRUE. : diffusion of dissolved variables in sediment 
!  choice_flxdiss_diffsed : choice for expression of dissolved fluxes at sediment-water interface
!             = 1 : Fick law : gradient between Cv_wat at dz(1)/2
!             = 2 : Fick law : gradient between Cv_wat at distance epdifi
!             = 3 : empirical formula proposed by Boudreau and Shaw-Hanratty, depending on shear velocity
          choice_flxdiss_diffsed=3 ! formulation Boudreau/Shaw-Hanratty
          xdifsi1 = 1.e-6
          xdifs1= 1.e-8
          epdifi = 0.01_rsh  ! not required if choice_flxdiss_diffsed=3

         IF(l_diffused .or. l_consolid .OR. l_bioturb)l_dyn_insed = .TRUE.

         ! initialization of dissolved variables, temperature and salinity in sediment
         !  IF(l_initsed_vardiss=.TRUE.) ==> initialization with initiales concentrations in water at bottom
         !  otherwise by default :  concentrations have already been initialized in cas_init_sedim 
         !  with sinit_sed, tinit_sed and cini_sed
         l_initsed_vardiss=.TRUE.
         IF(l_initsed_vardiss) THEN
            DO k=ksmiuni,ksdmax
               cv_sed(-1,k,:,:)=temp(1,:,:)
               cv_sed(0,k,:,:)=sal(1,:,:)
               DO iv=nvp+1,nv_adv
                 cv_sed(iv,k,:,:)=cv_wat(iv,1,:,:)
               ENDDO
            ENDDO
         ENDIF

         csedmin=30.0_rsh
         csed_mud_ini=400.0_rsh   !mud concentration in inital sediment (normally = cfreshmud if no consolidation)
         xperm1=4.0e-12        ! parameter for permeability : best value for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)
         xperm2=-6.0_rsh       ! parameter for permeability : best value for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)
         xsigma1=6.0e5         ! parameter for effective stress : best value for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)
         xsigma2=6.0_rsh       ! parameter for effective stress : best value for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)
       !!!!!!!!!!!!!!!!!!! dtcal_sed=10.0_rsh
         activlayer=0.02_rsh    ! active layer thickness (m)
         frmudcr2=0.7_rsh      ! threshold above which muddy behaviour
         coef_frmudcr1=1000._rsh       ! coef. for threshold below which sandy behaviour
         x1toce_mud=0.7_rsh            ! coef. for shear strength computation
         x2toce_mud=0._rsh             ! coef. for shear strength computation
         E0_sand_option=0      ! choice of formulation for E0_sand evaluation :
                         !  E0_sand_option= 0 E0_sand= E0_sand read in this namelist 
                         !  E0_sand_option= 1 E0_sand evaluated with Van Rijn (1984) formulation 
                         !  E0_sand_option= 2 E0_sand evaluated with erodimetry formulation (min(0.27,1000*d50-0.01)*(tenfo-toce)**n_eros_sand)
         E0_sand_para=0.0002_rsh            
         n_eros_sand=1.0_rsh           ! power of excess adim. shear stress in sand erosion law
         E0_mud=0.0002_rsh
         n_eros_mud=1._rsh             ! power of excess adim. shear stress in mud erosion law
         ero_option=3                  ! choice of erosion formulation for sand-mud mixing 
           !  ****these formulations are debatable and must be considered carefully by the user.
           !  **** Other laws are possible and have to be programmed. *****
           !  **** Ces formulations sont discutables et doivent etre considerees avec soin par l utlisateur. ***
           !  **** Bien d autres lois sont possibles et peuvent etre programmees *****
                           !     ero_option= 0 : pure mud behavior (for all particles and whatever the mixture)
                           !     ero_option= 1 : linear interpolation between sand and mud behavior, depend on proportions of the mixture
                           !     ero_option= 2 : formulation derived from that of J. Vareilles (2013)
                           !     ero_option= 3 : formulations proposed by B. Mengual (2015) with exponential coefficients depend on proportions of the mixture
         xexp_ero= 2                   ! used only if ero_option=3 : adjustment on exponential variation  (more brutal when xexp_ero high)
         corfluer1=0.0_rsh             ! correction of E0 for very low concentrations (default : 0.0)
         corfluer2=350._rsh            ! concentration defining "low concentrations" in previous correction
         cmudcr=2000._rsh               ! relative mud concentration in surficial layer below which fresh mud deposits can be mixed
         cvolmaxsort=0.58_rsh          ! maximum volumic concentration for sorted sediment (reference value= 0.58)
         cvolmaxmel=0.67_rsh           ! maximum volumic concentration for mixed sediment  (reference value= 0.67)
          xsigma1sg=xsigma1/9.81_rsh
!*****************************************************************  
          fws2=fricwav*0.5_rsh  
!  if key_flocmod : caracterisation of flocs (see cas_init_sedim) 
   ! specificity for morphosedimentary computation
   ! --------------------------------------------- 
 

   ! specificity for morphosedimentary computation
   ! --------------------------------------------- 
#if defined key_castestmorpho 
          !l_morphocoupl=.TRUE.  (deja specifie au debut)
          MF=200.0_rsh                         ! morphological factor (not relevant for 1DV computation)
          !slopecr=.01_rsh (n est pas utilise)
          slopefac=10.0_rsh
          coef_erolat=0.002_rsh
          l_erolat_wet_cell=.FALSE.
          l_bathy_smoothing=.FALSE.
          htncrit_eros=0.05_rsh                 ! water height below which no erosion flux is allowed
#else
          !l_morphocoupl=.FALSE. (deja specifie au debut)
          MF=1.0_rsh          ! si impose, pourrait eviter test dans effdep...
          !slopecr=1.0_rsh     ! valeur par defaut pour ne pas tenir compte effet de pente dans effdep (n est pas utilise)
          slopefac=0.0_rsh    ! valeur par defaut pour ne pas tenir compte effet de pente dans effdep
          coef_erolat=0.0_rsh
          l_erolat_wet_cell=.FALSE.
          htncrit_eros=0.0_rsh               ! water height below which no erosion flux is allowed
#endif 
   
      !  IF (l_repsed) CALL sedinit_fromfile
        IF (l_morphocoupl.and.(.NOT.l_initfromfile)) THEN
           CALL sed_alloc_morpho
           PRINT*,'!!!!!!!!! CALL sed_morphoinit dans casestuar !!!!!!!!'
           CALL sed_morphoinit
        ENDIF

     !  **********choice output in sediment *******************************
     ! output saleinity and temperture in sediment ?
       l_outsed_saltemp=.FALSE.
     !  sediment output over all layers                     --- option 1, default ---
     !  choice_nivsed_out=1  
     !
     !  sediment output over a fixed number of layers from the surface ----option 2---
     !  choice_nivsed_out=2
     !   nk_nivsed_out=10
     !   ep_nivsed_out(:)=0.0_rsh
     ! 
     !  sediment output : max thickness (in m)                       ---- option 3---
     !  choice_nivsed_out=3
     !   thickness max MUST be > dzsmaxuni=dzsmaxinit
     !   nk_nivsed_out=0 
     !   epmax_nivsed_out=0.030_rsh
     !
     !  sediment output : 5 layers max - thicknesses in m            ---- option 4---
     !   donnera une couche supplmentaire avec resultats integres jusqu au fond
     !  choice_nivsed_out=4
     !   nk_nivsed_out=5
     !   ep_nivsed_out(1)=0.025_rsh
     !   ep_nivsed_out(2)=0.025_rsh
     !   ep_nivsed_out(3)=0.025_rsh
     !   ep_nivsed_out(4)=0.050_rsh
     !   ep_nivsed_out(5)=0.100_rsh
     !
   CALL MUSTANG_compatibility
#endif
!  end of step 13 (sediment parameters)
  
  
   ! 14- adding river and outflows (water and variable input water column)
    !   ==================================================================

#if defined key_castestriver                                          
   nbfl=1
   CALL cas_init_river(nbfl)
        ! over-assignment
        iriv(1)=imax-1
        jriv(1)=jmilieu
        kriv(1)=kmax
        !river_thickness(1)=100.0_rsh
        ! for an estuary :
        river_thickness(1)=-1.0_rsh
        debmoy(1)=400.0_rsh
        deb_ampli_var=0.0_rsh
        concssmriver=0.05_rsh
        tmfle(1)=10.0_rsh
        amtfle(1)=0.0_rsh

#if defined key_substance                                                                                               
#ifdef key_sedim_MUSTANG
#ifdef key_nodiss_insed
   nbrej=3
#else
   nbrej=4
#endif

#else
   nbrej=1
#endif

   CALL cas_init_outflow
        ! over-assignment
        numriv(:)=1 
        l_rejflux(:)=.FALSE.                                                                                                  
                              
#ifdef key_sedim_MUSTANG
        nucon(1)=1                                                       
        nucon(2)=2  
        nucon(3)=3                                                                                                           
        ! sables, vase                                     
        valrej(:,1)=0.0_rsh                                             
        valrej(:,2)=0.0_rsh 
!        valrej(:,3)=0.5_rsh   
        !valrej(:,3)=concssmriver*debmoy(1)/kmax 
        valrej(:,3)=concssmriver
#ifdef key_nodiss_insed
        nbrej=3                                                          
#else
        nbrej=4                                                          
        nucon(4)=4                                                                                                           
        valrej(:,4)=10.0_rsh 
#endif
#else
        nbrej=1                                                          
        nucon(1)=1                                                       
        valrej(:,1)=10.0_rsh 
#endif

#endif
! fin substance

#endif
! fin castestriver

   ! for specific information on screen
   ! ==================================
#ifdef key_substance
   WRITE(*,*) 'Number of substances :',nv_tot
  ! DO iv=1,nv_tot
  !   cv_wat(iv,:,:,:) = cini_wat(iv)
  !   WRITE(*,*) 'initialisation :',cini_wat(iv)
  ! ENDDO
#if defined key_sedim_MUSTANG
     WRITE(*,*)'Number of substances :',nv_tot
     DO iv=1,nv_tot
       WRITE(*,*) iv,': "',TRIM(name_var(iv))
       WRITE(*,*) iv,': "',TRIM(name_var(iv))
     END DO
#endif
#endif



   PRINT_DBG*, 'exit CAS_castest_estuar routine'
#endif
  END SUBROUTINE cas_estuar
