MODULE thermcell_mod
      
      IMPLICIT NONE
      
      
! Flags for computations
                                                                  !  default
      INTEGER,PARAMETER :: dqimpl                     = 1         !  1        flag for thermcell_dq version (1 : implicit scheme || 0 : explicit scheme)
      
! Physical parameters
      
      REAL,PARAMETER :: r_aspect_thermals             = 1.0       !           Aspect ratio of the thermals (width / height)
      REAL,PARAMETER :: tau_thermals                  = 0.        !  0.       Relaxation time
      REAL,PARAMETER :: betalpha                      = 1.0       !  0.9       - factor between 0 (e=d) and 1 (rho*fraca=cst)
      REAL,PARAMETER :: afact                         = 1.        !  2./3.    Buoyancy contribution - factor between 0 and 1
      REAL,PARAMETER :: fact_epsilon                  = 1.e-4     !  2.e-3    Friction at plume borders - exponential decrease
      REAL,PARAMETER :: nu                            = 0.000     !           Geometrical contributions to entrainment and detrainment
      REAL,PARAMETER :: alpha_max                     = 0.7       !  0.7      Maximal permitted updraft fraction
      REAL,PARAMETER :: fomass_max                    = 0.5       !  0.5      Maximal permitted outgoing layer mass fraction
      REAL,PARAMETER :: pres_limit                    = 2.e5      !  1.e5     Minimal permitted pressure to trigger a thermal plume
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : linf is used to set the lowest possible first level because we allow it
!      to begin higher than the surface. It is set to 2 in order to remove the
!      first layer for gas giant.
!      If there is a surface, it has to be set to 1.
!      If someone want to call more than once the thermal plume model in some
!      grid points, this variable may become a saved array of INTEGER with size
!      ngrid. 
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      INTEGER,PARAMETER :: linf                       = 1         !     1
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : d_temp is an artificial virtual potential temperature offset added in
!      layer linf which can be used to force convection to begin in it.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      REAL,PARAMETER :: d_temp                        = 0.        !     0.
      
! Physical constants
      
      REAL,SAVE :: RTT
      REAL,SAVE :: RG
      REAL,SAVE :: RKAPPA
      REAL,SAVE :: RPI
      REAL,SAVE :: RD
      
!$OMP THREADPRIVATE(RTT, RG, RKAPPA, RPI, RD)
      
      
      CONTAINS
      
      SUBROUTINE init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV)
         
         IMPLICIT NONE
         
         REAL g
         REAL rcp
         REAL r
         REAL pi
         REAL T_h2o_ice_liq
         REAL RV
         
         RTT = T_h2o_ice_liq
         RG = g
         RKAPPA = rcp
         RPI = pi
         RD = r
         
         RETURN
      END SUBROUTINE init_thermcell_mod
      
END MODULE thermcell_mod