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

Last change on this file since 4066 was 4065, checked in by jbclement, 2 weeks ago

PEM:
Major refactor following the previous ones (r3989 and r3991) completing the large structural reorganization and cleanup of the PEM codebase. This revision introduces newly designed modules, standardizes interfaces with explicit ini/end APIs and adds native NetCDF I/O together with explicit PCM/PEM adapters. In detail:

  • Some PEM models were corrected or improved:
    • Frost/perennial ice semantics are clarified via renaming;
    • Soil temperature remapping clarified, notably by removing the rescaling of temperature deviation;
    • Geothermal flux for the PCM is computed based on the PEM state;
  • New explicit PEM/PCM adapters ("set_*"/"build4PCM_*") to decouple PEM internal representation from PCM file layouts and reconstruct consistent fields returned to the PCM;
  • New explicit build/teardown routines that centralize allocation and initialization ordering, reducing accidental use of uninitialized data and making the model lifecycle explicit;
  • Add native read/write helpers for NetCDF that centralize all low-level NetCDF interactions with major improvements (and more simplicity) compared to legacy PEM/PCM I/O (see the modules "io_netcdf" and "output"). They support reading, creation and writing of "diagevol.nc" (renamed from "diagpem.nc") and start/restart files;
  • Provide well-focused modules ("numerics"/"maths"/"utility"/"display") to host commonly-used primitives:
    • "numerics" defines numerical types and constants for reproducibility, portability across compilers and future transitions (e.g. quadruple precision experiments);
    • "display" provides a single controlled interface for runtime messages, status output and diagnostics, avoiding direct 'print'/'write' to enable silent mode, log redirection, and MPI-safe output in the future.
    • "utility" (new module) hosts generic helpers used throughout the code (e.g. "int2str" or "real2str");
  • Add modules "clim_state_init"/"clim_state_rec" which provide robust read/write logic for "start/startfi/startpem", including 1D fallbacks, mesh conversions and dimension checks. NetCDF file creation is centralized and explicit. Restart files are now self-consistent and future-proof, requiring changes only to affected variables;
  • Add module "atmosphere" which computes pressure fields, reconstructs potential temperature and air mass. It also holds the whole logic to define sigma or hybrid coordinates for altitudes;
  • Add module "geometry" to centrilize dimensions logic and grid conversions routines (including 2 new ones "dyngrd2vect"/"vect2dyngrd");
  • Add module "slopes" to isolate slopes handling;
  • Add module "surface" to isolate surface management. Notably, albedo and emissivity are now fully reconstructed following the PCM settings;
  • Add module "allocation" to check the dimension initialization and centrilize allocation/deallocation;
  • Finalize module decomposition and renaming to consolidate domain-based modules, purpose-based routines and physics/process-based variables;
  • The main program now drives a clearer sequence of conceptual steps (initialization / reading / evolution / update / build / writing) and fails explicitly instead of silently defaulting;
  • Ice table logic is made restart-safe;
  • 'Open'/'read' intrinsic logic is made safe and automatic;
  • Improve discoverability and standardize the data handling (private vs protected vs public);
  • Apply consistent documentation/header style already introduced;
  • Update deftank scripts to reflect new names and launch wrappers;

