  SUBROUTINE cas_tidalflat2dv

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_tidalflat2dv  ***
   !&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 option key_sedim_mixsed ) : 1 mud and 1 dissolved variable
   !&E       key_sedim_mudsand is no more operationnal since r2276 with MUSTANG : 1 mud and 1 dissolved variable
   !&E                                                + key_castestsable : 2 different sand 
   !&E                                         (or) key_castestsablovase : 2 sand and 1 mud variable
   !&E                                               (or) key_castestsub : 1 sand and 2 mud variable and 1 part (non constitutive) and 1 dissolved
   !&E                                             (or) key_castestconta : 2 mud variable and 3 part (non constitutive) and 3 dissolved
   !&E
   !&E                      + key_castestriver                           : with a river upstream and outflow 
   !&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
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : rsh,rlg,riolg,valmanq,kmax,imin,imax,lchain,&
                          ljminm1,ljmaxp2,liminm1,limaxp2,liminp1,limax
   USE comvars2d,  ONLY : fisac,gwac,dtmin,dtmax,dt,halfdt,suffix,fvisc,kx,ky,qmax,iscreenlog
   USE comvars3d,  ONLY : nzinit,kzinit,qmaxz,l_stateeq_lin,saliref_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
   USE toolgeom,  ONLY : f_zu
   USE comvarp2d, ONLY : temp,sal
   USE comvars3d, ONLY : tetaref,saliref
   USE comvars2d, ONLY : icon,fu,fv,h0,tobs,out_dir,hminfrot
   USE obccombine,ONLY : l_obc_mars,l_obc_west,l_obc_tide
   USE comvars2d, ONLY : h0_g,hx_g,hy_g,h0fond,hm,dhj,hrdef,vis2d,vish_xe,vish_phi
   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 sedimorpho,  ONLY : sed_morphoinit
#endif
#if defined key_castestriver
   USE comfluxdo,  ONLY : nbfl,iriv,jriv,kriv,debmoy,tmfle,amtfle,l_tempriv,l_flestat,river_thickness                                     
#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)    :: dxval,dyval,x,dzsmaxinit
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst,zr
   REAL(KIND=rsh)    :: hshallow,hdeep 
   REAL(KIND=rlg)    :: tool_datosec
   LOGICAL           :: l_filesubs


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

   PRINT_DBG*, 'enter CAS_TIDALFALT2DV routine'

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

   ! dynamics
   CALL cas_init_param

   ! over-assignement

   ! tracers : among others, saliref=saliref_lin=35.5_rsh, tetaref=tetaref_lin=10.0_rsh
   CALL cas_init_tspar

   ! over-assignement
   kx=1.0_rsh
   ky=1.0_rsh
   qmax=8.0_rsh
   qmaxz=10.0_rsh

   saliref=35.50_rsh
   l_stateeq_lin=.FALSE.
   saliref_lin=saliref

   ! 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

#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

   ! tidal flat bottom, plus walls
   ! S.Petton : hshallow=-1.0_rsh
   ! S.Petton : hdeep=3.0_rsh

   h0fond=0.01_rsh
   hshallow=-4.0_rsh
   hdeep=16.0_rsh

   hm=2.00_rsh
   hrdef=20.0
   dhj=0.1

   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

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

   ! over_assignement
   ! if using morphodynamic 
#if defined key_castestmorpho && defined key_sedim_MUSTANG
   l_morphocoupl=.TRUE.
#endif

   ! Grid parametrization
   ! --------------------

   ! assignement : input variables for cas_init_grid routine
   dxval=500.0_rsh
   dyval=500.0_rsh
   fisac=0.0_riolg
   gwac=0.0_riolg
   CALL cas_init_grid(dxval,dyval)

   ! over-assignement
   fu(:,:)=0.0_rsh
   fv(:,:)=0.0_rsh


   ! IO parametrization
   ! ------------------

   ! assignement : input variables for cas_init_io routine
   date_start=   '01/01/1900 00:00:00'
#if defined key_castestmorpho
   date_end  =   '05/06/1900 06:00:00'
#else
   date_end  =   '05/01/1900 06:00:00'
#endif
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'
   dtmin=30.0_rlg
   dtmax=50.0_rlg
   dt   =30.0_rlg
   halfdt=dt*0.5_rlg
   suffix='tidalflat2dv'

