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

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

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