Context Navigation

← Previous Change
Next Change →

libf

Timestamp:

Aug 28, 2025, 11:39:38 PM (8 months ago)

Author:

llange

Message:

Mars PCM
Cleaning of the subroutine 'soilwater'
LL

File:

: 1 edited

trunk/LMDZ.MARS/libf/phymars/soilwater.F90 (modified) (2 diffs)

Legend:

: Unmodified
: Added
: Removed

trunk/LMDZ.MARS/libf/phymars/soilwater.F90

-                      r3726
+                      r3905
+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, saturation_water_ice)
+      use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg, &
+                           ads_const_D, ads_massive_ice
+      use comcstfi_h, only: pi, r
+      use tracer_mod, only: igcm_h2o_vap
+      use surfdat_h, only: watercaptag ! use mis par AP15 essai
+      use geometry_mod, only: cell_area, latitude_deg
+      use write_output_mod, only: write_output
+implicit none
+! =================================================================================================  =
+! Description:
+!
+! This subroutine calculates the profile of water vapor, adsorbed water and ice in the regolith,
+! down to depth(nsoil) ~ 18m. It calculates the flux of water vapor (zdqsdifrego) between the
+! suburface and the atmosphere.
+!
+! Water vapor and adsorbed water are treated as two separate subsurface tracers that are connected
+! by adsorption / desorption coefficients, including an adsorption / desorption kinetic factor.
+!
+! The water adsorption isotherm is linear, with a saturation plateau corresponding to one monolayer
+! 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
+! compared to the saturation vapor pressure, and if the latter is reached, ice is formed, the water vapor
+! pressure is fixed to its saturation value and the adsorption - diffusion equation is solved again.
+! If ice was already present, its sublimation or growth is calculated.
+!
+! This subroutine is called by vdifc.F if the parameters "regolithe" and "water" are set to true in
+! callphys.def. It returns the flux of water vapor that is injected into or out of the regolith,
+! and which serves as a boundary condition to calculate the vertical atmospheric profile of water vapor
+! in vdifc. It also calculates the exchange between the subsurface and surface ice.
+!
+! It requires three subsurface tracers to be defined in traceur.def :
+!    -  h2o_vap_soil (subsurface water vapor)
+!    -  h2o_ice_soil (subsurface ice)
+!    -  h2o_ads_vap  (adsorbed water)
+!
+! Author : Pierre - Yves Meslin (pmeslin@irap.omp.eu)
+!
+! =================================================================================================  =
+! Arguments :
+! ======================================================================
+! Inputs :
+! ----------------------------------------------------------------------
+integer, intent(in) :: ngrid                    ! number of points in the model (all lat and long point in one 1D array)
+integer, intent(in) :: nlayer                   ! number of layers
+integer, intent(in) :: nq                       ! number of tracers
+integer, intent(in) :: nsoil
+integer, intent(in) :: nqsoil
+logical, save :: firstcall_soil = .true.             ! triggers the initialization
+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 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
+real, intent(in) :: pb(ngrid, nlayer)           ! Coefficients zb
+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) :: 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
+! 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)
+! Outputs :
+! ----------------------------------------------------------------------
+real, intent(out) :: zdqsdifrego(ngrid)                     ! Flux from subsurface (positive pointing outwards)
+real, intent(out) :: zq1temp2(ngrid)                        ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg)
+real*8, intent(out) :: saturation_water_ice(ngrid, nsoil)   ! Water pore ice saturation level (formerly Sw)
+! Outputs for the output files
+real*8 :: preduite(ngrid, nsoil)                ! how close the water vapor density is to adsorption saturation
+integer :: exch(ngrid)                          ! translates the logical variable exchange into a -1 or 1
+real*8 :: mass_h2o(ngrid)                       ! Mass of subsurface water column at each point (kg.m - 2) (formerly mh2otot)
+real*8 :: mass_ice(ngrid)                       ! Mass of subsurface ice at each point (kg.m - 2) (formerly micetot)
+real*8 :: mass_h2o_tot                          ! Mass of subsurface water over the whole planet
+real*8 :: mass_ice_tot                          ! Mass of subsurface ice over the whole planet
+real*8 :: nsurf(ngrid)                          ! surface tracer density in kg/m^3
+real*8 :: close_out(ngrid, nsoil)           ! output for close_top and close_bottom
+! Local (saved) variables :
+! ======================================================================
+real*8 :: P_sat_soil(ngrid, nsoil)              ! water saturation pressure of subsurface cells (at miD-layers) (formerly Psatsoil)
+real*8 :: nsatsoil(ngrid, nsoil)                ! gas number density at saturation pressure
+real*8, allocatable, save :: znsoil(:, :)       ! Water vapor soil concentration (kg.m - 3 of pore air)
+real*8 :: znsoilprev(ngrid, nsoil)              ! value of znsoil in the previous timestep
+real*8 :: znsoilprev2(ngrid, nsoil)             ! second variable for the value of znsoil in the previous timestep
+real*8 :: zdqsoil(ngrid, nsoil)                 ! Increase of pqsoil if sublimation of ice
+real*8, allocatable, save :: ice(:, :)          ! Ice content of subsurface cells (at miD-layers) (kg / m3)
+real*8 :: iceprev(ngrid, nsoil)
+logical :: ice_index_flag(nsoil)                ! flag for ice presence
+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, 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 :: 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)
+real*8 :: beta0(ngrid)
+! Adsorbtion/Desorption variables and parameters
+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)
+real*8 :: dztot1(nsoil)
+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)
+real*8 :: rhosurf(ngrid)            ! Surface air density
+! Porosity and tortuosity variables
+real*8, allocatable, save :: porosity_ice_free(:, :)  ! Porosity with no ice present (formerly eps0)
+real*8, allocatable, save :: porosity(:, :)           ! Porosity (formerly eps)
+real*8 :: porosity_prev(ngrid, nsoil)                 ! Porosity from previous timestep
+real*8, allocatable, save :: tortuosity(:, :)         ! Tortuosity factor (formerly tortuo)
+real*8 :: saturation_water_ice_inter(ngrid, nsoil)          ! Water pore ice saturation level at the interlayer
+! Diffussion coefficients
+real*8 :: D_mid(ngrid, nsoil)       ! Midlayer diffusion coefficients
+real*8 :: D_inter(ngrid, nsoil)     ! Interlayer diffusion coefficients (formerly D)
+real*8, allocatable, save :: D0(:, :)     ! Diffusion coefficient prefactor :
+                                          ! If (medium = 1) : D0 = value of D_mid for saturation_water_ice = 0, ie. poro / tortuo * Dm (in m2 / s)
+                                          ! If (medium = 2) : D0 = porosity_tortuosity factor (dimensionless)
+real*8 :: omega(ngrid, nsoil)       ! H2O - CO2 collision integral
+real*8 :: vth(ngrid, nsoil)         ! H2O thermal speed
+real*8, parameter :: Dk0 = 0.459D0  ! Knudsen diffusion coefficient (for saturation_water_ice = 0) Value for a (5 / 1) bidispersed random assembly of spheres
+                                    ! (dimensionless, ie. divided by thermal speed and characteristic meshsize of the porous medium)
+real*8 :: Dk(ngrid, nsoil)          ! Knudsen diffusion coefficient (divided by porosity / tortuosity factor)
+real*8 :: Dk_inter(ngrid, nsoil)          ! Knudsen diffusion coefficient at the interlayer
+real*8 :: Dm(ngrid, nsoil)          ! Molecular diffusion coefficient
+real*8 :: Dm_inter(ngrid, nsoil)          ! Molecular diffusion coefficient at the interlayer
+real*8, parameter :: choke_fraction = 0.8D0  ! fraction of ice that prevents further diffusion
+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 = .true.
+! Coefficients for the Diffusion calculations
+real*8 :: A(ngrid, nsoil)           ! Coefficient in the diffusion formula
+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 (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
+real*8 :: zdelta(ngrid, nsoil - 1)        ! Delta coefficients
+! Distances between layers
+real*8, allocatable, save :: interlayer_dz(:, :)      ! Distance between the interlayer points in m (formerly interdz)
+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 (adsorption) monolayer saturation ! modified 2020
+real*8, allocatable, save :: meshsize(:, :)     ! scaling/dimension factor for the por size
+real*8, allocatable, save :: rho_soil(:, :)     ! Soil density (bulk -  kg / m3) (formerly rhosoil)
+real*8, allocatable, save :: cste_ads(:, :)     ! Prefactor of adsorption coeff. k
+! general constants
+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 = 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
+real*8, parameter:: Sm = 10.6D-20         ! Surface of the water molecule (m2) (only needed in the theoretical formula which is not used right now)
+! 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 ! 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)
+! variables and parameters for the H2O map
+integer, parameter :: n_long_H2O = 66!180          ! number of longitudes of the new H2O map
+integer, parameter :: n_lat_H2O = 50 !87            ! number of latitudes of the new H2O map
+real*8, parameter :: rho_H2O_ice = 920.0D0      ! Ice density (formerly rhoice)
+real :: latH2O(n_lat_H2O*n_long_H2O)            ! Latitude at H2O map points
+real :: longH2O(n_lat_H2O*n_long_H2O)           ! Longitude at H2O map points
+real :: H2O_depth_data(n_lat_H2O*n_long_H2O)    ! depth of the ice layer in in g/cm^2 at H2O map points
+real :: H2O_cover_data(n_lat_H2O*n_long_H2O)    ! H2O content of the cover layer at H2O map points (not used yet)
+real :: dataH2O(n_lat_H2O*n_long_H2O)           ! H2O content of the ice layer at H2O map points
+real :: latH2O_temp(n_lat_H2O*n_long_H2O)       ! intermediate variable
+real :: longH2O_temp(n_lat_H2O*n_long_H2O)      ! intermediate variable
+real :: dataH2O_temp(n_lat_H2O*n_long_H2O)      ! intermediate variable
+real :: H2O_depth_data_temp(n_lat_H2O*n_long_H2O)     ! intermediate variable
+real, allocatable, save :: H2O(:)                            ! H2O map interpolated at GCM grid points (in wt%)
+real, allocatable, save :: H2O_depth(:)                      ! H2O map ice depth interpolated at GCM in g/cm^2
+! variables for the 1D case
+real*8, parameter :: mmr_h2o = 0.6D-4     ! Water vapor mass mixing ratio (for initialization only)
+real*8 :: diff(ngrid, nsoil)              ! difference between total water content and ice + vapor content
+                                          ! (only used for output, inconstistent: should just be adswater)
+real :: var_flux_soil(ngrid, nsoil)       ! Output of the flux between soil layers (kg/m^3/s)
+real :: default_diffcoeff = 4e-4          ! Diffusion coefficient by default (no variations with Temperature and pressure (m/s^2)
+real :: tol_massiveice = 50.              ! Tolerance to account for massive ice (kg/m^3)
+! Loop variables and counters
+integer :: ig, ik, i, j                   ! loop variables
+logical :: output_trigger                 ! used to limit amount of written output
+integer, save :: n                        ! number of runs of this subroutine
+integer :: sublim_n                       ! number of sublimation loop runs
+integer :: satur_mono_n                   ! number of monolayer saturation loop runs ! added 2020
+if (.not. allocated(znsoil)) then
+      allocate( znsoil(ngrid, nsoil) )
+      allocate( ice(ngrid, nsoil) )
+      allocate( adswater(ngrid, nsoil) )
+      allocate( ztot1(ngrid, nsoil) )
+      allocate( porosity_ice_free(ngrid, nsoil) )
+      allocate( porosity(ngrid, nsoil) )
+      allocate( tortuosity(ngrid, nsoil) )
+      allocate( D0(ngrid, nsoil) )
+      allocate( interlayer_dz(ngrid, nsoil) )
+      allocate( midlayer_dz(ngrid, 0:nsoil) )
+!      allocate( zalpha(ngrid, nsoil-1) )
+!      allocate( zbeta(ngrid, nsoil-1) )
+      allocate( zalpha(ngrid, nsoil) )    ! extra entry to match output format
+      allocate( zbeta(ngrid, nsoil) )     ! extra entry to match output format
+      allocate( meshsize(ngrid, nsoil) )
+      allocate( rho_soil(ngrid, nsoil) )
+      allocate( cste_ads(ngrid, nsoil) )
+      allocate( H2O(ngrid) )
+      allocate( H2O_depth(ngrid) )
+      allocate( close_top(ngrid, nsoil) )
+      allocate( close_bottom(ngrid, nsoil) )
+      allocate( over_mono_sat_flag(ngrid, nsoil) )
+endif
+! Used commons
+! mugaz ! Molar mass of the atmosphere (g.mol-1) ~43.49
+! Initialisation :
+! =================
+if (firstcall_soil) then
+      n = 0
+      close_top = .false.
+      close_bottom = .false.
+      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
+      do ig = 1, ngrid
+!            if(.not.watercaptag(ig)) then
+            ! Initialization of soil parameters
+            ! ===================================
+                  ! initialize the midlayer distances
+                  midlayer_dz(ig, 0) = mlayer(0)
+                  do ik = 1, nsoil - 1
+                        midlayer_dz(ig, ik) = mlayer(ik) - mlayer(ik - 1)
+                  enddo
+                  ! initialize the interlayer distances
+                  interlayer_dz(ig, 1) = layer(1)
+                  do ik = 2, nsoil
+                        interlayer_dz(ig, ik) = layer(ik) - layer(ik - 1)
+                  enddo
+                  ! Choice of porous medium and D0
+                  ! ===============================  =
+                  do ik = 1, nsoil
+                        ! These properties are defined here in order to enable custom profiles
+                        if(ads_massive_ice) then
+                              if(pqsoil(ig, ik, igcm_h2o_ice_soil).gt.tol_massiveice) then
+                                    porosity_ice_free(ig, ik) = 0.999999
+                              else
+                                    porosity_ice_free(ig, ik) = porosity_reg
+                              endif
+                        else
+                           porosity_ice_free(ig, ik) = porosity_reg
+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, saturation_water_ice)
+    use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, &
+                         porosity_reg, ads_const_D, ads_massive_ice
+    use comcstfi_h, only: pi, r
+    use tracer_mod, only: igcm_h2o_vap
+    use surfdat_h, only: watercaptag
+    use geometry_mod, only: cell_area, latitude_deg
+    use write_output_mod, only: write_output
+    implicit none
+    ! =================================================================================================
+    ! Description:
+    !
+    ! This subroutine calculates the profile of water vapor, adsorbed water and ice in the regolith,
+    ! down to depth(nsoil) ~ 18m. It calculates the flux of water vapor (zdqsdifrego) between the
+    ! subsurface and the atmosphere.
+    !
+    ! Water vapor and adsorbed water are treated as two separate subsurface tracers that are connected
+    ! by adsorption / desorption coefficients, including an adsorption / desorption kinetic factor.
+    !
+    ! The water adsorption isotherm is linear, with a saturation plateau corresponding to one monolayer
+    ! 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
+    ! compared to the saturation vapor pressure, and if the latter is reached, ice is formed, the water vapor
+    ! pressure is fixed to its saturation value and the adsorption-diffusion equation is solved again.
+    ! If ice was already present, its sublimation or growth is calculated.
+    !
+    ! This subroutine is called by vdifc.F if the parameters "adsorption_soil" and "water" are set to true in
+    ! callphys.def. It returns the flux of water vapor that is injected into or out of the regolith,
+    ! and which serves as a boundary condition to calculate the vertical atmospheric profile of water vapor
+    ! in vdifc. It also calculates the exchange between the subsurface and surface ice.
+    !
+    ! It requires three subsurface tracers to be defined in traceur.def:
+    !    - h2o_vap_soil (subsurface water vapor)
+    !    - h2o_ice_soil (subsurface ice)
+    !    - h2o_ads_vap  (adsorbed water)
+    !
+    ! Authors: Pierre-Yves Meslin (pmeslin@irap.omp.eu), cleaned by Lucas Lange
+    ! For previous version (with numerous commented lines), see -r 3904
+    ! =================================================================================================
+    ! Arguments
+    ! ======================================================================
+    ! Inputs
+    integer, intent(in) :: ngrid                    ! number of points in the model (all lat and long point in one 1D array)
+    integer, intent(in) :: nlayer                   ! number of layers
+    integer, intent(in) :: nq                       ! number of tracers
+    integer, intent(in) :: nsoil
+    integer, intent(in) :: nqsoil
+    logical, save :: firstcall_soil = .true.        ! triggers the initialization
+    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 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
+    real, intent(in) :: pb(ngrid, nlayer)           ! coefficients zb
+    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) :: 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              ! check we are at the last subtimestep for output
+    ! 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)
+    ! Outputs
+    real, intent(out) :: zdqsdifrego(ngrid)         ! flux from subsurface (positive pointing outwards)
+    real, intent(out) :: zq1temp2(ngrid)            ! temporary atmospheric mixing ratio after exchange with subsurface (kg/kg)
+    real*8, intent(out) :: saturation_water_ice(ngrid, nsoil) ! water pore ice saturation level (formerly Sw)
+    ! Outputs for the output files
+    real*8 :: preduite(ngrid, nsoil)                ! how close the water vapor density is to adsorption saturation
+    integer :: exch(ngrid)                          ! translates the logical variable exchange into a -1 or 1
+    real*8 :: mass_h2o(ngrid)                       ! mass of subsurface water column at each point (kg.m-2) (formerly mh2otot)
+    real*8 :: mass_ice(ngrid)                       ! mass of subsurface ice at each point (kg.m-2) (formerly micetot)
+    real*8 :: mass_h2o_tot                          ! mass of subsurface water over the whole planet
+    real*8 :: mass_ice_tot                          ! mass of subsurface ice over the whole planet
+    real*8 :: nsurf(ngrid)                          ! surface tracer density in kg/m^3
+    real*8 :: close_out(ngrid, nsoil)               ! output for close_top and close_bottom
+    ! Local (saved) variables
+    ! ======================================================================
+    real*8 :: P_sat_soil(ngrid, nsoil)              ! water saturation pressure of subsurface cells (at mid-layers) (formerly Psatsoil)
+    real*8 :: nsatsoil(ngrid, nsoil)                ! gas number density at saturation pressure
+    real*8, allocatable, save :: znsoil(:, :)       ! water vapor soil concentration (kg.m-3 of pore air)
+    real*8 :: znsoilprev(ngrid, nsoil)              ! value of znsoil in the previous timestep
+    real*8 :: znsoilprev2(ngrid, nsoil)             ! second variable for the value of znsoil in the previous timestep
+    real*8 :: zdqsoil(ngrid, nsoil)                 ! increase of pqsoil if sublimation of ice
+    real*8, allocatable, save :: ice(:, :)          ! ice content of subsurface cells (at mid-layers) (kg/m3)
+    real*8 :: iceprev(ngrid, nsoil)
+    logical :: ice_index_flag(nsoil)                ! flag for ice presence
+    real*8, allocatable, save :: adswater(:, :)     ! adsorbed water in subsurface cells (at mid-layers)
+    real*8 :: adswater_temp(ngrid, nsoil)           ! temporary variable for adsorbed water
+    logical, allocatable, save :: over_mono_sat_flag(:, :) ! flag for cells above monolayer saturation
+    real*8 :: adswprev(ngrid, nsoil)
+    logical :: recompute_cell_ads_flag(nsoil)       ! flag to check whether coefficients have changed and need recomputing
+    real*8, save :: adswater_sat_mono                ! adsorption monolayer saturation value (kg.m-3 of regolith)
+    real*8 :: delta0(ngrid)                         ! coefficient delta(0) modified
+    real*8 :: alpha0(ngrid)
+    real*8 :: beta0(ngrid)
+    ! Adsorption/Desorption variables and parameters
+    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
+    real*8 :: Ka2(ngrid, nsoil)                     ! adsorption time constant (s-1) after monolayer saturation (second segment of the bilinear function)
+    real*8 :: Kd2(ngrid, nsoil)                     ! desorption time constant (s-1) after monolayer saturation (second segment of the bilinear function)
+    real*8 :: k_ads_eq2(ngrid, nsoil)               ! equilibrium adsorption coefficient after monolayer saturation (second segment of the bilinear function)
+    real*8 :: c0(ngrid, nsoil)                      ! intercept of the second line in the bilinear function
+    real*8, parameter :: kinetic_factor = 0.01      ! (inverse of) characteristic time to reach adsorption equilibrium (s-1):
+                                                     ! fixed arbitrarily when kinetics factors are unknown
+    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 :: flux(ngrid, 0:nsoil-1)                ! fluxes at interfaces (kg.m-2.s-1) (positive = upward)
+    real*8 :: rho(ngrid)                            ! air density (first subsurface layer)
+    real*8 :: rhosurf(ngrid)                        ! surface air density
+    ! Porosity and tortuosity variables
+    real*8, allocatable, save :: porosity_ice_free(:, :)  ! porosity with no ice present (formerly eps0)
+    real*8, allocatable, save :: porosity(:, :)           ! porosity (formerly eps)
+    real*8 :: porosity_prev(ngrid, nsoil)                 ! porosity from previous timestep
+    real*8, allocatable, save :: tortuosity(:, :)         ! tortuosity factor (formerly tortuo)
+    real*8 :: saturation_water_ice_inter(ngrid, nsoil)    ! water pore ice saturation level at the interlayer
+    ! Diffusion coefficients
+    real*8 :: D_mid(ngrid, nsoil)                   ! midlayer diffusion coefficients
+    real*8 :: D_inter(ngrid, nsoil)                 ! interlayer diffusion coefficients (formerly D)
+    real*8, allocatable, save :: D0(:, :)           ! diffusion coefficient prefactor
+    real*8 :: omega(ngrid, nsoil)                   ! H2O-CO2 collision integral
+    real*8 :: vth(ngrid, nsoil)                     ! H2O thermal speed
+    real*8, parameter :: Dk0 = 0.459D0              ! Knudsen diffusion coefficient (for saturation_water_ice = 0)
+    real*8 :: Dk(ngrid, nsoil)                      ! Knudsen diffusion coefficient (divided by porosity/tortuosity factor)
+    real*8 :: Dk_inter(ngrid, nsoil)                ! Knudsen diffusion coefficient at the interlayer
+    real*8 :: Dm(ngrid, nsoil)                      ! molecular diffusion coefficient
+    real*8 :: Dm_inter(ngrid, nsoil)                ! molecular diffusion coefficient at the interlayer
+    real*8, parameter :: choke_fraction = 0.8D0     ! fraction of ice that prevents further diffusion
+    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 rises over saturation
+    logical, parameter :: print_closure_warnings = .true.
+    ! Coefficients for the diffusion calculations
+    real*8 :: A(ngrid, nsoil)                       ! coefficient in the diffusion formula
+    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 (before monolayer saturation)
+    real*8 :: F(ngrid, nsoil)                       ! coefficient in the diffusion formula (before monolayer saturation)
+    real*8 :: E2(ngrid, nsoil)                      ! coefficient in the diffusion formula (after monolayer saturation)
+    real*8 :: F2(ngrid, nsoil)                      ! coefficient in the diffusion formula (after monolayer saturation)
+    real*8, allocatable, save :: zalpha(:, :)       ! alpha coefficients
+    real*8, allocatable, save :: zbeta(:, :)        ! beta coefficients
+    real*8 :: zdelta(ngrid, nsoil-1)                ! delta coefficients
+    ! Distances between layers
+    real*8, allocatable, save :: interlayer_dz(:, :)      ! distance between the interlayer points in m (formerly interdz)
+    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 (adsorption) monolayer saturation
+    real*8, allocatable, save :: meshsize(:, :)     ! scaling/dimension factor for the pore size
+    real*8, allocatable, save :: rho_soil(:, :)     ! soil density (bulk - kg/m3) (formerly rhosoil)
+    real*8, allocatable, save :: cste_ads(:, :)     ! prefactor of adsorption coeff. k
+    ! General constants
+    real*8, parameter :: kb = 1.38065D-23           ! Boltzmann constant
+    real*8, parameter :: mw = 2.988D-26             ! water molecular weight
+    ! Adsorption coefficients
+    real*8, parameter :: enthalpy_ads = 35.D3       ! enthalpy of adsorption (J.mol-1) options 21.D3, 35.D3, and 60.D3
+    real*8, parameter :: enthalpy_ads2 = 21.D3      ! enthalpy of adsorption (J.mol-1) for the second segment
+    real*8, parameter :: DeltaQ_ads = 21.D3         ! heat of adsorption - enthalpy of liquefaction/sublimation for the first segment (J.mol-1)
+    real*8, parameter :: DeltaQ_ads2 = 21.D3        ! heat of adsorption - enthalpy of liquefaction/sublimation for the second segment (J.mol-1)
+    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
+    real*8, parameter:: Sm = 10.6D-20               ! surface of the water molecule (m2)
+    ! Reference values for choice_ads = 2
+    real*8, parameter :: Tref = 243.15D0
+    real*8, parameter :: nref = 2.31505631D-6       ! not used anymore (for the time being)
+    real*8, parameter :: Kd_ref = 0.0206D0          ! not used for the time being
+    real*8, parameter :: Ka_ref = 2.4D-4            ! not used for the time being
+    real*8, parameter :: Kref = 0.205D-6            ! value obtained from the fit of all adsorption data from Pommerol (2009)
+    real*8, parameter :: Kref2 = 0.108D-7           ! value obtained from the fit of all adsorption data from Pommerol (2009)
+    logical :: recompute_all_cells_ads_flag          ! flag to check whether there is a cell of a column that requires recomputing
+    logical :: sublimation_flag
+    logical :: condensation_flag(nsoil)
+    ! Variables and parameters for the H2O map
+    integer, parameter :: n_long_H2O = 66           ! number of longitudes of the new H2O map
+    integer, parameter :: n_lat_H2O = 50            ! number of latitudes of the new H2O map
+    real*8, parameter :: rho_H2O_ice = 920.0D0      ! ice density (formerly rhoice)
+    real :: latH2O(n_lat_H2O*n_long_H2O)            ! latitude at H2O map points
+    real :: longH2O(n_lat_H2O*n_long_H2O)           ! longitude at H2O map points
+    real :: H2O_depth_data(n_lat_H2O*n_long_H2O)    ! depth of the ice layer in g/cm^2 at H2O map points
+    real :: H2O_cover_data(n_lat_H2O*n_long_H2O)    ! H2O content of the cover layer at H2O map points (not used yet)
+    real :: dataH2O(n_lat_H2O*n_long_H2O)           ! H2O content of the ice layer at H2O map points
+    real :: latH2O_temp(n_lat_H2O*n_long_H2O)       ! intermediate variable
+    real :: longH2O_temp(n_lat_H2O*n_long_H2O)      ! intermediate variable
+    real :: dataH2O_temp(n_lat_H2O*n_long_H2O)      ! intermediate variable
+    real :: H2O_depth_data_temp(n_lat_H2O*n_long_H2O) ! intermediate variable
+    real, allocatable, save :: H2O(:)               ! H2O map interpolated at GCM grid points (in wt%)
+    real, allocatable, save :: H2O_depth(:)         ! H2O map ice depth interpolated at GCM in g/cm^2
+    ! Variables for the 1D case
+    real*8, parameter :: mmr_h2o = 0.6D-4           ! water vapor mass mixing ratio (for initialization only)
+    real*8 :: diff(ngrid, nsoil)                    ! difference between total water content and ice + vapor content
+    real :: var_flux_soil(ngrid, nsoil)             ! output of the flux between soil layers (kg/m^3/s)
+    real :: default_diffcoeff = 4e-4                ! diffusion coefficient by default (no variations with Temperature and pressure (m/s^2)
+    real :: tol_massiveice = 50.                    ! tolerance to account for massive ice (kg/m^3)
+    ! Loop variables and counters
+    integer :: ig, ik, i, j                         ! loop variables
+    logical :: output_trigger                       ! used to limit amount of written output
+    integer, save :: n                              ! number of runs of this subroutine
+    integer :: sublim_n                             ! number of sublimation loop runs
+    integer :: satur_mono_n                         ! number of monolayer saturation loop runs
+    ! Allocate arrays if not already done
+    if (.not. allocated(znsoil)) then
+        allocate(znsoil(ngrid, nsoil))
+        allocate(ice(ngrid, nsoil))
+        allocate(adswater(ngrid, nsoil))
+        allocate(ztot1(ngrid, nsoil))
+        allocate(porosity_ice_free(ngrid, nsoil))
+        allocate(porosity(ngrid, nsoil))
+        allocate(tortuosity(ngrid, nsoil))
+        allocate(D0(ngrid, nsoil))
+        allocate(interlayer_dz(ngrid, nsoil))
+        allocate(midlayer_dz(ngrid, 0:nsoil))
+        allocate(zalpha(ngrid, nsoil))       ! extra entry to match output format
+        allocate(zbeta(ngrid, nsoil))        ! extra entry to match output format
+        allocate(meshsize(ngrid, nsoil))
+        allocate(rho_soil(ngrid, nsoil))
+        allocate(cste_ads(ngrid, nsoil))
+        allocate(H2O(ngrid))
+        allocate(H2O_depth(ngrid))
+        allocate(close_top(ngrid, nsoil))
+        allocate(close_bottom(ngrid, nsoil))
+        allocate(over_mono_sat_flag(ngrid, nsoil))
+    endif
+    ! Initialization
+    ! ================
+    if (firstcall_soil) then
+        n = 0
+        close_top = .false.
+        close_bottom = .false.
+        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
+        do ig = 1, ngrid
+            ! Initialize soil parameters
+            ! ===========================
+            ! Initialize the midlayer distances
+            midlayer_dz(ig, 0) = mlayer(0)
+            do ik = 1, nsoil - 1
+                midlayer_dz(ig, ik) = mlayer(ik) - mlayer(ik - 1)
+            enddo
+            ! Initialize the interlayer distances
+            interlayer_dz(ig, 1) = layer(1)
+            do ik = 2, nsoil
+                interlayer_dz(ig, ik) = layer(ik) - layer(ik - 1)
+            enddo
+            ! Choice of porous medium and D0
+            ! ===============================
+            do ik = 1, nsoil
+                ! These properties are defined here in order to enable custom profiles
+                if (ads_massive_ice) then
+                    if (pqsoil(ig, ik, igcm_h2o_ice_soil) .gt. tol_massiveice) then
+                        porosity_ice_free(ig, ik) = 0.999999
+                    else
+                        porosity_ice_free(ig, ik) = porosity_reg
+                    endif
+                else
+                    porosity_ice_free(ig, ik) = porosity_reg
+                endif
+                tortuosity(ig, ik) = 1.5D0
+                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
+                D0(ig, ik) = porosity_ice_free(ig, ik) / tortuosity(ig, ik)
+            enddo
+            ! Choice of adsorption coefficients
+            ! ==================================
+            ! Unit = kg/m3; From best fit of all adsorption data from Pommerol et al. (2009)
+            adswater_sat_mono = 2.6998D-7 * S * rho_soil(ig, 1)
+            ! Initialization of ice, water vapor and adsorbed water profiles
+            ! ===============================================================
+            do ik = 1, nsoil
+                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)
+                saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
+                if (saturation_water_ice(ig, ik) .lt. 0.D0) then
+                    print *, "WARNING!! soil water ice negative at ", ig, ik
+                    print *, "saturation value: ", saturation_water_ice(ig, ik)
+                    print *, "setting saturation to 0"
+                    saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0)
+                endif
+                porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
+                ! Initialize soil as if all points are below the monolayer saturation level
+                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
+        enddo
+        print *, "initializing H2O data"
+        do ig = 1, ngrid
+            output_trigger = .true.
+            do ik = 1, nsoil
+                saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
+                saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0)
+                porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
+            enddo
+        enddo
+        print *, "Kinetic factor: ", kinetic_factor
+        if (kinetic_factor .eq. 0) then
+            print *, "WARNING: adsorption is switched off"
+        endif
+        firstcall_soil = .false.
+    endif  ! of "if firstcall_soil"
+    ! Update properties in case of the sublimation of massive ice
+    if (ads_massive_ice) then
+        do ig = 1, ngrid
+            do ik = 1, nsoil
+                if (pqsoil(ig, ik, igcm_h2o_ice_soil) .gt. tol_massiveice) then
+                    porosity_ice_free(ig, ik) = 0.999999
+                else
+                    porosity_ice_free(ig, ik) = porosity_reg
+                endif
+            enddo
+        enddo
+    endif
+! Computation of new (new step) transport coefficients
+! ===================================================
+do ig = 1, ngrid  ! loop over all gridpoints
+    if (.not. watercaptag(ig)) then  ! if there is no polar cap
+        do ik = 1, nsoil
+            ! Calculate water saturation level (pore volume filled by ice / total pore volume)
+            saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D0)
+            if (saturation_water_ice(ig, ik) < 0.D0) then
+                print *, "WARNING!! soil water ice negative at ", ig, ik
+                print *, "saturation value: ", saturation_water_ice(ig, ik)
+                print *, "setting saturation to 0"
+                saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0)
+            endif
+            ! Save the old porosity
+            porosity_prev(ig, ik) = porosity(ig, ik)
+            ! Calculate the new porosity
+            porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
+            ! Note: saturation_water_ice and porosity are computed from the ice content of the previous timestep
+        enddo
+        ! Calculate the air density in the first subsurface layer
+        rho(ig) = pplev(ig, 1) / (r * ptsoil(ig, 1))
+        ! Calculate diffusion coefficients at mid- and interlayers (with ice content from the previous timestep)
+        ! Dependence on saturation_water_ice taken from Hudson et al. (2009) and Meslin et al. (SSSAJ, 2010)
+        do ik = 1, nsoil  ! loop over all soil levels
+            ! Calculate H2O mean thermal velocity
+            vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, ik) / (pi * 18.D-3))
+            ! Calculate the H2O - CO2 collision integral (after Mellon and Jakosky, JGR, 1993)
+            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 * ptsoil(ig, ik) ** (3.D0 / 2.D0) &
+                       / (pplev(ig, 1) * omega(ig, ik)) * 1.D-4
+            ! Calculate the Knudsen diffusion coefficient (divided by D0 = porosity / tortuosity)
+            Dk(ig, ik) = Dk0 / D0(ig, ik) * meshsize(ig, ik) * vth(ig, ik)
+            ! Calculate midlayer diffusion coefficient (with dependence on saturation from Meslin et al., 2010)
+            D_mid(ig, ik) = D0(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) ** 2.D0 &
+                          * 1.D0 / (1.D0 / Dm(ig, ik) + 1.D0 / Dk(ig, ik))
+            if (ads_const_D) D_mid(ig, ik) = default_diffcoeff
+        enddo
+        ! Calculate interlayer diffusion coefficient as interpolation of the midlayer coefficients
+        do ik = 1, nsoil - 1
+            Dm_inter(ig, ik) = ((layer(ik) - mlayer(ik - 1)) * Dm(ig, ik + 1) &
+                              + (mlayer(ik) - layer(ik)) * Dm(ig, ik)) / (mlayer(ik) - mlayer(ik - 1))
+            Dk_inter(ig, ik) = ((layer(ik) - mlayer(ik - 1)) * Dk(ig, ik + 1) &
+                              + (mlayer(ik) - layer(ik)) * Dk(ig, ik)) / (mlayer(ik) - mlayer(ik - 1))
+            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))
+            endif
+            saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))
+            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))
+            if (ads_const_D) D_inter(ig, ik) = default_diffcoeff
+        enddo
+    endif
+enddo
+! Computation of new adsorption/desorption constants (Ka, Kd coefficients), and C, E, F coefficients
+! ==================================================================================================
+do ig = 1, ngrid  ! loop over all gridpoints
+    if (.not. watercaptag(ig)) then  ! if there is no polar cap
+        do ik = 1, nsoil  ! loop over all soil levels
+            if (choice_ads == 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 * ptsoil(ig, ik))) / tau0)
+                ! Calculate the desorption coefficient
+                Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
+                ! Calculate the absorption coefficient
+                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 (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)
+                ! 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))
+            elseif (choice_ads == 2) then  ! Modified 2020 (with DeltaQ instead of Q)
+                ! First segment of the bilinear function (before monolayer saturation)
+                ! Calculate the equilibrium adsorption coefficient
+                k_ads_eq(ig, ik) = 8.314D0 * rho_soil(ig, ik) * S * ptsoil(ig, ik) / 18.D-3 * Kref &
+                                 * dexp(DeltaQ_ads / 8.314D0 * (1.D0 / ptsoil(ig, ik) - 1.D0 / Tref))
+                ! Calculate the desorption coefficient
+                Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
+                ! Calculate the absorption coefficient
+                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 (after monolayer saturation) - added 2020
+                ! Calculate the equilibrium adsorption coefficient
+                k_ads_eq2(ig, ik) = 8.314D0 * rho_soil(ig, ik) * S * ptsoil(ig, ik) / 18.D-3 * Kref2 &
+                                  * dexp(DeltaQ_ads2 / 8.314D0 * (1.D0 / ptsoil(ig, ik) - 1.D0 / Tref))
+                ! 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))
+            elseif (choice_ads == 0) then  ! No adsorption/desorption
+                Ka(ig, ik) = 0.D0
+                Kd(ig, ik) = 0.D0
+                Ka2(ig, ik) = 0.D0
+                Kd2(ig, ik) = 0.D0
+            endif
+            ! Calculate the amount of water vapor at monolayer saturation - modified 2020
+            if (choice_ads /= 0) 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)
+            znsoilprev2(ig, ik) = znsoil(ig, ik)  ! needed in "special case" loop
+            adswprev(ig, ik) = adswater(ig, ik)
+            iceprev(ig, ik) = ice(ig, ik)
+            ! Calculate the total ice + water vapor content
+            ztot1(ig, ik) = porosity_prev(ig, ik) * znsoil(ig, ik) + ice(ig, ik)
+        enddo
+    endif
+enddo
+! Computation of delta, alpha and beta coefficients
+! =================================================
+do ig = 1, ngrid  ! loop over all gridpoints
+    ! Initialize the flux to zero to avoid undefined values
+    zdqsdifrego(ig) = 0.D0
+    do ik = 0, nsoil - 1
+        flux(ig, ik) = 0.D0
+    enddo
+    if (.not. watercaptag(ig)) then  ! if there is no polar cap
+        do ik = 1, nsoil  ! loop over all soil levels
+            if (ice(ig, ik) == 0.D0) then
+                ice_index_flag(ik) = .false.
+            else
+                ice_index_flag(ik) = .true.
+            endif
+        enddo
+        do ik = 1, nsoil  ! loop over all soil levels
+            ! 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))
+            endif
+            ! Calculate the saturation pressure
+            P_sat_soil(ig, ik) = 611.D0 * dexp(22.5D0 * (1.D0 - 273.16D0 / ptsoil(ig, ik)))
+            ! Calculate the gas number density at saturation pressure
+            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)
+        beta0(ig) = (pb(ig, 2) * pc(ig, 2) + pq(ig, 1, igcm_h2o_vap) * pa(ig, 1) + pqsurf(ig)) / delta0(ig)
+        ! Set loop flags
+        do ik = 1, nsoil
+            condensation_flag = .false.
+        enddo
+        sublimation_flag = .true.
+        sublim_n = 0
+        do while (sublimation_flag)  ! while there is sublimation
+            sublim_n = sublim_n + 1
+            if (sublim_n > 100) then
+                print *, "sublimation not converging in call ", n
+                print *, "might as well stop now"
+                stop
+            endif
+            ! Reset flag
+            sublimation_flag = .false.
+            recompute_all_cells_ads_flag = .true.
+            satur_mono_n = 0
+            do while (recompute_all_cells_ads_flag)
+                satur_mono_n = satur_mono_n + 1
+                ! Reset loop flag
+                recompute_all_cells_ads_flag = .false.
+                ! Calculate the surface air density
+                rhosurf(ig) = pplev(ig, 1) / (r * ptsrf(ig))
+                ! Calculate the coefficients in the uppermost layer
+                if (exchange(ig)) then  ! if there is surface exchange
+                    ! Calculate the delta, alpha, and beta coefficients
+                    zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) &
+                                  + D_inter(ig, 1) / midlayer_dz(ig, 1) &
+                                  + porosity_ice_free(ig, 1) * pb(ig, 1) / (rho(ig) * ptimestep) * (1.D0 - alpha0(ig))
+                    zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1)
+                    zbeta(ig, 1) = (porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * beta0(ig) &
+                                 + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1)) / zdelta(ig, 1)
+                else
+                    ! Calculate the delta, alpha, and beta coefficients without exchange
+                    zdelta(ig, 1) = A(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1) &
+                                  + D_mid(ig, 1) / midlayer_dz(ig, 0) + porosity(ig, 1) * C(ig, 1)
+                    zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1)
+                    zbeta(ig, 1) = (D_mid(ig, 1) / midlayer_dz(ig, 0) * qsat_surf(ig) * rhosurf(ig) &
+                                 + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1)) / zdelta(ig, 1)
+                endif
+                ! Check for ice
+                if (ice_index_flag(1)) then
+                    ! Set the alpha coefficient to zero
+                    zalpha(ig, 1) = 0.D0
+                    ! Set the beta coefficient to the number density at saturation pressure
+                    zbeta(ig, 1) = nsatsoil(ig, 1)
+                endif
+                do ik = 2, nsoil - 1  ! go through all other layers
+                    ! Calculate delta, alpha, and beta coefficients
+                    zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) &
+                                   + D_inter(ig, ik) / midlayer_dz(ig, ik) &
+                                   + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1))
+                    zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik)
+                    zbeta(ig, ik) = (D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) &
+                                   + B(ig, ik) + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev(ig, ik)) / zdelta(ig, ik)
+                    ! Check for ice
+                    if (ice_index_flag(ik)) then
+                        ! Set the alpha coefficient to zero
+                        zalpha(ig, ik) = 0.D0
+                        ! Set the beta coefficient to the number density at saturation pressure
+                        zbeta(ig, ik) = nsatsoil(ig, ik)
+                    endif
+                enddo
+                ! Computation of pqsoil, ztot1, zq1temp2 and zdqsdifrego
+                ! =======================================================
+                ! Calculate the preliminary amount of water vapor in the bottom layer
+                if (ice_index_flag(nsoil)) then  ! if there is ice
+                    ! Set the vapor number density to saturation
+                    znsoil(ig, nsoil) = nsatsoil(ig, nsoil)
+                else
+                    ! Calculate the vapor number density in the lowest layer
+                    znsoil(ig, nsoil) = (D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) &
+                                       + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev(ig, nsoil) + B(ig, nsoil)) &
+                                      / (porosity(ig, nsoil) * C(ig, nsoil) + A(ig, nsoil) &
+                                       + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)))
+                endif
+                ! 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))
+                    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 if any cell is over monolayer saturation
+                do ik = 1, nsoil  ! loop over all soil levels
+                    ! Change of coefficients
+                    recompute_cell_ads_flag(ik) = .false.  ! Assume there is no change of coefficients
+                    if (adswater_temp(ig, ik) >= 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
+                            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)
+                            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))
                         endif