#ifdef key_substance
   date_startsub='01/01/1900 01:00:00'
   tdebsubs=tool_datosec(date_startsub)
#ifdef key_sedim_MUSTANG
   ! consolidation only if  l_consolid=.true.  or l_diffused or l_bioturb
   date_start_sed='30/01/2000 00:00:00'
   tstart_sed=tool_datosec(date_start_sed)
   tstart_sed=MAX(tstart_sed,tdebsubs)
   t_sed=tstart_sed
   dt_sed=500.0_rlg
#endif
#endif

   CALL cas_init_io(date_start,date_end,date_start3d,date_startagrif,suffix)

   ! over-assignement
   iextmin(1)=imin
   iextmax(1)=imax
   jextmin(1)=2
   jextmax(1)=2
#if defined key_castestmorpho
   pasor(1)=10800.0_rlg 
#else
   pasor(1)=600.0_rlg
#endif
   tobs=2.0_rlg*3600.0_rlg
   out_dir='./'

! nom du fichier resultat :
! =========================================
   fileresu(1)='results.nc'
! =========================================
   fileresu(1)= TRIM(out_dir)//TRIM(fileresu(1))
   date_startout='01/01/1900 00:00:00'
   date_endout=date_end
   tdebsor(1)=tool_datosec(date_startout)
   tfinsor(1)=tool_datosec(date_endout)
   !pasor(1)=3.0_rsh*86400.0_rlg
   l_diag=.FALSE.

   ! Open boundaries
   ! ---------------

   ! no open boundary
   CALL cas_init_obc

   ! over-assignement
   ! geometry (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.

   ! Viscosity  (none by default)
   !-----------

   ! none
   CALL cas_init_visc

   ! over-assignement
   fvisc=0.1_rsh
   vis2d(:,:) = fvisc*0.01_rsh*dyval**1.15_rsh
   vish_xe(:,:,:) = fvisc*0.01_rsh*dyval**1.15_rsh
   vish_phi(:,:,:)= fvisc*0.01_rsh*dyval**1.15_rsh

   ! Turbulence
   ! ----------

   z0bottom=0.0035
   z0surface=1.0_rsh
   nzinit=0.0001_rsh
   kzinit=0.0001_rsh
! ancien tidalflat
   z0bottom=0.00001
   z0surface=0.0001_rsh
   nzinit=0.001_rsh
   kzinit=0.00001_rsh
   
   CALL cas_init_turb(kzinit,nzinit,z0surface,z0bottom)

   ! over-assignement
   bgdiff=0.000001_rsh
   turb_nbeq=2
   turb_0eq_option=1
   turb_2eq_option=2
   !OptionSF=5
   hminfrot=0.001_rsh


   ! Meteorology (none by default)
   ! -----------

   ! no forcing
   CALL cas_init_meteo

   ! Initialization of the pronostic variables
   ! -----------------------------------------

   ! set to 0,tetaref,saliref
   CALL cas_init_progvar

   ! over-assignement
!!!!!!!     N0 = 10e-2 s-1    !!!!!!!
!! EOS lineaire avec -g.dboz/dz = -g * alpha * T | alpha=0.0016)
   DO j=ljminm1,ljmaxp2
     DO i=liminm1,limaxp2
       IF (h0(i,j) > -valmanq) THEN
         DO k=1,kmax
           zr=-f_zu(h0(i,j),0.0_rsh,k,i,j)
           temp(k,i,j)= tetaref - 0.001_rsh/(9.81_rsh*0.0016)*zr
         END DO
       ELSE
         temp(:,i,j) = valmanq
       END IF
     END DO
   END DO

   DO i=liminp1,limax
     x=(i-1)*dxval
     IF (x < 30000.0_rsh ) THEN
       sal(:,i,:)=saliref
     ELSE IF (x > 80000.0_rsh ) THEN
       sal(:,i,:)= 0.0_rsh
     ELSE
       sal(:,i,:)=saliref*(1.0_rsh-(x-30000.0_rsh)/50000.0_rsh)
     END IF
   END DO


   ! Miscellaneous
   ! -------------

   ! Domain of calculation
   ! Percent of wet surface on each cell
   CALL cas_init_misc


   ! main initialization relative to the substances
   ! ----------------------------------------------

   ! assignement : input variables for cas_init_subs routine



   ! Substances
   !--------------
