source: trunk/LMDZ.COMMON/libf/evolution/soil.F90 @ 4175

MODULE soil
!-----------------------------------------------------------------------
! NAME
!     soil
!
! DESCRIPTION
!     Soil state and properties for PEM.
!
! AUTHORS & DATE
!     L. Lange, 04/2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics, only: dp, di, k4

! DECLARATION
! -----------
implicit none

! PARAMETERS
! ----------
real(dp), parameter :: rho_regolith      = 2000._dp ! Density of the martian regolith [kg/m^3], Zent et al. 1995; Buhler and Piqueux 2021
real(dp), parameter :: regolith_porosity = 0.4_dp   ! Porosity of the martian regolith for a random loose packing of monodisperse sphere, Scott 1960
real(dp), parameter :: breccia_porosity  = 0.1_dp   ! Porosity of breccia

logical(k4),                             protected          :: do_soil            ! To run with subsurface physics
logical(k4),                             protected          :: reg_thprop_dependp ! Regolith thermal properties depend on surface pressure
real(dp),                                protected          :: flux_geo           ! Geothermal flux [W/m^2]
real(dp),                                protected          :: depth_breccia      ! Depth of breccia [m]
real(dp),                                protected          :: depth_bedrock      ! Depth of bedrock [m]
real(dp), allocatable, dimension(:,:,:), protected, private :: TI_PCM             ! Soil thermal inertia in the "startfi.nc" at the beginning [SI]

! VARIABLES
! ---------
real(dp), allocatable, dimension(:)     :: layer              ! Soil layer depths [m]
real(dp), allocatable, dimension(:)     :: mlayer             ! Soil mid-layer depths [m]
real(dp), allocatable, dimension(:,:,:) :: TI                 ! Soil thermal inertia [SI]
real(dp), allocatable, dimension(:,:)   :: inertiedat         ! Soil thermal inertia saved as reference for current climate [SI]
real(dp), allocatable, dimension(:,:)   :: inertiedat_PCM     ! Soil thermal inertia saved as reference for current climate in the PCM [SI]
real(dp), allocatable, dimension(:,:)   :: mthermdiff         ! Mid-layer thermal diffusivity [SI]
real(dp), allocatable, dimension(:,:)   :: thermdiff          ! Inter-layer thermal diffusivity [SI]
real(dp), allocatable, dimension(:)     :: coefq              ! q_{k+1/2} coefficients [SI]
real(dp), allocatable, dimension(:,:)   :: coefd              ! d_k coefficients [SI]
real(dp), allocatable, dimension(:,:)   :: alph               ! alpha_k coefficients [SI]
real(dp), allocatable, dimension(:,:)   :: beta               ! beta_k coefficients [SI]
real(dp)                                :: mu                 ! mu coefficient [SI]
integer(di)                             :: index_breccia      ! Last index of the depth grid before having breccia
integer(di)                             :: index_bedrock      ! Last index of the depth grid before having bedrock
real(dp)                                :: volcapa            ! Soil volumetric heat capacity

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE ini_soil()
!-----------------------------------------------------------------------
! NAME
!     ini_soil
!
! DESCRIPTION
!     Allocate soil arrays based on grid dimensions.
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil, nsoil_PCM

! DECLARATION
! -----------
implicit none

! CODE
! ----
allocate(layer(nsoil))
allocate(mlayer(0:nsoil - 1))
allocate(TI(ngrid,nsoil,nslope))
allocate(mthermdiff(ngrid,0:nsoil - 1))
allocate(thermdiff(ngrid,nsoil - 1))
allocate(coefq(0:nsoil - 1))
allocate(coefd(ngrid,nsoil - 1))
allocate(alph(ngrid,nsoil - 1))
allocate(beta(ngrid,nsoil - 1))
allocate(inertiedat(ngrid,nsoil))
allocate(inertiedat_PCM(ngrid,nsoil_PCM))
allocate(TI_PCM(ngrid,nsoil_PCM,nslope))
inertiedat_PCM(:,:) = 0._dp
TI_PCM(:,:,:) = 0._dp

END SUBROUTINE ini_soil
!=======================================================================

!=======================================================================
SUBROUTINE end_soil
!-----------------------------------------------------------------------
! NAME
!     end_soil
!
! DESCRIPTION
!     Deallocate all soil arrays.
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! CODE
! ----
if (allocated(layer)) deallocate(layer)
if (allocated(mlayer)) deallocate(mlayer)
if (allocated(TI)) deallocate(TI)
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(inertiedat)) deallocate(inertiedat)
if (allocated(inertiedat_PCM)) deallocate(inertiedat_PCM)
if (allocated(TI_PCM)) deallocate(TI_PCM)

END SUBROUTINE end_soil
!=======================================================================

!=======================================================================
SUBROUTINE set_soil_config(do_soil_in,reg_thprop_dependp_in,flux_geo_in,depth_breccia_in,depth_bedrock_in)
!-----------------------------------------------------------------------
! NAME
!     set_soil_config
!
! DESCRIPTION
!     Setter for "soil" configuration parameters.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility,  only: real2str, bool2str
use display,  only: print_msg, LVL_NFO
use stoppage, only: stop_clean

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
logical(k4), intent(in) :: do_soil_in, reg_thprop_dependp_in
real(dp),    intent(in) :: flux_geo_in, depth_breccia_in, depth_bedrock_in

