Index: trunk/LMDZ.MARS/libf/phymars/soilwater.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/soilwater.F90	(revision 3219)
+++ trunk/LMDZ.MARS/libf/phymars/soilwater.F90	(revision 3223)
@@ -1,5 +1,5 @@
-subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf,tsoil, ptimestep,  & 
+subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep,  & 
       exchange, qsat_surf, pq, pa, pb, pc, pd, pdqsdifpot, pqsurf,  & 
-      pqsoil, pplev, rhoatmo,writeoutput,zdqsdifrego, zq1temp2) 
+      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2) 
 
 
@@ -10,5 +10,4 @@
       use geometry_mod, only: cell_area, latitude_deg
       use write_output_mod, only: write_output
-
 implicit none 
 
@@ -20,5 +19,5 @@
 ! suburface and the atmosphere.
 ! 
-! Water vapor and adsorbed water are treated as two separate subsurface tracers that are connected
+! Water vapor and nsoiladsorbed water are treated as two separate subsurface tracers that are connected
 ! by adsorption / desorption coefficients, including an adsorption / desorption kinetic factor.
 ! 
@@ -26,4 +25,5 @@
 ! of adsorbed water molecules (i.e. an approximation of a Langmuir - type isotherm). The slope of the
 ! isotherm is defined by the enthalpy of adsorption (in kJ / mol).
+! New since 2020 : the adsorption isotherm is now linear by part, to better simulate Type II or Type III isotherms (and not only Type I)
 ! 
 ! The linearized adsorption - diffusion equation is solved first, then the water vapor pressure is
@@ -68,10 +68,10 @@
 integer, intent(in) :: nqsoil
 logical, save :: firstcall_soil = .true.             ! triggers the initialization 
-real, intent(in) :: tsoil(ngrid, nsoil)         ! Subsurface temperatures
+real, intent(in) :: ptsoil(ngrid, nsoil)         ! Subsurface temperatures
 real, intent(in) :: ptsrf(ngrid)                ! Surface temperature
 real, intent(in) :: ptimestep                   ! length of the timestep (unit depends on run.def file)
 logical, intent(in) :: exchange(ngrid)          ! logical :: array describing whether there is exchange with the atmosphere at the current timestep
 
-real, intent(in) :: qsat_surf(ngrid)            ! Saturation mixing ratio at the surface
+real, intent(in) :: qsat_surf(ngrid)           ! Saturation vapor pressure at the current surface temperature (formerly qsat)
 real, intent(in) :: pq(ngrid, nlayer, nq)       ! Tracer atmospheric mixing ratio
 real, intent(in) :: pa(ngrid, nlayer)           ! Coefficients za
@@ -79,16 +79,15 @@
 real, intent(in) :: pc(ngrid, nlayer)           ! Coefficients zc
 real, intent(in) :: pd(ngrid, nlayer)           ! Coefficients zd
-real, intent(in) :: pdqsdifpot(ngrid)       ! Potential pdqsdif (without subsurface - atmosphere exchange)
+real, intent(in) :: pdqsdifpot(ngrid, nq)       ! Potential pdqsdif (without subsurface - atmosphere exchange)
 
 real, intent(in) :: pplev(ngrid, nlayer+1)      ! Atmospheric pressure levels
 real, intent(in) :: rhoatmo(ngrid)              ! Atmospheric air density (1st layer) (not used right now)
-
-logical, intent(in) :: writeoutput              ! just to check we are at the last subtimestep and we can write
+logical, intent(in) :: writeoutput              ! just to check we are at the last subtimestep and we
 
 ! Variables modified :
 ! ----------------------------------------------------------------------
-real, intent(inout) :: pqsoil(ngrid, nsoil, nqsoil) ! Subsurface tracers (water vapor and ice) 
-real, intent(in) :: pqsurf(ngrid)           ! Water ice on the surface
-                                                      ! (in kg.m - 3 : m - 3 of pore air for water vapor and m - 3 of regolith for water ice and adsorbed water)
+real, intent(inout) :: pqsoil(ngrid, nsoil, 3*nq) ! Subsurface tracers (water vapor and ice) 
+real, intent(in) :: pqsurf(ngrid)                 ! Water ice on the surface
+                                                  ! (in kg.m - 3 : m - 3 of pore air for water vapor and m - 3 of regolith for water ice and adsorbed water)
 ! Outputs :
 ! ---------------------------------------------------------------------- 
@@ -120,8 +119,17 @@
 real*8, allocatable, save :: adswater(:, :)     ! Adsorbed water in subsurface cells (at miD-layers) (...)
 real*8 :: adswater_temp(ngrid, nsoil)           ! temprory variable for adsorbed water
