  SUBROUTINE cas_obsthydroflume

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_obsthydroflume  ***
!&E
!&E ** Purpose : hydrodynamic experimental flume for fine calibration
!&E              of hydrodynamics in presence of obstructions (vegetation
!&E              or other).
!&E
!&E ** Description : The flume is by default a 2DV flume, 4 m long and 0.215 m
!&E                  deep. It is forced through a river discharge allowing to
!&E                  prescribe a constant velocity without feedbacks with bottom
!&E                  friction.
!&E                  This test-case only required the following CPP keys:
!&E                  -Dkey_castestriver , -Dkey_obstructions
!&E                  using -Dkey_dyn_adv_quick and -Dkey_tssub_adv_ultimatequickestmacho
!&E                  prevent generation of model instabilities at the edges of obstructions patch
!&E                  You also can use the CPP key : -Dkey_castestobstructions3D
!&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       !  2014-07 (F. Ganthy) main code
!&E       !  2016-09 (F. Ganthy) optimization on Abdelrhman 2007 method
!&E                              for computation on obstruction height :
!&E                              add number of obstruction segment to user-defined
!&E                              parameters
!&E       ! 2017-10 (F. Ganthy) Add new cpp keys related to new test-case options
!&E       ! 2017-10 (F. Ganthy) Some modifications:
!&E                             - Allowing multiple obstructions type in a single grid cell
!&E                             - Allowing 3D obstructions
!&E                             - Add simplified obstructions kind
!&E       ! 2017-11 (F. Ganthy) Remove useless variables related to turbulence coefficients
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters, ONLY : rsh,rlg,riolg,valmanq,kmax,imin,imax,jmin,jmax,lchain,&
                          ljminm1,ljmaxp2,liminm1,limaxp2,liminp1,limin,limax,ljmin,ljmax,&
                          liminm2,ljminm2
   USE comvars2d,  ONLY : fisac,gwac,dtmin,dtmax,dt,halfdt,suffix,fvisc,kx,ky,qmax,iscreenlog,z0b,l_modele2d,l_smagor
   USE comvars3d,  ONLY : nzinit,kzinit,qmaxz,l_stateeq_lin,saliref_lin

   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 : pasor,l_out_sal,l_out_temp,fileresu,tdebsor,tfinsor, &
                         iextmin,iextmax,jextmin,jextmax
   USE toolgeom,  ONLY : f_zu
   USE comvarp2d, ONLY : temp,sal,ssh,u,v
   USE comvars3d, ONLY : tetaref,saliref
   USE comvars2d, ONLY : icon,fu,fv,h0,tobs,out_dir,hminfrot
   USE comvarp3d, ONLY : uz,vz,ect,psi
   USE obccombine,ONLY : l_obc_mars,l_obc_east
   USE comvars2d, ONLY : h0_g,hx_g,hy_g,h0fond,hm,dhj,hrdef,vis2d,vish_xe,vish_phi
   USE comvars3d, ONLY : turb_nbeq,turb_2eq_option,bgdiff
   USE comturb, ONLY   : eps
   USE comdiag, ONLY : l_diag,file_diag,tdebdiag,tfindiag
   USE comfluxdo,  ONLY : nbfl,iriv,jriv,kriv,debmoy,tmfle,amtfle,l_tempriv,l_flestat,river_thickness

   USE comsubstance
   USE casinit,   ONLY : cas_init_subs,cas_init_subs_disspart,cas_init_outflow,cas_init_subs_budget

   USE parameters, ONLY : ksdmax
   USE comMUSTANG
   USE casinit,  ONLY : cas_init_sedim
   USE sedimorpho,  ONLY : sed_morphoinit



   USE comvars2d,ONLY : uforcing

   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,iobst
   INTEGER           :: ni_wall_e,ni_wall_w,nj_wall_n,nj_wall_s
   INTEGER           :: IERR_MPI
   REAL(KIND=rsh)    :: dxval,dyval,dzsmaxinit
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst,zr
   REAL(KIND=rlg)    :: tool_datosec
   LOGICAL           :: l_filesubs,l_nc4par







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

     !PRINT*, 'enter CAS_OBSTHYDROFLUME3D routine'







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

! dynamics
   CALL cas_init_param

! over-assignement
   l_modele2d = .FALSE.

! 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
   IF(l_modele2d) THEN
     icon=0
   ELSE
     icon=2
   ENDIF
   saliref = 35.50_rsh
   tetaref = 10.0_rsh
   l_stateeq_lin=.FALSE.
   saliref_lin=saliref

! Bathymetry
! ----------

! assignement : input variables for cas_init_bathy routine
   ni_wall_w = 2  ! land over ni_wall_w cells
   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

   h0fond   = 0.01_rsh
   hm       = 0.01_rsh


   depthcst = 1.0_rsh


   dhj      = 0.1

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

! over_assignement


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

! assignement : input variables for cas_init_grid routine


   dxval = 0.05_rsh
   dyval = 0.05_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'

   date_end        = '01/01/1900 01:00:00'

   date_start3d    = '01/01/1900 00:00:00'
   date_startagrif = '01/01/1900 00:00:00'

   dtmin  = 0.002_rlg
   dtmax  = 0.002_rlg
   dt     = 0.002_rlg

   suffix = 'obsthydroflume2dv'
! If specific test-cases, overwrite suffix









   halfdt=dt*0.5_rlg


   date_startsub='01/01/1900 00:05:00'
   tdebsubs=tool_datosec(date_startsub)

