| 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 | |
|---|
| 6 | use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg, & |
|---|
| 7 | ads_const_D, ads_massive_ice |
|---|
| 8 | use comcstfi_h |
|---|
| 9 | use tracer_mod |
|---|
| 10 | use surfdat_h, only: watercaptag ! use mis par AP15 essai |
|---|
| 11 | use geometry_mod, only: cell_area, latitude_deg |
|---|
| 12 | use write_output_mod, only: write_output |
|---|
| 13 | implicit none |
|---|
| 14 | |
|---|
| 15 | ! ================================================================================================= = |
|---|
| 16 | ! Description: |
|---|
| 17 | ! |
|---|
| 18 | ! This subroutine calculates the profile of water vapor, adsorbed water and ice in the regolith, |
|---|
| 19 | ! down to depth(nsoil) ~ 18m. It calculates the flux of water vapor (zdqsdifrego) between the |
|---|
| 20 | ! suburface and the atmosphere. |
|---|
| 21 | ! |
|---|
| 22 | ! Water vapor and adsorbed water are treated as two separate subsurface tracers that are connected |
|---|
| 23 | ! by adsorption / desorption coefficients, including an adsorption / desorption kinetic factor. |
|---|
| 24 | ! |
|---|
| 25 | ! The water adsorption isotherm is linear, with a saturation plateau corresponding to one monolayer |
|---|
| 26 | ! of adsorbed water molecules (i.e. an approximation of a Langmuir - type isotherm). The slope of the |
|---|
| 27 | ! isotherm is defined by the enthalpy of adsorption (in kJ / mol). |
|---|
| 28 | ! 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) |
|---|
| 29 | ! |
|---|
| 30 | ! The linearized adsorption - diffusion equation is solved first, then the water vapor pressure is |
|---|
| 31 | ! compared to the saturation vapor pressure, and if the latter is reached, ice is formed, the water vapor |
|---|
| 32 | ! pressure is fixed to its saturation value and the adsorption - diffusion equation is solved again. |
|---|
| 33 | ! If ice was already present, its sublimation or growth is calculated. |
|---|
| 34 | ! |
|---|
| 35 | ! This subroutine is called by vdifc.F if the parameters "regolithe" and "water" are set to true in |
|---|
| 36 | ! callphys.def. It returns the flux of water vapor that is injected into or out of the regolith, |
|---|
| 37 | ! and which serves as a boundary condition to calculate the vertical atmospheric profile of water vapor |
|---|
| 38 | ! in vdifc. It also calculates the exchange between the subsurface and surface ice. |
|---|
| 39 | ! |
|---|
| 40 | ! It requires three subsurface tracers to be defined in traceur.def : |
|---|
| 41 | ! - h2o_vap_soil (subsurface water vapor) |
|---|
| 42 | ! - h2o_ice_soil (subsurface ice) |
|---|
| 43 | ! - h2o_ads_vap (adsorbed water) |
|---|
| 44 | ! |
|---|
| 45 | ! Author : Pierre - Yves Meslin (pmeslin@irap.omp.eu) |
|---|
| 46 | ! |
|---|
| 47 | ! ================================================================================================= = |
|---|
| 48 | |
|---|
| 49 | ! Libraries : |
|---|
| 50 | ! =========== = |
|---|
| 51 | !#include "dimensions.h" |
|---|
| 52 | !#include "dimphys.h" |
|---|
| 53 | !#include "comsoil.h" |
|---|
| 54 | #include "callkeys.h" |
|---|
| 55 | !#include "comcstfi.h" |
|---|
| 56 | !#include "tracer.h" |
|---|
| 57 | !#include "watercap.h" |
|---|
| 58 | |
|---|
| 59 | |
|---|
| 60 | ! Arguments : |
|---|
| 61 | ! ====================================================================== |
|---|
| 62 | |
|---|
| 63 | ! Inputs : |
|---|
| 64 | ! ---------------------------------------------------------------------- |
|---|
| 65 | integer, intent(in) :: ngrid ! number of points in the model (all lat and long point in one 1D array) |
|---|
| 66 | integer, intent(in) :: nlayer ! number of layers |
|---|
| 67 | integer, intent(in) :: nq ! number of tracers |
|---|
| 68 | integer, intent(in) :: nsoil |
|---|
| 69 | integer, intent(in) :: nqsoil |
|---|
| 70 | logical, save :: firstcall_soil = .true. ! triggers the initialization |
|---|
| 71 | real, intent(in) :: ptsoil(ngrid, nsoil) ! Subsurface temperatures |
|---|
| 72 | real, intent(in) :: ptsrf(ngrid) ! Surface temperature |
|---|
| 73 | real, intent(in) :: ptimestep ! length of the timestep (unit depends on run.def file) |
|---|
| 74 | logical, intent(in) :: exchange(ngrid) ! logical :: array describing whether there is exchange with the atmosphere at the current timestep |
|---|
| 75 | |
|---|
| 76 | real, intent(in) :: qsat_surf(ngrid) ! Saturation mixing ratio at the surface at the current surface temperature (formerly qsat) |
|---|
| 77 | real, intent(in) :: pq(ngrid, nlayer, nq) ! Tracer atmospheric mixing ratio |
|---|
| 78 | real, intent(in) :: pa(ngrid, nlayer) ! Coefficients za |
|---|
| 79 | real, intent(in) :: pb(ngrid, nlayer) ! Coefficients zb |
|---|
| 80 | real, intent(in) :: pc(ngrid, nlayer) ! Coefficients zc |
|---|
| 81 | real, intent(in) :: pd(ngrid, nlayer) ! Coefficients zd |
|---|
| 82 | real, intent(in) :: pdqsdifpot(ngrid) ! Potential pdqsdif (without subsurface - atmosphere exchange) |
|---|
| 83 | |
|---|
| 84 | real, intent(in) :: pplev(ngrid, nlayer+1) ! Atmospheric pressure levels |
|---|
| 85 | real, intent(in) :: rhoatmo(ngrid) ! Atmospheric air density (1st layer) (not used right now) |
|---|
| 86 | logical, intent(in) :: writeoutput ! just to check we are at the last subtimestep and we |
|---|
| 87 | |
|---|
| 88 | ! Variables modified : |
|---|
| 89 | ! ---------------------------------------------------------------------- |
|---|
| 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) |
|---|
| 93 | ! Outputs : |
|---|
| 94 | ! ---------------------------------------------------------------------- |
|---|
| 95 | real, intent(out) :: zdqsdifrego(ngrid) ! Flux from subsurface (positive pointing outwards) |
|---|
| 96 | real, intent(out) :: zq1temp2(ngrid) ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg) |
|---|
| 97 | real*8, intent(out) :: saturation_water_ice(ngrid, nsoil) ! Water pore ice saturation level (formerly Sw) |
|---|
| 98 | |
|---|
| 99 | ! Outputs for the output files |
|---|
| 100 | real*8 :: preduite(ngrid, nsoil) ! how close the water vapor density is to adsorption saturation |
|---|
| 101 | integer :: exch(ngrid) ! translates the logical variable exchange into a -1 or 1 |
|---|
| 102 | real*8 :: mass_h2o(ngrid) ! Mass of subsurface water column at each point (kg.m - 2) (formerly mh2otot) |
|---|
| 103 | real*8 :: mass_ice(ngrid) ! Mass of subsurface ice at each point (kg.m - 2) (formerly micetot) |
|---|
| 104 | real*8 :: mass_h2o_tot ! Mass of subsurface water over the whole planet |
|---|
| 105 | real*8 :: mass_ice_tot ! Mass of subsurface ice over the whole planet |
|---|
| 106 | real*8 :: nsurf(ngrid) ! surface tracer density in kg/m^3 |
|---|
| 107 | real*8 :: close_out(ngrid, nsoil) ! output for close_top and close_bottom |
|---|
| 108 | |
|---|
| 109 | ! Local (saved) variables : |
|---|
| 110 | ! ====================================================================== |
|---|
| 111 | |
|---|
| 112 | real*8 :: P_sat_soil(ngrid, nsoil) ! water saturation pressure of subsurface cells (at miD-layers) (formerly Psatsoil) |
|---|
| 113 | real*8 :: nsatsoil(ngrid, nsoil) ! gas number density at saturation pressure |
|---|
| 114 | real*8, allocatable, save :: znsoil(:, :) ! Water vapor soil concentration (kg.m - 3 of pore air) |
|---|
| 115 | real*8 :: znsoilprev(ngrid, nsoil) ! value of znsoil in the previous timestep |
|---|
| 116 | real*8 :: znsoilprev2(ngrid, nsoil) ! second variable for the value of znsoil in the previous timestep |
|---|
| 117 | real*8 :: zdqsoil(ngrid, nsoil) ! Increase of pqsoil if sublimation of ice |
|---|
| 118 | real*8, allocatable, save :: ice(:, :) ! Ice content of subsurface cells (at miD-layers) (kg / m3) |
|---|
| 119 | real*8 :: iceprev(ngrid, nsoil) |
|---|
| 120 | logical :: ice_index_flag(nsoil) ! flag for ice presence |
|---|
| 121 | real*8, allocatable, save :: adswater(:, :) ! Adsorbed water in subsurface cells (at miD-layers) (...) |
|---|
| 122 | real*8 :: adswater_temp(ngrid, nsoil) ! temprory variable for adsorbed water |
|---|
| 123 | logical, allocatable, save :: over_mono_sat_flag(:, :) ! Formerly ads_water_saturation_flag_2(nsoil) (see descritption of the variable recompute_cell_ads_flag for an explanation ! ARNAU |
|---|
| 124 | real*8 :: adswprev(ngrid, nsoil) |
|---|
| 125 | logical :: recompute_cell_ads_flag(nsoil) ! ARNAU |
|---|
| 126 | ! Formerly ads_water_saturation_flag_1 but with a twist. This variable now |
|---|
| 127 | ! checks whether coefficients have changed and need recomputing. If the cell |
|---|
| 128 | ! is seen to be over the monolayer saturation level (i.e. the cell fulfills the |
|---|
| 129 | ! condition adswater_temp(ig, ik) > adswater_sat_mono) but the coefficients |
|---|
| 130 | ! have been computed assuming that the cell is below the monolayer saturation |
|---|
| 131 | ! layer (i.e. the variable over_mono_sat_flag(ig, ik) had been set to .false.), |
|---|
| 132 | ! then the cell needs to have its coefficients recomputed according to the |
|---|
| 133 | ! previous condition: adswater_temp(ig, ik) > adswater_sat_mono. Then, |
|---|
| 134 | ! the variable recompute_cell_ads_flag(ik) becomes .true.. ! ARNAU |
|---|
| 135 | real*8, save :: adswater_sat_mono ! Adsorption monolayer saturation value (kg.m - 3 of regolith) ! ARNAU |
|---|
| 136 | real*8 :: delta0(ngrid) ! Coefficient delta(0) modified |
|---|
| 137 | real*8 :: alpha0(ngrid) |
|---|
| 138 | real*8 :: beta0(ngrid) |
|---|
| 139 | |
|---|
| 140 | ! Adsorbtion/Desorption variables and parameters |
|---|
| 141 | real*8 :: Ka(ngrid, nsoil) ! Adsorption time constant (s - 1) before monolayer saturation (first segment of the bilinear function) |
|---|
| 142 | real*8 :: Kd(ngrid, nsoil) ! Desorption time constant (s - 1) before monolayer saturation (first segment of the bilinear function) |
|---|
| 143 | real*8 :: k_ads_eq(ngrid, nsoil) ! Equilibrium adsorption coefficient (formerly kads) before monolayer saturation (first segment of the bilinear function); unitless (converts kg/m3 into kg/m3) |
|---|
| 144 | real*8 :: Ka2(ngrid, nsoil) ! Adsorption time constant (s - 1) after monolayer saturation (second segment of the bilinear function) ! modified 2020 |
|---|
| 145 | real*8 :: Kd2(ngrid, nsoil) ! Desorption time constant (s - 1) after monolayer saturation (second segment of the bilinear function) ! modified 2020 |
|---|
| 146 | real*8 :: k_ads_eq2(ngrid, nsoil) ! Equilibrium adsorption coefficient (formerly kads) after monolayer saturation (second segment of the bilinear function); unitless ! modified 2020 |
|---|
| 147 | real*8 :: c0(ngrid, nsoil) ! Intercept of the second line in the bilinear function ! modified 2020 |
|---|
| 148 | real*8, parameter :: kinetic_factor = 0.01 ! (inverse of) Characteristic time to reach adsorption equilibrium (s - 1): |
|---|
| 149 | ! fixed arbitrarily when kinetics factors are unknown |
|---|
| 150 | ! possible values: ! = 1 / 1800 s - 1 ! / 1.16D-6 / !( = 10 earth days)! / 1.D0 / ! / 1.2D-5 / ! |
|---|
| 151 | |
|---|
| 152 | real*8, allocatable, save :: ztot1(:, :) ! Total (water vapor + ice) content (kg.m - 3 of soil) |
|---|
| 153 | real*8 :: dztot1(nsoil) |
|---|
| 154 | real*8 :: h2otot(ngrid, nsoil) ! Total water content (ice + water vapor + adsorbed water) |
|---|
| 155 | real*8 :: flux(ngrid, 0:nsoil - 1) ! Fluxes at interfaces (kg.m - 2.s - 1) (positive = upward) |
|---|
| 156 | real*8 :: rho(ngrid) ! Air density (first subsurface layer) |
|---|
| 157 | real*8 :: rhosurf(ngrid) ! Surface air density |
|---|
| 158 | |
|---|
| 159 | |
|---|
| 160 | ! Porosity and tortuosity variables |
|---|
| 161 | real*8, allocatable, save :: porosity_ice_free(:, :) ! Porosity with no ice present (formerly eps0) |
|---|
| 162 | real*8, allocatable, save :: porosity(:, :) ! Porosity (formerly eps) |
|---|
| 163 | real*8 :: porosity_prev(ngrid, nsoil) ! Porosity from previous timestep |
|---|
| 164 | real*8, allocatable, save :: tortuosity(:, :) ! Tortuosity factor (formerly tortuo) |
|---|
| 165 | |
|---|
| 166 | real*8 :: saturation_water_ice_inter(ngrid, nsoil) ! Water pore ice saturation level at the interlayer |
|---|
| 167 | |
|---|
| 168 | ! Diffussion coefficients |
|---|
| 169 | real*8 :: D_mid(ngrid, nsoil) ! Midlayer diffusion coefficients |
|---|
| 170 | real*8 :: D_inter(ngrid, nsoil) ! Interlayer diffusion coefficients (formerly D) |
|---|
| 171 | real*8, allocatable, save :: D0(:, :) ! Diffusion coefficient prefactor : |
|---|
| 172 | ! If (medium = 1) : D0 = value of D_mid for saturation_water_ice = 0, ie. poro / tortuo * Dm (in m2 / s) |
|---|
| 173 | ! If (medium = 2) : D0 = porosity_tortuosity factor (dimensionless) |
|---|
| 174 | real*8 :: omega(ngrid, nsoil) ! H2O - CO2 collision integral |
|---|
| 175 | real*8 :: vth(ngrid, nsoil) ! H2O thermal speed |
|---|
| 176 | real*8, parameter :: Dk0 = 0.459D0 ! Knudsen diffusion coefficient (for saturation_water_ice = 0) Value for a (5 / 1) bidispersed random assembly of spheres |
|---|
| 177 | ! (dimensionless, ie. divided by thermal speed and characteristic meshsize of the porous medium) |
|---|
| 178 | real*8 :: Dk(ngrid, nsoil) ! Knudsen diffusion coefficient (divided by porosity / tortuosity factor) |
|---|
| 179 | real*8 :: Dk_inter(ngrid, nsoil) ! Knudsen diffusion coefficient at the interlayer |
|---|
| 180 | real*8 :: Dm(ngrid, nsoil) ! Molecular diffusion coefficient |
|---|
| 181 | real*8 :: Dm_inter(ngrid, nsoil) ! Molecular diffusion coefficient at the interlayer |
|---|
| 182 | |
|---|
| 183 | real*8, parameter :: choke_fraction = 0.8D0 ! fraction of ice that prevents further diffusion |
|---|
| 184 | logical, allocatable, save :: close_top(:, :) ! closing diffusion at the top of a layer if ice rises over saturation |
|---|
| 185 | logical, allocatable, save :: close_bottom(:, :) ! closing diffusion at the bottom of a layer if ice risies over saturation |
|---|
| 186 | logical, parameter :: print_closure_warnings = .true. |
|---|
| 187 | |
|---|
| 188 | ! Coefficients for the Diffusion calculations |
|---|
| 189 | real*8 :: A(ngrid, nsoil) ! Coefficient in the diffusion formula |
|---|
| 190 | real*8 :: B(ngrid, nsoil) ! Coefficient in the diffusion formula |
|---|
| 191 | real*8 :: C(ngrid, nsoil) ! Coefficient in the diffusion formula |
|---|
| 192 | real*8 :: E(ngrid, nsoil) ! Coefficient in the diffusion formula (before monolayer saturation) ! added 2020 |
|---|
| 193 | real*8 :: F(ngrid, nsoil) ! Coefficient in the diffusion formula (before monolayer saturation) ! added 2020 |
|---|
| 194 | real*8 :: E2(ngrid, nsoil) ! Coefficient in the diffusion formula (after monolayer saturation) ! added 2020 |
|---|
| 195 | real*8 :: F2(ngrid, nsoil) ! Coefficient in the diffusion formula (after monolayer saturation) ! added 2020 |
|---|
| 196 | real*8, allocatable, save :: zalpha(:, :) ! Alpha coefficients |
|---|
| 197 | real*8, allocatable, save :: zbeta(:, :) ! Beta coefficients |
|---|
| 198 | real*8 :: zdelta(ngrid, nsoil - 1) ! Delta coefficients |
|---|
| 199 | |
|---|
| 200 | ! Distances between layers |
|---|
| 201 | real*8, allocatable, save :: interlayer_dz(:, :) ! Distance between the interlayer points in m (formerly interdz) |
|---|
| 202 | real*8, allocatable, save :: midlayer_dz(:, :) ! Distance between the midcell points in m (formerly middz) |
|---|
| 203 | |
|---|
| 204 | real*8 :: nsat(ngrid, nsoil) ! amount of water vapor at (adsorption) monolayer saturation ! modified 2020 |
|---|
| 205 | |
|---|
| 206 | real*8, allocatable, save :: meshsize(:, :) ! scaling/dimension factor for the por size |
|---|
| 207 | real*8, allocatable, save :: rho_soil(:, :) ! Soil density (bulk - kg / m3) (formerly rhosoil) |
|---|
| 208 | real*8, allocatable, save :: cste_ads(:, :) ! Prefactor of adsorption coeff. k |
|---|
| 209 | |
|---|
| 210 | ! general constants |
|---|
| 211 | real*8, parameter :: kb = 1.38065D-23 ! Boltzmann constant |
|---|
| 212 | real*8, parameter :: mw = 2.988D-26 ! Water molecular weight |
|---|
| 213 | !real*8, parameter :: rho_H2O_ice = 920.D0 ! Ice density |
|---|
| 214 | |
|---|
| 215 | ! adsorption coefficients |
|---|
| 216 | real*8, parameter :: enthalpy_ads = 35.D3 ! Enthalpy of adsorption (J.mole - 1) options 21.D3, 35.D3, and 60.D3 |
|---|
| 217 | real*8, parameter :: enthalpy_ads2 = 21.D3 ! Enthalpy of adsorption (J.mole - 1) options 21.D3, 35.D3, and 60.D3 for the second segment ! added 2020 |
|---|
| 218 | real*8, parameter :: DeltaQ_ads = 21.D3 ! 10.D3 ! This is the DeltaQ_ads = heat of adsorption - enthalpy of liquefaction/sublimation = E1 - EL and may be equal to RT*ln(c), where c is the BET constant (from BET1938). This is for the first segment (J.mole - 1) ! added 2020 |
|---|
| 219 | real*8, parameter :: DeltaQ_ads2 = 21.D3 ! 10.D3 ! This is the DeltaQ_ads = heat of adsorption - enthalpy of liquefaction/sublimation = E1 - EL and may be equal to RT*ln(c), where c is the BET constant (from BET1938). This is for the second segment (J.mole - 1) ! added 2020 |
|---|
| 220 | real*8, parameter :: tau0 = 1D-14 |
|---|
| 221 | real*8, parameter :: S = 17.D3 ! Soil specific surface area (m2.kg - 1 of solid) options: 17.D3 and 106.D3 |
|---|
| 222 | real*8, parameter:: Sm = 10.6D-20 ! Surface of the water molecule (m2) (only needed in the theoretical formula which is not used right now) |
|---|
| 223 | |
|---|
| 224 | |
|---|
| 225 | ! Reference values for choice_ads = 2 |
|---|
| 226 | real*8, parameter :: Tref = 243.15D0 |
|---|
| 227 | real*8, parameter :: nref = 2.31505631D-6 ! calculated as 18.D-3 / (8.314D0 * Tref) * 0.26D0 ! not used anymore (for the time being) |
|---|
| 228 | real*8, parameter :: Kd_ref = 0.0206D0 ! Not used for the time being (would require specific measurements to be known, but it is rarely measured) |
|---|
| 229 | real*8, parameter :: Ka_ref = 2.4D-4 ! Not used for the time being |
|---|
| 230 | ! 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 |
|---|
| 231 | ! real*8, parameter :: Kref2 = 5.95D-7 ! Value from data analysis Pommerol2009 ! ARNAU |
|---|
| 232 | real*8, parameter :: Kref = 0.205D-6 ! Value obtained from the fit of all adsorption data from Pommmerol (2009) (see Arnau's report - p.28: = yp/xp = 2.6998E-7/3.5562E-2, divided by psat(243.15K=37 Pa; will need to be multiplied by RT/M to be unitless before multiplying by znsoil, which in kg(water)/m3(air) and not in pascals) |
|---|
| 233 | real*8, parameter :: Kref2 = 0.108D-7 ! Value obtained from the fit of all adsorption data from Pommmerol (2009) (see Arnau's report - p.28 = m2; divided by psat(243.15K)=37 Pa) |
|---|
| 234 | |
|---|
| 235 | |
|---|
| 236 | logical :: recompute_all_cells_ads_flag ! Old saturation_flag but with a behaviour change ! ARNAU |
|---|
| 237 | ! The old saturation_flag was used to check whether the cell was saturated or |
|---|
| 238 | ! not, and therefore to assign the saturation value adswater(ig, ik) = |
|---|
| 239 | ! adswater_sat instead of the usual adswater(ig, ik) = adswater_temp(ig, ik) |
|---|
| 240 | ! The old routine using saturation_flag did not require that the A, B, C,... |
|---|
| 241 | ! adsorption coefficients be computed, but the new soilwater subroutine |
|---|
| 242 | ! does. Therefore, the variable recompute_all_cells_ads_flag checks whether |
|---|
| 243 | ! there is a cell of a column that requires recomputing. If at least one cell |
|---|
| 244 | ! requires recomputing (i.e. recompute_cell_ads_flag(ik) is .true.), then |
|---|
| 245 | ! recompute_all_cells_ads_flag becomes true, and the adsorption coefficients, |
|---|
| 246 | ! as well as the recursive equation (appearing in this code as `backward |
|---|
| 247 | ! loop over all layers`) ! ARNAU |
|---|
| 248 | logical :: sublimation_flag |
|---|
| 249 | logical :: condensation_flag(nsoil) |
|---|
| 250 | |
|---|
| 251 | ! variables and parameters for the H2O map |
|---|
| 252 | integer, parameter :: n_long_H2O = 66!180 ! number of longitudes of the new H2O map |
|---|
| 253 | integer, parameter :: n_lat_H2O = 50 !87 ! number of latitudes of the new H2O map |
|---|
| 254 | real*8, parameter :: rho_H2O_ice = 920.0D0 ! Ice density (formerly rhoice) |
|---|
| 255 | real :: latH2O(n_lat_H2O*n_long_H2O) ! Latitude at H2O map points |
|---|
| 256 | real :: longH2O(n_lat_H2O*n_long_H2O) ! Longitude at H2O map points |
|---|
| 257 | real :: H2O_depth_data(n_lat_H2O*n_long_H2O) ! depth of the ice layer in in g/cm^2 at H2O map points |
|---|
| 258 | real :: H2O_cover_data(n_lat_H2O*n_long_H2O) ! H2O content of the cover layer at H2O map points (not used yet) |
|---|
| 259 | real :: dataH2O(n_lat_H2O*n_long_H2O) ! H2O content of the ice layer at H2O map points |
|---|
| 260 | real :: latH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable |
|---|
| 261 | real :: longH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable |
|---|
| 262 | real :: dataH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable |
|---|
| 263 | real :: H2O_depth_data_temp(n_lat_H2O*n_long_H2O) ! intermediate variable |
|---|
| 264 | real, allocatable, save :: H2O(:) ! H2O map interpolated at GCM grid points (in wt%) |
|---|
| 265 | real, allocatable, save :: H2O_depth(:) ! H2O map ice depth interpolated at GCM in g/cm^2 |
|---|
| 266 | |
|---|
| 267 | ! variables for the 1D case |
|---|
| 268 | real*8, parameter :: mmr_h2o = 0.6D-4 ! Water vapor mass mixing ratio (for initialization only) |
|---|
| 269 | real*8 :: diff(ngrid, nsoil) ! difference between total water content and ice + vapor content |
|---|
| 270 | ! (only used for output, inconstistent: should just be adswater) |
|---|
| 271 | real :: var_flux_soil(ngrid, nsoil) ! Output of the flux between soil layers (kg/m^3/s) |
|---|
| 272 | real :: default_diffcoeff = 4e-4 ! Diffusion coefficient by default (no variations with Temperature and pressure (m/s^2) |
|---|
| 273 | real :: tol_massiveice = 50. ! Tolerance to account for massive ice (kg/m^3) |
|---|
| 274 | ! Loop variables and counters |
|---|
| 275 | integer :: ig, ik, i, j ! loop variables |
|---|
| 276 | logical :: output_trigger ! used to limit amount of written output |
|---|
| 277 | integer, save :: n ! number of runs of this subroutine |
|---|
| 278 | integer :: sublim_n ! number of sublimation loop runs |
|---|
| 279 | integer :: satur_mono_n ! number of monolayer saturation loop runs ! added 2020 |
|---|
| 280 | |
|---|
| 281 | |
|---|
| 282 | if (.not. allocated(znsoil)) then |
|---|
| 283 | allocate( znsoil(ngrid, nsoil) ) |
|---|
| 284 | allocate( ice(ngrid, nsoil) ) |
|---|
| 285 | allocate( adswater(ngrid, nsoil) ) |
|---|
| 286 | allocate( ztot1(ngrid, nsoil) ) |
|---|
| 287 | allocate( porosity_ice_free(ngrid, nsoil) ) |
|---|
| 288 | allocate( porosity(ngrid, nsoil) ) |
|---|
| 289 | allocate( tortuosity(ngrid, nsoil) ) |
|---|
| 290 | allocate( D0(ngrid, nsoil) ) |
|---|
| 291 | allocate( interlayer_dz(ngrid, nsoil) ) |
|---|
| 292 | allocate( midlayer_dz(ngrid, 0:nsoil) ) |
|---|
| 293 | ! allocate( zalpha(ngrid, nsoil-1) ) |
|---|
| 294 | ! allocate( zbeta(ngrid, nsoil-1) ) |
|---|
| 295 | allocate( zalpha(ngrid, nsoil) ) ! extra entry to match output format |
|---|
| 296 | allocate( zbeta(ngrid, nsoil) ) ! extra entry to match output format |
|---|
| 297 | allocate( meshsize(ngrid, nsoil) ) |
|---|
| 298 | allocate( rho_soil(ngrid, nsoil) ) |
|---|
| 299 | allocate( cste_ads(ngrid, nsoil) ) |
|---|
| 300 | allocate( H2O(ngrid) ) |
|---|
| 301 | allocate( H2O_depth(ngrid) ) |
|---|
| 302 | allocate( close_top(ngrid, nsoil) ) |
|---|
| 303 | allocate( close_bottom(ngrid, nsoil) ) |
|---|
| 304 | allocate( over_mono_sat_flag(ngrid, nsoil) ) |
|---|
| 305 | endif |
|---|
| 306 | |
|---|
| 307 | ! Used commons |
|---|
| 308 | ! mugaz ! Molar mass of the atmosphere (g.mol-1) ~43.49 |
|---|
| 309 | |
|---|
| 310 | |
|---|
| 311 | |
|---|
| 312 | ! Initialisation : |
|---|
| 313 | ! ================= |
|---|
| 314 | |
|---|
| 315 | if (firstcall_soil) then |
|---|
| 316 | n = 0 |
|---|
| 317 | close_top = .false. |
|---|
| 318 | close_bottom = .false. |
|---|
| 319 | print *, "adsorption enthalpy first segment: ", enthalpy_ads |
|---|
| 320 | print *, "adsorption enthalpy second segment: ", enthalpy_ads2 |
|---|
| 321 | print *, "adsorption BET constant C first segment: ", DeltaQ_ads |
|---|
| 322 | print *, "adsorption BET constant C second segment: ", DeltaQ_ads2 |
|---|
| 323 | print *, "specific surface area: ", S |
|---|
| 324 | |
|---|
| 325 | do ig = 1, ngrid |
|---|
| 326 | ! if(.not.watercaptag(ig)) then |
|---|
| 327 | |
|---|
| 328 | ! Initialization of soil parameters |
|---|
| 329 | ! =================================== |
|---|
| 330 | |
|---|
| 331 | ! initialize the midlayer distances |
|---|
| 332 | midlayer_dz(ig, 0) = mlayer(0) |
|---|
| 333 | |
|---|
| 334 | do ik = 1, nsoil - 1 |
|---|
| 335 | midlayer_dz(ig, ik) = mlayer(ik) - mlayer(ik - 1) |
|---|
| 336 | enddo |
|---|
| 337 | |
|---|
| 338 | ! initialize the interlayer distances |
|---|
| 339 | interlayer_dz(ig, 1) = layer(1) |
|---|
| 340 | do ik = 2, nsoil |
|---|
| 341 | interlayer_dz(ig, ik) = layer(ik) - layer(ik - 1) |
|---|
| 342 | enddo |
|---|
| 343 | |
|---|
| 344 | ! Choice of porous medium and D0 |
|---|
| 345 | ! =============================== = |
|---|
| 346 | do ik = 1, nsoil |
|---|
| 347 | |
|---|
| 348 | ! These properties are defined here in order to enable custom profiles |
|---|
| 349 | if(ads_massive_ice) then |
|---|
| 350 | if(pqsoil(ig, ik, igcm_h2o_ice_soil).gt.tol_massiveice) then |
|---|
| 351 | porosity_ice_free(ig, ik) = 0.999999 |
|---|
| 352 | else |
|---|
| 353 | porosity_ice_free(ig, ik) = porosity_reg |
|---|
| 354 | endif |
|---|
| 355 | else |
|---|
| 356 | porosity_ice_free(ig, ik) = porosity_reg |
|---|
| 357 | endif |
|---|
| 358 | tortuosity(ig, ik) = 1.5D0 |
|---|
| 359 | rho_soil(ig, ik) = 1.3D3 ! in kg/m3 of regolith (incl. porosity) |
|---|
| 360 | |
|---|
| 361 | meshsize(ig, ik) = 5.D-6 ! Characteristic size 1e - 6 = 1 micron : diameter(mode 5) = 10 microns, diameter(mode 1) = 1 microns |
|---|
| 362 | ! Example : with meshsize = 5.E - 6 : grain sizes = 50 and 10 microns |
|---|
| 363 | D0(ig, ik) = porosity_ice_free(ig, ik) / tortuosity(ig, ik) |
|---|
| 364 | |
|---|
| 365 | enddo |
|---|
| 366 | |
|---|
| 367 | ! Choice of adsorption coefficients |
|---|
| 368 | ! =================================== |
|---|
| 369 | 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 |
|---|
| 370 | ! 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) |
|---|
| 371 | ! theoretical formula is = rho_soil(ig, 1) * S / Sm * mw |
|---|
| 372 | |
|---|
| 373 | ! Numerical values above are for SSA = 13.7 m2 / g and plateau = 0.8wt% |
|---|
| 374 | ! adswater_sat_mono = 6D-2 * rho_soil(ig, 1) ! Experimental value for JSC Mars - 1 (Pommerol et al., 2009) |
|---|
| 375 | ! with SSA = 106 m2 / g and plateau (net) = 6wt% |
|---|
| 376 | |
|---|
| 377 | ! Initialisation of ice, water vapor and adsorbed water profiles |
|---|
| 378 | ! =============================================================== = |
|---|
| 379 | do ik = 1, nsoil |
|---|
| 380 | ! Initialisation of pqsoil(t = 0) |
|---|
| 381 | ! in 1D simulations the initialization happens here |
|---|
| 382 | ! if (ngrid.eq.1) then |
|---|
| 383 | ! znsoil(ig, ik) = mmr_h2o * mugaz * 1.D-3 * pplev(ig, 1) & |
|---|
| 384 | ! / (8.314D0 * ptsoil(ig, nsoil - 4)) |
|---|
| 385 | ! ! znsoil(ig, ik) = 0. |
|---|
| 386 | ! ice(ig, ik) = 0.D0 ! 0.1D0 !414.D0 |
|---|
| 387 | |
|---|
| 388 | ! saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) |
|---|
| 389 | ! porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) |
|---|
| 390 | |
|---|
| 391 | ! if (choice_ads.eq.1) then |
|---|
| 392 | ! vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, nsoil - 4) & |
|---|
| 393 | ! / (pi * 18.D-3)) |
|---|
| 394 | |
|---|
| 395 | ! first segment of the bilinear function |
|---|
| 396 | ! k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & |
|---|
| 397 | ! / (dexp(-enthalpy_ads & |
|---|
| 398 | ! / (8.314D0 * ptsoil(ig, nsoil - 4))) / tau0) |
|---|
| 399 | ! Kd(ig, ik) = kinetic_factor & |
|---|
| 400 | ! / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 401 | ! Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / & |
|---|
| 402 | ! (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 403 | |
|---|
| 404 | ! second segment of the bilinear function ! added 2020 |
|---|
| 405 | ! k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & |
|---|
| 406 | ! / (dexp(-enthalpy_ads2 & |
|---|
| 407 | ! / (8.314D0 * ptsoil(ig, nsoil - 4))) / tau0) |
|---|
| 408 | ! Kd2(ig, ik) = kinetic_factor & |
|---|
| 409 | ! / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 410 | ! Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / & |
|---|
| 411 | ! (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 412 | ! |
|---|
| 413 | ! elseif (choice_ads.eq.2) then ! par rapport \E0 valeur experimentale de reference |
|---|
| 414 | ! ! Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 * |
|---|
| 415 | ! ! & (1. / ptsoil(ig, nsoil - 4) - 1. / Tref)) |
|---|
| 416 | ! ! Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref * |
|---|
| 417 | ! ! & sqrt(ptsoil(ig, nsoil - 4) / Tref) / nref! |
|---|
| 418 | ! |
|---|
| 419 | ! ! first segment of the bilinear function |
|---|
| 420 | ! 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 |
|---|
| 421 | ! |
|---|
| 422 | ! Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 423 | ! |
|---|
| 424 | ! Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 425 | ! |
|---|
| 426 | ! ! second segment of the bilinear function ! added 2020 |
|---|
| 427 | ! 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 |
|---|
| 428 | ! |
|---|
| 429 | ! Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 430 | ! |
|---|
| 431 | ! Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 432 | ! endif |
|---|
| 433 | |
|---|
| 434 | ! if (k_ads_eq(ig,ik) * znsoil(ig, ik) .gt. adswater_sat_mono) then ! modified 2024, after correction of c0 expression |
|---|
| 435 | ! c0(ig, ik) = ( 1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono |
|---|
| 436 | ! adswater(ig, ik) = c0(ig,ik) + k_ads_eq2(ig, ik) * znsoil(ig, ik) |
|---|
| 437 | ! else |
|---|
| 438 | ! adswater(ig, ik) = k_ads_eq(ig, ik) * znsoil(ig, ik) |
|---|
| 439 | ! endif |
|---|
| 440 | |
|---|
| 441 | ! else ! in 3D simulations initialisation happens with newstart.F |
|---|
| 442 | znsoil(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_soil) |
|---|
| 443 | ice(ig, ik) = pqsoil(ig, ik, igcm_h2o_ice_soil) |
|---|
| 444 | adswater(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_ads) |
|---|
| 445 | ! endif |
|---|
| 446 | |
|---|
| 447 | saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) |
|---|
| 448 | |
|---|
| 449 | if (saturation_water_ice(ig, ik).lt.0.D0) then |
|---|
| 450 | print *, "WARNING!! soil water ice negative at ", ig, ik |
|---|
| 451 | print *, "saturation value: ", saturation_water_ice(ig, ik) |
|---|
| 452 | print *, "setting saturation to 0" |
|---|
| 453 | saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) |
|---|
| 454 | endif |
|---|
| 455 | |
|---|
| 456 | porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) |
|---|
| 457 | |
|---|
| 458 | ! Initialise soil as if all points where below the monolayer saturation level (added 2020) => now depends on value of adswater (modified in 2024) |
|---|
| 459 | if (adswater(ig,ik).gt.adswater_sat_mono) then |
|---|
| 460 | over_mono_sat_flag(ig, ik) = .true. ! |
|---|
| 461 | else |
|---|
| 462 | over_mono_sat_flag(ig, ik) = .false. |
|---|
| 463 | endif |
|---|
| 464 | |
|---|
| 465 | enddo |
|---|
| 466 | ! endif |
|---|
| 467 | enddo |
|---|
| 468 | |
|---|
| 469 | print *, "initializing H2O data" |
|---|
| 470 | ! |
|---|
| 471 | ! ! 1.6 intitalization of the water content |
|---|
| 472 | ! open(40,file='H2O_data') |
|---|
| 473 | ! do i = 1, n_lat_H2O*n_long_H2O |
|---|
| 474 | ! read(40,*) latH2O_temp(i), longH2O_temp(i), H2O_cover_data(i), dataH2O_temp(i), H2O_depth_data_temp(i) |
|---|
| 475 | ! ! write(*, *) 'depth data ', H2O_depth_data(i) |
|---|
| 476 | ! enddo |
|---|
| 477 | ! close(40) |
|---|
| 478 | ! |
|---|
| 479 | ! ! 1.6.2. put latH2O, longH2O and dataH2O in the right format to pass on to mapTh |
|---|
| 480 | ! ! in the datafile the latitudes are listed from negative to positive, but mapTh needs them the other way around |
|---|
| 481 | ! do i = 0, n_lat_H2O - 1 |
|---|
| 482 | ! do j = 1, n_long_H2O |
|---|
| 483 | ! latH2O( n_long_H2O * i + j ) = latH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j) |
|---|
| 484 | ! longH2O( n_long_H2O * i + j ) = longH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j) |
|---|
| 485 | ! dataH2O( n_long_H2O * i + j ) = dataH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j) |
|---|
| 486 | ! H2O_depth_data( n_long_H2O * i + j ) = H2O_depth_data_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j) |
|---|
| 487 | ! enddo |
|---|
| 488 | ! enddo |
|---|
| 489 | ! |
|---|
| 490 | ! call mapTh(latH2O, longH2O, dataH2O, n_long_H2O, n_lat_H2O, ngrid, H2O) |
|---|
| 491 | |
|---|
| 492 | ! call mapTh(latH2O, longH2O, H2O_depth_data, n_long_H2O, n_lat_H2O, ngrid, H2O_depth) |
|---|
| 493 | ! |
|---|
| 494 | do ig = 1, ngrid |
|---|
| 495 | ! convert depth from g/cm^2 to m |
|---|
| 496 | ! print *, H2O_depth(ig), rho_soil(ig, 1) |
|---|
| 497 | ! H2O_depth(ig) = H2O_depth(ig) * 10 / rho_soil(ig, 1) |
|---|
| 498 | ! if ( (latitude_deg(ig).lt.40).and.(latitude_deg(ig).gt.-40)) then |
|---|
| 499 | ! H2O(ig) = 0. |
|---|
| 500 | ! endif |
|---|
| 501 | |
|---|
| 502 | output_trigger = .true. |
|---|
| 503 | do ik = 1, nsoil |
|---|
| 504 | ! if (H2O_depth(ig).le.layer(ik)) then |
|---|
| 505 | ! ice(ig, ik) = min(H2O(ig), rho_H2O_ice * porosity_ice_free(ig, ik)) |
|---|
| 506 | ! if (output_trigger) then |
|---|
| 507 | ! print *, "ice set at: ", ig, ik, "to :", ice(ig,ik), "depth: ", H2O_depth(ig) |
|---|
| 508 | ! output_trigger = .false. |
|---|
| 509 | ! endif |
|---|
| 510 | ! endif |
|---|
| 511 | |
|---|
| 512 | saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) |
|---|
| 513 | saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) |
|---|
| 514 | porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) |
|---|
| 515 | enddo |
|---|
| 516 | enddo |
|---|
| 517 | |
|---|
| 518 | print *, "Kinetic factor: ", kinetic_factor |
|---|
| 519 | if (kinetic_factor.eq.0) then |
|---|
| 520 | print *, "WARNING: adsorption is switched off" |
|---|
| 521 | endif |
|---|
| 522 | firstcall_soil = .false. |
|---|
| 523 | endif ! of "if firstcall_soil" |
|---|
| 524 | |
|---|
| 525 | |
|---|
| 526 | |
|---|
| 527 | |
|---|
| 528 | ! Update properties in case of the sublimation of massive ice |
|---|
| 529 | |
|---|
| 530 | if(ads_massive_ice) then |
|---|
| 531 | do ig = 1, ngrid |
|---|
| 532 | do ik = 1, nsoil |
|---|
| 533 | if(pqsoil(ig, ik, igcm_h2o_ice_soil).gt.tol_massiveice) then |
|---|
| 534 | porosity_ice_free(ig, ik) = 0.999999 |
|---|
| 535 | else |
|---|
| 536 | porosity_ice_free(ig, ik) = porosity_reg |
|---|
| 537 | endif |
|---|
| 538 | enddo |
|---|
| 539 | enddo |
|---|
| 540 | endif |
|---|
| 541 | |
|---|
| 542 | ! Computation of new (new step) transport coefficients : |
|---|
| 543 | ! ======================================================= = |
|---|
| 544 | |
|---|
| 545 | do ig = 1, ngrid ! loop over all gridpoints |
|---|
| 546 | if(.not.watercaptag(ig)) then ! if there is no polar cap |
|---|
| 547 | |
|---|
| 548 | do ik = 1, nsoil |
|---|
| 549 | |
|---|
| 550 | ! calculate Water saturation level (pore volume filled by ice / total pore volume) |
|---|
| 551 | saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) |
|---|
| 552 | |
|---|
| 553 | if (saturation_water_ice(ig, ik).lt.0.D0) then |
|---|
| 554 | print *, "WARNING!! soil water ice negative at ", ig, ik |
|---|
| 555 | print *, "saturation value: ", saturation_water_ice(ig, ik) |
|---|
| 556 | print *, "setting saturation to 0" |
|---|
| 557 | saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) |
|---|
| 558 | endif |
|---|
| 559 | |
|---|
| 560 | ! save the old porosity |
|---|
| 561 | porosity_prev(ig, ik) = porosity(ig, ik) |
|---|
| 562 | |
|---|
| 563 | ! calculate the new porosity |
|---|
| 564 | porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) |
|---|
| 565 | ! Note : saturation_water_ice and porosity are computed from the ice content of the previous timestep |
|---|
| 566 | enddo |
|---|
| 567 | |
|---|
| 568 | ! calculate the air density in the first subsurface layer |
|---|
| 569 | rho(ig) = pplev(ig, 1) / (r * ptsoil(ig, 1)) |
|---|
| 570 | |
|---|
| 571 | ! calculate diffusion coefficients at mid- and interlayers (with ice content from the previous timestep) |
|---|
| 572 | ! Dependence on saturation_water_ice taken from the relation obtained by Hudson et al. (2009) and Meslin et al. (SSSAJ, 2010) |
|---|
| 573 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 574 | |
|---|
| 575 | ! calculate H2O mean thermal velocity |
|---|
| 576 | vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, ik) / (pi * 18.D-3)) |
|---|
| 577 | |
|---|
| 578 | ! The diffusion coefficient is calculated |
|---|
| 579 | |
|---|
| 580 | ! calculate the H2O - CO2 collision integral (after Mellon and Jakosky, JGR, 1993) |
|---|
| 581 | omega(ig, ik) = 1.442D0 - 0.673D0 * dlog(2.516D-3 * ptsoil(ig, ik)) & |
|---|
| 582 | + 0.252D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 2.D0 & |
|---|
| 583 | - 0.047D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 3.D0 & |
|---|
| 584 | + 0.003D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 4.D0 |
|---|
| 585 | |
|---|
| 586 | ! calculate the molecular diffusion coefficient |
|---|
| 587 | Dm(ig, ik) = 4.867D0 * ptsoil(ig, ik) ** (3.D0 / 2.D0) & |
|---|
| 588 | / (pplev(ig, 1) * omega(ig, ik)) * 1.D-4 |
|---|
| 589 | |
|---|
| 590 | ! calculate the Knudsen diffusion coefficient (divided by D0 = porosity / tortuosity in this case) |
|---|
| 591 | Dk(ig, ik) = Dk0 / D0(ig, ik) * meshsize(ig, ik) * vth(ig, ik) |
|---|
| 592 | |
|---|
| 593 | ! calculate the midlayer diffusion coeff. (with dependence on saturation_water_ice from Meslin et al., 2010 - only exact for normal diffusion) |
|---|
| 594 | D_mid(ig, ik) = D0(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) ** 2.D0 * 1.D0 / (1.D0 / Dm(ig, ik) + 1.D0 / Dk(ig, ik)) |
|---|
| 595 | |
|---|
| 596 | if(ads_const_D) D_mid(ig, ik) = default_diffcoeff |
|---|
| 597 | |
|---|
| 598 | ! Midlayer diffusion coeff. (correlation by Rogers and Nielson, 1991) |
|---|
| 599 | ! D_mid(ig, ik) = D0(ig, ik) * (1. - saturation_water_ice(ig, ik)) * exp(-6. * saturation_water_ice(ig, ik) & |
|---|
| 600 | ! * porosity_ice_free(ig, ik) - 6. * saturation_water_ice(ig, ik) ** (14. * porosity_ice_free(ig, ik))) & |
|---|
| 601 | ! * 1. / (1. / Dm(ig, ik) + 1. / Dk(ig, ik)) |
|---|
| 602 | enddo |
|---|
| 603 | |
|---|
| 604 | ! calculate the interlayer diffusion coefficient as interpolation of the midlayer coefficients |
|---|
| 605 | do ik = 1, nsoil - 1 |
|---|
| 606 | Dm_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * Dm(ig, ik + 1) & |
|---|
| 607 | + (mlayer(ik) - layer(ik)) * Dm(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1)) |
|---|
| 608 | |
|---|
| 609 | Dk_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * Dk(ig, ik + 1) & |
|---|
| 610 | + (mlayer(ik) - layer(ik)) * Dk(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1)) |
|---|
| 611 | |
|---|
| 612 | ! saturation_water_ice_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * saturation_water_ice(ig, ik + 1) & |
|---|
| 613 | ! + (mlayer(ik) - layer(ik)) * saturation_water_ice(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1)) |
|---|
| 614 | |
|---|
| 615 | if (close_bottom(ig, ik).or.close_top(ig, ik+1)) then |
|---|
| 616 | saturation_water_ice_inter(ig, ik) = 0.999D0 |
|---|
| 617 | else |
|---|
| 618 | saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik)) ! new diffusion interaction |
|---|
| 619 | endif |
|---|
| 620 | |
|---|
| 621 | saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik)) ! new diffusion interaction |
|---|
| 622 | |
|---|
| 623 | 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)) |
|---|
| 624 | |
|---|
| 625 | if(ads_const_D) D_inter(ig, ik) = default_diffcoeff |
|---|
| 626 | |
|---|
| 627 | ! D_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * D_mid(ig, ik + 1) & ! old implementation with direct interpolation |
|---|
| 628 | ! + (mlayer(ik) - layer(ik)) * D_mid(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1)) |
|---|
| 629 | enddo |
|---|
| 630 | endif |
|---|
| 631 | enddo |
|---|
| 632 | |
|---|
| 633 | ! Computation of new adsorption / desorption constants (Ka, Kd coefficients), and C, E, F coefficients |
|---|
| 634 | ! =================================================================================================== |
|---|
| 635 | |
|---|
| 636 | do ig = 1, ngrid ! loop over all gridpoints |
|---|
| 637 | if(.not.watercaptag(ig)) then ! if there is no polar cap |
|---|
| 638 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 639 | |
|---|
| 640 | if (choice_ads.eq.1) then ! Constant, uniform diffusion coefficient D0 (prescribed) |
|---|
| 641 | |
|---|
| 642 | ! first segment of the bilinear function (before monolayer saturation) |
|---|
| 643 | ! calculate the equilibrium adsorption coefficient |
|---|
| 644 | k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & |
|---|
| 645 | / (dexp(-enthalpy_ads / (8.314D0 * ptsoil(ig, ik))) / tau0) |
|---|
| 646 | |
|---|
| 647 | ! calculate the desorption coefficient |
|---|
| 648 | Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 649 | |
|---|
| 650 | ! calculate the absorption coefficient |
|---|
| 651 | Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 652 | |
|---|
| 653 | |
|---|
| 654 | ! second segment of the bilinear function (after monolayer saturation) ! added 2020 |
|---|
| 655 | ! calculate the equilibrium adsorption coefficient |
|---|
| 656 | k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & |
|---|
| 657 | / (dexp(-enthalpy_ads2 / (8.314D0 * ptsoil(ig, ik))) / tau0) ! added 2020 |
|---|
| 658 | |
|---|
| 659 | ! calculate the desorption coefficient |
|---|
| 660 | Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! added 2020 |
|---|
| 661 | |
|---|
| 662 | ! calculate the absorption coefficient |
|---|
| 663 | Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! added 2020 |
|---|
| 664 | |
|---|
| 665 | ! Kd(ig, ik) = exp(-enthalpy_ads / (8.314 * ptsoil(ig, ik))) |
|---|
| 666 | ! & / tau0 ! * 1.D-18 |
|---|
| 667 | ! Ka(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4. ! * 1.D-18 |
|---|
| 668 | |
|---|
| 669 | ! if ((ngrid.eq.1).and.(ik.eq.18)) then ! if 1D simulation and uppermost layer |
|---|
| 670 | ! print * , 'Ka, Kd, Ka / Kd = ', Ka(ig, ik), Kd(ig, ik), Ka(ig, ik) / Kd(ig, ik) |
|---|
| 671 | ! print * , 'k_ads_eq = ', k_ads_eq(ig, ik) |
|---|
| 672 | ! print * , 'porosity = ', porosity(ig, ik) |
|---|
| 673 | ! endif |
|---|
| 674 | |
|---|
| 675 | elseif (choice_ads.eq.2) then ! modified 2020 (with DeltaQ instead of Q) |
|---|
| 676 | ! Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 * |
|---|
| 677 | ! & (1. / ptsoil(ig, ik) - 1. / Tref)) |
|---|
| 678 | ! Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref * sqrt(ptsoil(ig, ik) / Tref) |
|---|
| 679 | ! & / nref |
|---|
| 680 | |
|---|
| 681 | ! first segment of the bilinear function (before monolayer saturation) |
|---|
| 682 | ! calculate the equilibrium adsorption coefficient |
|---|
| 683 | k_ads_eq(ig, ik) = 8.314D0 * rho_soil(ig, ik) * S * ptsoil(ig,ik)/18.D-3 *Kref * dexp(DeltaQ_ads / 8.314D0 * & |
|---|
| 684 | (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 |
|---|
| 685 | |
|---|
| 686 | ! calculate the desorption coefficient |
|---|
| 687 | Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 688 | |
|---|
| 689 | ! calculate the absorption coefficient |
|---|
| 690 | Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) |
|---|
| 691 | |
|---|
| 692 | ! second segment of the bilinear function (after monolayer saturation) ! added 2020 |
|---|
| 693 | ! calculate the equilibrium adsorption coefficient |
|---|
| 694 | k_ads_eq2(ig, ik) = 8.314D0 * rho_soil(ig, ik) * S * ptsoil(ig,ik)/18.D-3 * Kref2 * dexp(DeltaQ_ads2 / 8.314D0 * & |
|---|
| 695 | (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 |
|---|
| 696 | |
|---|
| 697 | ! calculate the desorption coefficient |
|---|
| 698 | Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 699 | |
|---|
| 700 | ! calculate the absorption coefficient |
|---|
| 701 | Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) |
|---|
| 702 | |
|---|
| 703 | elseif (choice_ads.eq.0) then!No adsorption/desorption |
|---|
| 704 | |
|---|
| 705 | Ka(ig, ik) = 0. |
|---|
| 706 | Kd(ig, ik) = 0. |
|---|
| 707 | Ka2(ig, ik) = 0. |
|---|
| 708 | Kd2(ig, ik) = 0. |
|---|
| 709 | endif |
|---|
| 710 | |
|---|
| 711 | ! calculate the amount of water vapor at monolayer saturation ! modified 2020 |
|---|
| 712 | if(choice_ads.ne.0) nsat(ig, ik) = adswater_sat_mono * Kd(ig, ik) / Ka(ig, ik) |
|---|
| 713 | |
|---|
| 714 | ! calculate the intercept of the second line in the bilinear function ! added 2020; corrected 2024 |
|---|
| 715 | c0(ig, ik) = ( 1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono |
|---|
| 716 | |
|---|
| 717 | ! calculate C, E, and F coefficients for later calculations |
|---|
| 718 | C(ig, ik) = interlayer_dz(ig, ik) / ptimestep |
|---|
| 719 | |
|---|
| 720 | ! first segment of the bilinear function (before monolayer saturation) |
|---|
| 721 | E(ig, ik) = Ka(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik)) |
|---|
| 722 | F(ig, ik) = Kd(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik)) |
|---|
| 723 | |
|---|
| 724 | ! second segment of the bilinear function (after monolayer saturation) ! added 2020 |
|---|
| 725 | E2(ig, ik) = Ka2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik)) |
|---|
| 726 | F2(ig, ik) = Kd2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik)) |
|---|
| 727 | |
|---|
| 728 | ! save the old water concentration |
|---|
| 729 | znsoilprev(ig, ik) = znsoil(ig, ik) |
|---|
| 730 | znsoilprev2(ig, ik) = znsoil(ig, ik) ! needed in "special case" loop because adswprev can be changed before this loop |
|---|
| 731 | adswprev(ig, ik) = adswater(ig, ik) !!!! Besoin de adswprev2 ??? |
|---|
| 732 | iceprev(ig, ik) = ice(ig, ik) ! needed in "special case" loop for the same reason |
|---|
| 733 | |
|---|
| 734 | ! calculate the total ice + water vapor content |
|---|
| 735 | ztot1(ig, ik) = porosity_prev(ig, ik) * znsoil(ig, ik) + ice(ig, ik) |
|---|
| 736 | enddo |
|---|
| 737 | endif |
|---|
| 738 | enddo |
|---|
| 739 | |
|---|
| 740 | ! Computation of delta, alpha and beta coefficients ETC !!!! |
|---|
| 741 | ! =========================================================== |
|---|
| 742 | do ig = 1, ngrid ! loop over all gridpoints |
|---|
| 743 | ! initialize the flux to zero to avoid undefined values |
|---|
| 744 | zdqsdifrego(ig) = 0.D0 |
|---|
| 745 | do ik = 0, nsoil - 1 |
|---|
| 746 | flux(ig, ik) = 0. |
|---|
| 747 | enddo |
|---|
| 748 | |
|---|
| 749 | if(.not.watercaptag(ig)) then ! if there is no polar cap |
|---|
| 750 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 751 | if (ice(ig, ik).eq.0.) then |
|---|
| 752 | ice_index_flag(ik) = .false. |
|---|
| 753 | else |
|---|
| 754 | ice_index_flag(ik) = .true. |
|---|
| 755 | endif |
|---|
| 756 | enddo |
|---|
| 757 | |
|---|
| 758 | |
|---|
| 759 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 760 | |
|---|
| 761 | ! calculate A and B coefficients ! modified 2020 |
|---|
| 762 | if ( .not.over_mono_sat_flag(ig, ik) ) then ! Assume cell below the monolayer saturation |
|---|
| 763 | ! calculate A and B coefficients (first segment of the bilinear function) |
|---|
| 764 | A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 765 | B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) |
|---|
| 766 | else ! Assume cell over the monolayer saturation ! added 2020 |
|---|
| 767 | ! calculate A and B coefficients (second segment of the bilinear function) |
|---|
| 768 | A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 769 | B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & |
|---|
| 770 | + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 |
|---|
| 771 | endif |
|---|
| 772 | |
|---|
| 773 | ! calculate the saturation pressure |
|---|
| 774 | 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) ? |
|---|
| 775 | |
|---|
| 776 | ! calculate the gas number density at saturation pressure |
|---|
| 777 | nsatsoil(ig, ik) = (P_sat_soil(ig, ik) * mw) / (kb * ptsoil(ig, ik)) |
|---|
| 778 | enddo |
|---|
| 779 | |
|---|
| 780 | ! calculate the alpha, beta, and delta coefficients for the surface layer |
|---|
| 781 | delta0(ig) = pa(ig, 1) + pb(ig, 2) * (1.D0 - pd(ig, 2)) + porosity_ice_free(ig, 1) * pb(ig, 1) |
|---|
| 782 | alpha0(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / delta0(ig) |
|---|
| 783 | |
|---|
| 784 | beta0(ig) = ( pb(ig, 2) * pc(ig, 2) + pq(ig, 1, igcm_h2o_vap) * pa(ig, 1) + pqsurf(ig) ) & |
|---|
| 785 | / delta0(ig) |
|---|
| 786 | |
|---|
| 787 | ! set loop flag |
|---|
| 788 | do ik = 1, nsoil |
|---|
| 789 | condensation_flag = .false. |
|---|
| 790 | enddo |
|---|
| 791 | |
|---|
| 792 | sublimation_flag = .true. |
|---|
| 793 | sublim_n = 0 |
|---|
| 794 | do while (sublimation_flag) ! while there is sublimation |
|---|
| 795 | ! print *, "sublimation loop: ", sublim_n |
|---|
| 796 | sublim_n = sublim_n + 1 |
|---|
| 797 | |
|---|
| 798 | if (sublim_n.gt.100) then |
|---|
| 799 | print *, "sublimation not converging in call ", n |
|---|
| 800 | print *, "might as well stop now" |
|---|
| 801 | stop |
|---|
| 802 | endif |
|---|
| 803 | |
|---|
| 804 | ! reset flag |
|---|
| 805 | sublimation_flag = .false. |
|---|
| 806 | |
|---|
| 807 | recompute_all_cells_ads_flag = .true. |
|---|
| 808 | satur_mono_n = 0 |
|---|
| 809 | do while (recompute_all_cells_ads_flag) |
|---|
| 810 | ! print *, satur_mono_n |
|---|
| 811 | satur_mono_n = satur_mono_n + 1 |
|---|
| 812 | |
|---|
| 813 | ! reset loop flag |
|---|
| 814 | recompute_all_cells_ads_flag = .false. |
|---|
| 815 | |
|---|
| 816 | ! calculate the surface air density |
|---|
| 817 | rhosurf(ig) = pplev(ig, 1) / (r * ptsrf(ig)) |
|---|
| 818 | |
|---|
| 819 | ! calculate the coefficients in the upermost layer |
|---|
| 820 | if (exchange(ig)) then ! if there is surface exchange |
|---|
| 821 | |
|---|
| 822 | ! calculate the delta, alpha, and beta coefficients |
|---|
| 823 | zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) & |
|---|
| 824 | + D_inter(ig, 1) / midlayer_dz(ig, 1) & !!! est - ce que le terme pb inclut le bon rho ? |
|---|
| 825 | + porosity_ice_free(ig, 1) * pb(ig, 1) / (rho(ig) * ptimestep) * (1.D0 - alpha0(ig)) |
|---|
| 826 | !pourrait remplacer pb/(ptime*rho(1/2)) par zcdv/ptime |
|---|
| 827 | |
|---|
| 828 | zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1) |
|---|
| 829 | |
|---|
| 830 | zbeta(ig, 1) = ( porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * beta0(ig) & |
|---|
| 831 | + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1) ) & |
|---|
| 832 | / zdelta(ig, 1) |
|---|
| 833 | else |
|---|
| 834 | ! calculate the delta, alpha, and beta coefficients without exchange |
|---|
| 835 | zdelta(ig, 1) = A(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1) & |
|---|
| 836 | + D_mid(ig, 1) / midlayer_dz(ig, 0) & |
|---|
| 837 | + porosity(ig, 1) * C(ig, 1) |
|---|
| 838 | |
|---|
| 839 | zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1. / zdelta(ig, 1) |
|---|
| 840 | |
|---|
| 841 | zbeta(ig, 1) = ( D_mid(ig, 1) / midlayer_dz(ig, 0) * qsat_surf(ig) * rhosurf(ig) & |
|---|
| 842 | + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1) ) & |
|---|
| 843 | / zdelta(ig, 1) |
|---|
| 844 | |
|---|
| 845 | endif |
|---|
| 846 | |
|---|
| 847 | ! check for ice |
|---|
| 848 | if (ice_index_flag(1)) then |
|---|
| 849 | ! set the alpha coefficient to zero |
|---|
| 850 | zalpha(ig, 1) = 0.D0 |
|---|
| 851 | ! set the beta coefficient to the number density at saturation pressure |
|---|
| 852 | zbeta(ig, 1) = nsatsoil(ig, 1) |
|---|
| 853 | endif |
|---|
| 854 | |
|---|
| 855 | do ik = 2, nsoil - 1 ! go through all other layers |
|---|
| 856 | |
|---|
| 857 | ! calculate delta, alpha, and beta coefficients |
|---|
| 858 | |
|---|
| 859 | zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) & |
|---|
| 860 | + D_inter(ig, ik) / midlayer_dz(ig, ik) & |
|---|
| 861 | + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1)) |
|---|
| 862 | |
|---|
| 863 | zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik) |
|---|
| 864 | |
|---|
| 865 | zbeta(ig, ik) = ( D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) & |
|---|
| 866 | + B(ig, ik) & |
|---|
| 867 | + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev(ig, ik) ) / zdelta(ig, ik) |
|---|
| 868 | |
|---|
| 869 | ! check for ice |
|---|
| 870 | if (ice_index_flag(ik)) then |
|---|
| 871 | ! set the alpha coefficient to zero |
|---|
| 872 | zalpha(ig, ik) = 0.D0 |
|---|
| 873 | ! set the beta coefficient to the number density at saturation pressure |
|---|
| 874 | zbeta(ig, ik) = nsatsoil(ig, ik) |
|---|
| 875 | endif |
|---|
| 876 | enddo |
|---|
| 877 | |
|---|
| 878 | ! Computation of pqsoil, ztot1, zq1temp2 and zdqsdifrego |
|---|
| 879 | ! ======================================================= = |
|---|
| 880 | |
|---|
| 881 | ! calculate the preliminary amount of water vapor in the bottom layer |
|---|
| 882 | if (ice_index_flag(nsoil)) then ! if there is ice |
|---|
| 883 | ! set the vapor number density to saturation |
|---|
| 884 | znsoil(ig, nsoil) = nsatsoil(ig, nsoil) |
|---|
| 885 | else |
|---|
| 886 | ! calculate the vapor number density in the lowest layer |
|---|
| 887 | znsoil(ig, nsoil) = ( D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) & |
|---|
| 888 | + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev(ig, nsoil) + B(ig, nsoil) ) & |
|---|
| 889 | / ( porosity(ig, nsoil) * C(ig, nsoil) + A(ig, nsoil) & |
|---|
| 890 | + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)) ) |
|---|
| 891 | endif |
|---|
| 892 | |
|---|
| 893 | ! calculate the preliminary amount of adsorbed water |
|---|
| 894 | if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024 |
|---|
| 895 | adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) & |
|---|
| 896 | / (1.D0 + ptimestep * Kd(ig, nsoil)) |
|---|
| 897 | else |
|---|
| 898 | adswater_temp(ig, nsoil) = ( Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) + ptimestep * c0(ig,nsoil) * Kd2(ig,nsoil)) & |
|---|
| 899 | / (1.D0 + ptimestep * Kd2(ig, nsoil)) |
|---|
| 900 | endif |
|---|
| 901 | |
|---|
| 902 | |
|---|
| 903 | |
|---|
| 904 | do ik = nsoil-1, 1, -1 ! backward loop over all layers |
|---|
| 905 | |
|---|
| 906 | znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik) |
|---|
| 907 | |
|---|
| 908 | if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024 |
|---|
| 909 | adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) & |
|---|
| 910 | / (1.D0 + ptimestep * Kd(ig, ik)) |
|---|
| 911 | else |
|---|
| 912 | adswater_temp(ig, ik) = ( Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) + ptimestep * c0(ig,ik) * Kd2(ig,ik)) & |
|---|
| 913 | / (1.D0 + ptimestep * Kd2(ig, ik)) |
|---|
| 914 | endif |
|---|
| 915 | enddo |
|---|
| 916 | |
|---|
| 917 | ! check if any cell is over monolayer saturation |
|---|
| 918 | do ik = 1, nsoil ! loop over all soil levels |
|---|
| 919 | |
|---|
| 920 | ! 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 |
|---|
| 921 | ! recompute_cell_ads_flag(ik) = .true. |
|---|
| 922 | ! cycle ! Jump loop |
|---|
| 923 | ! endif |
|---|
| 924 | |
|---|
| 925 | ! Change of coefficients ! ARNAU |
|---|
| 926 | recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients |
|---|
| 927 | if ( adswater_temp(ig, ik).ge.adswater_sat_mono ) then |
|---|
| 928 | |
|---|
| 929 | print *, "NOTICE: over monolayer saturation" |
|---|
| 930 | |
|---|
| 931 | if ( .not.over_mono_sat_flag(ig, ik) ) then |
|---|
| 932 | ! If the cell enters this scope, it |
|---|
| 933 | ! means that the cell is over the monolayer |
|---|
| 934 | ! saturation after calculation, but the |
|---|
| 935 | ! coefficients assume it is below. Therefore, |
|---|
| 936 | ! the cell needs to be recomputed |
|---|
| 937 | |
|---|
| 938 | over_mono_sat_flag(ig, ik) = .true. |
|---|
| 939 | recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients |
|---|
| 940 | |
|---|
| 941 | ! 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) |
|---|
| 942 | A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 943 | B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & |
|---|
| 944 | + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 |
|---|
| 945 | endif |
|---|
| 946 | else |
|---|
| 947 | if ( over_mono_sat_flag(ig, ik) ) then |
|---|
| 948 | ! If the cell enters this scope, it |
|---|
| 949 | ! means that the cell is below the monolayer |
|---|
| 950 | ! saturation after calculation, but the |
|---|
| 951 | ! coefficients assume it is above. Therefore, |
|---|
| 952 | ! the cell needs to be recomputed |
|---|
| 953 | |
|---|
| 954 | over_mono_sat_flag(ig, ik) = .false. |
|---|
| 955 | recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients |
|---|
| 956 | |
|---|
| 957 | ! calculate A and B coefficients (reset to first segment of the bilinear function) |
|---|
| 958 | A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 959 | B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) |
|---|
| 960 | endif |
|---|
| 961 | endif |
|---|
| 962 | enddo |
|---|
| 963 | |
|---|
| 964 | ! if one cell needs to be recomputed, then all the column is to be recomputed too |
|---|
| 965 | do ik = 1, nsoil |
|---|
| 966 | if ( recompute_cell_ads_flag(ik) ) then |
|---|
| 967 | recompute_all_cells_ads_flag = .true. |
|---|
| 968 | else |
|---|
| 969 | 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) |
|---|
| 970 | endif |
|---|
| 971 | enddo |
|---|
| 972 | enddo ! end loop while recompute_all_cells_ads_flag (monolayer saturation) |
|---|
| 973 | |
|---|
| 974 | !? I'm not sure if this should be calculated here again. I have a feeling that ztot1 is meant |
|---|
| 975 | ! as the value from the previous timestep |
|---|
| 976 | !do ik = 1, nsoil |
|---|
| 977 | ! ztot1(ig, ik) = porosity_prev(ig, ik) * znsoil(ig, ik) + ice(ig, ik) |
|---|
| 978 | !enddo |
|---|
| 979 | |
|---|
| 980 | ! Calculation of Fnet + Fads |
|---|
| 981 | ! ============================= |
|---|
| 982 | |
|---|
| 983 | ! calulate the flux variable |
|---|
| 984 | |
|---|
| 985 | ! calculate the flux in the top layer |
|---|
| 986 | if (exchange(ig)) then ! if there is exchange |
|---|
| 987 | ! calculate the flux taking diffusion to the atmosphere into account |
|---|
| 988 | flux(ig, 0) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep & |
|---|
| 989 | * ( znsoil(ig, 1) / rho(ig) - beta0(ig) - alpha0(ig) * znsoil(ig, 1) / rho(ig) ) |
|---|
| 990 | elseif (pqsurf(ig).gt.0.) then |
|---|
| 991 | ! assume that the surface is covered by waterice (if it is co2ice it should not call this subroutine at all) |
|---|
| 992 | flux(ig, 0) = D_mid(ig, 1) / midlayer_dz(ig, 0) * ( znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig) ) |
|---|
| 993 | endif |
|---|
| 994 | |
|---|
| 995 | ! check if the ground would take up water from the surface but there is none |
|---|
| 996 | if ((.not.exchange(ig)).and.(pqsurf(ig).eq.0.).and.(flux(ig, 0).lt.0.)) then |
|---|
| 997 | flux(ig, 0) = 0. |
|---|
| 998 | endif |
|---|
| 999 | |
|---|
| 1000 | ! calculate the flux in all other layers |
|---|
| 1001 | do ik = 1, nsoil - 1 |
|---|
| 1002 | flux(ig, ik) = D_inter(ig, ik) * ( znsoil(ig, ik + 1) - znsoil(ig, ik) ) / midlayer_dz(ig, ik) |
|---|
| 1003 | enddo |
|---|
| 1004 | |
|---|
| 1005 | ! calculate dztot1 which descibes the change in ztot1 (water vapor and ice). It is only used for output (directly and indirectly) |
|---|
| 1006 | |
|---|
| 1007 | ! calculate the change in ztot1 |
|---|
| 1008 | |
|---|
| 1009 | do ik = 1, nsoil - 1 |
|---|
| 1010 | dztot1(ik) = ( flux(ig, ik) - flux(ig, ik - 1) ) / interlayer_dz(ig, ik) & |
|---|
| 1011 | - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) & |
|---|
| 1012 | + B(ig, ik) / interlayer_dz(ig, ik) |
|---|
| 1013 | enddo |
|---|
| 1014 | |
|---|
| 1015 | dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) & |
|---|
| 1016 | - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) & |
|---|
| 1017 | + B(ig, nsoil) / interlayer_dz(ig, nsoil) |
|---|
| 1018 | |
|---|
| 1019 | ! Condensation / sublimation |
|---|
| 1020 | do ik = 1, nsoil ! loop over all layers |
|---|
| 1021 | ! print *, "ice in layer ", ik, ": ", ice(ig, ik) |
|---|
| 1022 | ! print *, "vapor in layer ", ik, ": ", znsoil(ig, ik) |
|---|
| 1023 | ! print *, "sat dens in layer ", ik, ": ", nsatsoil(ig, ik) |
|---|
| 1024 | ! print *, "" |
|---|
| 1025 | if (ice_index_flag(ik)) then ! if there is ice |
|---|
| 1026 | |
|---|
| 1027 | ! Compute ice content |
|---|
| 1028 | ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik) |
|---|
| 1029 | |
|---|
| 1030 | if (ice(ig, ik).lt.0.D0) then ! if all the ice is used up |
|---|
| 1031 | ! print *, "########complete sublimation in layer", ik, " cell:", ig |
|---|
| 1032 | ! print *, "ztot1: ", ztot1(ig, ik) |
|---|
| 1033 | ! print *, "dztot1*timestep: ", dztot1(ik) * ptimestep |
|---|
| 1034 | ! print *, "vapour: ", porosity(ig, ik) * nsatsoil(ig, ik) |
|---|
| 1035 | ! print *, "znsoil: ", znsoil(ig, ik) |
|---|
| 1036 | ! print *, "nsatsoil: ", nsatsoil(ig, ik) |
|---|
| 1037 | ! print *, "porosity: ", porosity(ig, ik) |
|---|
| 1038 | ! print *, "ice: ", ice(ig, ik) |
|---|
| 1039 | ! print *, "exchange: ", exchange(ig) |
|---|
| 1040 | ! print *, "co2ice: ", co2ice(ig) |
|---|
| 1041 | ! print *, "" |
|---|
| 1042 | ! print *, "zalpha: ", zalpha(ig, ik) |
|---|
| 1043 | ! print *, "zbeta: ", zbeta(ig, ik) |
|---|
| 1044 | ! print *, "" |
|---|
| 1045 | |
|---|
| 1046 | ! set the ice content to zero |
|---|
| 1047 | ice(ig, ik) = 0.D0 |
|---|
| 1048 | |
|---|
| 1049 | ! change the water vapor from the previous timestep. (Watch out! could go wrong) |
|---|
| 1050 | znsoilprev(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik) |
|---|
| 1051 | ! print *, "new znsoilprev", znsoilprev(ig, ik) |
|---|
| 1052 | |
|---|
| 1053 | ! remove the ice flag and raise the sublimation flag |
|---|
| 1054 | ice_index_flag(ik) = .false. |
|---|
| 1055 | sublimation_flag = .true. |
|---|
| 1056 | endif |
|---|
| 1057 | |
|---|
| 1058 | elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then ! if there is condenstation |
|---|
| 1059 | !ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik) |
|---|
| 1060 | |
|---|
| 1061 | |
|---|
| 1062 | ! print *, "=========== new condensation in layer", ik, " cell:", ig |
|---|
| 1063 | ! print *, "znsoil, nsatsoil: ", znsoil(ig, ik), nsatsoil(ig, ik) |
|---|
| 1064 | ! print *, "ztot1: ", ztot1(ig, ik) |
|---|
| 1065 | ! print *, "dztot1: ", dztot1(ik) |
|---|
| 1066 | ! print *, "ice: ", ice(ig,ik) |
|---|
| 1067 | ! print *, "" |
|---|
| 1068 | ! print *, "zalpha: ", zalpha(ig, ik) |
|---|
| 1069 | ! print *, "zbeta: ", zbeta(ig, ik) |
|---|
| 1070 | ! print *, "" |
|---|
| 1071 | |
|---|
| 1072 | !if (ice(ig, ik).lt.0.D0) then |
|---|
| 1073 | ! set the ice content to zero |
|---|
| 1074 | ! ice(ig, ik) = 0.D0 |
|---|
| 1075 | if (.not.condensation_flag(ik)) then |
|---|
| 1076 | ! raise the ice and sublimation flags |
|---|
| 1077 | ice_index_flag(ik) = .true. |
|---|
| 1078 | sublimation_flag = .true. |
|---|
| 1079 | ! print *, condensation_flag(ik) |
|---|
| 1080 | condensation_flag(ik) = .true. |
|---|
| 1081 | ! print *, condensation_flag(ik) |
|---|
| 1082 | ! print *, "set condensation flag" |
|---|
| 1083 | ! else |
|---|
| 1084 | ! print *, "condensation loop supressed" |
|---|
| 1085 | endif |
|---|
| 1086 | endif |
|---|
| 1087 | |
|---|
| 1088 | ! calculate the difference between total water content and ice + vapor content (only used for output) |
|---|
| 1089 | diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) & |
|---|
| 1090 | + adswater(ig, ik) - ztot1(ig, ik) - dztot1(ik) * ptimestep |
|---|
| 1091 | enddo ! loop over all layer |
|---|
| 1092 | enddo ! end first loop while sublimation_flag (condensation / sublimation) |
|---|
| 1093 | |
|---|
| 1094 | if (exchange(ig)) then ! if there is exchange |
|---|
| 1095 | |
|---|
| 1096 | ! calculate the temporty mixing ratio above the surface |
|---|
| 1097 | zq1temp2(ig) = beta0(ig) + alpha0(ig) * znsoil(ig, 1) / rho(ig) |
|---|
| 1098 | ! calculate the flux from the subsurface |
|---|
| 1099 | zdqsdifrego(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * (zq1temp2(ig) - znsoil(ig, 1) / rho(ig) ) |
|---|
| 1100 | |
|---|
| 1101 | else |
|---|
| 1102 | ! calculate the flux from the subsurface |
|---|
| 1103 | zdqsdifrego(ig) = D_mid(ig, 1) / midlayer_dz(ig, 0) * (znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig)) ! faut - il faire intervenir la porosite ? |
|---|
| 1104 | if ((pqsurf(ig).eq.0.).and.(zdqsdifrego(ig).lt.0.)) then |
|---|
| 1105 | zdqsdifrego(ig) = 0. |
|---|
| 1106 | endif |
|---|
| 1107 | endif |
|---|
| 1108 | |
|---|
| 1109 | ! Special case where there is not enough ice on the surface to supply the subsurface for the whole timestep |
|---|
| 1110 | ! (when exchange with the atmosphere is disabled): the upper boundary condition becomes a flux condition |
|---|
| 1111 | ! (and not qsat_surf) and all the remaining surface ice is sublimated and transferred to the subsurface. |
|---|
| 1112 | ! the actual change in surface ice happens in a higher subroutine |
|---|
| 1113 | ! ========================================================================================================= = |
|---|
| 1114 | |
|---|
| 1115 | if ( (.not.exchange(ig)) & |
|---|
| 1116 | .and. ( (-zdqsdifrego(ig) * ptimestep) & |
|---|
| 1117 | .gt.( pqsurf(ig) + pdqsdifpot(ig) * ptimestep) ) & |
|---|
| 1118 | .and.( (pqsurf(ig) + pdqsdifpot(ig) * ptimestep).gt.0. ) ) then |
|---|
| 1119 | |
|---|
| 1120 | ! calculate a new flux from the subsurface |
|---|
| 1121 | zdqsdifrego(ig) = -( pqsurf(ig) + pdqsdifpot(ig) * ptimestep ) / ptimestep |
|---|
| 1122 | |
|---|
| 1123 | ! ! check case where there is CO2 ice on the surface but qsurf = 0 |
|---|
| 1124 | ! if (co2ice(ig).gt.0.) then |
|---|
| 1125 | ! zdqsdifrego(ig) = 0.D0 |
|---|
| 1126 | ! endif |
|---|
| 1127 | do ik = 1, nsoil |
|---|
| 1128 | |
|---|
| 1129 | ! calculate A and B coefficients ! modified 2020 |
|---|
| 1130 | if ( .not.over_mono_sat_flag(ig, ik) ) then ! Assume cell below the monolayer saturation |
|---|
| 1131 | ! calculate A and B coefficients (first segment of the bilinear function) |
|---|
| 1132 | A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1133 | B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) |
|---|
| 1134 | else ! Assume cell over the monolayer saturation ! added 2020 |
|---|
| 1135 | ! calculate A and B coefficients (second segment of the bilinear function) |
|---|
| 1136 | A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1137 | B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & |
|---|
| 1138 | + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 |
|---|
| 1139 | endif |
|---|
| 1140 | |
|---|
| 1141 | ! initialize all flags for the loop |
|---|
| 1142 | |
|---|
| 1143 | |
|---|
| 1144 | ! initialize the ice |
|---|
| 1145 | ice(ig, ik) = iceprev(ig, ik) |
|---|
| 1146 | if (iceprev(ig, ik).eq.0.) then |
|---|
| 1147 | ice_index_flag(ik) = .false. |
|---|
| 1148 | else |
|---|
| 1149 | ice_index_flag(ik) = .true. |
|---|
| 1150 | endif |
|---|
| 1151 | enddo |
|---|
| 1152 | |
|---|
| 1153 | ! loop while there is new sublimation |
|---|
| 1154 | |
|---|
| 1155 | sublimation_flag = .true. |
|---|
| 1156 | sublim_n = 0 |
|---|
| 1157 | do while (sublimation_flag) |
|---|
| 1158 | ! print *, "special case sublimation loop: ", sublim_n |
|---|
| 1159 | sublim_n = sublim_n + 1 |
|---|
| 1160 | if (sublim_n.gt.100) then |
|---|
| 1161 | print *, "special case sublimation not converging in call ", n |
|---|
| 1162 | print *, "might as well stop now" |
|---|
| 1163 | stop |
|---|
| 1164 | endif |
|---|
| 1165 | |
|---|
| 1166 | ! reset the sublimation flag |
|---|
| 1167 | sublimation_flag = .false. |
|---|
| 1168 | |
|---|
| 1169 | ! loop until good values accounting for monolayer saturation are found |
|---|
| 1170 | recompute_all_cells_ads_flag = .true. |
|---|
| 1171 | do while (recompute_all_cells_ads_flag) |
|---|
| 1172 | |
|---|
| 1173 | ! reset loop flag |
|---|
| 1174 | recompute_all_cells_ads_flag = .false. |
|---|
| 1175 | |
|---|
| 1176 | ! calculate the Delta, Alpha, and Beta coefficients in the top layer |
|---|
| 1177 | zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1) |
|---|
| 1178 | |
|---|
| 1179 | zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1) |
|---|
| 1180 | |
|---|
| 1181 | zbeta(ig, 1) = ( porosity_prev(ig, 1) * C(ig, 1) * znsoilprev2(ig, 1) + B(ig, 1) - zdqsdifrego(ig) ) & |
|---|
| 1182 | / zdelta(ig, 1) |
|---|
| 1183 | |
|---|
| 1184 | ! check for ice |
|---|
| 1185 | if (ice_index_flag(1)) then |
|---|
| 1186 | ! set the alpha coefficient to zero |
|---|
| 1187 | zalpha(ig, 1) = 0.D0 |
|---|
| 1188 | ! set the beta coefficient to the number density at saturation pressure |
|---|
| 1189 | zbeta(ig, 1) = nsatsoil(ig, 1) |
|---|
| 1190 | endif |
|---|
| 1191 | |
|---|
| 1192 | do ik = 2, nsoil - 1 ! loop over all other layers |
|---|
| 1193 | |
|---|
| 1194 | ! calculate the Delta, Alpha, and Beta coefficients in the layer |
|---|
| 1195 | zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) + D_inter(ig, ik) / midlayer_dz(ig, ik) & |
|---|
| 1196 | + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1)) |
|---|
| 1197 | |
|---|
| 1198 | zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik) |
|---|
| 1199 | |
|---|
| 1200 | zbeta(ig, ik) = ( D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) & |
|---|
| 1201 | + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev2(ig, ik) + B(ig, ik) ) / zdelta(ig, ik) |
|---|
| 1202 | |
|---|
| 1203 | ! check for ice |
|---|
| 1204 | if (ice_index_flag(ik)) then |
|---|
| 1205 | ! set the alpha coefficient to zero |
|---|
| 1206 | zalpha(ig, ik) = 0.D0 |
|---|
| 1207 | ! set the beta coefficient to the number density at saturation pressure |
|---|
| 1208 | zbeta(ig, ik) = nsatsoil(ig, ik) |
|---|
| 1209 | endif |
|---|
| 1210 | enddo |
|---|
| 1211 | |
|---|
| 1212 | ! calculate the preliminary amount of water vapor in the bottom layer |
|---|
| 1213 | if (ice_index_flag(nsoil)) then |
|---|
| 1214 | ! set the vapor number density to saturation |
|---|
| 1215 | znsoil(ig, nsoil) = nsatsoil(ig, nsoil) |
|---|
| 1216 | else |
|---|
| 1217 | ! calculate the vapor number density in the lowest layer |
|---|
| 1218 | znsoil(ig, nsoil) = ( D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) & |
|---|
| 1219 | + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev2(ig, nsoil) + B(ig, nsoil) ) & |
|---|
| 1220 | / ( porosity(ig, nsoil) * C(ig, nsoil) & |
|---|
| 1221 | + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)) & |
|---|
| 1222 | + A(ig, nsoil) ) |
|---|
| 1223 | endif |
|---|
| 1224 | |
|---|
| 1225 | |
|---|
| 1226 | ! calculate the preliminary amount of adsorbed water |
|---|
| 1227 | if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024 |
|---|
| 1228 | adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) & |
|---|
| 1229 | / (1.D0 + ptimestep * Kd(ig, nsoil)) |
|---|
| 1230 | else |
|---|
| 1231 | adswater_temp(ig, nsoil) = ( Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) + ptimestep * c0(ig,nsoil) * Kd2(ig,nsoil)) & |
|---|
| 1232 | / (1.D0 + ptimestep * Kd2(ig, nsoil)) |
|---|
| 1233 | endif |
|---|
| 1234 | |
|---|
| 1235 | |
|---|
| 1236 | |
|---|
| 1237 | do ik = nsoil-1, 1, -1 ! backward loop over all layers |
|---|
| 1238 | |
|---|
| 1239 | znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik) |
|---|
| 1240 | |
|---|
| 1241 | if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024 |
|---|
| 1242 | adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) & |
|---|
| 1243 | / (1.D0 + ptimestep * Kd(ig, ik)) |
|---|
| 1244 | else |
|---|
| 1245 | adswater_temp(ig, ik) = ( Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) + ptimestep * c0(ig,ik) * Kd2(ig,ik)) & |
|---|
| 1246 | / (1.D0 + ptimestep * Kd2(ig, ik)) |
|---|
| 1247 | endif |
|---|
| 1248 | enddo |
|---|
| 1249 | |
|---|
| 1250 | |
|---|
| 1251 | ! check for any saturation |
|---|
| 1252 | do ik = 1, nsoil |
|---|
| 1253 | ! Change of coefficients ! ARNAU |
|---|
| 1254 | recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients ! ARNAU |
|---|
| 1255 | if ( adswater_temp(ig, ik).gt.adswater_sat_mono ) then |
|---|
| 1256 | |
|---|
| 1257 | print *, "NOTICE: over monolayer saturation" |
|---|
| 1258 | |
|---|
| 1259 | if ( .not.over_mono_sat_flag(ig, ik) ) then |
|---|
| 1260 | ! If the cell enters this scope, it |
|---|
| 1261 | ! means that the cell is over the monolayer |
|---|
| 1262 | ! saturation after calculation, but the |
|---|
| 1263 | ! coefficients assume it is below. Therefore, |
|---|
| 1264 | ! the cell needs to be recomputed |
|---|
| 1265 | |
|---|
| 1266 | over_mono_sat_flag(ig, ik) = .true. |
|---|
| 1267 | recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients |
|---|
| 1268 | |
|---|
| 1269 | ! 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) |
|---|
| 1270 | A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1271 | B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & |
|---|
| 1272 | + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 |
|---|
| 1273 | endif |
|---|
| 1274 | else |
|---|
| 1275 | if ( over_mono_sat_flag(ig, ik) ) then |
|---|
| 1276 | ! If the cell enters this scope, it |
|---|
| 1277 | ! means that the cell is below the monolayer |
|---|
| 1278 | ! saturation after calculation, but the |
|---|
| 1279 | ! coefficients assume it is above. Therefore, |
|---|
| 1280 | ! the cell needs to be recomputed |
|---|
| 1281 | |
|---|
| 1282 | over_mono_sat_flag(ig, ik) = .false. |
|---|
| 1283 | recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients |
|---|
| 1284 | |
|---|
| 1285 | ! calculate A and B coefficients (reset to first segment of the bilinear function) |
|---|
| 1286 | A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1287 | B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) |
|---|
| 1288 | endif |
|---|
| 1289 | endif |
|---|
| 1290 | enddo |
|---|
| 1291 | |
|---|
| 1292 | ! raise the saturation flag if any layer has saturated and reset the first layer saturation flag |
|---|
| 1293 | do ik = 1, nsoil ! modified 2020 |
|---|
| 1294 | if ( recompute_cell_ads_flag(ik) ) then |
|---|
| 1295 | recompute_all_cells_ads_flag = .true. |
|---|
| 1296 | else |
|---|
| 1297 | 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) |
|---|
| 1298 | endif |
|---|
| 1299 | enddo |
|---|
| 1300 | enddo ! end while loop (adsorption saturation) |
|---|
| 1301 | |
|---|
| 1302 | ! set the flux to the previously calculated value for the top layer |
|---|
| 1303 | flux(ig, 0) = zdqsdifrego(ig) |
|---|
| 1304 | |
|---|
| 1305 | ! calculate the flux for all other layers |
|---|
| 1306 | do ik = 1, nsoil - 1 |
|---|
| 1307 | flux(ig, ik) = D_inter(ig, ik) * (znsoil(ig, ik + 1) - znsoil(ig, ik)) / midlayer_dz(ig, ik) |
|---|
| 1308 | enddo |
|---|
| 1309 | |
|---|
| 1310 | ! calculate the change in ztot1 |
|---|
| 1311 | do ik = 1, nsoil - 1 |
|---|
| 1312 | dztot1(ik) = (flux(ig, ik) - flux(ig, ik - 1)) / interlayer_dz(ig, ik) & |
|---|
| 1313 | - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) & |
|---|
| 1314 | + B(ig, ik) / interlayer_dz(ig, ik) |
|---|
| 1315 | enddo |
|---|
| 1316 | |
|---|
| 1317 | dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) & |
|---|
| 1318 | - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) & |
|---|
| 1319 | + B(ig, nsoil) / interlayer_dz(ig, nsoil) |
|---|
| 1320 | |
|---|
| 1321 | |
|---|
| 1322 | ! Condensation / sublimation |
|---|
| 1323 | do ik = 1, nsoil |
|---|
| 1324 | if (ice_index_flag(ik)) then |
|---|
| 1325 | |
|---|
| 1326 | ! Compute ice content |
|---|
| 1327 | ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep & |
|---|
| 1328 | - porosity(ig, ik) * nsatsoil(ig, ik) |
|---|
| 1329 | |
|---|
| 1330 | if (ice(ig, ik).lt.0.D0) then ! if all the ice is used up |
|---|
| 1331 | |
|---|
| 1332 | ! print *, "############complete sublimation in layer ", ik |
|---|
| 1333 | ! print *, "ztot1: ", ztot1(ig, ik) |
|---|
| 1334 | ! print *, "dztot1*timestep: ", dztot1(ik) * ptimestep |
|---|
| 1335 | ! print *, "vapour: ", porosity(ig, ik) * nsatsoil(ig, ik) |
|---|
| 1336 | ! print *, "znsoil: ", znsoil(ig, ik) |
|---|
| 1337 | ! print *, "nsatsoil: ", nsatsoil(ig, ik) |
|---|
| 1338 | ! print *, "porosity: ", porosity(ig, ik) |
|---|
| 1339 | ! print *, "ice: ", ice(ig, ik) |
|---|
| 1340 | ! print *, "exchange: ", exchange(ig) |
|---|
| 1341 | ! print *, "co2ice: ", co2ice(ig) |
|---|
| 1342 | ! print *, "" |
|---|
| 1343 | ! print *, "zalpha: ", zalpha(ig, ik) |
|---|
| 1344 | ! print *, "zbeta: ", zbeta(ig, ik) |
|---|
| 1345 | ! print *, "" |
|---|
| 1346 | |
|---|
| 1347 | ! set the ice content to zero |
|---|
| 1348 | ice(ig, ik) = 0.D0 |
|---|
| 1349 | |
|---|
| 1350 | |
|---|
| 1351 | if (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then |
|---|
| 1352 | print *, "WARNING! complete sublimation, but vapor is oversaturated" |
|---|
| 1353 | print *, "special case loop in cell", ig, ik ,"will be supressed" |
|---|
| 1354 | else |
|---|
| 1355 | ! change the water vapor from the previous timestep. (Watch out! could go wrong) |
|---|
| 1356 | znsoilprev2(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik) |
|---|
| 1357 | ! print *, "new znsoilprev", znsoilprev(ig, ik) |
|---|
| 1358 | |
|---|
| 1359 | ! remove the ice flag and raise the sublimation flag |
|---|
| 1360 | ice_index_flag(ik) = .false. |
|---|
| 1361 | sublimation_flag = .true. |
|---|
| 1362 | endif |
|---|
| 1363 | endif |
|---|
| 1364 | |
|---|
| 1365 | elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then ! if there is new ice through condensation |
|---|
| 1366 | if (.not.condensation_flag(ik)) then |
|---|
| 1367 | ! raise the ice and sublimation flags |
|---|
| 1368 | ice_index_flag(ik) = .true. |
|---|
| 1369 | sublimation_flag = .true. |
|---|
| 1370 | ! print *, condensation_flag(ik) |
|---|
| 1371 | condensation_flag(ik) = .true. |
|---|
| 1372 | ! print *, condensation_flag(ik) |
|---|
| 1373 | ! print *, "set condensation flag" |
|---|
| 1374 | ! else |
|---|
| 1375 | ! print *, "special case condensation loop supressed" |
|---|
| 1376 | endif |
|---|
| 1377 | endif |
|---|
| 1378 | enddo |
|---|
| 1379 | enddo ! end first loop while (condensation / sublimation) |
|---|
| 1380 | endif ! Special case for all ice on the surface is used up |
|---|
| 1381 | |
|---|
| 1382 | do ik = 1, nsoil |
|---|
| 1383 | |
|---|
| 1384 | diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) & ! only used for output |
|---|
| 1385 | + adswater(ig, ik) - adswprev(ig, ik) & |
|---|
| 1386 | - ztot1(ig, ik) - dztot1(ik) * ptimestep !? This is inconsistent, in the not special case the previous adsorbed water is not counted |
|---|
| 1387 | !? also this just overwrites the previous calculation if I see that correctly |
|---|
| 1388 | |
|---|
| 1389 | ! calculate the total amount of water |
|---|
| 1390 | ztot1(ig, ik) = porosity(ig, ik) * znsoil(ig, ik) + ice(ig, ik) |
|---|
| 1391 | h2otot(ig, ik) = adswater(ig, ik) + ztot1(ig, ik) |
|---|
| 1392 | enddo |
|---|
| 1393 | endif ! if there is no polar cap |
|---|
| 1394 | enddo ! loop over all gridpoints |
|---|
| 1395 | |
|---|
| 1396 | ! check for choking condition |
|---|
| 1397 | do ig = 1, ngrid |
|---|
| 1398 | if(.not.watercaptag(ig)) then |
|---|
| 1399 | do ik = 1, nsoil - 1 |
|---|
| 1400 | 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 |
|---|
| 1401 | if ( flux(ig, ik - 1).gt.0.) then |
|---|
| 1402 | if (.not.close_top(ig, ik).and.print_closure_warnings) then |
|---|
| 1403 | print *, "NOTICE: closing top of soil layer due to inward flux", ig, ik |
|---|
| 1404 | endif |
|---|
| 1405 | close_top(ig, ik) = .true. |
|---|
| 1406 | elseif (flux(ig, ik - 1).lt.0.) then |
|---|
| 1407 | if (close_top(ig, ik).and.print_closure_warnings) then |
|---|
| 1408 | print *, "NOTICE: opening top of soil layer due to outward flux", ig, ik |
|---|
| 1409 | endif |
|---|
| 1410 | close_top(ig, ik) = .false. |
|---|
| 1411 | endif |
|---|
| 1412 | |
|---|
| 1413 | if ( flux(ig, ik).lt.0.) then |
|---|
| 1414 | if (.not.close_bottom(ig, ik).and.print_closure_warnings) then |
|---|
| 1415 | print *, "NOTICE: closing bottom of soil layer due to inward flux", ig, ik |
|---|
| 1416 | endif |
|---|
| 1417 | close_bottom(ig, ik) = .true. |
|---|
| 1418 | elseif (flux(ig, ik).gt.0.) then |
|---|
| 1419 | if (close_bottom(ig, ik).and.print_closure_warnings) then |
|---|
| 1420 | print *, "NOTICE: opening bottom of soil layer due to outward flux", ig, ik |
|---|
| 1421 | endif |
|---|
| 1422 | close_bottom(ig, ik) = .false. |
|---|
| 1423 | endif |
|---|
| 1424 | |
|---|
| 1425 | ! if ( (flux(ig, ik).lt.flux(ig, ik - 1)).and.(flux(ig, ik).gt.0D0) ) then |
|---|
| 1426 | ! if (.not.close_top(ig, ik).and.print_closure_warnings) then |
|---|
| 1427 | ! print *, "NOTICE: closing top of soil layer due to ice", ig, ik |
|---|
| 1428 | ! endif |
|---|
| 1429 | ! close_top(ig, ik) = .true. |
|---|
| 1430 | ! |
|---|
| 1431 | ! elseif ( (flux(ig, ik).lt.flux(ig, ik - 1)).and.(flux(ig, ik - 1).lt.0D0) ) then |
|---|
| 1432 | ! if (.not.close_bottom(ig, ik).and.print_closure_warnings) then |
|---|
| 1433 | ! print *, "NOTICE: closing bottom of soil layer due to ice", ig, ik |
|---|
| 1434 | ! endif |
|---|
| 1435 | ! close_bottom(ig, ik) = .true. |
|---|
| 1436 | ! endif |
|---|
| 1437 | ! |
|---|
| 1438 | ! if (close_top(ig, ik).and.close_bottom(ig, ik)) then |
|---|
| 1439 | ! close_top(ig, ik) = .false. |
|---|
| 1440 | ! close_bottom(ig, ik) = .false. |
|---|
| 1441 | ! if (print_closure_warnings) then |
|---|
| 1442 | ! print *, "WARNING: Reopening soil layer after complete closure:", ig, ik |
|---|
| 1443 | ! endif |
|---|
| 1444 | ! elseif (close_top(ig, ik).or.close_bottom(ig, ik)) then |
|---|
| 1445 | ! if ((flux(ig, ik) - flux(ig, ik - 1)).gt.0D0) then |
|---|
| 1446 | ! close_top(ig, ik) = .false. |
|---|
| 1447 | ! close_bottom(ig, ik) = .false. |
|---|
| 1448 | ! if (print_closure_warnings) then |
|---|
| 1449 | ! print *, "WARNING: Reopening soil layer due to rising ice:", ig, ik |
|---|
| 1450 | ! endif |
|---|
| 1451 | ! endif |
|---|
| 1452 | ! |
|---|
| 1453 | ! if (close_top(ig, ik).and.(flux(ig, ik - 1).lt.0D0)) then |
|---|
| 1454 | ! close_top = .false. |
|---|
| 1455 | ! if (print_closure_warnings) then |
|---|
| 1456 | ! print *, "NOTICE: Reopening top of soil layer due to upward flux:", ig, ik |
|---|
| 1457 | ! endif |
|---|
| 1458 | ! endif |
|---|
| 1459 | ! |
|---|
| 1460 | ! if (close_bottom(ig, ik).and.(flux(ig, ik).gt.0D0)) then |
|---|
| 1461 | ! close_bottom = .false. |
|---|
| 1462 | ! if (print_closure_warnings) then |
|---|
| 1463 | ! print *, "NOTICE: Reopening bottom of soil layer due to downward flux:", ig, ik |
|---|
| 1464 | ! endif |
|---|
| 1465 | ! endif |
|---|
| 1466 | ! endif |
|---|
| 1467 | else ! opening condition that ice content has fallen sufficiently |
|---|
| 1468 | if (close_top(ig, ik).or.close_bottom(ig, ik).and.print_closure_warnings) then |
|---|
| 1469 | print *, "NOTICE: Reopening soillayer due to decrease in ice", ig, ik |
|---|
| 1470 | endif |
|---|
| 1471 | close_top(ig, ik) = .false. |
|---|
| 1472 | close_bottom(ig, ik) = .false. |
|---|
| 1473 | endif |
|---|
| 1474 | enddo |
|---|
| 1475 | endif |
|---|
| 1476 | enddo |
|---|
| 1477 | |
|---|
| 1478 | ! Computation of total mass |
|---|
| 1479 | ! =========================== |
|---|
| 1480 | |
|---|
| 1481 | do ig = 1, ngrid |
|---|
| 1482 | ! initialize the masses to zero |
|---|
| 1483 | mass_h2o(ig) = 0.D0 |
|---|
| 1484 | mass_ice(ig) = 0.D0 |
|---|
| 1485 | |
|---|
| 1486 | if(.not.watercaptag(ig)) then |
|---|
| 1487 | do ik = 1, nsoil |
|---|
| 1488 | ! calculate the total water and ice mass |
|---|
| 1489 | mass_h2o(ig) = mass_h2o(ig) + h2otot(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1490 | mass_ice(ig) = mass_ice(ig) + ice(ig, ik) * interlayer_dz(ig, ik) |
|---|
| 1491 | |
|---|
| 1492 | ! calculate how close the water vapor content is to saturizing the adsorbed water |
|---|
| 1493 | if(choice_ads.ne.0) preduite(ig, ik) = znsoil(ig, ik) / nsat(ig, ik) |
|---|
| 1494 | |
|---|
| 1495 | ! write the results to the return variable |
|---|
| 1496 | pqsoil(ig, ik, igcm_h2o_vap_soil) = znsoil(ig, ik) |
|---|
| 1497 | pqsoil(ig, ik, igcm_h2o_ice_soil) = ice(ig, ik) |
|---|
| 1498 | pqsoil(ig, ik, igcm_h2o_vap_ads) = adswater(ig, ik) |
|---|
| 1499 | |
|---|
| 1500 | |
|---|
| 1501 | if (close_top(ig, ik)) then |
|---|
| 1502 | close_out(ig, ik) = 1 |
|---|
| 1503 | elseif (close_bottom(ig, ik)) then |
|---|
| 1504 | close_out(ig, ik) = -1 |
|---|
| 1505 | else |
|---|
| 1506 | close_out(ig, ik) = 0 |
|---|
| 1507 | endif |
|---|
| 1508 | enddo |
|---|
| 1509 | endif |
|---|
| 1510 | |
|---|
| 1511 | if (watercaptag(ig)) then |
|---|
| 1512 | do ik = 1, nsoil |
|---|
| 1513 | saturation_water_ice(ig, ik) = -1 |
|---|
| 1514 | enddo |
|---|
| 1515 | endif |
|---|
| 1516 | |
|---|
| 1517 | ! convert the logical :: flag exchange to a numeric output |
|---|
| 1518 | if (watercaptag(ig)) then |
|---|
| 1519 | exch(ig) = -2. |
|---|
| 1520 | elseif (exchange(ig)) then |
|---|
| 1521 | exch(ig) = 1. |
|---|
| 1522 | else |
|---|
| 1523 | exch(ig) = -1. |
|---|
| 1524 | endif |
|---|
| 1525 | |
|---|
| 1526 | enddo |
|---|
| 1527 | |
|---|
| 1528 | ! calculate the global total value |
|---|
| 1529 | mass_h2o_tot = 0.D0 |
|---|
| 1530 | mass_ice_tot = 0.D0 |
|---|
| 1531 | do ig = 1, ngrid |
|---|
| 1532 | mass_h2o_tot = mass_h2o_tot + mass_h2o(ig) * cell_area(ig) |
|---|
| 1533 | mass_ice_tot = mass_ice_tot + mass_ice(ig) * cell_area(ig) |
|---|
| 1534 | enddo |
|---|
| 1535 | |
|---|
| 1536 | ! print and iterate the call number |
|---|
| 1537 | ! print * , 'Subsurface call n = ', n |
|---|
| 1538 | n = n + 1 |
|---|
| 1539 | |
|---|
| 1540 | ! ----------------------------------------------------------------------- |
|---|
| 1541 | ! Output in diagfi and diagsoil |
|---|
| 1542 | ! ----------------------------------------------------------------------- |
|---|
| 1543 | |
|---|
| 1544 | ! reformat flux, because it has a unusual shape |
|---|
| 1545 | do ig = 1, ngrid |
|---|
| 1546 | nsurf(ig) = rhosurf(ig) * pq(ig, 1, igcm_h2o_vap) |
|---|
| 1547 | |
|---|
| 1548 | |
|---|
| 1549 | do ik = 1, nsoil - 1 |
|---|
| 1550 | var_flux_soil(ig, ik) = flux(ig, ik) |
|---|
| 1551 | enddo |
|---|
| 1552 | var_flux_soil(ig, nsoil) = 0. |
|---|
| 1553 | enddo |
|---|
| 1554 | |
|---|
| 1555 | if(writeoutput) then |
|---|
| 1556 | !print *, "flux shape: ", shape(flux) |
|---|
| 1557 | !print *, "adswater shape ", shape(adswater) |
|---|
| 1558 | |
|---|
| 1559 | !print *, "ngrid= ", ngrid |
|---|
| 1560 | |
|---|
| 1561 | !if (ngrid.eq.1) then ! if in 1D |
|---|
| 1562 | |
|---|
| 1563 | ! print *, "writing 1D data" |
|---|
| 1564 | |
|---|
| 1565 | ! print *, real(adswater(:, :), 4) |
|---|
| 1566 | ! print *, shape(adswater(:, :)) |
|---|
| 1567 | |
|---|
| 1568 | zalpha(1, nsoil) = 0 |
|---|
| 1569 | zbeta(1, nsoil) = 0 |
|---|
| 1570 | |
|---|
| 1571 | ! call write_output("zalpha","diffusion alpha coefficient", "unknown",real(zalpha(1, :))) |
|---|
| 1572 | |
|---|
| 1573 | ! call write_output("zbeta","diffusion beta coefficient", "unknown",real(zbeta(1, :))) |
|---|
| 1574 | |
|---|
| 1575 | ! call write_output("adswater","adswater", "kg / m^3",real(adswater(1, :))) |
|---|
| 1576 | |
|---|
| 1577 | ! call write_output("znsoil","znsoil", "kg / m^3",real(znsoil(1, :))) |
|---|
| 1578 | |
|---|
| 1579 | ! call write_output("ice","ice", "kg / m^3",real(ice(1, :))) |
|---|
| 1580 | |
|---|
| 1581 | ! call write_output("h2otot","total water content", "kg / m^3",real(h2otot(1, :))) |
|---|
| 1582 | |
|---|
| 1583 | ! call write_output("flux_soil","flux_soil", "kg / m^2",real(flux(1, :))) |
|---|
| 1584 | |
|---|
| 1585 | ! call write_output("flux","flux soil", "kg / m^2",real(zdqsdifrego(:))) |
|---|
| 1586 | |
|---|
| 1587 | ! call write_output("flux_surf","surface flux", "kg / m^2",real(zdqsdifrego(1))) |
|---|
| 1588 | |
|---|
| 1589 | ! call write_output("exchange","exchange", "boolean",real(exch(1))) |
|---|
| 1590 | |
|---|
| 1591 | !else |
|---|
| 1592 | |
|---|
| 1593 | ! print *, mass_h2o_tot |
|---|
| 1594 | ! print *, real(mass_h2o_tot, 4) |
|---|
| 1595 | |
|---|
| 1596 | ! call write_output("dztot1","change in dztot", "unknown",real(dztot1)) |
|---|
| 1597 | |
|---|
| 1598 | call write_output("flux_soillayer","flux of water between the soil layers", "kg m-2 s-1",var_flux_soil) |
|---|
| 1599 | |
|---|
| 1600 | ! call write_output("A_coef","A coefficient of dztot1", "unknown",real(B)) |
|---|
| 1601 | |
|---|
| 1602 | ! call write_output("B_coef","B coefficient of dztot1", "unknown",real(B)) |
|---|
| 1603 | |
|---|
| 1604 | ! call write_output("H2O_init","initialized H2O", "kg/m^2",real(H2O)) |
|---|
| 1605 | |
|---|
| 1606 | ! call write_output("H2O_depth_init","initialized H2O depth ", "m",real(H2O_depth)) |
|---|
| 1607 | |
|---|
| 1608 | call write_output("ice_saturation_soil","Water ice saturation in the soil layers", "Percent",saturation_water_ice) |
|---|
| 1609 | |
|---|
| 1610 | ! call write_output("mass_h2o_tot","total mass of subsurface water over the planet", "kg",real(mass_h2o_tot)) |
|---|
| 1611 | |
|---|
| 1612 | ! call write_output("mass_ice_tot","total mass of subsurface ice over the planet", "kg",real(mass_ice_tot)) |
|---|
| 1613 | |
|---|
| 1614 | call write_output("mass_h2o_soil","Mass of subsurface water column at each point", "kg m-2",(mass_h2o)) |
|---|
| 1615 | |
|---|
| 1616 | call write_output("mass_ice_soil","Mass of subsurface ice at each point", "kg m-2",(mass_ice)) |
|---|
| 1617 | |
|---|
| 1618 | call write_output("znsoil","Water vapor soil concentration ", "kg.m - 3 of pore air",znsoil) |
|---|
| 1619 | |
|---|
| 1620 | call write_output("nsatsoil","subsurface water vapor saturation density", "kg/m^3",real(nsatsoil)) |
|---|
| 1621 | |
|---|
| 1622 | call write_output("nsurf","surface water vapor density", "kg/m^3",real(nsurf)) |
|---|
| 1623 | |
|---|
| 1624 | call write_output("adswater","subsurface adsorbed water", "kg/m^3",real(adswater)) |
|---|
| 1625 | |
|---|
| 1626 | ! call write_output("subsurfice","subsurface ice", "kg/m^3",real(ice)) |
|---|
| 1627 | |
|---|
| 1628 | call write_output("flux_rego",'flux of water from the regolith','kg/m^2',zdqsdifrego) |
|---|
| 1629 | |
|---|
| 1630 | ! call write_output("exchange",'exchange','no unit',real(exch)) |
|---|
| 1631 | |
|---|
| 1632 | ! call write_output("close",'close top, bottom, or none (1, -1, 0)','no unit',real(close_out)) |
|---|
| 1633 | |
|---|
| 1634 | ! call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_inter) |
|---|
| 1635 | |
|---|
| 1636 | !endif |
|---|
| 1637 | endif |
|---|
| 1638 | RETURN |
|---|
| 1639 | END |
|---|
| 1640 | |
|---|
| 1641 | |
|---|