-logical :: ads_water_saturation_flag_2(nsoil) 
+logical, allocatable, save  :: over_mono_sat_flag(:, :) ! Formerly ads_water_saturation_flag_2(nsoil) (see descritption of the variable recompute_cell_ads_flag for an explanation ! ARNAU
 real*8 :: adswprev(ngrid, nsoil) 
-logical :: ads_water_saturation_flag_1(nsoil) 
-real*8, save :: adswater_sat                    ! Adsorption saturation value (kg.m - 3 of regolith) 
+logical :: recompute_cell_ads_flag(nsoil) ! ARNAU
+! Formerly ads_water_saturation_flag_1 but with a twist. This variable now
+! checks whether coefficients have changed and need recomputing. If the cell
+! is seen to be over the monolayer saturation level (i.e. the cell fulfills the
+! condition adswater_temp(ig, ik) > adswater_sat_mono) but the coefficients
+! have been computed assuming that the cell is below the monolayer saturation
+! layer (i.e. the variable over_mono_sat_flag(ig, ik) had been set to .false.),
+! then the cell needs to have its coefficients recomputed according to the
+! previous condition: adswater_temp(ig, ik) > adswater_sat_mono. Then,
+! the variable recompute_cell_ads_flag(ik) becomes .true.. ! ARNAU
+real*8, save :: adswater_sat_mono                    ! Adsorption monolayer saturation value (kg.m - 3 of regolith) ! ARNAU
 real*8 :: delta0(ngrid)                         ! Coefficient delta(0) modified
 real*8 :: alpha0(ngrid) 
@@ -129,14 +137,18 @@
 
 ! Adsorbtion/Desorption variables and parameters
-real*8 :: Ka(ngrid, nsoil)            ! Adsorption time constant (s - 1)
-real*8 :: Kd(ngrid, nsoil)            ! Desorption time constant (s - 1)
-real*8 :: k_ads_eq(ngrid, nsoil)      ! Equilibrium adsorption coefficient (unitless) (formerly kads)
-real*8, parameter :: kinetic_factor = 1      ! (inverse of) Characteristic time to reach adsorption equilibrium (s - 1): 
+real*8 :: Ka(ngrid, nsoil)            ! Adsorption time constant (s - 1) before monolayer saturation (first segment of the bilinear function)
+real*8 :: Kd(ngrid, nsoil)            ! Desorption time constant (s - 1) before monolayer saturation (first segment of the bilinear function)
+real*8 :: k_ads_eq(ngrid, nsoil)      ! Equilibrium adsorption coefficient (formerly kads) before monolayer saturation (first segment of the bilinear function); unitless (converts kg/m3 into kg/m3)
+real*8 :: Ka2(ngrid, nsoil)           ! Adsorption time constant (s - 1) after monolayer saturation (second segment of the bilinear function) ! modified 2020 
+real*8 :: Kd2(ngrid, nsoil)           ! Desorption time constant (s - 1) after monolayer saturation (second segment of the bilinear function) ! modified 2020
+real*8 :: k_ads_eq2(ngrid, nsoil)     ! Equilibrium adsorption coefficient (formerly kads) after monolayer saturation (second segment of the bilinear function); unitless ! modified 2020
+real*8 :: c0(ngrid, nsoil)            ! Intercept of the second line in the bilinear function ! modified 2020
+real*8, parameter :: kinetic_factor = 0.01      ! (inverse of) Characteristic time to reach adsorption equilibrium (s - 1): 
                                                 ! fixed arbitrarily when kinetics factors are unknown
                                                 ! possible values: ! = 1 / 1800 s - 1 ! / 1.16D-6 /  !( =  10 earth days)! / 1.D0 /  ! / 1.2D-5 /  !
                                           
-real*8, allocatable, save :: ztot1(:, :)  ! Total (water vapor +  ice) content (kg.m - 3 of soil) !? why does this not include adsorbed water?
+real*8, allocatable, save :: ztot1(:, :)  ! Total (water vapor +  ice) content (kg.m - 3 of soil) 
 real*8 :: dztot1(nsoil) 
-real*8 :: h2otot(ngrid, nsoil)      ! Total water content (ice +  water vapor +  adsorbed water) (....)
+real*8 :: h2otot(ngrid, nsoil)      ! Total water content (ice +  water vapor +  adsorbed water) 
 real*8 :: flux(ngrid, 0:nsoil - 1)  ! Fluxes at interfaces (kg.m - 2.s - 1) (positive = upward)
 real*8 :: rho(ngrid)                ! Air density (first subsurface layer)
@@ -171,5 +183,5 @@
 logical, allocatable, save :: close_top(:, :)         ! closing diffusion at the top of a layer if ice rises over saturation
 logical, allocatable, save :: close_bottom(:, :)      ! closing diffusion at the bottom of a layer if ice risies over saturation
-logical, parameter :: print_closure_warnings = .false.
+logical, parameter :: print_closure_warnings = .true.
 
 ! Coefficients for the Diffusion calculations
@@ -177,6 +189,8 @@
 real*8 :: B(ngrid, nsoil)           ! Coefficient in the diffusion formula
 real*8 :: C(ngrid, nsoil)           ! Coefficient in the diffusion formula
-real*8 :: E(ngrid, nsoil)           ! Coefficient in the diffusion formula
-real*8 :: F(ngrid, nsoil)           ! Coefficient in the diffusion formula
+real*8 :: E(ngrid, nsoil)           ! Coefficient in the diffusion formula (before monolayer saturation) ! added 2020 
+real*8 :: F(ngrid, nsoil)           ! Coefficient in the diffusion formula (before monolayer saturation) ! added 2020
+real*8 :: E2(ngrid, nsoil)           ! Coefficient in the diffusion formula (after monolayer saturation) ! added 2020
+real*8 :: F2(ngrid, nsoil)           ! Coefficient in the diffusion formula (after monolayer saturation) ! added 2020
 real*8, allocatable, save :: zalpha(:, :) ! Alpha coefficients
 real*8, allocatable, save :: zbeta(:, :)  ! Beta coefficients
@@ -187,5 +201,5 @@
 real*8, allocatable, save :: midlayer_dz(:, :)        ! Distance between the midcell points in m (formerly middz)
 
-real*8 :: nsat(ngrid, nsoil)                    ! amount of water vapor at adsoption saturation
+real*8 :: nsat(ngrid, nsoil)                    ! amount of water vapor at (adsorption) monolayer saturation ! modified 2020
 
 real*8, allocatable, save :: meshsize(:, :)     ! scaling/dimension factor for the por size
@@ -196,7 +210,11 @@
 real*8, parameter :: kb = 1.38065D-23     ! Boltzmann constant
 real*8, parameter :: mw = 2.988D-26		! Water molecular weight
+!real*8, parameter :: rho_H2O_ice = 920.D0    ! Ice density
 
 ! adsorption coefficients
-real*8, parameter :: enthalpy_ads = 21.D3 ! Enthalpy of adsorption (J.mole - 1) options 21.D3, 35.D3, and 60.D3
+real*8, parameter :: enthalpy_ads = 35.D3 ! Enthalpy of adsorption (J.mole - 1) options 21.D3, 35.D3, and 60.D3
+real*8, parameter :: enthalpy_ads2 = 21.D3 ! Enthalpy of adsorption (J.mole - 1) options 21.D3, 35.D3, and 60.D3 for the second segment ! added 2020
+real*8, parameter :: DeltaQ_ads = 21.D3 ! 10.D3 ! This is the DeltaQ_ads = heat of adsorption - enthalpy of liquefaction/sublimation = E1 - EL and may be equal to RT*ln(c), where c is the BET constant (from BET1938). This is for the first segment (J.mole - 1) ! added 2020
+real*8, parameter :: DeltaQ_ads2 = 21.D3 ! 10.D3 ! This is the DeltaQ_ads = heat of adsorption - enthalpy of liquefaction/sublimation = E1 - EL and may be equal to RT*ln(c), where c is the BET constant (from BET1938). This is for the second segment (J.mole - 1) ! added 2020 
 real*8, parameter :: tau0 = 1D-14 
 real*8, parameter :: S = 17.D3            ! Soil specific surface area (m2.kg - 1 of solid) options: 17.D3 and 106.D3
@@ -206,10 +224,25 @@
 ! Reference values for choice_ads = 2
 real*8, parameter :: Tref = 243.15D0
-real*8, parameter :: nref = 2.31505631D-6 ! calculated as 18.D-3 / (8.314D0 * Tref) * 0.26D0
-real*8, parameter :: Kd_ref = 0.0206D0
-real*8, parameter :: Ka_ref = 2.4D-4
-real*8, parameter :: Kref = 1.23D7        ! Volcanic tuff @ 243.15 K (obtained at low P / Psat)
-
-logical :: saturation_flag 
+real*8, parameter :: nref = 2.31505631D-6 ! calculated as 18.D-3 / (8.314D0 * Tref) * 0.26D0 ! not used anymore (for the time being)
+real*8, parameter :: Kd_ref = 0.0206D0   ! Not used for the time being (would require specific measurements to be known, but it is rarely measured)
+real*8, parameter :: Ka_ref = 2.4D-4     ! Not used for the time being
+! real*8, parameter :: Kref = 6.27D6 ! Value from data analysis from Pommerol2009 ! Old value 1.23D7        ! Volcanic tuff @ 243.15 K (obtained at low P / Psat) ! ARNAU 
+! real*8, parameter :: Kref2 = 5.95D-7 ! Value from data analysis Pommerol2009 ! ARNAU 
+real*8, parameter :: Kref = 0.205D-6 ! Value obtained from the fit of all adsorption data from Pommmerol (2009) (see Arnau's report - p.28: = yp/xp = 2.6998E-7/3.5562E-2, divided by psat(243.15K=37 Pa; will need to be multiplied by RT/M to be unitless before multiplying by znsoil, which in kg(water)/m3(air) and not in pascals)
+real*8, parameter :: Kref2 = 0.108D-7 ! Value obtained from the fit of all adsorption data from Pommmerol (2009) (see Arnau's report - p.28 = m2; divided by psat(243.15K)=37 Pa) 
+
+
+logical :: recompute_all_cells_ads_flag ! Old saturation_flag but with a behaviour change ! ARNAU
+! The old saturation_flag was used to check whether the cell was saturated or 
+! not, and therefore to assign the saturation value adswater(ig, ik) = 
+! adswater_sat instead of the usual adswater(ig, ik) = adswater_temp(ig, ik) 
+! The old routine using saturation_flag did not require that the A, B, C,... 
+! adsorption coefficients be computed, but the new soilwater subroutine 
+! does. Therefore, the variable recompute_all_cells_ads_flag checks whether
+! there is a cell of a column that requires recomputing. If at least one cell
+! requires recomputing (i.e. recompute_cell_ads_flag(ik) is .true.), then 
+! recompute_all_cells_ads_flag becomes true, and the adsorption coefficients,
+! as well as the recursive equation (appearing in this code as `backward 
+! loop over all layers`) ! ARNAU
 logical :: sublimation_flag
 logical :: condensation_flag(nsoil) 
@@ -243,11 +276,8 @@
 integer, save :: n                        ! number of runs of this subroutine
 integer :: sublim_n                       ! number of sublimation loop runs
-integer :: satur_n                        ! number of saturation loop runs
-
-
+integer :: satur_mono_n                   ! number of monolayer saturation loop runs ! added 2020
 
 
 if (.not. allocated(znsoil)) then
-      print*, 'Flag A'
       allocate( znsoil(ngrid, nsoil) )
       allocate( ice(ngrid, nsoil) )
@@ -271,4 +301,5 @@
       allocate( close_top(ngrid, nsoil) )
       allocate( close_bottom(ngrid, nsoil) )
+      allocate( over_mono_sat_flag(ngrid, nsoil) )
 endif
 
@@ -276,12 +307,8 @@
 ! mugaz ! Molar mass of the atmosphere (g.mol-1) ~43.49
 
-! TODOs
-! right now the subroutine gets called even if there is co2 ice on the surface
-! this should not be the case!! (Note it should be called, but here the assumption is that every cover is water ice)
 
 
 ! Initialisation :
 ! ================= 
-
 
 if (firstcall_soil) then
@@ -289,5 +316,8 @@
       close_top = .false.
       close_bottom = .false.
-      print *, "adsorption enthalpy: ", enthalpy_ads
+      print *, "adsorption enthalpy first segment: ", enthalpy_ads
+      print *, "adsorption enthalpy second segment: ", enthalpy_ads2
+      print *, "adsorption BET constant C first segment: ", DeltaQ_ads
+      print *, "adsorption BET constant C second segment: ", DeltaQ_ads2
       print *, "specific surface area: ", S
 
@@ -315,8 +345,8 @@
                   do ik = 1, nsoil
                         
-                        ! These properties are defined here in order to enable custum profiles
+                        ! These properties are defined here in order to enable custom profiles
                         porosity_ice_free(ig, ik) = 0.45D0
                         tortuosity(ig, ik) = 1.5D0
-                        rho_soil(ig, ik) = 1.3D3
+                        rho_soil(ig, ik) = 1.3D3 ! in kg/m3 of regolith (incl. porosity)
 
                         meshsize(ig, ik) = 5.D-6  ! Characteristic size 1e - 6 = 1 micron : diameter(mode 5) = 10 microns, diameter(mode 1) = 1 microns
@@ -327,10 +357,11 @@
 
                   ! Choice of adsorption coefficients
-                  ! ===================================                   
-                  adswater_sat = 0.8D-2 * S / 13.7D3 * rho_soil(ig, 1)  ! Experimental value for volcanic tuff (Pommerol et al., 2009)
+                  ! ===================================
+                    adswater_sat_mono = 2.6998D-7 * S * rho_soil(ig,1)  ! Unit = kg/m3; From best fit of all adsoprtion data from Pommerol et al. (2009) - See Arnau's report                   
+                  ! adswater_sat_mono = 0.8D-2 * S / 13.7D3 * rho_soil(ig, 1)  ! Unit = kg/m3 ; Experimental value for volcanic tuff (Pommerol et al., 2009) 
                   ! theoretical formula is = rho_soil(ig, 1) * S / Sm * mw
 
-                  ! with SSA = 13.7 m2 / g and plateau = 0.8wt%
-                  !                       adswater_sat = 6D-2 * rho_soil(ig, 1)            ! Experimental value for JSC Mars - 1 (Pommerol et al., 2009)
+                  ! Numerical values above are for SSA = 13.7 m2 / g and plateau = 0.8wt%
+                  !                       adswater_sat_mono = 6D-2 * rho_soil(ig, 1)            ! Experimental value for JSC Mars - 1 (Pommerol et al., 2009)
                   ! with SSA = 106 m2 / g and plateau (net) = 6wt%
 
@@ -340,53 +371,68 @@
                         ! Initialisation of pqsoil(t = 0)
                         ! in 1D simulations the initialization happens here
-                        if (ngrid.eq.1) then
-                              znsoil(ig, ik) = mmr_h2o * mugaz * 1.D-3 * pplev(ig, 1) &
-                              / (8.314D0 * tsoil(ig, nsoil - 4))
-                              !                                   znsoil(ig, ik) = 0.
-                              ice(ig, ik) = 0.D0  ! 0.1D0 !414.D0
-                              
-                              saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
-                              porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
-
-                              if (choice_ads.eq.1) then
-                                    vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * tsoil(ig, nsoil - 4) &
-                                          / (pi * 18.D-3))
-                                    k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
-                                          / (dexp(-enthalpy_ads &
-                                          / (8.314D0 * tsoil(ig, nsoil - 4))) / tau0)
-                                    Kd(ig, ik) = kinetic_factor  &
-                                          / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
-                                    Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) /  &
-                                          (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
-
-                              elseif (choice_ads.eq.2) then  ! par rapport \E0 valeur experimentale de reference
-                                    !                                         Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 *
-                                    !     &                                         (1. / tsoil(ig, nsoil - 4) - 1. / Tref))
-                                    !                                         Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref *
-                                    !     &                                         sqrt(tsoil(ig, nsoil - 4) / Tref) / nref
-
-                                    k_ads_eq(ig, ik) = Kref * dexp(enthalpy_ads / 8.314D0 * (1.D0 / tsoil(ig, ik) - 1.D0 / Tref))
-                                          
-                                    Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
-                                          
-                                    Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
-        
-                              endif
-
-                              if(choice_ads .ne. 0) adswater(ig, ik) = min(Ka(ig, ik) / Kd(ig, ik) * znsoil(ig, ik), adswater_sat)
-                              
-                        else  ! in 3D simulations initialisation happens with newstart.F
+!                       if (ngrid.eq.1) then
+!                             znsoil(ig, ik) = mmr_h2o * mugaz * 1.D-3 * pplev(ig, 1) &
+!                              / (8.314D0 * ptsoil(ig, nsoil - 4))
+!                              !                                   znsoil(ig, ik) = 0.
+!                              ice(ig, ik) = 0.D0  ! 0.1D0 !414.D0
+                              
+!                              saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
+!                              porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
+                              
+!                              if (choice_ads.eq.1) then
+!                                    vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, nsoil - 4) &
+!                                          / (pi * 18.D-3))
+
+                                    ! first segment of the bilinear function 
+!                                    k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
+!                                          / (dexp(-enthalpy_ads &
+!                                          / (8.314D0 * ptsoil(ig, nsoil - 4))) / tau0)
+!                                    Kd(ig, ik) = kinetic_factor  &
+!                                          / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
+!                                    Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) /  &
+!                                          (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
+
+                                    ! second segment of the bilinear function ! added 2020
+!                                    k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
+!                                          / (dexp(-enthalpy_ads2 &
+!                                          / (8.314D0 * ptsoil(ig, nsoil - 4))) / tau0) 
+!                                    Kd2(ig, ik) = kinetic_factor  &
+!                                          / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) 
+!                                    Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) /  &
+!                                          (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) 
+!
+!                              elseif (choice_ads.eq.2) then  ! par rapport \E0 valeur experimentale de reference
+!                                    !                                         Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 *
+!                                    !     &                                         (1. / ptsoil(ig, nsoil - 4) - 1. / Tref))
+!                                    !                                         Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref *
+!                                    !     &                                         sqrt(ptsoil(ig, nsoil - 4) / Tref) / nref!
+!
+!                                    ! first segment of the bilinear function  
+!                                    k_ads_eq(ig, ik) = 8.314D0 * ptsoil(ig,ik)/18.D-3 * Kref * dexp(DeltaQ_ads / 8.314D0 * (1.D0 / ptsoil(ig, ik) - 1.D0 / Tref)) ! Changed enthalpy_ads to DeltaQ_ads to ensure correct dependance/behaviour of k_ads_eq with temperature ; prefactor RT/M is to convert Kref in proper units to use with znsoil
+!                                          
+!                                    Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
+!                                          
+!                                    Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
+!
+!                                    ! second segment of the bilinear function ! added 2020
+!                                    k_ads_eq2(ig, ik) = 8.314D0 * ptsoil(ig,ik)/18.D-3 * Kref2 * dexp(DeltaQ_ads2 / 8.314D0 * (1.D0 / ptsoil(ig, ik) - 1.D0 / Tref)) ! Changed enthalpy_ads2 to DeltaQ_ads2 to ensure correct dependance/behaviour of k_ads_eq with temperature
+!                                          
+!                                    Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) 
+!                                          
+!                                    Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) 
+!                              endif
+                               
+!                              if (k_ads_eq(ig,ik) * znsoil(ig, ik) .gt. adswater_sat_mono) then ! modified 2024, after correction of c0 expression
+!                                 c0(ig, ik) = ( 1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono 
+!                                 adswater(ig, ik) = c0(ig,ik) + k_ads_eq2(ig, ik) * znsoil(ig, ik)
+!                              else
+!                                 adswater(ig, ik) = k_ads_eq(ig, ik) * znsoil(ig, ik)
+!                              endif 
+                              
+!                        else  ! in 3D simulations initialisation happens with newstart.F
                               znsoil(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_soil)
                               ice(ig, ik) = pqsoil(ig, ik, igcm_h2o_ice_soil)
                               adswater(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_ads)
-                        endif
-
-                        if (choice_ads.eq.0) then ! no adsorption
-
-                               Ka(:, :) = 0.
-                               Kd(:, :) = 0.
-                               adswater(:,:) = 0.
-
-                        endif
+!                        endif
 
                         saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
@@ -400,4 +446,12 @@
                         
                         porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
+                                                
+                        ! Initialise soil as if all points where below the monolayer saturation level (added 2020) => now depends on value of adswater (modified in 2024) 
+                        if (adswater(ig,ik).gt.adswater_sat_mono) then
+                           over_mono_sat_flag(ig, ik) = .true. !  
+                        else 
+                            over_mono_sat_flag(ig, ik) = .false.   
+                        endif      
+
                   enddo
 !            endif
@@ -455,5 +509,5 @@
       print *, "Kinetic factor: ", kinetic_factor
       if (kinetic_factor.eq.0) then
-            print *, "WARNING: adsorption is switched of"
+            print *, "WARNING: adsorption is switched off"
       endif
       firstcall_soil = .false.
@@ -488,5 +542,5 @@
 
             ! calculate the air density in the first subsurface layer
-            rho(ig) = pplev(ig, 1) / (r * tsoil(ig, 1))
+            rho(ig) = pplev(ig, 1) / (r * ptsoil(ig, 1))
 
             ! calculate diffusion coefficients at mid- and interlayers (with ice content from the previous timestep)
@@ -495,16 +549,16 @@
 
                   ! calculate H2O mean thermal velocity
-                  vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * tsoil(ig, ik) / (pi * 18.D-3))
+                  vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, ik) / (pi * 18.D-3))
 
                   ! The diffusion coefficient is calculated
 
                   ! calculate the H2O - CO2 collision integral (after Mellon and Jakosky, JGR, 1993)
