MODULE comsoil_h_PEM

implicit none

integer, parameter                  :: nsoilmx_PEM = 68   ! number of layers in the PEM
real, allocatable, dimension(:)     :: layer_PEM          ! soil layer depths [m]
real, allocatable, dimension(:)     :: mlayer_PEM         ! soil mid-layer depths [m]
real, allocatable, dimension(:,:,:) :: TI_PEM             ! soil thermal inertia [SI]
real, allocatable, dimension(:,:)   :: inertiedat_PEM     ! soil thermal inertia saved as reference for current climate [SI]
! Variables (FC: built in firstcall in soil.F)
real, allocatable, dimension(:,:,:) :: tsoil_PEM          ! sub-surface temperatures [K]
real, allocatable, dimension(:,:)   :: mthermdiff_PEM     ! (FC) mid-layer thermal diffusivity [SI]
real, allocatable, dimension(:,:)   :: thermdiff_PEM      ! (FC) inter-layer thermal diffusivity [SI]
real, allocatable, dimension(:)     :: coefq_PEM          ! (FC) q_{k+1/2} coefficients [SI]
real, allocatable, dimension(:,:)   :: coefd_PEM          ! (FC) d_k coefficients [SI]
real, allocatable, dimension(:,:)   :: alph_PEM           ! (FC) alpha_k coefficients [SI]
real, allocatable, dimension(:,:)   :: beta_PEM           ! beta_k coefficients [SI]
real                                :: mu_PEM             ! mu coefficient [SI]
real                                :: fluxgeo            ! Geothermal flux [W/m^2]
real                                :: depth_breccia      ! Depth at which we have breccia [m]
real                                :: depth_bedrock      ! Depth at which we have bedrock [m]
integer                             :: index_breccia      ! last index of the depth grid before having breccia
integer                             :: index_bedrock      ! last index of the depth grid before having bedrock
logical                             :: soil_pem           ! True by default, to run with the subsurface physic. Read in pem.def
real                                :: ini_huge_h2oice    ! Initial value for huge reservoirs of H2O ice [kg/m^2]
logical                             :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure

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

SUBROUTINE ini_comsoil_h_PEM(ngrid,nslope)

implicit none

integer, intent(in) :: ngrid  ! number of atmospheric columns
integer, intent(in) :: nslope ! number of slope within a mesh

allocate(layer_PEM(nsoilmx_PEM))
allocate(mlayer_PEM(0:nsoilmx_PEM - 1))
allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM - 1))
allocate(thermdiff_PEM(ngrid,nsoilmx_PEM - 1))
allocate(coefq_PEM(0:nsoilmx_PEM - 1))
allocate(coefd_PEM(ngrid,nsoilmx_PEM - 1))
allocate(alph_PEM(ngrid,nsoilmx_PEM - 1))
allocate(beta_PEM(ngrid,nsoilmx_PEM - 1))
allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))

END SUBROUTINE ini_comsoil_h_PEM

!=======================================================================

SUBROUTINE end_comsoil_h_PEM

implicit none

if (allocated(layer_PEM)) deallocate(layer_PEM)
if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
if (allocated(TI_PEM)) deallocate(TI_PEM)
if (allocated(tsoil_PEM)) deallocate(tsoil_PEM)
if (allocated(mthermdiff_PEM)) deallocate(mthermdiff_PEM)
if (allocated(thermdiff_PEM)) deallocate(thermdiff_PEM)
if (allocated(coefq_PEM)) deallocate(coefq_PEM)
if (allocated(coefd_PEM)) deallocate(coefd_PEM)
if (allocated(alph_PEM)) deallocate(alph_PEM)
if (allocated(beta_PEM)) deallocate(beta_PEM)
if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)

END SUBROUTINE end_comsoil_h_PEM

END MODULE comsoil_h_PEM