source: trunk/LMDZ.COMMON/libf/evolution/sorption.F90 @ 4065

MODULE sorption
!-----------------------------------------------------------------------
! NAME
!     sorption
!
! DESCRIPTION
!     CO2 and H2O adsorption/desorption in regolith following Zent & Quinn
!     (1995) and Jackosky et al. (1997).
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! PARAMETERS
! ----------
logical(k4), protected :: do_sorption ! Flag to compute adsorption/desorption

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

!=======================================================================
SUBROUTINE set_sorption_config(do_sorption_in)
!-----------------------------------------------------------------------
! NAME
!     set_sorption_config
!
! DESCRIPTION
!     Setter for 'sorption' configuration parameters.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use utility, only: bool2str
use display, only: print_msg

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

! ARGUMENTS
! ---------
logical(k4), intent(in) :: do_sorption_in

! CODE
! ----
do_sorption = do_sorption_in
call print_msg('do_sorption = '//bool2str(do_sorption))

END SUBROUTINE set_sorption_config
!=======================================================================

!=======================================================================
SUBROUTINE evolve_regolith_adsorption(d_h2oice,d_co2ice,h2oice,co2ice,tsoil,TI,ps,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
!-----------------------------------------------------------------------
! NAME
!     evolve_regolith_adsorption
!
! DESCRIPTION
!     Compute both H2O and CO2 adsorption in regolith.
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 12/2025
!     C. Metz, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

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

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! tendencies on the glaciers [1]
real(dp), dimension(:,:),   intent(in)    :: h2oice             ! H2O ice at the surface [kg/m^2]
real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Soil Thermal inertia [J/K/^2/s^1/2]
real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of co2 in the first layer [kg/kg]
real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2o in the first layer [kg/kg]
real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of h2o adsorbed [kg/m^3]
real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of co2 adsorbed [kg/m^3]
real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg

! LOCAL VARIABLES
! ---------------
real(dp), dimension(ngrid,nsoil,nslope) :: theta_h2o_ads ! Fraction of the pores occupied by H2O molecules

! CODE
! ----
call evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,theta_h2o_ads,h2o_ads_reg,delta_h2o_ads)
call evolve_co2_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads)

END SUBROUTINE evolve_regolith_adsorption
!=======================================================================

!=======================================================================
SUBROUTINE evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,theta_h2o_ads,h2o_ads_reg,delta_h2o_ads)
!-----------------------------------------------------------------------
! NAME
!     evolve_h2o_ads
!
! DESCRIPTION
!     Compute H2O adsorption in regolith following Jackosky et al. (1997).
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,   only: ngrid, nslope, nsoil, nday
use soil,       only: layer, index_breccia
use slopes,     only: subslope_dist, def_slope_mean
use atmosphere, only: ap, bp
use planet,     only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith
use maths,      only: pi

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! Tendencies on the glaciers ()
real(dp), dimension(:,:),   intent(in)    :: h2oice             ! H2O ice at the surface [kg/m^2]
real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
real(dp), dimension(:,:),   intent(in)    :: ps                 ! Surface pressure [Pa]
real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of co2 in the first layer [kg/kg]
real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2o in the first layer [kg/kg]
real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Soil Thermal inertia [J/K/^2/s^1/2]
real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg        ! Density of h2o adsorbed [kg/m^3]
real(dp), dimension(:,:,:), intent(out)   :: theta_h2o_ads      ! Fraction of the pores occupied by H2O molecules
real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of h2o adsorbed [kg/m^3]

! LOCAL VARIABLES
! ---------------
! Constants:
real(dp) :: Ko =  1.57e-8_dp            ! Jackosky et al. 1997
real(dp) :: e = 2573.9_dp               ! Jackosky et al. 1997
real(dp) :: mu = 0.48_dp                ! Jackosky et al. 1997
real(dp) :: m_theta = 2.84e-7_dp        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
! real(dp) :: as = 18.9e3_dp              ! Specific area, Buhler & Piqueux 2021
real(dp) :: as = 9.48e4_dp              ! Specific area, Zent
real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
real(dp), dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
real(dp)                                        :: K                       ! Used to compute theta
integer(di)                                     :: ig, iloop, islope, it   ! For loops
logical(k4),  dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
logical(k4),  dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
real(dp), dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference density of h2o adsorbed per slope [kg/m^3]
real(dp), dimension(ngrid,nsoil,nslope)         :: dm_h2o_regolith_slope   ! Elementary h2o mass adsorded per mesh per slope
real(dp)                                        :: A, B                    ! Used to compute the mean mass above the surface
real(dp)                                        :: p_sat                   ! Saturated vapor pressure of ice
real(dp), dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
real(dp), dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
real(dp), dimension(:,:), allocatable           :: pvapor                  ! Partial pressure above the surface
real(dp), dimension(:)  , allocatable           :: pvapor_avg              ! Yearly averaged

