Index: /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 3455)
+++ /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 3456)
@@ -1,11 +1,13 @@
-module adsorption_mod
-
-  implicit none
-
-  LOGICAL adsorption_pem ! True by default, to compute adsorption/desorption. Read in  pem.def
-  real, save, allocatable :: co2_adsorbded_phys(:,:,:)  ! co2 that is in the regolith [kg/m^2]
-  real, save, allocatable :: h2o_adsorbded_phys(:,:,:)  ! h2o that is in the regolith [kg/m^2]
-
-  contains
+MODULE adsorption_mod
+
+implicit none
+
+logical                                   :: adsorption_pem     ! True by default, to compute adsorption/desorption. Read in pem.def
+real, save, allocatable, dimension(:,:,:) :: co2_adsorbded_phys ! co2 that is in the regolith [kg/m^2]
+real, save, allocatable, dimension(:,:,:) :: h2o_adsorbded_phys ! h2o that is in the regolith [kg/m^2]
+
+!=======================================================================
+contains
+!=======================================================================
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -16,67 +18,62 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  subroutine ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
-
-  implicit none
-  integer,intent(in) :: ngrid ! number of atmospheric columns
-  integer,intent(in) :: nslope ! number of slope within a mesh 
-  integer,intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM
-    allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
-    allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
-  end subroutine ini_adsorption_h_PEM
-
-!!! -----------------------------------------------
-
-  subroutine end_adsorption_h_PEM
-
-  implicit none
-    if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
-    if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys)
-  end subroutine end_adsorption_h_PEM
-
-!!! -----------------------------------------------
-
-  subroutine regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
-                                m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg)
-
-! inputs
- INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! size dimension: physics x subslope x soil x timeseries
- REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) !tendancies on the glaciers [1]
- REAL,INTENT(IN) :: waterice(ngrid,nslope)              ! water ice at the surface [kg/m^2]
- REAL,INTENT(IN) :: co2ice(ngrid,nslope)                ! co2 ice at the surface [kg/m^2]
- REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)      ! Soil Thermal inertia (J/K/^2/s^1/2)
- REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)   ! Soil temperature (K)
- REAL,INTENT(IN) :: ps(ngrid,timelen)                   ! Average surface pressure [Pa]
- REAL,INTENT(IN) :: q_co2(ngrid,timelen)                ! Mass mixing ratio of co2 in the first layer (kg/kg)
- REAL,INTENT(IN) :: q_h2o(ngrid,timelen)                ! Mass mixing ratio of H2o in the first layer (kg/kg)
-
-! outputs
- REAL,INTENT(INOUT) :: m_h2o_completesoil(ngrid,nsoil_PEM,nslope) ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)     
- REAL,INTENT(OUT) :: delta_mh2oreg(ngrid)                         ! Difference density of h2o adsorbed (kg/m^3)
-
- REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)   ! Density of co2 adsorbed (kg/m^3)
- REAL,INTENT(OUT) :: delta_mco2reg(ngrid)                            ! Difference density of co2 adsorbed (kg/m^3)
-  
-! local variables
- REAL :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope) ! Fraction of the pores occupied by H2O molecules
+SUBROUTINE ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
+
+implicit none
+
+integer, intent(in) :: ngrid       ! number of atmospheric columns
+integer, intent(in) :: nslope      ! number of slope within a mesh
+integer, intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM
+
+allocate(co2_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
+allocate(h2o_adsorbded_phys(ngrid,nsoilmx_PEM,nslope))
+
+END SUBROUTINE ini_adsorption_h_PEM
+
+!=======================================================================
+SUBROUTINE end_adsorption_h_PEM
+
+if (allocated(co2_adsorbded_phys)) deallocate(co2_adsorbded_phys)
+if (allocated(h2o_adsorbded_phys)) deallocate(h2o_adsorbded_phys)
+
+END SUBROUTINE end_adsorption_h_PEM
+
+!=======================================================================
+SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
+                               m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg)
+
+implicit none
+
+! Inputs
+integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM, timelen  ! size dimension: physics x subslope x soil x timeseries
+real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! tendancies on the glaciers [1]
+real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! water ice at the surface [kg/m^2]
+real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! co2 ice at the surface [kg/m^2]
+real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Soil Thermal inertia (J/K/^2/s^1/2)
+real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Soil temperature (K)
+real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Average surface pressure [Pa]
+real,    dimension(ngrid,timelen),          intent(in) :: q_co2                              ! Mass mixing ratio of co2 in the first layer (kg/kg)
+real,    dimension(ngrid,timelen),          intent(in) :: q_h2o                              ! Mass mixing ratio of H2o in the first layer (kg/kg)
+
+! Outputs
+real, dimension(ngrid), intent(out) :: delta_mh2oreg ! Difference density of h2o adsorbed (kg/m^3)
+real, dimension(ngrid), intent(out) :: delta_mco2reg ! Difference density of co2 adsorbed (kg/m^3)
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil ! Density of co2 adsorbed (kg/m^3)
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil ! Density of h2o adsorbed (kg/m^3)
+
+! Local variables
+real, dimension(ngrid,nsoil_PEM,nslope) :: theta_h2o_adsorbded ! Fraction of the pores occupied by H2O molecules
 ! -------------
