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

Last change on this file since 4092 was 4074, checked in by jbclement, 13 days ago

PEM:

  • Correct management of H2O ice tendency in 1D when there is not enough ice anymore.
  • Clean initialization of allocatable module arrays (especially needed when no slope)
  • One more renaming for consistency + few small updates thoughout the code.

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.