! CODE
! ----
! 0. Some initializations
allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pvapor(ngrid,nday),pvapor_avg(ngrid))
A = 1._dp/m_co2 - 1._dp/m_noco2
B = 1._dp/m_noco2
theta_h2o_ads = 0._dp
dm_h2o_regolith_slope = 0._dp
ispermanent_h2oglaciers = .false.
ispermanent_co2glaciers = .false.

! 0.1 Look at perennial ice
where (abs(d_h2oice(:,:)) > 1.e-5 .and. h2oice(:,:) > 0._dp) ispermanent_h2oglaciers(:,:) = .true.
where (abs(d_co2ice(:,:)) > 1.e-5 .and. co2ice(:,:) > 0._dp) ispermanent_co2glaciers(:,:) = .true.

! 0.2 Compute the partial pressure of vapor
! a. the molecular mass into the column
do ig = 1,ngrid
    mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B)
end do

! b. pressure level
do it = 1,nday
    do ig = 1,ngrid
        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
    end do
end do

! c. Vapor pressure
pvapor = mass_mean/m_h2o*q_h2o_ts*zplev_mean
pvapor_avg = sum(pvapor,2)/nday
deallocate(pvapor,zplev_mean,mass_mean)

! 1. we compute the mass of H2O adsorded in each layer of the meshes
do ig = 1,ngrid
    do islope = 1,nslope
        do iloop = 1,index_breccia
            K = Ko*exp(e/tsoil(ig,iloop,islope))
            if (TI(ig,iloop,islope) < inertie_thresold) then
                theta_h2o_ads(ig,iloop,islope) = (K*pvapor_avg(ig)/(1._dp + K*pvapor_avg(ig)))**mu
            else
                p_sat = exp(beta_clap_h2o/tsoil(ig,iloop,islope) + alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ...
                theta_h2o_ads(ig,iloop,islope) = (K*p_sat/(1._dp + K*p_sat))**mu
            end if
            dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_ads(ig,iloop,islope)*m_theta*rho_regolith
        end do
    end do
end do

! 2. Check the exchange between the atmosphere and the regolith
do ig = 1,ngrid
    delta_h2o_ads(ig) = 0._dp
    do islope = 1,nslope
        deltam_reg_slope(ig,islope) = 0._dp
        do iloop = 1,index_breccia
            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
                if (iloop == 1) then
                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop))
                else
                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1))
                end if
            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
                deltam_reg_complete(ig,iloop,islope) = 0._dp
            end if
            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
        end do
        delta_h2o_ads(ig) = delta_h2o_ads(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180._dp)
    end do
end do
h2o_ads_reg = dm_h2o_regolith_slope

END SUBROUTINE evolve_h2o_ads
!=======================================================================

!=======================================================================
SUBROUTINE evolve_co2_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads)

!-----------------------------------------------------------------------
! NAME
!     evolve_co2_ads
!
! DESCRIPTION
!     Compute CO2 adsorption in regolith following Zent & Quinn (1995).
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry,   only: ngrid, nslope, nsoil, nday
use soil,       only: layer, index_breccia
use slopes,     only: subslope_dist, def_slope_mean
use atmosphere, only: ap, bp
use planet,     only: m_co2, m_noco2, rho_regolith
use maths,      only: pi

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Mean Soil Temperature [K]
real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Mean Thermal Inertia [USI]
real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! Tendencies on the glaciers ()
real(dp), dimension(:,:),   intent(in)    :: q_co2_ts, q_h2o_ts ! Mass mixing ratio of co2 and h2o in the first layer [kg/kg]
real(dp), dimension(:,:),   intent(in)    :: h2oice             ! Water ice at the surface [kg/m^2]
real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
real(dp), dimension(:,:,:), intent(inout) :: co2_ads_reg        ! Density of co2 adsorbed [kg/m^3]
real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of co2 adsorbed [kg/m^3]

! LOCAL VARIABLES
! ---------------
! Constants:
real(dp) :: alpha = 7.512e-6_dp        ! Zent & Quinn 1995
real(dp) :: beta = -1541.5_dp          ! Zent & Quinn 1995
real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
real(dp) :: m_theta = 4.27e-7_dp       ! Mass of co2 per m^2 absorbed
! real(dp) :: as = 18.9e3_dp           ! Specific area, Buhler & Piqueux 2021
real(dp) :: as = 9.48e4_dp             ! Same as previous but from zent
real(dp)                                           :: A, B                    ! Used to compute the mean mass above the surface
integer(di)                                        :: ig, islope, iloop, it   ! For loops
real(dp),    dimension(ngrid,nsoil,nslope)         :: dm_co2_regolith_slope   ! Elementary mass adsorded per mesh per slope
logical(k4), dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
logical(k4), dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
real(dp),    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
real(dp),    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference in the mass per slope  [kg/m^3]
real(dp),    dimension(ngrid,nsoil,nslope)         :: m_h2o_adsorbed          ! Density of CO2 adsorbed [kg/m^3]
real(dp),    dimension(ngrid,nsoil,nslope)         :: theta_h2o_ads           ! Fraction of the pores occupied by H2O molecules
real(dp),    dimension(ngrid)                      :: delta_mh2o              ! Difference density of h2o adsorbed [kg/m^3]
! nday array are allocated because heavy...
real(dp),    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
real(dp),    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
real(dp),    dimension(:,:), allocatable           :: pco2                    ! Partial pressure above the surface
real(dp),    dimension(:),   allocatable           :: pco2_avg                ! Yearly averaged