-
-
 ! Compute H2O adsorption, then CO2 adsorption
-
-  call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
-                                   theta_h2o_adsorbded,m_h2o_completesoil,delta_mh2oreg)
-
-
-  call 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_mco2reg)
-
-   RETURN
-  end subroutine
-
-!------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-
-  subroutine regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
-                                   theta_h2o_adsorbded,m_h2o_completesoil,delta_mreg)
+call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
+                            theta_h2o_adsorbded,m_h2o_completesoil,delta_mh2oreg)
+call 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_mco2reg)
+
+END SUBROUTINE
+
+!=======================================================================
+SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
+                                  theta_h2o_adsorbded,m_h2o_completesoil,delta_mreg)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -90,5 +87,5 @@
 use comsoil_h_PEM,         only: layer_PEM, index_breccia
 use comslope_mod,          only: subslope_dist, def_slope_mean
-use vertical_layers_mod,   only: ap,bp
+use vertical_layers_mod,   only: ap, bp
 use constants_marspem_mod, only: alpha_clap_h2o, beta_clap_h2o, m_h2o, m_co2,m_noco2, rho_regolith
 
@@ -99,145 +96,125 @@
 #endif
 
- implicit none
-! inputs
- INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension
- REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) ! Tendencies on the glaciers ()
- REAL,INTENT(IN) :: waterice(ngrid,nslope)              ! Water ice at the surface [kg/m^2]
- REAL,INTENT(IN) :: co2ice(ngrid,nslope)                ! CO2 ice at the surface [kg/m^2]
- REAL,INTENT(IN) :: ps(ngrid,timelen)                   ! Surface pressure (Pa)          
- REAL,INTENT(IN) :: q_co2(ngrid,timelen)                ! Mass mixing ratio of co2 in the first layer (kg/kg)
- REAL,INTENT(IN) :: q_h2o(ngrid,timelen)                ! Mass mixing ratio of H2o in the first layer (kg/kg)
- REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)      ! Soil Thermal inertia (J/K/^2/s^1/2)
- REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)   ! Soil temperature (K)
-
-! outputs
- REAL,INTENT(INOUT) :: m_h2o_completesoil(ngrid,nsoil_PEM,nslope) ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)     
- REAL,INTENT(OUT) :: theta_h2o_adsorbded(ngrid,nsoil_PEM,nslope)  ! Fraction of the pores occupied by H2O molecules
- REAL,INTENT(OUT) :: delta_mreg(ngrid)                            ! Difference density of h2o adsorbed (kg/m^3)
-
-! constants
- REAL :: Ko =  1.57e-8            ! Jackosky et al. 1997
- REAL :: e = 2573.9               ! Jackosky et al. 1997
- REAL :: mu = 0.48                ! Jackosky et al. 1997
- real :: m_theta = 2.84e-7        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
-! real :: as = 18.9e3             ! Specific area, Buhler & Piqueux 2021
- real :: as = 9.48e4              ! Specific area, Zent 
- real ::  inertie_thresold = 800. ! TI > 800 means cementation
-
- ! local variables
- REAL :: deltam_reg_complete(ngrid,index_breccia,nslope) ! Difference in the mass per slope and soil layer (kg/m^3)
- real :: K                        ! Used to compute theta
- integer ig, iloop, islope, it    ! For loops
- INTEGER :: ispermanent_co2glaciers(ngrid,nslope)      ! Check if the co2 glacier is permanent
- INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)      ! Check if the h2o glacier is permanent
- REAL :: deltam_reg_slope(ngrid,nslope)                ! Difference density of h2o adsorbed per slope (kg/m^3)
- REAL :: dm_h2o_regolith_slope(ngrid,nsoil_PEM,nslope) ! Elementary h2o mass adsorded per mesh per slope
- real :: A,B                                           ! Used to compute the mean mass above the surface
- real :: p_sat                                         ! Saturated vapor pressure of ice
- real,allocatable :: mass_mean(:,:)                    ! Mean mass above the surface
- real,allocatable :: zplev_mean(:,:)                   ! Pressure above the surface
- real,allocatable :: pvapor(:,:)                       ! Partial pressure above the surface
- real, allocatable :: pvapor_avg(:)                    ! Yearly averaged
+implicit none
+
+! Inputs
+integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM,timelen   ! Size dimension
+real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! Tendencies on the glaciers ()
+real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
+real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! CO2 ice at the surface [kg/m^2]
+real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Surface pressure (Pa)
+real,    dimension(ngrid,timelen),          intent(in) :: q_co2                              ! Mass mixing ratio of co2 in the first layer (kg/kg)
+real,    dimension(ngrid,timelen),          intent(in) :: q_h2o                              ! Mass mixing ratio of H2o in the first layer (kg/kg)
+real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Soil Thermal inertia (J/K/^2/s^1/2)
+real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Soil temperature (K)
+
+! Outputs
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)
+real, dimension(ngrid,nsoil_PEM,nslope), intent(out) :: theta_h2o_adsorbded ! Fraction of the pores occupied by H2O molecules
+real, dimension(ngrid),                  intent(out) :: delta_mreg          ! Difference density of h2o adsorbed (kg/m^3)
+
+! Constants
+real :: Ko =  1.57e-8            ! Jackosky et al. 1997
+real :: e = 2573.9               ! Jackosky et al. 1997
+real :: mu = 0.48                ! Jackosky et al. 1997
+real :: m_theta = 2.84e-7        ! Mass of h2o per m^2 absorbed Jackosky et al. 1997
+! real :: as = 18.9e3              ! Specific area, Buhler & Piqueux 2021
+real :: as = 9.48e4              ! Specific area, Zent
+real ::  inertie_thresold = 800. ! TI > 800 means cementation
+
+! Local variables
+real,    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer (kg/m^3)
+real                                           :: K                       ! Used to compute theta
+integer                                        :: ig, iloop, islope, it   ! For loops
+logical, dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
+logical, dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
+real,    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference density of h2o adsorbed per slope (kg/m^3)
+real,    dimension(ngrid,nsoil_PEM,nslope)     :: dm_h2o_regolith_slope   ! Elementary h2o mass adsorded per mesh per slope
+real                                           :: A, B                    ! Used to compute the mean mass above the surface
+real                                           :: p_sat                   ! Saturated vapor pressure of ice
+real,    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
+real,    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
+real,    dimension(:,:), allocatable           :: pvapor                  ! Partial pressure above the surface
+real,    dimension(:)  , allocatable           :: pvapor_avg              ! Yearly averaged
 
 ! 0. Some initializations
