  SUBROUTINE cas_obstsedflumecycl

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_obstsedflumecycl  ***
   !&E
   !&E ** Purpose : test case for flume experiments with cyclic boundary conditions
   !&E              sediment and obstructions
   !&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 force applied on momentum equation
   !&E                  This test-case only required the following CPP keys:
   !&E                  -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-08-11  Original code (F. Ganthy) adapted from caskelvin
   !&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---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : rsh,rlg,riolg,kmax,grav,pi,imin,imax,jmin,jmax, &
                          limin,liminp1,limax,ljmin,ljminp1,ljminp2,ljmax, &
                          ljmaxm1,ljminm2,ljmaxp2,liminm2,limaxp2, &
                          limaxp1,ljmaxp1,valmanq, &
                          ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north
   USE comvars2d, ONLY : fisac,gwac,fu,fv,dtmin,dtmax,dt,halfdt,suffix,l_obc_cycl_x, &
                         l_obc_cycl_y,icon,qmax,t,l_modele2d,h0,hx,h0fond, &
                         jh,jb,ig,id,hminfrot,kx,ky,hm,dhj,fvisc,out_dir,vis2d,vish_xe,vish_phi
   USE comvarp2d, ONLY : ssh,u,v
   USE comvarp3d, ONLY : uz,vz,ect,psi
   USE comvars3d, ONLY : kzinit,nzinit,qmaxz,turb_nbeq,turb_2eq_option,l_stateeq_lin,saliref_lin, &
                         saliref,bgdiff
   USE comturb, ONLY   : eps
   USE obccombine, ONLY : l_obc_cycl
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh
   USE output,    ONLY : pasor,l_out_temp,fileresu,tdebsor,tfinsor, &
                         iextmin,iextmax,jextmin,jextmax
#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,               &
#ifdef key_obstructions
                         cas_init_obstructions,                        &
#endif
                         cas_init_river
   USE comfluxdo, ONLY : phieau
   USE comvars2d, ONLY : uforcing
#ifdef key_obstructions
  USE obstructions
#endif


   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif,date_startout,date_endout
   INTEGER           :: i,j,iobst,k
   INTEGER           :: nj_wall_n,nj_wall_s
   REAL(KIND=rsh)    :: dxval,dyval
   REAL(KIND=rsh)    :: fuval
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst
   REAL(KIND=rlg)    :: tool_datosec
#ifdef key_castestobstructions3D
   ! FOR PATCHY DISTRIBUTION
   INTEGER                                 :: nb_patch
   REAL(KIND=rsh)                          :: diam_patch,distp
   REAL(KIND=rsh),DIMENSION(limin:limax)   :: xxp
   REAL(KIND=rsh),DIMENSION(ljmin:ljmax)   :: yyp
   REAL(KIND=rsh),DIMENSION(:),ALLOCATABLE :: xpatch,ypatch
#endif

   !!----------------------------------------------------------------------
   !! * Executable part
#ifdef key_castestobstructions3D
   PRINT_DBG*, 'enter CAS_OBSTSEDFLUMECYCL3D routine'
#else
     PRINT_DBG*, 'enter CAS_OBSTSEDFLUMECYCL2DV routine'
#endif

   ! 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
   icon=2
   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
#ifdef key_castestobstructions3D
   depthcst = 1.0_rsh
#else
   depthcst = 0.215_rsh
#endif
   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
   h0fond   = 0.01_rsh
   hm       = 0.01_rsh
   dhj      = 0.1_rsh


   ! flat bottom, plus walls
   CALL cas_init_bathy(depthcst,nj_wall_n=nj_wall_n,nj_wall_s=nj_wall_s)

   ! over-assignement

   ! grid parametrization
   ! --------------------

   ! assignement : input variables for cas_init_grid routine
#ifdef key_castestobstructions3D
   dxval=0.1_rsh
   dyval=0.1_rsh
#else
   dxval=0.05_rsh
   dyval=0.05_rsh
#endif
   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 00:30:00'
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'
#ifdef key_castestobstructions3D
   dtmin=0.05_rlg
   dtmax=0.05_rlg
   dt=0.05_rlg
   suffix='obstsedflumecycl3d'
#else
   dtmin=0.0025_rlg
   dtmax=0.0025_rlg
   dt=0.0025_rlg
   suffix='obstsedflumecycl2dv'
#endif
   halfdt=dt*0.5_rlg

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

   ! over-assignement
   iextmin(1) = imin
   iextmax(1) = imax
#ifdef key_castestobstructions3D
   jextmin(1) = jmin
   jextmax(1) = jmax
#else
   jextmin(1) = 2
   jextmax(1) = 2
#endif
   pasor(1)=10.0_rlg
