  SUBROUTINE cas_1DV
#if defined key_1DV

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_1DV  ***
   !&E
   !&E ** Purpose : definition of a test case for 1DV general  routine test
   !&E              For developpers
   !&E
   !&E ** Description :
   !&E                  1DV test case based on:
   !&E                  - 1DV hydrodynamics test case "cas_turbkatophil"
   !&E                  - 2DV sediment test case "cas_tidalflat2DV"
   !&E
   !&E   MUST be USED with : key_1DV, key_castest, key_dyn_adv_linear
   !&E          forcing by :
   !&E              either : pente de la surface libre                                          : key_1DV_force_gradssh 
   !&E              either : courant barotrope oscillant                                        : key_1DV_force_ubarotrope   
   !&E              either : courant barotrope lu dans un fichier                               : key_1DV_force_ufile  
   !&E              either : data of tenfon due to wave in an ASCII file issued from a 3D model : key_1DV_force_towfile
   !&E
   !&E           possibility to have wind (define wind_veloc)   
   !&E           possibility to have a oscillating barotropic courrent or surface slope (to be programmed in forcingKtest.F90)
   !&E           possibility to introduce an outflow (apport de variable)
   !&E           possibility to introduce bioturbation (key_castest_1DVbioturb)) (not operational)
   !&E
   !&E
   !&E   *****PARAMETERS AND DATE THE USER HAS TO CHANGE OR VERIFY ARE BETWEEN STARS LINES***********
   !&E   
   !&E   *****LES PARAMETRES ET LES DONNEES A MODIFIER OU A VERIFIER PAR L UTILISATEUR POUR LE CAS1dv**********
   !&E          **********SONT NOTEES ENTRE DES LIGNES D ETOILES********************
   !&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       !  2014-05  (F. Grasso) new test case for sediment consolidation
   !&E       !  2015-04  (B.Thouvenin & P. Le Hir) creating new test case cas1DV general
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : rsh,rlg,riolg,valmanq,kmax,imin,imax,lchain,&
                          ljminm1,ljmaxp2,liminm1,limaxp2
   USE comvarp3d, ONLY : ect
#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
   USE output,    ONLY : fileresu,tdebsor,tfinsor,nb_out,pasor,   &
                         l_out_ssh,l_out_ubt,l_out_u3d, &
                         l_out_sal,l_out_temp,l_out_wind,l_out_buoy,l_out_wave,&
                         l_out_visc,l_out_turb,l_out_off,l_out_atmsflx,        &
                         iextmin,iextmax,stepk,igextmin,igextmax,          &
                         jextmin,jextmax,kextmin,kextmax,l_out_zaxis_w,   &
                         jgextmin,jgextmax,l_filebydate,file_output
   USE toolgeom,  ONLY : f_zu
   USE comvarp2d, ONLY : temp,sal
   USE comvars3d, ONLY : tetaref,saliref,nzinit,kzinit,qmaxz,l_stateeq_lin,rhoref,  &
                         turb_0eq_option,turb_nbeq,turb_2eq_option,bgdiff,OptionSF, &
                         tetaref_lin,saliref_lin
   USE comvars2d, ONLY : fisac,gwac,dtmin,dtmax,dt,suffix,fvisc,kx,ky,qmax,halfdt, &
                         icon,fu,fv,h0,tobs,hminfrot,iscreenlog,ierrorlog,icon, &
                         h0_g,hx_g,hy_g,h0fond,hm,vis2d,vish_xe,vish_phi,name_out_xe, &
                         botstressmax,hminfrot,hmaxfrot
   USE comdiag,   ONLY : l_diag,file_diag,tdebdiag,tfindiag
   USE comturb,   ONLY : eps
   USE cometdyn,   ONLY : wind_veloc,wind_dir
#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
#endif
   USE cometdyn    , ONLY : gradxe_force_1DV_cste, gradxe_force_1DV_ampli,ubaro_force_1DV_cste,  &
                            ubaro_force_1DV_ampli, gradxe_force_1DV, ubaro_force_1DV,  &
                            l_force1DV_variable,cds,modulation_ampli_tide