! CODE
! ----
! 0. Some initializations
allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pco2(ngrid,nday),pco2_avg(ngrid))
m_h2o_adsorbed = 0._dp
A = 1._dp/m_co2 - 1._dp/m_noco2
B = 1._dp/m_noco2
dm_co2_regolith_slope = 0._dp
delta_co2_ads = 0._dp
ispermanent_h2oglaciers = .false.
ispermanent_co2glaciers = .false.

! 0.1 Look at perennial ice
where (abs(d_h2oice(:,:)) > 1.e-5 .and. h2oice(:,:) > 0._dp) ispermanent_h2oglaciers(:,:) = .true.
where (abs(d_co2ice(:,:)) > 1.e-5 .and. co2ice(:,:) > 0._dp) ispermanent_co2glaciers(:,:) = .true.

! 0.2  Compute the partial pressure of CO2
! a. the molecular mass into the column
do ig = 1,ngrid
    mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B)
end do

! b. pressure level
do it = 1,nday
    do ig = 1,ngrid
        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
    end do
end do

! c. Vapor pressure
pco2 = mass_mean/m_co2*q_co2_ts*zplev_mean
pco2_avg(:) = sum(pco2(:,:),2)/nday

deallocate(zplev_mean,mass_mean,pco2)

! 1. Compute the fraction of the pores occupied by H2O
call evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_co2_ts,q_h2o_ts,tsoil,TI,theta_h2o_ads,m_h2o_adsorbed,delta_mh2o)

! 2. we compute the mass of co2 adsorded in each layer of the meshes
do ig = 1,ngrid
    do islope = 1,nslope
        do iloop = 1,index_breccia
            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
                dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
                                                         (alpha*pco2_avg(ig) + sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
            else
                if (abs(co2_ads_reg(ig,iloop,islope)) < minieps) then !!! we are at first call
                    dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig) &
                                                             /(alpha*pco2_avg(ig)+sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
                else ! no change: permanent ice stick the atoms of CO2
                    dm_co2_regolith_slope(ig,iloop,islope) = co2_ads_reg(ig,iloop,islope)
                end if
            end if
        end do
    end do
end do

! 3. Check the exchange between the atmosphere and the regolith
do ig = 1,ngrid
    delta_co2_ads(ig) = 0._dp
    do islope = 1,nslope
        deltam_reg_slope(ig,islope) = 0._dp
        do iloop = 1,index_breccia
            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
                if (iloop == 1) then
                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop))
                else
                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1))
                end if
            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
                deltam_reg_complete(ig,iloop,islope) = 0._dp
            end if
            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
        end do
        delta_co2_ads(ig) = delta_co2_ads(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180._dp)
    end do
end do
co2_ads_reg = dm_co2_regolith_slope

END SUBROUTINE evolve_co2_ads
!=======================================================================

!=======================================================================
SUBROUTINE compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o)
!-----------------------------------------------------------------------
! NAME
!     compute_totmass_adsorbed
!
! DESCRIPTION
!     Compute the total mass of CO2 and H2O adsorbed in the regolith.
!
! AUTHORS & DATE
!     L. Lange, 2023
!     JB Clement, 12/2025
!     C. Metz, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use geometry, only: ngrid, nslope, nsoil, cell_area
use slopes,   only: subslope_dist, def_slope_mean
use soil,     only: layer
use maths,    only: pi
use display,  only: print_msg
use utility,  only: real2str

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

! ARGUMENTS
! ---------
real(dp), dimension(:,:,:), intent(in)  :: h2o_ads_reg, co2_ads_reg
real(qp),                   intent(out) :: totmass_adsco2, totmass_adsh2o

! LOCAL VARIABLES
! ---------------
integer(di) :: i, islope, l
real(dp)    :: thickness, slope_corr

! CODE
! ----
totmass_adsco2 = 0._qp
totmass_adsh2o = 0._qp
do i = 1,ngrid
    do islope = 1,nslope
        slope_corr = subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
        do l = 1,nsoil
            if (l == 1) then
                thickness = layer(l)
            else
                thickness = layer(l) - layer(l-1)
            end if
            totmass_adsco2 = totmass_adsco2 + co2_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i)
            totmass_adsh2o = totmass_adsh2o + h2o_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i)
        end do
    end do
end do
call print_msg("Total mass of CO2 adsorbed in the regolith [kg] = "//real2str(totmass_adsco2))
call print_msg("Total mass of H2O adsorbed in the regolith [kg] = "//real2str(totmass_adsh2o))

END SUBROUTINE compute_totmass_adsorbed
!=======================================================================

END MODULE sorption