+                        tortuosity(ig, ik) = 1.5D0
+                        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
+                                                  ! Example : with meshsize = 5.E - 6 : grain sizes = 50 and 10 microns
+                        D0(ig, ik) = porosity_ice_free(ig, ik) / tortuosity(ig, ik)
+                  enddo
+                  ! Choice of adsorption coefficients
+                  ! ===================================
+                    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
+                  ! 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%
+                  ! Initialisation of ice, water vapor and adsorbed water profiles
+                  ! ===============================================================  =
+                  do ik = 1, nsoil
+                        ! 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 * 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
+                        saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
+                        if (saturation_water_ice(ig, ik).lt.0.D0) then
+                              print *, "WARNING!! soil water ice negative at ", ig, ik
+                              print *, "saturation value: ", saturation_water_ice(ig, ik)
+                              print *, "setting saturation to 0"
+                              saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0)
+                    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
+                            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
+                        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
+      enddo
+      print *, "initializing H2O data"
+!
+!      ! 1.6 intitalization of the water content
+!      open(40,file='H2O_data')
+!      do i = 1, n_lat_H2O*n_long_H2O
+!            read(40,*) latH2O_temp(i), longH2O_temp(i), H2O_cover_data(i), dataH2O_temp(i), H2O_depth_data_temp(i)
+!            ! write(*, *) 'depth data ', H2O_depth_data(i)
+!      enddo
+!      close(40)
+!
+!      ! 1.6.2. put latH2O, longH2O and dataH2O in the right format to pass on to mapTh
+!      ! in the datafile the latitudes are listed from negative to positive, but mapTh needs them the other way around
+!      do i = 0, n_lat_H2O - 1
+!            do j = 1, n_long_H2O
+!                  latH2O( n_long_H2O * i + j ) = latH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j)
+!                  longH2O( n_long_H2O * i + j ) = longH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j)
+!                  dataH2O( n_long_H2O * i + j ) = dataH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j)
+!                  H2O_depth_data( n_long_H2O * i + j ) = H2O_depth_data_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j)
+!            enddo
+!      enddo
+!
+!      call mapTh(latH2O, longH2O, dataH2O, n_long_H2O, n_lat_H2O, ngrid, H2O)
+!      call mapTh(latH2O, longH2O, H2O_depth_data, n_long_H2O, n_lat_H2O, ngrid, H2O_depth)
+!
+      do ig = 1, ngrid
+            ! convert depth from g/cm^2 to m
+!            print *, H2O_depth(ig), rho_soil(ig, 1)
+!            H2O_depth(ig) = H2O_depth(ig) * 10 / rho_soil(ig, 1)
+!            if ( (latitude_deg(ig).lt.40).and.(latitude_deg(ig).gt.-40)) then
+!                  H2O(ig) = 0.
+!            endif
+            output_trigger = .true.
+            do ik = 1, nsoil
+!                  if (H2O_depth(ig).le.layer(ik)) then
+!                        ice(ig, ik) = min(H2O(ig), rho_H2O_ice * porosity_ice_free(ig, ik))
+!                        if (output_trigger) then
+!                              print *, "ice set at: ", ig, ik, "to :", ice(ig,ik), "depth: ", H2O_depth(ig)
+!                              output_trigger = .false.
+!                        endif
+!                  endif
+                  saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
+                  saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0)
+                  porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
+                    endif
+                enddo
+                ! If one cell needs to be recomputed, then all the column is to be recomputed too
+                do ik = 1, nsoil
+                    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
+                    endif
+                enddo
+            enddo  ! end loop while recompute_all_cells_ads_flag (monolayer saturation)
+            ! Calculation of Fnet + Fads
+            ! ===========================
+            ! Calculate the flux variable
+            ! Calculate the flux in the top layer
+            if (exchange(ig)) then  ! if there is exchange
+                ! Calculate the flux taking diffusion to the atmosphere into account
+                flux(ig, 0) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep &
+                            * (znsoil(ig, 1) / rho(ig) - beta0(ig) - alpha0(ig) * znsoil(ig, 1) / rho(ig))
+            elseif (pqsurf(ig) > 0.D0) then
+                ! Assume that the surface is covered by water ice
+                flux(ig, 0) = D_mid(ig, 1) / midlayer_dz(ig, 0) * (znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig))
+            endif
+            ! Check if the ground would take up water from the surface but there is none
+            if ((.not. exchange(ig)) .and. (pqsurf(ig) == 0.D0) .and. (flux(ig, 0) < 0.D0)) then
+                flux(ig, 0) = 0.D0
+            endif
+            ! Calculate the flux in all other layers
+            do ik = 1, nsoil - 1
+                flux(ig, ik) = D_inter(ig, ik) * (znsoil(ig, ik + 1) - znsoil(ig, ik)) / midlayer_dz(ig, ik)
             enddo