#ifdef key_biolo
   USE cometemp,   ONLY : sflx_solarcst,sflx_rad_maxdepth
   USE bioloinit
   USE subreaddat, ONLY : sub_read_vardiag
#endif
#if defined key_1DV_force_towfile
   USE readforcingk1DV, ONLY : file_tenfonw_1DV,cas1dv_read_forcing
#elif defined key_1DV_force_ufile
   USE readforcingk1DV, ONLY : file_ubarotrop_1DV,cas1dv_read_forcing
#endif 

   IMPLICIT NONE

   !! * Local declarations
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   CHARACTER(LEN=19) :: date_startout,date_endout
   INTEGER           :: i,j,k,isubs,iv,ig
   INTEGER           :: ni_wall_e,ni_wall_w,nj_wall_n,nj_wall_s
   INTEGER           :: IERR_MPI
   REAL(KIND=rsh)    :: dxval,dyval,z0bottom,z0surface,depthcst,zr
   REAL(KIND=rlg)    :: tool_datosec
   LOGICAL           :: l_filesubs
#ifdef key_sedim_MUSTANG
   REAL(KIND=rsh)    :: sinit_sed,tinit_sed,z0sedinit,dzsmaxinit
   LOGICAL           :: l_csedinit_varz,l_dzsinit_varz,l_casinit_hsed
#endif

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

   PRINT_DBG*, 'enter CAS_1DV_general routine'


   ! 1- general parametrization (each parameter can be reassigned)
   ! ==========================
   ! define h0fond, hm, hminim, fwetmin, dhj, rhoair, chp, coext
   ! define l_tssubadv_over_smalldepth, l_genpot, l_stability
   CALL cas_init_param
   

   ! 2- Water depth  (Bathymetry)
   ! ===========================
   ! assignement : input variables for cas_init_bathy routine
   ! the domain is seen as a WE channel (flat bottom, plus walls)
   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
   hc= 0.1_rsh
#endif
   hm=0.5_rsh
!  ************************* height of water column (metres)************
#ifdef key_biolo
   depthcst = 100.0_rsh
#else
   depthcst = 30.0_rsh
#endif
!  ********************************************************************* 
   CALL cas_init_bathy(depthcst,nj_wall_n=nj_wall_n,nj_wall_s=nj_wall_s)



   ! 3- Grid parametrization in the water column
   ! ===========================================
   ! latitude - longitude (Coriolis non operational), dx & dy
   fisac=1.0_riolg
   gwac=1.0_riolg
   dxval = 200.0_rsh
   dyval = 200.0_rsh
!  *********************** if generalised sigma : parameters selection**
#if defined key_siggen || defined key_gencoord
   theta_sig= 6.0_rsh ! entre 0 et 20
   b_sig= 0.0_rsh     ! entre 0 et 1
#endif
!  *********************************************************************
   CALL cas_init_grid(dxval,dyval)
   !  no Coriolis effect here in 1DV (computation in one direction only)
   fu(:,:)=0.0_rsh
   fv(:,:)=0.0_rsh


   ! 4- Open boundaries and horizontal viscosity  (none by default)
   ! ===========================================
   CALL cas_init_obc  
   CALL cas_init_visc
 
   
   ! 5- Vertical Turbulence
   ! ======================
!  **** selection of roughness and vertical viscosity and diffusivity ******
   z0bottom=0.0001_rsh
   z0surface  =0.0001_rsh
   nzinit=0.00001_rsh       !initial vertical viscosity   (.001)
   kzinit=0.00001_rsh     !initial vertical diffusivity (.00001)
#ifdef key_biolo
   kzinit=0.001_rsh
   z0surface  =1.0_rsh
#endif
   CALL cas_init_turb(kzinit,nzinit,z0surface,z0bottom)
   ! over-assignment
