source: trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90 @ 3983

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

PEM:
Clearing all the tags and code related to the Generic PCM.
JBC
File size: 18.9 KB
Line 
1MODULE adsorption_mod
2
3implicit none
4
5logical                             :: adsorption_pem     ! True by default, to compute adsorption/desorption. Read in pem.def
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 comsoil_h_PEM,         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 comcstfi_h,            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 comsoil_h_PEM,         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 comcstfi_h,            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 adsorption_mod
352
Note: See TracBrowser for help on using the repository browser.