Changeset 3223


Ignore:
Timestamp:
Feb 16, 2024, 4:07:39 PM (9 months ago)
Author:
llange
Message:

Mars PCM
Adsorption: commit the last version developed by Pierre Yves Meslin
LL

Location:
trunk/LMDZ.MARS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3219 r3223  
    45004500Small correction for the surface layer scheme: the value of the critical Richardson varys if one used the classic yamada scheme (0.2 in this case) or can be a tunning parameter if one uses the atke scheme
    45014501Also add the possibily to write the tke in the diagfi.nc when calling pbl_parameter
     4502
     4503== 16/02/2024 == LL
     4504Adsorption: commit the last version developed by Pierre Yves Meslin
     4505

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

    r3167 r3223  
    1 subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf,tsoil, ptimestep,  &
     1subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep,  &
    22      exchange, qsat_surf, pq, pa, pb, pc, pd, pdqsdifpot, pqsurf,  &
    3       pqsoil, pplev, rhoatmo,writeoutput,zdqsdifrego, zq1temp2)
     3      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2)
    44
    55
    … …  
    1010      use geometry_mod, only: cell_area, latitude_deg
    1111      use write_output_mod, only: write_output
    12 
    1312implicit none
    1413
    … …  
    2019! suburface and the atmosphere.
    2120!
    22 ! Water vapor and adsorbed water are treated as two separate subsurface tracers that are connected
     21! Water vapor and nsoiladsorbed water are treated as two separate subsurface tracers that are connected
    2322! by adsorption / desorption coefficients, including an adsorption / desorption kinetic factor.
    2423!
    … …  
    2625! of adsorbed water molecules (i.e. an approximation of a Langmuir - type isotherm). The slope of the
    2726! isotherm is defined by the enthalpy of adsorption (in kJ / mol).
     27! 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)
    2828!
    2929! The linearized adsorption - diffusion equation is solved first, then the water vapor pressure is
    … …  
    6868integer, intent(in) :: nqsoil
    6969logical, save :: firstcall_soil = .true.             ! triggers the initialization
    70 real, intent(in) :: tsoil(ngrid, nsoil)         ! Subsurface temperatures
     70real, intent(in) :: ptsoil(ngrid, nsoil)         ! Subsurface temperatures
    7171real, intent(in) :: ptsrf(ngrid)                ! Surface temperature
    7272real, intent(in) :: ptimestep                   ! length of the timestep (unit depends on run.def file)
    7373logical, intent(in) :: exchange(ngrid)          ! logical :: array describing whether there is exchange with the atmosphere at the current timestep
    7474
    75 real, intent(in) :: qsat_surf(ngrid)            ! Saturation mixing ratio at the surface
     75real, intent(in) :: qsat_surf(ngrid)           ! Saturation vapor pressure at the current surface temperature (formerly qsat)
    7676real, intent(in) :: pq(ngrid, nlayer, nq)       ! Tracer atmospheric mixing ratio
    7777real, intent(in) :: pa(ngrid, nlayer)           ! Coefficients za
    … …  
    7979real, intent(in) :: pc(ngrid, nlayer)           ! Coefficients zc
    8080real, intent(in) :: pd(ngrid, nlayer)           ! Coefficients zd
    81 real, intent(in) :: pdqsdifpot(ngrid)       ! Potential pdqsdif (without subsurface - atmosphere exchange)
     81real, intent(in) :: pdqsdifpot(ngrid, nq)       ! Potential pdqsdif (without subsurface - atmosphere exchange)
    8282
    8383real, intent(in) :: pplev(ngrid, nlayer+1)      ! Atmospheric pressure levels
    8484real, intent(in) :: rhoatmo(ngrid)              ! Atmospheric air density (1st layer) (not used right now)
    85 
    86 logical, intent(in) :: writeoutput              ! just to check we are at the last subtimestep and we can write
     85logical, intent(in) :: writeoutput              ! just to check we are at the last subtimestep and we
    8786
    8887! Variables modified :
    8988! ----------------------------------------------------------------------
    90 real, intent(inout) :: pqsoil(ngrid, nsoil, nqsoil) ! Subsurface tracers (water vapor and ice)
    91 real, intent(in) :: pqsurf(ngrid)           ! Water ice on the surface
    92                                                       ! (in kg.m - 3 : m - 3 of pore air for water vapor and m - 3 of regolith for water ice and adsorbed water)
     89real, intent(inout) :: pqsoil(ngrid, nsoil, 3*nq) ! Subsurface tracers (water vapor and ice)
     90real, intent(in) :: pqsurf(ngrid)                 ! Water ice on the surface
     91                                                  ! (in kg.m - 3 : m - 3 of pore air for water vapor and m - 3 of regolith for water ice and adsorbed water)
    9392! Outputs :
    9493! ----------------------------------------------------------------------
    … …  
    120119real*8, allocatable, save :: adswater(:, :)     ! Adsorbed water in subsurface cells (at miD-layers) (...)
    121120real*8 :: adswater_temp(ngrid, nsoil)           ! temprory variable for adsorbed water
    122 logical :: ads_water_saturation_flag_2(nsoil)
     121logical, 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
    123122real*8 :: adswprev(ngrid, nsoil)
    124 logical :: ads_water_saturation_flag_1(nsoil)
    125 real*8, save :: adswater_sat                    ! Adsorption saturation value (kg.m - 3 of regolith)
     123logical :: recompute_cell_ads_flag(nsoil) ! ARNAU
     124! Formerly ads_water_saturation_flag_1 but with a twist. This variable now
     125! checks whether coefficients have changed and need recomputing. If the cell
     126! is seen to be over the monolayer saturation level (i.e. the cell fulfills the
     127! condition adswater_temp(ig, ik) > adswater_sat_mono) but the coefficients
     128! have been computed assuming that the cell is below the monolayer saturation
     129! layer (i.e. the variable over_mono_sat_flag(ig, ik) had been set to .false.),
     130! then the cell needs to have its coefficients recomputed according to the
     131! previous condition: adswater_temp(ig, ik) > adswater_sat_mono. Then,
     132! the variable recompute_cell_ads_flag(ik) becomes .true.. ! ARNAU
     133real*8, save :: adswater_sat_mono                    ! Adsorption monolayer saturation value (kg.m - 3 of regolith) ! ARNAU
    126134real*8 :: delta0(ngrid)                         ! Coefficient delta(0) modified
    127135real*8 :: alpha0(ngrid)
    … …  
    129137
    130138! Adsorbtion/Desorption variables and parameters
    131 real*8 :: Ka(ngrid, nsoil)            ! Adsorption time constant (s - 1)
    132 real*8 :: Kd(ngrid, nsoil)            ! Desorption time constant (s - 1)
    133 real*8 :: k_ads_eq(ngrid, nsoil)      ! Equilibrium adsorption coefficient (unitless) (formerly kads)
    134 real*8, parameter :: kinetic_factor = 1      ! (inverse of) Characteristic time to reach adsorption equilibrium (s - 1):
     139real*8 :: Ka(ngrid, nsoil)            ! Adsorption time constant (s - 1) before monolayer saturation (first segment of the bilinear function)
     140real*8 :: Kd(ngrid, nsoil)            ! Desorption time constant (s - 1) before monolayer saturation (first segment of the bilinear function)
     141real*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)
     142real*8 :: Ka2(ngrid, nsoil)           ! Adsorption time constant (s - 1) after monolayer saturation (second segment of the bilinear function) ! modified 2020
     143real*8 :: Kd2(ngrid, nsoil)           ! Desorption time constant (s - 1) after monolayer saturation (second segment of the bilinear function) ! modified 2020
     144real*8 :: k_ads_eq2(ngrid, nsoil)     ! Equilibrium adsorption coefficient (formerly kads) after monolayer saturation (second segment of the bilinear function); unitless ! modified 2020
     145real*8 :: c0(ngrid, nsoil)            ! Intercept of the second line in the bilinear function ! modified 2020
     146real*8, parameter :: kinetic_factor = 0.01      ! (inverse of) Characteristic time to reach adsorption equilibrium (s - 1):
    135147                                                ! fixed arbitrarily when kinetics factors are unknown
    136148                                                ! possible values: ! = 1 / 1800 s - 1 ! / 1.16D-6 /  !( =  10 earth days)! / 1.D0 /  ! / 1.2D-5 /  !
    137149                                         
    138 real*8, allocatable, save :: ztot1(:, :)  ! Total (water vapor +  ice) content (kg.m - 3 of soil) !? why does this not include adsorbed water?
     150real*8, allocatable, save :: ztot1(:, :)  ! Total (water vapor +  ice) content (kg.m - 3 of soil)
    139151real*8 :: dztot1(nsoil)
    140 real*8 :: h2otot(ngrid, nsoil)      ! Total water content (ice +  water vapor +  adsorbed water) (....)
     152real*8 :: h2otot(ngrid, nsoil)      ! Total water content (ice +  water vapor +  adsorbed water)
    141153real*8 :: flux(ngrid, 0:nsoil - 1)  ! Fluxes at interfaces (kg.m - 2.s - 1) (positive = upward)
    142154real*8 :: rho(ngrid)                ! Air density (first subsurface layer)
    … …  
    171183logical, allocatable, save :: close_top(:, :)         ! closing diffusion at the top of a layer if ice rises over saturation
    172184logical, allocatable, save :: close_bottom(:, :)      ! closing diffusion at the bottom of a layer if ice risies over saturation
    173 logical, parameter :: print_closure_warnings = .false.
     185logical, parameter :: print_closure_warnings = .true.
    174186
    175187! Coefficients for the Diffusion calculations
    … …  
    177189real*8 :: B(ngrid, nsoil)           ! Coefficient in the diffusion formula
    178190real*8 :: C(ngrid, nsoil)           ! Coefficient in the diffusion formula
    179 real*8 :: E(ngrid, nsoil)           ! Coefficient in the diffusion formula
    180 real*8 :: F(ngrid, nsoil)           ! Coefficient in the diffusion formula
     191real*8 :: E(ngrid, nsoil)           ! Coefficient in the diffusion formula (before monolayer saturation) ! added 2020
     192real*8 :: F(ngrid, nsoil)           ! Coefficient in the diffusion formula (before monolayer saturation) ! added 2020
     193real*8 :: E2(ngrid, nsoil)           ! Coefficient in the diffusion formula (after monolayer saturation) ! added 2020
     194real*8 :: F2(ngrid, nsoil)           ! Coefficient in the diffusion formula (after monolayer saturation) ! added 2020
    181195real*8, allocatable, save :: zalpha(:, :) ! Alpha coefficients
    182196real*8, allocatable, save :: zbeta(:, :)  ! Beta coefficients
    … …  
    187201real*8, allocatable, save :: midlayer_dz(:, :)        ! Distance between the midcell points in m (formerly middz)
    188202
    189 real*8 :: nsat(ngrid, nsoil)                    ! amount of water vapor at adsoption saturation
     203real*8 :: nsat(ngrid, nsoil)                    ! amount of water vapor at (adsorption) monolayer saturation ! modified 2020
    190204
    191205real*8, allocatable, save :: meshsize(:, :)     ! scaling/dimension factor for the por size
    … …  
    196210real*8, parameter :: kb = 1.38065D-23     ! Boltzmann constant
    197211real*8, parameter :: mw = 2.988D-26             ! Water molecular weight
     212!real*8, parameter :: rho_H2O_ice = 920.D0    ! Ice density
    198213
    199214! adsorption coefficients
    200 real*8, parameter :: enthalpy_ads = 21.D3 ! Enthalpy of adsorption (J.mole - 1) options 21.D3, 35.D3, and 60.D3
     215real*8, parameter :: enthalpy_ads = 35.D3 ! Enthalpy of adsorption (J.mole - 1) options 21.D3, 35.D3, and 60.D3
     216real*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
     217real*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
     218real*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
    201219real*8, parameter :: tau0 = 1D-14
    202220real*8, parameter :: S = 17.D3            ! Soil specific surface area (m2.kg - 1 of solid) options: 17.D3 and 106.D3
    … …  
    206224! Reference values for choice_ads = 2
    207225real*8, parameter :: Tref = 243.15D0
    208 real*8, parameter :: nref = 2.31505631D-6 ! calculated as 18.D-3 / (8.314D0 * Tref) * 0.26D0
    209 real*8, parameter :: Kd_ref = 0.0206D0
    210 real*8, parameter :: Ka_ref = 2.4D-4
    211 real*8, parameter :: Kref = 1.23D7        ! Volcanic tuff @ 243.15 K (obtained at low P / Psat)
    212 
    213 logical :: saturation_flag
     226real*8, parameter :: nref = 2.31505631D-6 ! calculated as 18.D-3 / (8.314D0 * Tref) * 0.26D0 ! not used anymore (for the time being)
     227real*8, parameter :: Kd_ref = 0.0206D0   ! Not used for the time being (would require specific measurements to be known, but it is rarely measured)
     228real*8, parameter :: Ka_ref = 2.4D-4     ! Not used for the time being
     229! 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
     230! real*8, parameter :: Kref2 = 5.95D-7 ! Value from data analysis Pommerol2009 ! ARNAU
     231real*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)
     232real*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)
     233
     234
     235logical :: recompute_all_cells_ads_flag ! Old saturation_flag but with a behaviour change ! ARNAU
     236! The old saturation_flag was used to check whether the cell was saturated or
     237! not, and therefore to assign the saturation value adswater(ig, ik) =
     238! adswater_sat instead of the usual adswater(ig, ik) = adswater_temp(ig, ik)
     239! The old routine using saturation_flag did not require that the A, B, C,...
     240! adsorption coefficients be computed, but the new soilwater subroutine
     241! does. Therefore, the variable recompute_all_cells_ads_flag checks whether
     242! there is a cell of a column that requires recomputing. If at least one cell
     243! requires recomputing (i.e. recompute_cell_ads_flag(ik) is .true.), then
     244! recompute_all_cells_ads_flag becomes true, and the adsorption coefficients,
     245! as well as the recursive equation (appearing in this code as `backward
     246! loop over all layers`) ! ARNAU
    214247logical :: sublimation_flag
    215248logical :: condensation_flag(nsoil)
    … …  
    243276integer, save :: n                        ! number of runs of this subroutine
    244277integer :: sublim_n                       ! number of sublimation loop runs
    245 integer :: satur_n                        ! number of saturation loop runs
    246 
    247 
     278integer :: satur_mono_n                   ! number of monolayer saturation loop runs ! added 2020
    248279
    249280
    250281if (.not. allocated(znsoil)) then
    251       print*, 'Flag A'
    252282      allocate( znsoil(ngrid, nsoil) )
    253283      allocate( ice(ngrid, nsoil) )
    … …  
    271301      allocate( close_top(ngrid, nsoil) )
    272302      allocate( close_bottom(ngrid, nsoil) )
     303      allocate( over_mono_sat_flag(ngrid, nsoil) )
    273304endif
    274305
    … …  
    276307! mugaz ! Molar mass of the atmosphere (g.mol-1) ~43.49
    277308
    278 ! TODOs
    279 ! right now the subroutine gets called even if there is co2 ice on the surface
    280 ! this should not be the case!! (Note it should be called, but here the assumption is that every cover is water ice)
    281309
    282310
    283311! Initialisation :
    284312! =================
    285 
    286313
    287314if (firstcall_soil) then
    … …  
    289316      close_top = .false.
    290317      close_bottom = .false.
    291       print *, "adsorption enthalpy: ", enthalpy_ads
     318      print *, "adsorption enthalpy first segment: ", enthalpy_ads
     319      print *, "adsorption enthalpy second segment: ", enthalpy_ads2
     320      print *, "adsorption BET constant C first segment: ", DeltaQ_ads
     321      print *, "adsorption BET constant C second segment: ", DeltaQ_ads2
    292322      print *, "specific surface area: ", S
    293323
    … …  
    315345                  do ik = 1, nsoil
    316346                       
    317                         ! These properties are defined here in order to enable custum profiles
     347                        ! These properties are defined here in order to enable custom profiles
    318348                        porosity_ice_free(ig, ik) = 0.45D0
    319349                        tortuosity(ig, ik) = 1.5D0
    320                         rho_soil(ig, ik) = 1.3D3
     350                        rho_soil(ig, ik) = 1.3D3 ! in kg/m3 of regolith (incl. porosity)
    321351
    322352                        meshsize(ig, ik) = 5.D-6  ! Characteristic size 1e - 6 = 1 micron : diameter(mode 5) = 10 microns, diameter(mode 1) = 1 microns
    … …  
    327357
    328358                  ! Choice of adsorption coefficients
    329                   ! ===================================                   
    330                   adswater_sat = 0.8D-2 * S / 13.7D3 * rho_soil(ig, 1)  ! Experimental value for volcanic tuff (Pommerol et al., 2009)
     359                  ! ===================================
     360                    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                   
     361                  ! 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)
    331362                  ! theoretical formula is = rho_soil(ig, 1) * S / Sm * mw
    332363
    333                   ! with SSA = 13.7 m2 / g and plateau = 0.8wt%
    334                   !                       adswater_sat = 6D-2 * rho_soil(ig, 1)            ! Experimental value for JSC Mars - 1 (Pommerol et al., 2009)
     364                  ! Numerical values above are for SSA = 13.7 m2 / g and plateau = 0.8wt%
     365                  !                       adswater_sat_mono = 6D-2 * rho_soil(ig, 1)            ! Experimental value for JSC Mars - 1 (Pommerol et al., 2009)
    335366                  ! with SSA = 106 m2 / g and plateau (net) = 6wt%
    336367
    … …  
    340371                        ! Initialisation of pqsoil(t = 0)
    341372                        ! in 1D simulations the initialization happens here
    342                         if (ngrid.eq.1) then
    343                               znsoil(ig, ik) = mmr_h2o * mugaz * 1.D-3 * pplev(ig, 1) &
    344                               / (8.314D0 * tsoil(ig, nsoil - 4))
    345                               !                                   znsoil(ig, ik) = 0.
    346                               ice(ig, ik) = 0.D0  ! 0.1D0 !414.D0
    347                              
    348                               saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
    349                               porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
    350 
    351                               if (choice_ads.eq.1) then
    352                                     vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * tsoil(ig, nsoil - 4) &
    353                                           / (pi * 18.D-3))
    354                                     k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
    355                                           / (dexp(-enthalpy_ads &
    356                                           / (8.314D0 * tsoil(ig, nsoil - 4))) / tau0)
    357                                     Kd(ig, ik) = kinetic_factor  &
    358                                           / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
    359                                     Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) /  &
    360                                           (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
    361 
    362                               elseif (choice_ads.eq.2) then  ! par rapport \E0 valeur experimentale de reference
    363                                     !                                         Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 *
    364                                     !     &                                         (1. / tsoil(ig, nsoil - 4) - 1. / Tref))
    365                                     !                                         Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref *
    366                                     !     &                                         sqrt(tsoil(ig, nsoil - 4) / Tref) / nref
    367 
    368                                     k_ads_eq(ig, ik) = Kref * dexp(enthalpy_ads / 8.314D0 * (1.D0 / tsoil(ig, ik) - 1.D0 / Tref))
    369                                          
    370                                     Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
    371                                          
    372                                     Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
    373        
    374                               endif
    375 
    376                               if(choice_ads .ne. 0) adswater(ig, ik) = min(Ka(ig, ik) / Kd(ig, ik) * znsoil(ig, ik), adswater_sat)
    377                              
    378                         else  ! in 3D simulations initialisation happens with newstart.F
     373!                       if (ngrid.eq.1) then
     374!                             znsoil(ig, ik) = mmr_h2o * mugaz * 1.D-3 * pplev(ig, 1) &
     375!                              / (8.314D0 * ptsoil(ig, nsoil - 4))
     376!                              !                                   znsoil(ig, ik) = 0.
     377!                              ice(ig, ik) = 0.D0  ! 0.1D0 !414.D0
     378                             
     379!                              saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
     380!                              porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
     381                             
     382!                              if (choice_ads.eq.1) then
     383!                                    vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, nsoil - 4) &
     384!                                          / (pi * 18.D-3))
     385
     386                                    ! first segment of the bilinear function
     387!                                    k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
     388!                                          / (dexp(-enthalpy_ads &
     389!                                          / (8.314D0 * ptsoil(ig, nsoil - 4))) / tau0)
     390!                                    Kd(ig, ik) = kinetic_factor  &
     391!                                          / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
     392!                                    Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) /  &
     393!                                          (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
     394
     395                                    ! second segment of the bilinear function ! added 2020
     396!                                    k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
     397!                                          / (dexp(-enthalpy_ads2 &
     398!                                          / (8.314D0 * ptsoil(ig, nsoil - 4))) / tau0)
     399!                                    Kd2(ig, ik) = kinetic_factor  &
     400!                                          / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik))
     401!                                    Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) /  &
     402!                                          (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik))
     403!
     404!                              elseif (choice_ads.eq.2) then  ! par rapport \E0 valeur experimentale de reference
     405!                                    !                                         Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 *
     406!                                    !     &                                         (1. / ptsoil(ig, nsoil - 4) - 1. / Tref))
     407!                                    !                                         Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref *
     408!                                    !     &                                         sqrt(ptsoil(ig, nsoil - 4) / Tref) / nref!
     409!
     410!                                    ! first segment of the bilinear function 
     411!                                    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
     412!                                         
     413!                                    Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
     414!                                         
     415!                                    Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
     416!
     417!                                    ! second segment of the bilinear function ! added 2020
     418!                                    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
     419!                                         
     420!                                    Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik))
     421!                                         
     422!                                    Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik))
     423!                              endif
     424                               
     425!                              if (k_ads_eq(ig,ik) * znsoil(ig, ik) .gt. adswater_sat_mono) then ! modified 2024, after correction of c0 expression
     426!                                 c0(ig, ik) = ( 1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono
     427!                                 adswater(ig, ik) = c0(ig,ik) + k_ads_eq2(ig, ik) * znsoil(ig, ik)
     428!                              else
     429!                                 adswater(ig, ik) = k_ads_eq(ig, ik) * znsoil(ig, ik)
     430!                              endif
     431                             
     432!                        else  ! in 3D simulations initialisation happens with newstart.F
    379433                              znsoil(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_soil)
    380434                              ice(ig, ik) = pqsoil(ig, ik, igcm_h2o_ice_soil)
    381435                              adswater(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_ads)
    382                         endif
    383 
    384                         if (choice_ads.eq.0) then ! no adsorption
    385 
    386                                Ka(:, :) = 0.
    387                                Kd(:, :) = 0.
    388                                adswater(:,:) = 0.
    389 
    390                         endif
     436!                        endif
    391437
    392438                        saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0)
    … …  
    400446                       
    401447                        porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik))
     448                                               
     449                        ! Initialise soil as if all points where below the monolayer saturation level (added 2020) => now depends on value of adswater (modified in 2024)
     450                        if (adswater(ig,ik).gt.adswater_sat_mono) then
     451                           over_mono_sat_flag(ig, ik) = .true. ! 
     452                        else
     453                            over_mono_sat_flag(ig, ik) = .false.   
     454                        endif     
     455
    402456                  enddo
    403457!            endif
    … …  
    455509      print *, "Kinetic factor: ", kinetic_factor
    456510      if (kinetic_factor.eq.0) then
    457             print *, "WARNING: adsorption is switched of"
     511            print *, "WARNING: adsorption is switched off"
    458512      endif
    459513      firstcall_soil = .false.
    … …  
    488542
    489543            ! calculate the air density in the first subsurface layer
    490             rho(ig) = pplev(ig, 1) / (r * tsoil(ig, 1))
     544            rho(ig) = pplev(ig, 1) / (r * ptsoil(ig, 1))
    491545
    492546            ! calculate diffusion coefficients at mid- and interlayers (with ice content from the previous timestep)
    … …  
    495549
    496550                  ! calculate H2O mean thermal velocity
    497                   vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * tsoil(ig, ik) / (pi * 18.D-3))
     551                  vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, ik) / (pi * 18.D-3))
    498552
    499553                  ! The diffusion coefficient is calculated
    500554
    501555                  ! calculate the H2O - CO2 collision integral (after Mellon and Jakosky, JGR, 1993)
    502                   omega(ig, ik) = 1.442D0 - 0.673D0 * dlog(2.516D-3 * tsoil(ig, ik)) &
    503                         + 0.252D0 * (dlog(2.516D-3 * tsoil(ig, ik))) ** 2.D0 &
    504                         - 0.047D0 * (dlog(2.516D-3 * tsoil(ig, ik))) ** 3.D0 &
    505                         + 0.003D0 * (dlog(2.516D-3 * tsoil(ig, ik))) ** 4.D0
     556                  omega(ig, ik) = 1.442D0 - 0.673D0 * dlog(2.516D-3 * ptsoil(ig, ik)) &
     557                        + 0.252D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 2.D0 &
     558                        - 0.047D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 3.D0 &
     559                        + 0.003D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 4.D0
    506560
    507561                  ! calculate the molecular diffusion coefficient
    508                   Dm(ig, ik) = 4.867D0 * tsoil(ig, ik) ** (3.D0 / 2.D0) &
     562                  Dm(ig, ik) = 4.867D0 * ptsoil(ig, ik) ** (3.D0 / 2.D0) &
    509563                        / (pplev(ig, 1) * omega(ig, ik)) * 1.D-4
    510564
    … …  
    532586!                       + (mlayer(ik) - layer(ik)) * saturation_water_ice(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1))
    533587                 
    534                   ! new diffusion interaction
    535588                  if (close_bottom(ig, ik).or.close_top(ig, ik+1)) then
    536589                        saturation_water_ice_inter(ig, ik) = 0.999D0
    537590                  else
    538                         saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik)) 
     591                        saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))  ! new diffusion interaction
    539592                  endif
     593                 
     594                  saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))  ! new diffusion interaction
    540595                 
    541596                  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))
    … …  
    556611                  if (choice_ads.eq.1) then  ! Constant, uniform diffusion coefficient D0 (prescribed)
    557612
     613                        ! first segment of the bilinear function (before monolayer saturation)
    558614                        ! calculate the equilibrium adsorption coefficient
    559615                        k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
    560                               / (dexp(-enthalpy_ads / (8.314D0 * tsoil(ig, ik))) / tau0)
     616                              / (dexp(-enthalpy_ads / (8.314D0 * ptsoil(ig, ik))) / tau0)
    561617
    562618                        ! calculate the desorption coefficient
    … …  
    566622                        Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
    567623
    568                         !                             Kd(ig, ik) = exp(-enthalpy_ads / (8.314 * tsoil(ig, ik)))
     624
     625                        ! second segment of the bilinear function (after monolayer saturation) ! added 2020
     626                        ! calculate the equilibrium adsorption coefficient
     627                        k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 &
     628                              / (dexp(-enthalpy_ads2 / (8.314D0 * ptsoil(ig, ik))) / tau0) ! added 2020
     629
     630                        ! calculate the desorption coefficient
     631                        Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! added 2020
     632
     633                        ! calculate the absorption coefficient
     634                        Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! added 2020
     635
     636                        !                             Kd(ig, ik) = exp(-enthalpy_ads / (8.314 * ptsoil(ig, ik)))
    569637                        !     &                             / tau0     ! * 1.D-18
    570638                        !                             Ka(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.   ! * 1.D-18
    … …  
    576644!                        endif
    577645
    578                   elseif (choice_ads.eq.2) then
     646                  elseif (choice_ads.eq.2) then ! modified 2020 (with DeltaQ instead of Q)
    579647                        !                                   Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 *
    580                         !     &                                   (1. / tsoil(ig, ik) - 1. / Tref))
    581                         !                                   Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref * sqrt(tsoil(ig, ik) / Tref)
     648                        !     &                                   (1. / ptsoil(ig, ik) - 1. / Tref))
     649                        !                                   Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref * sqrt(ptsoil(ig, ik) / Tref)
    582650                        !     &                                   / nref
    583651
     652                        ! first segment of the bilinear function (before monolayer saturation)
    584653                        ! calculate the equilibrium adsorption coefficient
    585                         k_ads_eq(ig, ik) = Kref * dexp(enthalpy_ads / 8.314D0 *  &
    586                               (1.D0 / tsoil(ig, ik) - 1.D0 / Tref))
     654                        k_ads_eq(ig, ik) = 8.314D0 * ptsoil(ig,ik)/18.D-3 *Kref * dexp(DeltaQ_ads / 8.314D0 *  &
     655                              (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
    587656
    588657                        ! calculate the desorption coefficient
    … …  
    591660                        ! calculate the absorption coefficient
    592661                        Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik))
    593                   else  ! no ads
    594 
    595                         Kd(ig, ik) = 0.
    596 
    597                         Ka(ig, ik) = 0.
     662
     663                        ! second segment of the bilinear function (after monolayer saturation) ! added 2020
     664                        ! calculate the equilibrium adsorption coefficient
     665                        k_ads_eq2(ig, ik) = 8.314D0 * ptsoil(ig,ik)/18.D-3 * Kref2 * dexp(DeltaQ_ads2 / 8.314D0 *  &
     666                              (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
     667
     668                        ! calculate the desorption coefficient
     669                        Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik))
     670
     671                        ! calculate the absorption coefficient
     672                        Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik))
    598673                  endif
    599674                 
    600                   ! calculate the amount of water vapor at adorption saturation
    601                  
    602                   if (choice_ads.ne.0) nsat(ig, ik) = adswater_sat * Kd(ig, ik) / Ka(ig, ik)
     675                  ! calculate the amount of water vapor at monolayer saturation ! modified 2020
     676                  nsat(ig, ik) = adswater_sat_mono * Kd(ig, ik) / Ka(ig, ik) !
     677
     678                  ! calculate the intercept of the second line in the bilinear function ! added 2020; corrected 2024
     679                  c0(ig, ik) = ( 1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono
    603680
    604681                  ! calculate C, E, and F coefficients for later calculations
    605682                  C(ig, ik) = interlayer_dz(ig, ik) / ptimestep
     683
     684                  ! first segment of the bilinear function (before monolayer saturation)
    606685                  E(ig, ik) = Ka(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik))
    607686                  F(ig, ik) = Kd(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik))
    608                  
     687
     688                  ! second segment of the bilinear function (after monolayer saturation) ! added 2020
     689                  E2(ig, ik) = Ka2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik))
     690                  F2(ig, ik) = Kd2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik))
     691
    609692                  ! save the old water concentration
    610693                  znsoilprev(ig, ik) = znsoil(ig, ik)
    … …  
    630713      if(.not.watercaptag(ig)) then  ! if there is no polar cap
    631714            do ik = 1, nsoil  ! loop over all soil levels
    632                  
    633                   ! reset loop variables
    634                   ads_water_saturation_flag_1(ik) = .false.
    635                   ads_water_saturation_flag_2(ik) = .false.
    636715                  if (ice(ig, ik).eq.0.) then
    637716                        ice_index_flag(ik) = .false.
    … …  
    641720            enddo
    642721
    643             ! No (adsorption) saturation assumed
     722
    644723            do ik = 1, nsoil  ! loop over all soil levels
    645724                 
    646                   ! calculate A and B coefficients
    647                   A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik)
    648                   B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik)
     725                  ! calculate A and B coefficients ! modified 2020
     726                  if ( .not.over_mono_sat_flag(ig, ik) ) then ! Assume cell below the monolayer saturation
     727                        ! calculate A and B coefficients (first segment of the bilinear function)
     728                        A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik)
     729                        B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik)
     730                  else ! Assume cell over the monolayer saturation ! added 2020
     731                        ! calculate A and B coefficients (second segment of the bilinear function)
     732                        A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik)
     733                        B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) &
     734                                  + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024
     735                  endif
    649736                 
    650737                  ! calculate the saturation pressure
    651                   P_sat_soil(ig, ik) = 611.D0 * dexp(22.5D0 * (1.D0 - 273.16D0 / tsoil(ig, ik))) ! maybe take a new expression (Hsublim = 51.1 kJ / mol) ?
     738                  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) ?
    652739                   
    653740                  ! calculate the gas number density at saturation pressure
    654                   nsatsoil(ig, ik) = (P_sat_soil(ig, ik) * mw) / (kb * tsoil(ig, ik))
     741                  nsatsoil(ig, ik) = (P_sat_soil(ig, ik) * mw) / (kb * ptsoil(ig, ik))
    655742            enddo
    656743
    657744            ! calculate the alpha, beta, and delta coefficients for the surface layer
    658745            delta0(ig) = pa(ig, 1) + pb(ig, 2) * (1.D0 - pd(ig, 2)) + porosity_ice_free(ig, 1) * pb(ig, 1)
    659            
    660746            alpha0(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / delta0(ig)
    661747           
    … …  
    683769                  sublimation_flag = .false.
    684770                 
    685                   saturation_flag = .true.
    686                   satur_n = 0
    687                   do while (saturation_flag)
    688 !                        print *, satur_n
    689                         satur_n = satur_n + 1
     771                  recompute_all_cells_ads_flag = .true.
     772                  satur_mono_n = 0
     773                  do while (recompute_all_cells_ads_flag)
     774!                        print *, satur_mono_n
     775                        satur_mono_n = satur_mono_n + 1
    690776                       
    691777                        ! reset loop flag
    692                         saturation_flag = .false.
     778                        recompute_all_cells_ads_flag = .false.
    693779                       
    694780                        ! calculate the surface air density
    … …  
    770856                       
    771857                        ! calculate the preliminary amount of adsorbed water
    772                         adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) &
    773                               / (1.D0 + ptimestep * Kd(ig, nsoil))
    774                        
     858                        if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024
     859                            adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) &
     860                                 / (1.D0 + ptimestep * Kd(ig, nsoil))
     861                        else
     862                             adswater_temp(ig, nsoil) = ( Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) + ptimestep * c0(ig,nsoil) * Kd2(ig,nsoil)) &
     863                                 / (1.D0 + ptimestep * Kd2(ig, nsoil))
     864                        endif                             
     865
     866
     867 
    775868                        do ik = nsoil-1, 1, -1  ! backward loop over all layers
    776869                                 
    777870                              znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik)
    778                              
    779                               adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) &
     871                             
     872                              if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024
     873                                adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) &
    780874                                    / (1.D0 + ptimestep * Kd(ig, ik))
     875                              else
     876                                 adswater_temp(ig, ik) = ( Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) + ptimestep * c0(ig,ik) * Kd2(ig,ik)) &
     877                                    / (1.D0 + ptimestep * Kd2(ig, ik))
     878                              endif 
    781879                        enddo
    782880                       
    783                         ! check for any saturation
     881                        ! check if any cell is over monolayer saturation
    784882                        do ik = 1, nsoil  ! loop over all soil levels
    785                               if ( (adswater_temp(ig, ik).ge.adswater_sat) &
    786                               .and.(.not.ads_water_saturation_flag_2(ik))) then  ! if saturation and there was no saturation previously
    787                              
    788                                     print *, "NOTICE: adsorption saturation"
     883
     884                              ! 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
     885                              !     recompute_cell_ads_flag(ik) = .true.
     886                              !     cycle ! Jump loop
     887                              ! endif
     888                             
     889                              ! Change of coefficients ! ARNAU
     890                              recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients
     891                              if ( adswater_temp(ig, ik).ge.adswater_sat_mono ) then
     892
     893                                    print *, "NOTICE: over monolayer saturation"
     894
     895                                    if ( .not.over_mono_sat_flag(ig, ik) ) then
     896                                          ! If the cell enters this scope, it
     897                                          ! means that the cell is over the monolayer
     898                                          ! saturation after calculation, but the
     899                                          ! coefficients assume it is below. Therefore,
     900                                          ! the cell needs to be recomputed
    789901                                   
    790                                     ! set the adsorbed water to saturation level
    791                                     adswater(ig, ik) = adswater_sat
     902                                          over_mono_sat_flag(ig, ik) = .true.
     903                                          recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients
     904
     905                                          ! 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)
     906                                          A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik)
     907                                          B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) &
     908                                                          + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024
     909                                    endif
     910                              else
     911                                    if ( over_mono_sat_flag(ig, ik) ) then
     912                                          ! If the cell enters this scope, it
     913                                          ! means that the cell is below the monolayer
     914                                          ! saturation after calculation, but the
     915                                          ! coefficients assume it is above. Therefore,
     916                                          ! the cell needs to be recomputed
    792917                                   
    793                                     ! raise the layer saturation flags 1
    794                                     ads_water_saturation_flag_1(ik) = .true.
    795                                     ads_water_saturation_flag_2(ik) = .true.
    796                                    
    797                                     ! change the A and B coefficients
    798                                     A(ig, ik) = 0.D0
    799                                     B(ig, ik) = (adswprev(ig, ik) - adswater_sat) * C(ig, ik)
    800                                    
    801                               ! if there has never been saturation keep the temporary value from earlier
    802                               elseif (.not.ads_water_saturation_flag_2(ik)) then
    803                                     adswater(ig, ik) = adswater_temp(ig, ik)
     918                                          over_mono_sat_flag(ig, ik) = .false.
     919                                          recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients
     920
     921                                          ! calculate A and B coefficients (reset to first segment of the bilinear function)
     922                                          A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik)
     923                                          B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik)
     924                                    endif
    804925                              endif
    805926                        enddo
    806927                       
    807                         ! raise the saturation flag if any layer has saturated and reset the first layer saturation flag
     928                        ! if one cell needs to be recomputed, then all the column is to be recomputed too
    808929                        do ik = 1, nsoil
    809                               if (ads_water_saturation_flag_1(ik)) then
    810                                     saturation_flag = .true.
    811                                     ads_water_saturation_flag_1(ik) = .false.
     930                              if ( recompute_cell_ads_flag(ik) ) then
     931                                    recompute_all_cells_ads_flag = .true.
     932                              else
     933                                    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)
    812934                              endif
    813935                        enddo
    814                   enddo  ! end loop while saturation_flag (adsorption saturation)
     936                  enddo  ! end loop while recompute_all_cells_ads_flag (monolayer saturation)
    815937                 
    816938                  !? I'm not sure if this should be calculated here again. I have a feeling that ztot1 is meant
    … …  
    9381060                  ! calculate the temporty mixing ratio above the surface
    9391061                  zq1temp2(ig) = beta0(ig) + alpha0(ig) * znsoil(ig, 1) / rho(ig)
    940 
    9411062                  ! calculate the flux from the subsurface
    9421063                  zdqsdifrego(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * (zq1temp2(ig) - znsoil(ig, 1) / rho(ig) )
    … …  
    9441065            else
    9451066                  ! calculate the flux from the subsurface
    946                   zdqsdifrego(ig) = -D_mid(ig, 1) / midlayer_dz(ig, 0) * (qsat_surf(ig) * rhosurf(ig) - znsoil(ig, 1)) ! faut - il faire intervenir la porosite ?
     1067                  zdqsdifrego(ig) = D_mid(ig, 1) / midlayer_dz(ig, 0) * (znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig)) ! faut - il faire intervenir la porosite ?
    9471068                  if ((pqsurf(ig).eq.0.).and.(zdqsdifrego(ig).lt.0.)) then
    9481069                        zdqsdifrego(ig) = 0.
    … …  
    9571078
    9581079            if ( (.not.exchange(ig)) &
    959             .and. ( -(zdqsdifrego(ig) * ptimestep) &
    960             .gt.( pqsurf(ig)  + pdqsdifpot(ig) * ptimestep) ) &
    961             .and.( (pqsurf(ig) + pdqsdifpot(ig) * ptimestep).gt.0. ) ) then
     1080            .and. ( (-zdqsdifrego(ig) * ptimestep) &
     1081            .gt.( pqsurf(ig) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep) ) &
     1082            .and.( (pqsurf(ig) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep).gt.0. ) ) then
    9621083
    9631084                  ! calculate a new flux from the subsurface
    964                   zdqsdifrego(ig) = -( pqsurf(ig) + pdqsdifpot(ig) * ptimestep ) / ptimestep
     1085                  zdqsdifrego(ig) = -( pqsurf(ig) + pdqsdifpot(ig, igcm_h2o_ice) * ptimestep ) / ptimestep
    9651086                 
    9661087!                  ! check case where there is CO2 ice on the surface but qsurf = 0
    … …  
    9681089!                        zdqsdifrego(ig) = 0.D0 
    9691090!                  endif
    970                  
    971                   ! calculate the inital A and B coefficients
    972                   do ik = 1, nsoil
     1091
     1092
     1093                  ! calculate A and B coefficients ! modified 2020
     1094                  if ( .not.over_mono_sat_flag(ig, ik) ) then ! Assume cell below the monolayer saturation
     1095                        ! calculate A and B coefficients (first segment of the bilinear function)
    9731096                        A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik)
    9741097                        B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik)
    975                   enddo
     1098                  else ! Assume cell over the monolayer saturation ! added 2020
     1099                        ! calculate A and B coefficients (second segment of the bilinear function)
     1100                        A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik)
     1101                        B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) &
     1102                                  + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024
     1103                  endif
    9761104                 
    9771105                  ! initialize all flags for the loop
    9781106                  do ik = 1, nsoil
    979                  
    980                         ! initialize the first and second layer saturation flag
    981                         ads_water_saturation_flag_1(ik) = .false.
    982                         ads_water_saturation_flag_2(ik) = .false.
    9831107                       
    9841108                        ! initialize the ice
    … …  
    10071131                        sublimation_flag = .false.
    10081132                                               
    1009                         ! loop while there is new saturation
    1010                         saturation_flag = .true.
    1011                         do while (saturation_flag)
    1012                              
    1013                               ! reset the saturation flag
    1014                               saturation_flag = .false.
     1133                        ! loop until good values accounting for monolayer saturation are found
     1134                        recompute_all_cells_ads_flag = .true.
     1135                        do while (recompute_all_cells_ads_flag)
     1136                             
     1137                              ! reset loop flag
     1138                              recompute_all_cells_ads_flag = .false.
    10151139                             
    10161140                              ! calculate the Delta, Alpha, and Beta coefficients in the top layer
    … …  
    10631187                              endif
    10641188                             
    1065                               ! calculate the preliminary amount of adsorbed water
    1066                               adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) &
    1067                                     / ( 1.D0 + ptimestep * Kd(ig, nsoil) )                                   
    1068                              
    1069                               do ik = nsoil - 1, 1, - 1  ! loop from the bottom up
    1070                              
    1071                                     znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik)
    1072                                    
     1189
     1190                             ! calculate the preliminary amount of adsorbed water
     1191                            if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024
     1192                                adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) &
     1193                                     / (1.D0 + ptimestep * Kd(ig, nsoil))
     1194                            else
     1195                                 adswater_temp(ig, nsoil) = ( Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) + ptimestep * c0(ig,nsoil) * Kd2(ig,nsoil)) &
     1196                                 / (1.D0 + ptimestep * Kd2(ig, nsoil))
     1197                            endif
     1198
     1199
     1200
     1201                            do ik = nsoil-1, 1, -1  ! backward loop over all layers
     1202   
     1203                                  znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik)
     1204   
     1205                                  if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024
    10731206                                    adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) &
    1074                                           / (1.D0 + ptimestep * Kd(ig, ik))
    1075                               enddo
     1207                                        / (1.D0 + ptimestep * Kd(ig, ik))
     1208                                  else
     1209                                     adswater_temp(ig, ik) = ( Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) + ptimestep * c0(ig,ik) * Kd2(ig,ik)) &
     1210                                        / (1.D0 + ptimestep * Kd2(ig, ik))
     1211                                  endif
     1212                            enddo
     1213
    10761214                             
    10771215                              ! check for any saturation
    10781216                              do ik = 1, nsoil
    1079                                     if (( adswater_temp(ig, ik).ge.adswater_sat ) &
    1080                                     .and.(.not.ads_water_saturation_flag_2(ik))) then  ! if there is saturation in a new layer
     1217                                    ! Change of coefficients ! ARNAU
     1218                                    recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients ! ARNAU
     1219                                    if ( adswater_temp(ig, ik).gt.adswater_sat_mono ) then
     1220                                   
     1221                                          print *, "NOTICE: over monolayer saturation"
     1222
     1223                                          if ( .not.over_mono_sat_flag(ig, ik) ) then
     1224                                               ! If the cell enters this scope, it
     1225                                               ! means that the cell is over the monolayer
     1226                                               ! saturation after calculation, but the
     1227                                               ! coefficients assume it is below. Therefore,
     1228                                               ! the cell needs to be recomputed
    10811229                                         
    1082                                           ! set the amount of adsorbed water to the saturation value
    1083                                           adswater(ig, ik) = adswater_sat
     1230                                                over_mono_sat_flag(ig, ik) = .true.
     1231                                                recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients
     1232
     1233                                                ! 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)
     1234                                                A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik)
     1235                                                B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) &
     1236                                                          + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024
     1237                                          endif
     1238                                    else
     1239                                          if ( over_mono_sat_flag(ig, ik) ) then
     1240                                               ! If the cell enters this scope, it
     1241                                               ! means that the cell is below the monolayer
     1242                                               ! saturation after calculation, but the
     1243                                               ! coefficients assume it is above. Therefore,
     1244                                               ! the cell needs to be recomputed
    10841245                                         
    1085                                           ! raise both saturation flags
    1086                                           ads_water_saturation_flag_1(ik) = .true.                                         
    1087                                           ads_water_saturation_flag_2(ik) = .true.
    1088 
    1089                                           ! set new A and B coefficients
    1090                                           A(ig, ik) = 0.D0
    1091                                           B(ig, ik) = (adswprev(ig, ik) - adswater_sat) * C(ig, ik)
    1092 
    1093                                     elseif (.not.ads_water_saturation_flag_2(ik)) then  ! if the layer has not saturated before
    1094                                          
    1095                                           ! adopt the preliminary value
    1096                                           adswater(ig, ik) = adswater_temp(ig, ik)
     1246                                                over_mono_sat_flag(ig, ik) = .false.
     1247                                                recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients
     1248
     1249                                                ! calculate A and B coefficients (reset to first segment of the bilinear function)
     1250                                                A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik)
     1251                                                B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik)
     1252                                          endif
    10971253                                    endif
    10981254                              enddo
    10991255                             
    11001256                              ! raise the saturation flag if any layer has saturated and reset the first layer saturation flag
    1101                               do ik = 1, nsoil
    1102                                     if (ads_water_saturation_flag_1(ik)) then
    1103                                           saturation_flag = .true.
    1104                                           ads_water_saturation_flag_1(ik) = .false.
     1257                              do ik = 1, nsoil ! modified 2020
     1258                                    if ( recompute_cell_ads_flag(ik) ) then
     1259                                          recompute_all_cells_ads_flag = .true.
     1260                                    else
     1261                                          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)
    11051262                                    endif
    11061263                              enddo
    … …  
    11701327                                    endif
    11711328
    1172                               elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then  ! if there is new ice through condesation
     1329                              elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then  ! if there is new ice through condensation
    11731330                                    if (.not.condensation_flag(ik)) then
    11741331                                    ! raise the ice and sublimation flags
    … …  
    12061363            do ik = 1, nsoil - 1
    12071364                  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
    1208                         if ( flux(ig, ik - 1).lt.0.) then
     1365                        if ( flux(ig, ik - 1).gt.0.) then
    12091366                              if (.not.close_top(ig, ik).and.print_closure_warnings) then
    12101367                                    print *, "NOTICE: closing top of soil layer due to inward flux", ig, ik     
    12111368                              endif
    12121369                              close_top(ig, ik) = .true.
    1213                         elseif (flux(ig, ik - 1).gt.0.) then
     1370                        elseif (flux(ig, ik - 1).lt.0.) then
    12141371                              if (close_top(ig, ik).and.print_closure_warnings) then
    12151372                                    print *, "NOTICE: opening top of soil layer due to outward flux", ig, ik     
    … …  
    12181375                        endif
    12191376                       
    1220                         if ( flux(ig, ik).gt.0.) then
     1377                        if ( flux(ig, ik).lt.0.) then
    12211378                              if (.not.close_bottom(ig, ik).and.print_closure_warnings) then
    12221379                                    print *, "NOTICE: closing bottom of soil layer due to inward flux", ig, ik     
    12231380                              endif
    12241381                              close_bottom(ig, ik) = .true.
    1225                         elseif (flux(ig, ik).lt.0.) then
     1382                        elseif (flux(ig, ik).gt.0.) then
    12261383                              if (close_bottom(ig, ik).and.print_closure_warnings) then
    12271384                                    print *, "NOTICE: opening bottom of soil layer due to outward flux", ig, ik     
    … …  
    12731430!                        endif
    12741431                  else  ! opening condition that ice content has fallen sufficiently
    1275                         if ((close_top(ig, ik).or.close_bottom(ig, ik)).and.print_closure_warnings) then
     1432                        if (close_top(ig, ik).or.close_bottom(ig, ik).and.print_closure_warnings) then
    12761433                              print *, "NOTICE: Reopening soillayer due to decrease in ice", ig, ik                             
    12771434                        endif
    … …  
    12981455                 
    12991456                  ! calculate how close the water vapor content is to saturizing the adsorbed water
    1300                   if (choice_ads.ne.0) preduite(ig, ik) = znsoil(ig, ik) / nsat(ig, ik)
     1457                  preduite(ig, ik) = znsoil(ig, ik) / nsat(ig, ik)
    13011458                 
    13021459                  ! write the results to the return variable
    … …  
    13061463                 
    13071464
    1308                   if (close_top(ig, ik).and.close_bottom(ig, ik)) then
    1309                         close_out(ig, ik) = 3
    1310                   elseif (close_top(ig, ik)) then
    1311                         close_out(ig, ik) = 2
     1465                  if (close_top(ig, ik)) then
     1466                        close_out(ig, ik) = 1
    13121467                  elseif (close_bottom(ig, ik)) then
    1313                         close_out(ig, ik) = 1
     1468                        close_out(ig, ik) = -1
    13141469                  else
    13151470                        close_out(ig, ik) = 0 
    … …  
    13621517enddo
    13631518
    1364 
    13651519if(writeoutput) then
    13661520!print *, "flux shape: ", shape(flux)
    … …  
    14461600!endif
    14471601endif
    1448 
    14491602RETURN
    14501603END
    14511604
     1605
Note: See TracChangeset for help on using the changeset viewer.