MODULE soil_settings_PEM_mod

implicit none

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

SUBROUTINE soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PCM,TI_PEM)

use comsoil_h_PEM, only: layer_PEM, mlayer_pem, depth_breccia, depth_bedrock, index_breccia, index_bedrock
use iostart,       only: inquire_field_ndims, get_var, get_field, inquire_field, inquire_dimension_length

implicit none

!=======================================================================
!  Author: LL, based on work by Ehouarn Millour (07/2006)
!
!  Purpose: Read and/or initialise soil depths and properties
!
! Modifications: Aug.2010 EM: use NetCDF90 to load variables (enables using
!                              r4 or r8 restarts independently of having compiled
!                              the GCM in r4 or r8)
!                June 2013 TN: Possibility to read files with a time axis
!
!  The various actions and variable read/initialized are:
!  1. Read/build layer (and midlayer) depth
!  2. Interpolate thermal inertia and temperature on the new grid, if necessary
!=======================================================================

!=======================================================================
!  arguments
!  ---------
!  inputs:
integer,                                 intent(in) :: ngrid     ! # of horizontal grid points
integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
integer,                                 intent(in) :: nsoil_PCM ! # of soil layers in the PCM
real, dimension(ngrid,nsoil_PCM,nslope), intent(in) :: TI_PCM    ! Thermal inertia  in the PCM [SI]

real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Thermal inertia   in the PEM [SI]
!=======================================================================
! local variables:
integer :: ig, iloop, islope ! loop counters
real    :: alpha, lay1       ! coefficients for building layers
real    :: index_powerlaw    ! coefficient to match the power low grid with the exponential one
!=======================================================================
! 1. Depth coordinate
! -------------------
! mlayer_PEM distribution: For the grid in common with the PCM: lay1*(1+k^(2.9)*(1-exp(-k/20))
! Then we use a power low : lay1*alpha**(k-1/2)
lay1 = 2.e-4
alpha = 2
do iloop = 0,nsoil_PCM - 1
    mlayer_PEM(iloop) = lay1*(1 + iloop**2.9*(1 - exp(-real(iloop)/20.)))
enddo

do iloop = 0,nsoil_PEM - 1
    if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then
        index_powerlaw = iloop
        exit
    endif
enddo

do iloop = nsoil_PCM,nsoil_PEM - 1
    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5))
enddo

! Build layer_PEM()
do iloop = 1,nsoil_PEM - 1
layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2.
enddo
layer_PEM(nsoil_PEM) = 2*mlayer_PEM(nsoil_PEM - 1) - mlayer_PEM(nsoil_PEM - 2)


! 2. Thermal inertia (note: it is declared in comsoil_h)
! ------------------
do ig = 1,ngrid
    do islope = 1,nslope
        do iloop = 1,nsoil_PCM
            TI_PEM(ig,iloop,islope) = TI_PCM(ig,iloop,islope)
        enddo
        if (nsoil_PEM > nsoil_PCM) then
            do iloop = nsoil_PCM + 1,nsoil_PEM
                TI_PEM(ig,iloop,islope) = TI_PCM(ig,nsoil_PCM,islope)
            enddo
        endif
    enddo
enddo

! 3. Index for subsurface layering
! ------------------
index_breccia = 1
do iloop = 1,nsoil_PEM - 1
    if (depth_breccia >= layer_PEM(iloop)) then
        index_breccia = iloop
    else
        exit
    endif
enddo

index_bedrock = 1
do iloop = 1,nsoil_PEM - 1
    if (depth_bedrock >= layer_PEM(iloop)) then
        index_bedrock = iloop
    else
        exit
    endif
enddo

END SUBROUTINE soil_settings_PEM

END MODULE soil_settings_PEM_mod