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

Last change on this file since 3989 was 3989, checked in by jbclement, 5 weeks ago

PEM:
Massive structural refactor of the PEM codebase for improved readability, consistency and maintainability. The goal is to modernize, standardize and consolidate the code while removing legacy complexity. In detail, this change:

  • Performs large-scale cleanup removing unused code, obsolete routines, duplicated functionality and deprecated initialization logic;
  • Removes "*_mod" wrappers;
  • Replaces mixed naming conventions with clearer variable names, domain-based file/module names and purpose-based routine names;
  • Adds native reading/writing capabilities to the PEM removing the cumbersome dependency on Mars PCM subroutines;
  • Centralizes soil/ice/orbital/PEM configuration logic into dedicated modules;
  • Simplifies routines and updates calls/interfaces to match the new structure.

This refactor significantly clarifies the codebase and provides a cleaner foundation for forthcoming developments.
JBC
File size: 18.9 KB
Line 
1MODULE sorption
2
3implicit none
4
5logical                             :: do_sorption       ! Flag to compute adsorption/desorption
6real, allocatable, dimension(:,:,:) :: co2_adsorbed_phys ! co2 that is in the regolith [kg/m^2]
7real, allocatable, dimension(:,:,:) :: h2o_adsorbed_phys ! h2o that is in the regolith [kg/m^2]
8
9!=======================================================================
10contains
11!=======================================================================
12
13!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
14!!!
15!!! Purpose: Compute CO2 and H2O adsorption, following the methods from Zent & Quinn 1995, Jackosky et al., 1997
16!!!
17!!! Author: LL, 01/2023
18!!!
19!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20SUBROUTINE ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
21
22implicit none
23
24integer, intent(in) :: ngrid       ! number of atmospheric columns
25integer, intent(in) :: nslope      ! number of slope within a mesh
26integer, intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM
27
28allocate(co2_adsorbed_phys(ngrid,nsoilmx_PEM,nslope))
29allocate(h2o_adsorbed_phys(ngrid,nsoilmx_PEM,nslope))
30
31END SUBROUTINE ini_adsorption_h_PEM
32
33!=======================================================================
34SUBROUTINE end_adsorption_h_PEM
35
36if (allocated(co2_adsorbed_phys)) deallocate(co2_adsorbed_phys)
37if (allocated(h2o_adsorbed_phys)) deallocate(h2o_adsorbed_phys)
38
39END SUBROUTINE end_adsorption_h_PEM
40
41!=======================================================================
42SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
43                               m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg)
44
45implicit none
46
47! Inputs
48integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM, timelen  ! size dimension: physics x subslope x soil x timeseries
49real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! tendancies on the glaciers [1]
50real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! water ice at the surface [kg/m^2]
51real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! co2 ice at the surface [kg/m^2]
52real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Soil Thermal inertia (J/K/^2/s^1/2)
53real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Soil temperature (K)
54real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Average surface pressure [Pa]
55real,    dimension(ngrid,timelen),          intent(in) :: q_co2                              ! Mass mixing ratio of co2 in the first layer (kg/kg)
56real,    dimension(ngrid,timelen),          intent(in) :: q_h2o                              ! Mass mixing ratio of H2o in the first layer (kg/kg)
57
58! Outputs
59real, dimension(ngrid), intent(out) :: delta_mh2oreg ! Difference density of h2o adsorbed (kg/m^3)
60real, dimension(ngrid), intent(out) :: delta_mco2reg ! Difference density of co2 adsorbed (kg/m^3)
61real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil ! Density of co2 adsorbed (kg/m^3)
62real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil ! Density of h2o adsorbed (kg/m^3)
63
64! Local variables
65real, dimension(ngrid,nsoil_PEM,nslope) :: theta_h2o_adsorbed ! Fraction of the pores occupied by H2O molecules
66! -------------
67! Compute H2O adsorption, then CO2 adsorption
68call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
69                            theta_h2o_adsorbed,m_h2o_completesoil,delta_mh2oreg)
70call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, &
71                            tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mco2reg)
72
73END SUBROUTINE regolith_adsorption
74
75!=======================================================================
76SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
77                                  theta_h2o_adsorbed,m_h2o_completesoil,delta_mreg)
78
79!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80!!!
81!!! Purpose: Compute H2O adsorption, following the methods from Jackosky et al., 1997
82!!!
83!!! Author: LL, 01/2023
84!!!
85!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
86
87use soil,         only: layer_PEM, index_breccia
88use comslope_mod,          only: subslope_dist, def_slope_mean
89use vertical_layers_mod,   only: ap, bp
90use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith
91use phys_constants,        only: pi
92
93implicit none
94
95! Inputs
96integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM,timelen   ! Size dimension
97real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! Tendencies on the glaciers ()
98real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
99real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! CO2 ice at the surface [kg/m^2]
100real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Surface pressure (Pa)
101real,    dimension(ngrid,timelen),          intent(in) :: q_co2                              ! Mass mixing ratio of co2 in the first layer (kg/kg)
102real,    dimension(ngrid,timelen),          intent(in) :: q_h2o                              ! Mass mixing ratio of H2o in the first layer (kg/kg)
103real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Soil Thermal inertia (J/K/^2/s^1/2)
104real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Soil temperature (K)
105
106! Outputs
107real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)
108real, dimension(ngrid,nsoil_PEM,nslope), intent(out) :: theta_h2o_adsorbed ! Fraction of the pores occupied by H2O molecules
109real, dimension(ngrid),                  intent(out) :: delta_mreg          ! Difference density of h2o adsorbed (kg/m^3)
110
111! Constants
112real :: Ko =  1.57e-8            ! Jackosky et al. 1997
113real :: e = 2573.9               ! Jackosky et al. 1997
114real :: mu = 0.48                ! Jackosky et al. 1997
115real :: m_theta = 2.84e-7        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
116! real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
117real :: as = 9.48e4              ! Specific area, Zent
118real ::  inertie_thresold = 800. ! TI > 800 means cementation
119
120! Local variables
121real,    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer (kg/m^3)
122real                                           :: K                       ! Used to compute theta
123integer                                        :: ig, iloop, islope, it   ! For loops
124logical, dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
125logical, dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
126real,    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference density of h2o adsorbed per slope (kg/m^3)
127real,    dimension(ngrid,nsoil_PEM,nslope)     :: dm_h2o_regolith_slope   ! Elementary h2o mass adsorded per mesh per slope
128real                                           :: A, B                    ! Used to compute the mean mass above the surface
129real                                           :: p_sat                   ! Saturated vapor pressure of ice
130real,    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
131real,    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
132real,    dimension(:,:), allocatable           :: pvapor                  ! Partial pressure above the surface
133real,    dimension(:)  , allocatable           :: pvapor_avg              ! Yearly averaged
134
135! 0. Some initializations
136allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pvapor(ngrid,timelen),pvapor_avg(ngrid))
137A = 1./m_co2 - 1./m_noco2
138B = 1./m_noco2
139theta_h2o_adsorbed = 0.
140dm_h2o_regolith_slope = 0.
141ispermanent_h2oglaciers = .false.
142ispermanent_co2glaciers = .false.
143
144! 0.1 Look at perennial ice
145do ig = 1,ngrid
146    do islope = 1,nslope
147        if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
148        if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
149    enddo
150enddo
151
152! 0.2 Compute the partial pressure of vapor
153! a. the molecular mass into the column
154do ig = 1,ngrid
155    mass_mean(ig,:) = 1/(A*q_co2(ig,:) + B)
156enddo
157
158! b. pressure level
159do it = 1,timelen
160    do ig = 1,ngrid
161        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
162    enddo
163enddo
164
165! c. Vapor pressure
166pvapor = mass_mean/m_h2o*q_h2o*zplev_mean
167pvapor_avg = sum(pvapor,2)/timelen
168deallocate(pvapor,zplev_mean,mass_mean)
169
170! 1. we compute the mass of H2O adsorded in each layer of the meshes
171do ig = 1,ngrid
172    do islope = 1,nslope
173        do iloop = 1,index_breccia
174            K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
175            if (TI_PEM(ig,iloop,islope) < inertie_thresold) then
176                theta_h2o_adsorbed(ig,iloop,islope) = (K*pvapor_avg(ig)/(1. + K*pvapor_avg(ig)))**mu
177            else
178                p_sat = exp(beta_clap_h2o/tsoil_PEM(ig,iloop,islope) + alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ...
179                theta_h2o_adsorbed(ig,iloop,islope) = (K*p_sat/(1. + K*p_sat))**mu
180            endif
181            dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbed(ig,iloop,islope)*m_theta*rho_regolith
182        enddo
183    enddo
184enddo
185
186! 2. Check the exchange between the atmosphere and the regolith
187do ig = 1,ngrid
188    delta_mreg(ig) = 0.
189    do islope = 1,nslope
190        deltam_reg_slope(ig,islope) = 0.
191        do iloop = 1,index_breccia
192            if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
193                if (iloop == 1) then
194                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop))
195                else
196                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop) - layer_PEM(iloop - 1))
197                endif
198            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
199                deltam_reg_complete(ig,iloop,islope) = 0.
200            endif
201            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
202        enddo
203        delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
204    enddo
205enddo
206m_h2o_completesoil = dm_h2o_regolith_slope
207
208END SUBROUTINE regolith_h2oadsorption
209
210!=======================================================================
211!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
212!!!
213!!! Purpose: Compute CO2  following the methods from Zent & Quinn 1995
214!!!
215!!! Author: LL, 01/2023
216!!!
217!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
218SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg)
219
220use soil,                  only: layer_PEM, index_breccia, index_breccia
221use comslope_mod,          only: subslope_dist, def_slope_mean
222use vertical_layers_mod,   only: ap, bp
223use constants_marspem_mod, only: m_co2, m_noco2, rho_regolith
224use phys_constants,        only: pi
225
226implicit none
227
228! Inputs:
229integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM,timelen   ! Size dimension
230real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Average surface pressure [Pa]
231real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Mean Soil Temperature [K]
232real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Mean Thermal Inertia [USI]
233real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! Tendencies on the glaciers ()
234real,    dimension(ngrid,timelen),          intent(in) :: q_co2, q_h2o                       ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
235real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
236real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! CO2 ice at the surface [kg/m^2]
237
238! Outputs:
239real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil ! Density of co2 adsorbed (kg/m^3)
240real, dimension(ngrid), intent(out) :: delta_mreg ! Difference density of co2 adsorbed (kg/m^3)
241
242! Constants:
243real :: alpha = 7.512e-6        ! Zent & Quinn 1995
244real :: beta = -1541.5          ! Zent & Quinn 1995
245real :: inertie_thresold = 800. ! TI > 800 means cementation
246real :: m_theta = 4.27e-7       ! Mass of co2 per m^2 absorbed
247! real :: as = 18.9e3             ! Specific area, Buhler & Piqueux 2021
248real :: as = 9.48e4             ! Same as previous but from zent
249
250! Local
251real                                           :: A, B                    ! Used to compute the mean mass above the surface
252integer                                        :: ig, islope, iloop, it   ! For loops
253real,    dimension(ngrid,nsoil_PEM,nslope)     :: dm_co2_regolith_slope   ! Elementary mass adsorded per mesh per slope
254logical, dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
255logical, dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
256real,    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer (kg/m^3)
257real,    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference in the mass per slope  (kg/m^3)
258real,    dimension(ngrid,nsoil_PEM,nslope)     :: m_h2o_adsorbed          ! Density of CO2 adsorbed (kg/m^3)
259real,    dimension(ngrid,nsoil_PEM,nslope)     :: theta_h2o_adsorbed      ! Fraction of the pores occupied by H2O molecules
260real,    dimension(ngrid)                      :: delta_mh2o              ! Difference density of h2o adsorbed (kg/m^3)
261! timelen array are allocated because heavy...
262real,    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
263real,    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
264real,    dimension(:,:), allocatable           :: pco2                    ! Partial pressure above the surface
265real,    dimension(:),   allocatable           :: pco2_avg                ! Yearly averaged
266
267! 0. Some initializations
268allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pco2(ngrid,timelen),pco2_avg(ngrid))
269m_h2o_adsorbed = 0.
270A = 1./m_co2 - 1./m_noco2
271B = 1./m_noco2
272dm_co2_regolith_slope = 0.
273delta_mreg = 0.
274ispermanent_h2oglaciers = .false.
275ispermanent_co2glaciers = .false.
276
277! 0.1 Look at perennial ice
278do ig = 1,ngrid
279    do islope = 1,nslope
280        if (abs(d_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
281        if (abs(d_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
282    enddo
283enddo
284
285! 0.2  Compute the partial pressure of CO2
286! a. the molecular mass into the column
287do ig = 1,ngrid
288    mass_mean(ig,:) = 1./(A*q_co2(ig,:) + B)
289enddo
290
291! b. pressure level
292do it = 1,timelen
293    do ig = 1,ngrid
294        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
295    enddo
296enddo
297
298! c. Vapor pressure
299pco2 = mass_mean/m_co2*q_co2*zplev_mean
300pco2_avg(:) = sum(pco2(:,:),2)/timelen
301
302deallocate(zplev_mean,mass_mean,pco2)
303
304! 1. Compute the fraction of the pores occupied by H2O
305call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
306                            theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
307
308! 2. we compute the mass of co2 adsorded in each layer of the meshes
309do ig = 1,ngrid
310    do islope = 1,nslope
311        do iloop = 1,index_breccia
312            if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
313                dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
314                                                         (alpha*pco2_avg(ig) + sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
315            else
316                if (abs(m_co2_completesoil(ig,iloop,islope)) < 1.e-10) then !!! we are at first call
317                    dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) &
318                                                             /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
319                else ! no change: permanent ice stick the atoms of CO2
320                    dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
321                endif
322            endif
323        enddo
324    enddo
325enddo
326
327! 3. Check the exchange between the atmosphere and the regolith
328do ig = 1,ngrid
329    delta_mreg(ig) = 0.
330    do islope = 1,nslope
331        deltam_reg_slope(ig,islope) = 0.
332        do iloop = 1,index_breccia
333            if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
334                if (iloop == 1) then
335                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop))
336                else
337                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop) - layer_PEM(iloop - 1))
338                endif
339            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
340                deltam_reg_complete(ig,iloop,islope) = 0.
341            endif
342            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
343        enddo
344        delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
345    enddo
346enddo
347m_co2_completesoil = dm_co2_regolith_slope
348
349END SUBROUTINE regolith_co2adsorption
350
351END MODULE sorption
Note: See TracBrowser for help on using the repository browser.