+      enddo
+      print *, "Kinetic factor: ", kinetic_factor
+      if (kinetic_factor.eq.0) then
+            print *, "WARNING: adsorption is switched off"
+      endif
+      firstcall_soil = .false.
+endif  ! of "if firstcall_soil"
+! Update properties in case of the sublimation of massive ice
+if(ads_massive_ice) then
+      do ig = 1, ngrid
+            do ik = 1, nsoil
+                  if(pqsoil(ig, ik, igcm_h2o_ice_soil).gt.tol_massiveice) then
+                        porosity_ice_free(ig, ik) = 0.999999
+                  else
+                        porosity_ice_free(ig, ik) = porosity_reg
+                  endif
+            ! Calculate dztot1 which describes the change in ztot1 (water vapor and ice)
+            do ik = 1, nsoil - 1
+                dztot1(ik) = (flux(ig, ik) - flux(ig, ik - 1)) / interlayer_dz(ig, ik) &
+                           - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) &
+                           + B(ig, ik) / interlayer_dz(ig, ik)
             enddo
+      enddo
+endif
+! Computation of new (new step) transport coefficients :
+! =======================================================  =
+do ig = 1, ngrid  ! loop over all gridpoints
+      if(.not.watercaptag(ig)) then  ! if there is no polar cap
+            dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) &
+                          - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) &
+                          + B(ig, nsoil) / interlayer_dz(ig, nsoil)
+            ! Condensation / sublimation
+            do ik = 1, nsoil  ! loop over all
+                if (ice_index_flag(ik)) then  ! if there is ice
+                    ! Compute ice content
+                    ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik)
+                    if (ice(ig, ik) < 0.D0) then  ! if all the ice is used up
+                        ! Set the ice content to zero
+                        ice(ig, ik) = 0.D0
+                        ! Change the water vapor from the previous timestep. (Watch out! could go wrong)
+                        znsoilprev(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik)
+                        ! Remove the ice flag and raise the sublimation flag
+                        ice_index_flag(ik) = .false.
+                        sublimation_flag = .true.
+                    endif
+                elseif (znsoil(ig, ik) > nsatsoil(ig, ik)) then  ! if there is condensation
+                    if (.not. condensation_flag(ik)) then
+                        ! Raise the ice and sublimation flags
+                        ice_index_flag(ik) = .true.
+                        sublimation_flag = .true.
+                        condensation_flag(ik) = .true.
+                    endif
+                endif
+                ! Calculate the difference between total water content and ice + vapor content (only used for output)
+                diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) &
+                             + adswater(ig, ik) - ztot1(ig, ik) - dztot1(ik) * ptimestep
+            enddo  ! loop over all layers
+        enddo  ! end first loop while sublimation_flag (condensation / sublimation)
+        if (exchange(ig)) then  ! if there is exchange
+            ! Calculate the temporary 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))
+        else
+            ! Calculate the flux from the subsurface
+            zdqsdifrego(ig) = D_mid(ig, 1) / midlayer_dz(ig, 0) * (znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig))
+            if ((pqsurf(ig) == 0.D0) .and. (zdqsdifrego(ig) < 0.D0)) then
+                zdqsdifrego(ig) = 0.D0
+            endif
+        endif
+        ! Special case where there is not enough ice on the surface to supply the subsurface for the whole timestep
+        ! (when exchange with the atmosphere is disabled): the upper boundary condition becomes a flux condition
+        ! (and not qsat_surf) and all the remaining surface ice is sublimated and transferred to the subsurface.
+        ! The actual change in surface ice happens in a higher subroutine
+        ! =====================================================================================================
+        if ((.not. exchange(ig)) &
+            .and. ((-zdqsdifrego(ig) * ptimestep) > (pqsurf(ig) + pdqsdifpot(ig) * ptimestep)) &
+            .and. ((pqsurf(ig) + pdqsdifpot(ig) * ptimestep) > 0.D0)) then
+            ! Calculate a new flux from the subsurface
+            zdqsdifrego(ig) = -(pqsurf(ig) + pdqsdifpot(ig) * ptimestep) / ptimestep
             do ik = 1, nsoil
