MODULE earth_const
  USE prec
  USE math_const
  
  REAL(rstd),SAVE :: radius=6.37122E6
  REAL(rstd),SAVE :: g=9.80616
  REAL(rstd),PARAMETER :: daysec=86400
  REAL(rstd),SAVE :: omega=7.292E-5
  REAL(rstd),SAVE :: kappa=0.2857143
  REAL(rstd),SAVE :: cpp=1004.70885
  REAL(rstd),SAVE :: cppv=1860.
  REAL(rstd),SAVE :: Rv=461.5
  REAL(rstd),SAVE :: Treff=273.
  REAL(rstd),SAVE :: preff=101325.
  REAL(rstd),SAVE :: pa=50000. ! default value set to preff/2 by disvert_std
  REAL(rstd),SAVE :: scale_height=8000. ! atmospheric scale height (m)
  REAL(rstd),SAVE :: scale_factor=1.
  REAL(rstd),SAVE :: gas_constant = 8.3144621 
  REAL(rstd),SAVE :: mu                 ! molar mass of the atmosphere

  INTEGER, PARAMETER,PUBLIC :: thermo_none=-99, thermo_theta=1, thermo_entropy=2, &
       thermo_moist=3, thermo_boussinesq=4, thermo_dry=10, thermo_fake_moist=11, thermo_moist_debug=100, &
       caldyn_vert_noncons=1, caldyn_vert_cons=2
  INTEGER, PUBLIC :: caldyn_thermo, caldyn_vert_variant, physics_thermo
  !$OMP THREADPRIVATE(caldyn_thermo) 
  !$OMP THREADPRIVATE(physics_thermo) 

  LOGICAL, SAVE :: boussinesq, hydrostatic
  !$OMP THREADPRIVATE(boussinesq)
  !$OMP THREADPRIVATE(hydrostatic) 
  LOGICAL :: dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro, dysl_caldyn_vert
  !$OMP THREADPRIVATE(dysl, dysl_geopot, dysl_pvort_only, dysl_caldyn_fast, dysl_caldyn_coriolis, dysl_slow_hydro, dysl_caldyn_vert)

CONTAINS
  
  SUBROUTINE init_earth_const
  USE getin_mod
  IMPLICIT NONE
  REAL(rstd) :: X=1
  
    CALL getin("radius",radius)
    CALL getin("g",g)
    CALL getin("scale_factor",scale_factor)
    CALL getin("omega",omega)  
    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)
    
    boussinesq=.FALSE.
    CALL getin("boussinesq",boussinesq)  
    PRINT *, "Note : key boussinesq is deprecated, use thermo=boussinesq instead"
    IF(boussinesq) STOP
    hydrostatic=.TRUE.
    CALL getin("hydrostatic",hydrostatic)  
    IF(boussinesq .AND. .NOT. hydrostatic) THEN
       PRINT *, 'boussinesq=.TRUE. and hydrostatic=.FALSE. : Non-hydrostatic boussinesq equations are not supported'
       STOP
    END IF
    
    mu=kappa/cpp
    radius=radius/scale_factor
    omega=omega*scale_factor
    
  END SUBROUTINE init_earth_const
  
  
END MODULE earth_const