#if defined key_substance
   unit_var(1)='kg.m-3'
#ifdef key_castestsub
   nv_grav=0
   nv_sand=1
   nv_mud=2
   nv_cp=0
   nv_ncp=0
   nv_sorb=1
   nv_dis=1
   nv_fix=0
   nv_dri=0
   nv_int=0
#elif key_castestconta
   nv_grav=0
   nv_sand=0
   nv_mud=2          ! MES_H, MES_L
   nv_cp=0
   nv_ncp=0
   nv_sorb=3          ! Cpcin_H, Cpeq_H, Cpeq_L
   nv_dis=3          ! Csubd,Cd,CdL
   nv_fix=0
   nv_dri=0 
   nv_int=0 
#elif key_castestsable
   nv_dis=0
   nvpc=2
   nv_grav=0
   nv_cp=0
   nv_ncp=0
   nv_sorb=0
   name_var(1)='sable_grossier'
   iv=2
   name_var(iv)='sable_fin'
   unit_var(iv)='kg.m-3'         
   nv_mud=0
   nv_sand=2
#elif key_castestsablovase 
   nv_mud=1
   nv_sand=2
   nv_dis=0
   nvpc=3
   nv_ncp=0
   nv_sorb=0
   name_var(1)='sable_grossier'
   iv=2
   name_var(iv)='sable_fin'
   unit_var(iv)='kg.m-3'
   iv=3         
   name_var(iv)='vase'
   unit_var(iv)='kg.m-3'         
#else
   name_var(1)='vase'
   nv_dis=1
   nvpc=1
   nv_grav=0
   nv_cp=0
   nv_ncp=0
   nv_sorb=0
   name_var(nb_var)='subs_DISS'
   nv_mud=1
   nv_sand=0
#endif

   nvpc=nv_cp+nv_mud+nv_sand+nv_grav
   nvp=nvpc+nv_ncp+nv_sorb
   name_filesubs = 'bidon'
   l_filesubs=.FALSE.

   IF(.NOT.l_filesubs) THEN
   ! allocation et initialisation des tableaux cv_wat,f90,t90...

     CALL cas_init_subs(date_start,date_startsub)   ! nvp must be defined before

   ELSE
     CALL cas_init_subs_budget(date_start)
   ENDIF

   ! over-assignement
   l_bilan=.TRUE.
   dtout_budget=1.00_rlg
 

   ! assignement : input variables for cas_init_subs_disspart routine

!===========================DEFINITION DES VARIABLES======================================
! first variable
!-----------------
#if defined key_castestsable || defined key_castestsablovase
   itypv(1)=2
#else
   name_var(1)='vase'
   itypv(1)=3
#endif
   typart(1)=1.
   ws_free_opt(1)=0
   ws_free_min(1)=.0005_rsh
   ws_free_max(1)=.0005_rsh
   ws_free_para(1,1)=0.0_rsh
   ws_free_para(2,1)=0.0_rsh
   ws_free_para(3,1)=0.0_rsh
   ws_free_para(4,1)=0.0_rsh
   ws_hind_opt(1)=0
   ws_hind_para(1,1)=.0_rsh
   ws_hind_para(2,1)=.0_rsh   

   tocd(1)=100.0_rsh
#ifdef key_sedim_MUSTANG
   cini_wat(1)=0.0_rsh
! cas test florence B2
!   cini_wat(1)=0.1_rsh
#else
   cini_wat(1)=0.05_rsh
#endif    
   t90(1)=0.0_rsh
   sub_flx_atm(1)=0.0_rsh
   cv_rain(1)=0.0_rsh
   cv_wat(1,:,80:150,:) = cini_wat(1)
   obc_cv_wat_w(1,:,:,:)=cini_wat(1)

!  second variable 
!-----------------
   isubs=2  
   unit_var(isubs)='kg.m-3'
#if defined key_castestsable || defined key_castestsablovase
!  variable sand
   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
   cv_wat(isubs,:,:,:) = cini_wat(2)
#else
!  variable dissoute
   name_var(isubs)='DISS'
   typart(isubs)=0.
   itypv(isubs)=6
   cini_wat(isubs)=100.0_rsh
   cv_wat(isubs,:,80:150,:) = cini_wat(2)