+                  ! calculate Water saturation level (pore volume filled by ice / total pore volume)
+                  saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
+                  if (saturation_water_ice(ig, ik).lt.0.D0) then
+                        print *, "WARNING!! soil water ice negative at ", ig, ik
+                        print *, "saturation value: ", saturation_water_ice(ig, ik)
+                        print *, "setting saturation to 0"
+                        saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0)
+                  endif
+                  ! save the old porosity
+                  porosity_prev(ig, ik) = porosity(ig, ik)
+                  ! calculate the new porosity
+                  porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
+                  ! Note : saturation_water_ice and porosity are computed from the ice content of the previous timestep
+                ! 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))
+                endif
+                ! Initialize the ice
+                ice(ig, ik) = iceprev(ig, ik)
+                if (iceprev(ig, ik) == 0.D0) then
+                    ice_index_flag(ik) = .false.
+                else
+                    ice_index_flag(ik) = .true.
+                endif
             enddo
+            ! calculate the air density in the first subsurface layer
+            rho(ig) = pplev(ig, 1) / (r * ptsoil(ig, 1))
+            ! calculate diffusion coefficients at mid- and interlayers (with ice content from the previous timestep)
+            ! Dependence on saturation_water_ice taken from the relation obtained by Hudson et al. (2009) and Meslin et al. (SSSAJ, 2010)
+            do ik = 1, nsoil  ! loop over all soil levels
+                  ! calculate H2O mean thermal velocity
+                  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 * 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 * ptsoil(ig, ik) ** (3.D0 / 2.D0) &
+                        / (pplev(ig, 1) * omega(ig, ik)) * 1.D-4
+                  ! calculate the Knudsen diffusion coefficient (divided by D0 = porosity / tortuosity in this case)
+                  Dk(ig, ik) = Dk0 / D0(ig, ik) * meshsize(ig, ik) * vth(ig, ik)
+                  ! calculate the midlayer diffusion coeff. (with dependence on saturation_water_ice from Meslin et al., 2010 -  only exact for normal diffusion)
+                  D_mid(ig, ik) = D0(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) ** 2.D0 * 1.D0 / (1.D0 / Dm(ig, ik) + 1.D0 / Dk(ig, ik))
+                  if(ads_const_D) D_mid(ig, ik) = default_diffcoeff
+                  ! Midlayer diffusion coeff. (correlation by Rogers and Nielson, 1991)
+                  !                             D_mid(ig, ik) = D0(ig, ik) * (1. - saturation_water_ice(ig, ik)) * exp(-6. * saturation_water_ice(ig, ik)  &
+                  !                                   * porosity_ice_free(ig, ik) - 6. * saturation_water_ice(ig, ik) ** (14. * porosity_ice_free(ig, ik))) &
+                  !                                   * 1. / (1. / Dm(ig, ik) + 1. / Dk(ig, ik))
+            enddo
+            ! calculate the interlayer diffusion coefficient as interpolation of the midlayer coefficients
+            do ik = 1, nsoil - 1
+                  Dm_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * Dm(ig, ik + 1) &
+                        + (mlayer(ik) - layer(ik)) * Dm(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1))
+                  Dk_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * Dk(ig, ik + 1) &
+                        + (mlayer(ik) - layer(ik)) * Dk(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1))
+!                 saturation_water_ice_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * saturation_water_ice(ig, ik + 1) &
+!                       + (mlayer(ik) - layer(ik)) * saturation_water_ice(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1))
+                  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))  ! 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))
+                  if(ads_const_D) D_inter(ig, ik) = default_diffcoeff
+!                  D_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * D_mid(ig, ik + 1) &  ! old implementation with direct interpolation
+!                        + (mlayer(ik) - layer(ik)) * D_mid(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1))
+            enddo
+      endif
+enddo
+! Computation of new adsorption / desorption constants (Ka, Kd coefficients), and C, E, F coefficients
+! ===================================================================================================
+do ig = 1, ngrid  ! loop over all gridpoints
+      if(.not.watercaptag(ig)) then  ! if there is no polar cap
+            do ik = 1, nsoil  ! loop over all soil levels
+                  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 * ptsoil(ig, ik))) / tau0)
+                        ! calculate the desorption coefficient
+                        Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
+                        ! calculate the absorption coefficient
+                        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 (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
+!                        if ((ngrid.eq.1).and.(ik.eq.18)) then  ! if 1D simulation and uppermost layer
+!                              print * , 'Ka, Kd, Ka / Kd = ', Ka(ig, ik), Kd(ig, ik), Ka(ig, ik) / Kd(ig, ik)
+!                              print * , 'k_ads_eq = ', k_ads_eq(ig, ik)
+!                              print * , 'porosity = ', porosity(ig, ik)
+!                        endif
+                  elseif (choice_ads.eq.2) then ! modified 2020 (with DeltaQ instead of Q)
+                        !                                   Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 *
+                        !     &                                   (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) = 8.314D0 * rho_soil(ig, ik) * S * 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
+                        Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
+                        ! calculate the absorption coefficient
+                        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 (after monolayer saturation) ! added 2020
+                        ! calculate the equilibrium adsorption coefficient
+                        k_ads_eq2(ig, ik) = 8.314D0 * rho_soil(ig, ik) * S * 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))
+                  elseif (choice_ads.eq.0)  then!No adsorption/desorption
+                        Ka(ig, ik) = 0.
+                        Kd(ig, ik) = 0.
+                        Ka2(ig, ik) = 0.
+                        Kd2(ig, ik) = 0.
+                  endif
+                  ! calculate the amount of water vapor at monolayer saturation ! modified 2020
+                  if(choice_ads.ne.0) 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)
+                  znsoilprev2(ig, ik) = znsoil(ig, ik)  ! needed in "special case" loop because adswprev can be changed before this loop
+                  adswprev(ig, ik) = adswater(ig, ik)  !!!! Besoin de adswprev2 ???
+                  iceprev(ig, ik) = ice(ig, ik)  ! needed in "special case" loop for the same reason
+                  ! calculate the total ice + water vapor content
+                  ztot1(ig, ik) = porosity_prev(ig, ik) * znsoil(ig, ik) + ice(ig, ik)
+            enddo
+      endif
+enddo
+! Computation of delta, alpha and beta coefficients ETC !!!!
+! ===========================================================
+do ig = 1, ngrid  ! loop over all gridpoints
+      ! initialize the flux to zero to avoid undefined values
+      zdqsdifrego(ig) = 0.D0
+      do ik = 0, nsoil - 1
+            flux(ig, ik) = 0.
+      enddo
+      if(.not.watercaptag(ig)) then  ! if there is no polar cap
+            do ik = 1, nsoil  ! loop over all soil levels
+                  if (ice(ig, ik).eq.0.) then
+                        ice_index_flag(ik) = .false.
+                  else
+                        ice_index_flag(ik) = .true.
+                  endif
+            enddo
+            do ik = 1, nsoil  ! loop over all soil levels
+                  ! 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 / 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 * 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)
+            beta0(ig) = ( pb(ig, 2) * pc(ig, 2) + pq(ig, 1, igcm_h2o_vap) * pa(ig, 1) + pqsurf(ig) ) &
+                  / delta0(ig)
+            ! set loop flag
+            do ik = 1, nsoil
+                  condensation_flag = .false.
+            enddo
+            ! Loop while there is new sublimation
             sublimation_flag = .true.
             sublim_n = 0
