!!======================================================================

  SUBROUTINE cas_innershelf

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_innershelf  ***
!&E
!&E ** Purpose : define the innershelf test case
!&E
!&E ** Description :
!&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       !  2005-10-13  Original code
!&E       !  2008-12    (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E
!&E---------------------------------------------------------------------





























































!! * 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,liminm1,ljminm1,valmanq
   USE comvars2d, ONLY : fisac,gwac,fu,fv,dtmin,dtmax,dt,halfdt,suffix,l_obc_cycl_x, &
                         l_obc_cycl_y,icon,qmax,t,h0,ig,l_modele2d,fu,dy_v, &
                         h0_g,hx_g,hy_g,vis2d,vish_phi,vish_xe,frot_lin,hmaxfrot, &
                         hminfrot,vismin,vismax,cosmag,l_smagor,l_spongw,sponge_nbcell
   USE comvarp2d, ONLY : ssh,temp
!USE comvarp3d, ONLY : uz,vz
   USE comvars3d, ONLY : kzinit,nzinit,qmaxz,saliref,rhoref,nz,l_stateeq_lin, &
                         turb_0eq_option,bgdiff,turb_nbeq,turb_2eq_option
   USE obccombine, ONLY : l_obc_cycl
   USE toolmpi, ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh
   USE output,    ONLY : pasor,l_out_sal,l_out_wind,l_out_turb,l_out_zaxis_w, &
                         l_out_visc,nbave,l_out_timeave

   USE toolgeom,  ONLY : f_zw,f_zu
   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 cometdyn,  ONLY : wind_veloc,windamp_duration,pa


   IMPLICIT NONE

!! * Arguments

!! * Local declarations
   LOGICAL           :: l_nc4par
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   INTEGER           :: i,j,k
   INTEGER           :: ni_wall_w,ni_wall_e
   REAL(KIND=rsh)    :: dxval,dyval
   REAL(KIND=rsh)    :: vismuloc,nu0
   REAL(KIND=rsh)    :: x,z,zsurf
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst
   REAL(KIND=rsh)    :: Length_ETA,Length_XI,cff
   REAL(KIND=rsh)    :: Ekdepth,zd0,zd1,strat,thknss
!   REAL(KIND=rsh), DIMENSION(imin:imax) :: Eknum


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

   !PRINT*, 'enter CAS_INNERSHELF 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
   saliref=35.0_rsh
   rhoref=1000.0_rsh
   icon=0


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

! assignement : input variables for cas_init_bathy routine
   depthcst = 0.0_rsh
   ni_wall_w = 2  ! land over ni_wall_w cells (from jmin row) at the southern part of the domain
   ni_wall_e = 1  ! land over ni_wall_e cells at the northern part of the domain


! As depthcst is null, specify your own bathymetry
   dxval = 1000.0_rsh

   DO i=imin,imax
     x=dxval*(REAL(i-2,rsh)+0.5_rsh)  ! -2 pour se caler avec Roms
!x=dxval*REAL(i-2,rsh)
     h0_g(i,:)=0.001_rsh*(200000.0_rsh-x)+4.0_rsh
     hy_g(i,:)=h0_g(i,:)
     x=dxval*REAL(i-1,rsh)  ! -2 pour se caler avec Roms
     hx_g(i,:)=0.001_rsh*(200000.0_rsh-x)+4.0_rsh
   END DO


! walls
   CALL cas_init_bathy(depthcst,ni_wall_w=ni_wall_w,ni_wall_e=ni_wall_e)

! over-assignement


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

! assignement : input variables for cas_init_grid routine
 
! parametres de forme
!---------------------
   dxval = 1000.0_rsh
   dyval = 1000.0_rsh
   fisac=-21.0_riolg
   gwac=1.0_riolg



   CALL cas_init_grid(dxval,dyval)

! over-assignement


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


! assignement : input variables for cas_init_io routine
   date_start=   '01/01/1900 00:00:00'
!date_end=   '02/01/1900 00:00:00'
   date_end=   '21/01/1900 00:00:00'
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'
   dtmin=100.0_rlg
   dtmax=100.0_rlg
   dt=100.0_rlg
   halfdt=dt*0.5_rlg

   suffix='innershelf'
   l_nc4par = .TRUE.

   CALL cas_init_io(date_start,date_end,date_start3d,date_startagrif,suffix,l_nc4par=l_nc4par,nb_file=2)

! over-assignement
   pasor(2)=3600.0_rlg
   pasor(2)=86400.0_rlg*20.0_rlg  ! 20 days  (instantaneous field)
!pasor(1)=86400.0_rlg
!pasor(1)=78720.0_rlg
   pasor(1)=86400.0_rlg*2.0_rlg  ! 2 days  (averaged field)
   l_out_sal(:)=.FALSE.
   l_out_wind(:)=.TRUE.
   l_out_visc(:)=.TRUE.
   l_out_turb(:)=.TRUE.
   l_out_zaxis_w(:)=.TRUE.
   l_out_timeave(1)=.TRUE.
   nbave=1


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

! no open boundary
   CALL cas_init_obc

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

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

! none
   CALL cas_init_visc

! over-assignement
!vish_xe(:,:,:) = 100.0_rsh
!vish_phi(:,:,:)= 100.0_rsh
!l_smagor=.TRUE.
!cosmag=0.2_rsh
!vismin=10.0_rsh
!vismax=700.0_rsh
!cosmag=0.1_rsh
!vismin=0.0_rsh
!vismax=700.0_rsh
   l_spongw=.TRUE.
   sponge_nbcell=10
   DO j=ljminp1,ljmax
     IF(l_spongw) THEN                        ! couche eponge a l ouest si ouverte
       IF(h0_g(ig(j),j).GT.-valmanq) THEN
         DO i=ig(j),MAX(MIN(ig(j)+sponge_nbcell,imax),imin)
           IF ((i >= limin).AND.(i <= limax)) THEN
!     vismuloc=(5.* exp(-real(i-1)/10.))
!     vismuloc=(i-(ig(j)+sponge_nbcell))*(1-vismul)/sponge_nbcell+1.0_rsh
           vismuloc=COS(pi/2.0_rsh*(i-ig(j))/sponge_nbcell)
           nu0=100.0_rsh
           vis2d(i,j)=MAX(nu0*vismuloc,vis2d(i,j))
           IF(.NOT.l_modele2d) THEN
             vish_phi(:,i,j)=MAX(nu0*vismuloc,vish_phi(:,i,j))
             vish_xe(:,i,j) =MAX(nu0*vismuloc,vish_xe(:,i,j))
           END IF
           ENDIF
         END DO
       END IF
     END IF
   END DO



! turbulence
! ----------

   z0bottom   =0.0035_rsh  ! not used because the bottom stress is linear and nz=cst, kz=cst
   z0surface  =0.0_rsh
!z0surface  =1.6_rsh
!z0surface  =0.6_rsh
   nzinit=1.0e-2
   kzinit=0.0_rsh

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

! over-assignement
   turb_0eq_option=1
   turb_nbeq=0

!hmaxfrot=200.0_rsh
   hminfrot=20.0_rsh
   hminfrot=0.0_rsh

! Ekman depth
!   Ekdepth = pi* SQRT( 2*nzinit/ABS(fu(1,1)) )
! Ekman number
!   Eknum = nzinit / ABS(fu(1,1)) / hx_g(:,1)**2


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

! no forcing
   CALL cas_init_meteo

! over-assignement
   wind_veloc = SQRT(0.07_rsh)  ! to get a wind stress of 0.07N/m2, assuming rhoair=1
   windamp_duration = 0.0_rlg ! None

!! Parallelisation pas prete ici
   DO j=ljminp1,ljmax
     pa(limin:limax,j) = pa(limin:limax,j-1) - fu(limin:limax,j)*rhoref*0.02_rsh*dy_v(limin:limax,j)
   END DO
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,pa(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,pa(liminm2:limaxp2,ljminm2:ljmaxp2))


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

! no river
   CALL cas_init_river


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

! set to 0,tetaref,saliref
   CALL cas_init_progvar

! over-assignement

! tracers
!--------




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

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

   !PRINT*, 'exit CAS_INNERSHELF routine'

  END SUBROUTINE cas_innershelf

!!======================================================================

  SUBROUTINE cas_innershelf_structure(h,D)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_innershelf_structure  ***
!&E
!&E ** Purpose : define the horizontally-varying structure which extends
!&E              1D Ekman s model to horizontal dimensions using Welander s method
!&E
!&E ** Description :
!&E
!&E ** Called by : cas_init
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :  Welander (1975), Brink (1983), Mitchum and Clarke (1986)
!&E
!&E ** History :
!&E       !  2012-07    (V. Garnier)  Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : imin,imax,pi,rsh
   USE comvars2d,  ONLY : fu

!! * Arguments
   REAL(rsh), DIMENSION(imin:imax), INTENT(in) :: h    ! topography
   REAL(rsh),                       INTENT(in) :: D    ! Ekman depth

!! * Local declarations
   REAL(KIND=rsh), DIMENSION(imin:imax) :: alpha,S1,S2,T1,T2   !  struture functions
!   S1p,S2p,T1p,T2p        !  derivative of struture functions
!Eknum                  !  Ekman number
   INTEGER                :: i
   REAL(rsh)              :: phd,hemis

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

   hemis = fu(1,1)/ABS(fu(1,1))

   DO i=imin,imax
     phd = pi*h(i)/D
     alpha(i) = ( cosh(phd)*cos(phd) )**2 + ( sinh(phd)*sin(phd) )**2
     S1(i) = cosh(phd)*cos(phd) / alpha(i)
     S2(i) = sinh(phd)*sin(hemis*phd) / alpha(i)
     T1(i) = cosh(phd)*sinh(phd) / alpha(i)
     T2(i) = cos(phd)*sin(hemis*phd) / alpha(i)
   END DO


  END SUBROUTINE cas_innershelf_structure
