MODULE lmdz_atke_turbulence_ini

implicit none

! declaration of constants and parameters
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 :: cinffac
  !$OMP THREADPRIVATE(cinffac)
  real :: l0,ric,ri0,ri1,cinf,cn,cepsilon,pr_slope,pr_asym,pr_neut,clmix,clmixshear,smmin,ctkes,cke
  !$OMP THREADPRIVATE(l0,ri0,ri1,ric,cinf,cn,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 corresponding to the Kolmogov dissipation scale 
                                 ! in planetary atmospheres (Chen et al 2016, JGR Atmos)
  !$OMP THREADPRIVATE(lmin)
  logical :: atke_ok_vdiff, atke_ok_virtual
  !$OMP THREADPRIVATE(atke_ok_vdiff,atke_ok_virtual)

CONTAINS

SUBROUTINE atke_ini(rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in)

  USE ioipsl_getin_p_mod, ONLY : getin_p

  real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in, rv_in, viscom_in, viscoh_in


  ! input arguments (universal constants for planet)
  !-------------------------------------------------
 
  ! gravity acceleration
  rg=rg_in
  ! dry gas constant (R/M, R=perfect gas constant and M is the molar mass of the fluid)
  rd=rd_in
  ! Pi number
  rpi=rpi_in
  ! cp per unit mass of the gas
  rcpd=rcpd_in
  ! water vapor constant (for simulations in Earth's atmosphere)
  rv=rv_in
  ! kinematic molecular viscosity for momentum
  viscom=viscom_in
  ! kinematic molecular viscosity for heat
  viscoh=viscoh_in


  ! Read flag values in .def files
  !-------------------------------


  ! 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=.false.
  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)


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


  ! calculation of cn = Sm value at Ri=0
  ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
  cn=(1./sqrt(cepsilon))**(2/3)

  ! asymptotic value of Sm for Ri=-Inf
  cinffac=2.0
  CALL getin_p('atke_cinffac',cinffac)
  cinf=cinffac*cn
  if (cinf .le. cn) then
        call abort_physic("atke_turbulence_ini", &
        'cinf cannot be lower than cn', 1)
  endif


 ! 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

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


 ! asymptotic turbulent prandt number value for Ri=-Inf
  pr_asym=0.4
  CALL getin_p('atke_pr_asym',pr_asym)
  if (pr_asym .ge. pr_neut) then
        call abort_physic("atke_turbulence_ini", &
        'pr_asym must be be lower than pr_neut', 1)
  endif



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


  ! calculation of Ri0 such that continuity in slope of Sm at Ri=0
  ri0=2./rpi*(cinf - cn)*ric/cn

  ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
  ri1 = -2./rpi * (pr_asym - pr_neut)


 RETURN

END SUBROUTINE atke_ini

END MODULE  lmdz_atke_turbulence_ini