!  ****************************** selection of turbulence closure **********
   turb_nbeq=2         ! number of equations in the turbulence model(0, 1 or 2)
   turb_0eq_option=1   ! option number for the zero equation model
   turb_2eq_option=2   !  option number for the two-equation model 
                       ! ( 1- k-kl ; 2- k-epsilon; 3- k-omega; 4- generic modele)
   ect(:,:,:)=0.000001_rsh   ! initial turbulent kinetic energy
   eps(:,:,:)=0.000001_rsh   ! initial disspation rate of ect
   bgdiff=0.000001_rsh  ! background vertical viscosity and diffusivity (in m2/s), 
                        ! with values such that nz >= bgdiff and kz >= bgdiff
!   hminfrot=0.5_rsh  ! minimum water column height used to compute bottom and surface shear stresses
!   hmaxfrot=5.0_rsh  ! maximimum thickness of the bottom boundary layer 
                      ! (used to compute bottomshear stresses)
#ifdef key_biolo
   turb_nbeq=2
   turb_0eq_option=0
   turb_2eq_option=2
   ect(:,:,:)=0.0001_rsh
   eps(:,:,:)=0.0001_rsh
   bgdiff=0.0001_rsh
   OptionSF=5
#endif
!***************************************************************************


   ! 6- Forcing Dynamics : Sea Surface Slope or barotropic current & wind
   ! ====================================================================
   ! cste : steady forcing
   ! ampli : mean tidal amplitude
   ! modulation_ampli_tide : modulation of tidal amplitude expressed in variation of tidal coefficient (70+-modul)
!  *************************************************************************
!  define if constant or variable forcing, and give the constant value
!  IF l_forceDV_variable=.TRUE. program the variation in forcingKtest.F90 with terms _cste and _ampli
! If forcing with a courrent readed in a file : choice l_force1DV_variable=.FALSE. and use CPP key :  key_1DV_force_ufile
   l_force1DV_variable=.TRUE.
   IF(l_force1DV_variable) THEN 
#if defined key_1DV_force_gradssh
     ! the sea surface slope is given : its order of magnitude can be
     ! deduced from the anticipated barotropic current U, 
     ! either through the approximation slope = 1.e-5 x U
     ! or slope =  1/g*[K.U/ln(h/ez0)]^2
     gradxe_force_1DV_ampli= 0.00001_rsh   ! (m /m linear)
     gradxe_force_1DV_cste= 0.0_rsh        ! (m /m linear)
     modulation_ampli_tide=30._rsh         ! simulaton a tidal coefficient between 70-30 to 70+30 
#elif key_1DV_force_ubarotrope
     ubaro_force_1DV_ampli=0.8_rsh         ! (m /s)
     ubaro_force_1DV_cste=0.1_rsh          ! (m /s)
     modulation_ampli_tide=30._rsh         ! simulaton a tidal coefficient between 70-30 to 70+30 
#endif 
   ELSE
     gradxe_force_1DV_cste= 0.00001_rsh
     ubaro_force_1DV_cste=0.25_rsh
   ENDIF
   gradxe_force_1DV=gradxe_force_1DV_cste
   ubaro_force_1DV=ubaro_force_1DV_cste
   
   !! verification
#if defined key_1DV_force_ufile && defined key_1DV_force_ubarotrope
        WRITE(ierrorlog,*)
        WRITE(ierrorlog,*)'****************************'
        WRITE(ierrorlog,*)' You have chosen 2 incompatible CPP keys'     
        WRITE(ierrorlog,*)' key_1DV_force_ufile and key_1DV_force_ubarotrope '     
        WRITE(ierrorlog,*)' choose between : '
        WRITE(ierrorlog,*)'key_1DV_force_ubarotrope and l_force1DV_variable=.TRUE. : oscillating U'
        WRITE(ierrorlog,*)'key_1DV_force_ufile and l_force1DV_variable=.FALSE. :  U in a file'
        STOP
#endif