-
-
-     allocate(mass_mean(ngrid,timelen))
-     allocate(zplev_mean(ngrid,timelen))
-     allocate(pvapor(ngrid,timelen))
-     allocate(pvapor_avg(ngrid))
-     A =(1/m_co2 - 1/m_noco2)
-     B=1/m_noco2
-     theta_h2o_adsorbded(:,:,:) = 0.
-     dm_h2o_regolith_slope(:,:,:) = 0.
-
-#ifndef CPP_STD
-
-!0.1 Look at perennial ice
-  do ig = 1,ngrid
-    do islope = 1,nslope
-     if((abs(tend_h2oglaciers(ig,islope)).gt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then
-        ispermanent_h2oglaciers(ig,islope) = 1
-     else
-        ispermanent_h2oglaciers(ig,islope) = 0
-     endif
-
-     if((abs(tend_co2glaciers(ig,islope)).gt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
-        ispermanent_co2glaciers(ig,islope) = 1
-     else
-        ispermanent_co2glaciers(ig,islope) = 0
-     endif
-    enddo
-   enddo
-
-!   0.2 Compute the partial pressure of vapor
-!a. the molecular mass into the column
-     do ig = 1,ngrid
-       mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B)
-     enddo
-
-
-! b. pressure level
-     do it = 1,timelen
-       do ig = 1,ngrid
-         zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
-       enddo
-     enddo
-! c. Vapor pressure
-     pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:)
-     pvapor_avg(:) = sum(pvapor(:,:),2)/timelen
-#endif
-     deallocate(pvapor)
-     deallocate(zplev_mean)
-     deallocate(mass_mean)
-#ifndef CPP_STD
-
-! 1. we compute the mass of H2O adsorded in each layer of the meshes  
-
- do ig = 1,ngrid
-  do islope = 1,nslope
-    do iloop = 1,index_breccia
-      K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
-      if(TI_PEM(ig,iloop,islope).lt.inertie_thresold)  then
-        theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu
-      else
-        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 ...
-        theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1+K*p_sat))**mu
-      endif
-        dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*m_theta*rho_regolith
-   enddo
-  enddo
- enddo
-
-! 2. Check the exchange between the atmosphere and the regolith
-
-  do ig = 1,ngrid
-   delta_mreg(ig) = 0.
-   do islope = 1,nslope
-    deltam_reg_slope(ig,islope) = 0.
-    do iloop = 1,index_breccia
-       if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
-             if(iloop==1) then
-                 deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope)) &
-                                                   *(layer_PEM(iloop))
-             else
-                 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))
-             endif
-       else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
-             deltam_reg_complete(ig,iloop,islope) = 0.
-       endif
-       deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
-    enddo
-   delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-   enddo
-  enddo
-   m_h2o_completesoil(:,:,:) = dm_h2o_regolith_slope(:,:,:)
-
- RETURN
-#endif
-  end subroutine
-
-!------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pvapor(ngrid,timelen),pvapor_avg(ngrid))
+A = 1./m_co2 - 1./m_noco2
+B = 1./m_noco2
+theta_h2o_adsorbded = 0.
+dm_h2o_regolith_slope = 0.
+ispermanent_h2oglaciers = .false.
+ispermanent_co2glaciers = .false.
+
+#ifndef CPP_STD
+    ! 0.1 Look at perennial ice
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            if (abs(tend_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
+            if (abs(tend_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
+        enddo
+    enddo
+
+    ! 0.2 Compute the partial pressure of vapor
+    ! a. the molecular mass into the column
+    do ig = 1,ngrid
+        mass_mean(ig,:) = 1/(A*q_co2(ig,:) + B)
+    enddo
+
+    ! b. pressure level
+    do it = 1,timelen
+        do ig = 1,ngrid
+            zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
+        enddo
+    enddo
+
+    ! c. Vapor pressure
+    pvapor = mass_mean/m_h2o*q_h2o*zplev_mean
+    pvapor_avg = sum(pvapor,2)/timelen
+#endif
+deallocate(pvapor,zplev_mean,mass_mean)
+
+#ifndef CPP_STD
+    ! 1. we compute the mass of H2O adsorded in each layer of the meshes
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            do iloop = 1,index_breccia
+                K = Ko*exp(e/tsoil_PEM(ig,iloop,islope))
+                if (TI_PEM(ig,iloop,islope) < inertie_thresold) then
+                    theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1. + K*pvapor_avg(ig)))**mu
+                else
+                    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 ...
+                    theta_h2o_adsorbded(ig,iloop,islope) = (K*p_sat/(1. + K*p_sat))**mu
+                endif
+                dm_h2o_regolith_slope(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*m_theta*rho_regolith
+            enddo
+        enddo
+    enddo
+
+    ! 2. Check the exchange between the atmosphere and the regolith
+    do ig = 1,ngrid
+        delta_mreg(ig) = 0.
+        do islope = 1,nslope
+            deltam_reg_slope(ig,islope) = 0.
+            do iloop = 1,index_breccia
+                if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
+                    if (iloop == 1) then
+                        deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope))*(layer_PEM(iloop))
+                    else
+                        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))
+                    endif
+                else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
+                    deltam_reg_complete(ig,iloop,islope) = 0.
+                endif
+                deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
+            enddo
+            delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+        enddo
+    enddo
+    m_h2o_completesoil = dm_h2o_regolith_slope
+#endif
+END SUBROUTINE
+
+!=======================================================================
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -247,7 +224,5 @@
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-  SUBROUTINE 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)
+SUBROUTINE 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)
 
 use comsoil_h_PEM,         only: layer_PEM, index_breccia, index_breccia