This revision is a structural milestone aiming to be behavior-preserving where possible. It has been tested via compilation and short integration runs. However, due to extensive renames, moves, and API changes, full validation is still ongoing.
Note: the revision includes one (possibly two) easter egg hidden in the code for future archaeologists of the PEM. No physics were harmed.
JBC
File size: 20.1 KB
Line 
1MODULE sorption
2!-----------------------------------------------------------------------
3! NAME
4!     sorption
5!
6! DESCRIPTION
7!     CO2 and H2O adsorption/desorption in regolith following Zent & Quinn
8!     (1995) and Jackosky et al. (1997).
9!
10! AUTHORS & DATE
11!     L. Lange, 2023
12!     JB Clement, 02/2026
13!
14! NOTES
15!
16!-----------------------------------------------------------------------
17
18! DEPENDENCIES
19! ------------
20use numerics, only: dp, qp, di, k4, minieps
21
22! DECLARATION
23! -----------
24implicit none
25
26! PARAMETERS
27! ----------
28logical(k4), protected :: do_sorption ! Flag to compute adsorption/desorption
29
30contains
31!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32
33!=======================================================================
34SUBROUTINE set_sorption_config(do_sorption_in)
35!-----------------------------------------------------------------------
36! NAME
37!     set_sorption_config
38!
39! DESCRIPTION
40!     Setter for 'sorption' configuration parameters.
41!
42! AUTHORS & DATE
43!     JB Clement, 02/2026
44!
45! NOTES
46!
47!-----------------------------------------------------------------------
48
49! DEPENDENCIES
50! ------------
51use utility, only: bool2str
52use display, only: print_msg
53
54! DECLARATION
55! -----------
56implicit none
57
58! ARGUMENTS
59! ---------
60logical(k4), intent(in) :: do_sorption_in
61
62! CODE
63! ----
64do_sorption = do_sorption_in
65call print_msg('do_sorption = '//bool2str(do_sorption))
66
67END SUBROUTINE set_sorption_config
68!=======================================================================
69
70!=======================================================================
71SUBROUTINE 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)
72!-----------------------------------------------------------------------
73! NAME
74!     evolve_regolith_adsorption
75!
76! DESCRIPTION
77!     Compute both H2O and CO2 adsorption in regolith.
78!
79! AUTHORS & DATE
80!     L. Lange, 2023
81!     JB Clement, 12/2025
82!     C. Metz, 02/2026
83!
84! NOTES
85!
86!-----------------------------------------------------------------------
87
88! DEPENDENCIES
89! ------------
90use geometry, only: ngrid, nslope, nsoil
91
92! DECLARATION
93! -----------
94implicit none
95
96! ARGUMENTS
97! ---------
98real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! tendencies on the glaciers [1]
99real(dp), dimension(:,:),   intent(in)    :: h2oice             ! H2O ice at the surface [kg/m^2]
100real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
101real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Soil Thermal inertia [J/K/^2/s^1/2]
102real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
103real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
104real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of co2 in the first layer [kg/kg]
105real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2o in the first layer [kg/kg]
106real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of h2o adsorbed [kg/m^3]
107real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of co2 adsorbed [kg/m^3]
108real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg
109
110! LOCAL VARIABLES
111! ---------------
112real(dp), dimension(ngrid,nsoil,nslope) :: theta_h2o_ads ! Fraction of the pores occupied by H2O molecules
113
114! CODE
115! ----
116call 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)
117call evolve_co2_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads)
118
119END SUBROUTINE evolve_regolith_adsorption
120!=======================================================================
121
122!=======================================================================
123SUBROUTINE 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)
124!-----------------------------------------------------------------------
125! NAME
126!     evolve_h2o_ads
127!
128! DESCRIPTION
129!     Compute H2O adsorption in regolith following Jackosky et al. (1997).
130!
131! AUTHORS & DATE
132!     L. Lange, 2023
133!     JB Clement, 2025
134!
135! NOTES
136!
137!-----------------------------------------------------------------------
138
139! DEPENDENCIES
140! ------------
141use geometry,   only: ngrid, nslope, nsoil, nday
142use soil,       only: layer, index_breccia
143use slopes,     only: subslope_dist, def_slope_mean
144use atmosphere, only: ap, bp
145use planet,     only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith
146use maths,      only: pi
147
148! DECLARATION
149! -----------
150implicit none
151
152! ARGUMENTS
153! ---------
154real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! Tendencies on the glaciers ()
155real(dp), dimension(:,:),   intent(in)    :: h2oice             ! H2O ice at the surface [kg/m^2]
156real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
157real(dp), dimension(:,:),   intent(in)    :: ps                 ! Surface pressure [Pa]
158real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of co2 in the first layer [kg/kg]
159real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2o in the first layer [kg/kg]
160real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Soil Thermal inertia [J/K/^2/s^1/2]
161real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
162real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg        ! Density of h2o adsorbed [kg/m^3]
163real(dp), dimension(:,:,:), intent(out)   :: theta_h2o_ads      ! Fraction of the pores occupied by H2O molecules
164real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of h2o adsorbed [kg/m^3]
165
166! LOCAL VARIABLES
167! ---------------
168! Constants:
169real(dp) :: Ko =  1.57e-8_dp            ! Jackosky et al. 1997
170real(dp) :: e = 2573.9_dp               ! Jackosky et al. 1997
171real(dp) :: mu = 0.48_dp                ! Jackosky et al. 1997
172real(dp) :: m_theta = 2.84e-7_dp        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
173! real(dp) :: as = 18.9e3_dp              ! Specific area, Buhler & Piqueux 2021
174real(dp) :: as = 9.48e4_dp              ! Specific area, Zent
175real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
176real(dp), dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
177real(dp)                                        :: K                       ! Used to compute theta
178integer(di)                                     :: ig, iloop, islope, it   ! For loops
179logical(k4),  dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
180logical(k4),  dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
181real(dp), dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference density of h2o adsorbed per slope [kg/m^3]
182real(dp), dimension(ngrid,nsoil,nslope)         :: dm_h2o_regolith_slope   ! Elementary h2o mass adsorded per mesh per slope
183real(dp)                                        :: A, B                    ! Used to compute the mean mass above the surface
184real(dp)                                        :: p_sat                   ! Saturated vapor pressure of ice
185real(dp), dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
186real(dp), dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
187real(dp), dimension(:,:), allocatable           :: pvapor                  ! Partial pressure above the surface
188real(dp), dimension(:)  , allocatable           :: pvapor_avg              ! Yearly averaged
189
190! CODE
191! ----
192! 0. Some initializations
193allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pvapor(ngrid,nday),pvapor_avg(ngrid))
194A = 1._dp/m_co2 - 1._dp/m_noco2
195B = 1._dp/m_noco2
196theta_h2o_ads = 0._dp
197dm_h2o_regolith_slope = 0._dp
198ispermanent_h2oglaciers = .false.
199ispermanent_co2glaciers = .false.
200
201! 0.1 Look at perennial ice
202where (abs(d_h2oice(:,:)) > 1.e-5 .and. h2oice(:,:) > 0._dp) ispermanent_h2oglaciers(:,:) = .true.
203where (abs(d_co2ice(:,:)) > 1.e-5 .and. co2ice(:,:) > 0._dp) ispermanent_co2glaciers(:,:) = .true.
204
205! 0.2 Compute the partial pressure of vapor
206! a. the molecular mass into the column
207do ig = 1,ngrid
208    mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B)
209end do
210
211! b. pressure level
212do it = 1,nday
213    do ig = 1,ngrid
214        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
215    end do
216end do
217
218! c. Vapor pressure
219pvapor = mass_mean/m_h2o*q_h2o_ts*zplev_mean
220pvapor_avg = sum(pvapor,2)/nday
221deallocate(pvapor,zplev_mean,mass_mean)
222
223! 1. we compute the mass of H2O adsorded in each layer of the meshes
224do ig = 1,ngrid
225    do islope = 1,nslope
226        do iloop = 1,index_breccia
227            K = Ko*exp(e/tsoil(ig,iloop,islope))
228            if (TI(ig,iloop,islope) < inertie_thresold) then
229                theta_h2o_ads(ig,iloop,islope) = (K*pvapor_avg(ig)/(1._dp + K*pvapor_avg(ig)))**mu
230            else
231                p_sat = exp(beta_clap_h2o/tsoil(ig,iloop,islope) + alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ...
232                theta_h2o_ads(ig,iloop,islope) = (K*p_sat/(1._dp + K*p_sat))**mu
233            end if
234            dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_ads(ig,iloop,islope)*m_theta*rho_regolith
235        end do
236    end do
237end do
238
239! 2. Check the exchange between the atmosphere and the regolith
240do ig = 1,ngrid
241    delta_h2o_ads(ig) = 0._dp
242    do islope = 1,nslope
243        deltam_reg_slope(ig,islope) = 0._dp
244        do iloop = 1,index_breccia
245            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
246                if (iloop == 1) then
247                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop))
248                else
249                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1))
250                end if
251            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
252                deltam_reg_complete(ig,iloop,islope) = 0._dp
253            end if
254            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
255        end do
256        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)
257    end do
258end do
259h2o_ads_reg = dm_h2o_regolith_slope
260
261END SUBROUTINE evolve_h2o_ads
262!=======================================================================
263
264!=======================================================================
265SUBROUTINE evolve_co2_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads)
266
267!-----------------------------------------------------------------------
268! NAME
269!     evolve_co2_ads
270!
271! DESCRIPTION
272!     Compute CO2 adsorption in regolith following Zent & Quinn (1995).
273!
274! AUTHORS & DATE
275!     L. Lange, 2023
276!     JB Clement, 2025
277!
278! NOTES
279!
280!-----------------------------------------------------------------------
281
282! DEPENDENCIES
283! ------------
284use geometry,   only: ngrid, nslope, nsoil, nday
285use soil,       only: layer, index_breccia
286use slopes,     only: subslope_dist, def_slope_mean
287use atmosphere, only: ap, bp
288use planet,     only: m_co2, m_noco2, rho_regolith
289use maths,      only: pi
290
291! DECLARATION
292! -----------
293implicit none
294
295! ARGUMENTS
296! ---------
297real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
298real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Mean Soil Temperature [K]
299real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Mean Thermal Inertia [USI]
300real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! Tendencies on the glaciers ()
301real(dp), dimension(:,:),   intent(in)    :: q_co2_ts, q_h2o_ts ! Mass mixing ratio of co2 and h2o in the first layer [kg/kg]
302real(dp), dimension(:,:),   intent(in)    :: h2oice             ! Water ice at the surface [kg/m^2]
303real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
304real(dp), dimension(:,:,:), intent(inout) :: co2_ads_reg        ! Density of co2 adsorbed [kg/m^3]
305real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of co2 adsorbed [kg/m^3]
306
307! LOCAL VARIABLES
308! ---------------
309! Constants:
310real(dp) :: alpha = 7.512e-6_dp        ! Zent & Quinn 1995
311real(dp) :: beta = -1541.5_dp          ! Zent & Quinn 1995
312real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
313real(dp) :: m_theta = 4.27e-7_dp       ! Mass of co2 per m^2 absorbed
314! real(dp) :: as = 18.9e3_dp           ! Specific area, Buhler & Piqueux 2021
315real(dp) :: as = 9.48e4_dp             ! Same as previous but from zent
316real(dp)                                           :: A, B                    ! Used to compute the mean mass above the surface
317integer(di)                                        :: ig, islope, iloop, it   ! For loops
318real(dp),    dimension(ngrid,nsoil,nslope)         :: dm_co2_regolith_slope   ! Elementary mass adsorded per mesh per slope
319logical(k4), dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
320logical(k4), dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
321real(dp),    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
322real(dp),    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference in the mass per slope  [kg/m^3]
323real(dp),    dimension(ngrid,nsoil,nslope)         :: m_h2o_adsorbed          ! Density of CO2 adsorbed [kg/m^3]
324real(dp),    dimension(ngrid,nsoil,nslope)         :: theta_h2o_ads           ! Fraction of the pores occupied by H2O molecules
325real(dp),    dimension(ngrid)                      :: delta_mh2o              ! Difference density of h2o adsorbed [kg/m^3]
326! nday array are allocated because heavy...
327real(dp),    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
328real(dp),    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
329real(dp),    dimension(:,:), allocatable           :: pco2                    ! Partial pressure above the surface
330real(dp),    dimension(:),   allocatable           :: pco2_avg                ! Yearly averaged
331
332! CODE
333! ----
334! 0. Some initializations
335allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pco2(ngrid,nday),pco2_avg(ngrid))
336m_h2o_adsorbed = 0._dp
337A = 1._dp/m_co2 - 1._dp/m_noco2
338B = 1._dp/m_noco2
339dm_co2_regolith_slope = 0._dp
340delta_co2_ads = 0._dp
341ispermanent_h2oglaciers = .false.
342ispermanent_co2glaciers = .false.
343
344! 0.1 Look at perennial ice
345where (abs(d_h2oice(:,:)) > 1.e-5 .and. h2oice(:,:) > 0._dp) ispermanent_h2oglaciers(:,:) = .true.
346where (abs(d_co2ice(:,:)) > 1.e-5 .and. co2ice(:,:) > 0._dp) ispermanent_co2glaciers(:,:) = .true.
347
348! 0.2  Compute the partial pressure of CO2
349! a. the molecular mass into the column
350do ig = 1,ngrid
351    mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B)
352end do
353
354! b. pressure level
355do it = 1,nday
356    do ig = 1,ngrid
357        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
358    end do
359end do
360
361! c. Vapor pressure
362pco2 = mass_mean/m_co2*q_co2_ts*zplev_mean
363pco2_avg(:) = sum(pco2(:,:),2)/nday
364
365deallocate(zplev_mean,mass_mean,pco2)
366
367! 1. Compute the fraction of the pores occupied by H2O
368call 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)
369
370! 2. we compute the mass of co2 adsorded in each layer of the meshes
371do ig = 1,ngrid
372    do islope = 1,nslope
373        do iloop = 1,index_breccia
374            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
375                dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
376                                                         (alpha*pco2_avg(ig) + sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
377            else
378                if (abs(co2_ads_reg(ig,iloop,islope)) < minieps) then !!! we are at first call
379                    dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig) &
380                                                             /(alpha*pco2_avg(ig)+sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
381                else ! no change: permanent ice stick the atoms of CO2
382                    dm_co2_regolith_slope(ig,iloop,islope) = co2_ads_reg(ig,iloop,islope)
383                end if
384            end if
385        end do
386    end do
387end do
388
389! 3. Check the exchange between the atmosphere and the regolith
390do ig = 1,ngrid
391    delta_co2_ads(ig) = 0._dp
392    do islope = 1,nslope
393        deltam_reg_slope(ig,islope) = 0._dp
394        do iloop = 1,index_breccia
395            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
396                if (iloop == 1) then
397                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop))
398                else
399                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1))
400                end if
401            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
402                deltam_reg_complete(ig,iloop,islope) = 0._dp
403            end if
404            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
405        end do
406        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)
407    end do
408end do
409co2_ads_reg = dm_co2_regolith_slope
410
411END SUBROUTINE evolve_co2_ads
412!=======================================================================
413
414!=======================================================================
415SUBROUTINE compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o)
416!-----------------------------------------------------------------------
417! NAME
418!     compute_totmass_adsorbed
419!
420! DESCRIPTION
421!     Compute the total mass of CO2 and H2O adsorbed in the regolith.
422!
423! AUTHORS & DATE
424!     L. Lange, 2023
425!     JB Clement, 12/2025
426!     C. Metz, 02/2026
427!
428! NOTES
429!
430!-----------------------------------------------------------------------
431
432! DEPENDENCIES
433! ------------
434use geometry, only: ngrid, nslope, nsoil, cell_area
435use slopes,   only: subslope_dist, def_slope_mean
436use soil,     only: layer
437use maths,    only: pi
438use display,  only: print_msg
439use utility,  only: real2str
440
441! DECLARATION
442! -----------
443implicit none
444
445! ARGUMENTS
446! ---------
447real(dp), dimension(:,:,:), intent(in)  :: h2o_ads_reg, co2_ads_reg
448real(qp),                   intent(out) :: totmass_adsco2, totmass_adsh2o
449
450! LOCAL VARIABLES
451! ---------------
452integer(di) :: i, islope, l
453real(dp)    :: thickness, slope_corr
454
455! CODE
456! ----
457totmass_adsco2 = 0._qp
458totmass_adsh2o = 0._qp
459do i = 1,ngrid
460    do islope = 1,nslope
461        slope_corr = subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
462        do l = 1,nsoil
463            if (l == 1) then
464                thickness = layer(l)
465            else
466                thickness = layer(l) - layer(l-1)
467            end if
468            totmass_adsco2 = totmass_adsco2 + co2_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i)
469            totmass_adsh2o = totmass_adsh2o + h2o_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i)
470        end do
471    end do
472end do
473call print_msg("Total mass of CO2 adsorbed in the regolith [kg] = "//real2str(totmass_adsco2))
474call print_msg("Total mass of H2O adsorbed in the regolith [kg] = "//real2str(totmass_adsh2o))
475
476END SUBROUTINE compute_totmass_adsorbed
477!=======================================================================
478
479END MODULE sorption
Note: See TracBrowser for help on using the repository browser.