! 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


   l_nc4par = .TRUE.
   CALL cas_init_io(date_start,date_end,date_start3d,date_startagrif,suffix,l_nc4par)

! over-assignement
   iextmin(1) = imin
   iextmax(1) = imax

   jextmin(1) = 2
   jextmax(1) = 2


   pasor(1)   = 10.0_rlg

   tobs       = 2.0_rlg*3600.0_rlg
   out_dir    = './'

! Name of results file :
! =========================================











   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)
   l_diag        = .FALSE.

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

! no open boundary
   CALL cas_init_obc

! over-assignement
! geometry (where are the open boundaries ?)
   l_obc_east = .TRUE. 
! used scheme to manage obc
   l_obc_mars = .TRUE.


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

! none
   CALL cas_init_visc

! over-assignement
   l_smagor        = .TRUE.
   fvisc           = 1.0_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.00005_rsh

   z0surface = 0.00005_rsh
   IF (l_modele2d) THEN
     nzinit    = 0.0001_rsh
     kzinit    = 0.0001_rsh
   ELSE
     nzinit    = 0.0001_rsh
     kzinit    = 0.0001_rsh
   ENDIF
   
   CALL cas_init_turb(kzinit,nzinit,z0surface,z0bottom)

! over-assignement
   IF(.NOT.l_modele2d) THEN
     bgdiff          = 0.000001_rsh
     turb_nbeq       = 2
     turb_2eq_option = 2
     hminfrot        = 0.001_rsh
     ect(:,:,:)      = 1.0e-8
     psi(:,:,:)      = 1.0e-12
     eps(:,:,:)      = 1.0e-12
   ENDIF

! 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
   IF(.NOT.l_modele2d) THEN
     sal(:,:,:)  = 0.0_rsh
     temp(:,:,:) = 15.0_rsh
     WHERE( h0(liminm2:limaxp2,ljminm2:ljmaxp2)==-valmanq ) ssh(liminm2:limaxp2,ljminm2:ljmaxp2)=valmanq
     DO k=1,kmax
       WHERE( h0(liminm1:limaxp2,ljminm1:ljmaxp2)==-valmanq ) sal(k,liminm1:limaxp2,ljminm1:ljmaxp2)=valmanq
       WHERE( h0(liminm1:limaxp2,ljminm1:ljmaxp2)==-valmanq ) temp(k,liminm1:limaxp2,ljminm1:ljmaxp2)=valmanq
     END DO 
   ENDIF

   uforcing    = 0.5_rsh ! (m.s-1) depth averaged velocity for initialization and used to compute river discharge


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

   nv_grav=0
   nv_sand=1
   nv_mud=1
   nv_cp=0
   nv_ncp=0
   nv_sorb=0
   nv_dis=0
   nv_fix=0
   nv_dri=0
   nv_int=0
   isubs=1
   name_var(isubs)='sable_fin'
   unit_var(isubs)='kg.m-3'
   isubs=2
   name_var(isubs)='vase'
   unit_var(isubs)='kg.m-3'

   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

!===========================DEFINITION DES VARIABLES======================================
! first variable (sable fin)
!-----------------
   itypv(1) = 2 
   typart(1)=1.
   ws_free_opt(1)=0
   ws_free_min(1)=.005_rsh
   ws_free_max(1)=.005_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
   t90(1)=0.0_rsh
   sub_flx_atm(1)=0.0_rsh
   cv_rain(1)=0.0_rsh
   cv_wat(1,:,:,:) = cini_wat(1)
   obc_cv_wat_w(1,:,:,:)=cini_wat(1)

! second variable (vase)
!-----------------------
   itypv(2) = 3
   typart(2)=1.
   ws_free_opt(2)=0
   ws_free_min(2)=.0005_rsh
   ws_free_max(2)=.0005_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(2)=100.0_rsh
   t90(2)=0.0_rsh
   sub_flx_atm(2)=0.0_rsh
   cv_rain(2)=0.0_rsh
   cv_wat(2,:,:,:) = cini_wat(2)
   obc_cv_wat_w(2,:,:,:)=cini_wat(2)
!================================ FIN DEFINITION DES VARIABLES =====================

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

      CALL cas_init_subs_disspart(l_filesubs)

! over-assignement


      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

!First : sable fin
     ros(1)=2600.0_rsh
     diam_sed(1)=0.0001_rsh
!Second : mud
     ros(2)=2600.0_rsh
     diam_sed(2)=0.0001_rsh

   ENDIF


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



! Initialization of sedimentary parameters
! ----------------------------------------
! pre-assignement
   l_consolid = .FALSE.
   hseduni = 0.1_rsh
   ksmiuni = 1
   ksmauni = 2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! 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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   cseduni = 1508.0_rsh
   c_sedtot(:,:,:)=cseduni
   cini_sed(1)=0.2_rsh
   cini_sed(2)=0.8_rsh
   dzsmaxinit=10.0_rsh

! sediment output sur toutes les couches ---option 1---
   nk_nivsed_out=ksdmax
   ep_nivsed_out(:)=0.0_rsh
! Initialization of sediment layers
! ---------------------------------
   CALL cas_init_sedim(l_filesubs,dzsmaxinit=dzsmaxinit)

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


   IF(l_morphocoupl) CALL sed_morphoinit

! test sortie des couches
!   dzsmaxuni = 0.02_rsh


! fin sedim




! Obstructions
!--------------


   !PRINT*, 'exit CAS_OBSTHYDROFLUME2DV routine'


  END SUBROUTINE cas_obsthydroflume
