| 1 | subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep, & |
|---|
| 2 | exchange, qsat_surf, pq, pa, pb, pc, pd, pdqsdifpot, pqsurf, & |
|---|
| 3 | pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2, saturation_water_ice) |
|---|
| 4 | |
|---|
| 5 | use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, & |
|---|
| 6 | porosity_reg, ads_const_D, ads_massive_ice |
|---|
| 7 | use comcstfi_h, only: pi, r |
|---|
| 8 | use tracer_mod, only: igcm_h2o_vap |
|---|
| 9 | use surfdat_h, only: watercaptag |
|---|
| 10 | use geometry_mod, only: cell_area, latitude_deg |
|---|
| 11 | #ifndef MESOSCALE |
|---|
| 12 | use write_output_mod, only: write_output |
|---|
| 13 | #endif |
|---|
| 14 | |
|---|
| 15 | implicit none |
|---|
| 16 | |
|---|
| 17 | ! ================================================================================================= |
|---|
| 18 | ! Description: |
|---|
| 19 | ! |
|---|
| 20 | ! This subroutine calculates the profile of water vapor, adsorbed water and ice in the regolith, |
|---|
| 21 | ! down to depth(nsoil) ~ 18m. It calculates the flux of water vapor (zdqsdifrego) between the |
|---|
| 22 | ! subsurface and the atmosphere. |
|---|
| 23 | ! |
|---|
| 24 | ! Water vapor and adsorbed water are treated as two separate subsurface tracers that are connected |
|---|
| 25 | ! by adsorption / desorption coefficients, including an adsorption / desorption kinetic factor. |
|---|
| 26 | ! |
|---|
| 27 | ! The water adsorption isotherm is linear, with a saturation plateau corresponding to one monolayer |
|---|
| 28 | ! of adsorbed water molecules (i.e. an approximation of a Langmuir-type isotherm). The slope of the |
|---|
| 29 | ! isotherm is defined by the enthalpy of adsorption (in kJ/mol). |
|---|
| 30 | ! New since 2020: the adsorption isotherm is now linear by part, to better simulate Type II or Type III |
|---|
| 31 | ! isotherms (and not only Type I) |
|---|
| 32 | ! |
|---|
| 33 | ! The linearized adsorption-diffusion equation is solved first, then the water vapor pressure is |
|---|
| 34 | ! compared to the saturation vapor pressure, and if the latter is reached, ice is formed, the water vapor |
|---|
| 35 | ! pressure is fixed to its saturation value and the adsorption-diffusion equation is solved again. |
|---|
| 36 | ! If ice was already present, its sublimation or growth is calculated. |
|---|
| 37 | ! |
|---|
| 38 | ! This subroutine is called by vdifc.F if the parameters "adsorption_soil" and "water" are set to true in |
|---|
| 39 | ! callphys.def. It returns the flux of water vapor that is injected into or out of the regolith, |
|---|
| 40 | ! and which serves as a boundary condition to calculate the vertical atmospheric profile of water vapor |
|---|
| 41 | ! in vdifc. It also calculates the exchange between the subsurface and surface ice. |
|---|
| 42 | ! |
|---|
| 43 | ! It requires three subsurface tracers to be defined in traceur.def: |
|---|
| 44 | ! - h2o_vap_soil (subsurface water vapor) |
|---|
| 45 | ! - h2o_ice_soil (subsurface ice) |
|---|
| 46 | ! - h2o_ads_vap (adsorbed water) |
|---|
| 47 | ! |
|---|
| 48 | ! Authors: Pierre-Yves Meslin (pmeslin@irap.omp.eu), cleaned by Lucas Lange |
|---|
| 49 | ! For previous version (with numerous commented lines), see -r 3904 |
|---|
| 50 | ! ================================================================================================= |
|---|
| 51 | |
|---|
| 52 | ! Arguments |
|---|
| 53 | ! ====================================================================== |
|---|
| 54 | |
|---|
| 55 | ! Inputs |
|---|
| 56 | integer, intent(in) :: ngrid ! number of points in the model (all lat and long point in one 1D array) |
|---|
| 57 | integer, intent(in) :: nlayer ! number of layers |
|---|
| 58 | integer, intent(in) :: nq ! number of tracers |
|---|
| 59 | integer, intent(in) :: nsoil |
|---|
| 60 | integer, intent(in) :: nqsoil |
|---|
| 61 | logical, save :: firstcall_soil = .true. ! triggers the initialization |
|---|
| 62 | real, intent(in) :: ptsoil(ngrid, nsoil) ! subsurface temperatures |
|---|
| 63 | real, intent(in) :: ptsrf(ngrid) ! surface temperature |
|---|
| 64 | real, intent(in) :: ptimestep ! length of the timestep (unit depends on run.def file) |
|---|
| 65 | logical, intent(in) :: exchange(ngrid) ! logical array describing whether there is exchange with the atmosphere at the current timestep |
|---|
| 66 | |
|---|
| 67 | real, intent(in) :: qsat_surf(ngrid) ! saturation mixing ratio at the surface at the current surface temperature (formerly qsat) |
|---|
| 68 | real, intent(in) :: pq(ngrid, nlayer, nq) ! tracer atmospheric mixing ratio |
|---|
| 69 | real, intent(in) :: pa(ngrid, nlayer) ! coefficients za |
|---|
| 70 | real, intent(in) :: pb(ngrid, nlayer) ! coefficients zb |
|---|
| 71 | real, intent(in) :: pc(ngrid, nlayer) ! coefficients zc |
|---|
| 72 | real, intent(in) :: pd(ngrid, nlayer) ! coefficients zd |
|---|
| 73 | real, intent(in) :: pdqsdifpot(ngrid) ! potential pdqsdif (without subsurface-atmosphere exchange) |
|---|
| 74 | |
|---|
| 75 | real, intent(in) :: pplev(ngrid, nlayer+1) ! atmospheric pressure levels |
|---|
| 76 | real, intent(in) :: rhoatmo(ngrid) ! atmospheric air density (1st layer) (not used right now) |
|---|
| 77 | logical, intent(in) :: writeoutput ! check we are at the last subtimestep for output |
|---|
| 78 | |
|---|
| 79 | ! Variables modified |
|---|
| 80 | real, intent(inout) :: pqsoil(ngrid, nsoil, nqsoil) ! subsurface tracers (water vapor and ice) |
|---|
| 81 | real, intent(in) :: pqsurf(ngrid) ! water ice on the surface |
|---|
| 82 | ! (in kg.m-3: m-3 of pore air for water vapor and m-3 of regolith for water ice and adsorbed water) |
|---|
| 83 | |
|---|
| 84 | ! Outputs |
|---|
| 85 | real, intent(out) :: zdqsdifrego(ngrid) ! flux from subsurface (positive pointing outwards) |
|---|
| 86 | real, intent(out) :: zq1temp2(ngrid) ! temporary atmospheric mixing ratio after exchange with subsurface (kg/kg) |
|---|
| 87 | real*8, intent(out) :: saturation_water_ice(ngrid, nsoil) ! water pore ice saturation level (formerly Sw) |
|---|
| 88 | |
|---|
| 89 | ! Outputs for the output files |
|---|
| 90 | real*8 :: preduite(ngrid, nsoil) ! how close the water vapor density is to adsorption saturation |
|---|
| 91 | integer :: exch(ngrid) ! translates the logical variable exchange into a -1 or 1 |
|---|
| 92 | real*8 :: mass_h2o(ngrid) ! mass of subsurface water column at each point (kg.m-2) (formerly mh2otot) |
|---|
| 93 | real*8 :: mass_ice(ngrid) ! mass of subsurface ice at each point (kg.m-2) (formerly micetot) |
|---|
| 94 | real*8 :: mass_h2o_tot ! mass of subsurface water over the whole planet |
|---|
| 95 | real*8 :: mass_ice_tot ! mass of subsurface ice over the whole planet |
|---|
| 96 | real*8 :: nsurf(ngrid) ! surface tracer density in kg/m^3 |
|---|
| 97 | real*8 :: close_out(ngrid, nsoil) ! output for close_top and close_bottom |
|---|
| 98 | |
|---|
| 99 | ! Local (saved) variables |
|---|
| 100 | ! ====================================================================== |
|---|
| 101 | |
|---|
| 102 | real*8 :: P_sat_soil(ngrid, nsoil) ! water saturation pressure of subsurface cells (at mid-layers) (formerly Psatsoil) |
|---|
| 103 | real*8 :: nsatsoil(ngrid, nsoil) ! gas number density at saturation pressure |
|---|
| 104 | real*8, allocatable, save :: znsoil(:, :) ! water vapor soil concentration (kg.m-3 of pore air) |
|---|
| 105 | real*8 :: znsoilprev(ngrid, nsoil) ! value of znsoil in the previous timestep |
|---|
| 106 | real*8 :: znsoilprev2(ngrid, nsoil) ! second variable for the value of znsoil in the previous timestep |
|---|
| 107 | real*8 :: zdqsoil(ngrid, nsoil) ! increase of pqsoil if sublimation of ice |
|---|
| 108 | real*8, allocatable, save :: ice(:, :) ! ice content of subsurface cells (at mid-layers) (kg/m3) |
|---|
| 109 | real*8 :: iceprev(ngrid, nsoil) |
|---|
| 110 | logical :: ice_index_flag(nsoil) ! flag for ice presence |
|---|
| 111 | real*8, allocatable, save :: adswater(:, :) ! adsorbed water in subsurface cells (at mid-layers) |
|---|
| 112 | real*8 :: adswater_temp(ngrid, nsoil) ! temporary variable for adsorbed water |
|---|
| 113 | logical, allocatable, save :: over_mono_sat_flag(:, :) ! flag for cells above monolayer saturation |
|---|
| 114 | real*8 :: adswprev(ngrid, nsoil) |
|---|
| 115 | logical :: recompute_cell_ads_flag(nsoil) ! flag to check whether coefficients have changed and need recomputing |
|---|
| 116 | |
|---|
| 117 | real*8, save :: adswater_sat_mono ! adsorption monolayer saturation value (kg.m-3 of regolith) |
|---|
| 118 | real*8 :: delta0(ngrid) ! coefficient delta(0) modified |
|---|
| 119 | real*8 :: alpha0(ngrid) |
|---|
| 120 | real*8 :: beta0(ngrid) |
|---|
| 121 | |
|---|
| 122 | ! Adsorption/Desorption variables and parameters |
|---|
| 123 | real*8 :: Ka(ngrid, nsoil) ! adsorption time constant (s-1) before monolayer saturation (first segment of the bilinear function) |
|---|
| 124 | real*8 :: Kd(ngrid, nsoil) ! desorption time constant (s-1) before monolayer saturation (first segment of the bilinear function) |
|---|
| 125 | real*8 :: k_ads_eq(ngrid, nsoil) ! equilibrium adsorption coefficient (formerly kads) before monolayer saturation |
|---|
| 126 | real*8 :: Ka2(ngrid, nsoil) ! adsorption time constant (s-1) after monolayer saturation (second segment of the bilinear function) |
|---|
| 127 | real*8 :: Kd2(ngrid, nsoil) ! desorption time constant (s-1) after monolayer saturation (second segment of the bilinear function) |
|---|
| 128 | real*8 :: k_ads_eq2(ngrid, nsoil) ! equilibrium adsorption coefficient after monolayer saturation (second segment of the bilinear function) |
|---|
| 129 | real*8 :: c0(ngrid, nsoil) ! intercept of the second line in the bilinear function |
|---|
| 130 | real*8, parameter :: kinetic_factor = 0.01 ! (inverse of) characteristic time to reach adsorption equilibrium (s-1): |
|---|
| 131 | ! fixed arbitrarily when kinetics factors are unknown |
|---|
| 132 | |
|---|
| 133 | real*8, allocatable, save :: ztot1(:, :) ! total (water vapor + ice) content (kg.m-3 of soil) |
|---|
| 134 | real*8 :: dztot1(nsoil) |
|---|
| 135 | real*8 :: h2otot(ngrid, nsoil) ! total water content (ice + water vapor + adsorbed water) |
|---|
| 136 | real*8 :: flux(ngrid, 0:nsoil-1) ! fluxes at interfaces (kg.m-2.s-1) (positive = upward) |
|---|
| 137 | real*8 :: rho(ngrid) ! air density (first subsurface layer) |
|---|
| 138 | real*8 :: rhosurf(ngrid) ! surface air density |
|---|
| 139 | |
|---|
| 140 | ! Porosity and tortuosity variables |
|---|
| 141 | real*8, allocatable, save :: porosity_ice_free(:, :) ! porosity with no ice present (formerly eps0) |
|---|
| 142 | real*8, allocatable, save :: porosity(:, :) ! porosity (formerly eps) |
|---|
| 143 | real*8 :: porosity_prev(ngrid, nsoil) ! porosity from previous timestep |
|---|
| 144 | real*8, allocatable, save :: tortuosity(:, :) ! tortuosity factor (formerly tortuo) |
|---|
| 145 | |
|---|
| 146 | real*8 :: saturation_water_ice_inter(ngrid, nsoil) ! water pore ice saturation level at the interlayer |
|---|
| 147 | |
|---|
| 148 | ! Diffusion coefficients |
|---|
| 149 | real*8 :: D_mid(ngrid, nsoil) ! midlayer diffusion coefficients |
|---|
| 150 | real*8 :: D_inter(ngrid, nsoil) ! interlayer diffusion coefficients (formerly D) |
|---|
| 151 | real*8, allocatable, save :: D0(:, :) ! diffusion coefficient prefactor |
|---|
| 152 | real*8 :: omega(ngrid, nsoil) ! H2O-CO2 collision integral |
|---|
| 153 | real*8 :: vth(ngrid, nsoil) ! H2O thermal speed |
|---|
| 154 | real*8, parameter :: Dk0 = 0.459D0 ! Knudsen diffusion coefficient (for saturation_water_ice = 0) |
|---|
| 155 | real*8 :: Dk(ngrid, nsoil) ! Knudsen diffusion coefficient (divided by porosity/tortuosity factor) |
|---|
| 156 | real*8 :: Dk_inter(ngrid, nsoil) ! Knudsen diffusion coefficient at the interlayer |
|---|
| 157 | real*8 :: Dm(ngrid, nsoil) ! molecular diffusion coefficient |
|---|
| 158 | real*8 :: Dm_inter(ngrid, nsoil) ! molecular diffusion coefficient at the interlayer |
|---|
| 159 | |
|---|
| 160 | real*8, parameter :: choke_fraction = 0.8D0 ! fraction of ice that prevents further diffusion |
|---|
| 161 | logical, allocatable, save :: close_top(:, :) ! closing diffusion at the top of a layer if ice rises over saturation |
|---|
| 162 | logical, allocatable, save :: close_bottom(:, :) ! closing diffusion at the bottom of a layer if ice rises over saturation |
|---|
| 163 | logical, parameter :: print_closure_warnings = .true. |
|---|
| 164 | |
|---|
| 165 | ! Coefficients for the diffusion calculations |
|---|
| 166 | real*8 :: A(ngrid, nsoil) ! coefficient in the diffusion formula |
|---|
| 167 | real*8 :: B(ngrid, nsoil) ! coefficient in the diffusion formula |
|---|
| 168 | real*8 :: C(ngrid, nsoil) ! coefficient in the diffusion formula |
|---|
| 169 | real*8 :: E(ngrid, nsoil) ! coefficient in the diffusion formula (before monolayer saturation) |
|---|
| 170 | real*8 :: F(ngrid, nsoil) ! coefficient in the diffusion formula (before monolayer saturation) |
|---|
| 171 | real*8 :: E2(ngrid, nsoil) ! coefficient in the diffusion formula (after monolayer saturation) |
|---|
| 172 | real*8 :: F2(ngrid, nsoil) ! coefficient in the diffusion formula (after monolayer saturation) |
|---|
| 173 | real*8, allocatable, save :: zalpha(:, :) ! alpha coefficients |
|---|
| 174 | real*8, allocatable, save :: zbeta(:, :) ! beta coefficients |
|---|
| 175 | real*8 :: zdelta(ngrid, nsoil-1) ! delta coefficients |
|---|
| 176 | |
|---|
| 177 | ! Distances between layers |
|---|
| 178 | real*8, allocatable, save :: interlayer_dz(:, :) ! distance between the interlayer points in m (formerly interdz) |
|---|
| 179 | real*8, allocatable, save :: midlayer_dz(:, :) ! distance between the midcell points in m (formerly middz) |
|---|
| 180 | |
|---|
| 181 | real*8 :: nsat(ngrid, nsoil) ! amount of water vapor at (adsorption) monolayer saturation |
|---|
| 182 | |
|---|
| 183 | real*8, allocatable, save :: meshsize(:, :) ! scaling/dimension factor for the pore size |
|---|
| 184 | real*8, allocatable, save :: rho_soil(:, :) ! soil density (bulk - kg/m3) (formerly rhosoil) |
|---|
| 185 | real*8, allocatable, save :: cste_ads(:, :) ! prefactor of adsorption coeff. k |
|---|
| 186 | |
|---|
| 187 | ! General constants |
|---|
| 188 | real*8, parameter :: kb = 1.38065D-23 ! Boltzmann constant |
|---|
| 189 | real*8, parameter :: mw = 2.988D-26 ! water molecular weight |
|---|
| 190 | |
|---|
| 191 | ! Adsorption coefficients |
|---|
| 192 | real*8, parameter :: enthalpy_ads = 35.D3 ! enthalpy of adsorption (J.mol-1) options 21.D3, 35.D3, and 60.D3 |
|---|
| 193 | real*8, parameter :: enthalpy_ads2 = 21.D3 ! enthalpy of adsorption (J.mol-1) for the second segment |
|---|
| 194 | real*8, parameter :: DeltaQ_ads = 21.D3 ! heat of adsorption - enthalpy of liquefaction/sublimation for the first segment (J.mol-1) |
|---|
| 195 | real*8, parameter :: DeltaQ_ads2 = 21.D3 ! heat of adsorption - enthalpy of liquefaction/sublimation for the second segment (J.mol-1) |
|---|
| 196 | real*8, parameter :: tau0 = 1D-14 |
|---|
| 197 | real*8, parameter :: S = 17.D3 ! soil specific surface area (m2.kg-1 of solid) options: 17.D3 and 106.D3 |
|---|
| 198 | real*8, parameter:: Sm = 10.6D-20 ! surface of the water molecule (m2) |
|---|
| 199 | |
|---|
| 200 | ! Reference values for choice_ads = 2 |
|---|
| 201 | real*8, parameter :: Tref = 243.15D0 |
|---|
| 202 | real*8, parameter :: nref = 2.31505631D-6 ! not used anymore (for the time being) |
|---|
| 203 | real*8, parameter :: Kd_ref = 0.0206D0 ! not used for the time being |
|---|
| 204 | real*8, parameter :: Ka_ref = 2.4D-4 ! not used for the time being |
|---|
| 205 | real*8, parameter :: Kref = 0.205D-6 ! value obtained from the fit of all adsorption data from Pommerol (2009) |
|---|
| 206 | real*8, parameter :: Kref2 = 0.108D-7 ! value obtained from the fit of all adsorption data from Pommerol (2009) |
|---|
| 207 | |
|---|
| 208 | logical :: recompute_all_cells_ads_flag ! flag to check whether there is a cell of a column that requires recomputing |
|---|
| 209 | logical :: sublimation_flag |
|---|
| 210 | logical :: condensation_flag(nsoil) |
|---|
| 211 | |
|---|
| 212 | ! Variables and parameters for the H2O map |
|---|
| 213 | integer, parameter :: n_long_H2O = 66 ! number of longitudes of the new H2O map |
|---|
| 214 | integer, parameter :: n_lat_H2O = 50 ! number of latitudes of the new H2O map |
|---|
| 215 | real*8, parameter :: rho_H2O_ice = 920.0D0 ! ice density (formerly rhoice) |
|---|
| 216 | real :: latH2O(n_lat_H2O*n_long_H2O) ! latitude at H2O map points |
|---|
| 217 | real :: longH2O(n_lat_H2O*n_long_H2O) ! longitude at H2O map points |
|---|
| 218 | real :: H2O_depth_data(n_lat_H2O*n_long_H2O) ! depth of the ice layer in g/cm^2 at H2O map points |
|---|
| 219 | real :: H2O_cover_data(n_lat_H2O*n_long_H2O) ! H2O content of the cover layer at H2O map points (not used yet) |
|---|
| 220 | real :: dataH2O(n_lat_H2O*n_long_H2O) ! H2O content of the ice layer at H2O map points |
|---|
| 221 | real :: latH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable |
|---|
| 222 | real :: longH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable |
|---|
| 223 | real :: dataH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable |
|---|
| 224 | real :: H2O_depth_data_temp(n_lat_H2O*n_long_H2O) ! intermediate variable |
|---|
| 225 | real, allocatable, save :: H2O(:) ! H2O map interpolated at GCM grid points (in wt%) |
|---|
| 226 | real, allocatable, save :: H2O_depth(:) ! H2O map ice depth interpolated at GCM in g/cm^2 |
|---|
| 227 | |
|---|
| 228 | ! Variables for the 1D case |
|---|
| 229 | real*8, parameter :: mmr_h2o = 0.6D-4 ! water vapor mass mixing ratio (for initialization only) |
|---|
| 230 | real*8 :: diff(ngrid, nsoil) ! difference between total water content and ice + vapor content |
|---|
| 231 | real :: var_flux_soil(ngrid, nsoil) ! output of the flux between soil layers (kg/m^3/s) |
|---|
| 232 | real :: default_diffcoeff = 4e-4 ! diffusion coefficient by default (no variations with Temperature and pressure (m/s^2) |
|---|
| 233 | real :: tol_massiveice = 50. ! tolerance to account for massive ice (kg/m^3) |
|---|
| 234 | |
|---|
| 235 | ! Loop variables and counters |
|---|
| 236 | integer :: ig, ik, i, j ! loop variables |
|---|
| 237 | logical :: output_trigger ! used to limit amount of written output |
|---|
| 238 | integer, save :: n ! number of runs of this subroutine |
|---|
| 239 | integer :: sublim_n ! number of sublimation loop runs |
|---|
| 240 | integer :: satur_mono_n ! number of monolayer saturation loop runs |
|---|
| 241 | |
|---|
| 242 | ! Allocate arrays if not already done |
|---|
| 243 | if (.not. allocated(znsoil)) then |
|---|
| 244 | allocate(znsoil(ngrid, nsoil)) |
|---|
| 245 | allocate(ice(ngrid, nsoil)) |
|---|
| 246 | allocate(adswater(ngrid, nsoil)) |
|---|
| 247 | allocate(ztot1(ngrid, nsoil)) |
|---|
| 248 | allocate(porosity_ice_free(ngrid, nsoil)) |
|---|
| 249 | allocate(porosity(ngrid, nsoil)) |
|---|
| 250 | allocate(tortuosity(ngrid, nsoil)) |
|---|
| 251 | allocate(D0(ngrid, nsoil)) |
|---|
| 252 | allocate(interlayer_dz(ngrid, nsoil)) |
|---|
| 253 | allocate(midlayer_dz(ngrid, 0:nsoil)) |
|---|
| 254 | allocate(zalpha(ngrid, nsoil)) ! extra entry to match output format |
|---|
| 255 | allocate(zbeta(ngrid, nsoil)) ! extra entry to match output format |
|---|
| 256 | allocate(meshsize(ngrid, nsoil)) |
|---|
| 257 | allocate(rho_soil(ngrid, nsoil)) |
|---|
| 258 | allocate(cste_ads(ngrid, nsoil)) |
|---|
| 259 | allocate(H2O(ngrid)) |
|---|
| 260 | allocate(H2O_depth(ngrid)) |
|---|
| 261 | allocate(close_top(ngrid, nsoil)) |
|---|
| 262 | allocate(close_bottom(ngrid, nsoil)) |
|---|
| 263 | allocate(over_mono_sat_flag(ngrid, nsoil)) |
|---|
| 264 | endif |
|---|
| 265 | |
|---|
| 266 | ! Initialization |
|---|
| 267 | ! ================ |
|---|
| 268 | |
|---|
| 269 | if (firstcall_soil) then |
|---|
| 270 | n = 0 |
|---|
| 271 | close_top = .false. |
|---|
| 272 | close_bottom = .false. |
|---|
| 273 | print *, "adsorption enthalpy first segment: ", enthalpy_ads |
|---|
| 274 | print *, "adsorption enthalpy second segment: ", enthalpy_ads2 |
|---|
| 275 | print *, "adsorption BET constant C first segment: ", DeltaQ_ads |
|---|
| 276 | print *, "adsorption BET constant C second segment: ", DeltaQ_ads2 |
|---|
| 277 | print *, "specific surface area: ", S |
|---|
| 278 | |
|---|
| 279 | do ig = 1, ngrid |
|---|
| 280 | ! Initialize soil parameters |
|---|
| 281 | ! =========================== |
|---|
| 282 | |
|---|
| 283 | ! Initialize the midlayer distances |
|---|
| 284 | midlayer_dz(ig, 0) = mlayer(0) |
|---|
| 285 | |
|---|
| 286 | do ik = 1, nsoil - 1 |
|---|
| 287 | midlayer_dz(ig, ik) = mlayer(ik) - mlayer(ik - 1) |
|---|
| 288 | enddo |
|---|
| 289 | |
|---|
| 290 | ! Initialize the interlayer distances |
|---|
| 291 | interlayer_dz(ig, 1) = layer(1) |
|---|
| 292 | do ik = 2, nsoil |
|---|
| 293 | interlayer_dz(ig, ik) = layer(ik) - layer(ik - 1) |
|---|
| 294 | enddo |
|---|
| 295 | |
|---|
| 296 | ! Choice of porous medium and D0 |
|---|
| 297 | ! =============================== |
|---|
| 298 | do ik = 1, nsoil |
|---|
| 299 | ! These properties are defined here in order to enable custom profiles |
|---|
| 300 | if (ads_massive_ice) then |
|---|
| 301 | if (pqsoil(ig, ik, igcm_h2o_ice_soil) .gt. tol_massiveice) then |
|---|
| 302 | porosity_ice_free(ig, ik) = 0.999999 |
|---|
| 303 | else |
|---|
| 304 | porosity_ice_free(ig, ik) = porosity_reg |
|---|
| 305 | endif |
|---|
| 306 | else |
|---|
| 307 | porosity_ice_free(ig, ik) = porosity_reg |
|---|
| 308 | endif |
|---|
| 309 | tortuosity(ig, ik) = 1.5D0 |
|---|
| 310 | rho_soil(ig, ik) = 1.3D3 ! in kg/m3 of regolith (incl. porosity) |
|---|
| 311 | |
|---|
| 312 | meshsize(ig, ik) = 5.D-6 ! characteristic size 1e-6 = 1 micron |
|---|
| 313 | D0(ig, ik) = porosity_ice_free(ig, ik) / tortuosity(ig, ik) |
|---|
| 314 | enddo |
|---|
| 315 | |
|---|
| 316 | ! Choice of adsorption coefficients |
|---|
| 317 | ! ================================== |
|---|
| 318 | ! Unit = kg/m3; From best fit of all adsorption data from Pommerol et al. (2009) |
|---|
| 319 | adswater_sat_mono = 2.6998D-7 * S * rho_soil(ig, 1) |
|---|
| 320 | |
|---|
| 321 | ! Initialization of ice, water vapor and adsorbed water profiles |
|---|
| 322 | ! =============================================================== |
|---|
| 323 | do ik = 1, nsoil |
|---|
| 324 | znsoil(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_soil) |
|---|
| 325 | ice(ig, ik) = pqsoil(ig, ik, igcm_h2o_ice_soil) |
|---|
| 326 | adswater(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_ads) |
|---|
| 327 | |
|---|
| 328 | saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) |
|---|
| 329 | |
|---|
| 330 | if (saturation_water_ice(ig, ik) .lt. 0.D0) then |
|---|
| 331 | print *, "WARNING!! soil water ice negative at ", ig, ik |
|---|
| 332 | print *, "saturation value: ", saturation_water_ice(ig, ik) |
|---|
| 333 | print *, "setting saturation to 0" |
|---|
| 334 | saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) |
|---|
| 335 | endif |
|---|
| 336 | |
|---|
| 337 | porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) |
|---|
| 338 | |
|---|
| 339 | ! Initialize soil as if all points are below the monolayer saturation level |
|---|
| 340 | if (adswater(ig, ik) .gt. adswater_sat_mono) then |
|---|
| 341 | over_mono_sat_flag(ig, ik) = .true. |
|---|
| 342 | else |
|---|
| 343 | over_mono_sat_flag(ig, ik) = .false. |
|---|
| 344 | endif |
|---|
| 345 | enddo |
|---|
| 346 | enddo |
|---|
| 347 | |
|---|
| 348 | print *, "initializing H2O data" |
|---|
| 349 | |
|---|
| 350 | do ig = 1, ngrid |
|---|
| 351 | output_trigger = .true. |
|---|
| 352 | do ik = 1, nsoil |
|---|
| 353 | saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) |
|---|
| 354 | saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) |
|---|
| 355 | porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) |
|---|
| 356 | enddo |
|---|
| 357 | enddo |
|---|
| 358 | |
|---|
| 359 | print *, "Kinetic factor: ", kinetic_factor |
|---|
| 360 | if (kinetic_factor .eq. 0) then |
|---|
| 361 | print *, "WARNING: adsorption is switched off" |
|---|
| 362 | endif |
|---|
| 363 | firstcall_soil = .false. |
|---|
| 364 | endif ! of "if firstcall_soil" |
|---|
| 365 | |
|---|
| 366 | ! Update properties in case of the sublimation of massive ice |
|---|
| 367 | if (ads_massive_ice) then |
|---|
| 368 | do ig = 1, ngrid |
|---|
| 369 | do ik = 1, nsoil |
|---|
| 370 | if (pqsoil(ig, ik, igcm_h2o_ice_soil) .gt. tol_massiveice) then |
|---|
| 371 | porosity_ice_free(ig, ik) = 0.999999 |
|---|
| 372 | else |
|---|
| 373 | porosity_ice_free(ig, ik) = porosity_reg |
|---|
| 374 | endif |
|---|
| 375 | enddo |
|---|
| 376 | enddo |
|---|
| 377 | endif |
|---|
| 378 | |
|---|
| 379 | ! Computation of new (new step) transport coefficients |
|---|
| 380 | ! =================================================== |
|---|
| 381 | |
|---|
| 382 | do ig = 1, ngrid ! loop over all gridpoints |
|---|
| 383 | if (.not. watercaptag(ig)) then ! if there is no polar cap |
|---|
| 384 | |
|---|
| 385 | do ik = 1, nsoil |
|---|
| 386 | ! Calculate water saturation level (pore volume filled by ice / total pore volume) |
|---|
| 387 | saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D0) |
|---|
| 388 | |
|---|
| 389 | if (saturation_water_ice(ig, ik) < 0.D0) then |
|---|
| 390 | print *, "WARNING!! soil water ice negative at ", ig, ik |
|---|
| 391 | print *, "saturation value: ", saturation_water_ice(ig, ik) |
|---|
| 392 | print *, "setting saturation to 0" |
|---|
| 393 | saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) |
|---|
| 394 | endif |
|---|
| 395 | |
|---|
| 396 | ! Save the old porosity |
|---|
| 397 | porosity_prev(ig, ik) = porosity(ig, ik) |
|---|
| 398 | |
|---|
| 399 | ! Calculate the new porosity |
|---|
| 400 | porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) |
|---|
| 401 | ! Note: saturation_water_ice and porosity are computed from the ice content of the previous timestep |
|---|
| 402 | enddo |
|---|
| 403 | |
|---|
| 404 | ! Calculate the air density in the first subsurface layer |
|---|
| 405 | rho(ig) = pplev(ig, 1) / (r * ptsoil(ig, 1)) |
|---|
| 406 | |
|---|
| 407 | ! Calculate diffusion coefficients at mid- and interlayers (with ice content from the previous timestep) |
|---|
| 408 | ! Dependence on saturation_water_ice taken from Hudson et al. (2009) and Meslin et al. (SSSAJ, 2010) |
|---|
| 409 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 410 | |
|---|
| 411 | ! Calculate H2O mean thermal velocity |
|---|
| 412 | vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, ik) / (pi * 18.D-3)) |
|---|
| 413 | |
|---|
| 414 | ! Calculate the H2O - CO2 collision integral (after Mellon and Jakosky, JGR, 1993) |
|---|
| 415 | omega(ig, ik) = 1.442D0 - 0.673D0 * dlog(2.516D-3 * ptsoil(ig, ik)) & |
|---|
| 416 | + 0.252D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 2.D0 & |
|---|
| 417 | - 0.047D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 3.D0 & |
|---|
| 418 | + 0.003D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 4.D0 |
|---|
| 419 | |
|---|
| 420 | ! Calculate the molecular diffusion coefficient |
|---|
| 421 | Dm(ig, ik) = 4.867D0 * ptsoil(ig, ik) ** (3.D0 / 2.D0) & |
|---|
| 422 | / (pplev(ig, 1) * omega(ig, ik)) * 1.D-4 |
|---|
| 423 | |
|---|
| 424 | ! Calculate the Knudsen diffusion coefficient (divided by D0 = porosity / tortuosity) |
|---|
| 425 | Dk(ig, ik) = Dk0 / D0(ig, ik) * meshsize(ig, ik) * vth(ig, ik) |
|---|
| 426 | |
|---|
| 427 | ! Calculate midlayer diffusion coefficient (with dependence on saturation from Meslin et al., 2010) |
|---|
| 428 | D_mid(ig, ik) = D0(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) ** 2.D0 & |
|---|
| 429 | * 1.D0 / (1.D0 / Dm(ig, ik) + 1.D0 / Dk(ig, ik)) |
|---|
| 430 | |
|---|
| 431 | if (ads_const_D) D_mid(ig, ik) = default_diffcoeff |
|---|
| 432 | |
|---|
| 433 | enddo |
|---|
| 434 | |
|---|
| 435 | ! Calculate interlayer diffusion coefficient as interpolation of the midlayer coefficients |
|---|
| 436 | do ik = 1, nsoil - 1 |
|---|
| 437 | Dm_inter(ig, ik) = ((layer(ik) - mlayer(ik - 1)) * Dm(ig, ik + 1) & |
|---|
| 438 | + (mlayer(ik) - layer(ik)) * Dm(ig, ik)) / (mlayer(ik) - mlayer(ik - 1)) |
|---|
| 439 | |
|---|
| 440 | Dk_inter(ig, ik) = ((layer(ik) - mlayer(ik - 1)) * Dk(ig, ik + 1) & |
|---|
| 441 | + (mlayer(ik) - layer(ik)) * Dk(ig, ik)) / (mlayer(ik) - mlayer(ik - 1)) |
|---|
| 442 | |
|---|
| 443 | if (close_bottom(ig, ik) .or. close_top(ig, ik + 1)) then |
|---|
| 444 | saturation_water_ice_inter(ig, ik) = 0.999D0 |
|---|
| 445 | else |
|---|
| 446 | saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik)) |
|---|
| 447 | endif |
|---|
| 448 | |
|---|
| 449 | saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik)) |
|---|
| 450 | |
|---|
| 451 | D_inter(ig, ik) = D0(ig, ik) * (1.D0 - saturation_water_ice_inter(ig, ik)) ** 2.D0 & |
|---|
| 452 | * 1.D0 / (1.D0 / Dm_inter(ig, ik) + 1.D0 / Dk_inter(ig, ik)) |
|---|
| 453 | |
|---|
| 454 | if (ads_const_D) D_inter(ig, ik) = default_diffcoeff |
|---|
| 455 | |
|---|
| 456 | enddo |
|---|
| 457 | endif |
|---|
| 458 | enddo |
|---|
| 459 | |
|---|
| 460 | ! Computation of new adsorption/desorption constants (Ka, Kd coefficients), and C, E, F coefficients |
|---|
| 461 | ! ================================================================================================== |
|---|
| 462 | |
|---|
| 463 | do ig = 1, ngrid ! loop over all gridpoints |
|---|
| 464 | if (.not. watercaptag(ig)) then ! if there is no polar cap |
|---|
| 465 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 466 | |
|---|
| 467 | if (choice_ads == 1) then ! Constant, uniform diffusion coefficient D0 (prescribed) |
|---|
| 468 | |
|---|
| 469 | ! First segment of the bilinear function (before monolayer saturation) |
|---|
| 470 | ! Calculate the equilibrium adsorption coefficient |
|---|
| 471 | k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & |
|---|
| 472 | / (dexp(-enthalpy_ads / (8.314D0 * ptsoil(ig, ik))) / tau0) |
|---|
| 473 | |
|---|
| 474 | ! Calculate the desorption coefficient |
|---|
| 475 | Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 476 | |
|---|
| 477 | ! Calculate the absorption coefficient |
|---|
| 478 | Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 479 | |
|---|
| 480 | ! Second segment of the bilinear function (after monolayer saturation) - added 2020 |
|---|
| 481 | ! Calculate the equilibrium adsorption coefficient |
|---|
| 482 | k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & |
|---|
| 483 | / (dexp(-enthalpy_ads2 / (8.314D0 * ptsoil(ig, ik))) / tau0) |
|---|
| 484 | |
|---|
| 485 | ! Calculate the desorption coefficient |
|---|
| 486 | Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 487 | |
|---|
| 488 | ! Calculate the absorption coefficient |
|---|
| 489 | Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 490 | |
|---|
| 491 | elseif (choice_ads == 2) then ! Modified 2020 (with DeltaQ instead of Q) |
|---|
| 492 | |
|---|
| 493 | ! First segment of the bilinear function (before monolayer saturation) |
|---|
| 494 | ! Calculate the equilibrium adsorption coefficient |
|---|
| 495 | k_ads_eq(ig, ik) = 8.314D0 * rho_soil(ig, ik) * S * ptsoil(ig, ik) / 18.D-3 * Kref & |
|---|
| 496 | * dexp(DeltaQ_ads / 8.314D0 * (1.D0 / ptsoil(ig, ik) - 1.D0 / Tref)) |
|---|
| 497 | |
|---|
| 498 | ! Calculate the desorption coefficient |
|---|
| 499 | Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 500 | |
|---|
| 501 | ! Calculate the absorption coefficient |
|---|
| 502 | Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 503 | |
|---|
| 504 | ! Second segment of the bilinear function (after monolayer saturation) - added 2020 |
|---|
| 505 | ! Calculate the equilibrium adsorption coefficient |
|---|
| 506 | k_ads_eq2(ig, ik) = 8.314D0 * rho_soil(ig, ik) * S * ptsoil(ig, ik) / 18.D-3 * Kref2 & |
|---|
| 507 | * dexp(DeltaQ_ads2 / 8.314D0 * (1.D0 / ptsoil(ig, ik) - 1.D0 / Tref)) |
|---|
| 508 | |
|---|
| 509 | ! Calculate the desorption coefficient |
|---|
| 510 | Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 511 | |
|---|
| 512 | ! Calculate the absorption coefficient |
|---|
| 513 | Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 514 | |
|---|
| 515 | elseif (choice_ads == 0) then ! No adsorption/desorption |
|---|
| 516 | Ka(ig, ik) = 0.D0 |
|---|
| 517 | Kd(ig, ik) = 0.D0 |
|---|
| 518 | Ka2(ig, ik) = 0.D0 |
|---|
| 519 | Kd2(ig, ik) = 0.D0 |
|---|
| 520 | endif |
|---|
| 521 | |
|---|
| 522 | ! Calculate the amount of water vapor at monolayer saturation - modified 2020 |
|---|
| 523 | if (choice_ads /= 0) nsat(ig, ik) = adswater_sat_mono * Kd(ig, ik) / Ka(ig, ik) |
|---|
| 524 | |
|---|
| 525 | ! Calculate the intercept of the second line in the bilinear function - added 2020; corrected 2024 |
|---|
| 526 | c0(ig, ik) = (1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono |
|---|
| 527 | |
|---|
| 528 | ! Calculate C, E, and F coefficients for later calculations |
|---|
| 529 | C(ig, ik) = interlayer_dz(ig, ik) / ptimestep |
|---|
| 530 | |
|---|
| 531 | ! First segment of the bilinear function (before monolayer saturation) |
|---|
| 532 | E(ig, ik) = Ka(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik)) |
|---|
| 533 | F(ig, ik) = Kd(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik)) |
|---|
| 534 | |
|---|
| 535 | ! Second segment of the bilinear function (after monolayer saturation) - added 2020 |
|---|
| 536 | E2(ig, ik) = Ka2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik)) |
|---|
| 537 | F2(ig, ik) = Kd2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik)) |
|---|
| 538 | |
|---|
| 539 | ! Save the old water concentration |
|---|
| 540 | znsoilprev(ig, ik) = znsoil(ig, ik) |
|---|
| 541 | znsoilprev2(ig, ik) = znsoil(ig, ik) ! needed in "special case" loop |
|---|
| 542 | adswprev(ig, ik) = adswater(ig, ik) |
|---|
| 543 | iceprev(ig, ik) = ice(ig, ik) |
|---|
| 544 | |
|---|
| 545 | ! Calculate the total ice + water vapor content |
|---|
| 546 | ztot1(ig, ik) = porosity_prev(ig, ik) * znsoil(ig, ik) + ice(ig, ik) |
|---|
| 547 | enddo |
|---|
| 548 | endif |
|---|
| 549 | enddo |
|---|
| 550 | |
|---|
| 551 | ! Computation of delta, alpha and beta coefficients |
|---|
| 552 | ! ================================================= |
|---|
| 553 | |
|---|
| 554 | do ig = 1, ngrid ! loop over all gridpoints |
|---|
| 555 | ! Initialize the flux to zero to avoid undefined values |
|---|
| 556 | zdqsdifrego(ig) = 0.D0 |
|---|
| 557 | do ik = 0, nsoil - 1 |
|---|
| 558 | flux(ig, ik) = 0.D0 |
|---|
| 559 | enddo |
|---|
| 560 | |
|---|
| 561 | if (.not. watercaptag(ig)) then ! if there is no polar cap |
|---|
| 562 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 563 | if (ice(ig, ik) == 0.D0) then |
|---|
| 564 | ice_index_flag(ik) = .false. |
|---|
| 565 | else |
|---|
| 566 | ice_index_flag(ik) = .true. |
|---|
| 567 | endif |
|---|
| 568 | enddo |
|---|
| 569 | |
|---|
| 570 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 571 | |
|---|
| 572 | ! Calculate A and B coefficients - modified 2020 |
|---|
| 573 | if (.not. over_mono_sat_flag(ig, ik)) then ! Assume cell below the monolayer saturation |
|---|
| 574 | ! Calculate A and B coefficients (first segment of the bilinear function) |
|---|
| 575 | A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 576 | B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) |
|---|
| 577 | else ! Assume cell over the monolayer saturation - added 2020 |
|---|
| 578 | ! Calculate A and B coefficients (second segment of the bilinear function) |
|---|
| 579 | A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 580 | B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & |
|---|
| 581 | + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig, ik) * c0(ig, ik)) |
|---|
| 582 | endif |
|---|
| 583 | |
|---|
| 584 | ! Calculate the saturation pressure |
|---|
| 585 | P_sat_soil(ig, ik) = 611.D0 * dexp(22.5D0 * (1.D0 - 273.16D0 / ptsoil(ig, ik))) |
|---|
| 586 | |
|---|
| 587 | ! Calculate the gas number density at saturation pressure |
|---|
| 588 | nsatsoil(ig, ik) = (P_sat_soil(ig, ik) * mw) / (kb * ptsoil(ig, ik)) |
|---|
| 589 | enddo |
|---|
| 590 | |
|---|
| 591 | ! Calculate the alpha, beta, and delta coefficients for the surface layer |
|---|
| 592 | delta0(ig) = pa(ig, 1) + pb(ig, 2) * (1.D0 - pd(ig, 2)) + porosity_ice_free(ig, 1) * pb(ig, 1) |
|---|
| 593 | alpha0(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / delta0(ig) |
|---|
| 594 | |
|---|
| 595 | beta0(ig) = (pb(ig, 2) * pc(ig, 2) + pq(ig, 1, igcm_h2o_vap) * pa(ig, 1) + pqsurf(ig)) / delta0(ig) |
|---|
| 596 | |
|---|
| 597 | ! Set loop flags |
|---|
| 598 | do ik = 1, nsoil |
|---|
| 599 | condensation_flag(ik) = .false. |
|---|
| 600 | enddo |
|---|
| 601 | |
|---|
| 602 | sublimation_flag = .true. |
|---|
| 603 | sublim_n = 0 |
|---|
| 604 | do while (sublimation_flag) ! while there is sublimation |
|---|
| 605 | sublim_n = sublim_n + 1 |
|---|
| 606 | |
|---|
| 607 | if (sublim_n > 100) then |
|---|
| 608 | print *, "sublimation not converging in call ", n |
|---|
| 609 | print *, "might as well stop now" |
|---|
| 610 | stop |
|---|
| 611 | endif |
|---|
| 612 | |
|---|
| 613 | ! Reset flag |
|---|
| 614 | sublimation_flag = .false. |
|---|
| 615 | |
|---|
| 616 | recompute_all_cells_ads_flag = .true. |
|---|
| 617 | satur_mono_n = 0 |
|---|
| 618 | do while (recompute_all_cells_ads_flag) |
|---|
| 619 | satur_mono_n = satur_mono_n + 1 |
|---|
| 620 | |
|---|
| 621 | ! Reset loop flag |
|---|
| 622 | recompute_all_cells_ads_flag = .false. |
|---|
| 623 | |
|---|
| 624 | ! Calculate the surface air density |
|---|
| 625 | rhosurf(ig) = pplev(ig, 1) / (r * ptsrf(ig)) |
|---|
| 626 | |
|---|
| 627 | ! Calculate the coefficients in the uppermost layer |
|---|
| 628 | if (exchange(ig)) then ! if there is surface exchange |
|---|
| 629 | |
|---|
| 630 | ! Calculate the delta, alpha, and beta coefficients |
|---|
| 631 | zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) & |
|---|
| 632 | + D_inter(ig, 1) / midlayer_dz(ig, 1) & |
|---|
| 633 | + porosity_ice_free(ig, 1) * pb(ig, 1) / (rho(ig) * ptimestep) * (1.D0 - alpha0(ig)) |
|---|
| 634 | |
|---|
| 635 | zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1) |
|---|
| 636 | |
|---|
| 637 | zbeta(ig, 1) = (porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * beta0(ig) & |
|---|
| 638 | + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1)) / zdelta(ig, 1) |
|---|
| 639 | else |
|---|
| 640 | ! Calculate the delta, alpha, and beta coefficients without exchange |
|---|
| 641 | zdelta(ig, 1) = A(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1) & |
|---|
| 642 | + D_mid(ig, 1) / midlayer_dz(ig, 0) + porosity(ig, 1) * C(ig, 1) |
|---|
| 643 | |
|---|
| 644 | zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1) |
|---|
| 645 | |
|---|
| 646 | zbeta(ig, 1) = (D_mid(ig, 1) / midlayer_dz(ig, 0) * qsat_surf(ig) * rhosurf(ig) & |
|---|
| 647 | + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1)) / zdelta(ig, 1) |
|---|
| 648 | endif |
|---|
| 649 | |
|---|
| 650 | ! Check for ice |
|---|
| 651 | if (ice_index_flag(1)) then |
|---|
| 652 | ! Set the alpha coefficient to zero |
|---|
| 653 | zalpha(ig, 1) = 0.D0 |
|---|
| 654 | ! Set the beta coefficient to the number density at saturation pressure |
|---|
| 655 | zbeta(ig, 1) = nsatsoil(ig, 1) |
|---|
| 656 | endif |
|---|
| 657 | |
|---|
| 658 | do ik = 2, nsoil - 1 ! go through all other layers |
|---|
| 659 | |
|---|
| 660 | ! Calculate delta, alpha, and beta coefficients |
|---|
| 661 | zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) & |
|---|
| 662 | + D_inter(ig, ik) / midlayer_dz(ig, ik) & |
|---|
| 663 | + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1)) |
|---|
| 664 | |
|---|
| 665 | zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik) |
|---|
| 666 | |
|---|
| 667 | zbeta(ig, ik) = (D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) & |
|---|
| 668 | + B(ig, ik) + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev(ig, ik)) / zdelta(ig, ik) |
|---|
| 669 | |
|---|
| 670 | ! Check for ice |
|---|
| 671 | if (ice_index_flag(ik)) then |
|---|
| 672 | ! Set the alpha coefficient to zero |
|---|
| 673 | zalpha(ig, ik) = 0.D0 |
|---|
| 674 | ! Set the beta coefficient to the number density at saturation pressure |
|---|
| 675 | zbeta(ig, ik) = nsatsoil(ig, ik) |
|---|
| 676 | endif |
|---|
| 677 | enddo |
|---|
| 678 | |
|---|
| 679 | ! Computation of pqsoil, ztot1, zq1temp2 and zdqsdifrego |
|---|
| 680 | ! ======================================================= |
|---|
| 681 | |
|---|
| 682 | ! Calculate the preliminary amount of water vapor in the bottom layer |
|---|
| 683 | if (ice_index_flag(nsoil)) then ! if there is ice |
|---|
| 684 | ! Set the vapor number density to saturation |
|---|
| 685 | znsoil(ig, nsoil) = nsatsoil(ig, nsoil) |
|---|
| 686 | else |
|---|
| 687 | ! Calculate the vapor number density in the lowest layer |
|---|
| 688 | znsoil(ig, nsoil) = (D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) & |
|---|
| 689 | + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev(ig, nsoil) + B(ig, nsoil)) & |
|---|
| 690 | / (porosity(ig, nsoil) * C(ig, nsoil) + A(ig, nsoil) & |
|---|
| 691 | + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1))) |
|---|
| 692 | endif |
|---|
| 693 | |
|---|
| 694 | ! Calculate the preliminary amount of adsorbed water |
|---|
| 695 | if (.not. over_mono_sat_flag(ig, nsoil)) then ! modified 2024 |
|---|
| 696 | adswater_temp(ig, nsoil) = (Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil)) & |
|---|
| 697 | / (1.D0 + ptimestep * Kd(ig, nsoil)) |
|---|
| 698 | else |
|---|
| 699 | adswater_temp(ig, nsoil) = (Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) & |
|---|
| 700 | + ptimestep * c0(ig, nsoil) * Kd2(ig, nsoil)) & |
|---|
| 701 | / (1.D0 + ptimestep * Kd2(ig, nsoil)) |
|---|
| 702 | endif |
|---|
| 703 | |
|---|
| 704 | do ik = nsoil - 1, 1, -1 ! backward loop over all layers |
|---|
| 705 | |
|---|
| 706 | znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik) |
|---|
| 707 | |
|---|
| 708 | if (.not. over_mono_sat_flag(ig, ik)) then ! modified 2024 |
|---|
| 709 | adswater_temp(ig, ik) = (Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik)) & |
|---|
| 710 | / (1.D0 + ptimestep * Kd(ig, ik)) |
|---|
| 711 | else |
|---|
| 712 | adswater_temp(ig, ik) = (Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) & |
|---|
| 713 | + ptimestep * c0(ig, ik) * Kd2(ig, ik)) & |
|---|
| 714 | / (1.D0 + ptimestep * Kd2(ig, ik)) |
|---|
| 715 | endif |
|---|
| 716 | enddo |
|---|
| 717 | |
|---|
| 718 | ! Check if any cell is over monolayer saturation |
|---|
| 719 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 720 | |
|---|
| 721 | ! Change of coefficients |
|---|
| 722 | recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients |
|---|
| 723 | |
|---|
| 724 | if (adswater_temp(ig, ik) >= adswater_sat_mono) then |
|---|
| 725 | |
|---|
| 726 | print *, "NOTICE: over monolayer saturation" |
|---|
| 727 | |
|---|
| 728 | if (.not. over_mono_sat_flag(ig, ik)) then |
|---|
| 729 | ! If the cell enters this scope, it means that the cell is over the monolayer |
|---|
| 730 | ! saturation after calculation, but the coefficients assume it is below. |
|---|
| 731 | ! Therefore, the cell needs to be recomputed |
|---|
| 732 | |
|---|
| 733 | over_mono_sat_flag(ig, ik) = .true. |
|---|
| 734 | recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients |
|---|
| 735 | |
|---|
| 736 | ! Change the A and B coefficients (change from first segment to second segment) |
|---|
| 737 | A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 738 | B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & |
|---|
| 739 | + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig, ik) * c0(ig, ik)) |
|---|
| 740 | endif |
|---|
| 741 | else |
|---|
| 742 | if (over_mono_sat_flag(ig, ik)) then |
|---|
| 743 | ! If the cell enters this scope, it means that the cell is below the monolayer |
|---|
| 744 | ! saturation after calculation, but the coefficients assume it is above. |
|---|
| 745 | ! Therefore, the cell needs to be recomputed |
|---|
| 746 | |
|---|
| 747 | over_mono_sat_flag(ig, ik) = .false. |
|---|
| 748 | recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients |
|---|
| 749 | |
|---|
| 750 | ! Calculate A and B coefficients (reset to first segment of the bilinear function) |
|---|
| 751 | A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 752 | B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) |
|---|
| 753 | endif |
|---|
| 754 | endif |
|---|
| 755 | enddo |
|---|
| 756 | |
|---|
| 757 | ! If one cell needs to be recomputed, then all the column is to be recomputed too |
|---|
| 758 | do ik = 1, nsoil |
|---|
| 759 | if (recompute_cell_ads_flag(ik)) then |
|---|
| 760 | recompute_all_cells_ads_flag = .true. |
|---|
| 761 | else |
|---|
| 762 | adswater(ig, ik) = adswater_temp(ig, ik) ! if no recomputing is needed, store the value |
|---|
| 763 | endif |
|---|
| 764 | enddo |
|---|
| 765 | enddo ! end loop while recompute_all_cells_ads_flag (monolayer saturation) |
|---|
| 766 | |
|---|
| 767 | ! Calculation of Fnet + Fads |
|---|
| 768 | ! =========================== |
|---|
| 769 | |
|---|
| 770 | ! Calculate the flux variable |
|---|
| 771 | |
|---|
| 772 | ! Calculate the flux in the top layer |
|---|
| 773 | if (exchange(ig)) then ! if there is exchange |
|---|
| 774 | ! Calculate the flux taking diffusion to the atmosphere into account |
|---|
| 775 | flux(ig, 0) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep & |
|---|
| 776 | * (znsoil(ig, 1) / rho(ig) - beta0(ig) - alpha0(ig) * znsoil(ig, 1) / rho(ig)) |
|---|
| 777 | elseif (pqsurf(ig) > 0.D0) then |
|---|
| 778 | ! Assume that the surface is covered by water ice |
|---|
| 779 | flux(ig, 0) = D_mid(ig, 1) / midlayer_dz(ig, 0) * (znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig)) |
|---|
| 780 | endif |
|---|
| 781 | |
|---|
| 782 | ! Check if the ground would take up water from the surface but there is none |
|---|
| 783 | if ((.not. exchange(ig)) .and. (pqsurf(ig) == 0.D0) .and. (flux(ig, 0) < 0.D0)) then |
|---|
| 784 | flux(ig, 0) = 0.D0 |
|---|
| 785 | endif |
|---|
| 786 | |
|---|
| 787 | ! Calculate the flux in all other layers |
|---|
| 788 | do ik = 1, nsoil - 1 |
|---|
| 789 | flux(ig, ik) = D_inter(ig, ik) * (znsoil(ig, ik + 1) - znsoil(ig, ik)) / midlayer_dz(ig, ik) |
|---|
| 790 | enddo |
|---|
| 791 | |
|---|
| 792 | ! Calculate dztot1 which describes the change in ztot1 (water vapor and ice) |
|---|
| 793 | do ik = 1, nsoil - 1 |
|---|
| 794 | dztot1(ik) = (flux(ig, ik) - flux(ig, ik - 1)) / interlayer_dz(ig, ik) & |
|---|
| 795 | - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) & |
|---|
| 796 | + B(ig, ik) / interlayer_dz(ig, ik) |
|---|
| 797 | enddo |
|---|
| 798 | |
|---|
| 799 | dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) & |
|---|
| 800 | - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) & |
|---|
| 801 | + B(ig, nsoil) / interlayer_dz(ig, nsoil) |
|---|
| 802 | |
|---|
| 803 | ! Condensation / sublimation |
|---|
| 804 | do ik = 1, nsoil ! loop over all |
|---|
| 805 | if (ice_index_flag(ik)) then ! if there is ice |
|---|
| 806 | |
|---|
| 807 | ! Compute ice content |
|---|
| 808 | ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik) |
|---|
| 809 | |
|---|
| 810 | if (ice(ig, ik) < 0.D0) then ! if all the ice is used up |
|---|
| 811 | |
|---|
| 812 | ! Set the ice content to zero |
|---|
| 813 | ice(ig, ik) = 0.D0 |
|---|
| 814 | |
|---|
| 815 | ! Change the water vapor from the previous timestep. (Watch out! could go wrong) |
|---|
| 816 | znsoilprev(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik) |
|---|
| 817 | |
|---|
| 818 | ! Remove the ice flag and raise the sublimation flag |
|---|
| 819 | ice_index_flag(ik) = .false. |
|---|
| 820 | sublimation_flag = .true. |
|---|
| 821 | endif |
|---|
| 822 | |
|---|
| 823 | elseif (znsoil(ig, ik) > nsatsoil(ig, ik)) then ! if there is condensation |
|---|
| 824 | |
|---|
| 825 | if (.not. condensation_flag(ik)) then |
|---|
| 826 | ! Raise the ice and sublimation flags |
|---|
| 827 | ice_index_flag(ik) = .true. |
|---|
| 828 | sublimation_flag = .true. |
|---|
| 829 | condensation_flag(ik) = .true. |
|---|
| 830 | endif |
|---|
| 831 | endif |
|---|
| 832 | |
|---|
| 833 | ! Calculate the difference between total water content and ice + vapor content (only used for output) |
|---|
| 834 | diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) & |
|---|
| 835 | + adswater(ig, ik) - ztot1(ig, ik) - dztot1(ik) * ptimestep |
|---|
| 836 | enddo ! loop over all layers |
|---|
| 837 | enddo ! end first loop while sublimation_flag (condensation / sublimation) |
|---|
| 838 | |
|---|
| 839 | if (exchange(ig)) then ! if there is exchange |
|---|
| 840 | |
|---|
| 841 | ! Calculate the temporary mixing ratio above the surface |
|---|
| 842 | zq1temp2(ig) = beta0(ig) + alpha0(ig) * znsoil(ig, 1) / rho(ig) |
|---|
| 843 | ! Calculate the flux from the subsurface |
|---|
| 844 | zdqsdifrego(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * (zq1temp2(ig) - znsoil(ig, 1) / rho(ig)) |
|---|
| 845 | |
|---|
| 846 | else |
|---|
| 847 | ! Calculate the flux from the subsurface |
|---|
| 848 | zdqsdifrego(ig) = D_mid(ig, 1) / midlayer_dz(ig, 0) * (znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig)) |
|---|
| 849 | if ((pqsurf(ig) == 0.D0) .and. (zdqsdifrego(ig) < 0.D0)) then |
|---|
| 850 | zdqsdifrego(ig) = 0.D0 |
|---|
| 851 | endif |
|---|
| 852 | endif |
|---|
| 853 | |
|---|
| 854 | ! Special case where there is not enough ice on the surface to supply the subsurface for the whole timestep |
|---|
| 855 | ! (when exchange with the atmosphere is disabled): the upper boundary condition becomes a flux condition |
|---|
| 856 | ! (and not qsat_surf) and all the remaining surface ice is sublimated and transferred to the subsurface. |
|---|
| 857 | ! The actual change in surface ice happens in a higher subroutine |
|---|
| 858 | ! ===================================================================================================== |
|---|
| 859 | |
|---|
| 860 | if ((.not. exchange(ig)) & |
|---|
| 861 | .and. ((-zdqsdifrego(ig) * ptimestep) > (pqsurf(ig) + pdqsdifpot(ig) * ptimestep)) & |
|---|
| 862 | .and. ((pqsurf(ig) + pdqsdifpot(ig) * ptimestep) > 0.D0)) then |
|---|
| 863 | |
|---|
| 864 | ! Calculate a new flux from the subsurface |
|---|
| 865 | zdqsdifrego(ig) = -(pqsurf(ig) + pdqsdifpot(ig) * ptimestep) / ptimestep |
|---|
| 866 | |
|---|
| 867 | do ik = 1, nsoil |
|---|
| 868 | |
|---|
| 869 | ! Calculate A and B coefficients - modified 2020 |
|---|
| 870 | if (.not. over_mono_sat_flag(ig, ik)) then ! Assume cell below the monolayer saturation |
|---|
| 871 | ! Calculate A and B coefficients (first segment of the bilinear function) |
|---|
| 872 | A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 873 | B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) |
|---|
| 874 | else ! Assume cell over the monolayer saturation - added 2020 |
|---|
| 875 | ! Calculate A and B coefficients (second segment of the bilinear function) |
|---|
| 876 | A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 877 | B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & |
|---|
| 878 | + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig, ik) * c0(ig, ik)) |
|---|
| 879 | endif |
|---|
| 880 | |
|---|
| 881 | ! Initialize the ice |
|---|
| 882 | ice(ig, ik) = iceprev(ig, ik) |
|---|
| 883 | if (iceprev(ig, ik) == 0.D0) then |
|---|
| 884 | ice_index_flag(ik) = .false. |
|---|
| 885 | else |
|---|
| 886 | ice_index_flag(ik) = .true. |
|---|
| 887 | endif |
|---|
| 888 | enddo |
|---|
| 889 | |
|---|
| 890 | ! Loop while there is new sublimation |
|---|
| 891 | sublimation_flag = .true. |
|---|
| 892 | sublim_n = 0 |
|---|
| 893 | do while (sublimation_flag) |
|---|
| 894 | sublim_n = sublim_n + 1 |
|---|
| 895 | if (sublim_n > 100) then |
|---|
| 896 | print *, "special case sublimation not converging in call ", n |
|---|
| 897 | print *, "might as well stop now" |
|---|
| 898 | stop |
|---|
| 899 | endif |
|---|
| 900 | |
|---|
| 901 | ! Reset the sublimation flag |
|---|
| 902 | sublimation_flag = .false. |
|---|
| 903 | |
|---|
| 904 | ! Loop until good values accounting for monolayer saturation are found |
|---|
| 905 | recompute_all_cells_ads_flag = .true. |
|---|
| 906 | do while (recompute_all_cells_ads_flag) |
|---|
| 907 | |
|---|
| 908 | ! Reset loop flag |
|---|
| 909 | recompute_all_cells_ads_flag = .false. |
|---|
| 910 | |
|---|
| 911 | ! Calculate the Delta, Alpha, and Beta coefficients in the top layer |
|---|
| 912 | zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1) |
|---|
| 913 | |
|---|
| 914 | zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1) |
|---|
| 915 | |
|---|
| 916 | zbeta(ig, 1) = (porosity_prev(ig, 1) * C(ig, 1) * znsoilprev2(ig, 1) + B(ig, 1) - zdqsdifrego(ig)) / zdelta(ig, 1) |
|---|
| 917 | |
|---|
| 918 | ! Check for ice |
|---|
| 919 | if (ice_index_flag(1)) then |
|---|
| 920 | ! Set the alpha coefficient to zero |
|---|
| 921 | zalpha(ig, 1) = 0.D0 |
|---|
| 922 | ! Set the beta coefficient to the number density at saturation pressure |
|---|
| 923 | zbeta(ig, 1) = nsatsoil(ig, 1) |
|---|
| 924 | endif |
|---|
| 925 | |
|---|
| 926 | do ik = 2, nsoil - 1 ! loop over all other layers |
|---|
| 927 | |
|---|
| 928 | ! Calculate the Delta, Alpha, and Beta coefficients in the layer |
|---|
| 929 | zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) + D_inter(ig, ik) / midlayer_dz(ig, ik) & |
|---|
| 930 | + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1)) |
|---|
| 931 | |
|---|
| 932 | zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik) |
|---|
| 933 | |
|---|
| 934 | zbeta(ig, ik) = (D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) & |
|---|
| 935 | + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev2(ig, ik) + B(ig, ik)) / zdelta(ig, ik) |
|---|
| 936 | |
|---|
| 937 | ! Check for ice |
|---|
| 938 | if (ice_index_flag(ik)) then |
|---|
| 939 | ! Set the alpha coefficient to zero |
|---|
| 940 | zalpha(ig, ik) = 0.D0 |
|---|
| 941 | ! Set the beta coefficient to the number density at saturation pressure |
|---|
| 942 | zbeta(ig, ik) = nsatsoil(ig, ik) |
|---|
| 943 | endif |
|---|
| 944 | enddo |
|---|
| 945 | |
|---|
| 946 | ! Calculate the preliminary amount of water vapor in the bottom layer |
|---|
| 947 | if (ice_index_flag(nsoil)) then |
|---|
| 948 | ! Set the vapor number density to saturation |
|---|
| 949 | znsoil(ig, nsoil) = nsatsoil(ig, nsoil) |
|---|
| 950 | else |
|---|
| 951 | ! Calculate the vapor number density in the lowest layer |
|---|
| 952 | znsoil(ig, nsoil) = (D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) & |
|---|
| 953 | + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev2(ig, nsoil) + B(ig, nsoil)) & |
|---|
| 954 | / (porosity(ig, nsoil) * C(ig, nsoil) & |
|---|
| 955 | + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)) & |
|---|
| 956 | + A(ig, nsoil)) |
|---|
| 957 | endif |
|---|
| 958 | |
|---|
| 959 | ! Calculate the preliminary amount of adsorbed water |
|---|
| 960 | if (.not. over_mono_sat_flag(ig, nsoil)) then ! modified 2024 |
|---|
| 961 | adswater_temp(ig, nsoil) = (Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil)) & |
|---|
| 962 | / (1.D0 + ptimestep * Kd(ig, nsoil)) |
|---|
| 963 | else |
|---|
| 964 | adswater_temp(ig, nsoil) = (Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) & |
|---|
| 965 | + ptimestep * c0(ig, nsoil) * Kd2(ig, nsoil)) & |
|---|
| 966 | / (1.D0 + ptimestep * Kd2(ig, nsoil)) |
|---|
| 967 | endif |
|---|
| 968 | |
|---|
| 969 | do ik = nsoil - 1, 1, -1 ! backward loop over all layers |
|---|
| 970 | |
|---|
| 971 | znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik) |
|---|
| 972 | |
|---|
| 973 | if (.not. over_mono_sat_flag(ig, ik)) then ! modified 2024 |
|---|
| 974 | adswater_temp(ig, ik) = (Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik)) & |
|---|
| 975 | / (1.D0 + ptimestep * Kd(ig, ik)) |
|---|
| 976 | else |
|---|
| 977 | adswater_temp(ig, ik) = (Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) & |
|---|
| 978 | + ptimestep * c0(ig, ik) * Kd2(ig, ik)) & |
|---|
| 979 | / (1.D0 + ptimestep * Kd2(ig, ik)) |
|---|
| 980 | endif |
|---|
| 981 | enddo |
|---|
| 982 | |
|---|
| 983 | ! Check for any saturation |
|---|
| 984 | do ik = 1, nsoil |
|---|
| 985 | ! Change of coefficients |
|---|
| 986 | recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients |
|---|
| 987 | |
|---|
| 988 | if (adswater_temp(ig, ik) > adswater_sat_mono) then |
|---|
| 989 | |
|---|
| 990 | print *, "NOTICE: over monolayer saturation" |
|---|
| 991 | |
|---|
| 992 | if (.not. over_mono_sat_flag(ig, ik)) then |
|---|
| 993 | ! If the cell enters this scope, it means that the cell is over the monolayer |
|---|
| 994 | ! saturation after calculation, but the coefficients assume it is below. |
|---|
| 995 | ! Therefore, the cell needs to be recomputed |
|---|
| 996 | |
|---|
| 997 | over_mono_sat_flag(ig, ik) = .true. |
|---|
| 998 | recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients |
|---|
| 999 | |
|---|
| 1000 | ! Change the A and B coefficients |
|---|
| 1001 | A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1002 | B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & |
|---|
| 1003 | + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig, ik) * c0(ig, ik)) |
|---|
| 1004 | endif |
|---|
| 1005 | else |
|---|
| 1006 | if (over_mono_sat_flag(ig, ik)) then |
|---|
| 1007 | ! If the cell enters this scope, it means that the cell is below the monolayer |
|---|
| 1008 | ! saturation after calculation, but the coefficients assume it is above. |
|---|
| 1009 | ! Therefore, the cell needs to be recomputed |
|---|
| 1010 | |
|---|
| 1011 | over_mono_sat_flag(ig, ik) = .false. |
|---|
| 1012 | recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients |
|---|
| 1013 | |
|---|
| 1014 | ! Calculate A and B coefficients (reset to first segment of the bilinear function) |
|---|
| 1015 | A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1016 | B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) |
|---|
| 1017 | endif |
|---|
| 1018 | endif |
|---|
| 1019 | enddo |
|---|
| 1020 | |
|---|
| 1021 | ! Raise the saturation flag if any layer has saturated |
|---|
| 1022 | do ik = 1, nsoil |
|---|
| 1023 | if (recompute_cell_ads_flag(ik)) then |
|---|
| 1024 | recompute_all_cells_ads_flag = .true. |
|---|
| 1025 | else |
|---|
| 1026 | adswater(ig, ik) = adswater_temp(ig, ik) ! if no recomputing is needed, store the value |
|---|
| 1027 | endif |
|---|
| 1028 | enddo |
|---|
| 1029 | enddo ! end while loop (adsorption saturation) |
|---|
| 1030 | |
|---|
| 1031 | ! Set the flux to the previously calculated value for the top layer |
|---|
| 1032 | flux(ig, 0) = zdqsdifrego(ig) |
|---|
| 1033 | |
|---|
| 1034 | ! Calculate the flux for all other layers |
|---|
| 1035 | do ik = 1, nsoil - 1 |
|---|
| 1036 | flux(ig, ik) = D_inter(ig, ik) * (znsoil(ig, ik + 1) - znsoil(ig, ik)) / midlayer_dz(ig, ik) |
|---|
| 1037 | enddo |
|---|
| 1038 | |
|---|
| 1039 | ! Calculate the change in ztot1 |
|---|
| 1040 | do ik = 1, nsoil - 1 |
|---|
| 1041 | dztot1(ik) = (flux(ig, ik) - flux(ig, ik - 1)) / interlayer_dz(ig, ik) & |
|---|
| 1042 | - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) & |
|---|
| 1043 | + B(ig, ik) / interlayer_dz(ig, ik) |
|---|
| 1044 | enddo |
|---|
| 1045 | |
|---|
| 1046 | dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) & |
|---|
| 1047 | - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) & |
|---|
| 1048 | + B(ig, nsoil) / interlayer_dz(ig, nsoil) |
|---|
| 1049 | |
|---|
| 1050 | ! Condensation / sublimation |
|---|
| 1051 | do ik = 1, nsoil |
|---|
| 1052 | if (ice_index_flag(ik)) then |
|---|
| 1053 | |
|---|
| 1054 | ! Compute ice content |
|---|
| 1055 | ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik) |
|---|
| 1056 | |
|---|
| 1057 | if (ice(ig, ik) < 0.D0) then ! if all the ice is used up |
|---|
| 1058 | |
|---|
| 1059 | ! Set the ice content to zero |
|---|
| 1060 | ice(ig, ik) = 0.D0 |
|---|
| 1061 | |
|---|
| 1062 | if (znsoil(ig, ik) > nsatsoil(ig, ik)) then |
|---|
| 1063 | print *, "WARNING! complete sublimation, but vapor is oversaturated" |
|---|
| 1064 | print *, "special case loop in cell", ig, ik, "will be suppressed" |
|---|
| 1065 | else |
|---|
| 1066 | ! Change the water vapor from the previous timestep. (Watch out! could go wrong) |
|---|
| 1067 | znsoilprev2(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik) |
|---|
| 1068 | |
|---|
| 1069 | ! Remove the ice flag and raise the sublimation flag |
|---|
| 1070 | ice_index_flag(ik) = .false. |
|---|
| 1071 | sublimation_flag = .true. |
|---|
| 1072 | endif |
|---|
| 1073 | endif |
|---|
| 1074 | |
|---|
| 1075 | elseif (znsoil(ig, ik) > nsatsoil(ig, ik)) then ! if there is new ice through condensation |
|---|
| 1076 | if (.not. condensation_flag(ik)) then |
|---|
| 1077 | ! Raise the ice and sublimation flags |
|---|
| 1078 | ice_index_flag(ik) = .true. |
|---|
| 1079 | sublimation_flag = .true. |
|---|
| 1080 | condensation_flag(ik) = .true. |
|---|
| 1081 | endif |
|---|
| 1082 | endif |
|---|
| 1083 | enddo |
|---|
| 1084 | enddo ! end first loop while (condensation / sublimation) |
|---|
| 1085 | endif ! Special case for all ice on the surface is used up |
|---|
| 1086 | |
|---|
| 1087 | do ik = 1, nsoil |
|---|
| 1088 | |
|---|
| 1089 | diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) & ! only used for output |
|---|
| 1090 | + adswater(ig, ik) - adswprev(ig, ik) & |
|---|
| 1091 | - ztot1(ig, ik) - dztot1(ik) * ptimestep |
|---|
| 1092 | |
|---|
| 1093 | ! Calculate the total amount of water |
|---|
| 1094 | ztot1(ig, ik) = porosity(ig, ik) * znsoil(ig, ik) + ice(ig, ik) |
|---|
| 1095 | h2otot(ig, ik) = adswater(ig, ik) + ztot1(ig, ik) |
|---|
| 1096 | enddo |
|---|
| 1097 | endif ! if there is no polar cap |
|---|
| 1098 | enddo ! loop over all gridpoints |
|---|
| 1099 | |
|---|
| 1100 | ! Check for choking condition |
|---|
| 1101 | do ig = 1, ngrid |
|---|
| 1102 | if (.not. watercaptag(ig)) then |
|---|
| 1103 | do ik = 1, nsoil - 1 |
|---|
| 1104 | if (ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)) > choke_fraction) then ! if the ice is over saturation or choke_fraction |
|---|
| 1105 | if (flux(ig, ik - 1) > 0.D0) then |
|---|
| 1106 | if (.not. close_top(ig, ik) .and. print_closure_warnings) then |
|---|
| 1107 | print *, "NOTICE: closing top of soil layer due to inward flux", ig, ik |
|---|
| 1108 | endif |
|---|
| 1109 | close_top(ig, ik) = .true. |
|---|
| 1110 | elseif (flux(ig, ik - 1) < 0.D0) then |
|---|
| 1111 | if (close_top(ig, ik) .and. print_closure_warnings) then |
|---|
| 1112 | print *, "NOTICE: opening top of soil layer due to outward flux", ig, ik |
|---|
| 1113 | endif |
|---|
| 1114 | close_top(ig, ik) = .false. |
|---|
| 1115 | endif |
|---|
| 1116 | |
|---|
| 1117 | if (flux(ig, ik) < 0.D0) then |
|---|
| 1118 | if (.not. close_bottom(ig, ik) .and. print_closure_warnings) then |
|---|
| 1119 | print *, "NOTICE: closing bottom of soil layer due to inward flux", ig, ik |
|---|
| 1120 | endif |
|---|
| 1121 | close_bottom(ig, ik) = .true. |
|---|
| 1122 | elseif (flux(ig, ik) > 0.D0) then |
|---|
| 1123 | if (close_bottom(ig, ik) .and. print_closure_warnings) then |
|---|
| 1124 | print *, "NOTICE: opening bottom of soil layer due to outward flux", ig, ik |
|---|
| 1125 | endif |
|---|
| 1126 | close_bottom(ig, ik) = .false. |
|---|
| 1127 | endif |
|---|
| 1128 | |
|---|
| 1129 | else ! opening condition that ice content has fallen sufficiently |
|---|
| 1130 | if ((close_top(ig, ik) .or. close_bottom(ig, ik)) .and. print_closure_warnings) then |
|---|
| 1131 | print *, "NOTICE: Reopening soil layer due to decrease in ice", ig, ik |
|---|
| 1132 | endif |
|---|
| 1133 | close_top(ig, ik) = .false. |
|---|
| 1134 | close_bottom(ig, ik) = .false. |
|---|
| 1135 | endif |
|---|
| 1136 | enddo |
|---|
| 1137 | endif |
|---|
| 1138 | enddo |
|---|
| 1139 | |
|---|
| 1140 | ! Computation of total mass |
|---|
| 1141 | ! ========================= |
|---|
| 1142 | |
|---|
| 1143 | do ig = 1, ngrid |
|---|
| 1144 | ! Initialize the masses to zero |
|---|
| 1145 | mass_h2o(ig) = 0.D0 |
|---|
| 1146 | mass_ice(ig) = 0.D0 |
|---|
| 1147 | |
|---|
| 1148 | if (.not. watercaptag(ig)) then |
|---|
| 1149 | do ik = 1, nsoil |
|---|
| 1150 | ! Calculate the total water and ice mass |
|---|
| 1151 | mass_h2o(ig) = mass_h2o(ig) + h2otot(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1152 | mass_ice(ig) = mass_ice(ig) + ice(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1153 | |
|---|
| 1154 | ! Calculate how close the water vapor content is to saturizing the adsorbed water |
|---|
| 1155 | if (choice_ads /= 0) preduite(ig, ik) = znsoil(ig, ik) / nsat(ig, ik) |
|---|
| 1156 | |
|---|
| 1157 | ! Write the results to the return variable |
|---|
| 1158 | pqsoil(ig, ik, igcm_h2o_vap_soil) = znsoil(ig, ik) |
|---|
| 1159 | pqsoil(ig, ik, igcm_h2o_ice_soil) = ice(ig, ik) |
|---|
| 1160 | pqsoil(ig, ik, igcm_h2o_vap_ads) = adswater(ig, ik) |
|---|
| 1161 | |
|---|
| 1162 | if (close_top(ig, ik)) then |
|---|
| 1163 | close_out(ig, ik) = 1 |
|---|
| 1164 | elseif (close_bottom(ig, ik)) then |
|---|
| 1165 | close_out(ig, ik) = -1 |
|---|
| 1166 | else |
|---|
| 1167 | close_out(ig, ik) = 0 |
|---|
| 1168 | endif |
|---|
| 1169 | enddo |
|---|
| 1170 | endif |
|---|
| 1171 | |
|---|
| 1172 | if (watercaptag(ig)) then |
|---|
| 1173 | do ik = 1, nsoil |
|---|
| 1174 | saturation_water_ice(ig, ik) = -1 |
|---|
| 1175 | enddo |
|---|
| 1176 | endif |
|---|
| 1177 | |
|---|
| 1178 | ! Convert the logical flag exchange to a numeric output |
|---|
| 1179 | if (watercaptag(ig)) then |
|---|
| 1180 | exch(ig) = -2.D0 |
|---|
| 1181 | elseif (exchange(ig)) then |
|---|
| 1182 | exch(ig) = 1.D0 |
|---|
| 1183 | else |
|---|
| 1184 | exch(ig) = -1.D0 |
|---|
| 1185 | endif |
|---|
| 1186 | enddo |
|---|
| 1187 | |
|---|
| 1188 | ! Calculate the global total value |
|---|
| 1189 | mass_h2o_tot = 0.D0 |
|---|
| 1190 | mass_ice_tot = 0.D0 |
|---|
| 1191 | do ig = 1, ngrid |
|---|
| 1192 | mass_h2o_tot = mass_h2o_tot + mass_h2o(ig) * cell_area(ig) |
|---|
| 1193 | mass_ice_tot = mass_ice_tot + mass_ice(ig) * cell_area(ig) |
|---|
| 1194 | enddo |
|---|
| 1195 | |
|---|
| 1196 | ! Increment the call number |
|---|
| 1197 | n = n + 1 |
|---|
| 1198 | |
|---|
| 1199 | ! ----------------------------------------------------------------------- |
|---|
| 1200 | ! Output in diagfi and diagsoil |
|---|
| 1201 | ! ----------------------------------------------------------------------- |
|---|
| 1202 | |
|---|
| 1203 | ! Reformat flux, because it has an unusual shape |
|---|
| 1204 | do ig = 1, ngrid |
|---|
| 1205 | nsurf(ig) = rhosurf(ig) * pq(ig, 1, igcm_h2o_vap) |
|---|
| 1206 | |
|---|
| 1207 | do ik = 1, nsoil - 1 |
|---|
| 1208 | var_flux_soil(ig, ik) = flux(ig, ik) |
|---|
| 1209 | enddo |
|---|
| 1210 | var_flux_soil(ig, nsoil) = 0.D0 |
|---|
| 1211 | enddo |
|---|
| 1212 | |
|---|
| 1213 | #ifndef MESOSCALE |
|---|
| 1214 | if (writeoutput) then |
|---|
| 1215 | zalpha(1, nsoil) = 0.D0 |
|---|
| 1216 | zbeta(1, nsoil) = 0.D0 |
|---|
| 1217 | |
|---|
| 1218 | call write_output("flux_soillayer", "flux of water between the soil layers", "kg m-2 s-1", var_flux_soil) |
|---|
| 1219 | call write_output("ice_saturation_soil", "Water ice saturation in the soil layers", "Percent", saturation_water_ice) |
|---|
| 1220 | call write_output("mass_h2o_soil", "Mass of subsurface water column at each point", "kg m-2", mass_h2o) |
|---|
| 1221 | call write_output("mass_ice_soil", "Mass of subsurface ice at each point", "kg m-2", mass_ice) |
|---|
| 1222 | call write_output("znsoil", "Water vapor soil concentration", "kg.m - 3 of pore air", znsoil) |
|---|
| 1223 | call write_output("nsatsoil", "subsurface water vapor saturation density", "kg/m^3", real(nsatsoil)) |
|---|
| 1224 | call write_output("nsurf", "surface water vapor density", "kg/m^3", real(nsurf)) |
|---|
| 1225 | call write_output("adswater", "subsurface adsorbed water", "kg/m^3", real(adswater)) |
|---|
| 1226 | call write_output("flux_rego", 'flux of water from the regolith', 'kg/m^2', zdqsdifrego) |
|---|
| 1227 | endif |
|---|
| 1228 | #endif |
|---|
| 1229 | END |
|---|