!  ************************************************************************   
   ! Meteorology (none by default)
   ! default values of casinit.F90 not changed :
   ! l_meteo_stat=.TRUE.   l_meteo_hom=.TRUE.
   ! windamp_duration = 0.0_rlg     l_sflx_solarcst=.TRUE.   sflx_solarcst=0.0_rsh   l_sflx_radlossbot = .FALSE.
   ! nbatm=0    imeteo_dragtype=0    imeteo_exchtype=0    cds=0.0016_rsh
   ! paref=101500.0  ==> calcul de pa(:,:)=paref
   CALL cas_init_meteo
   wind_dir=90.0_rsh
   ! ****************** wind velocity and surface drag coefficient ********   
   wind_veloc=0.0_rsh
   cds=0.0011_rsh
   ! **********************************************************************   

#ifdef key_biolo
   sflx_solarcst=100.0_rsh
   sflx_rad_maxdepth=300.0_rsh
#endif

   ! 7- Temporal parameters
   ! ======================
   
!  **************** dates - time steps - start/end of simulaton**********  
   date_start=     '01/01/1900 00:00:00'
#ifdef key_biolo
   date_end=       '05/12/1901 12:00:00'
#else
   date_end=       '05/03/1900 12:00:00'
#endif
   date_start3d=date_start
   dtmin=60.0_rlg
   dtmax=600.0_rlg
   dt=120.0_rlg
   halfdt=dt*0.5_rlg
   suffix='1DV'
#ifdef key_biolo
#ifdef key_benthos
   suffix='1DV_biobenth'
#elif key_castest_sedbiomudb
   suffix='1DV_biosedb'
#else
   suffix='1DV_biosed'
#endif
#endif
#ifdef key_castest_1DVbioturb
   suffix='1DV_bioturb'
#endif
#ifdef key_castest_1DVconso
   suffix='1DV_conso'
#endif
#if defined key_castest_1DVbioturb && defined key_castest_1DVconso
   suffix='1DV_consoturb'
#endif

#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)
   dt_sed=60.0_rlg
   t_sed=tstart_sed+dt_sed
#endif
#endif
!  **********************************************************************   
   date_startagrif=date_start3d
   ! 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)
   ! over-assignment:
   iextmin(1)=2
   iextmax(1)=2
   jextmin(1)=2
   jextmax(1)=2
   fileresu(1)='results_'//TRIM(ADJUSTL(ADJUSTR(suffix)))//'.nc'
   

   ! 8- Saving parameters (results, and possible restart)
   ! ====================================================
   !    
   ! ******************************************************************* 
   pasor(1)=3600.0_rlg      ! time step for saving 
#ifdef key_biolo
   pasor(1)=7200.0_rlg
#endif
   date_startout=date_start3d
   date_endout=date_end
   tdebsor(1)=tool_datosec(date_startout)
   tfinsor(1)=tool_datosec(date_endout)
#ifdef key_biolo
   l_diag=.TRUE.
#else
   l_diag=.FALSE.
#endif
   l_out_buoy(1)=.TRUE.
   l_out_turb(1)=.TRUE.
   l_out_wind(1)=.FALSE.
   l_out_zaxis_w(1)=.TRUE.
   name_out_xe='XE'
   ! ***************************************************************** 
   ! possible restart from a saving file :
   !l_initfromfile=.FALSE.
   !file_init = 'save.nc'
   ! time step for saving file
   ! l_saverestart_1file=.FALSE.   
   ! l_saverestart_bydate=.FALSE.
   ! saverestart_step=5.0_rlg*3600.0_rlg   ! in hours
 

   ! 9- Initialiation of temperature and salinity  
   ! ============================================
   CALL cas_init_tspar
   ! over-assignment:   
   icon=1  ! no internal pressure gradient
!  ********************************************************************   
   qmaxz=10.0_rsh ! TVD scheme limiting parameter (modulate between upstream and quick advection schemes) vertical 
   l_stateeq_lin=.TRUE. 
   rhoref=1027.34_rsh
   saliref=35.0_rsh
   tetaref=10.0_rsh
   saliref_lin=saliref
   tetaref_lin=tetaref