@@ -262,156 +237,131 @@
 #endif
 
-      IMPLICIT NONE
-! Inputs:  
- INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM,timelen             ! Size dimension
- REAL,INTENT(IN) :: ps(ngrid,timelen)                               ! Average surface pressure [Pa]
- REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)               ! Mean Soil Temperature [K]
- REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)                  ! Mean Thermal Inertia [USI]
- REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) ! Tendencies on the glaciers ()
- REAL,INTENT(IN) :: q_co2(ngrid,timelen),q_h2o(ngrid,timelen)       ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
- REAL,INTENT(IN) :: waterice(ngrid,nslope)                          ! Water ice at the surface [kg/m^2]
- REAL,INTENT(IN) :: co2ice(ngrid,nslope)                            ! CO2 ice at the surface [kg/m^2] 
+implicit none
+
+! Inputs:
+integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM,timelen   ! Size dimension
+real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Average surface pressure [Pa]
+real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Mean Soil Temperature [K]
+real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Mean Thermal Inertia [USI]
+real,    dimension(ngrid,nslope),           intent(in) :: tend_h2oglaciers, tend_co2glaciers ! Tendencies on the glaciers ()
+real,    dimension(ngrid,timelen),          intent(in) :: q_co2, q_h2o                       ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
+real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
+real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! CO2 ice at the surface [kg/m^2]
 
 ! Outputs:
- REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)   ! Density of co2 adsorbed (kg/m^3)
- REAL,INTENT(OUT) :: delta_mreg(ngrid)                              ! Difference density of co2 adsorbed (kg/m^3)
-  
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil ! Density of co2 adsorbed (kg/m^3)
+real, dimension(ngrid), intent(out) :: delta_mreg ! Difference density of co2 adsorbed (kg/m^3)
+
 ! Constants:
-
- REAL :: alpha = 7.512e-6         ! Zent & Quinn 1995
- REAL :: beta =  -1541.5          ! Zent & Quinn 1995
- REAL ::  inertie_thresold = 800. ! TI > 800 means cementation
- real :: m_theta = 4.27e-7        ! Mass of co2 per m^2 absorbed 
+real :: alpha = 7.512e-6        ! Zent & Quinn 1995
+real :: beta = -1541.5          ! Zent & Quinn 1995
+real :: inertie_thresold = 800. ! TI > 800 means cementation
+real :: m_theta = 4.27e-7       ! Mass of co2 per m^2 absorbed
 ! real :: as = 18.9e3             ! Specific area, Buhler & Piqueux 2021
- real :: as = 9.48e4              ! Same as previous but from zent 
-! Local          
- real :: A,B                                             ! Used to compute the mean mass above the surface
- INTEGER :: ig,islope,iloop,it                           ! For loops
- REAL :: dm_co2_regolith_slope(ngrid,nsoil_PEM,nslope)   ! Elementary mass adsorded per mesh per slope
- INTEGER :: ispermanent_co2glaciers(ngrid,nslope)        ! Check if the co2 glacier is permanent
- INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)        ! Check if the h2o glacier is permanent
- REAL :: deltam_reg_complete(ngrid,index_breccia,nslope) ! Difference in the mass per slope and soil layer (kg/m^3)
- REAL :: deltam_reg_slope(ngrid,nslope)                  ! Difference in the mass per slope  (kg/m^3)
- REAL :: m_h2o_adsorbed(ngrid,nsoil_PEM,nslope)          ! Density of CO2 adsorbed (kg/m^3)
- REAL :: theta_h2o_adsorbed(ngrid,nsoil_PEM,nslope)      ! Fraction of the pores occupied by H2O molecules
- REAL :: delta_mh2o(ngrid)                               ! Difference density of h2o adsorbed (kg/m^3)
-!timelen array are allocated because heavy ...
- real,allocatable :: mass_mean(:,:)                      ! Mean mass above the surface
- real,allocatable :: zplev_mean(:,:)                     ! Pressure above the surface
- real,allocatable :: pco2(:,:)                           ! Partial pressure above the surface
- real, allocatable :: pco2_avg(:)                        ! Yearly averaged
+real :: as = 9.48e4             ! Same as previous but from zent
+
+! Local
+real                                           :: A, B                    ! Used to compute the mean mass above the surface
+integer                                        :: ig, islope, iloop, it   ! For loops
+real,    dimension(ngrid,nsoil_PEM,nslope)     :: dm_co2_regolith_slope   ! Elementary mass adsorded per mesh per slope
+logical, dimension(ngrid,nslope)               :: ispermanent_co2glaciers ! Check if the co2 glacier is permanent
+logical, dimension(ngrid,nslope)               :: ispermanent_h2oglaciers ! Check if the h2o glacier is permanent
+real,    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer (kg/m^3)
+real,    dimension(ngrid,nslope)               :: deltam_reg_slope        ! Difference in the mass per slope  (kg/m^3)
+real,    dimension(ngrid,nsoil_PEM,nslope)     :: m_h2o_adsorbed          ! Density of CO2 adsorbed (kg/m^3)
+real,    dimension(ngrid,nsoil_PEM,nslope)     :: theta_h2o_adsorbed      ! Fraction of the pores occupied by H2O molecules
+real,    dimension(ngrid)                      :: delta_mh2o              ! Difference density of h2o adsorbed (kg/m^3)
+! timelen array are allocated because heavy...
+real,    dimension(:,:), allocatable           :: mass_mean               ! Mean mass above the surface
+real,    dimension(:,:), allocatable           :: zplev_mean              ! Pressure above the surface
+real,    dimension(:,:), allocatable           :: pco2                    ! Partial pressure above the surface
+real,    dimension(:),   allocatable           :: pco2_avg                ! Yearly averaged
 
 ! 0. Some initializations