-                  omega(ig, ik) = 1.442D0 - 0.673D0 * dlog(2.516D-3 * tsoil(ig, ik)) &
-                        + 0.252D0 * (dlog(2.516D-3 * tsoil(ig, ik))) ** 2.D0 &
-                        - 0.047D0 * (dlog(2.516D-3 * tsoil(ig, ik))) ** 3.D0 &
-                        + 0.003D0 * (dlog(2.516D-3 * tsoil(ig, ik))) ** 4.D0
+                  omega(ig, ik) = 1.442D0 - 0.673D0 * dlog(2.516D-3 * ptsoil(ig, ik)) &
+                        + 0.252D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 2.D0 &
+                        - 0.047D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 3.D0 &
+                        + 0.003D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 4.D0
 
                   ! calculate the molecular diffusion coefficient
-                  Dm(ig, ik) = 4.867D0 * tsoil(ig, ik) ** (3.D0 / 2.D0) &
+                  Dm(ig, ik) = 4.867D0 * ptsoil(ig, ik) ** (3.D0 / 2.D0) &
                         / (pplev(ig, 1) * omega(ig, ik)) * 1.D-4
 
@@ -532,10 +586,11 @@
 !                       + (mlayer(ik) - layer(ik)) * saturation_water_ice(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1))
                   
