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

Last change on this file since 4138 was 4135, checked in by jbclement, 13 days ago

PEM:
Relocate Mars-specific parameters out of "planet" module and into the modules that own their physics domain.
JBC
File size: 18.3 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, LVL_NFO
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),LVL_NFO)
66
67END SUBROUTINE set_sorption_config
68!=======================================================================
69
70!=======================================================================
71SUBROUTINE evolve_regolith_adsorption(h2o_ice,co2_ice,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)    :: h2o_ice            ! H2O ice at the surface [kg/m^2]
99real(dp), dimension(:,:),   intent(in)    :: co2_ice            ! CO2 ice at the surface [kg/m^2]
100real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Soil Thermal inertia [J/K/^2/s^1/2]
101real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
102real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
103real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of CO2 in the first layer [kg/kg]
104real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2O in the first layer [kg/kg]
105real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of H2O adsorbed [kg/m^3]
106real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of CO2 adsorbed [kg/m^3]
107real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg
108
109! LOCAL VARIABLES
110! ---------------
111real(dp), dimension(ngrid,nsoil,nslope) :: theta_h2o_ads ! Fraction of the pores occupied by H2O molecules
112
113! CODE
114! ----
115call evolve_h2o_ads(h2o_ice,co2_ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,theta_h2o_ads,h2o_ads_reg,delta_h2o_ads)
116call evolve_co2_ads(h2o_ice,co2_ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads)
117
118END SUBROUTINE evolve_regolith_adsorption
119!=======================================================================
120
121!=======================================================================
122SUBROUTINE evolve_h2o_ads(h2o_ice,co2_ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,theta_h2o_ads,h2o_ads_reg,delta_h2o_ads)
123!-----------------------------------------------------------------------
124! NAME
125!     evolve_h2o_ads
126!
127! DESCRIPTION
128!     Compute H2O adsorption in regolith following Jackosky et al. (1997).
129!
130! AUTHORS & DATE
131!     L. Lange, 2023
132!     JB Clement, 2025
133!
134! NOTES
135!
136!-----------------------------------------------------------------------
137
138! DEPENDENCIES
139! ------------
140use geometry,   only: ngrid, nslope, nsoil, nday
141use soil,       only: layer, index_breccia, rho_regolith
142use slopes,     only: subslope_dist, def_slope_mean
143use atmosphere, only: ap, bp
144use physics,    only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2, A, B
145use maths,      only: pi
146
147! DECLARATION
148! -----------
149implicit none
150
151! ARGUMENTS
152! ---------
153real(dp), dimension(:,:),   intent(in)    :: h2o_ice       ! H2O ice at the surface [kg/m^2]
154real(dp), dimension(:,:),   intent(in)    :: co2_ice       ! CO2 ice at the surface [kg/m^2]
155real(dp), dimension(:,:),   intent(in)    :: ps            ! Surface pressure [Pa]
156real(dp), dimension(:,:),   intent(in)    :: q_co2_ts      ! Mass mixing ratio of CO2 in the first layer [kg/kg]
157real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts      ! Mass mixing ratio of H2O in the first layer [kg/kg]
158real(dp), dimension(:,:,:), intent(in)    :: TI            ! Soil Thermal inertia [J/K/^2/s^1/2]
159real(dp), dimension(:,:,:), intent(in)    :: tsoil         ! Soil temperature [K]
160real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg   ! Density of H2O adsorbed [kg/m^3]
161real(dp), dimension(:,:,:), intent(out)   :: theta_h2o_ads ! Fraction of the pores occupied by H2O molecules
162real(dp), dimension(:),     intent(out)   :: delta_h2o_ads ! Difference density of H2O adsorbed [kg/m^3]
163
164! LOCAL VARIABLES
165! ---------------
166! Constants:
167real(dp) :: Ko =  1.57e-8_dp            ! Jackosky et al. 1997
168real(dp) :: e = 2573.9_dp               ! Jackosky et al. 1997
169real(dp) :: mu = 0.48_dp                ! Jackosky et al. 1997
170real(dp) :: m_theta = 2.84e-7_dp        ! Mass of H2O per m^2 absorbed Jackosky et al. 1997
171! real(dp) :: as = 18.9e3_dp              ! Specific area, Buhler & Piqueux 2021
172real(dp) :: as = 9.48e4_dp              ! Specific area, Zent
173real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
174real(dp), dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
175real(dp)                                        :: K                       ! Used to compute theta
176integer(di)                                     :: ig, iloop, islope, it   ! For loops
177real(dp),    dimension(ngrid,nslope)            :: deltam_reg_slope        ! Difference density of H2O adsorbed per slope [kg/m^3]
178real(dp),    dimension(ngrid,nsoil,nslope)      :: dm_h2o_regolith_slope   ! Elementary H2O mass adsorded per mesh per slope
179real(dp)                                        :: p_sat                   ! Saturated vapor pressure of ice
180real(dp),    dimension(:,:), allocatable        :: mass_mean               ! Mean mass above the surface
181real(dp),    dimension(:,:), allocatable        :: zplev_mean              ! Pressure above the surface
182real(dp),    dimension(:,:), allocatable        :: pvapor                  ! Partial pressure above the surface
183real(dp),    dimension(:)  , allocatable        :: pvapor_avg              ! Yearly averaged
184
185! CODE
186! ----
187! 0. Some initializations
188allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pvapor(ngrid,nday),pvapor_avg(ngrid))
189theta_h2o_ads = 0._dp
190dm_h2o_regolith_slope = 0._dp
191
192! 0. Compute the partial pressure of vapor
193! a. the molecular mass into the column
194do ig = 1,ngrid
195    mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B)
196end do
197
198! b. pressure level
199do it = 1,nday
200    do ig = 1,ngrid
201        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
202    end do
203end do
204
205! c. Vapor pressure
206pvapor = mass_mean/m_h2o*q_h2o_ts*zplev_mean
207pvapor_avg = sum(pvapor,2)/nday
208deallocate(pvapor,zplev_mean,mass_mean)
209
210! 1. we compute the mass of H2O adsorded in each layer of the meshes
211do ig = 1,ngrid
212    do islope = 1,nslope
213        do iloop = 1,index_breccia
214            K = Ko*exp(e/tsoil(ig,iloop,islope))
215            if (TI(ig,iloop,islope) < inertie_thresold) then
216                theta_h2o_ads(ig,iloop,islope) = (K*pvapor_avg(ig)/(1._dp + K*pvapor_avg(ig)))**mu
217            else
218                p_sat = exp(beta_clap_h2o/tsoil(ig,iloop,islope) + alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ...
219                theta_h2o_ads(ig,iloop,islope) = (K*p_sat/(1._dp + K*p_sat))**mu
220            end if
221            dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_ads(ig,iloop,islope)*m_theta*rho_regolith
222        end do
223    end do
224end do
225
226! 2. Check the exchange between the atmosphere and the regolith
227do ig = 1,ngrid
228    delta_h2o_ads(ig) = 0._dp
229    do islope = 1,nslope
230        deltam_reg_slope(ig,islope) = 0._dp
231        do iloop = 1,index_breccia
232            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < minieps .and. abs(co2_ice(ig,islope)) < minieps) then
233                if (iloop == 1) then
234                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop))
235                else
236                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1))
237                end if
238            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
239                deltam_reg_complete(ig,iloop,islope) = 0._dp
240            end if
241            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
242        end do
243        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)
244    end do
245end do
246h2o_ads_reg = dm_h2o_regolith_slope
247
248END SUBROUTINE evolve_h2o_ads
249!=======================================================================
250
251!=======================================================================
252SUBROUTINE evolve_co2_ads(h2o_ice,co2_ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads)
253
254!-----------------------------------------------------------------------
255! NAME
256!     evolve_co2_ads
257!
258! DESCRIPTION
259!     Compute CO2 adsorption in regolith following Zent & Quinn (1995).
260!
261! AUTHORS & DATE
262!     L. Lange, 2023
263!     JB Clement, 2025
264!
265! NOTES
266!
267!-----------------------------------------------------------------------
268
269! DEPENDENCIES
270! ------------
271use geometry,   only: ngrid, nslope, nsoil, nday
272use soil,       only: layer, index_breccia, rho_regolith
273use slopes,     only: subslope_dist, def_slope_mean
274use atmosphere, only: ap, bp
275use physics,    only: m_co2, A, B
276use maths,      only: pi
277
278! DECLARATION
279! -----------
280implicit none
281
282! ARGUMENTS
283! ---------
284real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
285real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Mean Soil Temperature [K]
286real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Mean Thermal Inertia [USI]
287real(dp), dimension(:,:),   intent(in)    :: q_co2_ts, q_h2o_ts ! Mass mixing ratio of CO2 and H2O in the first layer [kg/kg]
288real(dp), dimension(:,:),   intent(in)    :: h2o_ice            ! Water ice at the surface [kg/m^2]
289real(dp), dimension(:,:),   intent(in)    :: co2_ice            ! CO2 ice at the surface [kg/m^2]
290real(dp), dimension(:,:,:), intent(inout) :: co2_ads_reg        ! Density of CO2 adsorbed [kg/m^3]
291real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of CO2 adsorbed [kg/m^3]
292
293! LOCAL VARIABLES
294! ---------------
295! Constants:
296real(dp) :: alpha = 7.512e-6_dp        ! Zent & Quinn 1995
297real(dp) :: beta = -1541.5_dp          ! Zent & Quinn 1995
298real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
299real(dp) :: m_theta = 4.27e-7_dp       ! Mass of CO2 per m^2 absorbed
300! real(dp) :: as = 18.9e3_dp           ! Specific area, Buhler & Piqueux 2021
301real(dp) :: as = 9.48e4_dp             ! Same as previous but from zent
302integer(di)                                        :: ig, islope, iloop, it   ! For loops
303real(dp),    dimension(ngrid,nsoil,nslope)         :: dm_co2_regolith_slope   ! Elementary mass adsorded per mesh per slope
304real(dp),    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
305real(dp),    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference in the mass per slope  [kg/m^3]
306real(dp),    dimension(ngrid,nsoil,nslope)         :: m_h2o_adsorbed          ! Density of CO2 adsorbed [kg/m^3]
307real(dp),    dimension(ngrid,nsoil,nslope)         :: theta_h2o_ads           ! Fraction of the pores occupied by H2O molecules
308real(dp),    dimension(ngrid)                      :: delta_mh2o              ! Difference density of H2O adsorbed [kg/m^3]
309! nday array are allocated because heavy...
310real(dp),    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
311real(dp),    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
312real(dp),    dimension(:,:), allocatable           :: pco2                    ! Partial pressure above the surface
313real(dp),    dimension(:),   allocatable           :: pco2_avg                ! Yearly averaged
314
315! CODE
316! ----
317! 0. Some initializations
318allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pco2(ngrid,nday),pco2_avg(ngrid))
319m_h2o_adsorbed = 0._dp
320dm_co2_regolith_slope = 0._dp
321delta_co2_ads = 0._dp
322
323! 0.  Compute the partial pressure of CO2
324! a. the molecular mass into the column
325do ig = 1,ngrid
326    mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B)
327end do
328
329! b. pressure level
330do it = 1,nday
331    do ig = 1,ngrid
332        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
333    end do
334end do
335
336! c. Vapor pressure
337pco2 = mass_mean/m_co2*q_co2_ts*zplev_mean
338pco2_avg(:) = sum(pco2(:,:),2)/nday
339
340deallocate(zplev_mean,mass_mean,pco2)
341
342! 1. Compute the fraction of the pores occupied by H2O
343call evolve_h2o_ads(h2o_ice,co2_ice,ps,q_co2_ts,q_h2o_ts,tsoil,TI,theta_h2o_ads,m_h2o_adsorbed,delta_mh2o)
344
345! 2. we compute the mass of CO2 adsorded in each layer of the meshes
346do ig = 1,ngrid
347    do islope = 1,nslope
348        do iloop = 1,index_breccia
349            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < minieps .and. abs(co2_ice(ig,islope)) < minieps) then
350                dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
351                                                         (alpha*pco2_avg(ig) + sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
352            else
353                if (abs(co2_ads_reg(ig,iloop,islope)) < minieps) then !!! we are at first call
354                    dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig) &
355                                                             /(alpha*pco2_avg(ig)+sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
356                else ! no change: permanent ice stick the atoms of CO2
357                    dm_co2_regolith_slope(ig,iloop,islope) = co2_ads_reg(ig,iloop,islope)
358                end if
359            end if
360        end do
361    end do
362end do
363
364! 3. Check the exchange between the atmosphere and the regolith
365do ig = 1,ngrid
366    delta_co2_ads(ig) = 0._dp
367    do islope = 1,nslope
368        deltam_reg_slope(ig,islope) = 0._dp
369        do iloop = 1,index_breccia
370            if (TI(ig,iloop,islope) < inertie_thresold .and. abs(h2o_ice(ig,islope)) < minieps .and. abs(co2_ice(ig,islope)) < minieps) then
371                if (iloop == 1) then
372                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop))
373                else
374                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1))
375                end if
376            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
377                deltam_reg_complete(ig,iloop,islope) = 0._dp
378            end if
379            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
380        end do
381        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)
382    end do
383end do
384co2_ads_reg = dm_co2_regolith_slope
385
386END SUBROUTINE evolve_co2_ads
387!=======================================================================
388
389!=======================================================================
390SUBROUTINE compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o)
391!-----------------------------------------------------------------------
392! NAME
393!     compute_totmass_adsorbed
394!
395! DESCRIPTION
396!     Compute the total mass of CO2 and H2O adsorbed in the regolith.
397!
398! AUTHORS & DATE
399!     L. Lange, 2023
400!     JB Clement, 12/2025
401!     C. Metz, 02/2026
402!
403! NOTES
404!
405!-----------------------------------------------------------------------
406
407! DEPENDENCIES
408! ------------
409use geometry, only: ngrid, nslope, nsoil, cell_area
410use slopes,   only: subslope_dist, def_slope_mean
411use soil,     only: layer
412use maths,    only: pi
413use display,  only: print_msg, LVL_NFO
414use utility,  only: real2str
415
416! DECLARATION
417! -----------
418implicit none
419
420! ARGUMENTS
421! ---------
422real(dp), dimension(:,:,:), intent(in)  :: h2o_ads_reg, co2_ads_reg
423real(qp),                   intent(out) :: totmass_adsco2, totmass_adsh2o
424
425! LOCAL VARIABLES
426! ---------------
427integer(di) :: i, islope, l
428real(dp)    :: thickness, slope_corr
429
430! CODE
431! ----
432totmass_adsco2 = 0._qp
433totmass_adsh2o = 0._qp
434do i = 1,ngrid
435    do islope = 1,nslope
436        slope_corr = subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
437        do l = 1,nsoil
438            if (l == 1) then
439                thickness = layer(l)
440            else
441                thickness = layer(l) - layer(l-1)
442            end if
443            totmass_adsco2 = totmass_adsco2 + co2_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i)
444            totmass_adsh2o = totmass_adsh2o + h2o_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i)
445        end do
446    end do
447end do
448call print_msg("Total mass of CO2 adsorbed in the regolith [kg] = "//real2str(totmass_adsco2),LVL_NFO)
449call print_msg("Total mass of H2O adsorbed in the regolith [kg] = "//real2str(totmass_adsh2o),LVL_NFO)
450
451END SUBROUTINE compute_totmass_adsorbed
452!=======================================================================
453
454END MODULE sorption
Note: See TracBrowser for help on using the repository browser.