!  ********************************************************************   


   ! 10- 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.

   ! ******************* possibility to read a file********************
   l_filesubs=.TRUE.
   IF (l_filesubs) THEN  ! we use "sub_read_var" instead of "cas_init_subs"
     name_filesubs = '../inputs/variable_test.dat'
#ifdef key_castest_1DVbioturb
     name_filesubs = '../inputs/variable_test_bioturb.dat'
#endif
#ifdef key_biolo
       ! par defaut : modele de base
#ifdef key_biolo_opt2
         name_filesubs = '../inputs/variable_sedbio_debase2.dat'
#else
       ! biolo option 1
#ifdef key_castest_sedbiomudb
         name_filesubs = '../inputs/variable_sedbio_debase1_mudb.dat'
#else
         name_filesubs = '../inputs/variable_sedbio_debase1.dat'
#endif
#ifdef key_benthos
         name_filesubs = '../inputs/variable_sedbio_benthos.dat'
#endif
#if defined key_karenia && defined key_psnz && defined key_phaeocystis && defined key_oxygen && defined key_benthos
         name_filesubs = '../inputs/variable_sedbio_previmer.dat'
#endif
#endif
#endif
!*****************************************************************   
   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
   
     ! 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
     
     ! give caracteristics of each variable 
     ! ------------------------------------
     ! typart=0: non-const. part., typart=1: const. part.
     ! itypv==1: GRAV, itypv==2: SAND, itypv==3: MUDS/PART, itypv==4: NCP,itypv==5: SORB,
     ! itypv==6: DISS, itypv==7: FIXE, itypv==8: DRIV, itypv==9: INTE, 

     !-------------------------------- beginning of substance definition
     !isubs=1   name_var(isubs)='Light Sand'
     !unit_var(isubs)='kg.m-3'  ! typart(isubs)=1     itypv(isubs)=2
     !ws_free_opt(2)=0     ws_free_min(2)=.0001_rsh     ws_free_max(2)=.0001_rsh     ws_free_para(1,2)=0.0_rsh
     ! ws_free_para(2,2)=0.0_rsh     ws_free_para(3,2)=0.0_rsh     ws_free_para(4,2)=0.0_rsh     ws_hind_opt(2)=0
     ! ws_hind_para(1,2)=.0_rsh     ws_hind_para(2,2)=.0_rsh 
     !tocd(isubs)=1.0_rsh     cini_wat(isubs)=0.0_rsh     cv_wat(isubs,:,:,:) = cini_wat(isubs)
     ! etc..
   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)
   ENDIF
   
#ifdef key_biolo
   !   VARIABLES DIAGNOSTIQUES
   ! ---------------------------------------
 ! lecture des parametres pour key_biolo ranges dans inputs comme le fichier de variables
   filevardiag='../inputs/vardiag_biodebase.dat'
#ifdef key_biolo_opt2
   filevardiag='../inputs/vardiag_biodebase_opt2.dat'
#else
#ifdef key_benthos
   filevardiag = '../inputs/vardiag_biobenthos.dat'
#endif
#if defined key_karenia && defined key_psnz && defined key_phaeocystis && defined key_oxygen && defined key_benthos
   filevardiag = '../inputs/vardiag_bioprevimer.dat'
#endif
#endif
   call sub_read_vardiag
   call biolo_userinit(cv_wat)
   call biolo_param('w')
#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

   ! 11- 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:
   ! introduction of an inital temperature vertical gradient 
   !! 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
   ! Miscellaneous  (ig,id...)
   CALL cas_init_misc
   
   
#ifdef key_sedim_MUSTANG

   ! 12- Initialization of sedimentary parameters
   ! ============================================
   
   !******************* before calling cas_init_sedim ********************************
#ifdef key_castest_1DVconso
   l_consolid = .TRUE.              ! option for accounting for consolidation
#else
   l_consolid = .FALSE.              ! option for accounting for consolidation
#endif
   hseduni = 0.10_rsh                 ! initial sediment height 
   !hseduni = 1.0_rsh                 ! initial sediment height 