+            do while (sublimation_flag)  ! while there is sublimation
+!                  print *, "sublimation loop: ", sublim_n
+                  sublim_n = sublim_n + 1
+                  if (sublim_n.gt.100) then
+                        print *, "sublimation not converging in call ", n
+                        print *, "might as well stop now"
+                        stop
+                  endif
+                  ! reset flag
+                  sublimation_flag = .false.
+                  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
+                        recompute_all_cells_ads_flag = .false.
+                        ! calculate the surface air density
+                        rhosurf(ig) = pplev(ig, 1) / (r * ptsrf(ig))
+                        ! calculate the coefficients in the upermost layer
+                        if (exchange(ig)) then  ! if there is surface exchange
+                              ! calculate the delta, alpha, and beta coefficients
+                              zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) &
+                                    + D_inter(ig, 1) / midlayer_dz(ig, 1) &  !!! est - ce que le terme pb inclut le bon rho ?
+                                    + porosity_ice_free(ig, 1) * pb(ig, 1) / (rho(ig) * ptimestep) * (1.D0 - alpha0(ig))
+                                    !pourrait remplacer pb/(ptime*rho(1/2)) par zcdv/ptime
+                              zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1)
+                              zbeta(ig, 1) = ( porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * beta0(ig) &
+                                    + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1) ) &
+                                    / zdelta(ig, 1)
+            do while (sublimation_flag)
+                sublim_n = sublim_n + 1
+                if (sublim_n > 100) then
+                    print *, "special case sublimation not converging in call ", n
+                    print *, "might as well stop now"
+                    stop
+                endif
+                ! Reset the sublimation flag
+                sublimation_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
+                    zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1)
+                    zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1)
+                    zbeta(ig, 1) = (porosity_prev(ig, 1) * C(ig, 1) * znsoilprev2(ig, 1) + B(ig, 1) - zdqsdifrego(ig)) / zdelta(ig, 1)
+                    ! Check for ice
+                    if (ice_index_flag(1)) then
+                        ! Set the alpha coefficient to zero
+                        zalpha(ig, 1) = 0.D0
+                        ! Set the beta coefficient to the number density at saturation pressure
+                        zbeta(ig, 1) = nsatsoil(ig, 1)
+                    endif
+                    do ik = 2, nsoil - 1  ! loop over all other layers
+                        ! Calculate the Delta, Alpha, and Beta coefficients in the layer
+                        zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) + D_inter(ig, ik) / midlayer_dz(ig, ik) &
+                                       + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1))
+                        zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik)
+                        zbeta(ig, ik) = (D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) &
+                                       + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev2(ig, ik) + B(ig, ik)) / zdelta(ig, ik)
+                        ! Check for ice
+                        if (ice_index_flag(ik)) then
+                            ! Set the alpha coefficient to zero
+                            zalpha(ig, ik) = 0.D0
+                            ! Set the beta coefficient to the number density at saturation pressure
+                            zbeta(ig, ik) = nsatsoil(ig, ik)
+                        endif
+                    enddo
+                    ! Calculate the preliminary amount of water vapor in the bottom layer
+                    if (ice_index_flag(nsoil)) then
+                        ! Set the vapor number density to saturation
+                        znsoil(ig, nsoil) = nsatsoil(ig, nsoil)
+                    else
+                        ! Calculate the vapor number density in the lowest layer
+                        znsoil(ig, nsoil) = (D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) &
+                                           + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev2(ig, nsoil) + B(ig, nsoil)) &
+                                          / (porosity(ig, nsoil) * C(ig, nsoil) &
+                                           + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)) &
+                                           + A(ig, nsoil))
+                    endif
+                    ! 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))
                         else
