MODULE earth_const
  USE prec, ONLY : rstd
  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

! Some variable are BIND(C) to be accessible from Python

  REAL(rstd), BIND(C) :: g=9.80616    ! acceleration of gravity assumed constant
  REAL(rstd):: radius=6.37122E6
  REAL(rstd) :: omega=7.292E-5
  REAL(rstd) :: scale_factor=1.
  REAL(rstd),PARAMETER :: daysec=86400

  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

  REAL(rstd) :: nu=0.35     ! exponent in variable-Cp law Cp=cpp*(T/Treff)^nu
  REAL(rstd) :: pa=50000.   ! default value set to preff/2 by disvert_std
  REAL(rstd) :: scale_height=8000. ! atmospheric scale height (m)
  REAL(rstd) :: gas_constant = 8.3144621 
  REAL(rstd) :: mu             ! molar mass (?)

  REAL(rstd) :: kappa=0.2857143

  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
  LOGICAL :: boussinesq
  INTEGER, BIND(C) :: caldyn_thermo, physics_thermo

CONTAINS
  
  SUBROUTINE init_earth_const
    USE getin_mod, ONLY : getin
    USE grid_param, ONLY : nqdyn
    USE omp_para, ONLY : is_master
    USE mpipara, ONLY : is_mpi_root
    CHARACTER(len=255) :: def

    CALL getin("g",g)
    CALL getin("radius",radius)
    CALL getin("omega",omega)  
    CALL getin("scale_factor",scale_factor)
    radius=radius/scale_factor
    omega=omega*scale_factor

    CALL getin("kappa",kappa)
    CALL getin("cpp",cpp)
    CALL getin("cppv",cppv)
    CALL getin("Rv",Rv)
    CALL getin("preff",preff)
    CALL getin("Treff",Treff)
    CALL getin("scale_height",scale_height)
    Rd = kappa*cpp ! kappa = Rd/Cp
    mu = kappa/cpp ! FIXME ??

    nqdyn=1
    physics_thermo = thermo_none

    def='theta'
    CALL getin('thermo',def)
    SELECT CASE(TRIM(def))
    CASE('boussinesq')
       boussinesq=.TRUE.
       caldyn_thermo=thermo_boussinesq
    CASE('theta')
       caldyn_thermo=thermo_theta
       physics_thermo=thermo_dry
    CASE('entropy')
       caldyn_thermo=thermo_entropy
       physics_thermo=thermo_dry
    CASE('variable_Cp')
       caldyn_thermo=thermo_variable_Cp
       physics_thermo=thermo_variable_Cp
       CALL getin("nu",nu)
    CASE('theta_fake_moist')
       caldyn_thermo=thermo_theta
       physics_thermo=thermo_fake_moist
    CASE('entropy_fake_moist')
       caldyn_thermo=thermo_entropy
       physics_thermo=thermo_fake_moist
    CASE('moist')
       caldyn_thermo=thermo_moist_debug
       physics_thermo=thermo_moist
       nqdyn = 2
    CASE DEFAULT
       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_thermo : <', &
            TRIM(def),'> options are <theta>, <entropy>'
       STOP
    END SELECT

    IF(is_master) THEN
       SELECT CASE(caldyn_thermo)
       CASE(thermo_theta)
          PRINT *, 'caldyn_thermo = thermo_theta'
       CASE(thermo_entropy)
          PRINT *, 'caldyn_thermo = thermo_entropy'
       CASE(thermo_moist_debug)
          PRINT *, 'caldyn_thermo = thermo_moist_debug'
       CASE DEFAULT
          STOP
       END SELECT

       SELECT CASE(physics_thermo)
       CASE(thermo_dry)
          PRINT *, 'physics_thermo = thermo_dry'
       CASE(thermo_fake_moist)
          PRINT *, 'physics_thermo = thermo_fake_moist'
       CASE(thermo_moist)
          PRINT *, 'physics_thermo = thermo_moist'
       END SELECT

    END IF

  END SUBROUTINE init_earth_const

END MODULE earth_const