! cseduni : concentration totale du sediment initial, a priori calculee en fonction d une concentration
!           relative de vase initiale; si on ne simule pas un culot de tassement, on part d un sediment deja consolide
!           typiquement : concentration relative de vase=600 kg/m3
!           on a cseduni=Ctot= crelmud/(1+f*(Crelmud/ros-1)) , 
!           ou f= fraction sableuse = [somme des conc. des sables et graviers]/ cseduni
! ici : f=0.4 et Crelvase=600 ==> Cseduni=867
   cseduni = 867.0_rsh               ! initial sediment concentration (ignored if l_csedinit_varz=.TRUE.)
   ksmiuni = 1                       ! default value of ksmi
   ksmauni = 75                      ! initial number of sediment layers
   dzsmaxinit = 0.5e-3               ! "theoretical" maximum layer thickness
   sinit_sed=35.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
   l_casinit_hsed=.TRUE.            ! .true. if we want to adjust the sediment thickness in order to be coherent with sediment parameters
                                     ! possible only if l_csedinit_varz=.FALSE.
   l_csedinit_varz=.FALSE.           ! true if csed varies within the sediment 
        ! if =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
        ! for example :
        ! sand
        !cv_sed(1,1:30,:,:)=500._rsh
        !cv_sed(2,1:30,:,:)=400._rsh
        !cv_sed(3,1:30,:,:)=500._rsh
        !cv_sed(1,31:60,:,:)=300._rsh
        !cv_sed(2,31:60,:,:)=250._rsh
        !cv_sed(3,31:60,:,:)=200._rsh
        !cv_sed(1,61:ksmauni,:,:)=200._rsh
        !cv_sed(2,61:ksmauni,:,:)=200._rsh
        !cv_sed(3,61:ksmauni,:,:)=50._rsh
        !cv_sed(-1,:,:,:)=10._rsh   ! temperature
        ! DO k=1,ksmauni
        !   cv_sed(-1,k,2,2)=(k+3)*0.5_rsh
        ! ENDDO
#ifdef key_biolo
#ifdef key_castest_sedmudb
 ! Ndet derniere variable constitutive
       ! cv_sed(nvpc,:,:,:)=0.0_rsh
#endif
       ! cv_sed(nvpc+1:nvp,:,:,:)=0.0_rsh
       ! cv_sed(iv_detr_N,:,:,:)=1000.0_rsh/hseduni
       ! cv_sed(iv_detr_Si,:,:,:)=1000.0_rsh/hseduni
       ! cv_sed(iv_detr_P,:,:,:)=200.0_rsh/hseduni
#endif
       ! DO iv=nvp+1,nv_state
       !         ! dissolved variables : cini_sed in conc/m3 of pore water
       !         ! + fixed variables : cini_sed in conc/m3 of sediment
       !   cv_sed(iv,:,:,:)=cini_sed(iv)
       ! ENDDO


       ! cv_sed( 0,:,:,:)=sinit_sed   ! salinite
        ! salinity and temperature in the sediment
       ! cv_sed(-1,:,:,:)=tinit_sed  ! temperature
        ! DO k=1,ksmauni
        !   cv_sed(-1,k,2,2)=(k+3)*0.5_rsh
        ! ENDDO
        ! cv_sed( 0,:,2,2)=35._rsh  ! salinite
   l_dzsinit_varz=.FALSE.            ! true if dzs (layer thickness) varies within the sediment
   IF(l_dzsinit_varz) THEN
        ! give  dzs of each initial sediment layer k (m) for i=j=2 [dzs(k,2,2)] 
        ! for example : dzs(1:10,:,2)=0.5_rsh
        !dzs(1:10,:,2)=0.015_rsh
        !dzs(11:19,:,2)=0.005_rsh
        !dzs(20,:,2)=0.002_rsh
     ! romaric
        !dzs(1:50,:,2)=0.002_rsh
        !dzs(51:75,:,2)=0.002_rsh
   ENDIF
   
    ! In case l_filesubs=.FALSE. define ros and diam of each particle 
    IF(.NOT.l_filesubs) THEN
       ros(1)=2600.0_rsh
       diam_sed(1)=0.000100_rsh
       ros(2)=2600.0_rsh
       diam_sed(2)=0.000020_rsh
      IF(.NOT. l_csedinit_varz) THEN
       ! and give uniform initial concentration of each variable in sediment
       ! ATTENTION cini_sed in fraction of total sediment for constitutive particulate variables
       !           cini_sed in mass/kg of total sediment for non constitutive ones
       !           cini_sed in mass/kg of sorbent sediment (associated constitutive particulaite variable) for SORB substances 
       !           cini_sed in mass/m3 of pore water for dissolved variables  
         cini_sed(1)=0.5       ! initial fraction [0-1] of total sediment concentration (cseduni)
         cini_sed(2)=0.5       ! initial fraction [0-1] of total sediment concentration (cseduni)
      ENDIF
    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)

   ! *************parameters to be checked (see parasedim.txt) ********
   dzsmin = 1.0e-5                  ! minimum layer thickness
   IF (l_consolid) THEN
     cfreshmud = 100.0_rsh            ! concentration of fresh muddy deposits
   ELSE
     cfreshmud = 500.0_rsh            ! concentration of fresh muddy deposits
     cfreshmud = 100.0_rsh      !test romaric      ! concentration of fresh muddy deposits
   ENDIF
