  SUBROUTINE cas_turbkatophil

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_turbkatophil  ***
!&E
!&E ** Purpose : definition of a test case to monitore the mixed layer depth
!&E              according to the vertical turbulence
!&E
!&E ** Description :
!&E                  1DV test case
!&E                  - Kato-Phillips (default) : wind-induced mixed layer
!&E                    The cooling of the water column is forced by a constant wind
!&E                    This test case involves the buoyancy and shear terms in the turbulence equations
!&E                  - Willis (-Dkey_castest_willis) : convective mixed layer
!&E                    The cooling of the water column is due to a negative solar flux
!&E                    that triggers an instability propagating in the column
!&E                  For both test cases, analytical solutions are available (threshold of tke)
!&E
!&E ** Called by : cas_init
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables : hx, hy, h0, sal, temp, dy, dhx, dhy, dsigu, dsigw
!&E                         surf, fisac, gwac, dgac, dfiac, fu, fv, sig, sigw
!&E                         tfiusr, tfivsr, goulet, govlet, ig, id, jb, jh
!&E                         ssh,u,v,uz,vz
!&E
!&E ** Reference :
!&E    Kato, H., Phillips, O.M., 1969. On the penetration of a turbulent layer into stratified fluid. Journal of Fluid Mechanics 37, 643-655
!&E    Willis, G. E., J. W. Deardorff, 1974: A Laboratory Model of the Unstable Planetary Boundary Layer. J. Atmos. Sci., 31, 1297-1307
!&E
!&E ** History :
!&E       !  2013-10  (S. Petton) Original code
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters, ONLY : rsh,rlg,riolg,valmanq,kmax,imin,imax,&
                          ljminm1,ljmaxp2,liminm1,limaxp2
   USE comvars2d,  ONLY : fisac,gwac,dtmin,dtmax,dt,halfdt,suffix
   USE comvars3d,  ONLY : nzinit,kzinit

   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,l_out_turb, &
                         l_out_buoy,l_out_wind,l_out_zaxis_w,   & 
                         iextmin,iextmax,jextmin,jextmax

   USE toolgeom,  ONLY : f_zu
   USE comvarp2d, ONLY : temp
   USE comvars3d, ONLY : tetaref
   USE comvars2d, ONLY : icon,fu,fv,h0
   USE cometdyn,  ONLY : wind_veloc ! ,wind_dir,rhoair,cds
   USE comvars2d, ONLY : l_obc_cycl_x
   USE obccombine,ONLY : l_obc_cycl
   USE comvars3d, ONLY : turb_0eq_option,turb_nbeq,turb_2eq_option,bgdiff,OptionSF
   USE cometemp,  ONLY : sflx_solarcst
   USE comdiag, ONLY : l_diag

   IMPLICIT NONE

!! * Arguments

!! * Local declarations
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   INTEGER           :: i,j,k
   INTEGER           :: ni_wall_e,ni_wall_w,nj_wall_n,nj_wall_s
   REAL(KIND=rsh)    :: dxval,dyval
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst,zr


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

   !PRINT*, 'enter CAS_TURBKATOPHIL 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
   icon=1


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

! assignement : input variables for cas_init_bathy routine
   depthcst = 50.0_rsh
   ni_wall_e = 1  ! land over ni_wall_e cells at the eastern part of the domain
   ni_wall_w = 2  ! land over ni_wall_w cells (from imin column) at the western part of the domain
   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



! flat bottom, plus walls

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




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

! assignement : input variables for cas_init_grid routine
 
! parametres de forme
!---------------------

   dxval = 340.0_rsh
   dyval = 340.0_rsh

   fisac=1.0_riolg
   gwac=1.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=   '02/01/1900 06:00:00'
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'
   dtmin=20.0_rlg
   dtmax=20.0_rlg
   dt=20.0_rlg
   suffix='kato_phillips'

   halfdt=dt*0.5_rlg

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

! over-assignement
   iextmin(1)=2
   iextmax(1)=2
   jextmin(1)=2
   jextmax(1)=2
   pasor(1)=3600.0_rlg
   l_out_sal(1)=.FALSE.
   l_out_temp(1)=.TRUE.
   l_out_buoy(1)=.TRUE.
   l_out_turb(1)=.TRUE.
   l_out_wind(1)=.TRUE.
   l_out_zaxis_w(1)=.TRUE.


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

! no open boundary
   CALL cas_init_obc

! over-assignement

   l_obc_cycl = .TRUE.
   l_obc_cycl_x=.TRUE. 



! viscosity   (none)
!-----------

! none
   CALL cas_init_visc

! over-assignement


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

   z0bottom=0.0035
   z0surface  =1.0_rsh
   nzinit=0.0001_rsh
   kzinit=0.0001_rsh

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

! over-assignement
   bgdiff=0.000001_rsh
   turb_nbeq=2
   turb_0eq_option=0
   turb_2eq_option=2
   OptionSF=5


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

! no forcing
   CALL cas_init_meteo

! over-assignement

   wind_veloc = SQRT(0.1_rsh)



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

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

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

   !PRINT*, 'exit CAS_TURBKATOPHIL routine'

  END SUBROUTINE cas_turbkatophil
