module param
      implicit none
      ! acceleration due to gravity [m/s^2]
      real, parameter :: g = 9.81
      ! molecular Prandtl number (air) [n/d]
      real, parameter :: Pr_m = 0.71
      ! kinematic viscosity of air [m^2/s]
      real, parameter :: nu_air = 0.000015e0
      ! von Karman constant [n/d]
      real, parameter :: kappa = 0.40
      ! inverse Prandtl (turbulent) number in neutral conditions [n/d]
      real, parameter :: Pr_t_0_inv = 1.15
      ! inverse Prandtl (turbulent) number in free convection [n/d]
      real, parameter :: Pr_t_inf_inv = 3.5

      ! stability function coeff. (unstable) [= g4 & g10 in deprecated code]
      real, parameter :: alpha_m = 16.0
      real, parameter :: alpha_h = 16.0
      real, parameter :: alpha_h_fix = 1.2

      ! stability function coeff. (stable)
      real, parameter :: beta_m = 4.7
      real, parameter :: beta_h = beta_m

      ! --- max Ri-bulk value in stable case ( < 1 / beta_m )
      real, parameter ::  Rib_max = 0.9 / beta_m

      ! Re fully roughness minimum value [n/d]
      real, parameter ::  Re_rough_min = 16.3

      ! roughness model coeff. [n/d]
      ! --- transitional mode
      !     B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
      real, parameter :: B1_rough = 5.0 / 6.0
      real, parameter :: B2_rough = 0.45
      real, parameter :: B3_rough = kappa * Pr_m
      ! --- fully rough mode (Re > Re_rough_min)
      !     B = B4 * Re^(B2)
      real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)

      real, parameter :: B_max_land = 2.0
      real, parameter :: B_max_ocean = 8.0

      ! Charnock parameters
      !     z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
      real, parameter :: gamma_c = 0.0144
      real, parameter :: Re_visc_min = 0.111

      real, parameter :: h_charnock = 10.0
      real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c))
      real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock

!      AN5=(A6/A0)**4
!      D1=(2.0E0*G10-AN5*G4-SQRT((AN5*G4)**2+4.0E0*AN5*G10*(G10-G4)))/(2.0E0*G10**2)
!      Y10=(1.0E0-G4*D1)**.25E0
!      X10=(1.0E0-G10*D1)**.5E0
!      P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
!      P0=ALOG((X10-1.0E0)/(X10+1.0E0))

end module PARAM