MODULE soil_settings_PEM_mod

implicit none

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

SUBROUTINE soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_GCM,TI_PEM)

!      use netcdf
use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, depth_breccia, depth_bedrock, index_breccia, index_bedrock
use comsoil_h,     only: inertiedat, layer, mlayer, volcapa

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_GCM ! # of soil layers in the GCM
real, dimension(ngrid,nsoil_GCM,nslope), intent(in) :: TI_GCM    ! Thermal inertia  in the GCM [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
logical :: found
real    :: alpha,lay1 ! coefficients for building layers
real    :: xmin, xmax ! to display min and max of a field
!=======================================================================
! 1. Depth coordinate
! -------------------
! 1.4 Build mlayer(), if necessary
! default mlayer distribution, following a power law:
!  mlayer(k)=lay1*alpha**(k-1/2)
lay1 = 2.e-4
alpha = 2
do iloop = 0,nsoil_PEM - 1
    mlayer_PEM(iloop) = lay1*(alpha**(iloop - 0.5))
enddo

! 1.5 Build layer(); following the same law as mlayer()
! Assuming layer distribution follows mid-layer law:
! layer(k)=lay1*alpha**(k-1)
lay1 = sqrt(mlayer_PEM(0)*mlayer_PEM(1))
alpha = mlayer_PEM(1)/mlayer_PEM(0)
do iloop = 1,nsoil_PEM
    layer_PEM(iloop) = lay1*(alpha**(iloop - 1))
enddo

! 2. Thermal inertia (note: it is declared in comsoil_h)
! ------------------
do ig = 1,ngrid
    do islope = 1,nslope
        do iloop = 1,nsoil_GCM
            TI_PEM(ig,iloop,islope) = TI_GCM(ig,iloop,islope)
        enddo
        if (nsoil_PEM > nsoil_GCM) then
            do iloop = nsoil_GCM + 1,nsoil_PEM
                TI_PEM(ig,iloop,islope) = TI_GCM(ig,nsoil_GCM,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