! CODE
! ----
do_soil = do_soil_in
reg_thprop_dependp = reg_thprop_dependp_in
flux_geo = flux_geo_in
depth_breccia = depth_breccia_in
depth_bedrock = depth_bedrock_in
call print_msg('do_soil            = '//bool2str(do_soil),LVL_NFO)
call print_msg('reg_thprop_dependp = '//bool2str(reg_thprop_dependp),LVL_NFO)
call print_msg('flux_geo           = '//real2str(flux_geo),LVL_NFO)
call print_msg('depth_breccia      = '//real2str(depth_breccia),LVL_NFO)
call print_msg('depth_bedrock      = '//real2str(depth_bedrock),LVL_NFO)
if (depth_breccia < 0._dp) call stop_clean(__FILE__,__LINE__,'''depth_breccia'' must be positive!',1)
if (depth_bedrock < 0._dp) call stop_clean(__FILE__,__LINE__,'''depth_bedrock'' must be positive!',1)
if (depth_bedrock <= depth_breccia) call stop_clean(__FILE__,__LINE__,'''depth_bedrock'' must be greater than ''depth_breccia''!',1)

END SUBROUTINE set_soil_config
!=======================================================================

!=======================================================================
SUBROUTINE set_TI_PCM(TI_PCM_in)
!-----------------------------------------------------------------------
! NAME
!     set_TI_PCM
!
! DESCRIPTION
!     Setter for 'TI_PCM'.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:,:,:), intent(in) :: TI_PCM_in

! CODE
! ----
TI_PCM(:,:,:) = TI_PCM_in(:,:,:)

END SUBROUTINE set_TI_PCM
!=======================================================================

!=======================================================================
SUBROUTINE set_inertiedat_PCM(inertiedat_PCM_in)
!-----------------------------------------------------------------------
! NAME
!     set_inertiedat_PCM
!
! DESCRIPTION
!     Setter for 'inertiedat_PCM'.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in) :: inertiedat_PCM_in

! CODE
! ----
inertiedat_PCM(:,:) = inertiedat_PCM_in(:,:)

END SUBROUTINE set_inertiedat_PCM
!=======================================================================

!=======================================================================
SUBROUTINE set_soil(therm_initia)
!-----------------------------------------------------------------------
! NAME
!     set_soil
!
! DESCRIPTION
!     Initialize soil depths and properties. Builds layer depths using
!     exponential then power-law distribution; interpolates thermal inertia
!     from PCM to PEM grid.
!
! AUTHORS & DATE
!     E. Millour, 07/2006
!     L. Lange, 2023
!     JB Clement, 2023-2025
!
! NOTES
!     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
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil, nsoil_PCM

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:,:,:), intent(inout) :: therm_initia ! Thermal inertia in the PEM [SI]

! LOCAL VARIABLES
! ---------------
integer(di) :: ig, iloop, islope ! loop counters
real(dp)    :: alpha, lay1       ! coefficients for building layers
real(dp)    :: index_powerlaw    ! coefficient to match the power low grid with the exponential one

! CODE
! ----
! 1. Depth coordinate
! -------------------
! mlayer 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(iloop) = lay1*(1 + iloop**2.9*(1._dp - exp(-real(iloop,dp)/20._dp)))
end do

do iloop = 0,nsoil - 1
    if (lay1*(alpha**(iloop - 0.5_dp)) > mlayer(nsoil_PCM - 1)) then
        index_powerlaw = iloop
        exit
    end if
end do

do iloop = nsoil_PCM,nsoil - 1
    mlayer(iloop) = lay1*(alpha**(index_powerlaw + (iloop - nsoil_PCM) - 0.5_dp))
end do

! Build layer()
do iloop = 1,nsoil - 1
    layer(iloop) = (mlayer(iloop) + mlayer(iloop - 1))/2._dp
end do
layer(nsoil) = 2._dp*mlayer(nsoil - 1) - mlayer(nsoil - 2)


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

! 3. Index for subsurface layering
! ------------------
index_breccia = 1
do iloop = 1,nsoil - 1
    if (depth_breccia >= layer(iloop)) then
        index_breccia = iloop
    else
        exit
    end if
end do

index_bedrock = 1
do iloop = 1,nsoil - 1
    if (depth_bedrock >= layer(iloop)) then
        index_bedrock = iloop
    else
        exit
    end if
end do

END SUBROUTINE set_soil
!=======================================================================

!=======================================================================
SUBROUTINE build4PCM_soil(tsoil_avg,tsoil_dev,inertiesoil4PCM,tsoil4PCM,flux_geo4PCM)
!-----------------------------------------------------------------------
! NAME
!     build4PCM_soil
!
! DESCRIPTION
!     Reconstructs the soil temperature profile, thermal inertia and
!     the geothermal flux for the PCM.
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!     C. Metz, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: nsoil_PCM
use display,  only: print_msg, LVL_NFO

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:,:,:), intent(in)  :: tsoil_avg, tsoil_dev
real(dp), dimension(:,:,:), intent(out) :: tsoil4PCM, inertiesoil4PCM
real(dp), dimension(:,:),   intent(out) :: flux_geo4PCM

! CODE
! ----
call print_msg('> Building soil thermal inertia for the PCM',LVL_NFO)
inertiesoil4PCM(:,:,:) = TI(:,1:nsoil_PCM,:)

call print_msg('> Building soil temperature profile for the PCM',LVL_NFO)
tsoil4PCM(:,:,:) = tsoil_avg(:,1:nsoil_PCM,:) + tsoil_dev(:,:,:)

call print_msg('> Building geothermal flux for the PCM',LVL_NFO)
flux_geo4PCM(:,:) = (tsoil_avg(:,nsoil_PCM + 1,:) - tsoil4PCM(:,nsoil_PCM,:))/(mlayer(nsoil_PCM + 1) - mlayer(nsoil_PCM))

END SUBROUTINE build4PCM_soil
!=======================================================================

END MODULE soil