source: trunk/LMDZ.COMMON/libf/evolution/soil_therm_inertia.F90

MODULE soil_therm_inertia
!-----------------------------------------------------------------------
! NAME
!     soil_therm_inertia
!
! DESCRIPTION
!     Compute and update soil thermal properties based on ice content,
!     pressure, and cementation state.
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! PARAMETERS
! ----------
real(dp), parameter :: reg_inertie_thresold = 370._dp ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
real(dp), parameter :: reg_inertie_minvalue = 50._dp  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
real(dp), parameter :: reg_inertie_maxvalue = 370._dp ! Maximum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
real(dp), parameter :: C = 0.0015_dp                  ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [unitless]
real(dp), parameter :: K = 8.1*1e4_dp                 ! Constant to derive TI as a function of P, from Presley and Christensen 1997 [Torr, or 133.3Pa]
real(dp), parameter :: Pa2Torr = 1./133.3_dp          ! Conversion from Pa to tor [Pa/Torr]

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

!=======================================================================
SUBROUTINE get_ice_TI(ispureice,pore_filling,surf_thermalinertia,ice_thermalinertia)
!-----------------------------------------------------------------------
! NAME
!     get_ice_TI
!
! DESCRIPTION
!     Compute ice thermal properties based on pore filling and purity.
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 2025
!
! NOTES
!     Uses Siegler et al. (2012) formula for pore-filling ice;
!     uses Mellon et al. (2004) value for pure water ice.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use planet, only: porosity

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

! ARGUMENTS
! ---------
logical(k4), intent(in)  :: ispureice           ! Boolean to check if ice is massive or just pore filling
real(dp),    intent(in)  :: pore_filling        ! ice pore filling in each layer (1)
real(dp),    intent(in)  :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)
real(dp),    intent(out) :: ice_thermalinertia  ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)

! LOCAL VARIABLES
! ---------------
real(dp) :: inertie_purewaterice = 2100._dp ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004)

! CODE
! ----
if (ispureice) then
    ice_thermalinertia = inertie_purewaterice
else
    ice_thermalinertia = sqrt(surf_thermalinertia**2 + porosity*pore_filling*inertie_purewaterice**2) ! Siegler et al., 2012
end if

END SUBROUTINE get_ice_TI
!=======================================================================

!=======================================================================
SUBROUTINE update_soil_TI(tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI)
!-----------------------------------------------------------------------
! NAME
!     update_soil_TI
!
! DESCRIPTION
!     Update soil thermal inertia based on ice table, regolith properties,
!     and pressure-dependent cementation.
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil
use soil,     only: volcapa, layer, inertiedat, depth_breccia, depth_bedrock, index_breccia, index_bedrock, reg_thprop_dependp
use planet,   only: TI_breccia, TI_bedrock, TI_regolith_avg, P610
use display,  only: print_msg

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

! ARGUMENTS
! ---------
real(dp),                   intent(in)    :: p_avg_new            ! Global average surface pressure [Pa]
real(dp), dimension(:,:),   intent(in)    :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
real(dp), dimension(:,:),   intent(in)    :: waterice             ! Surface Water ice [kg/m^2]
real(dp), dimension(:,:),   intent(in)    :: icetable_depth       ! Ice table depth [m]
real(dp), dimension(:,:),   intent(in)    :: icetable_thickness   ! Ice table thickness [m]
real(dp), dimension(:,:,:), intent(in)    :: ice_porefilling      ! Ice porefilling [m^3/m^3]
logical(k4),                intent(in)    :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
real(dp), dimension(:,:,:), intent(inout) :: TI                   ! Soil Thermal Inertia [J/m^2/K/s^1/2]

! LOCAL VARIABLES
! ---------------
integer(di)                       :: ig, islope, isoil, iref, iend ! Loop variables
real(dp), dimension(ngrid,nslope) :: regolith_inertia              ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
real(dp)                          :: delta                         ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
real(dp)                          :: ice_bottom_depth              ! Bottom depth of the subsurface ice [m]
real(dp)                          :: d_part                        ! Regolith particle size [micrometer]
real(dp)                          :: ice_inertia                   ! Inertia of water ice [SI]

! CODE
! ----
call print_msg("> Updating soil properties")

! 1. Modification of the regolith thermal inertia.
do islope = 1,nslope
    regolith_inertia(:,islope) = inertiedat(:,1)
    do ig = 1,ngrid
        if (tendencies_waterice(ig,islope) < -1.e-5_dp .and. abs(waterice(ig,islope)) < minieps) regolith_inertia(ig,islope) = TI_regolith_avg
        if (reg_thprop_dependp) then
            if (TI(ig,1,islope) < reg_inertie_thresold) then
                d_part = (regolith_inertia(ig,islope)**2/(volcapa*C*(P610*Pa2Torr)**0.6))**(-1./(0.11_dp*log10(P610*Pa2Torr/K))) ! compute particle size, in micrometer
                regolith_inertia(ig,islope) = sqrt(volcapa*C*(p_avg_new*Pa2Torr)**0.6*d_part**(-0.11_dp*log10(p_avg_new*Pa2Torr/K)))
                if (regolith_inertia(ig,islope) > reg_inertie_maxvalue) regolith_inertia(ig,islope) = reg_inertie_maxvalue
                if (regolith_inertia(ig,islope) < reg_inertie_minvalue) regolith_inertia(ig,islope) = reg_inertie_minvalue
            end if
        end if
    end do
end do

! 2. Build new Thermal Inertia
do islope = 1,nslope
    do ig = 1,ngrid
        do isoil = 1,index_breccia
            TI(ig,isoil,islope) = regolith_inertia(ig,islope)
        end do
        if (regolith_inertia(ig,islope) < TI_breccia) then
