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