-                  ! new diffusion interaction
                   if (close_bottom(ig, ik).or.close_top(ig, ik+1)) then
                         saturation_water_ice_inter(ig, ik) = 0.999D0
                   else
-                        saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))  
+                        saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))  ! new diffusion interaction
                   endif
+                  
+                  saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))  ! new diffusion interaction
                   
                   D_inter(ig, ik) = D0(ig, ik) * (1.D0 - saturation_water_ice_inter(ig, ik)) ** 2.D0 * 1.D0 / (1.D0 / Dm_inter(ig, ik) + 1.D0 / Dk_inter(ig, ik))
@@ -556,7 +611,8 @@
                   if (choice_ads.eq.1) then  ! Constant, uniform diffusion coefficient D0 (prescribed)
 
+                        ! first segment of the bilinear function (before monolayer saturation) 
                         ! calculate the equilibrium adsorption coefficient
                         k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
-                              / (dexp(-enthalpy_ads / (8.314D0 * tsoil(ig, ik))) / tau0)
+                              / (dexp(-enthalpy_ads / (8.314D0 * ptsoil(ig, ik))) / tau0)
 
                         ! calculate the desorption coefficient
@@ -566,5 +622,17 @@
                         Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
 
-                        !                             Kd(ig, ik) = exp(-enthalpy_ads / (8.314 * tsoil(ig, ik)))
+
+                        ! second segment of the bilinear function (after monolayer saturation) ! added 2020
+                        ! calculate the equilibrium adsorption coefficient
+                        k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
+                              / (dexp(-enthalpy_ads2 / (8.314D0 * ptsoil(ig, ik))) / tau0) ! added 2020 
+
+                        ! calculate the desorption coefficient
+                        Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! added 2020
+
+                        ! calculate the absorption coefficient
+                        Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! added 2020
+
+                        !                             Kd(ig, ik) = exp(-enthalpy_ads / (8.314 * ptsoil(ig, ik)))
                         !     &                             / tau0     ! * 1.D-18
                         !                             Ka(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.   ! * 1.D-18
@@ -576,13 +644,14 @@
 !                        endif
 
-                  elseif (choice_ads.eq.2) then
+                  elseif (choice_ads.eq.2) then ! modified 2020 (with DeltaQ instead of Q)
                         !                                   Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 *
-                        !     &                                   (1. / tsoil(ig, ik) - 1. / Tref))
-                        !                                   Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref * sqrt(tsoil(ig, ik) / Tref)
+                        !     &                                   (1. / ptsoil(ig, ik) - 1. / Tref))
+                        !                                   Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref * sqrt(ptsoil(ig, ik) / Tref)
                         !     &                                   / nref
 
