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
#ifndef MESOSCALE
    use write_output_mod, only: write_output
#endif

    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(ik) = .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, ik)) 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
                    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
                    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
            
            ! 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
            
            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 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
            
            ! Loop while there is new sublimation
            sublimation_flag = .true.
            sublim_n = 0
            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, ik)) 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
                        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
                    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
                            adswater(ig, ik) = adswater_temp(ig, ik)  ! if no recomputing is needed, store the value
                        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) &
                              - 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) < 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
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)) > 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 /= 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
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)
enddo

! Increment the call number
n = n + 1

! -----------------------------------------------------------------------
!      Output in diagfi and diagsoil
! -----------------------------------------------------------------------

! 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.D0
enddo

#ifndef MESOSCALE
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
#endif
END