subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep, & exchange, qsat_surf, pq, pa, pb, pc, pd, pdqsdifpot, pqsurf, & pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2, saturation_water_ice) use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg, & ads_const_D, ads_massive_ice use comcstfi_h use tracer_mod use surfdat_h, only: watercaptag ! use mis par AP15 essai use geometry_mod, only: cell_area, latitude_deg use write_output_mod, only: write_output implicit none ! ================================================================================================= = ! Description: ! ! This subroutine calculates the profile of water vapor, adsorbed water and ice in the regolith, ! down to depth(nsoil) ~ 18m. It calculates the flux of water vapor (zdqsdifrego) between the ! suburface and the atmosphere. ! ! Water vapor and adsorbed water are treated as two separate subsurface tracers that are connected ! by adsorption / desorption coefficients, including an adsorption / desorption kinetic factor. ! ! The water adsorption isotherm is linear, with a saturation plateau corresponding to one monolayer ! of adsorbed water molecules (i.e. an approximation of a Langmuir - type isotherm). The slope of the ! isotherm is defined by the enthalpy of adsorption (in kJ / mol). ! 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) ! ! The linearized adsorption - diffusion equation is solved first, then the water vapor pressure is ! compared to the saturation vapor pressure, and if the latter is reached, ice is formed, the water vapor ! pressure is fixed to its saturation value and the adsorption - diffusion equation is solved again. ! If ice was already present, its sublimation or growth is calculated. ! ! This subroutine is called by vdifc.F if the parameters "regolithe" and "water" are set to true in ! callphys.def. It returns the flux of water vapor that is injected into or out of the regolith, ! and which serves as a boundary condition to calculate the vertical atmospheric profile of water vapor ! in vdifc. It also calculates the exchange between the subsurface and surface ice. ! ! It requires three subsurface tracers to be defined in traceur.def : ! - h2o_vap_soil (subsurface water vapor) ! - h2o_ice_soil (subsurface ice) ! - h2o_ads_vap (adsorbed water) ! ! Author : Pierre - Yves Meslin (pmeslin@irap.omp.eu) ! ! ================================================================================================= = ! Libraries : ! =========== = !#include "dimensions.h" !#include "dimphys.h" !#include "comsoil.h" #include "callkeys.h" !#include "comcstfi.h" !#include "tracer.h" !#include "watercap.h" ! Arguments : ! ====================================================================== ! Inputs : ! ---------------------------------------------------------------------- integer, intent(in) :: ngrid ! number of points in the model (all lat and long point in one 1D array) integer, intent(in) :: nlayer ! number of layers integer, intent(in) :: nq ! number of tracers integer, intent(in) :: nsoil integer, intent(in) :: nqsoil logical, save :: firstcall_soil = .true. ! triggers the initialization real, intent(in) :: ptsoil(ngrid, nsoil) ! Subsurface temperatures real, intent(in) :: ptsrf(ngrid) ! Surface temperature real, intent(in) :: ptimestep ! length of the timestep (unit depends on run.def file) logical, intent(in) :: exchange(ngrid) ! logical :: array describing whether there is exchange with the atmosphere at the current timestep real, intent(in) :: qsat_surf(ngrid) ! Saturation mixing ratio at the surface at the current surface temperature (formerly qsat) real, intent(in) :: pq(ngrid, nlayer, nq) ! Tracer atmospheric mixing ratio real, intent(in) :: pa(ngrid, nlayer) ! Coefficients za real, intent(in) :: pb(ngrid, nlayer) ! Coefficients zb real, intent(in) :: pc(ngrid, nlayer) ! Coefficients zc real, intent(in) :: pd(ngrid, nlayer) ! Coefficients zd real, intent(in) :: pdqsdifpot(ngrid) ! Potential pdqsdif (without subsurface - atmosphere exchange) real, intent(in) :: pplev(ngrid, nlayer+1) ! Atmospheric pressure levels real, intent(in) :: rhoatmo(ngrid) ! Atmospheric air density (1st layer) (not used right now) logical, intent(in) :: writeoutput ! just to check we are at the last subtimestep and we ! Variables modified : ! ---------------------------------------------------------------------- real, intent(inout) :: pqsoil(ngrid, nsoil, nqsoil) ! Subsurface tracers (water vapor and ice) real, intent(in) :: pqsurf(ngrid) ! Water ice on the surface ! (in kg.m - 3 : m - 3 of pore air for water vapor and m - 3 of regolith for water ice and adsorbed water) ! Outputs : ! ---------------------------------------------------------------------- real, intent(out) :: zdqsdifrego(ngrid) ! Flux from subsurface (positive pointing outwards) real, intent(out) :: zq1temp2(ngrid) ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg) real*8, intent(out) :: saturation_water_ice(ngrid, nsoil) ! Water pore ice saturation level (formerly Sw) ! Outputs for the output files real*8 :: preduite(ngrid, nsoil) ! how close the water vapor density is to adsorption saturation integer :: exch(ngrid) ! translates the logical variable exchange into a -1 or 1 real*8 :: mass_h2o(ngrid) ! Mass of subsurface water column at each point (kg.m - 2) (formerly mh2otot) real*8 :: mass_ice(ngrid) ! Mass of subsurface ice at each point (kg.m - 2) (formerly micetot) real*8 :: mass_h2o_tot ! Mass of subsurface water over the whole planet real*8 :: mass_ice_tot ! Mass of subsurface ice over the whole planet real*8 :: nsurf(ngrid) ! surface tracer density in kg/m^3 real*8 :: close_out(ngrid, nsoil) ! output for close_top and close_bottom ! Local (saved) variables : ! ====================================================================== real*8 :: P_sat_soil(ngrid, nsoil) ! water saturation pressure of subsurface cells (at miD-layers) (formerly Psatsoil) real*8 :: nsatsoil(ngrid, nsoil) ! gas number density at saturation pressure real*8, allocatable, save :: znsoil(:, :) ! Water vapor soil concentration (kg.m - 3 of pore air) real*8 :: znsoilprev(ngrid, nsoil) ! value of znsoil in the previous timestep real*8 :: znsoilprev2(ngrid, nsoil) ! second variable for the value of znsoil in the previous timestep real*8 :: zdqsoil(ngrid, nsoil) ! Increase of pqsoil if sublimation of ice real*8, allocatable, save :: ice(:, :) ! Ice content of subsurface cells (at miD-layers) (kg / m3) real*8 :: iceprev(ngrid, nsoil) logical :: ice_index_flag(nsoil) ! flag for ice presence real*8, allocatable, save :: adswater(:, :) ! Adsorbed water in subsurface cells (at miD-layers) (...) real*8 :: adswater_temp(ngrid, nsoil) ! temprory variable for adsorbed water 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 real*8 :: adswprev(ngrid, nsoil) logical :: recompute_cell_ads_flag(nsoil) ! ARNAU ! Formerly ads_water_saturation_flag_1 but with a twist. This variable now ! checks whether coefficients have changed and need recomputing. If the cell ! is seen to be over the monolayer saturation level (i.e. the cell fulfills the ! condition adswater_temp(ig, ik) > adswater_sat_mono) but the coefficients ! have been computed assuming that the cell is below the monolayer saturation ! layer (i.e. the variable over_mono_sat_flag(ig, ik) had been set to .false.), ! then the cell needs to have its coefficients recomputed according to the ! previous condition: adswater_temp(ig, ik) > adswater_sat_mono. Then, ! the variable recompute_cell_ads_flag(ik) becomes .true.. ! ARNAU real*8, save :: adswater_sat_mono ! Adsorption monolayer saturation value (kg.m - 3 of regolith) ! ARNAU real*8 :: delta0(ngrid) ! Coefficient delta(0) modified real*8 :: alpha0(ngrid) real*8 :: beta0(ngrid) ! Adsorbtion/Desorption variables and parameters real*8 :: Ka(ngrid, nsoil) ! Adsorption time constant (s - 1) before monolayer saturation (first segment of the bilinear function) real*8 :: Kd(ngrid, nsoil) ! Desorption time constant (s - 1) before monolayer saturation (first segment of the bilinear function) 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) real*8 :: Ka2(ngrid, nsoil) ! Adsorption time constant (s - 1) after monolayer saturation (second segment of the bilinear function) ! modified 2020 real*8 :: Kd2(ngrid, nsoil) ! Desorption time constant (s - 1) after monolayer saturation (second segment of the bilinear function) ! modified 2020 real*8 :: k_ads_eq2(ngrid, nsoil) ! Equilibrium adsorption coefficient (formerly kads) after monolayer saturation (second segment of the bilinear function); unitless ! modified 2020 real*8 :: c0(ngrid, nsoil) ! Intercept of the second line in the bilinear function ! modified 2020 real*8, parameter :: kinetic_factor = 0.01 ! (inverse of) Characteristic time to reach adsorption equilibrium (s - 1): ! fixed arbitrarily when kinetics factors are unknown ! possible values: ! = 1 / 1800 s - 1 ! / 1.16D-6 / !( = 10 earth days)! / 1.D0 / ! / 1.2D-5 / ! real*8, allocatable, save :: ztot1(:, :) ! Total (water vapor + ice) content (kg.m - 3 of soil) real*8 :: dztot1(nsoil) real*8 :: h2otot(ngrid, nsoil) ! Total water content (ice + water vapor + adsorbed water) real*8 :: flux(ngrid, 0:nsoil - 1) ! Fluxes at interfaces (kg.m - 2.s - 1) (positive = upward) real*8 :: rho(ngrid) ! Air density (first subsurface layer) real*8 :: rhosurf(ngrid) ! Surface air density ! Porosity and tortuosity variables real*8, allocatable, save :: porosity_ice_free(:, :) ! Porosity with no ice present (formerly eps0) real*8, allocatable, save :: porosity(:, :) ! Porosity (formerly eps) real*8 :: porosity_prev(ngrid, nsoil) ! Porosity from previous timestep real*8, allocatable, save :: tortuosity(:, :) ! Tortuosity factor (formerly tortuo) real*8 :: saturation_water_ice_inter(ngrid, nsoil) ! Water pore ice saturation level at the interlayer ! Diffussion coefficients real*8 :: D_mid(ngrid, nsoil) ! Midlayer diffusion coefficients real*8 :: D_inter(ngrid, nsoil) ! Interlayer diffusion coefficients (formerly D) real*8, allocatable, save :: D0(:, :) ! Diffusion coefficient prefactor : ! If (medium = 1) : D0 = value of D_mid for saturation_water_ice = 0, ie. poro / tortuo * Dm (in m2 / s) ! If (medium = 2) : D0 = porosity_tortuosity factor (dimensionless) real*8 :: omega(ngrid, nsoil) ! H2O - CO2 collision integral real*8 :: vth(ngrid, nsoil) ! H2O thermal speed real*8, parameter :: Dk0 = 0.459D0 ! Knudsen diffusion coefficient (for saturation_water_ice = 0) Value for a (5 / 1) bidispersed random assembly of spheres ! (dimensionless, ie. divided by thermal speed and characteristic meshsize of the porous medium) real*8 :: Dk(ngrid, nsoil) ! Knudsen diffusion coefficient (divided by porosity / tortuosity factor) real*8 :: Dk_inter(ngrid, nsoil) ! Knudsen diffusion coefficient at the interlayer real*8 :: Dm(ngrid, nsoil) ! Molecular diffusion coefficient real*8 :: Dm_inter(ngrid, nsoil) ! Molecular diffusion coefficient at the interlayer real*8, parameter :: choke_fraction = 0.8D0 ! fraction of ice that prevents further diffusion logical, allocatable, save :: close_top(:, :) ! closing diffusion at the top of a layer if ice rises over saturation logical, allocatable, save :: close_bottom(:, :) ! closing diffusion at the bottom of a layer if ice risies over saturation logical, parameter :: print_closure_warnings = .true. ! Coefficients for the Diffusion calculations real*8 :: A(ngrid, nsoil) ! Coefficient in the diffusion formula real*8 :: B(ngrid, nsoil) ! Coefficient in the diffusion formula real*8 :: C(ngrid, nsoil) ! Coefficient in the diffusion formula real*8 :: E(ngrid, nsoil) ! Coefficient in the diffusion formula (before monolayer saturation) ! added 2020 real*8 :: F(ngrid, nsoil) ! Coefficient in the diffusion formula (before monolayer saturation) ! added 2020 real*8 :: E2(ngrid, nsoil) ! Coefficient in the diffusion formula (after monolayer saturation) ! added 2020 real*8 :: F2(ngrid, nsoil) ! Coefficient in the diffusion formula (after monolayer saturation) ! added 2020 real*8, allocatable, save :: zalpha(:, :) ! Alpha coefficients real*8, allocatable, save :: zbeta(:, :) ! Beta coefficients real*8 :: zdelta(ngrid, nsoil - 1) ! Delta coefficients ! Distances between layers real*8, allocatable, save :: interlayer_dz(:, :) ! Distance between the interlayer points in m (formerly interdz) real*8, allocatable, save :: midlayer_dz(:, :) ! Distance between the midcell points in m (formerly middz) real*8 :: nsat(ngrid, nsoil) ! amount of water vapor at (adsorption) monolayer saturation ! modified 2020 real*8, allocatable, save :: meshsize(:, :) ! scaling/dimension factor for the por size real*8, allocatable, save :: rho_soil(:, :) ! Soil density (bulk - kg / m3) (formerly rhosoil) real*8, allocatable, save :: cste_ads(:, :) ! Prefactor of adsorption coeff. k ! general constants real*8, parameter :: kb = 1.38065D-23 ! Boltzmann constant real*8, parameter :: mw = 2.988D-26 ! Water molecular weight !real*8, parameter :: rho_H2O_ice = 920.D0 ! Ice density ! adsorption coefficients real*8, parameter :: enthalpy_ads = 35.D3 ! Enthalpy of adsorption (J.mole - 1) options 21.D3, 35.D3, and 60.D3 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 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 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 real*8, parameter :: tau0 = 1D-14 real*8, parameter :: S = 17.D3 ! Soil specific surface area (m2.kg - 1 of solid) options: 17.D3 and 106.D3 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) ! Reference values for choice_ads = 2 real*8, parameter :: Tref = 243.15D0 real*8, parameter :: nref = 2.31505631D-6 ! calculated as 18.D-3 / (8.314D0 * Tref) * 0.26D0 ! not used anymore (for the time being) 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) real*8, parameter :: Ka_ref = 2.4D-4 ! Not used for the time being ! 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 ! real*8, parameter :: Kref2 = 5.95D-7 ! Value from data analysis Pommerol2009 ! ARNAU 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) 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) logical :: recompute_all_cells_ads_flag ! Old saturation_flag but with a behaviour change ! ARNAU ! The old saturation_flag was used to check whether the cell was saturated or ! not, and therefore to assign the saturation value adswater(ig, ik) = ! adswater_sat instead of the usual adswater(ig, ik) = adswater_temp(ig, ik) ! The old routine using saturation_flag did not require that the A, B, C,... ! adsorption coefficients be computed, but the new soilwater subroutine ! does. Therefore, the variable recompute_all_cells_ads_flag checks whether ! there is a cell of a column that requires recomputing. If at least one cell ! requires recomputing (i.e. recompute_cell_ads_flag(ik) is .true.), then ! recompute_all_cells_ads_flag becomes true, and the adsorption coefficients, ! as well as the recursive equation (appearing in this code as `backward ! loop over all layers`) ! ARNAU logical :: sublimation_flag logical :: condensation_flag(nsoil) ! variables and parameters for the H2O map integer, parameter :: n_long_H2O = 66!180 ! number of longitudes of the new H2O map integer, parameter :: n_lat_H2O = 50 !87 ! number of latitudes of the new H2O map real*8, parameter :: rho_H2O_ice = 920.0D0 ! Ice density (formerly rhoice) real :: latH2O(n_lat_H2O*n_long_H2O) ! Latitude at H2O map points real :: longH2O(n_lat_H2O*n_long_H2O) ! Longitude at H2O map points real :: H2O_depth_data(n_lat_H2O*n_long_H2O) ! depth of the ice layer in in g/cm^2 at H2O map points real :: H2O_cover_data(n_lat_H2O*n_long_H2O) ! H2O content of the cover layer at H2O map points (not used yet) real :: dataH2O(n_lat_H2O*n_long_H2O) ! H2O content of the ice layer at H2O map points real :: latH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable real :: longH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable real :: dataH2O_temp(n_lat_H2O*n_long_H2O) ! intermediate variable real :: H2O_depth_data_temp(n_lat_H2O*n_long_H2O) ! intermediate variable real, allocatable, save :: H2O(:) ! H2O map interpolated at GCM grid points (in wt%) real, allocatable, save :: H2O_depth(:) ! H2O map ice depth interpolated at GCM in g/cm^2 ! variables for the 1D case real*8, parameter :: mmr_h2o = 0.6D-4 ! Water vapor mass mixing ratio (for initialization only) real*8 :: diff(ngrid, nsoil) ! difference between total water content and ice + vapor content ! (only used for output, inconstistent: should just be adswater) real :: var_flux_soil(ngrid, nsoil) ! Output of the flux between soil layers (kg/m^3/s) real :: default_diffcoeff = 4e-4 ! Diffusion coefficient by default (no variations with Temperature and pressure (m/s^2) real :: tol_massiveice = 50. ! Tolerance to account for massive ice (kg/m^3) ! Loop variables and counters integer :: ig, ik, i, j ! loop variables logical :: output_trigger ! used to limit amount of written output integer, save :: n ! number of runs of this subroutine integer :: sublim_n ! number of sublimation loop runs integer :: satur_mono_n ! number of monolayer saturation loop runs ! added 2020 if (.not. allocated(znsoil)) then allocate( znsoil(ngrid, nsoil) ) allocate( ice(ngrid, nsoil) ) allocate( adswater(ngrid, nsoil) ) allocate( ztot1(ngrid, nsoil) ) allocate( porosity_ice_free(ngrid, nsoil) ) allocate( porosity(ngrid, nsoil) ) allocate( tortuosity(ngrid, nsoil) ) allocate( D0(ngrid, nsoil) ) allocate( interlayer_dz(ngrid, nsoil) ) allocate( midlayer_dz(ngrid, 0:nsoil) ) ! allocate( zalpha(ngrid, nsoil-1) ) ! allocate( zbeta(ngrid, nsoil-1) ) allocate( zalpha(ngrid, nsoil) ) ! extra entry to match output format allocate( zbeta(ngrid, nsoil) ) ! extra entry to match output format allocate( meshsize(ngrid, nsoil) ) allocate( rho_soil(ngrid, nsoil) ) allocate( cste_ads(ngrid, nsoil) ) allocate( H2O(ngrid) ) allocate( H2O_depth(ngrid) ) allocate( close_top(ngrid, nsoil) ) allocate( close_bottom(ngrid, nsoil) ) allocate( over_mono_sat_flag(ngrid, nsoil) ) endif ! Used commons ! mugaz ! Molar mass of the atmosphere (g.mol-1) ~43.49 ! Initialisation : ! ================= if (firstcall_soil) then n = 0 close_top = .false. close_bottom = .false. print *, "adsorption enthalpy first segment: ", enthalpy_ads print *, "adsorption enthalpy second segment: ", enthalpy_ads2 print *, "adsorption BET constant C first segment: ", DeltaQ_ads print *, "adsorption BET constant C second segment: ", DeltaQ_ads2 print *, "specific surface area: ", S do ig = 1, ngrid ! if(.not.watercaptag(ig)) then ! Initialization of soil parameters ! =================================== ! initialize the midlayer distances midlayer_dz(ig, 0) = mlayer(0) do ik = 1, nsoil - 1 midlayer_dz(ig, ik) = mlayer(ik) - mlayer(ik - 1) enddo ! initialize the interlayer distances interlayer_dz(ig, 1) = layer(1) do ik = 2, nsoil interlayer_dz(ig, ik) = layer(ik) - layer(ik - 1) enddo ! Choice of porous medium and D0 ! =============================== = do ik = 1, nsoil ! These properties are defined here in order to enable custom profiles if(ads_massive_ice) then if(pqsoil(ig, ik, igcm_h2o_ice_soil).gt.tol_massiveice) then porosity_ice_free(ig, ik) = 0.999999 else porosity_ice_free(ig, ik) = porosity_reg endif else porosity_ice_free(ig, ik) = porosity_reg endif tortuosity(ig, ik) = 1.5D0 rho_soil(ig, ik) = 1.3D3 ! in kg/m3 of regolith (incl. porosity) meshsize(ig, ik) = 5.D-6 ! Characteristic size 1e - 6 = 1 micron : diameter(mode 5) = 10 microns, diameter(mode 1) = 1 microns ! Example : with meshsize = 5.E - 6 : grain sizes = 50 and 10 microns D0(ig, ik) = porosity_ice_free(ig, ik) / tortuosity(ig, ik) enddo ! Choice of adsorption coefficients ! =================================== 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 ! 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) ! theoretical formula is = rho_soil(ig, 1) * S / Sm * mw ! Numerical values above are for SSA = 13.7 m2 / g and plateau = 0.8wt% ! adswater_sat_mono = 6D-2 * rho_soil(ig, 1) ! Experimental value for JSC Mars - 1 (Pommerol et al., 2009) ! with SSA = 106 m2 / g and plateau (net) = 6wt% ! Initialisation of ice, water vapor and adsorbed water profiles ! =============================================================== = do ik = 1, nsoil ! Initialisation of pqsoil(t = 0) ! in 1D simulations the initialization happens here ! if (ngrid.eq.1) then ! znsoil(ig, ik) = mmr_h2o * mugaz * 1.D-3 * pplev(ig, 1) & ! / (8.314D0 * ptsoil(ig, nsoil - 4)) ! ! znsoil(ig, ik) = 0. ! ice(ig, ik) = 0.D0 ! 0.1D0 !414.D0 ! saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) ! porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) ! if (choice_ads.eq.1) then ! vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, nsoil - 4) & ! / (pi * 18.D-3)) ! first segment of the bilinear function ! k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & ! / (dexp(-enthalpy_ads & ! / (8.314D0 * ptsoil(ig, nsoil - 4))) / tau0) ! Kd(ig, ik) = kinetic_factor & ! / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) ! Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / & ! (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) ! second segment of the bilinear function ! added 2020 ! k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & ! / (dexp(-enthalpy_ads2 & ! / (8.314D0 * ptsoil(ig, nsoil - 4))) / tau0) ! Kd2(ig, ik) = kinetic_factor & ! / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / & ! (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! ! elseif (choice_ads.eq.2) then ! par rapport \E0 valeur experimentale de reference ! ! Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 * ! ! & (1. / ptsoil(ig, nsoil - 4) - 1. / Tref)) ! ! Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref * ! ! & sqrt(ptsoil(ig, nsoil - 4) / Tref) / nref! ! ! ! first segment of the bilinear function ! 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 ! ! Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) ! ! Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) ! ! ! second segment of the bilinear function ! added 2020 ! 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 ! ! Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! ! Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! endif ! if (k_ads_eq(ig,ik) * znsoil(ig, ik) .gt. adswater_sat_mono) then ! modified 2024, after correction of c0 expression ! c0(ig, ik) = ( 1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono ! adswater(ig, ik) = c0(ig,ik) + k_ads_eq2(ig, ik) * znsoil(ig, ik) ! else ! adswater(ig, ik) = k_ads_eq(ig, ik) * znsoil(ig, ik) ! endif ! else ! in 3D simulations initialisation happens with newstart.F znsoil(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_soil) ice(ig, ik) = pqsoil(ig, ik, igcm_h2o_ice_soil) adswater(ig, ik) = pqsoil(ig, ik, igcm_h2o_vap_ads) ! endif saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) if (saturation_water_ice(ig, ik).lt.0.D0) then print *, "WARNING!! soil water ice negative at ", ig, ik print *, "saturation value: ", saturation_water_ice(ig, ik) print *, "setting saturation to 0" saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) endif porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) ! Initialise soil as if all points where below the monolayer saturation level (added 2020) => now depends on value of adswater (modified in 2024) if (adswater(ig,ik).gt.adswater_sat_mono) then over_mono_sat_flag(ig, ik) = .true. ! else over_mono_sat_flag(ig, ik) = .false. endif enddo ! endif enddo print *, "initializing H2O data" ! ! ! 1.6 intitalization of the water content ! open(40,file='H2O_data') ! do i = 1, n_lat_H2O*n_long_H2O ! read(40,*) latH2O_temp(i), longH2O_temp(i), H2O_cover_data(i), dataH2O_temp(i), H2O_depth_data_temp(i) ! ! write(*, *) 'depth data ', H2O_depth_data(i) ! enddo ! close(40) ! ! ! 1.6.2. put latH2O, longH2O and dataH2O in the right format to pass on to mapTh ! ! in the datafile the latitudes are listed from negative to positive, but mapTh needs them the other way around ! do i = 0, n_lat_H2O - 1 ! do j = 1, n_long_H2O ! latH2O( n_long_H2O * i + j ) = latH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j) ! longH2O( n_long_H2O * i + j ) = longH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j) ! dataH2O( n_long_H2O * i + j ) = dataH2O_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j) ! H2O_depth_data( n_long_H2O * i + j ) = H2O_depth_data_temp( n_long_H2O * (n_lat_H2O - (i + 1)) + j) ! enddo ! enddo ! ! call mapTh(latH2O, longH2O, dataH2O, n_long_H2O, n_lat_H2O, ngrid, H2O) ! call mapTh(latH2O, longH2O, H2O_depth_data, n_long_H2O, n_lat_H2O, ngrid, H2O_depth) ! do ig = 1, ngrid ! convert depth from g/cm^2 to m ! print *, H2O_depth(ig), rho_soil(ig, 1) ! H2O_depth(ig) = H2O_depth(ig) * 10 / rho_soil(ig, 1) ! if ( (latitude_deg(ig).lt.40).and.(latitude_deg(ig).gt.-40)) then ! H2O(ig) = 0. ! endif output_trigger = .true. do ik = 1, nsoil ! if (H2O_depth(ig).le.layer(ik)) then ! ice(ig, ik) = min(H2O(ig), rho_H2O_ice * porosity_ice_free(ig, ik)) ! if (output_trigger) then ! print *, "ice set at: ", ig, ik, "to :", ice(ig,ik), "depth: ", H2O_depth(ig) ! output_trigger = .false. ! endif ! endif saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) enddo enddo print *, "Kinetic factor: ", kinetic_factor if (kinetic_factor.eq.0) then print *, "WARNING: adsorption is switched off" endif firstcall_soil = .false. endif ! of "if firstcall_soil" ! Update properties in case of the sublimation of massive ice if(ads_massive_ice) then do ig = 1, ngrid do ik = 1, nsoil if(pqsoil(ig, ik, igcm_h2o_ice_soil).gt.tol_massiveice) then porosity_ice_free(ig, ik) = 0.999999 else porosity_ice_free(ig, ik) = porosity_reg endif enddo enddo endif ! Computation of new (new step) transport coefficients : ! ======================================================= = do ig = 1, ngrid ! loop over all gridpoints if(.not.watercaptag(ig)) then ! if there is no polar cap do ik = 1, nsoil ! calculate Water saturation level (pore volume filled by ice / total pore volume) saturation_water_ice(ig, ik) = min(ice(ig, ik) / (rho_H2O_ice * porosity_ice_free(ig, ik)), 0.999D-0) if (saturation_water_ice(ig, ik).lt.0.D0) then print *, "WARNING!! soil water ice negative at ", ig, ik print *, "saturation value: ", saturation_water_ice(ig, ik) print *, "setting saturation to 0" saturation_water_ice(ig, ik) = max(saturation_water_ice(ig, ik), 0.D0) endif ! save the old porosity porosity_prev(ig, ik) = porosity(ig, ik) ! calculate the new porosity porosity(ig, ik) = porosity_ice_free(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) ! Note : saturation_water_ice and porosity are computed from the ice content of the previous timestep enddo ! calculate the air density in the first subsurface layer rho(ig) = pplev(ig, 1) / (r * ptsoil(ig, 1)) ! calculate diffusion coefficients at mid- and interlayers (with ice content from the previous timestep) ! Dependence on saturation_water_ice taken from the relation obtained by Hudson et al. (2009) and Meslin et al. (SSSAJ, 2010) do ik = 1, nsoil ! loop over all soil levels ! calculate H2O mean thermal velocity vth(ig, ik) = dsqrt(8.D0 * 8.314D0 * ptsoil(ig, ik) / (pi * 18.D-3)) ! The diffusion coefficient is calculated ! calculate the H2O - CO2 collision integral (after Mellon and Jakosky, JGR, 1993) omega(ig, ik) = 1.442D0 - 0.673D0 * dlog(2.516D-3 * ptsoil(ig, ik)) & + 0.252D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 2.D0 & - 0.047D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 3.D0 & + 0.003D0 * (dlog(2.516D-3 * ptsoil(ig, ik))) ** 4.D0 ! calculate the molecular diffusion coefficient Dm(ig, ik) = 4.867D0 * ptsoil(ig, ik) ** (3.D0 / 2.D0) & / (pplev(ig, 1) * omega(ig, ik)) * 1.D-4 ! calculate the Knudsen diffusion coefficient (divided by D0 = porosity / tortuosity in this case) Dk(ig, ik) = Dk0 / D0(ig, ik) * meshsize(ig, ik) * vth(ig, ik) ! calculate the midlayer diffusion coeff. (with dependence on saturation_water_ice from Meslin et al., 2010 - only exact for normal diffusion) 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)) if(ads_const_D) D_mid(ig, ik) = default_diffcoeff ! Midlayer diffusion coeff. (correlation by Rogers and Nielson, 1991) ! D_mid(ig, ik) = D0(ig, ik) * (1. - saturation_water_ice(ig, ik)) * exp(-6. * saturation_water_ice(ig, ik) & ! * porosity_ice_free(ig, ik) - 6. * saturation_water_ice(ig, ik) ** (14. * porosity_ice_free(ig, ik))) & ! * 1. / (1. / Dm(ig, ik) + 1. / Dk(ig, ik)) enddo ! calculate the interlayer diffusion coefficient as interpolation of the midlayer coefficients do ik = 1, nsoil - 1 Dm_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * Dm(ig, ik + 1) & + (mlayer(ik) - layer(ik)) * Dm(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1)) Dk_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * Dk(ig, ik + 1) & + (mlayer(ik) - layer(ik)) * Dk(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1)) ! saturation_water_ice_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * saturation_water_ice(ig, ik + 1) & ! + (mlayer(ik) - layer(ik)) * saturation_water_ice(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1)) if (close_bottom(ig, ik).or.close_top(ig, ik+1)) then saturation_water_ice_inter(ig, ik) = 0.999D0 else saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik)) ! new diffusion interaction endif saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik)) ! new diffusion interaction 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)) if(ads_const_D) D_inter(ig, ik) = default_diffcoeff ! D_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * D_mid(ig, ik + 1) & ! old implementation with direct interpolation ! + (mlayer(ik) - layer(ik)) * D_mid(ig, ik) ) / (mlayer(ik) - mlayer(ik - 1)) enddo endif enddo ! Computation of new adsorption / desorption constants (Ka, Kd coefficients), and C, E, F coefficients ! =================================================================================================== do ig = 1, ngrid ! loop over all gridpoints if(.not.watercaptag(ig)) then ! if there is no polar cap do ik = 1, nsoil ! loop over all soil levels if (choice_ads.eq.1) then ! Constant, uniform diffusion coefficient D0 (prescribed) ! first segment of the bilinear function (before monolayer saturation) ! calculate the equilibrium adsorption coefficient k_ads_eq(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & / (dexp(-enthalpy_ads / (8.314D0 * ptsoil(ig, ik))) / tau0) ! calculate the desorption coefficient Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) ! calculate the absorption coefficient Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) ! second segment of the bilinear function (after monolayer saturation) ! added 2020 ! calculate the equilibrium adsorption coefficient k_ads_eq2(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4.D0 & / (dexp(-enthalpy_ads2 / (8.314D0 * ptsoil(ig, ik))) / tau0) ! added 2020 ! calculate the desorption coefficient Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! added 2020 ! calculate the absorption coefficient Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! added 2020 ! Kd(ig, ik) = exp(-enthalpy_ads / (8.314 * ptsoil(ig, ik))) ! & / tau0 ! * 1.D-18 ! Ka(ig, ik) = rho_soil(ig, ik) * S * vth(ig, ik) / 4. ! * 1.D-18 ! if ((ngrid.eq.1).and.(ik.eq.18)) then ! if 1D simulation and uppermost layer ! print * , 'Ka, Kd, Ka / Kd = ', Ka(ig, ik), Kd(ig, ik), Ka(ig, ik) / Kd(ig, ik) ! print * , 'k_ads_eq = ', k_ads_eq(ig, ik) ! print * , 'porosity = ', porosity(ig, ik) ! endif elseif (choice_ads.eq.2) then ! modified 2020 (with DeltaQ instead of Q) ! Kd(ig, ik) = Kd_ref * exp(-enthalpy_ads / 8.314 * ! & (1. / ptsoil(ig, ik) - 1. / Tref)) ! Ka(ig, ik) = rho_soil(ig, ik) * Ka_ref * sqrt(ptsoil(ig, ik) / Tref) ! & / nref ! first segment of the bilinear function (before monolayer saturation) ! calculate the equilibrium adsorption coefficient 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 ! calculate the desorption coefficient Kd(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) ! calculate the absorption coefficient Ka(ig, ik) = kinetic_factor * k_ads_eq(ig, ik) / (1.D0 + k_ads_eq(ig, ik) / porosity(ig, ik)) ! second segment of the bilinear function (after monolayer saturation) ! added 2020 ! calculate the equilibrium adsorption coefficient 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 ! calculate the desorption coefficient Kd2(ig, ik) = kinetic_factor / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) ! calculate the absorption coefficient Ka2(ig, ik) = kinetic_factor * k_ads_eq2(ig, ik) / (1.D0 + k_ads_eq2(ig, ik) / porosity(ig, ik)) elseif (choice_ads.eq.0) then!No adsorption/desorption Ka(ig, ik) = 0. Kd(ig, ik) = 0. Ka2(ig, ik) = 0. Kd2(ig, ik) = 0. endif ! calculate the amount of water vapor at monolayer saturation ! modified 2020 if(choice_ads.ne.0) nsat(ig, ik) = adswater_sat_mono * Kd(ig, ik) / Ka(ig, ik) ! calculate the intercept of the second line in the bilinear function ! added 2020; corrected 2024 c0(ig, ik) = ( 1.D0 - (k_ads_eq2(ig, ik) / k_ads_eq(ig, ik))) * adswater_sat_mono ! calculate C, E, and F coefficients for later calculations C(ig, ik) = interlayer_dz(ig, ik) / ptimestep ! first segment of the bilinear function (before monolayer saturation) E(ig, ik) = Ka(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik)) F(ig, ik) = Kd(ig, ik) / (1.D0 + ptimestep * Kd(ig, ik)) ! second segment of the bilinear function (after monolayer saturation) ! added 2020 E2(ig, ik) = Ka2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik)) F2(ig, ik) = Kd2(ig, ik) / (1.D0 + ptimestep * Kd2(ig, ik)) ! save the old water concentration znsoilprev(ig, ik) = znsoil(ig, ik) znsoilprev2(ig, ik) = znsoil(ig, ik) ! needed in "special case" loop because adswprev can be changed before this loop adswprev(ig, ik) = adswater(ig, ik) !!!! Besoin de adswprev2 ??? iceprev(ig, ik) = ice(ig, ik) ! needed in "special case" loop for the same reason ! calculate the total ice + water vapor content ztot1(ig, ik) = porosity_prev(ig, ik) * znsoil(ig, ik) + ice(ig, ik) enddo endif enddo ! Computation of delta, alpha and beta coefficients ETC !!!! ! =========================================================== do ig = 1, ngrid ! loop over all gridpoints ! initialize the flux to zero to avoid undefined values zdqsdifrego(ig) = 0.D0 do ik = 0, nsoil - 1 flux(ig, ik) = 0. enddo if(.not.watercaptag(ig)) then ! if there is no polar cap do ik = 1, nsoil ! loop over all soil levels if (ice(ig, ik).eq.0.) then ice_index_flag(ik) = .false. else ice_index_flag(ik) = .true. endif enddo do ik = 1, nsoil ! loop over all soil levels ! calculate A and B coefficients ! modified 2020 if ( .not.over_mono_sat_flag(ig, ik) ) then ! Assume cell below the monolayer saturation ! calculate A and B coefficients (first segment of the bilinear function) A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) else ! Assume cell over the monolayer saturation ! added 2020 ! calculate A and B coefficients (second segment of the bilinear function) A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 endif ! calculate the saturation pressure 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) ? ! calculate the gas number density at saturation pressure nsatsoil(ig, ik) = (P_sat_soil(ig, ik) * mw) / (kb * ptsoil(ig, ik)) enddo ! calculate the alpha, beta, and delta coefficients for the surface layer delta0(ig) = pa(ig, 1) + pb(ig, 2) * (1.D0 - pd(ig, 2)) + porosity_ice_free(ig, 1) * pb(ig, 1) alpha0(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / delta0(ig) beta0(ig) = ( pb(ig, 2) * pc(ig, 2) + pq(ig, 1, igcm_h2o_vap) * pa(ig, 1) + pqsurf(ig) ) & / delta0(ig) ! set loop flag do ik = 1, nsoil condensation_flag = .false. enddo sublimation_flag = .true. sublim_n = 0 do while (sublimation_flag) ! while there is sublimation ! print *, "sublimation loop: ", sublim_n sublim_n = sublim_n + 1 if (sublim_n.gt.100) then print *, "sublimation not converging in call ", n print *, "might as well stop now" stop endif ! reset flag sublimation_flag = .false. recompute_all_cells_ads_flag = .true. satur_mono_n = 0 do while (recompute_all_cells_ads_flag) ! print *, satur_mono_n satur_mono_n = satur_mono_n + 1 ! reset loop flag recompute_all_cells_ads_flag = .false. ! calculate the surface air density rhosurf(ig) = pplev(ig, 1) / (r * ptsrf(ig)) ! calculate the coefficients in the upermost layer if (exchange(ig)) then ! if there is surface exchange ! calculate the delta, alpha, and beta coefficients zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) & + D_inter(ig, 1) / midlayer_dz(ig, 1) & !!! est - ce que le terme pb inclut le bon rho ? + porosity_ice_free(ig, 1) * pb(ig, 1) / (rho(ig) * ptimestep) * (1.D0 - alpha0(ig)) !pourrait remplacer pb/(ptime*rho(1/2)) par zcdv/ptime zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1) zbeta(ig, 1) = ( porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * beta0(ig) & + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1) ) & / zdelta(ig, 1) else ! calculate the delta, alpha, and beta coefficients without exchange zdelta(ig, 1) = A(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1) & + D_mid(ig, 1) / midlayer_dz(ig, 0) & + porosity(ig, 1) * C(ig, 1) zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1. / zdelta(ig, 1) zbeta(ig, 1) = ( D_mid(ig, 1) / midlayer_dz(ig, 0) * qsat_surf(ig) * rhosurf(ig) & + porosity_prev(ig, 1) * C(ig, 1) * znsoilprev(ig, 1) + B(ig, 1) ) & / zdelta(ig, 1) endif ! check for ice if (ice_index_flag(1)) then ! set the alpha coefficient to zero zalpha(ig, 1) = 0.D0 ! set the beta coefficient to the number density at saturation pressure zbeta(ig, 1) = nsatsoil(ig, 1) endif do ik = 2, nsoil - 1 ! go through all other layers ! calculate delta, alpha, and beta coefficients zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) & + D_inter(ig, ik) / midlayer_dz(ig, ik) & + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1)) zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik) zbeta(ig, ik) = ( D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) & + B(ig, ik) & + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev(ig, ik) ) / zdelta(ig, ik) ! check for ice if (ice_index_flag(ik)) then ! set the alpha coefficient to zero zalpha(ig, ik) = 0.D0 ! set the beta coefficient to the number density at saturation pressure zbeta(ig, ik) = nsatsoil(ig, ik) endif enddo ! Computation of pqsoil, ztot1, zq1temp2 and zdqsdifrego ! ======================================================= = ! calculate the preliminary amount of water vapor in the bottom layer if (ice_index_flag(nsoil)) then ! if there is ice ! set the vapor number density to saturation znsoil(ig, nsoil) = nsatsoil(ig, nsoil) else ! calculate the vapor number density in the lowest layer znsoil(ig, nsoil) = ( D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) & + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev(ig, nsoil) + B(ig, nsoil) ) & / ( porosity(ig, nsoil) * C(ig, nsoil) + A(ig, nsoil) & + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)) ) endif ! calculate the preliminary amount of adsorbed water if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024 adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) & / (1.D0 + ptimestep * Kd(ig, nsoil)) else adswater_temp(ig, nsoil) = ( Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) + ptimestep * c0(ig,nsoil) * Kd2(ig,nsoil)) & / (1.D0 + ptimestep * Kd2(ig, nsoil)) endif do ik = nsoil-1, 1, -1 ! backward loop over all layers znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik) if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024 adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) & / (1.D0 + ptimestep * Kd(ig, ik)) else adswater_temp(ig, ik) = ( Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) + ptimestep * c0(ig,ik) * Kd2(ig,ik)) & / (1.D0 + ptimestep * Kd2(ig, ik)) endif enddo ! check if any cell is over monolayer saturation do ik = 1, nsoil ! loop over all soil levels ! 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 ! recompute_cell_ads_flag(ik) = .true. ! cycle ! Jump loop ! endif ! Change of coefficients ! ARNAU recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients if ( adswater_temp(ig, ik).ge.adswater_sat_mono ) then print *, "NOTICE: over monolayer saturation" if ( .not.over_mono_sat_flag(ig, ik) ) then ! If the cell enters this scope, it ! means that the cell is over the monolayer ! saturation after calculation, but the ! coefficients assume it is below. Therefore, ! the cell needs to be recomputed over_mono_sat_flag(ig, ik) = .true. recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients ! 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) A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 endif else if ( over_mono_sat_flag(ig, ik) ) then ! If the cell enters this scope, it ! means that the cell is below the monolayer ! saturation after calculation, but the ! coefficients assume it is above. Therefore, ! the cell needs to be recomputed over_mono_sat_flag(ig, ik) = .false. recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients ! calculate A and B coefficients (reset to first segment of the bilinear function) A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) endif endif enddo ! if one cell needs to be recomputed, then all the column is to be recomputed too do ik = 1, nsoil if ( recompute_cell_ads_flag(ik) ) then recompute_all_cells_ads_flag = .true. else 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) endif enddo enddo ! end loop while recompute_all_cells_ads_flag (monolayer saturation) !? I'm not sure if this should be calculated here again. I have a feeling that ztot1 is meant ! as the value from the previous timestep !do ik = 1, nsoil ! ztot1(ig, ik) = porosity_prev(ig, ik) * znsoil(ig, ik) + ice(ig, ik) !enddo ! Calculation of Fnet + Fads ! ============================= ! calulate the flux variable ! calculate the flux in the top layer if (exchange(ig)) then ! if there is exchange ! calculate the flux taking diffusion to the atmosphere into account flux(ig, 0) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep & * ( znsoil(ig, 1) / rho(ig) - beta0(ig) - alpha0(ig) * znsoil(ig, 1) / rho(ig) ) elseif (pqsurf(ig).gt.0.) then ! assume that the surface is covered by waterice (if it is co2ice it should not call this subroutine at all) flux(ig, 0) = D_mid(ig, 1) / midlayer_dz(ig, 0) * ( znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig) ) endif ! check if the ground would take up water from the surface but there is none if ((.not.exchange(ig)).and.(pqsurf(ig).eq.0.).and.(flux(ig, 0).lt.0.)) then flux(ig, 0) = 0. endif ! calculate the flux in all other layers do ik = 1, nsoil - 1 flux(ig, ik) = D_inter(ig, ik) * ( znsoil(ig, ik + 1) - znsoil(ig, ik) ) / midlayer_dz(ig, ik) enddo ! calculate dztot1 which descibes the change in ztot1 (water vapor and ice). It is only used for output (directly and indirectly) ! calculate the change in ztot1 do ik = 1, nsoil - 1 dztot1(ik) = ( flux(ig, ik) - flux(ig, ik - 1) ) / interlayer_dz(ig, ik) & - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) & + B(ig, ik) / interlayer_dz(ig, ik) enddo dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) & - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) & + B(ig, nsoil) / interlayer_dz(ig, nsoil) ! Condensation / sublimation do ik = 1, nsoil ! loop over all layers ! print *, "ice in layer ", ik, ": ", ice(ig, ik) ! print *, "vapor in layer ", ik, ": ", znsoil(ig, ik) ! print *, "sat dens in layer ", ik, ": ", nsatsoil(ig, ik) ! print *, "" if (ice_index_flag(ik)) then ! if there is ice ! Compute ice content ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik) if (ice(ig, ik).lt.0.D0) then ! if all the ice is used up ! print *, "########complete sublimation in layer", ik, " cell:", ig ! print *, "ztot1: ", ztot1(ig, ik) ! print *, "dztot1*timestep: ", dztot1(ik) * ptimestep ! print *, "vapour: ", porosity(ig, ik) * nsatsoil(ig, ik) ! print *, "znsoil: ", znsoil(ig, ik) ! print *, "nsatsoil: ", nsatsoil(ig, ik) ! print *, "porosity: ", porosity(ig, ik) ! print *, "ice: ", ice(ig, ik) ! print *, "exchange: ", exchange(ig) ! print *, "co2ice: ", co2ice(ig) ! print *, "" ! print *, "zalpha: ", zalpha(ig, ik) ! print *, "zbeta: ", zbeta(ig, ik) ! print *, "" ! set the ice content to zero ice(ig, ik) = 0.D0 ! change the water vapor from the previous timestep. (Watch out! could go wrong) znsoilprev(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik) ! print *, "new znsoilprev", znsoilprev(ig, ik) ! remove the ice flag and raise the sublimation flag ice_index_flag(ik) = .false. sublimation_flag = .true. endif elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then ! if there is condenstation !ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep - porosity(ig, ik) * nsatsoil(ig, ik) ! print *, "=========== new condensation in layer", ik, " cell:", ig ! print *, "znsoil, nsatsoil: ", znsoil(ig, ik), nsatsoil(ig, ik) ! print *, "ztot1: ", ztot1(ig, ik) ! print *, "dztot1: ", dztot1(ik) ! print *, "ice: ", ice(ig,ik) ! print *, "" ! print *, "zalpha: ", zalpha(ig, ik) ! print *, "zbeta: ", zbeta(ig, ik) ! print *, "" !if (ice(ig, ik).lt.0.D0) then ! set the ice content to zero ! ice(ig, ik) = 0.D0 if (.not.condensation_flag(ik)) then ! raise the ice and sublimation flags ice_index_flag(ik) = .true. sublimation_flag = .true. ! print *, condensation_flag(ik) condensation_flag(ik) = .true. ! print *, condensation_flag(ik) ! print *, "set condensation flag" ! else ! print *, "condensation loop supressed" endif endif ! calculate the difference between total water content and ice + vapor content (only used for output) diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) & + adswater(ig, ik) - ztot1(ig, ik) - dztot1(ik) * ptimestep enddo ! loop over all layer enddo ! end first loop while sublimation_flag (condensation / sublimation) if (exchange(ig)) then ! if there is exchange ! calculate the temporty mixing ratio above the surface zq1temp2(ig) = beta0(ig) + alpha0(ig) * znsoil(ig, 1) / rho(ig) ! calculate the flux from the subsurface zdqsdifrego(ig) = porosity_ice_free(ig, 1) * pb(ig, 1) / ptimestep * (zq1temp2(ig) - znsoil(ig, 1) / rho(ig) ) else ! calculate the flux from the subsurface zdqsdifrego(ig) = D_mid(ig, 1) / midlayer_dz(ig, 0) * (znsoil(ig, 1) - qsat_surf(ig) * rhosurf(ig)) ! faut - il faire intervenir la porosite ? if ((pqsurf(ig).eq.0.).and.(zdqsdifrego(ig).lt.0.)) then zdqsdifrego(ig) = 0. endif endif ! Special case where there is not enough ice on the surface to supply the subsurface for the whole timestep ! (when exchange with the atmosphere is disabled): the upper boundary condition becomes a flux condition ! (and not qsat_surf) and all the remaining surface ice is sublimated and transferred to the subsurface. ! the actual change in surface ice happens in a higher subroutine ! ========================================================================================================= = if ( (.not.exchange(ig)) & .and. ( (-zdqsdifrego(ig) * ptimestep) & .gt.( pqsurf(ig) + pdqsdifpot(ig) * ptimestep) ) & .and.( (pqsurf(ig) + pdqsdifpot(ig) * ptimestep).gt.0. ) ) then ! calculate a new flux from the subsurface zdqsdifrego(ig) = -( pqsurf(ig) + pdqsdifpot(ig) * ptimestep ) / ptimestep ! ! check case where there is CO2 ice on the surface but qsurf = 0 ! if (co2ice(ig).gt.0.) then ! zdqsdifrego(ig) = 0.D0 ! endif do ik = 1, nsoil ! calculate A and B coefficients ! modified 2020 if ( .not.over_mono_sat_flag(ig, ik) ) then ! Assume cell below the monolayer saturation ! calculate A and B coefficients (first segment of the bilinear function) A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) else ! Assume cell over the monolayer saturation ! added 2020 ! calculate A and B coefficients (second segment of the bilinear function) A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 endif ! initialize all flags for the loop ! initialize the ice ice(ig, ik) = iceprev(ig, ik) if (iceprev(ig, ik).eq.0.) then ice_index_flag(ik) = .false. else ice_index_flag(ik) = .true. endif enddo ! loop while there is new sublimation sublimation_flag = .true. sublim_n = 0 do while (sublimation_flag) ! print *, "special case sublimation loop: ", sublim_n sublim_n = sublim_n + 1 if (sublim_n.gt.100) then print *, "special case sublimation not converging in call ", n print *, "might as well stop now" stop endif ! reset the sublimation flag sublimation_flag = .false. ! loop until good values accounting for monolayer saturation are found recompute_all_cells_ads_flag = .true. do while (recompute_all_cells_ads_flag) ! reset loop flag recompute_all_cells_ads_flag = .false. ! calculate the Delta, Alpha, and Beta coefficients in the top layer zdelta(ig, 1) = A(ig, 1) + porosity(ig, 1) * C(ig, 1) + D_inter(ig, 1) / midlayer_dz(ig, 1) zalpha(ig, 1) = D_inter(ig, 1) / midlayer_dz(ig, 1) * 1.D0 / zdelta(ig, 1) zbeta(ig, 1) = ( porosity_prev(ig, 1) * C(ig, 1) * znsoilprev2(ig, 1) + B(ig, 1) - zdqsdifrego(ig) ) & / zdelta(ig, 1) ! check for ice if (ice_index_flag(1)) then ! set the alpha coefficient to zero zalpha(ig, 1) = 0.D0 ! set the beta coefficient to the number density at saturation pressure zbeta(ig, 1) = nsatsoil(ig, 1) endif do ik = 2, nsoil - 1 ! loop over all other layers ! calculate the Delta, Alpha, and Beta coefficients in the layer zdelta(ig, ik) = A(ig, ik) + porosity(ig, ik) * C(ig, ik) + D_inter(ig, ik) / midlayer_dz(ig, ik) & + D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * (1.D0 - zalpha(ig, ik - 1)) zalpha(ig, ik) = D_inter(ig, ik) / midlayer_dz(ig, ik) * 1.D0 / zdelta(ig, ik) zbeta(ig, ik) = ( D_inter(ig, ik - 1) / midlayer_dz(ig, ik - 1) * zbeta(ig, ik - 1) & + porosity_prev(ig, ik) * C(ig, ik) * znsoilprev2(ig, ik) + B(ig, ik) ) / zdelta(ig, ik) ! check for ice if (ice_index_flag(ik)) then ! set the alpha coefficient to zero zalpha(ig, ik) = 0.D0 ! set the beta coefficient to the number density at saturation pressure zbeta(ig, ik) = nsatsoil(ig, ik) endif enddo ! calculate the preliminary amount of water vapor in the bottom layer if (ice_index_flag(nsoil)) then ! set the vapor number density to saturation znsoil(ig, nsoil) = nsatsoil(ig, nsoil) else ! calculate the vapor number density in the lowest layer znsoil(ig, nsoil) = ( D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * zbeta(ig, nsoil - 1) & + porosity_prev(ig, nsoil) * C(ig, nsoil) * znsoilprev2(ig, nsoil) + B(ig, nsoil) ) & / ( porosity(ig, nsoil) * C(ig, nsoil) & + D_inter(ig, nsoil - 1) / midlayer_dz(ig, nsoil - 1) * (1.D0 - zalpha(ig, nsoil - 1)) & + A(ig, nsoil) ) endif ! calculate the preliminary amount of adsorbed water if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024 adswater_temp(ig, nsoil) = ( Ka(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) ) & / (1.D0 + ptimestep * Kd(ig, nsoil)) else adswater_temp(ig, nsoil) = ( Ka2(ig, nsoil) * ptimestep * znsoil(ig, nsoil) + adswprev(ig, nsoil) + ptimestep * c0(ig,nsoil) * Kd2(ig,nsoil)) & / (1.D0 + ptimestep * Kd2(ig, nsoil)) endif do ik = nsoil-1, 1, -1 ! backward loop over all layers znsoil(ig, ik) = zalpha(ig, ik) * znsoil(ig, ik + 1) + zbeta(ig, ik) if (.not.over_mono_sat_flag(ig, nsoil)) then ! modified 2024 adswater_temp(ig, ik) = ( Ka(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) ) & / (1.D0 + ptimestep * Kd(ig, ik)) else adswater_temp(ig, ik) = ( Ka2(ig, ik) * ptimestep * znsoil(ig, ik) + adswprev(ig, ik) + ptimestep * c0(ig,ik) * Kd2(ig,ik)) & / (1.D0 + ptimestep * Kd2(ig, ik)) endif enddo ! check for any saturation do ik = 1, nsoil ! Change of coefficients ! ARNAU recompute_cell_ads_flag(ik) = .false. ! Assume there is no change of coefficients ! ARNAU if ( adswater_temp(ig, ik).gt.adswater_sat_mono ) then print *, "NOTICE: over monolayer saturation" if ( .not.over_mono_sat_flag(ig, ik) ) then ! If the cell enters this scope, it ! means that the cell is over the monolayer ! saturation after calculation, but the ! coefficients assume it is below. Therefore, ! the cell needs to be recomputed over_mono_sat_flag(ig, ik) = .true. recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients ! 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) A(ig, ik) = E2(ig, ik) * interlayer_dz(ig, ik) B(ig, ik) = interlayer_dz(ig, ik) * F2(ig, ik) * adswprev(ig, ik) & + interlayer_dz(ig, ik) * (F2(ig, ik) * Kd2(ig, ik) * c0(ig, ik) * ptimestep - Kd2(ig,ik)*c0(ig,ik)) ! Corrected 2024 endif else if ( over_mono_sat_flag(ig, ik) ) then ! If the cell enters this scope, it ! means that the cell is below the monolayer ! saturation after calculation, but the ! coefficients assume it is above. Therefore, ! the cell needs to be recomputed over_mono_sat_flag(ig, ik) = .false. recompute_cell_ads_flag(ik) = .true. ! There is a change of coefficients ! calculate A and B coefficients (reset to first segment of the bilinear function) A(ig, ik) = E(ig, ik) * interlayer_dz(ig, ik) B(ig, ik) = interlayer_dz(ig, ik) * F(ig, ik) * adswprev(ig, ik) endif endif enddo ! raise the saturation flag if any layer has saturated and reset the first layer saturation flag do ik = 1, nsoil ! modified 2020 if ( recompute_cell_ads_flag(ik) ) then recompute_all_cells_ads_flag = .true. else 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) endif enddo enddo ! end while loop (adsorption saturation) ! set the flux to the previously calculated value for the top layer flux(ig, 0) = zdqsdifrego(ig) ! calculate the flux for all other layers do ik = 1, nsoil - 1 flux(ig, ik) = D_inter(ig, ik) * (znsoil(ig, ik + 1) - znsoil(ig, ik)) / midlayer_dz(ig, ik) enddo ! calculate the change in ztot1 do ik = 1, nsoil - 1 dztot1(ik) = (flux(ig, ik) - flux(ig, ik - 1)) / interlayer_dz(ig, ik) & - A(ig, ik) / interlayer_dz(ig, ik) * znsoil(ig, ik) & + B(ig, ik) / interlayer_dz(ig, ik) enddo dztot1(nsoil) = -flux(ig, nsoil - 1) / interlayer_dz(ig, nsoil) & - A(ig, nsoil) / interlayer_dz(ig, nsoil) * znsoil(ig, nsoil) & + B(ig, nsoil) / interlayer_dz(ig, nsoil) ! Condensation / sublimation do ik = 1, nsoil if (ice_index_flag(ik)) then ! Compute ice content ice(ig, ik) = ztot1(ig, ik) + dztot1(ik) * ptimestep & - porosity(ig, ik) * nsatsoil(ig, ik) if (ice(ig, ik).lt.0.D0) then ! if all the ice is used up ! print *, "############complete sublimation in layer ", ik ! print *, "ztot1: ", ztot1(ig, ik) ! print *, "dztot1*timestep: ", dztot1(ik) * ptimestep ! print *, "vapour: ", porosity(ig, ik) * nsatsoil(ig, ik) ! print *, "znsoil: ", znsoil(ig, ik) ! print *, "nsatsoil: ", nsatsoil(ig, ik) ! print *, "porosity: ", porosity(ig, ik) ! print *, "ice: ", ice(ig, ik) ! print *, "exchange: ", exchange(ig) ! print *, "co2ice: ", co2ice(ig) ! print *, "" ! print *, "zalpha: ", zalpha(ig, ik) ! print *, "zbeta: ", zbeta(ig, ik) ! print *, "" ! set the ice content to zero ice(ig, ik) = 0.D0 if (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then print *, "WARNING! complete sublimation, but vapor is oversaturated" print *, "special case loop in cell", ig, ik ,"will be supressed" else ! change the water vapor from the previous timestep. (Watch out! could go wrong) znsoilprev2(ig, ik) = ztot1(ig, ik) / porosity_prev(ig, ik) ! print *, "new znsoilprev", znsoilprev(ig, ik) ! remove the ice flag and raise the sublimation flag ice_index_flag(ik) = .false. sublimation_flag = .true. endif endif elseif (znsoil(ig, ik).gt.nsatsoil(ig, ik)) then ! if there is new ice through condensation if (.not.condensation_flag(ik)) then ! raise the ice and sublimation flags ice_index_flag(ik) = .true. sublimation_flag = .true. ! print *, condensation_flag(ik) condensation_flag(ik) = .true. ! print *, condensation_flag(ik) ! print *, "set condensation flag" ! else ! print *, "special case condensation loop supressed" endif endif enddo enddo ! end first loop while (condensation / sublimation) endif ! Special case for all ice on the surface is used up do ik = 1, nsoil diff(ig, ik) = ice(ig, ik) + porosity(ig, ik) * znsoil(ig, ik) & ! only used for output + adswater(ig, ik) - adswprev(ig, ik) & - ztot1(ig, ik) - dztot1(ik) * ptimestep !? This is inconsistent, in the not special case the previous adsorbed water is not counted !? also this just overwrites the previous calculation if I see that correctly ! calculate the total amount of water ztot1(ig, ik) = porosity(ig, ik) * znsoil(ig, ik) + ice(ig, ik) h2otot(ig, ik) = adswater(ig, ik) + ztot1(ig, ik) enddo endif ! if there is no polar cap enddo ! loop over all gridpoints ! check for choking condition do ig = 1, ngrid if(.not.watercaptag(ig)) then do ik = 1, nsoil - 1 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 if ( flux(ig, ik - 1).gt.0.) then if (.not.close_top(ig, ik).and.print_closure_warnings) then print *, "NOTICE: closing top of soil layer due to inward flux", ig, ik endif close_top(ig, ik) = .true. elseif (flux(ig, ik - 1).lt.0.) then if (close_top(ig, ik).and.print_closure_warnings) then print *, "NOTICE: opening top of soil layer due to outward flux", ig, ik endif close_top(ig, ik) = .false. endif if ( flux(ig, ik).lt.0.) then if (.not.close_bottom(ig, ik).and.print_closure_warnings) then print *, "NOTICE: closing bottom of soil layer due to inward flux", ig, ik endif close_bottom(ig, ik) = .true. elseif (flux(ig, ik).gt.0.) then if (close_bottom(ig, ik).and.print_closure_warnings) then print *, "NOTICE: opening bottom of soil layer due to outward flux", ig, ik endif close_bottom(ig, ik) = .false. endif ! if ( (flux(ig, ik).lt.flux(ig, ik - 1)).and.(flux(ig, ik).gt.0D0) ) then ! if (.not.close_top(ig, ik).and.print_closure_warnings) then ! print *, "NOTICE: closing top of soil layer due to ice", ig, ik ! endif ! close_top(ig, ik) = .true. ! ! elseif ( (flux(ig, ik).lt.flux(ig, ik - 1)).and.(flux(ig, ik - 1).lt.0D0) ) then ! if (.not.close_bottom(ig, ik).and.print_closure_warnings) then ! print *, "NOTICE: closing bottom of soil layer due to ice", ig, ik ! endif ! close_bottom(ig, ik) = .true. ! endif ! ! if (close_top(ig, ik).and.close_bottom(ig, ik)) then ! close_top(ig, ik) = .false. ! close_bottom(ig, ik) = .false. ! if (print_closure_warnings) then ! print *, "WARNING: Reopening soil layer after complete closure:", ig, ik ! endif ! elseif (close_top(ig, ik).or.close_bottom(ig, ik)) then ! if ((flux(ig, ik) - flux(ig, ik - 1)).gt.0D0) then ! close_top(ig, ik) = .false. ! close_bottom(ig, ik) = .false. ! if (print_closure_warnings) then ! print *, "WARNING: Reopening soil layer due to rising ice:", ig, ik ! endif ! endif ! ! if (close_top(ig, ik).and.(flux(ig, ik - 1).lt.0D0)) then ! close_top = .false. ! if (print_closure_warnings) then ! print *, "NOTICE: Reopening top of soil layer due to upward flux:", ig, ik ! endif ! endif ! ! if (close_bottom(ig, ik).and.(flux(ig, ik).gt.0D0)) then ! close_bottom = .false. ! if (print_closure_warnings) then ! print *, "NOTICE: Reopening bottom of soil layer due to downward flux:", ig, ik ! endif ! endif ! endif else ! opening condition that ice content has fallen sufficiently if (close_top(ig, ik).or.close_bottom(ig, ik).and.print_closure_warnings) then print *, "NOTICE: Reopening soillayer due to decrease in ice", ig, ik endif close_top(ig, ik) = .false. close_bottom(ig, ik) = .false. endif enddo endif enddo ! Computation of total mass ! =========================== do ig = 1, ngrid ! initialize the masses to zero mass_h2o(ig) = 0.D0 mass_ice(ig) = 0.D0 if(.not.watercaptag(ig)) then do ik = 1, nsoil ! calculate the total water and ice mass mass_h2o(ig) = mass_h2o(ig) + h2otot(ig, ik) * interlayer_dz(ig, ik) mass_ice(ig) = mass_ice(ig) + ice(ig, ik) * interlayer_dz(ig, ik) ! calculate how close the water vapor content is to saturizing the adsorbed water if(choice_ads.ne.0) preduite(ig, ik) = znsoil(ig, ik) / nsat(ig, ik) ! write the results to the return variable pqsoil(ig, ik, igcm_h2o_vap_soil) = znsoil(ig, ik) pqsoil(ig, ik, igcm_h2o_ice_soil) = ice(ig, ik) pqsoil(ig, ik, igcm_h2o_vap_ads) = adswater(ig, ik) if (close_top(ig, ik)) then close_out(ig, ik) = 1 elseif (close_bottom(ig, ik)) then close_out(ig, ik) = -1 else close_out(ig, ik) = 0 endif enddo endif if (watercaptag(ig)) then do ik = 1, nsoil saturation_water_ice(ig, ik) = -1 enddo endif ! convert the logical :: flag exchange to a numeric output if (watercaptag(ig)) then exch(ig) = -2. elseif (exchange(ig)) then exch(ig) = 1. else exch(ig) = -1. endif enddo ! calculate the global total value mass_h2o_tot = 0.D0 mass_ice_tot = 0.D0 do ig = 1, ngrid mass_h2o_tot = mass_h2o_tot + mass_h2o(ig) * cell_area(ig) mass_ice_tot = mass_ice_tot + mass_ice(ig) * cell_area(ig) enddo ! print and iterate the call number ! print * , 'Subsurface call n = ', n n = n + 1 ! ----------------------------------------------------------------------- ! Output in diagfi and diagsoil ! ----------------------------------------------------------------------- ! reformat flux, because it has a unusual shape do ig = 1, ngrid nsurf(ig) = rhosurf(ig) * pq(ig, 1, igcm_h2o_vap) do ik = 1, nsoil - 1 var_flux_soil(ig, ik) = flux(ig, ik) enddo var_flux_soil(ig, nsoil) = 0. enddo if(writeoutput) then !print *, "flux shape: ", shape(flux) !print *, "adswater shape ", shape(adswater) !print *, "ngrid= ", ngrid !if (ngrid.eq.1) then ! if in 1D ! print *, "writing 1D data" ! print *, real(adswater(:, :), 4) ! print *, shape(adswater(:, :)) zalpha(1, nsoil) = 0 zbeta(1, nsoil) = 0 ! call write_output("zalpha","diffusion alpha coefficient", "unknown",real(zalpha(1, :))) ! call write_output("zbeta","diffusion beta coefficient", "unknown",real(zbeta(1, :))) ! call write_output("adswater","adswater", "kg / m^3",real(adswater(1, :))) ! call write_output("znsoil","znsoil", "kg / m^3",real(znsoil(1, :))) ! call write_output("ice","ice", "kg / m^3",real(ice(1, :))) ! call write_output("h2otot","total water content", "kg / m^3",real(h2otot(1, :))) ! call write_output("flux_soil","flux_soil", "kg / m^2",real(flux(1, :))) ! call write_output("flux","flux soil", "kg / m^2",real(zdqsdifrego(:))) ! call write_output("flux_surf","surface flux", "kg / m^2",real(zdqsdifrego(1))) ! call write_output("exchange","exchange", "boolean",real(exch(1))) !else ! print *, mass_h2o_tot ! print *, real(mass_h2o_tot, 4) ! call write_output("dztot1","change in dztot", "unknown",real(dztot1)) call write_output("flux_soillayer","flux of water between the soil layers", "kg m-2 s-1",var_flux_soil) ! call write_output("A_coef","A coefficient of dztot1", "unknown",real(B)) ! call write_output("B_coef","B coefficient of dztot1", "unknown",real(B)) ! call write_output("H2O_init","initialized H2O", "kg/m^2",real(H2O)) ! call write_output("H2O_depth_init","initialized H2O depth ", "m",real(H2O_depth)) call write_output("ice_saturation_soil","Water ice saturation in the soil layers", "Percent",saturation_water_ice) ! call write_output("mass_h2o_tot","total mass of subsurface water over the planet", "kg",real(mass_h2o_tot)) ! call write_output("mass_ice_tot","total mass of subsurface ice over the planet", "kg",real(mass_ice_tot)) call write_output("mass_h2o_soil","Mass of subsurface water column at each point", "kg m-2",(mass_h2o)) call write_output("mass_ice_soil","Mass of subsurface ice at each point", "kg m-2",(mass_ice)) call write_output("znsoil","Water vapor soil concentration ", "kg.m - 3 of pore air",znsoil) call write_output("nsatsoil","subsurface water vapor saturation density", "kg/m^3",real(nsatsoil)) call write_output("nsurf","surface water vapor density", "kg/m^3",real(nsurf)) call write_output("adswater","subsurface adsorbed water", "kg/m^3",real(adswater)) ! call write_output("subsurfice","subsurface ice", "kg/m^3",real(ice)) call write_output("flux_rego",'flux of water from the regolith','kg/m^2',zdqsdifrego) ! call write_output("exchange",'exchange','no unit',real(exch)) ! call write_output("close",'close top, bottom, or none (1, -1, 0)','no unit',real(close_out)) ! call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_inter) !endif endif RETURN END