+                        ! first segment of the bilinear function (before monolayer saturation)
                         ! calculate the equilibrium adsorption coefficient
-                        k_ads_eq(ig, ik) = Kref * dexp(enthalpy_ads / 8.314D0 *  &
-                              (1.D0 / tsoil(ig, ik) - 1.D0 / Tref))
+                        k_ads_eq(ig, ik) = 8.314D0 * ptsoil(ig,ik)/18.D-3 *Kref * dexp(DeltaQ_ads / 8.314D0 *  &
+                              (1.D0 / ptsoil(ig, ik) - 1.D0 / Tref)) !  Changed enthalpy_ads to DeltaQ_ads to ensure correct dependance/behaviour of k_ads_eq with temperature
 
                         ! calculate the desorption coefficient
@@ -591,20 +660,34 @@
                         ! calculate the absorption coefficient
                         Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
-                  else  ! no ads
-
-                        Kd(ig, ik) = 0.
-
-                        Ka(ig, ik) = 0.
+
+                        ! second segment of the bilinear function (after monolayer saturation) ! added 2020
+                        ! calculate the equilibrium adsorption coefficient
+                        k_ads_eq2(ig, ik) = 8.314D0 * ptsoil(ig,ik)/18.D-3 * Kref2 * dexp(DeltaQ_ads2 / 8.314D0 *  &
+                              (1.D0 / ptsoil(ig, ik) - 1.D0 / Tref)) ! Changed enthalpy_ads2 to DeltaQ_ads2 to ensure correct dependance/behaviour of k_ads_eq with temperature
+
+                        ! calculate the desorption coefficient
+                        Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) 
+
+                        ! calculate the absorption coefficient
+                        Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) 
                   endif
                   
-                  ! calculate the amount of water vapor at adorption saturation
-                 
-                  if (choice_ads.ne.0) nsat(ig, ik) = adswater_sat * Kd(ig, ik) / Ka(ig, ik)
+                  ! calculate the amount of water vapor at monolayer saturation ! modified 2020
+                  nsat(ig, ik) = adswater_sat_mono * Kd(ig, ik) / Ka(ig, ik) ! 
+
+                  ! calculate the intercept of the second line in the bilinear function ! added 2020; corrected 2024 
+                  c0(ig, ik) = ( 1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono
 
                   ! calculate C, E, and F coefficients for later calculations
                   C(ig, ik) = interlayer_dz(ig, ik) / ptimestep
+
+                  ! first segment of the bilinear function (before monolayer saturation)
                   E(ig, ik) = Ka(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik))
                   F(ig, ik) = Kd(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik))
