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

Last change on this file since 3495 was 3456, checked in by jbclement, 6 weeks ago

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