-
-     allocate(mass_mean(ngrid,timelen))
-     allocate(zplev_mean(ngrid,timelen))
-     allocate(pco2(ngrid,timelen))
-     allocate(pco2_avg(ngrid))
-
-
-  
-      m_h2o_adsorbed(:,:,:) = 0.
-      A =(1/m_co2 - 1/m_noco2)
-      B=1/m_noco2
-
-     dm_co2_regolith_slope(:,:,:) = 0
-     delta_mreg(:) = 0.
-
-#ifndef CPP_STD
-
-!0.1 Look at perennial ice
-  do ig = 1,ngrid
-    do islope = 1,nslope
-     if((abs(tend_h2oglaciers(ig,islope)).gt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then
-        ispermanent_h2oglaciers(ig,islope) = 1
-     else
-        ispermanent_h2oglaciers(ig,islope) = 0
-     endif
-
-     if((abs(tend_co2glaciers(ig,islope)).gt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
-        ispermanent_co2glaciers(ig,islope) = 1
-     else
-        ispermanent_co2glaciers(ig,islope) = 0
-     endif
-    enddo
-   enddo
-
-!   0.2  Compute the partial pressure of CO2
-!a. the molecular mass into the column
-     do ig = 1,ngrid
-       mass_mean(ig,:) = 1./(A*q_co2(ig,:) +B)
-     enddo
-
-! b. pressure level
-     do it = 1,timelen
-       do ig = 1,ngrid
-         zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
-       enddo
-     enddo
-
-! c. Vapor pressure
-     pco2(:,:) = mass_mean(:,:)/m_co2*q_co2(:,:)*zplev_mean(:,:)
-     pco2_avg(:) = sum(pco2(:,:),2)/timelen
-
-     deallocate(zplev_mean)
-     deallocate(mass_mean)
-     deallocate(pco2)
-
-
-! 1. Compute the fraction of the pores occupied by H2O
-
- call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
-                                   theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
-
-! 2.  we compute the mass of co2 adsorded in each layer of the meshes  
-
- do ig = 1,ngrid
-  do islope = 1,nslope
-    do iloop = 1,index_breccia
-     if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
-     dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
-                                             (alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
-     else
-        if(abs(m_co2_completesoil(ig,iloop,islope)).lt.(1e-10)) then !!! we are at first call
-          dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1-theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) &
-                                                 /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope))) 
-        else ! no change: permanent ice stick the atoms of CO2
-          dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
-        endif
-     endif
-    enddo
- enddo
-enddo
-
-! 3. Check the exchange between the atmosphere and the regolith
-
-  do ig = 1,ngrid
-   delta_mreg(ig) = 0.
-   do islope = 1,nslope
-    deltam_reg_slope(ig,islope) = 0.
-    do iloop = 1,index_breccia
-       if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
-             if(iloop == 1) then 
-                 deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) &
-                                                   *(layer_PEM(iloop))               
-             else
-                 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))
-             endif
-       else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
-             deltam_reg_complete(ig,iloop,islope) = 0.
-       endif
-       deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
-    enddo
-   delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
-   enddo
-  enddo
-  m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:)
-
-!=======================================================================
-      RETURN
-#endif
-      END
-
-
-   end module
+allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pco2(ngrid,timelen),pco2_avg(ngrid))
+m_h2o_adsorbed = 0.
+A = 1./m_co2 - 1./m_noco2
+B = 1./m_noco2
+dm_co2_regolith_slope = 0.
+delta_mreg = 0.
+ispermanent_h2oglaciers = .false.
+ispermanent_co2glaciers = .false.
+
+#ifndef CPP_STD
+    ! 0.1 Look at perennial ice
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            if (abs(tend_h2oglaciers(ig,islope)) > 1.e-5 .and. abs(waterice(ig,islope)) > 0.) ispermanent_h2oglaciers(ig,islope) = .true.
+            if (abs(tend_co2glaciers(ig,islope)) > 1.e-5 .and. abs(co2ice(ig,islope)) > 0.) ispermanent_co2glaciers(ig,islope) = .true.
+        enddo
+    enddo
+
+    ! 0.2  Compute the partial pressure of CO2
+    ! a. the molecular mass into the column
+    do ig = 1,ngrid
+        mass_mean(ig,:) = 1./(A*q_co2(ig,:) + B)
+    enddo
+
+    ! b. pressure level
+    do it = 1,timelen
+        do ig = 1,ngrid
+            zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it)
+        enddo
+    enddo
+
+    ! c. Vapor pressure
+    pco2 = mass_mean/m_co2*q_co2*zplev_mean
+    pco2_avg(:) = sum(pco2(:,:),2)/timelen
+
+    deallocate(zplev_mean,mass_mean,pco2)
+
+    ! 1. Compute the fraction of the pores occupied by H2O
+    call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
+                                theta_h2o_adsorbed, m_h2o_adsorbed,delta_mh2o)
+
+    ! 2. we compute the mass of co2 adsorded in each layer of the meshes
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            do iloop = 1,index_breccia
+                if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
+                    dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig)/ &
+                                                             (alpha*pco2_avg(ig) + sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
+                else
+                    if (abs(m_co2_completesoil(ig,iloop,islope)) < 1.e-10) then !!! we are at first call
+                        dm_co2_regolith_slope(ig,iloop,islope) = as*rho_regolith*m_theta*(1. - theta_h2o_adsorbed(ig,iloop,islope))*alpha*pco2_avg(ig) &
+                                                                 /(alpha*pco2_avg(ig)+sqrt(tsoil_PEM(ig,iloop,islope))*exp(beta/tsoil_PEM(ig,iloop,islope)))
+                    else ! no change: permanent ice stick the atoms of CO2
+                        dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
+                    endif
+                endif
+            enddo
+        enddo
+    enddo
+
+    ! 3. Check the exchange between the atmosphere and the regolith
+    do ig = 1,ngrid
+        delta_mreg(ig) = 0.
+        do islope = 1,nslope
+            deltam_reg_slope(ig,islope) = 0.
+            do iloop = 1,index_breccia
+                if (TI_PEM(ig,iloop,islope) < inertie_thresold .and. .not. ispermanent_h2oglaciers(ig,islope) .and. .not. ispermanent_co2glaciers(ig,islope)) then
+                    if (iloop == 1) then
+                        deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope))*(layer_PEM(iloop))
+                    else
+                        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))
+                    endif
+                else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
+                    deltam_reg_complete(ig,iloop,islope) = 0.
+                endif
+                deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
+            enddo
+            delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+        enddo
+    enddo
+    m_co2_completesoil = dm_co2_regolith_slope
+#endif
+
+END SUBROUTINE regolith_co2adsorption
+
+END MODULE adsorption_mod
+
Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3455)
+++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3456)
@@ -436,2 +436,5 @@
 == 30/09/2024 == JBC
 Correction of the launching script for a relaunch situation.
