module comsoil_h_PEM

implicit none
  integer, parameter :: nsoilmx_PEM = 27               ! number of layers in the PEM
  real,save,allocatable,dimension(:) :: layer_PEM      ! soil layer depths   [m]
  real,save,allocatable,dimension(:) :: mlayer_PEM     ! soil mid-layer depths [m]
  real,save,allocatable,dimension(:,:,:) :: TI_PEM     ! soil thermal inertia [SI]
  real,save,allocatable,dimension(:,:) :: inertiedat_PEM ! soil thermal inertia saved as reference for current climate [SI]
  ! variables (FC: built in firstcall in soil.F)
  REAL,SAVE,ALLOCATABLE :: tsoil_PEM(:,:,:)       ! sub-surface temperatures [K]
  real,save,allocatable :: mthermdiff_PEM(:,:)  ! (FC) mid-layer thermal diffusivity [SI]
  real,save,allocatable :: thermdiff_PEM(:,:)   ! (FC) inter-layer thermal diffusivity [SI]
  real,save,allocatable :: coefq_PEM(:)         ! (FC) q_{k+1/2} coefficients [SI]
  real,save,allocatable :: coefd_PEM(:,:)       ! (FC) d_k coefficients [SI]
  real,save,allocatable :: alph_PEM(:,:)      ! (FC) alpha_k coefficients [SI]
  real,save,allocatable :: beta_PEM(:,:)      ! beta_k coefficients [SI]
  real,save :: mu_PEM                           ! mu coefficient [SI]
  real,save :: fluxgeo                               ! Geothermal flux [W/m^2]
  real,save :: depth_breccia                         ! Depth at which we have breccia [m]
  real,save :: depth_bedrock                         ! Depth at which we have bedrock [m]
  integer,save :: index_breccia                         ! last index of the depth grid before having breccia
  integer,save :: 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,save,allocatable,dimension(:) :: water_reservoir      ! Large reserve of water   [kg/m^2]
  real,save :: water_reservoir_nom
  logical, save :: 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)) 
    allocate(water_reservoir(ngrid))
  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)
    if (allocated(water_reservoir)) deallocate(water_reservoir)
  end subroutine end_comsoil_h_PEM

end module comsoil_h_PEM