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