-                  
+
+                  ! second segment of the bilinear function (after monolayer saturation) ! added 2020
+                  E2(ig, ik) = Ka2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik)) 
+                  F2(ig, ik) = Kd2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik)) 
+
                   ! save the old water concentration
                   znsoilprev(ig, ik) = znsoil(ig, ik)
@@ -630,8 +713,4 @@
       if(.not.watercaptag(ig)) then  ! if there is no polar cap
             do ik = 1, nsoil  ! loop over all soil levels
-                  
-                  ! reset loop variables
-                  ads_water_saturation_flag_1(ik) = .false.
-                  ads_water_saturation_flag_2(ik) = .false.
                   if (ice(ig, ik).eq.0.) then
                         ice_index_flag(ik) = .false.
@@ -641,21 +720,28 @@
             enddo
 
-            ! No (adsorption) saturation assumed
+
             do ik = 1, nsoil  ! loop over all soil levels
                   
-                  ! calculate A and B coefficients
-                  A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik)
-                  B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik)
+                  ! calculate A and B coefficients ! modified 2020
+                  if ( .not.over_mono_sat_flag(ig, ik) ) then ! Assume cell below the monolayer saturation
+                        ! calculate A and B coefficients (first segment of the bilinear function)
+                        A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) 
+                        B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) 
+                  else ! Assume cell over the monolayer saturation ! added 2020
+                        ! calculate A and B coefficients (second segment of the bilinear function)
+                        A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) 
+                        B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) &
+                                  + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 
+                  endif
                   
                   ! calculate the saturation pressure
-                  P_sat_soil(ig, ik) = 611.D0 * dexp(22.5D0 * (1.D0 - 273.16D0 / tsoil(ig, ik))) ! maybe take a new expression (Hsublim = 51.1 kJ / mol) ?
+                  P_sat_soil(ig, ik) = 611.D0 * dexp(22.5D0 * (1.D0 - 273.16D0 / ptsoil(ig, ik))) ! maybe take a new expression (Hsublim = 51.1 kJ / mol) ?
                    
                   ! calculate the gas number density at saturation pressure
-                  nsatsoil(ig, ik) = (P_sat_soil(ig, ik) * mw) / (kb * tsoil(ig, ik))
+                  nsatsoil(ig, ik) = (P_sat_soil(ig, ik) * mw) / (kb * ptsoil(ig, ik))
             enddo
 
             ! calculate the alpha, beta, and delta coefficients for the surface layer
             delta0(ig) = pa(ig, 1) + pb(ig, 2) * (1.D0 - pd(ig, 2)) + porosity_ice_free(ig, 1) * pb(ig, 1)
-            
             alpha0(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / delta0(ig)
             
@@ -683,12 +769,12 @@
                   sublimation_flag = .false.
                   
-                  saturation_flag = .true.
-                  satur_n = 0
-                  do while (saturation_flag)
-!                        print *, satur_n
-                        satur_n = satur_n + 1
+                  recompute_all_cells_ads_flag = .true.
+                  satur_mono_n = 0
+                  do while (recompute_all_cells_ads_flag)
+!                        print *, satur_mono_n
+                        satur_mono_n = satur_mono_n + 1
                         
                         ! reset loop flag
-                        saturation_flag = .false.
+                        recompute_all_cells_ads_flag = .false.
                         
                         ! calculate the surface air density
@@ -770,47 +856,83 @@
                         
                         ! calculate the preliminary amount of adsorbed water
-                        adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) &
-                              / (1.D0 + ptimestep * Kd(ig, nsoil))
-                        
+                        if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024
+                            adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) &
+                                 / (1.D0 + ptimestep * Kd(ig, nsoil))
+                        else
+                             adswater_temp(ig, nsoil) = ( Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) + ptimestep * c0(ig,nsoil) * Kd2(ig,nsoil)) &
+                                 / (1.D0 + ptimestep * Kd2(ig, nsoil))
+                        endif                             
+
+
+ 
                         do ik = nsoil-1, 1, -1  ! backward loop over all layers
                                   
                               znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik)
-                              
-                              adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) &
+                             
+                              if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024
+                                adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) &
                                     / (1.D0 + ptimestep * Kd(ig, ik))
+                              else
+                                 adswater_temp(ig, ik) = ( Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) + ptimestep * c0(ig,ik) * Kd2(ig,ik)) &
+                                    / (1.D0 + ptimestep * Kd2(ig, ik))
+                              endif  
                         enddo
                         
-                        ! check for any saturation
+                        ! check if any cell is over monolayer saturation
                         do ik = 1, nsoil  ! loop over all soil levels
-                              if ( (adswater_temp(ig, ik).ge.adswater_sat) &
-                              .and.(.not.ads_water_saturation_flag_2(ik))) then  ! if saturation and there was no saturation previously
-                              
-                                    print *, "NOTICE: adsorption saturation"
+
+                              ! if( (ik.lt.nsoil) .and. (recompute_cell_ads_flag(ik+1) = .true.) ) then ! Make this loop faster as all cells also need to be recomputed if the cell below needs to be recomputed ! ARNAU
+                              !     recompute_cell_ads_flag(ik) = .true.
+                              !     cycle ! Jump loop
+                              ! endif
+                              
+                              ! Change of coefficients ! ARNAU 
+                              recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients 
+                              if ( adswater_temp(ig, ik).ge.adswater_sat_mono ) then 
+
+                                    print *, "NOTICE: over monolayer saturation"
+
+                                    if ( .not.over_mono_sat_flag(ig, ik) ) then
+                                          ! If the cell enters this scope, it 
+                                          ! means that the cell is over the monolayer
+                                          ! saturation after calculation, but the 
+                                          ! coefficients assume it is below. Therefore, 
+                                          ! the cell needs to be recomputed
                                     
-                                    ! set the adsorbed water to saturation level
-                                    adswater(ig, ik) = adswater_sat
+                                          over_mono_sat_flag(ig, ik) = .true. 
+                                          recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients 
+
+                                          ! change the A and B coefficients (change from first segment to second segment of the bilinear function, as we are over the monolayer saturation is_cell_over_monolayer_sat)
+                                          A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) 
+                                          B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) &
+                                                          + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024
+                                    endif
+                              else
+                                    if ( over_mono_sat_flag(ig, ik) ) then
+                                          ! If the cell enters this scope, it 
+                                          ! means that the cell is below the monolayer
+                                          ! saturation after calculation, but the 
+                                          ! coefficients assume it is above. Therefore, 
+                                          ! the cell needs to be recomputed
                                     
