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