source: trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90 @ 3136

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
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_GCM - 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_GCM-1)) then
        index_powerlaw = iloop
        exit
    endif
enddo

do iloop = nsoil_GCM, nsoil_PEM-1 
    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop-nsoil_GCM)-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_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