MODULE atke_turbulence_ini_mod

implicit none

save

integer :: iflag_atke, iflag_num_atke, iflag_atke_lmix
  !$OMP THREADPRIVATE(iflag_atke, iflag_num_atke, iflag_atke_lmix)
  real :: kappa = 0.4 ! Von Karman constant
  !$OMP THREADPRIVATE(kappa)
  real :: l0,ric,ri0,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke
  !$OMP THREADPRIVATE(l0,ric,cinf,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke)
  integer :: lunout,prt_level
  !$OMP THREADPRIVATE(lunout,prt_level)
  real :: rg, rd, rpi, rcpd, rv
  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd, rv)

  real :: viscom, viscoh
  !$OMP THREADPRIVATE(viscom,viscoh)

  real :: lmin=0.01              ! minimum mixing length
  !$OMP THREADPRIVATE(lmin)

  logical :: atke_ok_vdiff, atke_ok_virtual
  !$OMP THREADPRIVATE(atke_ok_vdiff,atke_ok_virtual)

CONTAINS

SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_in, rv_in)

   USE ioipsl_getin_p_mod, ONLY : getin_p

  integer, intent(in) :: lunout_in,prt_level_in
  real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in


  lunout=lunout_in
  prt_level=prt_level_in
  rd=rd_in
  rg=rg_in
  rpi=rpi_in
  rcpd=rcpd_in
  rv=rv_in

  viscom=1.46E-5
  viscoh=2.06E-5

  ! flag that controls options in atke_compute_km_kh
  iflag_atke=0
  CALL getin_p('iflag_atke',iflag_atke)

  ! flag that controls the calculation of mixing length in atke
  iflag_atke_lmix=0
  CALL getin_p('iflag_atke_lmix',iflag_atke_lmix)

  if (iflag_atke .eq. 0 .and. iflag_atke_lmix>0) then
        call abort_physic("atke_turbulence_ini", &
        'stationary scheme must use mixing length formulation not depending on tke', 1)
  endif

  ! activate vertical diffusion of TKE or not
  atke_ok_vdiff=.false.
  CALL getin_p('atke_ok_vdiff',atke_ok_vdiff)


  ! account for vapor for flottability
  atke_ok_virtual=.true.
  CALL getin_p('atke_ok_virtual',atke_ok_virtual)


  ! flag that controls the numerical treatment of diffusion coeffiient calculation
  iflag_num_atke=0
  CALL getin_p('iflag_num_atke',iflag_num_atke)

  ! asymptotic mixing length in neutral conditions [m]
  ! Sun et al 2011, JAMC 
  ! between 10 and 40

  l0=15.0 
  CALL getin_p('atke_l0',l0)

  ! critical Richardson number 
  ric=0.25
  CALL getin_p('atke_ric',ric)

  ! asymptotic value of Sm for Ri=-Inf
  cinf=1.5
  CALL getin_p('atke_cinf',cinf)

  ! constant for tke dissipation calculation
  cepsilon=5.87 ! default value as in yamada4
  CALL getin_p('atke_cepsilon',cepsilon)


 ! coefficient for surface TKE 
 ! following Lenderink & Holtslag 2004, ctkes=(cepsilon)**(2/3)
 ! (provided by limit condition in neutral conditions)
  ctkes=(cepsilon)**(2./3.)

  ! slope of Pr=f(Ri) for stable conditions
  pr_slope=5.0 ! default value from Zilitinkevich et al. 2005
  CALL getin_p('atke_pr_slope',pr_slope)
  if (pr_slope .le. 1) then
        call abort_physic("atke_turbulence_ini", &
        'pr_slope has to be greater than 1 for consistency of the tke scheme', 1)
  endif

  ! asymptotic turbulent prandt number value for Ri=-Inf
  pr_asym=0.4
  CALL getin_p('atke_pr_asym',pr_asym)

  ! value of turbulent prandtl number in neutral conditions (Ri=0)
  pr_neut=0.8
  CALL getin_p('atke_pr_neut',pr_neut)

  ! coefficient for mixing length depending on local stratification
  clmix=0.5
  CALL getin_p('atke_clmix',clmix)
 
  ! coefficient for mixing length depending on local wind shear
  clmixshear=0.5
  CALL getin_p('atke_clmixshear',clmixshear)


  ! minimum anisotropy coefficient (defined here as minsqrt(Ez/Ek)) at large Ri. 
  ! From Zilitinkevich et al. 2013, it equals sqrt(0.03)~0.17
  
  smmin=0.17
  CALL getin_p('atke_smmin',smmin)


  ! ratio between the eddy diffusivity coeff for tke wrt that for momentum
  ! default value from Lenderink et al. 2004
  cke=2.
  CALL getin_p('atke_cke',cke)

 RETURN

END SUBROUTINE atke_ini

END MODULE  atke_turbulence_ini_mod