MODULE earth_const
   USE prec, ONLY: rstd
   USE thermodynamics_mod, ONLY: thermodynamics_t
   IMPLICIT NONE
   SAVE

! init_earth_const is called from icosa_init outside of an OMP PARALLEL section
! hence the variables set here are not THREADPRIVATE

   REAL(rstd) :: radius = 6.37122E6
   REAL(rstd) :: daysec = 86400 ! day length (s)
   REAL(rstd) :: omega = 7.292E-5
   REAL(rstd) :: scale_height = 8000. ! atmospheric scale height (m)
   REAL(rstd) :: scale_factor = 1.
   REAL(rstd) :: pa = 50000.   ! default value set to preff/2 by disvert_std
   REAL(rstd), BIND(C) :: g = 9.80616    ! acceleration of gravity assumed constant

   TYPE(thermodynamics_t) :: thermo ! should eventually replace kappa, Treff, preff, cpp, cppv, Rd, Rv

   REAL(rstd) :: kappa = 0.2857143
   REAL(rstd), BIND(C) :: Treff = 273., &     ! reference temperature, used in definition of entropy
                          preff = 101325., &  ! reference pressure, used in definition of entropy and theta
                          cpp = 1004.70885, & ! specific heat at constant pressure at t=Treff for dry air
                          cppv = 1860., &     ! specific heat at constant pressure at t=Treff for water vapor
                          Rd, &               ! specific perfect gas constant for dry air
                          Rv = 461.5, &       ! specific perfect gas constant for water vapor
                          nu = 0.35           ! exponent in variable-Cp law Cp=cpp*(T/Treff)^nu

   INTEGER, PARAMETER :: thermo_none = -99, thermo_theta = 1, thermo_entropy = 2, thermo_variable_Cp = 3, &
                         thermo_moist = 4, thermo_boussinesq = 5, thermo_dry = 10, thermo_fake_moist = 11, thermo_moist_debug = 100

   ! variables below are set in init_caldyn_gcm => need THREADPRIVATE

   INTEGER, BIND(C) :: caldyn_thermo
   !$OMP THREADPRIVATE(caldyn_thermo)

   INTEGER, BIND(C) :: physics_thermo
   !$OMP THREADPRIVATE(physics_thermo)

   LOGICAL :: boussinesq
   !$OMP THREADPRIVATE(boussinesq)

CONTAINS

   SUBROUTINE init_earth_const
      USE getin_mod, ONLY: getin

      ! @getin_doc g realm=planet type=real default=9.80616
      ! Acceleration of gravity, assumed constant within the traditional
      ! shallow-atmosphere approximation, in m^2/s.
      CALL getin("g", g)

      ! @getin_doc radius realm=planet type=real default=6371220
      ! Planetary radius, assumed constant within the traditional
      ! shallow-atmosphere approximation, in m.
      CALL getin("radius", radius)

      ! @getin_doc omega realm=planet type=real default=7.292e-5
      ! Planetary rotation rate, in radian/second.
      CALL getin("omega", omega)

      ! @getin_doc scale_factor realm=planet type=real default=1
      ! Scaling factor for DCMIP 2012 small-planet experiments.
      ! When set, the planetary radius is `radius/scale_factor`
      ! and the planetary rotation rate is `omega*scale_factor`.
      ! To also rescale time step `dt` and other time scales,
      ! see `time_style`.
      CALL getin("scale_factor", scale_factor)
      radius = radius/scale_factor
      omega = omega*scale_factor

      ! @getin_doc scale_height realm=planet type=real default=8000
      ! Atmospheric scale height in meters. For an isothermal atmosphere
      ! and a perfect gas, the scale height is equal to `RT/g`.
      ! This length scale is used by certain choices of the hybrid
      ! coefficients for the vertical pressure-based coordinate, see `disvert`.
      CALL getin("scale_height", scale_height)

      ! this is redundant with the thermodynamics module
      ! we want to get rid of this when all thermodynamics computation go through that module
      CALL getin("kappa", kappa)
      CALL getin("cpp", cpp)
      CALL getin("cppv", cppv)
      CALL getin("Rv", Rv)
      CALL getin("preff", preff)
      CALL getin("Treff", Treff)
      Rd = kappa*cpp ! kappa = Rd/Cp

      ! check for obsolete option ; global 'boussinesq' is set by init_caldyn_gcm
      boussinesq = .FALSE.
      CALL getin("boussinesq", boussinesq)
      IF (boussinesq) THEN
         PRINT *, "Note : key boussinesq is deprecated, use thermo=boussinesq instead"
         STOP
      END IF
   END SUBROUTINE init_earth_const

END MODULE earth_const