-                                    ! raise the layer saturation flags 1
-                                    ads_water_saturation_flag_1(ik) = .true.
-                                    ads_water_saturation_flag_2(ik) = .true.
-                                    
-                                    ! change the A and B coefficients
-                                    A(ig, ik) = 0.D0
-                                    B(ig, ik) = (adswprev(ig, ik) - adswater_sat) * C(ig, ik)
-                                    
-                              ! if there has never been saturation keep the temporary value from earlier
-                              elseif (.not.ads_water_saturation_flag_2(ik)) then
-                                    adswater(ig, ik) = adswater_temp(ig, ik)
+                                          over_mono_sat_flag(ig, ik) = .false. 
+                                          recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients 
+
+                                          ! calculate A and B coefficients (reset to first segment of the bilinear function)
+                                          A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) 
+                                          B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) 
+                                    endif
                               endif
                         enddo
                         
-                        ! raise the saturation flag if any layer has saturated and reset the first layer saturation flag
+                        ! if one cell needs to be recomputed, then all the column is to be recomputed too
                         do ik = 1, nsoil
-                              if (ads_water_saturation_flag_1(ik)) then
-                                    saturation_flag = .true.
-                                    ads_water_saturation_flag_1(ik) = .false.
+                              if ( recompute_cell_ads_flag(ik) ) then
+                                    recompute_all_cells_ads_flag = .true.
+                              else
+                                    adswater(ig, ik) = adswater_temp(ig, ik) ! if no recomputing is needed, store the value (it may be wrong if the cell below needs to be recomputed. It will be correct in the next loop iterations)
                               endif
                         enddo 
-                  enddo  ! end loop while saturation_flag (adsorption saturation)
+                  enddo  ! end loop while recompute_all_cells_ads_flag (monolayer saturation)
                   
                   !? I'm not sure if this should be calculated here again. I have a feeling that ztot1 is meant
@@ -938,5 +1060,4 @@
                   ! calculate the temporty mixing ratio above the surface
                   zq1temp2(ig) = beta0(ig) + alpha0(ig) * znsoil(ig, 1) / rho(ig)
-
                   ! calculate the flux from the subsurface
                   zdqsdifrego(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * (zq1temp2(ig) - znsoil(ig, 1) / rho(ig) )
@@ -944,5 +1065,5 @@
             else
                   ! calculate the flux from the subsurface
-                  zdqsdifrego(ig) = -D_mid(ig, 1) / midlayer_dz(ig, 0) * (qsat_surf(ig) * rhosurf(ig) - znsoil(ig, 1)) ! faut - il faire intervenir la porosite ?
+                  zdqsdifrego(ig) = D_mid(ig, 1) / midlayer_dz(ig, 0) * (znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig)) ! faut - il faire intervenir la porosite ?
                   if ((pqsurf(ig).eq.0.).and.(zdqsdifrego(ig).lt.0.)) then
                         zdqsdifrego(ig) = 0.
@@ -957,10 +1078,10 @@
 
             if ( (.not.exchange(ig)) &
-            .and. ( -(zdqsdifrego(ig) * ptimestep) &
-            .gt.( pqsurf(ig)  + pdqsdifpot(ig) * ptimestep) ) &
-            .and.( (pqsurf(ig) + pdqsdifpot(ig) * ptimestep).gt.0. ) ) then
+            .and. ( (-zdqsdifrego(ig) * ptimestep) &
+            .gt.( pqsurf(ig) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep) ) &
+            .and.( (pqsurf(ig) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep).gt.0. ) ) then
 
                   ! calculate a new flux from the subsurface
-                  zdqsdifrego(ig) = -( pqsurf(ig) + pdqsdifpot(ig) * ptimestep ) / ptimestep
+                  zdqsdifrego(ig) = -( pqsurf(ig) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep ) / ptimestep
                   
 !                  ! check case where there is CO2 ice on the surface but qsurf = 0
@@ -968,17 +1089,20 @@
 !                        zdqsdifrego(ig) = 0.D0  
 !                  endif
-                  
-                  ! calculate the inital A and B coefficients
-                  do ik = 1, nsoil
+
+
+                  ! calculate A and B coefficients ! modified 2020
+                  if ( .not.over_mono_sat_flag(ig, ik) ) then ! Assume cell below the monolayer saturation
+                        ! calculate A and B coefficients (first segment of the bilinear function)
                         A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik)
                         B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik)
-                  enddo
+                  else ! Assume cell over the monolayer saturation ! added 2020
+                        ! calculate A and B coefficients (second segment of the bilinear function)
+                        A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik)
+                        B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) &
+                                  + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024
+                  endif
                   
                   ! initialize all flags for the loop
                   do ik = 1, nsoil
-                  
-                        ! initialize the first and second layer saturation flag
-                        ads_water_saturation_flag_1(ik) = .false.
-                        ads_water_saturation_flag_2(ik) = .false.
                         
                         ! initialize the ice
@@ -1007,10 +1131,10 @@
                         sublimation_flag = .false.
                                                 
-                        ! loop while there is new saturation
-                        saturation_flag = .true.
-                        do while (saturation_flag)
-                              
-                              ! reset the saturation flag
-                              saturation_flag = .false.
+                        ! loop until good values accounting for monolayer saturation are found
+                        recompute_all_cells_ads_flag = .true.
+                        do while (recompute_all_cells_ads_flag)
+                              
+                              ! reset loop flag
+                              recompute_all_cells_ads_flag = .false.
                               
                               ! calculate the Delta, Alpha, and Beta coefficients in the top layer
@@ -1063,44 +1187,77 @@
                               endif
                               
-                              ! calculate the preliminary amount of adsorbed water
-                              adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) &
-                                    / ( 1.D0 + ptimestep * Kd(ig, nsoil) )                                    
-                              
-                              do ik = nsoil - 1, 1, - 1  ! loop from the bottom up
-                              
-                                    znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik)
-                                    
+
+                             ! calculate the preliminary amount of adsorbed water
+                            if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024
+                                adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) &
+                                     / (1.D0 + ptimestep * Kd(ig, nsoil))
+                            else
+                                 adswater_temp(ig, nsoil) = ( Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) + ptimestep * c0(ig,nsoil) * Kd2(ig,nsoil)) &
+                                 / (1.D0 + ptimestep * Kd2(ig, nsoil))
+                            endif
+
+
+
+                            do ik = nsoil-1, 1, -1  ! backward loop over all layers
+    
+                                  znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik)
+    
+                                  if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024
                                     adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) &
-                                          / (1.D0 + ptimestep * Kd(ig, ik))
-                              enddo
+                                        / (1.D0 + ptimestep * Kd(ig, ik))
+                                  else
+                                     adswater_temp(ig, ik) = ( Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) + ptimestep * c0(ig,ik) * Kd2(ig,ik)) &
+                                        / (1.D0 + ptimestep * Kd2(ig, ik))
+                                  endif
+                            enddo
+
                               
                               ! check for any saturation
                               do ik = 1, nsoil