+                              ! calculate the delta, alpha, and beta coefficients without exchange
+                              zdelta(ig, 1) = A(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1)  &
+                                    + D_mid(ig, 1) / midlayer_dz(ig, 0) &
+                                    + porosity(ig, 1) * C(ig, 1)
+                              zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1. / zdelta(ig, 1)
+                              zbeta(ig, 1) = ( D_mid(ig, 1) / midlayer_dz(ig, 0) * qsat_surf(ig) * rhosurf(ig) &
+                                    + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1) ) &
+                                    / zdelta(ig, 1)
+                            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
+                        ! check for ice
+                        if (ice_index_flag(1)) then
+                              ! set the alpha coefficient to zero
+                              zalpha(ig, 1) = 0.D0
+                              ! set the beta coefficient to the number density at saturation pressure
+                              zbeta(ig, 1) = nsatsoil(ig, 1)
+                    enddo
+                    ! Check for any saturation
+                    do ik = 1, nsoil
+                        ! Change of coefficients
+                        recompute_cell_ads_flag(ik) = .false.  ! Assume there is no change of coefficients
+                        if (adswater_temp(ig, ik) > 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
+                                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
+                                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))
+                            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
+                                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
+                        do ik = 2, nsoil - 1  ! go through all other layers
+                              ! calculate delta, alpha, and beta coefficients
+                              zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) &
+                                    + D_inter(ig, ik) / midlayer_dz(ig, ik) &
+                                    + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1))
+                              zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik)
+                              zbeta(ig, ik) = ( D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) &
+                                    + B(ig, ik) &
+                                    + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev(ig, ik) ) / zdelta(ig, ik)
+                              ! check for ice
+                              if (ice_index_flag(ik)) then
+                                    ! set the alpha coefficient to zero
+                                    zalpha(ig, ik) = 0.D0
+                                    ! set the beta coefficient to the number density at saturation pressure
+                                    zbeta(ig, ik) = nsatsoil(ig, ik)
+                              endif
+                        enddo
+                        ! Computation of pqsoil, ztot1, zq1temp2 and zdqsdifrego
+                        ! =======================================================  =
+                        ! calculate the preliminary amount of water vapor in the bottom layer
+                        if (ice_index_flag(nsoil)) then  ! if there is ice
+                              ! set the vapor number density to saturation
+                              znsoil(ig, nsoil) = nsatsoil(ig, nsoil)
+                    enddo
+                    ! Raise the saturation flag if any layer has saturated
+                    do ik = 1, nsoil
+                        if (recompute_cell_ads_flag(ik)) then
+                            recompute_all_cells_ads_flag = .true.
                         else
+                              ! calculate the vapor number density in the lowest layer
+                              znsoil(ig, nsoil) = ( D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) &
+                                    + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev(ig, nsoil) + B(ig, nsoil) ) &
+                                    / ( porosity(ig, nsoil) * C(ig, nsoil) + A(ig, nsoil) &
+                                    + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)) )
+                            adswater(ig, ik) = adswater_temp(ig, ik)  ! if no recomputing is needed, store the value
                         endif
+                        ! 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))
+                              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 if any cell is over monolayer saturation
+                        do ik = 1, nsoil  ! loop over all soil levels
+                              ! 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
+                                          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
+                                          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
+                        ! if one cell needs to be recomputed, then all the column is to be recomputed too
+                        do ik = 1, nsoil
+                              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 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
+                  ! as the value from the previous timestep
+                  !do ik = 1, nsoil
+                  !      ztot1(ig, ik) = porosity_prev(ig, ik) * znsoil(ig, ik) + ice(ig, ik)
+                  !enddo
+                  ! Calculation of Fnet + Fads
+                  ! =============================
+                  ! calulate the flux variable
+                  ! calculate the flux in the top layer
+                  if (exchange(ig)) then  ! if there is exchange
+                        ! calculate the flux taking diffusion to the atmosphere into account
+                        flux(ig, 0) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep &
+                              * ( znsoil(ig, 1) / rho(ig) - beta0(ig) - alpha0(ig) * znsoil(ig, 1) / rho(ig) )
+                  elseif (pqsurf(ig).gt.0.) then
+                        ! assume that the surface is covered by waterice (if it is co2ice it should not call this subroutine at all)
+                        flux(ig, 0) = D_mid(ig, 1) / midlayer_dz(ig, 0) * ( znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig) )
+                  endif
+                  ! check if the ground would take up water from the surface but there is none
+                  if ((.not.exchange(ig)).and.(pqsurf(ig).eq.0.).and.(flux(ig, 0).lt.0.)) then
+                        flux(ig, 0) = 0.
+                  endif
+                  ! calculate the flux in all other layers
+                  do ik = 1, nsoil - 1
+                        flux(ig, ik) = D_inter(ig, ik) * ( znsoil(ig, ik + 1) - znsoil(ig, ik) ) / midlayer_dz(ig, ik)
+                  enddo
+                  ! calculate dztot1 which descibes the change in ztot1 (water vapor and ice). It is only used for output (directly and indirectly)
+                  ! calculate the change in ztot1
+                  do ik = 1, nsoil - 1
+                        dztot1(ik) = ( flux(ig, ik) - flux(ig, ik - 1) ) / interlayer_dz(ig, ik) &
+                              - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) &
+                              + B(ig, ik) / interlayer_dz(ig, ik)
+                  enddo
+                  dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) &
+                        - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) &
+                        + B(ig, nsoil) / interlayer_dz(ig, nsoil)
+                  ! Condensation / sublimation
+                  do ik = 1, nsoil  ! loop over all layers
+!                        print *, "ice in layer      ", ik, ": ", ice(ig, ik)
+!                        print *, "vapor in layer    ", ik, ": ", znsoil(ig, ik)
+!                        print *, "sat dens in layer ", ik, ": ", nsatsoil(ig, ik)
+!                        print *, ""
+                        if (ice_index_flag(ik)) then  ! if there is ice
+                              ! Compute ice content
+                              ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik)
+                              if (ice(ig, ik).lt.0.D0) then  ! if all the ice is used up
+!                                    print *, "########complete sublimation in layer", ik, "  cell:", ig
+!                                    print *, "ztot1: ", ztot1(ig, ik)
+!                                    print *, "dztot1*timestep: ", dztot1(ik) * ptimestep
+!                                    print *, "vapour: ", porosity(ig, ik) * nsatsoil(ig, ik)
+!                                    print *, "znsoil: ", znsoil(ig, ik)
+!                                    print *, "nsatsoil: ", nsatsoil(ig, ik)
+!                                    print *, "porosity: ", porosity(ig, ik)
+!                                    print *, "ice: ", ice(ig, ik)
+!                                    print *, "exchange: ", exchange(ig)
+!                                    print *, "co2ice: ", co2ice(ig)
+!                                    print *, ""
+!                                    print *, "zalpha: ", zalpha(ig, ik)
+!                                    print *, "zbeta: ", zbeta(ig, ik)
+!                                    print *, ""
+                                    ! set the ice content to zero
+                                    ice(ig, ik) = 0.D0
+                                    ! change the water vapor from the previous timestep. (Watch out! could go wrong)
+                                    znsoilprev(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik)
+!                                    print *, "new znsoilprev", znsoilprev(ig, ik)
+                                    ! remove the ice flag and raise the sublimation flag
+                                    ice_index_flag(ik) = .false.
+                                    sublimation_flag = .true.
+                              endif
+                        elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then  ! if there is condenstation
+                              !ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik)
+!                              print *, "=========== new condensation in layer", ik, "  cell:", ig
+!                              print *, "znsoil, nsatsoil: ", znsoil(ig, ik), nsatsoil(ig, ik)
+!                              print *, "ztot1: ", ztot1(ig, ik)
+!                              print *, "dztot1: ", dztot1(ik)
+!                              print *, "ice: ", ice(ig,ik)
+!                              print *, ""
+!                              print *, "zalpha: ", zalpha(ig, ik)
+!                              print *, "zbeta: ", zbeta(ig, ik)
+!                              print *, ""
+                              !if (ice(ig, ik).lt.0.D0) then
+                                    ! set the ice content to zero
+                              !      ice(ig, ik) = 0.D0
+                              if (.not.condensation_flag(ik)) then
+                                    ! raise the ice and sublimation flags
+                                    ice_index_flag(ik) = .true.
+                                    sublimation_flag = .true.
+!                                    print *, condensation_flag(ik)
+                                    condensation_flag(ik) = .true.
+!                                    print *, condensation_flag(ik)
+!                                    print *, "set condensation flag"
+!                              else
+!                                    print *, "condensation loop supressed"
+                              endif
+                        endif
+                        ! calculate the difference between total water content and ice + vapor content (only used for output)
+                        diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) &
+                              + adswater(ig, ik) - ztot1(ig, ik) - dztot1(ik) * ptimestep
+                  enddo  ! loop over all layer
+            enddo  ! end first loop while sublimation_flag (condensation / sublimation)
+            if (exchange(ig)) then  ! if there is exchange
+                  ! 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) )
+            else
+                  ! calculate the flux from the subsurface
+                  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.
+                  endif
+            endif
+! Special case where there is not enough ice on the surface to supply the subsurface for the whole timestep
+! (when exchange with the atmosphere is disabled): the upper boundary condition becomes a flux condition
+! (and not qsat_surf) and all the remaining surface ice is sublimated and transferred to the subsurface.
+! the actual change in surface ice happens in a higher subroutine
+! =========================================================================================================  =
+            if ( (.not.exchange(ig)) &
+            .and. ( (-zdqsdifrego(ig) * ptimestep) &
+            .gt.( pqsurf(ig) + pdqsdifpot(ig) * ptimestep) ) &
+            .and.( (pqsurf(ig) + pdqsdifpot(ig) * ptimestep).gt.0. ) ) then
+                  ! calculate a new flux from the subsurface
+                  zdqsdifrego(ig) = -( pqsurf(ig) + pdqsdifpot(ig) * ptimestep ) / ptimestep
+!                  ! check case where there is CO2 ice on the surface but qsurf = 0
+!                  if (co2ice(ig).gt.0.) then
+!                        zdqsdifrego(ig) = 0.D0
+!                  endif
+                  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)
+                  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
+                        ! initialize the ice
+                        ice(ig, ik) = iceprev(ig, ik)
+                        if (iceprev(ig, ik).eq.0.) then
+                              ice_index_flag(ik) = .false.
+                        else
+                              ice_index_flag(ik) = .true.
+                        endif
+                  enddo
+                  ! loop while there is new sublimation
+                  sublimation_flag = .true.
+                  sublim_n = 0
+                  do while (sublimation_flag)
+!                        print *, "special case sublimation loop: ", sublim_n
+                        sublim_n = sublim_n + 1
+                        if (sublim_n.gt.100) then
+                              print *, "special case sublimation not converging in call ", n
+                              print *, "might as well stop now"
+                              stop
+                        endif
+                        ! reset the sublimation flag
+                        sublimation_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
+                              zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1)
+                              zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1)
+                              zbeta(ig, 1) = ( porosity_prev(ig, 1) * C(ig, 1) * znsoilprev2(ig, 1) + B(ig, 1) - zdqsdifrego(ig) ) &
+                                    / zdelta(ig, 1)
+                              ! check for ice
+                              if (ice_index_flag(1)) then
+                                          ! set the alpha coefficient to zero
+                                          zalpha(ig, 1) = 0.D0
+                                          ! set the beta coefficient to the number density at saturation pressure
+                                          zbeta(ig, 1) = nsatsoil(ig, 1)
+                              endif
+                              do ik = 2, nsoil - 1  ! loop over all other layers
+                                    ! calculate the Delta, Alpha, and Beta coefficients in the layer
+                                    zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) + D_inter(ig, ik) / midlayer_dz(ig, ik) &
+                                          + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1))
+                                    zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik)
+                                    zbeta(ig, ik) = ( D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) &
+                                          + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev2(ig, ik) + B(ig, ik) ) / zdelta(ig, ik)
+                                    ! check for ice
+                                    if (ice_index_flag(ik)) then
+                                          ! set the alpha coefficient to zero
+                                          zalpha(ig, ik) = 0.D0
+                                          ! set the beta coefficient to the number density at saturation pressure
+                                          zbeta(ig, ik) = nsatsoil(ig, ik)
+                                    endif
+                              enddo
+                              ! calculate the preliminary amount of water vapor in the bottom layer
+                              if (ice_index_flag(nsoil)) then
+                                    ! set the vapor number density to saturation
+                                    znsoil(ig, nsoil) = nsatsoil(ig, nsoil)
+                              else
+                                    ! calculate the vapor number density in the lowest layer
+                                    znsoil(ig, nsoil) = ( D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) &
+                                          + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev2(ig, nsoil) + B(ig, nsoil) ) &
+                                          / ( porosity(ig, nsoil) * C(ig, nsoil) &
+                                          + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)) &
+                                          + A(ig, nsoil) )
+                              endif
+                             ! 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))
+                                  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
+                                    ! 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
+                                                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
+                                                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 ! 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
+                        enddo  ! end while loop (adsorption saturation)
+                        ! set the flux to the previously calculated value for the top layer
+                        flux(ig, 0) = zdqsdifrego(ig)
+                        ! calculate the flux for all other layers
+                        do ik = 1, nsoil - 1
+                              flux(ig, ik) = D_inter(ig, ik) * (znsoil(ig, ik + 1) - znsoil(ig, ik)) / midlayer_dz(ig, ik)
+                        enddo
+                        ! calculate the change in ztot1
+                        do ik = 1, nsoil - 1
+                              dztot1(ik) = (flux(ig, ik) - flux(ig, ik - 1)) / interlayer_dz(ig, ik) &
+                                    - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) &
+                                    + B(ig, ik) / interlayer_dz(ig, ik)
+                        enddo
+                        dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) &
+                    enddo
+                enddo  ! end while loop (adsorption saturation)
+                ! Set the flux to the previously calculated value for the top layer
+                flux(ig, 0) = zdqsdifrego(ig)
+                ! Calculate the flux for all other layers
+                do ik = 1, nsoil - 1
+                    flux(ig, ik) = D_inter(ig, ik) * (znsoil(ig, ik + 1) - znsoil(ig, ik)) / midlayer_dz(ig, ik)
+                enddo
+                ! Calculate the change in ztot1
+                do ik = 1, nsoil - 1
+                    dztot1(ik) = (flux(ig, ik) - flux(ig, ik - 1)) / interlayer_dz(ig, ik) &
+                               - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) &
+                               + B(ig, ik) / interlayer_dz(ig, ik)
+                enddo
+                dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) &
                               - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) &
                               + B(ig, nsoil) / interlayer_dz(ig, nsoil)