#endif
   t90(isubs)=0.0_rsh
   sub_flx_atm(isubs)=0.0_rsh
   cv_rain(isubs)=0.0_rsh
   DO isubs=1,2
       obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs) 
   ENDDO

!  third variable 
!-----------------  

#if defined key_castestsablovase
   unit_var(3)='kg.m-3'
!  variable mud
   name_var(3)='vase'
   typart(3)=1.
   itypv(3)=3
   ws_free_opt(3)=0
   ws_free_min(3)=.0005_rsh
   ws_free_max(3)=.0005_rsh
   ws_free_para(1,3)=0.0_rsh
   ws_free_para(2,3)=0.0_rsh
   ws_free_para(3,3)=0.0_rsh
   ws_free_para(4,3)=0.0_rsh
   ws_hind_opt(3)=0
   ws_hind_para(1,3)=.0_rsh
   ws_hind_para(2,3)=.0_rsh 
   tocd(3)=1.0_rsh
   cini_wat(3)=0.02_rsh
   cv_wat(3,:,:,:) = cini_wat(3)
   t90(3)=0.0_rsh
   sub_flx_atm(3)=0.0_rsh
   cv_rain(3)=0.0_rsh
   obc_cv_wat_w(3,:,:,:)=cini_wat(3) 
#endif

! other tests
!-------------    
#if defined key_castestsub
   name_var(1)='sable'
   typart(1)=1.
   itypv(1)=2
   ws_free_opt(1)=0
   ws_free_min(1)=.0005_rsh
   ws_free_max(1)=.0005_rsh
   ws_free_para(1,1)=0.0_rsh
   ws_free_para(2,1)=0.0_rsh
   ws_free_para(3,1)=0.0_rsh
   ws_free_para(4,1)=0.0_rsh
   ws_hind_opt(1)=0
   ws_hind_para(1,1)=.0_rsh
   ws_hind_para(2,1)=.0_rsh 
   tocd(1)=100.0_rsh
   cini_wat(1)=0.0_rsh
#endif
#if defined key_castestsub || defined key_castestconta
#if defined key_castestsub
   isubs=2
#elif key_castestconta
   isubs=1
#endif
   name_var(isubs)='vase_H'
   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.02_rsh     
   t90(isubs)=0.0_rsh
   sub_flx_atm(isubs)=0.0_rsh
   cv_rain(isubs)=0.0_rsh
   cv_wat(isubs,:,80:150,:) = cini_wat(isubs)
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)

#if defined key_castestsub
   isubs=3
#elif key_castestconta
   isubs=2
#endif
   name_var(isubs)='vase_L'
   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.05_rsh     
   t90(isubs)=0.0_rsh
   sub_flx_atm(isubs)=0.0_rsh
   cv_rain(isubs)=0.0_rsh
   cv_wat(isubs,:,80:150,:) = cini_wat(isubs)
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
#endif

#ifdef key_castestsub
! one non constitutive particulate substance
!  sorbed on the first constitutive particulate substance
   isubs=4
   name_var(isubs)='part_H'
   unit_var(isubs)='g.m-3'
   typart(isubs)=0.
   itypv(isubs)=5
   name_varpc_assoc(isubs)='vase_H'
   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)
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
 
! one dissolved substance
   isubs=5
   name_var(isubs)='diss1'
   unit_var(isubs)='kg.m-3'
   typart(isubs)=0
   itypv(isubs)=6
   cini_wat(isubs)=100._rsh
   t90(isubs)=0._rsh
   !sub_flx_atm(isubs)=0.0000001_rsh
   sub_flx_atm(isubs)=0.0_rsh
   cv_rain(isubs)=0._rsh
   cv_wat(isubs,:,10:100,:) = cini_wat(isubs)
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
#endif

#if defined key_castestconta
! one first non constitutive particulate substance
!  sorbed on the second constitutive particulate substance
   isubs=3
   name_var(isubs)='Cp_Leq'
   unit_var(isubs)='g.m-3'
   typart(isubs)=0
   itypv(isubs)=5
   name_varpc_assoc(isubs)='vase_L'
   tocd(isubs)=1.0_rsh
   cini_wat(isubs)=0.1_rsh*cini_wat(2)   
   t90(isubs)=0.0_rsh
   sub_flx_atm(isubs)=0.0_rsh
   cv_rain(isubs)=0.0_rsh
   cv_wat(isubs,:,80:150,:) = cini_wat(isubs)
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
  