+
+== 11/10/2024 == JBC
+Cleaning of the adsorption module to make the debugging easier.
Index: /trunk/LMDZ.COMMON/libf/evolution/deftank/multiple_exec.sh
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/deftank/multiple_exec.sh	(revision 3455)
+++ /trunk/LMDZ.COMMON/libf/evolution/deftank/multiple_exec.sh	(revision 3456)
@@ -1,6 +1,6 @@
 #!/bin/bash
-#############################################################
-### Script to exectute mutliple scripts in subdirectories ###
-#############################################################
+############################################################
+### Script to execute multiple scripts in subdirectories ###
+############################################################
 
 # Name of the file to execute in all subdirectories
Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3455)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3456)
@@ -217,5 +217,5 @@
 ! b. Special case for inertiedat, inertiedat_PEM
         call get_field("inertiedat_PEM",inertiedat_PEM,found)
-        if (.not.found) then
+        if (.not. found) then
             do iloop = 1,nsoil_PCM
                 inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
@@ -259,5 +259,5 @@
             write(num,fmt='(i2.2)') islope
             call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
-            if (.not.found) then
+            if (.not. found) then
                 write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
                 write(*,*)'will reconstruct the values of Tsoil'
@@ -346,5 +346,5 @@
 
             if (.not. found) deltam_co2_regolith_phys = 0.
-            if (.not.found2) deltam_h2o_regolith_phys = 0.
+            if (.not. found2) deltam_h2o_regolith_phys = 0.
 
             write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