#ifdef key_biolo
   l_diffused=.TRUE.                ! option for diffusion in sediment   
#else
   l_diffused=.FALSE.               ! option for diffusion in sediment   
#endif 
#ifdef key_castest_1DVbioturb
   l_diffused = .TRUE.              ! option for accounting for diffusion
   l_bioturb = .TRUE.              ! option for accounting for bioturbation
#endif  
   IF(l_diffused .or. l_consolid .OR. l_bioturb)l_dyn_insed = .TRUE.

   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
   xdifs1 = 1.e-8           ! diffusion coef wihtin the sediment
   xdifs2 = 0.0_rsh        ! diffusion coef wihtin the sediment
   xdifsi1 = 1.e-6          ! diffusion coef wihtin the water sediment interface
   xdifsi2 = 0.0_rsh        ! diffusion coef wihtin the water sediment interface
   epdifi = 0.01_rsh        ! diffusion thickness in the water at the sediment-water interface
   ! required if l_bioturb=.TRUE.
   xbioturbmax=1.157e-12_rsh  ! max diffusion coeffient by bioturbation Db (in surface)
   xbioturbk=6.0_rsh          ! coef (slope) for bioturbation coefficient betwwen max Db at sediemnt surface and 0 at bottom
   dbiotu0=0.01_rsh           ! max depth beneath the sediment surface below which there is no bioturbation
   dbiotum=0.002_rsh          ! sediment thickness wherein the diffusion-bioturbation coefficient Db is constant (max)
   frmud_db_min=0.6_rsh       ! mud fraction limit (min) below which there is no Bioturbation
   frmud_db_max=0.8_rsh       ! mud fraction limit (max)above which the bioturbation coefficient Db is maximum (muddy sediment)

! 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=.FALSE.
   IF(l_initsed_vardiss) THEN
        DO k=ksmiuni,ksmauni
             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


#ifdef key_sedim_MUSTANG
   csedmin=30.0_rsh
! csed_mud_ini : 
   csed_mud_ini=300.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.001_rsh            ! coef. for shear strength computation
   x2toce_mud=1._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.006_rsh            
   n_eros_sand=1.6_rsh           ! power of excess adim. shear stress in sand erosion law
   E0_mud=0.00015_rsh
   n_eros_mud=0.5_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=600._rsh               ! relative mud concentration in surficial layer below which fresh mud deposits can be mixed
   cvolmaxsort=0.61_rsh          ! maximum volumic concentration for sorted sediment (reference value= 0.58)
   cvolmaxmel=0.61_rsh           ! maximum volumic concentration for mixed sediment  (reference value= 0.67)
   MF=1.0_rsh                    ! morphological factor (not relevant for 1DV computation)
   htncrit_eros=0.0_rsh          ! water height below which no erosion flux is allowed
   xsigma1sg=xsigma1/9.81_rsh
   coef_erolat=0.0_rsh
   l_erolat_wet_cell=.FALSE.