! a second non constitutive particulate substance
!  sorbed on the first constitutive particulate substance
! nulle dans l eau initialement, presente dans sediment sur la deuxieme variable contitutive
   isubs=4
   name_var(isubs)='Cp_Hcin'
   unit_var(isubs)='g.m-3'
   typart(isubs)=0
   itypv(isubs)=5
   name_varpc_assoc(isubs)='vase_H'
   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)
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
 
! a third non constitutive particulate substance
  ! nulle dans l eau initialement, presente dans sediment sur la deuxieme variable contitutive
!  sorbed on the first constitutive particulate substance
   isubs=5
   name_var(isubs)='Cp_Heq'
   unit_var(isubs)='g.m-3'
   typart(isubs)=0
   itypv(isubs)=5
   name_varpc_assoc(isubs)='vase_H'
   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)
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
 
! one first dissolved substance
   isubs=6
   name_var(isubs)='Csubd'
   unit_var(isubs)='kg.m-3'
   typart(isubs)=0
   itypv(isubs)=6
   cini_wat(isubs)=0.005_rsh
   t90(isubs)=0.0_rsh
   sub_flx_atm(isubs)=0.0_rsh
   cv_rain(isubs)=0.0_rsh
   DO i=imin,imax
        cv_wat(isubs,:,i,:) = cini_wat(isubs)*(0.1+0.9*(imax-i)/imax)
   END DO
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
          
! one second dissolved substance
   isubs=7
   name_var(isubs)='Cd'
   unit_var(isubs)='kg.m-3'
   typart(isubs)=0
   itypv(isubs)=6
   cini_wat(isubs)=0.005_rsh
   t90(isubs)=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)
   
! one thrid dissolved substance
   isubs=8
   name_var(isubs)='CdL'
   unit_var(isubs)='kg.m-3'
   typart(isubs)=0
   itypv(isubs)=6
   cini_wat(isubs)=0.00_rsh
   t90(isubs)=0.0_rsh
   sub_flx_atm(isubs)=0.0001_rsh
   cv_rain(isubs)=0.0_rsh
   cv_wat(isubs,:,:,:) = cini_wat(isubs)
   obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs)
#endif
!================================ FIN DEFINITION DES VARIABLES =====================

   ! si il existe des variables particulaires 
   IF (nvp>0) THEN

      CALL cas_init_subs_disspart(l_filesubs)

   ! over-assignement
#if defined key_castestsable || defined key_castestsablovase || defined key_castestsub
#ifdef key_sand2D
      l_subs2D(1)=.true.
      l_outsandrouse(1)=.true.
#endif
#endif

      DO j=ljminm1,ljmaxp2
        DO i=liminm1,limaxp2
          IF (h0(i,j) > -valmanq) THEN
    !        DO k=1,kmax
    !          zr=-f_zu(h0(i,j),0.0_rsh,k,i,j)
    !          cv_wat(1,k,i,j)= 0.1 - 0.00001_rsh/(9.81_rsh*0.0016)*zr
    !          cv_wat(2,k,i,j)= 50. - 0.001_rsh/(9.81_rsh*0.0016)*zr
     !        cv_wat(:,k,i,j) = cini_wat(:)
     !       END DO

          ELSE
             cv_wat(:,:,i,j) = valmanq
          END IF
        END DO
       END DO
    
    ! obligation de definir les ros et diam des particules apres l appel a cas_init_subs_disspart dans lasuelle ces variable sont allouees
    ! vase ne se sert pas du diametre mais il ne doit pas etre nul et il se sert de ros
#ifdef key_sedim_MUSTANG
! first variable : sable_grossier
       ros(1)=2600.0_rsh
       diam_sed(1)=0.0002_rsh
       !diam_sed(1)=0.0001_rsh
#if defined key_castestsable || defined key_castestsablovase
       ros(2:nvp)=2600.0_rsh
       diam_sed(2:nvp)=0.0001_rsh
#elif key_castestsub
       ros(2:nvp)=2600.0_rsh
       diam_sed(2:nvp)=0.0001_rsh
#elif key_castestconta
       ros(1:nvp)=2600.0_rsh
       diam_sed(1:nvp)=0.0001_rsh
#endif
#endif

   ENDIF
#endif
       ! fin substance


   ! River and outflows
   ! -------------------