-                                    if (( adswater_temp(ig, ik).ge.adswater_sat ) &
-                                    .and.(.not.ads_water_saturation_flag_2(ik))) then  ! if there is saturation in a new layer
+                                    ! Change of coefficients ! ARNAU
+                                    recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients ! ARNAU
+                                    if ( adswater_temp(ig, ik).gt.adswater_sat_mono ) then 
+                                    
+                                          print *, "NOTICE: over monolayer saturation"
+
+                                          if ( .not.over_mono_sat_flag(ig, ik) ) then
+                                               ! If the cell enters this scope, it 
+                                               ! means that the cell is over the monolayer
+                                               ! saturation after calculation, but the 
+                                               ! coefficients assume it is below. Therefore, 
+                                               ! the cell needs to be recomputed
                                           
-                                          ! set the amount of adsorbed water to the saturation value
-                                          adswater(ig, ik) = adswater_sat
+                                                over_mono_sat_flag(ig, ik) = .true. 
+                                                recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients 
+
+                                                ! change the A and B coefficients (change from first segment to second segment of the bilinear function, as we are over the monolayer saturation is_cell_over_monolayer_sat)
+                                                A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) 
+                                                B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) &
+                                                          + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024
+                                          endif
+                                    else
+                                          if ( over_mono_sat_flag(ig, ik) ) then
+                                               ! If the cell enters this scope, it 
+                                               ! means that the cell is below the monolayer
+                                               ! saturation after calculation, but the 
+                                               ! coefficients assume it is above. Therefore, 
+                                               ! the cell needs to be recomputed
                                           
-                                          ! raise both saturation flags
-                                          ads_water_saturation_flag_1(ik) = .true.                                          
-                                          ads_water_saturation_flag_2(ik) = .true.
-
-                                          ! set new A and B coefficients
-                                          A(ig, ik) = 0.D0
-                                          B(ig, ik) = (adswprev(ig, ik) - adswater_sat) * C(ig, ik)
-
-                                    elseif (.not.ads_water_saturation_flag_2(ik)) then  ! if the layer has not saturated before
-                                          
-                                          ! adopt the preliminary value
-                                          adswater(ig, ik) = adswater_temp(ig, ik)
+                                                over_mono_sat_flag(ig, ik) = .false. 
+                                                recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients 
+
+                                                ! calculate A and B coefficients (reset to first segment of the bilinear function)
+                                                A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) 
+                                                B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) 
+                                          endif
                                     endif
                               enddo
                               
                               ! raise the saturation flag if any layer has saturated and reset the first layer saturation flag
-                              do ik = 1, nsoil
-                                    if (ads_water_saturation_flag_1(ik)) then
-                                          saturation_flag = .true.
-                                          ads_water_saturation_flag_1(ik) = .false.
+                              do ik = 1, nsoil ! modified 2020
+                                    if ( recompute_cell_ads_flag(ik) ) then
+                                          recompute_all_cells_ads_flag = .true.
+                                    else
+                                          adswater(ig, ik) = adswater_temp(ig, ik) ! if no recomputing is needed, store the value (it may be wrong if the cell below needs to be recomputed. It will be correct in the next loop iterations)
                                     endif
                               enddo
@@ -1170,5 +1327,5 @@
                                     endif
 
-                              elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then  ! if there is new ice through condesation
+                              elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then  ! if there is new ice through condensation
                                     if (.not.condensation_flag(ik)) then
                                     ! raise the ice and sublimation flags
@@ -1206,10 +1363,10 @@
             do ik = 1, nsoil - 1
                   if (ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)).gt.choke_fraction) then  ! if the ice is over saturation or choke_fraction
-                        if ( flux(ig, ik - 1).lt.0.) then
+                        if ( flux(ig, ik - 1).gt.0.) then
                               if (.not.close_top(ig, ik).and.print_closure_warnings) then
                                     print *, "NOTICE: closing top of soil layer due to inward flux", ig, ik      
                               endif
                               close_top(ig, ik) = .true.
-                        elseif (flux(ig, ik - 1).gt.0.) then
+                        elseif (flux(ig, ik - 1).lt.0.) then
                               if (close_top(ig, ik).and.print_closure_warnings) then
                                     print *, "NOTICE: opening top of soil layer due to outward flux", ig, ik      
@@ -1218,10 +1375,10 @@
                         endif
                         
-                        if ( flux(ig, ik).gt.0.) then
+                        if ( flux(ig, ik).lt.0.) then
                               if (.not.close_bottom(ig, ik).and.print_closure_warnings) then
                                     print *, "NOTICE: closing bottom of soil layer due to inward flux", ig, ik      
                               endif
                               close_bottom(ig, ik) = .true.
-                        elseif (flux(ig, ik).lt.0.) then
+                        elseif (flux(ig, ik).gt.0.) then
                               if (close_bottom(ig, ik).and.print_closure_warnings) then
                                     print *, "NOTICE: opening bottom of soil layer due to outward flux", ig, ik      
@@ -1273,5 +1430,5 @@
 !                        endif
                   else  ! opening condition that ice content has fallen sufficiently
-                        if ((close_top(ig, ik).or.close_bottom(ig, ik)).and.print_closure_warnings) then
+                        if (close_top(ig, ik).or.close_bottom(ig, ik).and.print_closure_warnings) then
                               print *, "NOTICE: Reopening soillayer due to decrease in ice", ig, ik                              
                         endif
@@ -1298,5 +1455,5 @@
                   
                   ! calculate how close the water vapor content is to saturizing the adsorbed water
-                  if (choice_ads.ne.0) preduite(ig, ik) = znsoil(ig, ik) / nsat(ig, ik)
+                  preduite(ig, ik) = znsoil(ig, ik) / nsat(ig, ik)
                   
                   ! write the results to the return variable
@@ -1306,10 +1463,8 @@
                   
 
-                  if (close_top(ig, ik).and.close_bottom(ig, ik)) then
-                        close_out(ig, ik) = 3
-                  elseif (close_top(ig, ik)) then
-                        close_out(ig, ik) = 2
+                  if (close_top(ig, ik)) then
+                        close_out(ig, ik) = 1
                   elseif (close_bottom(ig, ik)) then
-                        close_out(ig, ik) = 1
+                        close_out(ig, ik) = -1
                   else
                         close_out(ig, ik) = 0  
@@ -1362,5 +1517,4 @@
 enddo
 
-
 if(writeoutput) then
 !print *, "flux shape: ", shape(flux)
@@ -1446,6 +1600,6 @@
 !endif
 endif
-
 RETURN
 END
 
+