#elif key_sedim_mudsand             
   csmaxmel = 1400.0_rsh
   frmudcr1=0.3_rsh
   frmudcr2=0.5_rsh
   E0sable=0.0002_rsh *0.3_rsh
   alphavase=1.0_rsh
   alphasable=1.0_rsh
   xeros = 0.0003_rsh *0.3_rsh   ! this formulation so as to have a similar erosion flux
                                 ! as for mixsed with x1toce_mud=0.3 and E0_mud=0.0002
   x1toce = 0.1_rsh
   x2toce = 0.0_rsh
   poro1 = 0.6_rsh
   xtas11 = 0.1e-3
   xtas12 = 50.0_rsh
   poro2 = 0.6_rsh
   xtas21 = 1.825e-5
   xtas22 = 57.0_rsh
   poro3 = 0.6_rsh
   xtas31 = 0.9e-08
   xtas32 = 13.0_rsh
   fexcs = 1.0_rsh
   cexcs=1._rsh-fexcs
   xcvgce=0.3_rsh
   fimpls=1.0_rsh
#endif
!*****************************************************************  
    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! initialisation forcage tenfon ou courant (lecture d un fichier ASCII resultats d un 3D en un point)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
#if defined key_1DV_force_towfile || defined key_1DV_force_ufile
#if defined key_1DV_force_towfile
   file_tenfonw_1DV='../inputs/tenfon_station_MANGAE2500/txt_files/122008_122009/tenfon_122008_122009_pt4_nj192_ni123_h0_27_8082.txt'
#elif defined key_1DV_force_ufile
   file_ubarotrop_1DV='../inputs/'
#endif 
   CALL cas1DV_read_forcing
#endif
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   fws2=fricwav*0.5_rsh  
!  if key_flocmod : caracterisation of flocs (see cas_init_sedim) 

!  **********choice output in sediment *******************************
!  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
!
#if defined key_sand2D
   l_outsandrouse(:)=.TRUE.
#endif
   CALL MUSTANG_compatibility
#endif
!  end of step 12 (sediment parameters)
  

!   13- adding a variable input water column  
!   ========================================
!
#ifdef key_substance
   ! **********************define number of inputs **********************
   ! one input (reject) per layer and per rejected variable 
   nbrej=0  
   IF(nbrej > 0) THEN
      CALL cas_init_outflow
      l_rejflux(:)=.TRUE.  ! no water input associated
                           ! then input as a flux (in quantity per second)
      ! example : here, mud input in surface ans sand input near bottom
      nucon(1)=2              ! Number of variable for input N1 (here SAND light)
      nucon(2)=3              ! Number of variable for input N1 (here mud)  
      krej(1)=1               ! layer of input 1
      krej(2)=kmax            ! layer of input 2
      ! fluxes rejected (in mass/s) :
      ! (first index : time variation, 2nd index : input number, see cas_init_outflow)                                                   
      valrej(:,1)=0.001_rsh    ! sable kg/s
      valrej(:,2)=0.005_rsh    ! vase  kg/s
   ! *********************************************************************
   ENDIF
#endif


   ! 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(*,*) 'initiatlisation :',cini_wat(iv)
   ENDDO
#if defined key_sedim_MUSTANG
     WRITE(*,*)'Number of substances :',nv_tot
     WRITE(*,*) 'cseduni (kg.m-3) =',cseduni
     DO iv=1,nv_tot
       WRITE(*,*) iv,': "',TRIM(name_var(iv)),'", cini_sed (fraction) =',cini_sed(iv)
       WRITE(*,*) iv,': "',TRIM(name_var(iv)),'", cini_wat (',TRIM(unit_var(iv)),') =',cini_wat(iv)
     END DO
#endif
#endif

   PRINT_DBG*, 'exit CAS_1DV general routine'
#endif
  END SUBROUTINE cas_1DV