!!! transition
            delta = depth_breccia
            TI(ig,index_breccia + 1,islope) = sqrt((layer(index_breccia + 1) - layer(index_breccia))/                  &
                                                  (((delta - layer(index_breccia))/(TI(ig,index_breccia,islope)**2)) + &
                                                  ((layer(index_breccia + 1) - delta)/(TI_breccia**2))))
            do isoil = index_breccia + 2,index_bedrock
                TI(ig,isoil,islope) = TI_breccia
            end do
        else ! we keep the high ti values
            do isoil = index_breccia + 1,index_bedrock
                TI(ig,isoil,islope) = TI(ig,index_breccia,islope)
            end do
        end if ! TI PEM and breccia comparison
!!! transition
        delta = depth_bedrock
        TI(ig,index_bedrock + 1,islope) = sqrt((layer(index_bedrock + 1) - layer(index_bedrock))/                  &
                                              (((delta - layer(index_bedrock))/(TI(ig,index_bedrock,islope)**2)) + &
                                              ((layer(index_bedrock + 1) - delta)/(TI_bedrock**2))))
        do isoil = index_bedrock + 2,nsoil
            TI(ig,isoil,islope) = TI_bedrock
        end do
    end do ! ig
end do ! islope

! 3. Build new TI in case of ice table
do ig = 1,ngrid
    do islope = 1,nslope
        if (icetable_depth(ig,islope) > -1.e-6_dp) then
        ! 3.0 Case where it is perennial ice
            if (icetable_depth(ig,islope) < minieps) then
                call get_ice_TI(.true.,1._dp,0._dp,ice_inertia)
                do isoil = 1,nsoil
                    TI(ig,isoil,islope) = ice_inertia
                end do
            else
                if (icetable_equilibrium) then
                    call get_ice_TI(.false.,1._dp,regolith_inertia(ig,islope),ice_inertia)
                    ! 3.1.1 find the index of the mixed layer
                    iref = 0 ! initialize iref
                    do isoil = 1,nsoil ! loop on layers to find the beginning of the ice table
                        if (icetable_depth(ig,islope) >= layer(isoil)) then
                            iref = isoil ! pure regolith layer up to here
                        else
                            exit ! correct iref was obtained in previous cycle
                        end if
                    end do
                    ! 3.1.2 find the index of the end of the ice table
                    iend = 0 ! initialize iend
                    ice_bottom_depth = icetable_depth(ig,islope) + icetable_thickness(ig,islope)
                    do isoil = 1,nsoil ! loop on layers to find the end of the ice table
                        if (ice_bottom_depth >= layer(isoil)) then
                            iend = isoil ! pure regolith layer up to here
                        else
                            exit ! correct iref was obtained in previous cycle
                        end if
                    end do
                    ! 3.2 Build the new ti
                    if (iref < nsoil) then
                        if (iref == iend) then
                            ! Ice table begins and end in the same mixture with three components
                            if (iref /= 0) then ! mixed layer
                                TI(ig,iref + 1,islope) = sqrt((layer(iref + 1) - layer(iref))/                                      &
                                                             (((icetable_depth(ig,islope) - layer(iref))/(TI(ig,iref,islope)**2)) + &
                                                             ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) +    &
                                                             ((layer(iref + 1) - ice_bottom_depth)/(TI(ig,iref + 1,islope)**2))))
                            else ! first layer is already a mixed layer
                                ! (ie: take layer(iref=0)=0)
                                TI(ig,1,islope) = sqrt((layer(1))/                                                        &
                                                      (((icetable_depth(ig,islope))/(TI(ig,1,islope)**2)) +               &
                                                      ((ice_bottom_depth - icetable_depth(ig,islope))/(ice_inertia**2)) + &
                                                      ((layer(2) - ice_bottom_depth)/(TI(ig,2,islope)**2))))
                            end if ! of if (iref /= 0)
                        else
                            if (iref /= 0) then ! mixed layer
                                TI(ig,iref + 1,islope) = sqrt((layer(iref + 1) - layer(iref))/                                      &
                                                             (((icetable_depth(ig,islope) - layer(iref))/(TI(ig,iref,islope)**2)) + &
                                                             ((layer(iref + 1) - icetable_depth(ig,islope))/(ice_inertia**2))))
                            else ! first layer is already a mixed layer
                                ! (ie: take layer(iref=0)=0)
                                TI(ig,1,islope) = sqrt((layer(1))/                                          &
                                                      (((icetable_depth(ig,islope))/(TI(ig,1,islope)**2)) + &
                                                      ((layer(1) - icetable_depth(ig,islope))/(ice_inertia**2))))
                            end if ! of if (iref /= 0)
                        end if ! iref == iend

                        TI(ig,iref + 2:iend,islope) = ice_inertia
                        if (iend < nsoil) then
                            TI(ig,iend + 1,islope) = sqrt((layer(iend + 1) - layer(iend))/                             &
                                                         (((ice_bottom_depth - layer(iend))/(TI(ig,iend,islope)**2)) + &
                                                         ((layer(iend + 1) - ice_bottom_depth)/(ice_inertia**2))))
                        end if
                    end if ! of if (iref < nsoil)
                else if (icetable_dynamic) then
                    do  isoil = 1,nsoil
                        call get_ice_TI(.false.,ice_porefilling(ig,isoil,islope),regolith_inertia(ig,islope),TI(ig,isoil,islope))
                    end do
                end if ! of if icetable_equilibrium
            end if ! permanent glaciers
        end if ! icetable_depth(ig,islope) > 0.
    end do !islope
end do !ig

END SUBROUTINE update_soil_TI
!=======================================================================

END MODULE soil_therm_inertia