+                        ! Condensation / sublimation
+                        do ik = 1, nsoil
+                              if (ice_index_flag(ik)) then
+                                    ! Compute ice content
+                                    ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep &
+                                          - porosity(ig, ik) * nsatsoil(ig, ik)
+                                    if (ice(ig, ik).lt.0.D0) then  ! if all the ice is used up
+!                                          print *, "############complete sublimation in layer ", ik
+!                                          print *, "ztot1: ", ztot1(ig, ik)
+!                                          print *, "dztot1*timestep: ", dztot1(ik) * ptimestep
+!                                          print *, "vapour: ", porosity(ig, ik) * nsatsoil(ig, ik)
+!                                          print *, "znsoil: ", znsoil(ig, ik)
+!                                          print *, "nsatsoil: ", nsatsoil(ig, ik)
+!                                          print *, "porosity: ", porosity(ig, ik)
+!                                          print *, "ice: ", ice(ig, ik)
+!                                          print *, "exchange: ", exchange(ig)
+!                                          print *, "co2ice: ", co2ice(ig)
+!                                          print *, ""
+!                                          print *, "zalpha: ", zalpha(ig, ik)
+!                                          print *, "zbeta: ", zbeta(ig, ik)
+!                                          print *, ""
+                                          ! set the ice content to zero
+                                          ice(ig, ik) = 0.D0
+                                          if (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then
+                                                print *, "WARNING! complete sublimation, but vapor is oversaturated"
+                                                print *, "special case loop in cell", ig, ik ,"will be supressed"
+                                          else
+                                                ! change the water vapor from the previous timestep. (Watch out! could go wrong)
+                                                znsoilprev2(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik)
+                                                ! print *, "new znsoilprev", znsoilprev(ig, ik)
+                                                ! remove the ice flag and raise the sublimation flag
+                                                ice_index_flag(ik) = .false.
+                                                sublimation_flag = .true.
+                                          endif
+                                    endif
+                              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
+                                    ice_index_flag(ik) = .true.
+                                    sublimation_flag = .true.
+!                                    print *, condensation_flag(ik)
+                                    condensation_flag(ik) = .true.
+!                                    print *, condensation_flag(ik)
+!                                    print *, "set condensation flag"
+!                              else
+!                                    print *, "special case condensation loop supressed"
+                              endif
+                              endif
+                        enddo
+                  enddo  ! end first loop while (condensation / sublimation)
+            endif  ! Special case for all ice on the surface is used up
+            do ik = 1, nsoil
+                  diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) &  ! only used for output
+                        + adswater(ig, ik) - adswprev(ig, ik) &
+                        - ztot1(ig, ik) - dztot1(ik) * ptimestep  !? This is inconsistent, in the not special case the previous adsorbed water is not counted
+                                                                  !? also this just overwrites the previous calculation if I see that correctly
+                  ! calculate the total amount of water
+                  ztot1(ig, ik) = porosity(ig, ik) * znsoil(ig, ik) + ice(ig, ik)
+                  h2otot(ig, ik) = adswater(ig, ik) + ztot1(ig, ik)
+            enddo
+      endif  ! if there is no polar cap
+                ! Condensation / sublimation
+                do ik = 1, nsoil
+                    if (ice_index_flag(ik)) then
+                        ! Compute ice content
+                        ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik)
+                        if (ice(ig, ik) < 0.D0) then  ! if all the ice is used up
+                            ! Set the ice content to zero
+                            ice(ig, ik) = 0.D0
+                            if (znsoil(ig, ik) > nsatsoil(ig, ik)) then
+                                print *, "WARNING! complete sublimation, but vapor is oversaturated"
+                                print *, "special case loop in cell", ig, ik, "will be suppressed"
+                            else
+                                ! Change the water vapor from the previous timestep. (Watch out! could go wrong)
+                                znsoilprev2(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik)
+                                ! Remove the ice flag and raise the sublimation flag
+                                ice_index_flag(ik) = .false.
+                                sublimation_flag = .true.
+                            endif
+                        endif
+                    elseif (znsoil(ig, ik) > nsatsoil(ig, ik)) then  ! if there is new ice through condensation
+                        if (.not. condensation_flag(ik)) then
+                            ! Raise the ice and sublimation flags
+                            ice_index_flag(ik) = .true.
+                            sublimation_flag = .true.
+                            condensation_flag(ik) = .true.
+                        endif
+                    endif
+                enddo
+            enddo  ! end first loop while (condensation / sublimation)
+        endif  ! Special case for all ice on the surface is used up
+        do ik = 1, nsoil
+            diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) &  ! only used for output
+                         + adswater(ig, ik) - adswprev(ig, ik) &
+                         - ztot1(ig, ik) - dztot1(ik) * ptimestep
+            ! Calculate the total amount of water
+            ztot1(ig, ik) = porosity(ig, ik) * znsoil(ig, ik) + ice(ig, ik)
+            h2otot(ig, ik) = adswater(ig, ik) + ztot1(ig, ik)
+        enddo
+    endif  ! if there is no polar cap
 enddo  ! loop over all gridpoints
 ! check for choking condition
+! Check for choking condition
 do ig = 1, ngrid