#if defined key_castestriver                                          
   nbfl=1
   CALL cas_init_river(nbfl)

   ! over-assignement

   iriv(1)=200
   jriv(1)=2
   kriv(1)=10
   river_thickness(1)=100.0_rsh
   debmoy(1)=25.0_rsh
   tmfle(1)=10.0_rsh
   amtfle(1)=0.0_rsh 

#if defined key_substance
#ifdef key_castestsub                                            
   nbrej=3                                                                                                      
#elif key_castestconta                                                
   nbrej=4                                                                                                    
#elif key_castestsablovase 
   nbrej=3
#else
   nbrej=2
#endif

   CALL cas_init_outflow

   ! over-assignement
   numriv(:)=1 
   l_rejflux(:)=.false.                                                                                                  
#ifdef key_castestsub                                            
   nucon(1)=2                                                       
   nucon(2)=4                                                       
   nucon(3)=5                                                       
   valrej(:,1)=1.0_rsh
   valrej(:,2)=10.0_rsh                                              
   valrej(:,3)=10.0_rsh                                              
#elif key_castestconta                                                
   nbrej=4                                                          
   nucon(1)=1                                                       
   nucon(2)=4                                                       
   nucon(3)=6                                                       
   nucon(4)=7                                                       
   ! vase, CP_HCIN, CSUBD, CD                                       
   valrej(:,1)=0.02_rsh                                             
   valrej(:,2)=0.003_rsh                                            
   valrej(:,3)=0.005_rsh                                            
   valrej(:,4)=0.0005_rsh                                           
#elif key_castestsablovase                                               
   nbrej=3                                                          
   nucon(1)=1                                                       
   nucon(2)=2  
   nucon(3)=3                                                                                                           
   ! sables, vase                                     
   valrej(:,1)=0.1_rsh                                             
   valrej(:,2)=0.2_rsh 
   valrej(:,3)=0.5_rsh   
#elif key_castestsable                                               
   nbrej=2                                                          
   nucon(1)=1                                                       
   nucon(2)=2                                                       
   ! sable, dissous                                     
   valrej(:,1)=0.1_rsh                                             
   valrej(:,2)=0.5_rsh                                                                                      
#endif
#endif
         ! fin substance
#endif
         ! fin castestriver

#ifdef key_sedim_MUSTANG

   ! Initialization of sedimentary parameters
   ! ----------------------------------------
   ! pre-assignement
   l_consolid = .FALSE.
   hseduni = 0.1_rsh
   ksmiuni = 1
   ksmauni =2

! test si on veut sortie des couches 
!   ksmauni=5
 
!! test avec florence B2
!    hseduni=0._rsh

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! ATTTENTION !!!
!!!!!!!!!!!!!!!!!!
!  the total concentration is given as an input, and the contribution
!                 of each variable (cini_sed(iv)) is given as a fraction of c_sedtot 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#if defined key_castestsablovase
   ! constitutives in fraction de 1 sable puis vase
   cseduni = 1508.0_rsh
   c_sedtot(:,:,:)=cseduni
   cini_sed(1)=0.4_rsh
   cini_sed(2)=0.4_rsh
   cini_sed(3)=0.2_rsh
#elif key_castestsable
   ! constitutives in fraction de 1 
   cseduni = 1508.0_rsh
   c_sedtot(:,:,:)=cseduni
   cini_sed(1)=0.40_rsh
   cini_sed(2)=0.6_rsh
#elif key_castestsub
   ! constitutives in fraction de 1 sable puis vase
   cseduni = 1508.0_rsh
   c_sedtot(:,:,:)=cseduni
   cini_sed(1)=0.98_rsh
   cini_sed(2)=0.02_rsh
   if(nb_var>2)cini_sed(3:nb_var)=0.0_rsh
#else
   cseduni = 400.0_rsh
   c_sedtot(:,:,:)=cseduni
   ! constitutives in fraction de 1 (vase seule)
   cini_sed(1)=1.0_rsh
   if(nb_var>1)cini_sed(2:nb_var)=0.0_rsh
#endif

#ifdef key_castestsub
   cini_sed(3)=0.0_rsh
   ! non constitutives in masse/kg of associated constitutive substance
   cini_sed(4)=10.0_rsh
   ! dissolved in masse/m3 of pore water 
   cini_sed(5)=0.0_rsh
#elif key_castestconta
   cini_sed(2)=0.0_rsh
   ! non constitutive in masse/kg of associated constitutive substance
   cini_sed(3)=0.0_rsh
   cini_sed(4)=1.0_rsh
   cini_sed(5)=0.01_rsh
   ! dissolved in masse/m3 of pore water 
   cini_sed(6)=0.0_rsh
   cini_sed(7)=0.0_rsh
   cini_sed(8)=0.0_rsh
#endif

! sediment output sur toutes les couches ---option 1---
   nk_nivsed_out=ksdmax
   ep_nivsed_out(:)=0.0_rsh

! sediment output sur un nombre fix de  couches en partant de la surface ----option 2---
!   nk_nivsed_out=3
 !  ep_nivsed_out(:)=0.0_rsh

!  sediment output : epaisseur max (en m) -- option 3--
!   nk_nivsed_out=0
!   epmax_nivsed_out=0.030_rsh

!  sediment output : 5 couches max - epaisseurs en m ----option 4----
!   nk_nivsed_out=4
!   ep_nivsed_out(1)=0.001_rsh
!   ep_nivsed_out(2)=0.005_rsh
!   ep_nivsed_out(3)=0.010_rsh
!   ep_nivsed_out(4)=0.050_rsh


   ! Initialization of sediment layers
   ! ---------------------------------
   dzsmaxinit = 10.0_rsh
   CALL cas_init_sedim(l_filesubs,dzsmaxinit=dzsmaxinit)

   ! over-assignement
   dzsmin = 1.0e-5
   cfreshmud = 400.0_rsh

#if defined key_castestmorpho && defined key_sedim_MUSTANG
! could be different from the default (equal to 0) if morpho using
   htncrit_eros=0.0_rsh
#endif
   IF(l_morphocoupl) CALL sed_morphoinit

! test sortie des couches 
!   dzsmaxuni = 0.02_rsh

#endif
        ! fin sedim
!
!OMPMPI barrier
!OMPMPI master
!   CALL_MPI ex_j_rsh(-1,0,1,limin,limax,ljminm1,ljmax,govlet(limin:limax,ljminm1:ljmax))
!   CALL_MPI ex_i_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,dx(liminm1:limaxp1,ljminm1:ljmaxp1))
!   CALL_MPI ex_j_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,dx(liminm1:limaxp1,ljminm1:ljmaxp1))
!   CALL_MPI ex_i_rsh(0,1,1,limin,limaxp1,ljminm1,ljmax,dx_v(limin:limaxp1,ljminm1:ljmax))
!   CALL_MPI ex_j_rsh(-1,0,1,limin,limaxp1,ljminm1,ljmax,dx_v(limin:limaxp1,ljminm1:ljmax))
!   CALL_MPI ex_i_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,dy(liminm1:limaxp1,ljminm1:ljmaxp1))
!   CALL_MPI ex_j_rsh(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,dy(liminm1:limaxp1,ljminm1:ljmaxp1))
!   CALL_MPI ex_i_rsh(-1,0,1,liminm1,limax,ljmin,ljmaxp1,dy_u(liminm1:limax,ljmin:ljmaxp1))
!   CALL_MPI ex_j_rsh(0,1,1,liminm1,limax,ljmin,ljmaxp1,dy_u(liminm1:limax,ljmin:ljmaxp1))
!   CALL_MPI ex_i_rsh(-1,0,1,liminm1,limax,ljmin,ljmax,dx_u(liminm1:limax,ljmin:ljmax))
!   CALL_MPI ex_i_rsh(-1,0,1,liminm1,limax,ljmin,ljmax,dx_f(liminm1:limax,ljmin:ljmax))
!   CALL_MPI ex_j_rsh(-1,0,1,limin,limax,ljminm1,ljmax,dy_v(limin:limax,ljminm1:ljmax))
!   CALL_MPI ex_j_rsh(-1,0,1,limin,limax,ljminm1,ljmax,dy_f(limin:limax,ljminm1:ljmax))
!OMPMPI end master
!OMPMPI barrier
!OMPMPI flush(govlet,dx,dx_v,dy,dy_u,dx_u,dx_f,dy_v,dy_f)



   PRINT_DBG*, 'exit CAS_TIDALFLAT2DV routine'

  END SUBROUTINE cas_tidalflat2dv
