source: trunk/LMDZ.MARS/libf/phymars/soilwater.F90

Last change on this file was 3726, checked in by emillour, 3 months ago

Mars PCM:
Turn "callkeys.h" into module "callkeys_mod.F90"
EM

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