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
  !$OMP THREADPRIVATE(l0, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, clmix)
  integer :: lunout,prt_level
  !$OMP THREADPRIVATE(lunout,prt_level)
  real :: rg, rd, rpi, rcpd
  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd)

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

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


CONTAINS

SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_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


  lunout=lunout_in
  prt_level=prt_level_in
  rd=rd_in
  rg=rg_in
  rpi=rpi_in
  rcpd=rcpd_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


  ! 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=16.6/2./sqrt(2.) ! default value as in yamada4
  CALL getin_p('atke_cepsilon',cepsilon)

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

END SUBROUTINE atke_ini

END MODULE  atke_turbulence_ini_mod