module comsoil_h

implicit none
! nsoilmx : number of subterranean layers
!EM: old soil routine:      integer, parameter :: nsoilmx = 10
  integer, parameter :: nsoilmx = 18 

  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 fromn control(35) from physicq start file
       !     in physdem (or set via tabfi, or initialized in
       !                 soil_settings.F)

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

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

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

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))
 
  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)
  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))
 
  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)

  end subroutine end_comsoil_h_slope_var

end module comsoil_h