+      if(.not.watercaptag(ig)) then
+            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).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).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
+                              endif
+                              close_top(ig, ik) = .false.
+                        endif
+                        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).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
+                              endif
+                              close_bottom(ig, ik) = .false.
+                        endif
+!                        if ( (flux(ig, ik).lt.flux(ig, ik - 1)).and.(flux(ig, ik).gt.0D0) ) then
+!                              if (.not.close_top(ig, ik).and.print_closure_warnings) then
+!                                    print *, "NOTICE: closing top of soil layer due to ice", ig, ik
+!                              endif
+!                              close_top(ig, ik) = .true.
+!
+!                        elseif ( (flux(ig, ik).lt.flux(ig, ik - 1)).and.(flux(ig, ik - 1).lt.0D0) ) then
+!                              if (.not.close_bottom(ig, ik).and.print_closure_warnings) then
+!                                    print *, "NOTICE: closing bottom of soil layer due to ice", ig, ik
+!                              endif
+!                              close_bottom(ig, ik) = .true.
+!                        endif
+!
+!                        if (close_top(ig, ik).and.close_bottom(ig, ik)) then
+!                              close_top(ig, ik) = .false.
+!                              close_bottom(ig, ik) = .false.
+!                              if (print_closure_warnings) then
+!                                    print *, "WARNING: Reopening soil layer after complete closure:", ig, ik
+!                              endif
+!                        elseif (close_top(ig, ik).or.close_bottom(ig, ik)) then
+!                              if ((flux(ig, ik) - flux(ig, ik - 1)).gt.0D0) then
+!                                    close_top(ig, ik) = .false.
+!                                    close_bottom(ig, ik) = .false.
+!                                    if (print_closure_warnings) then
+!                                          print *, "WARNING: Reopening soil layer due to rising ice:", ig, ik
+!                                    endif
+!                              endif
+!
+!                              if (close_top(ig, ik).and.(flux(ig, ik - 1).lt.0D0)) then
+!                                    close_top = .false.
+!                                    if (print_closure_warnings) then
+!                                          print *, "NOTICE: Reopening top of soil layer due to upward flux:", ig, ik
+!                                    endif
+!                              endif
+!
+!                              if (close_bottom(ig, ik).and.(flux(ig, ik).gt.0D0)) then
+!                                    close_bottom = .false.
+!                                    if (print_closure_warnings) then
+!                                          print *, "NOTICE: Reopening bottom of soil layer due to downward flux:", ig, ik
+!                                    endif
+!                              endif
+!                        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
+                              print *, "NOTICE: Reopening soillayer due to decrease in ice", ig, ik
+                        endif
+                        close_top(ig, ik) = .false.
+                        close_bottom(ig, ik) = .false.
+                  endif
+            enddo
+      endif
+    if (.not. watercaptag(ig)) then
+        do ik = 1, nsoil - 1
+            if (ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)) > choke_fraction) then  ! if the ice is over saturation or choke_fraction
+                if (flux(ig, ik - 1) > 0.D0) 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) < 0.D0) then
+                    if (close_top(ig, ik) .and. print_closure_warnings) then
+                        print *, "NOTICE: opening top of soil layer due to outward flux", ig, ik
+                    endif
+                    close_top(ig, ik) = .false.
+                endif
+                if (flux(ig, ik) < 0.D0) 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) > 0.D0) then
+                    if (close_bottom(ig, ik) .and. print_closure_warnings) then
+                        print *, "NOTICE: opening bottom of soil layer due to outward flux", ig, ik
+                    endif
+                    close_bottom(ig, ik) = .false.
+                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
+                    print *, "NOTICE: Reopening soil layer due to decrease in ice", ig, ik
+                endif
+                close_top(ig, ik) = .false.
+                close_bottom(ig, ik) = .false.
+            endif
+        enddo
+    endif
 enddo
 ! Computation of total mass
 ! ===========================
+! =========================
 do ig = 1, ngrid
+      ! initialize the masses to zero
+      mass_h2o(ig) = 0.D0
+      mass_ice(ig) = 0.D0
+      if(.not.watercaptag(ig)) then
+            do ik = 1, nsoil
+                  ! calculate the total water and ice mass
+                  mass_h2o(ig) = mass_h2o(ig) +  h2otot(ig, ik) * interlayer_dz(ig, ik)
+                  mass_ice(ig) = mass_ice(ig) +  ice(ig, ik) * interlayer_dz(ig, ik)
+                  ! 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)
+                  ! write the results to the return variable
+                  pqsoil(ig, ik, igcm_h2o_vap_soil) = znsoil(ig, ik)
+                  pqsoil(ig, ik, igcm_h2o_ice_soil) = ice(ig, ik)
+                  pqsoil(ig, ik, igcm_h2o_vap_ads) = adswater(ig, ik)
+                  if (close_top(ig, ik)) then
+                        close_out(ig, ik) = 1
+                  elseif (close_bottom(ig, ik)) then
+                        close_out(ig, ik) = -1
+                  else
+                        close_out(ig, ik) = 0
+                  endif
+            enddo
+      endif
+      if (watercaptag(ig)) then
+            do ik = 1, nsoil
+                  saturation_water_ice(ig, ik) = -1
+            enddo
+      endif
+      ! convert the logical :: flag exchange to a numeric output
+      if (watercaptag(ig)) then
+            exch(ig) = -2.
+      elseif (exchange(ig)) then
+            exch(ig) = 1.
+      else
+            exch(ig) = -1.
+      endif
+    ! Initialize the masses to zero
+    mass_h2o(ig) = 0.D0
+    mass_ice(ig) = 0.D0
+    if (.not. watercaptag(ig)) then
+        do ik = 1, nsoil
+            ! Calculate the total water and ice mass
+            mass_h2o(ig) = mass_h2o(ig) + h2otot(ig, ik) * interlayer_dz(ig, ik)
+            mass_ice(ig) = mass_ice(ig) + ice(ig, ik) * interlayer_dz(ig, ik)
+            ! Calculate how close the water vapor content is to saturizing the adsorbed water
+            if (choice_ads /= 0) preduite(ig, ik) = znsoil(ig, ik) / nsat(ig, ik)
+            ! Write the results to the return variable
+            pqsoil(ig, ik, igcm_h2o_vap_soil) = znsoil(ig, ik)
+            pqsoil(ig, ik, igcm_h2o_ice_soil) = ice(ig, ik)
+            pqsoil(ig, ik, igcm_h2o_vap_ads) = adswater(ig, ik)
+            if (close_top(ig, ik)) then
+                close_out(ig, ik) = 1
+            elseif (close_bottom(ig, ik)) then
+                close_out(ig, ik) = -1
+            else
+                close_out(ig, ik) = 0
+            endif
+        enddo
+    endif
+    if (watercaptag(ig)) then
+        do ik = 1, nsoil
+            saturation_water_ice(ig, ik) = -1
+        enddo
+    endif
+    ! Convert the logical flag exchange to a numeric output
+    if (watercaptag(ig)) then
+        exch(ig) = -2.D0
+    elseif (exchange(ig)) then
+        exch(ig) = 1.D0
+    else
+        exch(ig) = -1.D0
+    endif
 enddo
 ! calculate the global total value
+! Calculate the global total value
 mass_h2o_tot = 0.D0
 mass_ice_tot = 0.D0
 do ig = 1, ngrid
       mass_h2o_tot = mass_h2o_tot + mass_h2o(ig) * cell_area(ig)
       mass_ice_tot = mass_ice_tot + mass_ice(ig) * cell_area(ig)
+    mass_h2o_tot = mass_h2o_tot + mass_h2o(ig) * cell_area(ig)
+    mass_ice_tot = mass_ice_tot + mass_ice(ig) * cell_area(ig)
 enddo
+! print and iterate the call number
+! print * , 'Subsurface call n = ', n
+! Increment the call number
 n = n + 1
 …
 ! -----------------------------------------------------------------------
 ! reformat flux, because it has a unusual shape
+! Reformat flux, because it has an unusual shape
 do ig = 1, ngrid
+      nsurf(ig) = rhosurf(ig) * pq(ig, 1, igcm_h2o_vap)
+      do ik = 1, nsoil - 1
+              var_flux_soil(ig, ik) = flux(ig, ik)
+      enddo
+      var_flux_soil(ig, nsoil) = 0.
+    nsurf(ig) = rhosurf(ig) * pq(ig, 1, igcm_h2o_vap)
+    do ik = 1, nsoil - 1
+        var_flux_soil(ig, ik) = flux(ig, ik)
+    enddo
+    var_flux_soil(ig, nsoil) = 0.D0
 enddo
+if(writeoutput) then
+!print *, "flux shape: ", shape(flux)
+!print *, "adswater shape ", shape(adswater)
+!print *, "ngrid= ", ngrid
+!if (ngrid.eq.1) then  ! if in 1D
+!      print *, "writing 1D data"
+!      print *, real(adswater(:, :), 4)
+!      print *, shape(adswater(:, :))
+      zalpha(1, nsoil) = 0
+      zbeta(1, nsoil) = 0
+!       call write_output("zalpha","diffusion alpha coefficient", "unknown",real(zalpha(1, :)))
+!       call write_output("zbeta","diffusion beta coefficient", "unknown",real(zbeta(1, :)))
+!       call write_output("adswater","adswater", "kg / m^3",real(adswater(1, :)))
+!       call write_output("znsoil","znsoil", "kg / m^3",real(znsoil(1, :)))
+!       call write_output("ice","ice", "kg / m^3",real(ice(1, :)))
+!       call write_output("h2otot","total water content", "kg / m^3",real(h2otot(1, :)))
+!       call write_output("flux_soil","flux_soil", "kg / m^2",real(flux(1, :)))
+!       call write_output("flux","flux soil", "kg / m^2",real(zdqsdifrego(:)))
+!       call write_output("flux_surf","surface flux", "kg / m^2",real(zdqsdifrego(1)))
+!       call write_output("exchange","exchange", "boolean",real(exch(1)))
+!else
+!      print *, mass_h2o_tot
+!      print *, real(mass_h2o_tot, 4)
+!       call write_output("dztot1","change in dztot", "unknown",real(dztot1))
+        call write_output("flux_soillayer","flux of water between the soil layers", "kg m-2 s-1",var_flux_soil)
+!       call write_output("A_coef","A coefficient of dztot1", "unknown",real(B))
+!       call write_output("B_coef","B coefficient of dztot1", "unknown",real(B))
+!       call write_output("H2O_init","initialized H2O", "kg/m^2",real(H2O))
+!      call write_output("H2O_depth_init","initialized H2O depth ", "m",real(H2O_depth))
+       call write_output("ice_saturation_soil","Water ice saturation in the soil layers", "Percent",saturation_water_ice)
+!       call write_output("mass_h2o_tot","total mass of subsurface water over the planet", "kg",real(mass_h2o_tot))
+!       call write_output("mass_ice_tot","total mass of subsurface ice over the planet", "kg",real(mass_ice_tot))
+        call write_output("mass_h2o_soil","Mass of subsurface water column at each point", "kg m-2",(mass_h2o))
+        call write_output("mass_ice_soil","Mass of subsurface ice at each point", "kg m-2",(mass_ice))
+        call write_output("znsoil","Water vapor soil concentration ", "kg.m - 3 of pore air",znsoil)
+       call write_output("nsatsoil","subsurface water vapor saturation density", "kg/m^3",real(nsatsoil))
+       call write_output("nsurf","surface water vapor density", "kg/m^3",real(nsurf))
+       call write_output("adswater","subsurface adsorbed water", "kg/m^3",real(adswater))
+!       call write_output("subsurfice","subsurface ice", "kg/m^3",real(ice))
+        call write_output("flux_rego",'flux of water from the regolith','kg/m^2',zdqsdifrego)
+!       call write_output("exchange",'exchange','no unit',real(exch))
+!       call write_output("close",'close top, bottom, or none (1, -1, 0)','no unit',real(close_out))
+!       call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_inter)
+!endif
+if (writeoutput) then
+    zalpha(1, nsoil) = 0.D0
+    zbeta(1, nsoil) = 0.D0
+    call write_output("flux_soillayer", "flux of water between the soil layers", "kg m-2 s-1", var_flux_soil)
+    call write_output("ice_saturation_soil", "Water ice saturation in the soil layers", "Percent", saturation_water_ice)
+    call write_output("mass_h2o_soil", "Mass of subsurface water column at each point", "kg m-2", mass_h2o)
+    call write_output("mass_ice_soil", "Mass of subsurface ice at each point", "kg m-2", mass_ice)
+    call write_output("znsoil", "Water vapor soil concentration", "kg.m - 3 of pore air", znsoil)
+    call write_output("nsatsoil", "subsurface water vapor saturation density", "kg/m^3", real(nsatsoil))
+    call write_output("nsurf", "surface water vapor density", "kg/m^3", real(nsurf))
+    call write_output("adswater", "subsurface adsorbed water", "kg/m^3", real(adswater))
+    call write_output("flux_rego", 'flux of water from the regolith', 'kg/m^2', zdqsdifrego)
 endif
 END

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 3905 for trunk/LMDZ.MARS/libf

Legend:

trunk/LMDZ.MARS/libf/phymars/soilwater.F90

Download in other formats: