Ignore:
Timestamp:
Feb 12, 2026, 9:09:12 AM (3 weeks ago)
Author:
jbclement
Message:

PEM:
Major refactor following the previous ones (r3989 and r3991) completing the large structural reorganization and cleanup of the PEM codebase. This revision introduces newly designed modules, standardizes interfaces with explicit ini/end APIs and adds native NetCDF I/O together with explicit PCM/PEM adapters. In detail:

  • Some PEM models were corrected or improved:
    • Frost/perennial ice semantics are clarified via renaming;
    • Soil temperature remapping clarified, notably by removing the rescaling of temperature deviation;
    • Geothermal flux for the PCM is computed based on the PEM state;
  • New explicit PEM/PCM adapters ("set_*"/"build4PCM_*") to decouple PEM internal representation from PCM file layouts and reconstruct consistent fields returned to the PCM;
  • New explicit build/teardown routines that centralize allocation and initialization ordering, reducing accidental use of uninitialized data and making the model lifecycle explicit;
  • Add native read/write helpers for NetCDF that centralize all low-level NetCDF interactions with major improvements (and more simplicity) compared to legacy PEM/PCM I/O (see the modules "io_netcdf" and "output"). They support reading, creation and writing of "diagevol.nc" (renamed from "diagpem.nc") and start/restart files;
  • Provide well-focused modules ("numerics"/"maths"/"utility"/"display") to host commonly-used primitives:
    • "numerics" defines numerical types and constants for reproducibility, portability across compilers and future transitions (e.g. quadruple precision experiments);
    • "display" provides a single controlled interface for runtime messages, status output and diagnostics, avoiding direct 'print'/'write' to enable silent mode, log redirection, and MPI-safe output in the future.
    • "utility" (new module) hosts generic helpers used throughout the code (e.g. "int2str" or "real2str");
  • Add modules "clim_state_init"/"clim_state_rec" which provide robust read/write logic for "start/startfi/startpem", including 1D fallbacks, mesh conversions and dimension checks. NetCDF file creation is centralized and explicit. Restart files are now self-consistent and future-proof, requiring changes only to affected variables;
  • Add module "atmosphere" which computes pressure fields, reconstructs potential temperature and air mass. It also holds the whole logic to define sigma or hybrid coordinates for altitudes;
  • Add module "geometry" to centrilize dimensions logic and grid conversions routines (including 2 new ones "dyngrd2vect"/"vect2dyngrd");
  • Add module "slopes" to isolate slopes handling;
  • Add module "surface" to isolate surface management. Notably, albedo and emissivity are now fully reconstructed following the PCM settings;
  • Add module "allocation" to check the dimension initialization and centrilize allocation/deallocation;
  • Finalize module decomposition and renaming to consolidate domain-based modules, purpose-based routines and physics/process-based variables;
  • The main program now drives a clearer sequence of conceptual steps (initialization / reading / evolution / update / build / writing) and fails explicitly instead of silently defaulting;
  • Ice table logic is made restart-safe;
  • 'Open'/'read' intrinsic logic is made safe and automatic;
  • Improve discoverability and standardize the data handling (private vs protected vs public);
  • Apply consistent documentation/header style already introduced;
  • Update deftank scripts to reflect new names and launch wrappers;

This revision is a structural milestone aiming to be behavior-preserving where possible. It has been tested via compilation and short integration runs. However, due to extensive renames, moves, and API changes, full validation is still ongoing.
Note: the revision includes one (possibly two) easter egg hidden in the code for future archaeologists of the PEM. No physics were harmed.
JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/sorption.F90

    r3991 r4065  
    1010! AUTHORS & DATE
    1111!     L. Lange, 2023
    12 !     JB Clement, 2025
    13 !
    14 ! NOTES
    15 !
    16 !-----------------------------------------------------------------------
    17 
    18 ! DECLARATION
    19 ! -----------
    20 implicit none
    21 
    22 ! MODULE VARIABLES
    23 ! ----------------
    24 logical                             :: do_sorption       ! Flag to compute adsorption/desorption
    25 real, allocatable, dimension(:,:,:) :: co2_adsorbed_phys ! co2 that is in the regolith [kg/m^2]
    26 real, allocatable, dimension(:,:,:) :: h2o_adsorbed_phys ! h2o that is in the regolith [kg/m^2]
     12!     JB Clement, 02/2026
     13!
     14! NOTES
     15!
     16!-----------------------------------------------------------------------
     17
     18! DEPENDENCIES
     19! ------------
     20use numerics, only: dp, qp, di, k4, minieps
     21
     22! DECLARATION
     23! -----------
     24implicit none
     25
     26! PARAMETERS
     27! ----------
     28logical(k4), protected :: do_sorption ! Flag to compute adsorption/desorption
    2729
    2830contains
    … …  
    3032
    3133!=======================================================================
    32 SUBROUTINE ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
    33 !-----------------------------------------------------------------------
    34 ! NAME
    35 !     ini_adsorption_h_PEM
    36 !
    37 ! DESCRIPTION
    38 !     Allocate CO2 and H2O adsorption arrays.
     34SUBROUTINE set_sorption_config(do_sorption_in)
     35!-----------------------------------------------------------------------
     36! NAME
     37!     set_sorption_config
     38!
     39! DESCRIPTION
     40!     Setter for 'sorption' configuration parameters.
     41!
     42! AUTHORS & DATE
     43!     JB Clement, 02/2026
     44!
     45! NOTES
     46!
     47!-----------------------------------------------------------------------
     48
     49! DEPENDENCIES
     50! ------------
     51use utility, only: bool2str
     52use display, only: print_msg
     53
     54! DECLARATION
     55! -----------
     56implicit none
     57
     58! ARGUMENTS
     59! ---------
     60logical(k4), intent(in) :: do_sorption_in
     61
     62! CODE
     63! ----
     64do_sorption = do_sorption_in
     65call print_msg('do_sorption = '//bool2str(do_sorption))
     66
     67END SUBROUTINE set_sorption_config
     68!=======================================================================
     69
     70!=======================================================================
     71SUBROUTINE evolve_regolith_adsorption(d_h2oice,d_co2ice,h2oice,co2ice,tsoil,TI,ps,q_h2o_ts,q_co2_ts,h2o_ads_reg,co2_ads_reg,delta_h2o_ads,delta_co2_ads)
     72!-----------------------------------------------------------------------
     73! NAME
     74!     evolve_regolith_adsorption
     75!
     76! DESCRIPTION
     77!     Compute both H2O and CO2 adsorption in regolith.
     78!
     79! AUTHORS & DATE
     80!     L. Lange, 2023
     81!     JB Clement, 12/2025
     82!     C. Metz, 02/2026
     83!
     84! NOTES
     85!
     86!-----------------------------------------------------------------------
     87
     88! DEPENDENCIES
     89! ------------
     90use geometry, only: ngrid, nslope, nsoil
     91
     92! DECLARATION
     93! -----------
     94implicit none
     95
     96! ARGUMENTS
     97! ---------
     98real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! tendencies on the glaciers [1]
     99real(dp), dimension(:,:),   intent(in)    :: h2oice             ! H2O ice at the surface [kg/m^2]
     100real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
     101real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Soil Thermal inertia [J/K/^2/s^1/2]
     102real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
     103real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
     104real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of co2 in the first layer [kg/kg]
     105real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2o in the first layer [kg/kg]
     106real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of h2o adsorbed [kg/m^3]
     107real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of co2 adsorbed [kg/m^3]
     108real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg, co2_ads_reg
     109
     110! LOCAL VARIABLES
     111! ---------------
     112real(dp), dimension(ngrid,nsoil,nslope) :: theta_h2o_ads ! Fraction of the pores occupied by H2O molecules
     113
     114! CODE
     115! ----
     116call evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,theta_h2o_ads,h2o_ads_reg,delta_h2o_ads)
     117call evolve_co2_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads)
     118
     119END SUBROUTINE evolve_regolith_adsorption
     120!=======================================================================
     121
     122!=======================================================================
     123SUBROUTINE evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,theta_h2o_ads,h2o_ads_reg,delta_h2o_ads)
     124!-----------------------------------------------------------------------
     125! NAME
     126!     evolve_h2o_ads
     127!
     128! DESCRIPTION
     129!     Compute H2O adsorption in regolith following Jackosky et al. (1997).
    39130!
    40131! AUTHORS & DATE
    … …  
    46137!-----------------------------------------------------------------------
    47138
     139! DEPENDENCIES
     140! ------------
     141use geometry,   only: ngrid, nslope, nsoil, nday
     142use soil,       only: layer, index_breccia
     143use slopes,     only: subslope_dist, def_slope_mean
     144use atmosphere, only: ap, bp
     145use planet,     only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith
     146use maths,      only: pi
     147
    48148! DECLARATION
    49149! -----------
    … …  
    52152! ARGUMENTS
    53153! ---------
    54 integer, intent(in) :: ngrid       ! number of atmospheric columns
    55 integer, intent(in) :: nslope      ! number of slope within a mesh
    56 integer, intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM
    57 
    58 ! CODE
    59 ! ----
    60 allocate(co2_adsorbed_phys(ngrid,nsoilmx_PEM,nslope))
    61 allocate(h2o_adsorbed_phys(ngrid,nsoilmx_PEM,nslope))
    62 
    63 END SUBROUTINE ini_adsorption_h_PEM
    64 !=======================================================================
    65 
    66 !=======================================================================
    67 SUBROUTINE end_adsorption_h_PEM
    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 !-----------------------------------------------------------------------
    82 
    83 ! DECLARATION
    84 ! -----------
    85 implicit none
    86 
    87 ! CODE
    88 ! ----
    89 if (allocated(co2_adsorbed_phys)) deallocate(co2_adsorbed_phys)
    90 if (allocated(h2o_adsorbed_phys)) deallocate(h2o_adsorbed_phys)
    91 
    92 END SUBROUTINE end_adsorption_h_PEM
    93 !=======================================================================
    94 
    95 !=======================================================================
    96 SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
    97                                m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg)
    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 !-----------------------------------------------------------------------
    112 
    113 ! DECLARATION
    114 ! -----------
    115 implicit none
    116 
    117 ! ARGUMENTS
    118 ! ---------
    119 integer,                                    intent(in)    :: ngrid, nslope, nsoil_PEM, timelen  ! size dimension: physics x subslope x soil x timeseries
    120 real,    dimension(ngrid,nslope),           intent(in)    :: d_h2oglaciers, d_co2glaciers ! tendancies on the glaciers [1]
    121 real,    dimension(ngrid,nslope),           intent(in)    :: waterice                           ! water ice at the surface [kg/m^2]
    122 real,    dimension(ngrid,nslope),           intent(in)    :: co2ice                             ! co2 ice at the surface [kg/m^2]
    123 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: TI_PEM                             ! Soil Thermal inertia (J/K/^2/s^1/2)
    124 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: tsoil_PEM                          ! Soil temperature (K)
    125 real,    dimension(ngrid,timelen),          intent(in)    :: ps                                 ! Average surface pressure [Pa]
    126 real,    dimension(ngrid,timelen),          intent(in)    :: q_co2                              ! Mass mixing ratio of co2 in the first layer (kg/kg)
    127 real,    dimension(ngrid,timelen),          intent(in)    :: q_h2o                              ! Mass mixing ratio of H2o in the first layer (kg/kg)
    128 real,    dimension(ngrid),                  intent(out)   :: delta_mh2oreg                      ! Difference density of h2o adsorbed (kg/m^3)
    129 real,    dimension(ngrid),                  intent(out)   :: delta_mco2reg                      ! Difference density of co2 adsorbed (kg/m^3)
    130 real,    dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil                 ! Density of co2 adsorbed (kg/m^3)
    131 real,    dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil                 ! Density of h2o adsorbed (kg/m^3)
    132 
    133 ! LOCAL VARIABLES
    134 ! ---------------
    135 real, dimension(ngrid,nsoil_PEM,nslope) :: theta_h2o_adsorbed ! Fraction of the pores occupied by H2O molecules
    136 
    137 ! CODE
    138 ! ----
    139 call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
    140                             theta_h2o_adsorbed,m_h2o_completesoil,delta_mh2oreg)
    141 call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o, &
    142                             tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mco2reg)
    143 
    144 END SUBROUTINE regolith_adsorption
    145 !=======================================================================
    146 
    147 !=======================================================================
    148 SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
    149                                   theta_h2o_adsorbed,m_h2o_completesoil,delta_mreg)
    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 !-----------------------------------------------------------------------
    164 
    165 ! DEPENDENCIES
    166 ! ------------
    167 use soil,                  only: layer_PEM, index_breccia
    168 use comslope_mod,          only: subslope_dist, def_slope_mean
    169 use vertical_layers_mod,   only: ap, bp
    170 use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith
    171 use phys_constants,        only: pi
    172 
    173 ! DECLARATION
    174 ! -----------
    175 implicit none
    176 
    177 ! ARGUMENTS
    178 ! ---------
    179 integer,                                    intent(in)    :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension
    180 real,    dimension(ngrid,nslope),           intent(in)    :: d_h2oglaciers, d_co2glaciers     ! Tendencies on the glaciers ()
    181 real,    dimension(ngrid,nslope),           intent(in)    :: waterice                         ! Water ice at the surface [kg/m^2]
    182 real,    dimension(ngrid,nslope),           intent(in)    :: co2ice                           ! CO2 ice at the surface [kg/m^2]
    183 real,    dimension(ngrid,timelen),          intent(in)    :: ps                               ! Surface pressure (Pa)
    184 real,    dimension(ngrid,timelen),          intent(in)    :: q_co2                            ! Mass mixing ratio of co2 in the first layer (kg/kg)
    185 real,    dimension(ngrid,timelen),          intent(in)    :: q_h2o                            ! Mass mixing ratio of H2o in the first layer (kg/kg)
    186 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: TI_PEM                           ! Soil Thermal inertia (J/K/^2/s^1/2)
    187 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: tsoil_PEM                        ! Soil temperature (K)
    188 real,    dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil               ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)
    189 real,    dimension(ngrid,nsoil_PEM,nslope), intent(out)   :: theta_h2o_adsorbed               ! Fraction of the pores occupied by H2O molecules
    190 real,    dimension(ngrid),                  intent(out)   :: delta_mreg                       ! Difference density of h2o adsorbed (kg/m^3)
     154real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! Tendencies on the glaciers ()
     155real(dp), dimension(:,:),   intent(in)    :: h2oice             ! H2O ice at the surface [kg/m^2]
     156real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
     157real(dp), dimension(:,:),   intent(in)    :: ps                 ! Surface pressure [Pa]
     158real(dp), dimension(:,:),   intent(in)    :: q_co2_ts           ! Mass mixing ratio of co2 in the first layer [kg/kg]
     159real(dp), dimension(:,:),   intent(in)    :: q_h2o_ts           ! Mass mixing ratio of H2o in the first layer [kg/kg]
     160real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Soil Thermal inertia [J/K/^2/s^1/2]
     161real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Soil temperature [K]
     162real(dp), dimension(:,:,:), intent(inout) :: h2o_ads_reg        ! Density of h2o adsorbed [kg/m^3]
     163real(dp), dimension(:,:,:), intent(out)   :: theta_h2o_ads      ! Fraction of the pores occupied by H2O molecules
     164real(dp), dimension(:),     intent(out)   :: delta_h2o_ads      ! Difference density of h2o adsorbed [kg/m^3]
    191165
    192166! LOCAL VARIABLES
    193167! ---------------
    194168! Constants:
    195 real :: Ko =  1.57e-8            ! Jackosky et al. 1997
    196 real :: e = 2573.9               ! Jackosky et al. 1997
    197 real :: mu = 0.48                ! Jackosky et al. 1997
    198 real :: 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
    200 real :: as = 9.48e4              ! Specific area, Zent
    201 real ::  inertie_thresold = 800. ! TI > 800 means cementation
    202 real,    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer (kg/m^3)
    203 real                                           :: K                       ! Used to compute theta
    204 integer                                        :: ig, iloop, islope, it   ! For loops
    205 logical, dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
    206 logical, dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
    207 real,    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference density of h2o adsorbed per slope (kg/m^3)
    208 real,    dimension(ngrid,nsoil_PEM,nslope)     :: dm_h2o_regolith_slope   ! Elementary h2o mass adsorded per mesh per slope
    209 real                                           :: A, B                    ! Used to compute the mean mass above the surface
    210 real                                           :: p_sat                   ! Saturated vapor pressure of ice
    211 real,    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
    212 real,    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
    213 real,    dimension(:,:), allocatable           :: pvapor                  ! Partial pressure above the surface
    214 real,    dimension(:)  , allocatable           :: pvapor_avg              ! Yearly averaged
     169real(dp) :: Ko =  1.57e-8_dp            ! Jackosky et al. 1997
     170real(dp) :: e = 2573.9_dp               ! Jackosky et al. 1997
     171real(dp) :: mu = 0.48_dp                ! Jackosky et al. 1997
     172real(dp) :: m_theta = 2.84e-7_dp        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
     173! real(dp) :: as = 18.9e3_dp              ! Specific area, Buhler & Piqueux 2021
     174real(dp) :: as = 9.48e4_dp              ! Specific area, Zent
     175real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
     176real(dp), dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
     177real(dp)                                        :: K                       ! Used to compute theta
     178integer(di)                                     :: ig, iloop, islope, it   ! For loops
     179logical(k4), dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
     180logical(k4), dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
     181real(dp), dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference density of h2o adsorbed per slope [kg/m^3]
     182real(dp), dimension(ngrid,nsoil,nslope)         :: dm_h2o_regolith_slope   ! Elementary h2o mass adsorded per mesh per slope
     183real(dp)                                        :: A, B                    ! Used to compute the mean mass above the surface
     184real(dp)                                        :: p_sat                   ! Saturated vapor pressure of ice
     185real(dp), dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
     186real(dp), dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
     187real(dp), dimension(:,:), allocatable           :: pvapor                  ! Partial pressure above the surface
     188real(dp), dimension(:)  , allocatable           :: pvapor_avg              ! Yearly averaged
    215189
    216190! CODE
    217191! ----
    218192! 0. Some initializations
    219 allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pvapor(ngrid,timelen),pvapor_avg(ngrid))
    220 A = 1./m_co2 - 1./m_noco2
    221 B = 1./m_noco2
    222 theta_h2o_adsorbed = 0.
    223 dm_h2o_regolith_slope = 0.
     193allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pvapor(ngrid,nday),pvapor_avg(ngrid))
     194A = 1._dp/m_co2 - 1._dp/m_noco2
     195B = 1._dp/m_noco2
     196theta_h2o_ads = 0._dp
     197dm_h2o_regolith_slope = 0._dp
    224198ispermanent_h2oglaciers = .false.
    225199ispermanent_co2glaciers = .false.
    226200
    227201! 0.1 Look at perennial ice
    228 do 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.
    232     enddo
    233 enddo
     202where (abs(d_h2oice(:,:)) > 1.e-5 .and. h2oice(:,:) > 0._dp) ispermanent_h2oglaciers(:,:) = .true.
     203where (abs(d_co2ice(:,:)) > 1.e-5 .and. co2ice(:,:) > 0._dp) ispermanent_co2glaciers(:,:) = .true.
    234204
    235205! 0.2 Compute the partial pressure of vapor
    236206! a. the molecular mass into the column
    237207do ig = 1,ngrid
    238     mass_mean(ig,:) = 1/(A*q_co2(ig,:) + B)
    239 enddo
     208    mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B)
     209end do
    240210
    241211! b. pressure level
    242 do it = 1,timelen
     212do it = 1,nday
    243213    do ig = 1,ngrid
    244214        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
    245     enddo
    246 enddo
     215    end do
     216end do
    247217
    248218! c. Vapor pressure
    249 pvapor = mass_mean/m_h2o*q_h2o*zplev_mean
    250 pvapor_avg = sum(pvapor,2)/timelen
     219pvapor = mass_mean/m_h2o*q_h2o_ts*zplev_mean
     220pvapor_avg = sum(pvapor,2)/nday
    251221deallocate(pvapor,zplev_mean,mass_mean)
    252222
    … …  
    255225    do islope = 1,nslope
    256226        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
     227            K = Ko*exp(e/tsoil(ig,iloop,islope))
     228            if (TI(ig,iloop,islope) < inertie_thresold) then
     229                theta_h2o_ads(ig,iloop,islope) = (K*pvapor_avg(ig)/(1._dp + K*pvapor_avg(ig)))**mu
    260230            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
    265         enddo
    266     enddo
    267 enddo
     231                p_sat = exp(beta_clap_h2o/tsoil(ig,iloop,islope) + alpha_clap_h2o) ! we assume fixed temperature in the ice ... not really good but ...
     232                theta_h2o_ads(ig,iloop,islope) = (K*p_sat/(1._dp + K*p_sat))**mu
     233            end if
     234            dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_ads(ig,iloop,islope)*m_theta*rho_regolith
     235        end do
     236    end do
     237end do
    268238
    269239! 2. Check the exchange between the atmosphere and the regolith
    270240do ig = 1,ngrid
    271     delta_mreg(ig) = 0.
     241    delta_h2o_ads(ig) = 0._dp
    272242    do islope = 1,nslope
    273         deltam_reg_slope(ig,islope) = 0.
     243        deltam_reg_slope(ig,islope) = 0._dp
    274244        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
     245            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
    276246                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))
     247                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop))
    278248                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))
    280                 endif
     249                    deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - h2o_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1))
     250                end if
    281251            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
    282                 deltam_reg_complete(ig,iloop,islope) = 0.
    283             endif
     252                deltam_reg_complete(ig,iloop,islope) = 0._dp
     253            end if
    284254            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
    285         enddo
    286         delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    287     enddo
    288 enddo
    289 m_h2o_completesoil = dm_h2o_regolith_slope
    290 
    291 END SUBROUTINE regolith_h2oadsorption
    292 !=======================================================================
    293 
    294 !=======================================================================
    295 SUBROUTINE 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)
    296 
    297 !-----------------------------------------------------------------------
    298 ! NAME
    299 !     regolith_co2adsorption
     255        end do
     256        delta_h2o_ads(ig) = delta_h2o_ads(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180._dp)
     257    end do
     258end do
     259h2o_ads_reg = dm_h2o_regolith_slope
     260
     261END SUBROUTINE evolve_h2o_ads
     262!=======================================================================
     263
     264!=======================================================================
     265SUBROUTINE evolve_co2_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_h2o_ts,q_co2_ts,tsoil,TI,co2_ads_reg,delta_co2_ads)
     266
     267!-----------------------------------------------------------------------
     268! NAME
     269!     evolve_co2_ads
    300270!
    301271! DESCRIPTION
    … …  
    312282! DEPENDENCIES
    313283! ------------
    314 use soil,                  only: layer_PEM, index_breccia
    315 use comslope_mod,          only: subslope_dist, def_slope_mean
    316 use vertical_layers_mod,   only: ap, bp
    317 use constants_marspem_mod, only: m_co2, m_noco2, rho_regolith
    318 use phys_constants,        only: pi
     284use geometry,   only: ngrid, nslope, nsoil, nday
     285use soil,       only: layer, index_breccia
     286use slopes,     only: subslope_dist, def_slope_mean
     287use atmosphere, only: ap, bp
     288use planet,     only: m_co2, m_noco2, rho_regolith
     289use maths,      only: pi
    319290
    320291! DECLARATION
    … …  
    324295! ARGUMENTS
    325296! ---------
    326 integer,                                    intent(in)    :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension
    327 real,    dimension(ngrid,timelen),          intent(in)    :: ps                               ! Average surface pressure [Pa]
    328 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: tsoil_PEM                        ! Mean Soil Temperature [K]
    329 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: TI_PEM                           ! Mean Thermal Inertia [USI]
    330 real,    dimension(ngrid,nslope),           intent(in)    :: d_h2oglaciers, d_co2glaciers     ! Tendencies on the glaciers ()
    331 real,    dimension(ngrid,timelen),          intent(in)    :: q_co2, q_h2o                     ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
    332 real,    dimension(ngrid,nslope),           intent(in)    :: waterice                         ! Water ice at the surface [kg/m^2]
    333 real,    dimension(ngrid,nslope),           intent(in)    :: co2ice                           ! CO2 ice at the surface [kg/m^2]
    334 real,    dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil               ! Density of co2 adsorbed (kg/m^3)
    335 real,    dimension(ngrid),                  intent(out)   :: delta_mreg                       ! Difference density of co2 adsorbed (kg/m^3)
     297real(dp), dimension(:,:),   intent(in)    :: ps                 ! Average surface pressure [Pa]
     298real(dp), dimension(:,:,:), intent(in)    :: tsoil              ! Mean Soil Temperature [K]
     299real(dp), dimension(:,:,:), intent(in)    :: TI                 ! Mean Thermal Inertia [USI]
     300real(dp), dimension(:,:),   intent(in)    :: d_h2oice, d_co2ice ! Tendencies on the glaciers ()
     301real(dp), dimension(:,:),   intent(in)    :: q_co2_ts, q_h2o_ts ! Mass mixing ratio of co2 and h2o in the first layer [kg/kg]
     302real(dp), dimension(:,:),   intent(in)    :: h2oice             ! Water ice at the surface [kg/m^2]
     303real(dp), dimension(:,:),   intent(in)    :: co2ice             ! CO2 ice at the surface [kg/m^2]
     304real(dp), dimension(:,:,:), intent(inout) :: co2_ads_reg        ! Density of co2 adsorbed [kg/m^3]
     305real(dp), dimension(:),     intent(out)   :: delta_co2_ads      ! Difference density of co2 adsorbed [kg/m^3]
    336306
    337307! LOCAL VARIABLES
    338308! ---------------
    339309! Constants:
    340 real :: alpha = 7.512e-6        ! Zent & Quinn 1995
    341 real :: beta = -1541.5          ! Zent & Quinn 1995
    342 real :: inertie_thresold = 800. ! TI > 800 means cementation
    343 real :: m_theta = 4.27e-7       ! Mass of co2 per m^2 absorbed
    344 ! real :: as = 18.9e3             ! Specific area, Buhler & Piqueux 2021
    345 real :: as = 9.48e4             ! Same as previous but from zent
    346 real                                           :: A, B                    ! Used to compute the mean mass above the surface
    347 integer                                        :: ig, islope, iloop, it   ! For loops
    348 real,    dimension(ngrid,nsoil_PEM,nslope)     :: dm_co2_regolith_slope   ! Elementary mass adsorded per mesh per slope
    349 logical, dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
    350 logical, dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
    351 real,    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer (kg/m^3)
    352 real,    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference in the mass per slope  (kg/m^3)
    353 real,    dimension(ngrid,nsoil_PEM,nslope)     :: m_h2o_adsorbed          ! Density of CO2 adsorbed (kg/m^3)
    354 real,    dimension(ngrid,nsoil_PEM,nslope)     :: theta_h2o_adsorbed      ! Fraction of the pores occupied by H2O molecules
    355 real,    dimension(ngrid)                      :: delta_mh2o              ! Difference density of h2o adsorbed (kg/m^3)
    356 ! timelen array are allocated because heavy...
    357 real,    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
    358 real,    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
    359 real,    dimension(:,:), allocatable           :: pco2                    ! Partial pressure above the surface
    360 real,    dimension(:),   allocatable           :: pco2_avg                ! Yearly averaged
     310real(dp) :: alpha = 7.512e-6_dp        ! Zent & Quinn 1995
     311real(dp) :: beta = -1541.5_dp          ! Zent & Quinn 1995
     312real(dp) :: inertie_thresold = 800._dp ! TI > 800 means cementation
     313real(dp) :: m_theta = 4.27e-7_dp       ! Mass of co2 per m^2 absorbed
     314! real(dp) :: as = 18.9e3_dp           ! Specific area, Buhler & Piqueux 2021
     315real(dp) :: as = 9.48e4_dp             ! Same as previous but from zent
     316real(dp)                                           :: A, B                    ! Used to compute the mean mass above the surface
     317integer(di)                                        :: ig, islope, iloop, it   ! For loops
     318real(dp),    dimension(ngrid,nsoil,nslope)         :: dm_co2_regolith_slope   ! Elementary mass adsorded per mesh per slope
     319logical(k4), dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
     320logical(k4), dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
     321real(dp),    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer [kg/m^3]
     322real(dp),    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference in the mass per slope  [kg/m^3]
     323real(dp),    dimension(ngrid,nsoil,nslope)         :: m_h2o_adsorbed          ! Density of CO2 adsorbed [kg/m^3]
     324real(dp),    dimension(ngrid,nsoil,nslope)         :: theta_h2o_ads           ! Fraction of the pores occupied by H2O molecules
     325real(dp),    dimension(ngrid)                      :: delta_mh2o              ! Difference density of h2o adsorbed [kg/m^3]
     326! nday array are allocated because heavy...
     327real(dp),    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
     328real(dp),    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
     329real(dp),    dimension(:,:), allocatable           :: pco2                    ! Partial pressure above the surface
     330real(dp),    dimension(:),   allocatable           :: pco2_avg                ! Yearly averaged
    361331
    362332! CODE
    363333! ----
    364334! 0. Some initializations
    365 allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pco2(ngrid,timelen),pco2_avg(ngrid))
    366 m_h2o_adsorbed = 0.
    367 A = 1./m_co2 - 1./m_noco2
    368 B = 1./m_noco2
    369 dm_co2_regolith_slope = 0.
    370 delta_mreg = 0.
     335allocate(mass_mean(ngrid,nday),zplev_mean(ngrid,nday),pco2(ngrid,nday),pco2_avg(ngrid))
     336m_h2o_adsorbed = 0._dp
     337A = 1._dp/m_co2 - 1._dp/m_noco2
     338B = 1._dp/m_noco2
     339dm_co2_regolith_slope = 0._dp
     340delta_co2_ads = 0._dp
    371341ispermanent_h2oglaciers = .false.
    372342ispermanent_co2glaciers = .false.
    373343
    374344! 0.1 Look at perennial ice
    375 do 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.
    379     enddo
    380 enddo
     345where (abs(d_h2oice(:,:)) > 1.e-5 .and. h2oice(:,:) > 0._dp) ispermanent_h2oglaciers(:,:) = .true.
     346where (abs(d_co2ice(:,:)) > 1.e-5 .and. co2ice(:,:) > 0._dp) ispermanent_co2glaciers(:,:) = .true.
    381347
    382348! 0.2  Compute the partial pressure of CO2
    383349! a. the molecular mass into the column
    384350do ig = 1,ngrid
    385     mass_mean(ig,:) = 1./(A*q_co2(ig,:) + B)
    386 enddo
     351    mass_mean(ig,:) = 1._dp/(A*q_co2_ts(ig,:) + B)
     352end do
    387353
    388354! b. pressure level
    389 do it = 1,timelen
     355do it = 1,nday
    390356    do ig = 1,ngrid
    391357        zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
    392     enddo
    393 enddo
     358    end do
     359end do
    394360
    395361! c. Vapor pressure
    396 pco2 = mass_mean/m_co2*q_co2*zplev_mean
    397 pco2_avg(:) = sum(pco2(:,:),2)/timelen
     362pco2 = mass_mean/m_co2*q_co2_ts*zplev_mean
     363pco2_avg(:) = sum(pco2(:,:),2)/nday
    398364
    399365deallocate(zplev_mean,mass_mean,pco2)
    400366
    401367! 1. Compute the fraction of the pores occupied by H2O
    402 call 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)
     368call evolve_h2o_ads(d_h2oice,d_co2ice,h2oice,co2ice,ps,q_co2_ts,q_h2o_ts,tsoil,TI,theta_h2o_ads,m_h2o_adsorbed,delta_mh2o)
    404369
    405370! 2. we compute the mass of co2 adsorded in each layer of the meshes
    … …  
    407372    do islope = 1,nslope
    408373        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)))
     374            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
     375                dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
     376                                                         (alpha*pco2_avg(ig) + sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
    412377            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)))
     378                if (abs(co2_ads_reg(ig,iloop,islope)) < minieps) then !!! we are at first call
     379                    dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1._dp - theta_h2o_ads(ig,iloop,islope))*alpha*pco2_avg(ig) &
     380                                                             /(alpha*pco2_avg(ig)+sqrt(tsoil(ig,iloop,islope))*exp(beta/tsoil(ig,iloop,islope)))
    416381                else ! no change: permanent ice stick the atoms of CO2
    417                     dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
    418                 endif
    419             endif
    420         enddo
    421     enddo
    422 enddo
     382                    dm_co2_regolith_slope(ig,iloop,islope) = co2_ads_reg(ig,iloop,islope)
     383                end if
     384            end if
     385        end do
     386    end do
     387end do
    423388
    424389! 3. Check the exchange between the atmosphere and the regolith
    425390do ig = 1,ngrid
    426     delta_mreg(ig) = 0.
     391    delta_co2_ads(ig) = 0._dp
    427392    do islope = 1,nslope
    428         deltam_reg_slope(ig,islope) = 0.
     393        deltam_reg_slope(ig,islope) = 0._dp
    429394        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
     395            if (TI(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
    431396                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))
     397                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop))
    433398                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))
    435                 endif
     399                    deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - co2_ads_reg(ig,iloop,islope))*(layer(iloop) - layer(iloop - 1))
     400                end if
    436401            else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
    437                 deltam_reg_complete(ig,iloop,islope) = 0.
    438             endif
     402                deltam_reg_complete(ig,iloop,islope) = 0._dp
     403            end if
    439404            deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
    440         enddo
    441         delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    442     enddo
    443 enddo
    444 m_co2_completesoil = dm_co2_regolith_slope
    445 
    446 END SUBROUTINE regolith_co2adsorption
     405        end do
     406        delta_co2_ads(ig) = delta_co2_ads(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180._dp)
     407    end do
     408end do
     409co2_ads_reg = dm_co2_regolith_slope
     410
     411END SUBROUTINE evolve_co2_ads
     412!=======================================================================
     413
     414!=======================================================================
     415SUBROUTINE compute_totmass_adsorbed(h2o_ads_reg,co2_ads_reg,totmass_adsco2,totmass_adsh2o)
     416!-----------------------------------------------------------------------
     417! NAME
     418!     compute_totmass_adsorbed
     419!
     420! DESCRIPTION
     421!     Compute the total mass of CO2 and H2O adsorbed in the regolith.
     422!
     423! AUTHORS & DATE
     424!     L. Lange, 2023
     425!     JB Clement, 12/2025
     426!     C. Metz, 02/2026
     427!
     428! NOTES
     429!
     430!-----------------------------------------------------------------------
     431
     432! DEPENDENCIES
     433! ------------
     434use geometry, only: ngrid, nslope, nsoil, cell_area
     435use slopes,   only: subslope_dist, def_slope_mean
     436use soil,     only: layer
     437use maths,    only: pi
     438use display,  only: print_msg
     439use utility,  only: real2str
     440
     441! DECLARATION
     442! -----------
     443implicit none
     444
     445! ARGUMENTS
     446! ---------
     447real(dp), dimension(:,:,:), intent(in)  :: h2o_ads_reg, co2_ads_reg
     448real(qp),                   intent(out) :: totmass_adsco2, totmass_adsh2o
     449
     450! LOCAL VARIABLES
     451! ---------------
     452integer(di) :: i, islope, l
     453real(dp)    :: thickness, slope_corr
     454
     455! CODE
     456! ----
     457totmass_adsco2 = 0._qp
     458totmass_adsh2o = 0._qp
     459do i = 1,ngrid
     460    do islope = 1,nslope
     461        slope_corr = subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180._dp)
     462        do l = 1,nsoil
     463            if (l == 1) then
     464                thickness = layer(l)
     465            else
     466                thickness = layer(l) - layer(l-1)
     467            end if
     468            totmass_adsco2 = totmass_adsco2 + co2_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i)
     469            totmass_adsh2o = totmass_adsh2o + h2o_ads_reg(i,l,islope)*thickness*slope_corr*cell_area(i)
     470        end do
     471    end do
     472end do
     473call print_msg("Total mass of CO2 adsorbed in the regolith [kg] = "//real2str(totmass_adsco2))
     474call print_msg("Total mass of H2O adsorbed in the regolith [kg] = "//real2str(totmass_adsh2o))
     475
     476END SUBROUTINE compute_totmass_adsorbed
    447477!=======================================================================
    448478
Note: See TracChangeset for help on using the changeset viewer.