MODULE comsoil_h

implicit none

! nsoilmx : number of subterranean layers
integer, parameter :: nsoilmx = 57

real, save, allocatable, dimension(:)     :: layer       ! soil layer depths
real, save, allocatable, dimension(:)     :: mlayer      ! soil mid-layer depths
real, save, allocatable, dimension(:,:)   :: inertiedat  ! soil thermal inertia for present climate
real, save, allocatable, dimension(:,:,:) :: inertiesoil ! soil thermal inertia
real, save                                :: volcapa     ! soil volumetric heat capacity
! NB: volcapa is read from control(35) from physics start file
!     in physdem (or set via tabfi, or initialized in soil_settings.F)

!$OMP THREADPRIVATE(layer,mlayer,inertiedat,inertiesoil,volcapa)

! Variables (FC: built in firstcall in soil.F)
real, save, allocatable, dimension(:,:,:) :: tsoil      ! sub-surface temperatures (K)
real, save, allocatable, dimension(:,:,:) :: mthermdiff ! (FC) mid-layer thermal diffusivity
real, save, allocatable, dimension(:,:,:) :: thermdiff  ! (FC) inter-layer thermal diffusivity
real, save, allocatable, dimension(:)     :: coefq      ! (FC) q_{k+1/2} coefficients
real, save, allocatable, dimension(:,:,:) :: coefd      ! (FC) d_k coefficients
real, save, allocatable, dimension(:,:,:) :: alph       ! (FC) alpha_k coefficients
real, save, allocatable, dimension(:,:,:) :: beta       ! beta_k coefficients
real, save, allocatable, dimension(:,:)   :: flux_geo   ! Geothermal Flux (W/m^2)
real, save :: mu

!$OMP THREADPRIVATE(tsoil,mthermdiff,thermdiff,coefq,coefd,alph,beta,mu,flux_geo)

! Subsurface tracers:
logical, save                                  :: adsorption_soil      ! boolean to call adosrption (or not)
real,    save                                  :: choice_ads           ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
logical, save                                  :: ads_const_D          ! boolean to have constant diffusion coefficient
logical, save                                  :: ads_massive_ice      ! boolean to have massive subsurface ice
real,    save, allocatable, dimension(:,:,:,:) :: qsoil                ! subsurface tracers (kg/m^3 of regol)
integer, parameter :: nqsoil = 3                                       ! number of subsurface tracers, only three when working with water
integer, parameter :: igcm_h2o_vap_soil = 1
integer, parameter :: igcm_h2o_ice_soil = 2
integer, parameter :: igcm_h2o_vap_ads  = 3
REAL, parameter :: porosity_reg  = 0.45
!$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads)
!$OMP THREADPRIVATE(ads_const_D,ads_massive_ice)

!=======================================================================
contains
!=======================================================================

subroutine ini_comsoil_h(ngrid,nslope)

implicit none

integer,intent(in) :: ngrid  ! number of atmospheric columns
integer,intent(in) :: nslope ! number of sub grid slopes

allocate(layer(nsoilmx)) !soil layer depths
allocate(mlayer(0:nsoilmx - 1)) ! soil mid-layer depths
allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia for present climate
allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
allocate(mthermdiff(ngrid,0:nsoilmx - 1,nslope))
allocate(thermdiff(ngrid,nsoilmx - 1,nslope))
allocate(coefq(0:nsoilmx - 1))
allocate(coefd(ngrid,nsoilmx - 1,nslope))
allocate(alph(ngrid,nsoilmx - 1,nslope))
allocate(beta(ngrid,nsoilmx - 1,nslope))
allocate(flux_geo(ngrid,nslope))
allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))

END SUBROUTINE ini_comsoil_h

!=======================================================================
SUBROUTINE end_comsoil_h

implicit none

if (allocated(layer)) deallocate(layer)
if (allocated(mlayer)) deallocate(mlayer)
if (allocated(inertiedat)) deallocate(inertiedat)
if (allocated(inertiesoil)) deallocate(inertiesoil)
if (allocated(tsoil)) deallocate(tsoil)
if (allocated(mthermdiff)) deallocate(mthermdiff)
if (allocated(thermdiff)) deallocate(thermdiff)
if (allocated(coefq)) deallocate(coefq)
if (allocated(coefd)) deallocate(coefd)
if (allocated(alph)) deallocate(alph)
if (allocated(beta)) deallocate(beta)
if (allocated(flux_geo)) deallocate(flux_geo)
if (allocated(qsoil))  deallocate(qsoil)

END SUBROUTINE end_comsoil_h

!=======================================================================
SUBROUTINE ini_comsoil_h_slope_var(ngrid,nslope)

implicit none

integer,intent(in) :: ngrid  ! number of atmospheric columns
integer,intent(in) :: nslope ! number of sub grid slopes

allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
allocate(mthermdiff(ngrid,0:nsoilmx - 1,nslope))
allocate(thermdiff(ngrid,nsoilmx - 1,nslope))
allocate(coefd(ngrid,nsoilmx - 1,nslope))
allocate(alph(ngrid,nsoilmx - 1,nslope))
allocate(beta(ngrid,nsoilmx - 1,nslope))
allocate(flux_geo(ngrid,nslope))
allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))

END SUBROUTINE ini_comsoil_h_slope_var

!=======================================================================
SUBROUTINE end_comsoil_h_slope_var

implicit none

if (allocated(tsoil)) deallocate(tsoil)
if (allocated(inertiesoil)) deallocate(inertiesoil)
if (allocated(mthermdiff)) deallocate(mthermdiff)
if (allocated(thermdiff)) deallocate(thermdiff)
if (allocated(coefd)) deallocate(coefd)
if (allocated(alph)) deallocate(alph)
if (allocated(beta)) deallocate(beta)
if (allocated(flux_geo)) deallocate(flux_geo)
if (allocated(qsoil))  deallocate(qsoil)

END SUBROUTINE end_comsoil_h_slope_var

END MODULE comsoil_h