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

Last change on this file since 2857 was 2849, checked in by llange, 3 years ago

PEM
Update the codes for subsurface evolution to run with XIOS
1) Water density at the surface and in the soil is now read in the XIOS file
2) Reshape routine created as XIOS output have one element less in the longitude grid
3) The ice table is now computed using these water densities
4) Minor fixs in the main, adsorption module, and tendencies evolutions

LL
File size: 11.3 KB
Line 
1  module adsorption_mod 
2  implicit none
3  contains
4
5  subroutine regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, theta_h2o_adsorbded,m_h2o_adsorbed)
6#ifndef CPP_STD
7 use vertical_layers_mod, ONLY: ap,bp
8 use comsoil_h_PEM, only: n_1km
9#endif
10
11 implicit none
12
13 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! size dimension
14 REAL,INTENT(IN) :: ps(ngrid,timelen)                   ! surface pressure (Pa)         
15 REAL,INTENT(IN) :: q_co2(ngrid,timelen)                ! Mass mixing ratio of co2 in the first layer (kg/kg)
16 REAL,INTENT(IN) :: q_h2o(ngrid,timelen)                ! Mass mixing ratio of H2o in the first layer (kg/kg)
17 REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)      ! Soil Thermal inertia (J/K/^2/s^1/2)
18 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)   ! Soil temperature (K)
19
20! output
21 REAL,INTENT(OUT) ::  m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)     ! Density of h2o adsorbed (kg/m^3)
22 REAL,INTENT(OUT) :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules
23
24! constant
25 REAL :: Ko =  1.57e-8            ! Jackosky et al. 1997
26 REAL :: e = 2573.9               ! Jackosky et al. 1997
27 REAL :: mu = 0.48                ! Jackosky et al. 1997
28 REAL ::  inertie_thresold = 800. ! TI > 800 means cementation
29 real :: m_h2o = 18.01528E-3      ! Molecular weight of h2o (kg/mol)
30 real :: m_co2 = 44.01E-3         ! Molecular weight of co2 (kg/mol)
31 real :: m_noco2 = 33.37E-3       ! Molecular weight of non co2 (kg/mol)
32 REAL ::  rho_regolith = 2000.    ! density of the reoglith, Buhler & Piqueux 2021
33 real :: alpha_clapeyron = -6143.7! eq. (2) in Murphy & Koop 2005
34 real :: beta_clapeyron = 28.9074 ! eq. (2) in Murphy & Koop 2005
35 real :: mi = 2.84e-7             ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
36 real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
37
38! local variable
39 real :: A,B                                                   ! Used to compute the mean mass above the surface
40 real :: K                                                     ! Used to compute theta
41 real :: p_sat                                                 ! saturated vapor pressure of ice
42 integer ig,iloop, islope,isoil,it                             ! for loops
43 real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
44 real,allocatable :: zplev_mean(:,:)                           ! pressure above the surface
45 real,allocatable :: pvapor(:,:)                               ! partial pressure above the surface
46 real, allocatable :: pvapor_avg(:)                            ! yearly average vapor pressure above the surface
47
48! 0. Some initializations
49
50     allocate(mass_mean(ngrid,timelen))
51     allocate(zplev_mean(ngrid,timelen))
52     allocate(pvapor(ngrid,timelen))
53     allocate(pvapor_avg(ngrid))
54
55
56
57      m_h2o_adsorbed(:,:,:) = 0.
58      theta_h2o_adsorbded(:,:,:) = 0.
59      A =(1/m_co2 - 1/m_noco2)
60      B=1/m_noco2
61#ifndef CPP_STD
62! 1. Compute rho surface yearly averaged
63
64!   1.1 Compute the partial pressure of vapor
65!a. the molecular mass into the column
66     do ig = 1,ngrid
67       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
68     enddo
69
70! b. pressure level
71     do it = 1,timelen
72       do ig = 1,ngrid
73         zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
74       enddo
75     enddo
76
77! c. Vapor pressure
78     pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:)
79     pvapor_avg(:) = sum(pvapor(:,:),2)/timelen
80
81     deallocate(pvapor)
82     deallocate(zplev_mean)
83     deallocate(mass_mean)
84
85! 2. we compute the mass of co2 adsorded in each layer of the meshes 
86
87 do ig = 1,ngrid
88  do islope = 1,nslope
89    do iloop = 1,n_1km
90       K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
91     if(TI_PEM(ig,iloop,islope).lt.inertie_thresold)  then
92
93       theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu
94       m_h2o_adsorbed(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
95     else
96        p_sat =exp(alpha_clapeyron/tsoil_PEM(ig,iloop,islope) +beta_clapeyron) ! we assume fixed temperature in the ice ... not really:q good but ...
97        theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu
98        m_h2o_adsorbed(ig,iloop,islope) =as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith
99     endif
100   enddo
101  enddo
102 enddo
103
104 RETURN
105#endif
106  end subroutine
107
108  SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_completesoil,delta_mreg)
109#ifndef CPP_STD
110      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km
111      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
112      use comslope_mod, only : subslope_dist,def_slope_mean
113      use vertical_layers_mod, ONLY: ap,bp
114#endif
115
116      IMPLICIT NONE
117! Input: 
118 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen             ! size dimension
119 REAL,INTENT(IN) :: ps(ngrid,timelen)                               ! Average surface pressure [Pa]
120 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)               ! Mean Soil Temperature [K]
121 REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)                  ! Mean Thermal Inertia [USI]
122 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) !tendancies on the glaciers ()
123 REAL,INTENT(IN) :: q_co2(ngrid,timelen),q_h2o(ngrid,timelen)       ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
124 REAL,INTENT(IN) :: waterice(ngrid,nslope)                          ! water ice at the surface [kg/m^2]
125 REAL,INTENT(IN) :: co2ice(ngrid,nslope)                            ! co2 ice at the surface [kg/m^2]
126
127! Outputs:
128 REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)   ! Density of co2 adsorbed (kg/m^3)
129 REAL,INTENT(INOUT) :: delta_mreg(ngrid)                            ! Difference density of co2 adsorbed (kg/m^3)
130 
131! Constants:
132
133 REAL :: alpha = 7.512e-6 ! Zent & Quinn 1995
134 REAL :: beta =  -1541.5  ! Zent & Quinn 1995
135 REAL ::  inertie_thresold = 800. ! TI > 800 means cementation
136 REAL ::  rho_regolith = 2000. ! density of the reoglith, buhler & piqueux 2021
137 real :: m_co2 = 44.01E-3      ! Molecular weight of co2 (kg/mol)
138 real :: m_noco2 = 33.37E-3    ! Molecular weight of h2o (kg/mol)
139 real :: m_theta = 4.27e-7     ! Mass of co2 per m^2 absorbed
140 real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
141
142! Local         
143 real :: A,B                                             ! Used to compute the mean mass above the surface
144 INTEGER :: ig,islope,iloop,it                           ! for loops
145 REAL :: dm_co2_regolith_slope(ngrid,nsoil_PEM,nslope)   ! elementary mass adsorded per mesh per slope
146 INTEGER :: ispermanent_co2glaciers(ngrid,nslope)        ! Check if the glacier is permanent
147 INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)        ! Check if the glacier is permanent
148#ifndef CPP_STD
149 REAL :: deltam_reg_complete(ngrid,n_1km,nslope)         ! Difference in the mass per slope and soil layer (kg/m^3)
150#endif
151 REAL :: deltam_reg_slope(ngrid,nslope)                  ! Difference in the mass per slope  (kg/m^3)
152 REAL :: m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)          ! Density of CO2 adsorbed (kg/m^3)
153 REAL :: theta_h2o_adsorbed(ngrid,nsoil_PEM,nslope)     ! Fraction of the pores occupied by H2O molecules
154!timelen array are allocated because heavy ...
155 real,allocatable :: mass_mean(:,:)                            ! mean mass above the surface
156 real,allocatable :: zplev_mean(:,:)                           ! pressure above the surface
157 real,allocatable :: pco2(:,:)                                  ! partial pressure above the surface
158 real, allocatable :: pco2_avg(:)                              ! yearly averaged
159
160! 0. Some initializations
161
162     allocate(mass_mean(ngrid,timelen))
163     allocate(zplev_mean(ngrid,timelen))
164     allocate(pco2(ngrid,timelen))
165     allocate(pco2_avg(ngrid))
166
167#ifndef CPP_STD
168 
169      m_h2o_adsorbed(:,:,:) = 0.
170      A =(1/m_co2 - 1/m_noco2)
171      B=1/m_noco2
172
173     dm_co2_regolith_slope(:,:,:) = 0
174     delta_mreg(:) = 0.
175
176!0.1 Look at perenial ice
177  do ig = 1,ngrid
178    do islope = 1,nslope
179     if((abs(tend_h2oglaciers(ig,islope)).gt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then
180        ispermanent_h2oglaciers(ig,islope) = 1
181     else
182        ispermanent_h2oglaciers(ig,islope) = 0
183     endif
184
185     if((abs(tend_co2glaciers(ig,islope)).gt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
186        ispermanent_co2glaciers(ig,islope) = 1
187     else
188        ispermanent_co2glaciers(ig,islope) = 0
189     endif
190    enddo
191   enddo
192
193!   0.2  Compute the partial pressure of CO2
194!a. the molecular mass into the column
195     do ig = 1,ngrid
196       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
197     enddo
198
199! b. pressure level
200     do it = 1,timelen
201       do ig = 1,ngrid
202         zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
203       enddo
204     enddo
205
206! c. Vapor pressure
207     pco2(:,:) = mass_mean(:,:)/m_co2*q_co2(:,:)*zplev_mean(:,:)
208     pco2_avg(:) = sum(pco2(:,:),2)/timelen
209
210     deallocate(zplev_mean)
211     deallocate(mass_mean)
212     deallocate(pco2)
213
214
215! 1. Compute the fraction of the pores occupied by H2O
216 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen, ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,theta_h2o_adsorbed, m_h2o_adsorbed)
217
218! 2.  we compute the mass of co2 adsorded in each layer of the meshes 
219
220 do ig = 1,ngrid
221  do islope = 1,nslope
222    do iloop = 1,n_1km
223     if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
224     dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
225                                             (alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
226     else
227        if(abs(m_co2_completesoil(ig,iloop,islope)).lt.1-10) then !!! we are at first call
228          dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) &
229                                                 /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
230        else ! no change: permanent ice stick the atoms of CO2
231          dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
232        endif
233     endif
234  enddo
235 enddo
236enddo
237
238! 2. Check the exchange between the atmosphere and the regolith
239
240  do ig = 1,ngrid
241   delta_mreg(ig) = 0.
242   do islope = 1,nslope
243    deltam_reg_slope(ig,islope) = 0.
244    do iloop = 1,n_1km
245       if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
246             deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) &
247                                                   *(layer_PEM(iloop+1) - layer_PEM(iloop))
248       else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
249             deltam_reg_complete(ig,iloop,islope) = 0.
250       endif
251       deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
252    enddo
253   delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
254   enddo
255  enddo
256   m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:)
257
258!=======================================================================
259      RETURN
260#endif
261      END
262
263
264   end module
Note: See TracBrowser for help on using the repository browser.