! 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)

   ! open boundaries
   ! ---------------

   ! no open boundary
   CALL cas_init_obc

   ! over-assignement
   ! used scheme to manage obc
   l_obc_cycl = .TRUE.
   l_obc_cycl_x=.TRUE.
   ! physics introduced into the simulated area (tide or ogcm solution)

   ! viscosity
   !-----------

   ! 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   =1.0e-5
   z0surface  =1.0e-5_rsh
   kzinit=0.000001_rsh
   nzinit=0.000001_rsh

   ! Gaspard parametrization
   CALL cas_init_turb(kzinit,nzinit,z0surface,z0bottom)

   ! over-assignement
   bgdiff          = 1.0e-6
   turb_nbeq       = 2
   turb_2eq_option = 2
   hminfrot        = 0.01_rsh
   ect(:,:,:)      = 1.0e-8
   psi(:,:,:)      = 1.0e-12
   eps(:,:,:)      = 1.0e-12
   phieau(:,:,:)   = 0.0_rlg


   ! meteorology
   ! -----------

   ! no forcing
   CALL cas_init_meteo


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

   ! no river
   CALL cas_init_river


   ! initialization of the pronostic variables
   ! -----------------------------------------

   ! set to 0,tetaref,saliref
   CALL cas_init_progvar


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

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

   ! over-assignement
   id(:)=id(:)-ishift_obc_east
   ig(:)=ig(:)+ishift_obc_west
   jh(:)=jh(:)-jshift_obc_north
   jb(:)=jb(:)+jshift_obc_south

   ! Dynamics : forcing
#ifdef key_casobstflume_cycle
      ! Velocity forcing for flume experiments
   uforcing = 0.12_rsh/depthcst
#endif  


    ! Obstructions
   !--------------
#ifdef key_obstructions
   ! Number of obstruction variables
   obst_nbvar = 1
   CALL obst_alloc_nbvar
   ! Initialisation of each obstruction variables (the following lines should be repeated for 
   ! each variable
   ! FIRST VARIABLE
   iobst=1
   obst_varnum(iobst)               = iobst             ! * (1)  VARIABLE NUMBER
   obst_varname(iobst)              = 'ZosteraNoltei'   ! * (2)  VARIABLE NAME
   l_obst_flexible(iobst)           = .TRUE.            ! * (3)  VARIABLE IS RIGID OR FLEXIBLE
   l_obst_downward(iobst)           = .FALSE.           ! * (4)  VARIABLE IS UPWARD OR DOWNWARD
   l_obst_cylindre(iobst)           = .FALSE.           ! * (5)  LOGICAL FOR SHAPE OF OBSTRUCTIONS
   l_obst_3dobst(iobst)             = .FALSE.           ! * (6)  LOGICAL FOR FULL 3D OBSTRUCTIONS
   l_obst_noturb(iobst)             = .FALSE.           ! * (7)  LOGICAL FOR USE SIMPLIFIED FORMULATION (Z0)
   l_obst_filechar(iobst)           = .FALSE.           ! * (8)  LOGICAL FOR USE OF TEMPORAL OBSTRUCTION CHARACTERISTICS
   obst_fn_char(iobst)              = './'              ! * (9)  FILENAME FOR TEMPORAL OBSTRUCTION CHARACTERISTICS
   l_obst_filedistri(iobst)         = .FALSE.           ! * (10) LOGICAL FOR USE A FILE FOR VERTICAL DISTRIBUTION OF OBSTRUCTIONS
   obst_fn_distrib(iobst)           = './'              ! * (11) FILENAME FOR VERTICAL DISTRIBUTION FILE
   l_obst_fracxy(iobst)             = .TRUE.            ! * (12) LOGICAL TO TAKE ACCOUNT FOR COVERAGE CORRECTION WITHIN GRID CELL
   obst_fracxy_type(iobst)          = 2                 ! * (13) INTEGER FOR THE TYPE OF COVERAGE CORRECTION
   obst_c_fracxy_k0(iobst)          = 6.0_rsh           ! * (14) FIRST PARAMETER FOR COVERAGE CORRECTION (TYPE=1 or TYPE=2)
   obst_c_fracxy_k1(iobst)          = 12.0_rsh          ! * (15) SECOND PARAMETER FOR COVERAGE CORRECTION (TYPE=2)
   obst_c_fracxy_l(iobst)           = 3.5_rsh           ! * (16) THIRD PARAMETER FOR COVERAGE CORRECTION (TYPE=2)
   obst_i_height(iobst)             = 0.15_rsh          ! * (17) INITIAL OBSTRUCTION HEIGHT
   obst_i_width(iobst)              = 0.001_rsh         ! * (18) INITIAL OBSTRUCTION WIDTH
   obst_i_thick(iobst)              = 0.0002_rsh        ! * (19) INITIAL OBSTRUCTION THICKNESS
   obst_i_dens(iobst)               = 50000.0_rsh       ! * (20) INITIAL OBSTRUCTION DENSITY
   l_obst_abdelposture(iobst)       = .FALSE.           ! * (21) USE ABDELRHMAN 2007 METHOD FOR OBSTRUCTION POSTURE
   obst_i_abdel_nmax(iobst)         = 3                 ! * (22) NUMBER OF OBSTRUCTION HEIGHT DISCRETIZATION FOR ABDELRHMAN 2007 METHOD
   obst_i_rho(iobst)                = 650.0_rsh         ! * (23) INITIAL OBSTRUCTION VOLUMETRIC MASS
   l_obst_param_height(iobst)       = .TRUE.            ! * (24) USE PARAMETERIZATION (FROM VELOCITY) FOR FLEXIBLE OBSTRUCTION HEIGHT
   obst_c_height_x0(iobst)          = 0.95              ! * (25) FIRST COEFFICIENT FOR OBSTRUCTION HEIGHT
   obst_c_height_x1(iobst)          = -3.6              ! * (26) SECOND COEFFICIENT FOR OBSTRUCTION HEIGHT
   l_obst_drag_cste(iobst)          = .TRUE.            ! * (27) USE CONSTANT DRAG COEFFICIENT
   obst_i_drag(iobst)               = 1.5               ! * (28) DRAG COEFFICIENT FOR OBSTRUCTION
   obst_i_lift(iobst)               = 0.5               ! * (29) LIFT COEFFICIENT FOR OBSTRUCTIONS
   obst_c_shelter(iobst)            = 100.0_rsh         ! * (30) SCALING FACTOR FOR SHELETERING COEFFICIENT
   obst_c_lz(iobst)                 = 0.1               ! * (32) FORM COEFFICIENT FOR DISSIPATION TIME-SCALE WITHIN OBSTRUCTIONS
   l_obst_z0bstress(iobst)          = .TRUE.            ! * (33) IMPACT OF OBSTRUCTIONS ON SEDIMENTOLOGICAl ROUGHNESS LENGTH
   l_obst_z0bstress_strick(iobst)   = .FALSE.           ! * (34) USE OF Z0 CORRECTION FROM STRICKLER COEFFICIENT
   l_obst_z0bstress_param3d(iobst)  = .TRUE.            ! * (35) USE OF Z0 CORRECTION FROM PARAMETERIZATION FOR 3D FORMULATIONS
   obst_i_z0bstress3d(iobst)        = 0.0012_rsh        ! * (36) INITIAL VALUE FOR SEDIMENTOLOGICAL ROUGHNESS LENGTH
   obst_c_z0bstress3d_x0(iobst)     = 0.0011_rsh        ! * (37) FIRST PARAMETER FOR ROUGHNESS LENGTH PARAMETERIZATION FOR 3D FORMULATIONS
   obst_c_z0bstress3d_x1(iobst)     = 0.2974_rsh        ! * (38) SECOND PARAMETER FOR ROUGHNESS LENGTH PARAMETERIZATION FOR 3D FORMULATIONS
   l_obst_z0bstress_param2dh(iobst) = .FALSE.           ! * (39) USE OF Z0 CORRECTION FROM PARAMETERIZATION FOR 2Donly FORMULATIONS
   obst_i_z0bstress2dh(iobst)       = 0.0042_rsh        ! * (40) INITIAL VALUE FOR SEDIMENTOLOGICAL ROUGHNESS LENGTH FOR 2D FORMULATIONS
   obst_c_z0bstress2dh_x0(iobst)    = 0.0001_rsh        ! * (41) FIRST PARAMETER FOR ROUGHNESS LENGTH PARAMETERIZATION FOR 2D FORMULATIONS
   obst_c_z0bstress2dh_x1(iobst)    = -1.3348_rsh       ! * (42) SECOND PARAMETER FOR ROUGHNESS LENGTH PARAMETERIZATION FOR 2D FORMULATIONS
   ! END of Initialization of first obstruction variable

   ! The roughness length for bottom shear stress (outside obstructions)
   obst_i_z0bstress = 0.00024_rsh

   ! Finalizing initialization
   CALL cas_init_obstructions

   ! over-assignement
   ! Position of the fist variable over the domain
#ifdef key_castestobstructions3D
   ! Position of the fist variable over the domain
   DO i=imin,imax
     xxp(i) = REAL(i,rsh)*dxval
   ENDDO
   DO j=jmin,jmax
     yyp(j) = REAL(j,rsh)*dyval
   ENDDO
   nb_patch = 1
   ALLOCATE(xpatch(1:nb_patch))
   ALLOCATE(ypatch(1:nb_patch))
   diam_patch = 0.4_rsh
   xpatch(1) = 2.0_rsh
   ypatch(1) = MAXVAL(yyp)/2.0_rsh
   DO k=1,nb_patch
     DO i=imin,imax
       DO j=jmin,jmax
         distp = sqrt((xxp(i)-xpatch(k))**2.0_rsh + (yyp(j)-ypatch(k))**2.0_rsh)
         IF (distp.LE.(diam_patch/2.0_rsh)) obst_position(iobst,i,j) = 1.0_rsh
       ENDDO
     ENDDO
   ENDDO

#else
   obst_position(iobst,36:53,:) = 1.0_rsh
#endif

#endif

#ifdef key_castestobstructions3D
   PRINT_DBG*, 'exit CAS_OBSTSEDFLUMECYCL3D routine'
#else
     PRINT_DBG*, 'exit CAS_OBSTSEDFLUMECYCL2DV routine'
#endif


  END SUBROUTINE cas_obstsedflumecycl
