diff --git a/CHANGELOG.md b/CHANGELOG.md index 38a972a5..17f70c5d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,13 @@ [#234](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/234) - Save water stress factor and water potential into csv files. [#229](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/229) +- Add an option to define soil layer settings through a csv file discussed in + [#237](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/237) and + [#241](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/241) and added in + [#243](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/243) +- Remove unnecessary function discussed in + [#244](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/244) and removed in + [#243](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/243) **Fixed:** @@ -22,8 +29,7 @@ - Defining the indices of the first four layers discussed in [#237](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/237) and fixed in [#238](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/238) - - + # [1.5.0](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.5.0) - 3 Jan 2024 diff --git a/docs/STEMMUS_SCOPE_on_CRIB.md b/docs/STEMMUS_SCOPE_on_CRIB.md index 0ff3628d..b35d0668 100644 --- a/docs/STEMMUS_SCOPE_on_CRIB.md +++ b/docs/STEMMUS_SCOPE_on_CRIB.md @@ -34,6 +34,14 @@ provides parameter inputs for PROSPECT, leaf_biochemical, fluorescence, soil, canopy, aerodynamic, angles, photosynthetic temperature dependence functional parameters, etc. + - input_soilLayThick.csv (optional): A file to change the discretization of + the soil layers of the STEMMUS model. An example of this file is in + [example_data folder](../example_data). This file (if needed) should be copied into the + `InputPath` folder. If this file is used, it will override the default settings of + the soil layers. The file has three columns: 1) layer number, 2) layer thickness, + and 3) maximum root depth. The user is free to change the values of the three columns. + Also, the number of rows determines the number of the soil layers and the total + thickness of the soil column (sum of soil layer thickness). 2. Config file: it is a text file that sets the paths **required** by the model. For example, see [config_file_crib.txt](../config_file_crib.txt) in this @@ -124,4 +132,4 @@ main script using MATLAB command line in a terminal: ```bash matlab -nodisplay -nosplash -nodesktop -r "run('STEMMUS_SCOPE.m');exit;" -``` \ No newline at end of file +``` diff --git a/docs/STEMMUS_SCOPE_on_Snellius.md b/docs/STEMMUS_SCOPE_on_Snellius.md index 3106037b..19525286 100644 --- a/docs/STEMMUS_SCOPE_on_Snellius.md +++ b/docs/STEMMUS_SCOPE_on_Snellius.md @@ -17,6 +17,7 @@ Dutch National supercomputer hosted at SURF. - radiationdata - soil_spectra - input_data.xlsx + - input_soilThick.csv (optional) For the explanation of the directories see [Dataflow of STEMMUS_SCOPE on CRIB](./STEMMUS_SCOPE_on_CRIB.md#dataflow-of-stemmus_scope-on-crib). @@ -150,4 +151,4 @@ mkdir -p slurm sbatch run_stemmus_scope_snellius.sh ``` -This creates a log file per each forcing file in the folder `slurm`. \ No newline at end of file +This creates a log file per each forcing file in the folder `slurm`. diff --git a/example_data/input_soilLayThick.csv b/example_data/input_soilLayThick.csv new file mode 100644 index 00000000..bca17029 --- /dev/null +++ b/example_data/input_soilLayThick.csv @@ -0,0 +1,61 @@ +NL,LayThick (cm),RootDepth (cm) +1,1,450 +2,1, +3,1, +4,2, +5,2, +6,2, +7,2, +8,2, +9,2, +10,2, +11,2, +12,2, +13,2, +14,2, +15,2.5, +16,2.5, +17,2.5, +18,2.5, +19,5, +20,5, +21,5, +22,5, +23,5, +24,10, +25,10, +26,10, +27,10, +28,10, +29,10, +30,10, +31,10, +32,10, +33,10, +34,10, +35,10, +36,10, +37,10, +38,10, +39,10, +40,10, +41,15, +42,15, +43,20, +44,20, +45,20, +46,20, +47,20, +48,20, +49,20, +50,20, +51,20, +52,20, +53,20, +54,20, +55,25, +56,25, +57,25, +58,25, +59,25, +60,25, diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE index 0cdc650c..8b070199 100755 Binary files a/run_model_on_snellius/exe/STEMMUS_SCOPE and b/run_model_on_snellius/exe/STEMMUS_SCOPE differ diff --git a/src/+conductivity/+hydraulicConductivity/calcuulateDTheta_LLh.m b/src/+conductivity/+hydraulicConductivity/calculateDTheta_LLh.m similarity index 93% rename from src/+conductivity/+hydraulicConductivity/calcuulateDTheta_LLh.m rename to src/+conductivity/+hydraulicConductivity/calculateDTheta_LLh.m index eaf543c6..6208e678 100644 --- a/src/+conductivity/+hydraulicConductivity/calcuulateDTheta_LLh.m +++ b/src/+conductivity/+hydraulicConductivity/calculateDTheta_LLh.m @@ -1,4 +1,4 @@ -function dtheta_llh = calcuulateDTheta_LLh(dtheta_uuh, theta_m, theta_uu, theta_ll, gamma_hh, SoilVariables, VanGenuchten) +function dtheta_llh = calculateDTheta_LLh(dtheta_uuh, theta_m, theta_uu, theta_ll, gamma_hh, SoilVariables, VanGenuchten, ModelSettings) theta_s = VanGenuchten.Theta_s; theta_r = VanGenuchten.Theta_r; alpha = VanGenuchten.Alpha; @@ -14,8 +14,6 @@ theta_l = SoilVariables.Theta_L; h_frez = SoilVariables.h_frez; - % get model settings - ModelSettings = io.getModelSettings(); heatTerm = hh + hh_frez; if ModelSettings.SWCC == 1 if ModelSettings.SFCC == 1 diff --git a/src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m b/src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m index ef5514dd..777df1f7 100644 --- a/src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m +++ b/src/+conductivity/+hydraulicConductivity/calculateDTheta_UUh.m @@ -1,4 +1,4 @@ -function dtheta_uuh = calculateDTheta_UUh(theta_uu, theta_m, theta_ll, gamma_hh, SoilVariables, VanGenuchten) +function dtheta_uuh = calculateDTheta_UUh(theta_uu, theta_m, theta_ll, gamma_hh, SoilVariables, VanGenuchten, ModelSettings) theta_s = VanGenuchten.Theta_s; theta_r = VanGenuchten.Theta_r; alpha = VanGenuchten.Alpha; @@ -14,9 +14,6 @@ phi_s = SoilVariables.Phi_s; lamda = SoilVariables.Lamda; - % get model settings - ModelSettings = io.getModelSettings(); - heatTerm = hh + hh_frez; if ModelSettings.SWCC == 1 if ModelSettings.SFCC == 1 diff --git a/src/+conductivity/+hydraulicConductivity/calculateSe.m b/src/+conductivity/+hydraulicConductivity/calculateSe.m index 08a98e9a..255963df 100644 --- a/src/+conductivity/+hydraulicConductivity/calculateSe.m +++ b/src/+conductivity/+hydraulicConductivity/calculateSe.m @@ -1,6 +1,4 @@ -function se = calculateSe(theta_ll, gamma_hh, SoilVariables) - % get model settings - ModelSettings = io.getModelSettings(); +function se = calculateSe(theta_ll, gamma_hh, SoilVariables, ModelSettings) % get soil constants SoilConstants = io.getSoilConstants(); diff --git a/src/+conductivity/+hydraulicConductivity/calculateTheta_II.m b/src/+conductivity/+hydraulicConductivity/calculateTheta_II.m index 125f191c..4e3a4fc8 100644 --- a/src/+conductivity/+hydraulicConductivity/calculateTheta_II.m +++ b/src/+conductivity/+hydraulicConductivity/calculateTheta_II.m @@ -1,7 +1,4 @@ -function theta_ii = calculateTheta_II(tt, xcap, hh, Theta_II) - - % get model settings - ModelSettings = io.getModelSettings(); +function theta_ii = calculateTheta_II(tt, xcap, hh, Theta_II, ModelSettings) Tf1 = 273.15 + 1; Tf2 = 273.15 - 3; diff --git a/src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m b/src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m index a580dfbd..6c5aa9dd 100644 --- a/src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m +++ b/src/+conductivity/+hydraulicConductivity/calculateTheta_LL.m @@ -1,6 +1,4 @@ -function theta_ll = calculateTheta_LL(theta_uu, theta_ii, theta_m, gamma_hh, SoilVariables, VanGenuchten) - % get model settings - ModelSettings = io.getModelSettings(); +function theta_ll = calculateTheta_LL(theta_uu, theta_ii, theta_m, gamma_hh, SoilVariables, VanGenuchten, ModelSettings) % load Constants Constants = io.define_constants(); diff --git a/src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m b/src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m index 8d4348a3..db468611 100644 --- a/src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m +++ b/src/+conductivity/+hydraulicConductivity/calculateTheta_UU.m @@ -1,4 +1,4 @@ -function theta_uu = calculateTheta_UU(theta_m, gamma_hh, SoilVariables, VanGenuchten) +function theta_uu = calculateTheta_UU(theta_m, gamma_hh, SoilVariables, VanGenuchten, ModelSettings) hh = SoilVariables.hh; phi_s = SoilVariables.Phi_s; @@ -10,9 +10,6 @@ n = VanGenuchten.n; m = VanGenuchten.m; - % get model settings - ModelSettings = io.getModelSettings(); - % calculate theta_uu if ModelSettings.SWCC == 1 if ModelSettings.SFCC == 1 diff --git a/src/+conductivity/+hydraulicConductivity/fixHeat.m b/src/+conductivity/+hydraulicConductivity/fixHeat.m index 4dcb8d97..f8626a56 100644 --- a/src/+conductivity/+hydraulicConductivity/fixHeat.m +++ b/src/+conductivity/+hydraulicConductivity/fixHeat.m @@ -1,7 +1,4 @@ -function [hh, hh_frez] = fixHeat(hh, hh_frez, Phi_s) - - % get model settings - ModelSettings = io.getModelSettings(); +function [hh, hh_frez] = fixHeat(hh, hh_frez, Phi_s, ModelSettings) if ModelSettings.SWCC == 1 if ModelSettings.SFCC ~= 1 diff --git a/src/+conductivity/calculateGasConductivity.m b/src/+conductivity/calculateGasConductivity.m index ccca06b3..0d7d3b82 100644 --- a/src/+conductivity/calculateGasConductivity.m +++ b/src/+conductivity/calculateGasConductivity.m @@ -1,4 +1,4 @@ -function k_g = calculateGasConductivity(InitialValues, TransportCoefficient, VanGenuchten, SoilVariables) +function k_g = calculateGasConductivity(InitialValues, TransportCoefficient, VanGenuchten, SoilVariables, ModelSettings) %{ This is to calculate the intrinsic permeability of soil for gas flow. Scanlon, B. R. (2000), Soil gas movement in unsaturated systems, in @@ -10,9 +10,6 @@ 20107, doi:10.1029/2011JD015835, 2011. %} - % get model settings - ModelSettings = io.getModelSettings(); - % load Constants Constants = io.define_constants(); diff --git a/src/+conductivity/calculateGasDispersivity.m b/src/+conductivity/calculateGasDispersivity.m index 4bbea955..b7d483f5 100644 --- a/src/+conductivity/calculateGasDispersivity.m +++ b/src/+conductivity/calculateGasDispersivity.m @@ -1,4 +1,4 @@ -function GasDispersivity = calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g) +function GasDispersivity = calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g, ModelSettings) %{ This is to calculate the gas phase longitudinal dispersivity. Zeng, Y., Su, Z., Wan, L. and Wen, i.: Numerical analysis of @@ -7,9 +7,6 @@ 20107, doi:10.1029/2011JD015835, 2011. %} - % get model settings - ModelSettings = io.getModelSettings(); - % get model constants Constants = io.define_constants(); diff --git a/src/+conductivity/calculateHydraulicConductivity.m b/src/+conductivity/calculateHydraulicConductivity.m index b6c8d2c4..2ed12a85 100644 --- a/src/+conductivity/calculateHydraulicConductivity.m +++ b/src/+conductivity/calculateHydraulicConductivity.m @@ -1,4 +1,4 @@ -function SoilVariables = calculateHydraulicConductivity(SoilVariables, VanGenuchten, KIT, L_f) +function SoilVariables = calculateHydraulicConductivity(SoilVariables, VanGenuchten, KIT, L_f, ModelSettings) %{ This is to calculate the hydraulic conductivity of soil, based on hydraulic conductivity models (like VG and others). @@ -11,9 +11,6 @@ hydraulic conductivity models (like VG and others). 20107, doi:10.1029/2011JD015835, 2011. %} - % get model settings - ModelSettings = io.getModelSettings(); - % load Constants Constants = io.define_constants(); @@ -73,26 +70,26 @@ hydraulic conductivity models (like VG and others). SV = sliceVector(SV, lengthX, MN); SV = sliceMatrix(SV, lengthX, ModelSettings.nD, MN, j); - [hh, hh_frez] = conductivity.hydraulicConductivity.fixHeat(SV.hh, SV.hh_frez, SV.Phi_s); + [hh, hh_frez] = conductivity.hydraulicConductivity.fixHeat(SV.hh, SV.hh_frez, SV.Phi_s, ModelSettings); SV.hh = hh; SV.hh_frez = hh_frez; Gamma_hh = conductivity.hydraulicConductivity.calculateGamma_hh(SV.hh); Theta_m = conductivity.hydraulicConductivity.calculateTheta_m(Gamma_hh, VG, SV.POR); - Theta_UU = conductivity.hydraulicConductivity.calculateTheta_UU(Theta_m, Gamma_hh, SV, VG); + Theta_UU = conductivity.hydraulicConductivity.calculateTheta_UU(Theta_m, Gamma_hh, SV, VG, ModelSettings); % circular calculation of Theta_II! See issue 181, item 3 % Theta_II is soil ice content, % Theta_LL is liquid water content, % Theta_UU is the total water content before soil freezing. The % 'Theta_UU' is set as saturation. - Theta_II = conductivity.hydraulicConductivity.calculateTheta_II(SV.TT, SV.XCAP, SV.hh, SV.Theta_II); - Theta_LL = conductivity.hydraulicConductivity.calculateTheta_LL(Theta_UU, Theta_II, Theta_m, Gamma_hh, SV, VG); + Theta_II = conductivity.hydraulicConductivity.calculateTheta_II(SV.TT, SV.XCAP, SV.hh, SV.Theta_II, ModelSettings); + Theta_LL = conductivity.hydraulicConductivity.calculateTheta_LL(Theta_UU, Theta_II, Theta_m, Gamma_hh, SV, VG, ModelSettings); Theta_II = (Theta_UU - Theta_LL) * Constants.RHOL / Constants.RHOI; % ice water contentTheta_II - DTheta_UUh = conductivity.hydraulicConductivity.calculateDTheta_UUh(Theta_UU, Theta_m, Theta_LL, Gamma_hh, SV, VG); - DTheta_LLh = conductivity.hydraulicConductivity.calcuulateDTheta_LLh(DTheta_UUh, Theta_m, Theta_UU, Theta_LL, Gamma_hh, SV, VG); - Se = conductivity.hydraulicConductivity.calculateSe(Theta_LL, Gamma_hh, SV); + DTheta_UUh = conductivity.hydraulicConductivity.calculateDTheta_UUh(Theta_UU, Theta_m, Theta_LL, Gamma_hh, SV, VG, ModelSettings); + DTheta_LLh = conductivity.hydraulicConductivity.calculateDTheta_LLh(DTheta_UUh, Theta_m, Theta_UU, Theta_LL, Gamma_hh, SV, VG, ModelSettings); + Se = conductivity.hydraulicConductivity.calculateSe(Theta_LL, Gamma_hh, SV, ModelSettings); % Ratio_ice used in Condg_k_g.m if Theta_UU ~= 0 diff --git a/src/+conductivity/calculateSoilThermalProperites.m b/src/+conductivity/calculateSoilThermalProperites.m index b6c35006..2162fb44 100644 --- a/src/+conductivity/calculateSoilThermalProperites.m +++ b/src/+conductivity/calculateSoilThermalProperites.m @@ -1,4 +1,4 @@ -function [ETCON, EHCAP, TETCON, EfTCON, ZETA] = calculateSoilThermalProperites(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV) +function [ETCON, EHCAP, TETCON, EfTCON, ZETA] = calculateSoilThermalProperites(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, ModelSettings, DRHOVT, L, RHOV) %{ This is used to calculate Heat Capacity and Thermal Conductivity. %} @@ -38,9 +38,6 @@ % load Constants Constants = io.define_constants(); - % get model settings - ModelSettings = io.getModelSettings(); - MN = 0; TCON(1) = 1.37e-3 * 4.182; for i = 1:ModelSettings.NL diff --git a/src/+conductivity/calculateThermalConductivityCapacity.m b/src/+conductivity/calculateThermalConductivityCapacity.m index cd9bf097..80d9bcd0 100644 --- a/src/+conductivity/calculateThermalConductivityCapacity.m +++ b/src/+conductivity/calculateThermalConductivityCapacity.m @@ -1,4 +1,4 @@ -function ThermalConductivityCapacity = calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV) +function ThermalConductivityCapacity = calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, ModelSettings, DRHOVT, L, RHOV) %{ This is to calculate thermal conductivity and thermal capacity. Chung, S. O., and R. Horton (1987), Soil heat and water flow with a @@ -9,9 +9,6 @@ 20107, doi:10.1029/2011JD015835, 2011. %} - % get model settings - ModelSettings = io.getModelSettings(); - % load Constants Constants = io.define_constants(); @@ -21,7 +18,7 @@ RHO_bulk = ThermalConductivity.RHO_bulk; if ModelSettings.ThmrlCondCap == 1 - [ETCON, EHCAP, TETCON, EfTCON, ZETA] = conductivity.calculateSoilThermalProperites(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV); + [ETCON, EHCAP, TETCON, EfTCON, ZETA] = conductivity.calculateSoilThermalProperites(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, ModelSettings, DRHOVT, L, RHOV); end for i = 1:ModelSettings.NL diff --git a/src/+conductivity/calculateTransportCoefficient.m b/src/+conductivity/calculateTransportCoefficient.m index 17b399d6..945c52d6 100644 --- a/src/+conductivity/calculateTransportCoefficient.m +++ b/src/+conductivity/calculateTransportCoefficient.m @@ -1,4 +1,4 @@ -function TransportCoefficient = calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t) +function TransportCoefficient = calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t, ModelSettings) %{ This is to calculate the transport coefficient for absorbed liquid flow due to temperature gradient. @@ -20,9 +20,6 @@ TransportCoefficient.MU_W = InitialValues.MU_W; TransportCoefficient.D_Ta = InitialValues.D_Ta; - % get model settings - ModelSettings = io.getModelSettings(); - % load Constants Constants = io.define_constants(); diff --git a/src/+conductivity/calculateVaporVariables.m b/src/+conductivity/calculateVaporVariables.m index 90254610..e4aa24a3 100644 --- a/src/+conductivity/calculateVaporVariables.m +++ b/src/+conductivity/calculateVaporVariables.m @@ -1,4 +1,4 @@ -function VaporVariables = calculateVaporVariables(InitialValues, SoilVariables, VanGenuchten, ThermalConductivityCapacity, TT) +function VaporVariables = calculateVaporVariables(InitialValues, SoilVariables, VanGenuchten, ModelSettings, ThermalConductivityCapacity, TT) %{ This is to calculate vapor diffusivity, vapor dispersivity, and vapor enhancement factor. @@ -23,9 +23,6 @@ Eta = InitialValues.Eta; D_A = InitialValues.D_A; - % get model settings - ModelSettings = io.getModelSettings(); - for i = 1:ModelSettings.NL for j = 1:ModelSettings.nD MN = i + j - 1; diff --git a/src/+dryair/assembleCoefficientMatrices.m b/src/+dryair/assembleCoefficientMatrices.m index fea11a5d..70c6f2fc 100644 --- a/src/+dryair/assembleCoefficientMatrices.m +++ b/src/+dryair/assembleCoefficientMatrices.m @@ -1,4 +1,4 @@ -function [RHS, AirMatrices, SAVE] = assembleCoefficientMatrices(AirMatrices, SoilVariables, Delt_t, P_g, GroundwaterSettings) +function [RHS, AirMatrices, SAVE] = assembleCoefficientMatrices(AirMatrices, SoilVariables, Delt_t, P_g, ModelSettings, GroundwaterSettings) %{ Assemble the coefficient matrices of Equation 4.32 STEMMUS Technical Notes, page 44, for dry air equation. @@ -14,7 +14,6 @@ C6 = AirMatrices.C6; C7 = AirMatrices.C7; - ModelSettings = io.getModelSettings(); n = ModelSettings.NN; % Alias of SoilVariables diff --git a/src/+dryair/calculateBoundaryConditions.m b/src/+dryair/calculateBoundaryConditions.m index 917a49d8..387cc3a0 100644 --- a/src/+dryair/calculateBoundaryConditions.m +++ b/src/+dryair/calculateBoundaryConditions.m @@ -1,10 +1,9 @@ -function [RHS, AirMatrices] = calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT, P_gg, GroundwaterSettings) +function [RHS, AirMatrices] = calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT, P_gg, ModelSettings, GroundwaterSettings) %{ Determine the boundary condition for solving the dry air equation. %} TopPg = 100 .* (ForcingData.Pg_msr); - ModelSettings = io.getModelSettings(); n = ModelSettings.NN; if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa diff --git a/src/+dryair/calculateDryAirParameters.m b/src/+dryair/calculateDryAirParameters.m index e12a40b8..5204d99a 100644 --- a/src/+dryair/calculateDryAirParameters.m +++ b/src/+dryair/calculateDryAirParameters.m @@ -1,10 +1,9 @@ function AirVariabes = calculateDryAirParameters(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, VaporVariables, ... - P_gg, Xah, XaT, Xaa, RHODA, GroundwaterSettings) + P_gg, Xah, XaT, Xaa, RHODA, ModelSettings, GroundwaterSettings) %{ Calculate all the parameters related to dry air equation e.g., Equation 3.59-3.64, STEMMUS Technical Notes, page 27-28. %} - ModelSettings = io.getModelSettings(); Constants = io.define_constants(); AirVariabes.Cah = InitialValues.Cah; diff --git a/src/+dryair/calculateMatricCoefficients.m b/src/+dryair/calculateMatricCoefficients.m index 6f09d3ce..e1ec4802 100644 --- a/src/+dryair/calculateMatricCoefficients.m +++ b/src/+dryair/calculateMatricCoefficients.m @@ -1,12 +1,10 @@ -function AirMatrices = calculateMatricCoefficients(AirVariabes, InitialValues, GroundwaterSettings) +function AirMatrices = calculateMatricCoefficients(AirVariabes, InitialValues, ModelSettings, GroundwaterSettings) %{ Calculate all the parameters related to matric coefficients e.g., c1-c7 as in Equation 4.32 STEMMUS Technical Notes, page 44, which is an example for soil moisture equation, but for dry air equation. %} - ModelSettings = io.getModelSettings(); - AirMatrices.C1 = InitialValues.C1; AirMatrices.C2 = InitialValues.C2; AirMatrices.C3 = InitialValues.C3; diff --git a/src/+dryair/solveDryAirEquations.m b/src/+dryair/solveDryAirEquations.m index 7693017f..dd5c93cb 100644 --- a/src/+dryair/solveDryAirEquations.m +++ b/src/+dryair/solveDryAirEquations.m @@ -1,5 +1,5 @@ function [AirVariabes, RHS, SAVE, P_gg] = solveDryAirEquations(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, VaporVariables, ... - BoundaryCondition, ForcingData, P_gg, P_g, Xah, XaT, Xaa, RHODA, KT, Delt_t, GroundwaterSettings) + BoundaryCondition, ForcingData, P_gg, P_g, Xah, XaT, Xaa, RHODA, KT, Delt_t, ModelSettings, GroundwaterSettings) %{ Solve the dry air equation with the Thomas algorithm to update the soil air pressure 'P_gg', the finite difference time-stepping scheme is @@ -7,14 +7,14 @@ Technical Notes' section 4, Equation 4.32. %} AirVariabes = dryair.calculateDryAirParameters(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, VaporVariables, ... - P_gg, Xah, XaT, Xaa, RHODA, GroundwaterSettings); + P_gg, Xah, XaT, Xaa, RHODA, ModelSettings, GroundwaterSettings); - AirMatrices = dryair.calculateMatricCoefficients(AirVariabes, InitialValues, GroundwaterSettings); + AirMatrices = dryair.calculateMatricCoefficients(AirVariabes, InitialValues, ModelSettings, GroundwaterSettings); - [RHS, AirMatrices, SAVE] = dryair.assembleCoefficientMatrices(AirMatrices, SoilVariables, Delt_t, P_g, GroundwaterSettings); + [RHS, AirMatrices, SAVE] = dryair.assembleCoefficientMatrices(AirMatrices, SoilVariables, Delt_t, P_g, ModelSettings, GroundwaterSettings); - [RHS, AirMatrices] = dryair.calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT, P_gg, GroundwaterSettings); + [RHS, AirMatrices] = dryair.calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT, P_gg, ModelSettings, GroundwaterSettings); - [AirMatrices, P_gg, RHS] = dryair.solveTridiagonalMatrixEquations(RHS, AirMatrices, GroundwaterSettings); + [AirMatrices, P_gg, RHS] = dryair.solveTridiagonalMatrixEquations(RHS, AirMatrices, ModelSettings, GroundwaterSettings); end diff --git a/src/+dryair/solveTridiagonalMatrixEquations.m b/src/+dryair/solveTridiagonalMatrixEquations.m index b966a4e6..3c4e09e7 100644 --- a/src/+dryair/solveTridiagonalMatrixEquations.m +++ b/src/+dryair/solveTridiagonalMatrixEquations.m @@ -1,11 +1,9 @@ -function [AirMatrices, P_gg, RHS] = solveTridiagonalMatrixEquations(RHS, AirMatrices, GroundwaterSettings) +function [AirMatrices, P_gg, RHS] = solveTridiagonalMatrixEquations(RHS, AirMatrices, ModelSettings, GroundwaterSettings) %{ Use Thomas algorithm to solve the tridiagonal matrix equations, which is in the form of Equation 4.25 STEMMUS Technical Notes, page 41. %} - ModelSettings = io.getModelSettings(); - if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top else % Groundwater Coupling is activated diff --git a/src/+energy/assembleCoefficientMatrices.m b/src/+energy/assembleCoefficientMatrices.m index d95caaa0..d0e50ff8 100644 --- a/src/+energy/assembleCoefficientMatrices.m +++ b/src/+energy/assembleCoefficientMatrices.m @@ -1,4 +1,4 @@ -function [RHS, EnergyMatrices, SAVE] = assembleCoefficientMatrices(InitialValues, EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, GroundwaterSettings) +function [RHS, EnergyMatrices, SAVE] = assembleCoefficientMatrices(InitialValues, EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, ModelSettings, GroundwaterSettings) %{ assembles the coefficient matrices of Equation 4.32, STEMMUS Technical Notes, page 44, the example was only shown for the soil moisture @@ -14,7 +14,6 @@ C6 = EnergyMatrices.C6; C7 = EnergyMatrices.C7; - ModelSettings = io.getModelSettings(); n = ModelSettings.NN; if ModelSettings.Soilairefc == 1 % see https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/227 diff --git a/src/+energy/calculateBoundaryConditions.m b/src/+energy/calculateBoundaryConditions.m index 46e86213..4e59856c 100644 --- a/src/+energy/calculateBoundaryConditions.m +++ b/src/+energy/calculateBoundaryConditions.m @@ -1,13 +1,11 @@ -function [RHS, EnergyMatrices] = calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, ... - SoilVariables, Precip, EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, GroundwaterSettings) +function [RHS, EnergyMatrices] = calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, SoilVariables, ... + Precip, EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings) %{ Determine the boundary condition for solving the energy equation, see STEMMUS Technical Notes. %} - ModelSettings = io.getModelSettings(); n = ModelSettings.NN; - Constants = io.define_constants(); if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa diff --git a/src/+energy/calculateEnergyFluxes.m b/src/+energy/calculateEnergyFluxes.m index d8c6de3e..cbcae312 100644 --- a/src/+energy/calculateEnergyFluxes.m +++ b/src/+energy/calculateEnergyFluxes.m @@ -1,12 +1,9 @@ -function [QET, QEB] = calculateEnergyFluxes(SAVE, TT) +function [QET, QEB] = calculateEnergyFluxes(SAVE, TT, ModelSettings) %{ Calculate the energy fluxes on the boundary nodes, see STEMMUS Technical Notes. %} - ModelSettings = io.getModelSettings(); - n = ModelSettings.NN; - if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top else % Groundwater Coupling is activated @@ -15,6 +12,6 @@ end % So far these are unused vars - QET = SAVE(2, 1, 2) - SAVE(2, 2, 2) * TT(n - 1) - SAVE(2, 3, 2) * TT(n); + QET = SAVE(2, 1, 2) - SAVE(2, 2, 2) * TT(ModelSettings.NN - 1) - SAVE(2, 3, 2) * TT(ModelSettings.NN); QEB = -SAVE(1, 1, 2) + SAVE(1, 2, 2) * TT(indxBotm) + SAVE(1, 3, 2) * TT(indxBotm + 1); end diff --git a/src/+energy/calculateEnergyParameters.m b/src/+energy/calculateEnergyParameters.m index e9c2cf51..9ae53494 100644 --- a/src/+energy/calculateEnergyParameters.m +++ b/src/+energy/calculateEnergyParameters.m @@ -1,12 +1,11 @@ function EnergyVariables = calculateEnergyParameters(InitialValues, SoilVariables, HeatVariables, TransportCoefficient, AirVariabes, ... - VaporVariables, GasDispersivity, ThermalConductivityCapacity, ... - DRHOVh, DRHOVT, KL_T, Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, GroundwaterSettings) + VaporVariables, GasDispersivity, ThermalConductivityCapacity, DRHOVh, DRHOVT, ... + KL_T, Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, ModelSettings, GroundwaterSettings) %{ Calculate all the parameters related to energy balance equation e.Constants.g., Equation 3.65-3.73, STEMMUS Technical Notes, page 29-32. %} - ModelSettings = io.getModelSettings(); Constants = io.define_constants(); if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa diff --git a/src/+energy/calculateMatricCoefficients.m b/src/+energy/calculateMatricCoefficients.m index ce21243d..80b54202 100644 --- a/src/+energy/calculateMatricCoefficients.m +++ b/src/+energy/calculateMatricCoefficients.m @@ -1,10 +1,9 @@ -function EnergyMatrices = calculateMatricCoefficients(EnergyVariables, InitialValues, GroundwaterSettings) +function EnergyMatrices = calculateMatricCoefficients(EnergyVariables, InitialValues, ModelSettings, GroundwaterSettings) %{ Calculate all the parameters related to matric coefficients e.g., c1-c7 as in Equation 4.32 STEMMUS Technical Notes, page 44, which is an example for soil moisture equation, but here it is for energy equation. %} - ModelSettings = io.getModelSettings(); EnergyMatrices.C1 = InitialValues.C1; EnergyMatrices.C2 = InitialValues.C2; diff --git a/src/+energy/solveEnergyBalanceEquations.m b/src/+energy/solveEnergyBalanceEquations.m index 54877fae..14ce36bb 100644 --- a/src/+energy/solveEnergyBalanceEquations.m +++ b/src/+energy/solveEnergyBalanceEquations.m @@ -2,7 +2,7 @@ AirVariabes, VaporVariables, GasDispersivity, ThermalConductivityCapacity, ... HBoundaryFlux, BoundaryCondition, ForcingData, DRHOVh, DRHOVT, KL_T, ... Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, Delt_t, P_g, P_gg, ... - TOLD, Precip, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, GroundwaterSettings) + TOLD, Precip, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings) %{ Solve the Energy balance equation with the Thomas algorithm to update the soil temperature 'SoilVariables.TT', the finite difference @@ -11,17 +11,17 @@ %} EnergyVariables = energy.calculateEnergyParameters(InitialValues, SoilVariables, HeatVariables, TransportCoefficient, AirVariabes, ... - VaporVariables, GasDispersivity, ThermalConductivityCapacity, ... - DRHOVh, DRHOVT, KL_T, Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, GroundwaterSettings); + VaporVariables, GasDispersivity, ThermalConductivityCapacity, DRHOVh, DRHOVT, ... + KL_T, Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, ModelSettings, GroundwaterSettings); - EnergyMatrices = energy.calculateMatricCoefficients(EnergyVariables, InitialValues, GroundwaterSettings); + EnergyMatrices = energy.calculateMatricCoefficients(EnergyVariables, InitialValues, ModelSettings, GroundwaterSettings); - [RHS, EnergyMatrices, SAVE] = energy.assembleCoefficientMatrices(InitialValues, EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, GroundwaterSettings); + [RHS, EnergyMatrices, SAVE] = energy.assembleCoefficientMatrices(InitialValues, EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, ModelSettings, GroundwaterSettings); - [RHS, EnergyMatrices] = energy.calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, ... - SoilVariables, Precip, EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, GroundwaterSettings); + [RHS, EnergyMatrices] = energy.calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, SoilVariables, ... + Precip, EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings); - [SoilVariables, CHK, RHS, EnergyMatrices] = energy.solveTridiagonalMatrixEquations(EnergyMatrices, SoilVariables, RHS, CHK, GroundwaterSettings); + [SoilVariables, CHK, RHS, EnergyMatrices] = energy.solveTridiagonalMatrixEquations(EnergyMatrices, SoilVariables, RHS, CHK, ModelSettings, GroundwaterSettings); if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top @@ -30,7 +30,6 @@ indxBotm = GroundwaterSettings.indxBotmLayer; end - ModelSettings = io.getModelSettings(); if any(isnan(SoilVariables.TT)) || any(SoilVariables.TT(1:ModelSettings.NN) < ForcingData.Tmin) for i = indxBotm:ModelSettings.NN SoilVariables.TT(i) = TOLD(i); diff --git a/src/+energy/solveTridiagonalMatrixEquations.m b/src/+energy/solveTridiagonalMatrixEquations.m index e3ace4a4..66e459a2 100644 --- a/src/+energy/solveTridiagonalMatrixEquations.m +++ b/src/+energy/solveTridiagonalMatrixEquations.m @@ -1,11 +1,9 @@ -function [SoilVariables, CHK, RHS, EnergyMatrices] = solveTridiagonalMatrixEquations(EnergyMatrices, SoilVariables, RHS, CHK, GroundwaterSettings) +function [SoilVariables, CHK, RHS, EnergyMatrices] = solveTridiagonalMatrixEquations(EnergyMatrices, SoilVariables, RHS, CHK, ModelSettings, GroundwaterSettings) %{ Use Thomas algorithm to solve the tridiagonal matrix equations, which is in the form of Equation 4.25, STEMMUS Technical Notes, page 41. %} - ModelSettings = io.getModelSettings(); - if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top else % Groundwater Coupling is activated diff --git a/src/+groundwater/calculateGroundwaterRecharge.m b/src/+groundwater/calculateGroundwaterRecharge.m index 48c090dc..2435756e 100644 --- a/src/+groundwater/calculateGroundwaterRecharge.m +++ b/src/+groundwater/calculateGroundwaterRecharge.m @@ -1,4 +1,4 @@ -function [depToGWT_end, indxGWLay_end, gwfluxes] = calculateGroundwaterRecharge(EnergyVariables, SoilVariables, depToGWT_strt, indxGWLay_strt, KT, GroundwaterSettings) +function gwfluxes = calculateGroundwaterRecharge(EnergyVariables, SoilVariables, KT, ModelSettings, GroundwaterSettings) %{ Added by Mostafa, modified after Lianyu The concept followed to calculate groundwater recharge can be found in: @@ -11,20 +11,20 @@ The index order for the soil layers in STEMMUS is from bottom to top (NN:1), where NN is the number of STEMMUS soil layers, which is the opposite of MODFLOW (top to bottom). So, when converting information between STEMMUS and MODFLOW, indices need to be flipped + The concept of the moving balancing domain, used in the HYDRUS-MODFLOW paper and also STEMMUS-MODFLOW, is not adopted here (its effect is minor) + Instead, a different simple approach is followed (extract recharge from the layer above groundwater table) + %%%%%%%%%% Variables definitions %%%%%%%%%% Equations of the water balance are implemented in this function (Equations 8-13 of HYDRUS-MODFLOW paper) as follows: recharge = recharge_init + (SY - sy) * DeltZ, where: gwfluxes structure that includes the outputs (groundwater recharge and it's individual components) + indxRchrg index of the soil layer where the recharge is calculated recharge groundwater recharge, which is the upper boundary flux of the top layer of the phreatic aquifer (after the correction of the specific yield) recharge_init upper boundary flux into the moving balancing domain (before the correction of the specific yield) SS specific storage of MODFLOW aquifers SY large-scale specific yield of the phreatic aquifer sy small-scale specific yield (dynamically changing water yield) caused by fluctuation of the water table - depToGWT_strt water table depth: depth from top of the soil up to the phreatic surface at the start of the current time step - indxGWLay_strt index of the soil layer that includes the phreatic surface at the start of the current time step - depToGWT_end water table depth: depth from top of the soil up to the phreatic surface at the end of the current time step - indxGWLay_end index of the soil layer that includes the phreatic surface at the end of the current time step Theta_L soil moisture at the start of the current time step (bottom to top) Theta_LL soil moisture at the end of the current time step (bottom to top) STheta_L soil moisture at the start of the current time step (top to bottom) @@ -34,19 +34,7 @@ soilThick cumulative soil layers thickness (from top to bottom) aqlevels elevation of top surface level and all bottom levels of aquifer layers, received from MODFLOW through BMI %} - % Start Recharge calculations - % (a) Define the upper and lower boundaries of the moving balancing domain - % The moving balancing domain is located between depToGWT_strt and depToGWT_end - [depToGWT_end, indxGWLay_end] = groundwater.findPhreaticSurface(SoilVariables.hh, KT, GroundwaterSettings.soilThick, GroundwaterSettings.indxBotmLayer_R); - % indxRchrg and indxRchrgMax are the indices of the upper and lower levels of the moving boundary - % Following the HYDRUS-MODFLOW paper and also STEMMUS-MODFLOW, indxRchrg and indxRchrgMax are defined as in the next two lines - indxRchrgMax = max(indxGWLay_strt, indxGWLay_end) + 2; % the positive 2 is a user-specified value to define lower boundary of the moving boundary - % indxRchrg = min(indxGWLay_strt, indxGWLay_end) - 3; % the negative 2 or 3 is a user-specified value to define upper boundary of the moving boundary - % However, I comment the line above and use a slightly different way (ignore the moving layer boundaries and extract recharge from the two layers above groundwater table) - indxRchrg = min(indxGWLay_strt, indxGWLay_end) - 1; - indxRchrg_above = min(indxGWLay_strt, indxGWLay_end) - 2; - - % (b) Call the fluxes to get the initial recharge + % (a) Call the fluxes that contribute to the initial recharge % flip the fluxes (because STEMMUS calculations are from bottom to top, and MODFLOW needs the recharge from top to bottom) QLh_flip = flip(EnergyVariables.QL_h(1, :)); % liquid flux due to matric potential gradient QLT_flip = flip(EnergyVariables.QL_T(1, :)); % liquid flux due to temperature gradient @@ -56,21 +44,27 @@ soilThick cumulative soil layers thickness (from top to bottom) QVa_flip = flip(EnergyVariables.QVa(1, :)); % vapor water flux due to air pressure gradient Q_flip = flip(EnergyVariables.Qtot(1, :)); % total flux (liquid + vapor) - % (c) Get the recharge_init (before the correction of the specific yield)) + % (b) Get the recharge_init (before the correction of the specific yield)) % to avoid a zero flux value at the layer which recharge will be exported - if Q_flip(indxRchrg) == 0 - recharge_init = Q_flip(indxRchrg_above); - else - recharge_init = Q_flip(indxRchrg); % mean([Q_flip(indxRchrg), Q_flip(indxRchrg_above)]) + for j = 1:ModelSettings.NL - 1 + rech_lay = Q_flip(j); + rech_nextlay = Q_flip(j + 1); + if (rech_nextlay == 0 || rech_nextlay == -0) && rech_lay ~= 0 + recharge_init = rech_lay; + indxRchrg = j; % index of the recharge layer + break + else + continue + end end %{ - % (d) Calculations of SY + % (c) Calculations of SY % Note: In the HYDRUS-MODFLOW paper, Sy (from MODFLOW) was used. In Lianyu STEMMUS_MODFLOW code, a combination of Sy and Ss was used aqlevels = GroundwaterSettings.aqlevels; % elevation of top surface level and all bottom levels of aquifer layers numAqL = GroundwaterSettings.numAqL; % number of MODFLOW aquifer layers soilThick = GroundwaterSettings.soilThick; % cumulative soil layer thickness (from top to bottom) - indxAqLay = groundwater.calculateIndexAquifer(aqlevels, numAqL, soilThick); % index of MODFLOW aquifer layers for each STEMMUS soil layer + indxAqLay = groundwater.calculateIndexAquifer(aqlevels, numAqL, soilThick, ModelSettings.NN); % index of MODFLOW aquifer layers for each STEMMUS soil layer K = indxAqLay(indxGWLay_end); Thk = aqlevels(1) - aqlevels(K) - depToGWT_end; @@ -78,24 +72,21 @@ soilThick cumulative soil layers thickness (from top to bottom) SS = GroundwaterSettings.SS; S = (SY(K) - SS(K) * Thk) * (depToGWT_strt - depToGWT_end); - % (e) Calculations of sy - ModelSettings = io.getModelSettings(); - NN = ModelSettings.NN; % Number of nodes - NL = ModelSettings.NL; % Number of layers + % (d) Calculations of sy Theta_L = SoilVariables.Theta_L; Theta_LL = SoilVariables.Theta_LL; % Flip the soil moisture from top layer to bottom (opposite of Theta_L) - STheta_L(1) = Theta_L(NL, 2); - STheta_L(2:1:NN) = Theta_L(NN - 1:-1:1, 1); - STheta_LL(1) = Theta_LL(NL, 2); - STheta_LL(2:1:NN) = Theta_LL(NN - 1:-1:1, 1); + STheta_L(1) = Theta_L(ModelSettings.NL, 2); + STheta_L(2:1:ModelSettings.NN) = Theta_L(ModelSettings.NN - 1:-1:1, 1); + STheta_LL(1) = Theta_LL(ModelSettings.NL, 2); + STheta_LL(2:1:ModelSettings.NN) = Theta_LL(ModelSettings.NN - 1:-1:1, 1); sy = 0; for i = indxRchrg:indxRchrgMax - 1 sy = sy + 0.5 * (soilThick(i + 1) - soilThick(i)) * (STheta_LL(i) + STheta_LL(i + 1) - STheta_L(i) - STheta_L(i + 1)); end - % (f) Aggregate c, d, and e to get recharge + % (e) Aggregate b, c, and d to get recharge % recharge = recharge_init + S - sy; %} @@ -106,6 +97,12 @@ soilThick cumulative soil layers thickness (from top to bottom) gwfluxes.recharge = 0; end + % check 1 + diff = abs(GroundwaterSettings.indxBotmLayer_R - indxRchrg); + if diff > 1 + warning('Index of the bottom layer boundary that includes groundwater head ~= index of the recharge layer + 1!'); + end + % Outputs to be exported in csv gwfluxes.QLh = QLh_flip(indxRchrg); gwfluxes.QLT = QLT_flip(indxRchrg); diff --git a/src/+groundwater/calculateIndexAquifer.m b/src/+groundwater/calculateIndexAquifer.m index 6b6e63c8..829668d2 100644 --- a/src/+groundwater/calculateIndexAquifer.m +++ b/src/+groundwater/calculateIndexAquifer.m @@ -1,13 +1,10 @@ -function indxAqLay = calculateIndexAquifer(aqlevels, numAqL, soilThick) +function indxAqLay = calculateIndexAquifer(aqlevels, numAqL, soilThick, ModelSettings) %{ Assign the index of the MODFLOW aquifer that corresponds to each STEMMUS soil layer. aqlevels elevation of top surface level and all bottom levels of aquifer layers numAqL number of MODFLOW aquifer layers, received from MODFLOW through BMI %} - % Load model settings - ModelSettings = io.getModelSettings(); - numAqN = numAqL + 1; % number of MODFLOW aquifer nodes indxAqLay = zeros(ModelSettings.NN, 1); indxAqLay(1) = 1; diff --git a/src/+groundwater/calculateIndexBottomLayer.m b/src/+groundwater/calculateIndexBottomLayer.m index 6e215404..eedd762a 100644 --- a/src/+groundwater/calculateIndexBottomLayer.m +++ b/src/+groundwater/calculateIndexBottomLayer.m @@ -1,4 +1,4 @@ -function [indxBotmLayer, indxBotmLayer_R] = calculateIndexBottomLayer(soilThick, gw_Dep) +function [indxBotmLayer, indxBotmLayer_R] = calculateIndexBottomLayer(soilThick, gw_Dep, ModelSettings) %{ Calculate the index of the bottom layer level using MODFLOW data @@ -9,9 +9,6 @@ indxBotmLayer index of the bottom layer that contains the current headBotm indxBotmLayer_R = []; - % Load model settings - ModelSettings = io.getModelSettings(); - for i = 1:ModelSettings.NL midThick = (soilThick(i) + soilThick(i + 1)) / 2; if gw_Dep >= soilThick(i) && gw_Dep < soilThick(i + 1) diff --git a/src/+groundwater/findPhreaticSurface.m b/src/+groundwater/findPhreaticSurface.m deleted file mode 100644 index d5077d3e..00000000 --- a/src/+groundwater/findPhreaticSurface.m +++ /dev/null @@ -1,69 +0,0 @@ -function [depToGWT, indxGWLay] = findPhreaticSurface(soilWaterPotential, KT, soilThick, indxBotmLayer_R) - %{ - added by Mostafa, modified after Lianyu - This function finds the soil layer that includes the phreatic surface (saturated zone) and is used in the groundwater recharge calculations - - %%%%%%%%%% Important notes %%%%%%%%%% - The index order for the soil layers in STEMMUS is from bottom to top (NN:1), where NN is the number of STEMMUS soil layers, - which is the opposite of MODFLOW (top to bottom). So, when converting information between STEMMUS and MODFLOW, indices need to be flipped - - Although the groundwater depth (variable -> "gw_Dep") is calculated already from the MODFLOW inputs, it is re-calculated again in this function for two reasons: - 1) In this function, the groundwater depth is re-calculated, but this time is from the STEMMUS matric potential, so to ensure that ... - STEMMUS understood the location of the groundwater depth correctly. The re-calculated groundwater depth is stored in the "depToGWT" variable. - 2) The depToGWT needs to be assigned to a certain soil layer. Because of that assignment, there will be a slight difference between gw_Dep ... - and depToGWT (e.g. difference = 1-10 cm based on the thickness of the soil layer that contains the depToGWT). - - %%%%%%%%%% Variables definitions %%%%%%%%%% - depToGWT water table depth: depth from top soil layer to groundwater level, calculated from STEMMUS variables - indxBotmLayer_R index of the bottom layer that contains the current headBotmLayer (top to bottom) - indxGWLay index of the soil layer that includes the phreatic surface, so the recharge is .... - extracted from that layer (index order is from top layer to bottom) - Note: indxGWLay must equal to indxBotmLayer_R - Shh soil matric potential (from top layer to bottom; opposite of hh) - soilThick cumulative soil layers thickness (from top to bottom) - %} - - % Load model settings - ModelSettings = io.getModelSettings(); - NN = ModelSettings.NN; % number of nodes - - % Call the matric potential - Shh(1:1:NN) = soilWaterPotential(NN:-1:1); - depToGWT = 0; % starting value for initialization, updated below - indxGWLay = NN - 2; % starting value for initialization, updated below - - % Find the phreatic surface (saturated zone) - for i = NN:-1:2 - hh_lay = Shh(i); - soilThick_lay = soilThick(i); - hh_nextlay = Shh(i - 1); - soilThick_nextlay = soilThick(i - 1); - % apply a condition to find the groundwater table from the matric potential by differentiating .... - % between a) first layer with positive or zero head value and b) last layer with negative head value - if hh_lay > -1e-5 && hh_nextlay <= -1e-5 - depToGWT = (hh_lay * soilThick_nextlay - hh_nextlay * soilThick_lay) / (hh_lay - hh_nextlay); - midThick = (soilThick(i) + soilThick(i - 1)) / 2; - if depToGWT >= midThick - indxGWLay = i; - elseif depToGWT < midThick - indxGWLay = i - 1; - end - break - end - end - - if KT > 1 % start the checks from the first time step (after initialization) - % check 1 - if depToGWT <= 0 - warning('The phreatic surface level is equal or higher than the land surface level!'); - % check 2 - elseif depToGWT > ModelSettings.Tot_Depth % total soil depth - warning('The phreatic surface level is lower than the end of the soil column!'); - end - % check 3 - diff = abs(indxGWLay - indxBotmLayer_R); - if diff > 1 - warning('Index of the bottom layer boundary that includes groundwater head ~= index of the layer that has zero matric potential!'); - end - end -end diff --git a/src/+groundwater/initializeGroundwaterSettings.m b/src/+groundwater/initializeGroundwaterSettings.m index 999a4f3c..052a64e9 100644 --- a/src/+groundwater/initializeGroundwaterSettings.m +++ b/src/+groundwater/initializeGroundwaterSettings.m @@ -24,7 +24,7 @@ tempBotm groundwater temperature (C) at the bottom layer, received fr GroundwaterSettings.GroundwaterCoupling = 0; % (value = 0 -> deactivate coupling, or = 1 -> activate coupling); default = 0, update value to = 1 -> through BMI % Initialize the variables (head, temperature) at the bottom boundary (start of saturated zone) - GroundwaterSettings.headBotmLayer = 1950.0; % groundwater head (cm) at bottom layer + GroundwaterSettings.headBotmLayer = 1750.0; % groundwater head (cm) at bottom layer GroundwaterSettings.tempBotm = NaN; % groundwater temperature at bottom layer (C) % Call MODFLOW layers information (number of aquifer layers and their elevations, etc) diff --git a/src/+init/applySoilHeteroEffect.m b/src/+init/applySoilHeteroEffect.m index 1a85a68b..b956efe9 100644 --- a/src/+init/applySoilHeteroEffect.m +++ b/src/+init/applySoilHeteroEffect.m @@ -1,4 +1,4 @@ -function [SoilVariables, VanGenuchten] = applySoilHeteroEffect(SoilProperties, SoilVariables, VanGenuchten) +function [SoilVariables, VanGenuchten] = applySoilHeteroEffect(SoilProperties, SoilVariables, VanGenuchten, ModelSettings) initX = SoilVariables.InitialValues.initX; initND = SoilVariables.InitialValues.initND; @@ -8,7 +8,6 @@ ImpedF = repelem(3, 6); % get model settings - ModelSettings = io.getModelSettings(); Eqlspace = ModelSettings.Eqlspace; % Get soil constants for StartInit @@ -34,14 +33,14 @@ InitLnth(i) = ModelSettings.Tot_Depth - SoilConstants.Elmn_Lnth; for subRoutine = 5:-1:1 if abs(InitLnth(i) - initND(subRoutine)) < 1e-10 - [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i); + [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ModelSettings, ImpedF, Dmark, i); SoilVariables.InitialValues.initH = initH; Dmark = i + 2; end end if abs(InitLnth(i)) < 1e-10 subRoutine = 0; - [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i); + [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ModelSettings, ImpedF, Dmark, i); SoilVariables.InitialValues.initH = initH; end end diff --git a/src/+init/applySoilHeteroWithInitialFreezing.m b/src/+init/applySoilHeteroWithInitialFreezing.m index 07e97ae5..88ecfd7e 100644 --- a/src/+init/applySoilHeteroWithInitialFreezing.m +++ b/src/+init/applySoilHeteroWithInitialFreezing.m @@ -1,12 +1,9 @@ -function [SoilVariables] = applySoilHeteroWithInitialFreezing(LatentHeatOfFreezing, SoilVariables) +function [SoilVariables] = applySoilHeteroWithInitialFreezing(LatentHeatOfFreezing, SoilVariables, ModelSettings) SoilVariables.ISFT = 0; - % get model settings - ModelSettings = io.getModelSettings(); - for i = 1:ModelSettings.NN - SoilVariables.h_frez = init.updateHfreez(i, LatentHeatOfFreezing, SoilVariables); + SoilVariables.h_frez = init.updateHfreez(i, LatentHeatOfFreezing, SoilVariables, ModelSettings); SoilVariables.hh_frez(i) = SoilVariables.h_frez(i); SoilVariables.h(i) = SoilVariables.h(i) - SoilVariables.h_frez(i); SoilVariables.hh(i) = SoilVariables.h(i); diff --git a/src/+init/calculateInitialThermal.m b/src/+init/calculateInitialThermal.m index 3406d893..29532c66 100644 --- a/src/+init/calculateInitialThermal.m +++ b/src/+init/calculateInitialThermal.m @@ -1,4 +1,4 @@ -function ThermalConductivity = calculateInitialThermal(SoilVariables, VanGenuchten) +function ThermalConductivity = calculateInitialThermal(SoilVariables, VanGenuchten, ModelSettings) FEHCAP = []; % see issue 139 HCAP(1) = 0.998 * 4.182; @@ -25,9 +25,6 @@ GA1 = 0.035; GA2 = 0.013; - % get model settings - ModelSettings = io.getModelSettings(); - % Sum over all phases of dry porous media to find the dry heat capacity for j = 1:ModelSettings.NL % and the sums in the dry thermal conductivity; diff --git a/src/+init/defineInitialValues.m b/src/+init/defineInitialValues.m index af769214..77ca0bf6 100644 --- a/src/+init/defineInitialValues.m +++ b/src/+init/defineInitialValues.m @@ -1,9 +1,8 @@ -function InitialValues = defineInitialValues(Dur_tot) +function InitialValues = defineInitialValues(Dur_tot, ModelSettings) %{ %} % get model settings - ModelSettings = io.getModelSettings(); NN = ModelSettings.NN; % Number of nodes; mN = ModelSettings.mN; mL = ModelSettings.mL; % Number of elements. Prevending the exceeds of size of arraies; diff --git a/src/+init/setBoundaryCondition.m b/src/+init/setBoundaryCondition.m index 5e83869d..f65b0c84 100644 --- a/src/+init/setBoundaryCondition.m +++ b/src/+init/setBoundaryCondition.m @@ -1,4 +1,4 @@ -function BoundaryCondition = setBoundaryCondition(SoilVariables, ForcingData, initialLandcoverClass) +function BoundaryCondition = setBoundaryCondition(SoilVariables, ForcingData, initialLandcoverClass, ModelSettings) % Variables information for boundary condition settings % NBCh Indicator for type of surface boundary condition on mass euqation to be applied; % "1"--Specified matric head; @@ -50,9 +50,6 @@ Ta_msr = ForcingData.Ta_msr; - % get model settings - ModelSettings = io.getModelSettings(); - % NBCh: Moisture Surface B.C.: % 1 --Specified matric head(BCh); % 2 --Specified flux(BCh); diff --git a/src/+init/soilHeteroSubroutine.m b/src/+init/soilHeteroSubroutine.m index f9d45937..042619ba 100644 --- a/src/+init/soilHeteroSubroutine.m +++ b/src/+init/soilHeteroSubroutine.m @@ -1,4 +1,4 @@ -function [SoilVariables, VanGenuchten, initH] = soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, ML) +function [SoilVariables, VanGenuchten, initH] = soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ModelSettings, ImpedF, Dmark, ML) %{ considering soil heterogeneity effect %} @@ -7,9 +7,6 @@ initT = SoilVariables.InitialValues.initT; initH = SoilVariables.InitialValues.initH; - % get model settings - ModelSettings = io.getModelSettings(); - switch subRoutine case 0 from_id = Dmark; @@ -56,13 +53,13 @@ end J = SoilVariables.IS(i); - [SoilVariables, VanGenuchten] = init.updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, j, J); + [SoilVariables, VanGenuchten] = init.updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, ModelSettings, j, J); SoilVariables.Imped(i) = ImpedF(J); - initH(indexOfInit) = init.updateInitH(initX(indexOfInit), VanGenuchten, SoilVariables, j); + initH(indexOfInit) = init.updateInitH(initX(indexOfInit), VanGenuchten, SoilVariables, j, ModelSettings); if subRoutine == 5 - Btmh = init.updateBtmh(VanGenuchten, SoilVariables, i); + Btmh = init.updateBtmh(VanGenuchten, SoilVariables, i, ModelSettings); SoilVariables.T(i) = SoilVariables.BtmT + (i - 1) * (initT(indexOfInit) - SoilVariables.BtmT) / ML; SoilVariables.h(i) = Btmh + (i - 1) * (initH(indexOfInit) - Btmh) / ML; SoilVariables.IH(i) = 1; % Index of wetting history of soil which would be assumed as dry at the first with the value of 1 diff --git a/src/+init/updateBtmh.m b/src/+init/updateBtmh.m index cb0c1182..d5e4430a 100644 --- a/src/+init/updateBtmh.m +++ b/src/+init/updateBtmh.m @@ -1,6 +1,4 @@ -function Btmh = updateBtmh(VanGenuchten, SoilVariables, i) - % get model settings - ModelSettings = io.getModelSettings(); +function Btmh = updateBtmh(VanGenuchten, SoilVariables, i, ModelSettings) if ModelSettings.SWCC == 1 % VG soil water retention model Btmh = init.calcInitH(VanGenuchten.Theta_s(i), VanGenuchten.Theta_r(i), SoilVariables.BtmX, VanGenuchten.n(i), VanGenuchten.m(i), VanGenuchten.Alpha(i)); diff --git a/src/+init/updateHfreez.m b/src/+init/updateHfreez.m index 0d64ef48..19955471 100644 --- a/src/+init/updateHfreez.m +++ b/src/+init/updateHfreez.m @@ -1,4 +1,4 @@ -function h_frez = updateHfreez(i, LatentHeatOfFreezing, SoilVariables) +function h_frez = updateHfreez(i, LatentHeatOfFreezing, SoilVariables, ModelSettings) % get Constants Constants = io.define_constants(); @@ -7,9 +7,6 @@ h = SoilVariables.h; Phi_s = SoilVariables.Phi_s; - % get model settings - ModelSettings = io.getModelSettings(); - if T(i) <= 0 h_frez(i) = LatentHeatOfFreezing * 1e4 * (T(i)) / Constants.g / ModelSettings.T0; else diff --git a/src/+init/updateInitH.m b/src/+init/updateInitH.m index 6c73c33a..31017d2c 100644 --- a/src/+init/updateInitH.m +++ b/src/+init/updateInitH.m @@ -1,6 +1,4 @@ -function initH = updateInitH(initX, VanGenuchten, SoilVariables, i) - % get model settings - ModelSettings = io.getModelSettings(); +function initH = updateInitH(initX, VanGenuchten, SoilVariables, i, ModelSettings) if ModelSettings.SWCC == 1 % VG soil water retention model initH = init.calcInitH(VanGenuchten.Theta_s(i), VanGenuchten.Theta_r(i), initX, VanGenuchten.n(i), VanGenuchten.m(i), VanGenuchten.Alpha(i)); diff --git a/src/+init/updateSoilVariables.m b/src/+init/updateSoilVariables.m index 9273812e..1e2e786b 100644 --- a/src/+init/updateSoilVariables.m +++ b/src/+init/updateSoilVariables.m @@ -1,4 +1,4 @@ -function [SoilVariables, VanGenuchten] = updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, i, j) +function [SoilVariables, VanGenuchten] = updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, ModelSettings, i, j) % Get soil constants for StartInit SoilConstants = io.getSoilConstants(); @@ -12,9 +12,6 @@ SoilVariables.XSOC(i) = SoilVariables.VPERSOC(j); SoilVariables.XK(i) = SoilConstants.XK; - % get model settings - ModelSettings = io.getModelSettings(); - if ModelSettings.SWCC == 1 % VG soil water retention model VanGenuchten.Theta_s(i) = SoilProperties.SaturatedMC(j); VanGenuchten.Theta_r(i) = SoilProperties.ResidualMC(j); diff --git a/src/+groundwater/calculateSoilLayerThickness.m b/src/+io/calculateSoilLayerThickness.m similarity index 81% rename from src/+groundwater/calculateSoilLayerThickness.m rename to src/+io/calculateSoilLayerThickness.m index 59ddf85d..dae1f9b0 100644 --- a/src/+groundwater/calculateSoilLayerThickness.m +++ b/src/+io/calculateSoilLayerThickness.m @@ -1,4 +1,4 @@ -function soilThick = calculateSoilLayerThickness() +function soilThick = calculateSoilLayerThickness(ModelSettings) %{ Calculate soil layers thickness (cumulative layers thickness; e.g. 1, 2, 3, 5, 10, ......., 480, total soil depth) @@ -7,9 +7,6 @@ soilThick cumulative soil layers thickness (from top to bottom) %} - % Load model settings - ModelSettings = io.getModelSettings(); - soilThick = zeros(ModelSettings.NN, 1); % cumulative soil layers thickness soilThick(1) = 0; diff --git a/src/+io/getModelSettings.m b/src/+io/getModelSettings.m index b39cc6a6..978feddd 100644 --- a/src/+io/getModelSettings.m +++ b/src/+io/getModelSettings.m @@ -73,9 +73,4 @@ ModelSettings.ML = ML; ModelSettings.DeltZ = DeltZ; ModelSettings.DeltZ_R = DeltZ_R; - - ModelSettings.NN = ModelSettings.NL + 1; % Number of nodes; - ModelSettings.mN = ModelSettings.NN + 1; - ModelSettings.mL = ModelSettings.NL + 1; % Number of elements. Prevending the exceeds of size of arraies; - ModelSettings.nD = 2; end diff --git a/src/+io/loadSoilInitialValues.m b/src/+io/loadSoilInitialValues.m index ebf04fdb..94aaa06f 100644 --- a/src/+io/loadSoilInitialValues.m +++ b/src/+io/loadSoilInitialValues.m @@ -1,4 +1,4 @@ -function [InitialValues, BtmX, BtmT, Tss] = loadSoilInitialValues(InputPath, TimeProperties, SoilProperties, ForcingData) +function [InitialValues, BtmX, BtmT, Tss] = loadSoilInitialValues(InputPath, TimeProperties, SoilProperties, ForcingData, ModelSettings) %{ Producing initial soil moisture and soil temperature profile %} @@ -9,7 +9,6 @@ Ta_msr = ForcingData.Ta_msr; % Get model settings - ModelSettings = io.getModelSettings(); Tot_Depth = ModelSettings.Tot_Depth; SWCC = ModelSettings.SWCC; diff --git a/src/+io/readSoilLayerSettings.m b/src/+io/readSoilLayerSettings.m new file mode 100644 index 00000000..daacaf14 --- /dev/null +++ b/src/+io/readSoilLayerSettings.m @@ -0,0 +1,12 @@ +function SoilLayerSettings = readSoilLayerSettings(soilLayersFile) + soildata = dlmread(soilLayersFile, ',', 1, 0); % skip column names + + SoilLayerSettings.NL = soildata(end, 1); + SoilLayerSettings.ML = SoilLayerSettings.NL; + + SoilLayerSettings.DeltZ_R = transpose(soildata(:, 2)); + SoilLayerSettings.DeltZ = flip(SoilLayerSettings.DeltZ_R); + + SoilLayerSettings.Tot_Depth = sum(SoilLayerSettings.DeltZ_R); + SoilLayerSettings.R_depth = soildata(1, 3); +end diff --git a/src/+soilmoisture/assembleCoefficientMatrices.m b/src/+soilmoisture/assembleCoefficientMatrices.m index 11955776..99a5de08 100644 --- a/src/+soilmoisture/assembleCoefficientMatrices.m +++ b/src/+soilmoisture/assembleCoefficientMatrices.m @@ -1,9 +1,8 @@ -function HeatMatrices = assembleCoefficientMatrices(HeatVariables, InitialValues, Srt, GroundwaterSettings) +function HeatMatrices = assembleCoefficientMatrices(HeatVariables, InitialValues, Srt, ModelSettings, GroundwaterSettings) %{ Assemble the coefficient matrices of Equation 4.32 (STEMMUS Technical Notes, page 44). %} - ModelSettings = io.getModelSettings(); % Define HeatMatrices structure HeatMatrices.C1 = InitialValues.C1; diff --git a/src/+soilmoisture/calculateBoundaryConditions.m b/src/+soilmoisture/calculateBoundaryConditions.m index 6d1fff71..960c98d4 100644 --- a/src/+soilmoisture/calculateBoundaryConditions.m +++ b/src/+soilmoisture/calculateBoundaryConditions.m @@ -1,9 +1,9 @@ function [AVAIL0, RHS, HeatMatrices, Precip, ForcingData] = calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ... - TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap, GroundwaterSettings) + TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap, ModelSettings, GroundwaterSettings) %{ Determine the boundary condition for solving the soil moisture equation. %} - ModelSettings = io.getModelSettings(); + % n is the index of n_th item n = ModelSettings.NN; @@ -33,7 +33,6 @@ soilThick = GroundwaterSettings.soilThick; % cumulative soil layers thickness topLevel = GroundwaterSettings.topLevel; headBotmLayer = GroundwaterSettings.headBotmLayer; % groundwater head at bottom layer, received from MODFLOW through BMI - RHS(indxBotm) = headBotmLayer - topLevel + soilThick(indxBotmLayer_R); % (RHS = zero at the end, need to check with Yijian and Lianyu) if BoundaryCondition.NBChB == 1 % Specify matric head at bottom to be ---BoundaryCondition.BChB; RHS(indxBotm) = headBotmLayer - topLevel + soilThick(indxBotmLayer_R); C4(indxBotm, 1) = 1; @@ -46,6 +45,7 @@ RHS(indxBotm) = RHS(indxBotm) - SoilVariables.KL_h(indxBotm, 1); end end + % Apply the surface boundary condition called for by BoundaryCondition.NBCh if BoundaryCondition.NBCh == 1 % Specified matric head at surface---equal to hN; % h_SUR: Observed matric potential at surface. This variable diff --git a/src/+soilmoisture/calculateEvapotranspiration.m b/src/+soilmoisture/calculateEvapotranspiration.m index d240af1e..450a5d5d 100644 --- a/src/+soilmoisture/calculateEvapotranspiration.m +++ b/src/+soilmoisture/calculateEvapotranspiration.m @@ -1,6 +1,4 @@ -function [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, RWUs, RWUg] = calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt, GroundwaterSettings) - - ModelSettings = io.getModelSettings(); +function [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, RWUs, RWUg] = calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt, ModelSettings, GroundwaterSettings) if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top diff --git a/src/+soilmoisture/calculateMatricCoefficients.m b/src/+soilmoisture/calculateMatricCoefficients.m index 0aef14bf..39e67789 100644 --- a/src/+soilmoisture/calculateMatricCoefficients.m +++ b/src/+soilmoisture/calculateMatricCoefficients.m @@ -1,11 +1,10 @@ function [HeatVariables, SoilVariables] = calculateMatricCoefficients(SoilVariables, VaporVariables, GasDispersivity, InitialValues, ... - RHOV, DRHOVh, DRHOVT, D_Ta, GroundwaterSettings) + RHOV, DRHOVh, DRHOVT, D_Ta, ModelSettings, GroundwaterSettings) %{ Calculate all the parameters related to matric coefficients (e.g., c1-c8) as in Equation 4.32 (STEMMUS Technical Notes, page 44). %} - ModelSettings = io.getModelSettings(); Constants = io.define_constants(); % Add a new variables to SoilVariables diff --git a/src/+soilmoisture/calculateTimeDerivatives.m b/src/+soilmoisture/calculateTimeDerivatives.m index a2544abb..c0eff4c9 100644 --- a/src/+soilmoisture/calculateTimeDerivatives.m +++ b/src/+soilmoisture/calculateTimeDerivatives.m @@ -1,4 +1,4 @@ -function [RHS, HeatMatrices, boundaryFluxArray] = calculateTimeDerivatives(InitialValues, SoilVariables, HeatMatrices, Delt_t, P_gg, GroundwaterSettings) +function [RHS, HeatMatrices, boundaryFluxArray] = calculateTimeDerivatives(InitialValues, SoilVariables, HeatMatrices, Delt_t, P_gg, ModelSettings, GroundwaterSettings) %{ Perform the finite difference of the time derivatives of the matrix equation as Equation 4.32 shows in 'STEMMUS Technical Notes', section 4. @@ -20,7 +20,6 @@ boundaryFluxArray = InitialValues.SAVE; RHS = InitialValues.RHS; - ModelSettings = io.getModelSettings(); % n is the index of n_th item n = ModelSettings.NN; diff --git a/src/+soilmoisture/calculatesSoilWaterFluxes.m b/src/+soilmoisture/calculatesSoilWaterFluxes.m index f5f97b91..86a5abf3 100644 --- a/src/+soilmoisture/calculatesSoilWaterFluxes.m +++ b/src/+soilmoisture/calculatesSoilWaterFluxes.m @@ -1,8 +1,7 @@ -function HBoundaryFlux = calculatesSoilWaterFluxes(SAVE, hh, GroundwaterSettings) +function HBoundaryFlux = calculatesSoilWaterFluxes(SAVE, hh, ModelSettings, GroundwaterSettings) %{ Calculate the soil water fluxes on the boundary node. %} - ModelSettings = io.getModelSettings(); if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top diff --git a/src/+soilmoisture/solveSoilMoistureBalance.m b/src/+soilmoisture/solveSoilMoistureBalance.m index 0d061b86..8d9a6bde 100644 --- a/src/+soilmoisture/solveSoilMoistureBalance.m +++ b/src/+soilmoisture/solveSoilMoistureBalance.m @@ -1,5 +1,5 @@ function [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, ... - BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, GroundwaterSettings) + BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings) %{ Solve the soil moisture balance equation with the Thomas algorithm to update the soil matric potential 'hh', the finite difference @@ -7,16 +7,16 @@ 'STEMMUS Technical Notes' section 4, Equation 4.32. %} [HeatVariables, SoilVariables] = soilmoisture.calculateMatricCoefficients(SoilVariables, VaporVariables, GasDispersivity, InitialValues, ... - RHOV, DRHOVh, DRHOVT, D_Ta, GroundwaterSettings); + RHOV, DRHOVh, DRHOVT, D_Ta, ModelSettings, GroundwaterSettings); - HeatMatrices = soilmoisture.assembleCoefficientMatrices(HeatVariables, InitialValues, Srt, GroundwaterSettings); + HeatMatrices = soilmoisture.assembleCoefficientMatrices(HeatVariables, InitialValues, Srt, ModelSettings, GroundwaterSettings); - [RHS, HeatMatrices, boundaryFluxArray] = soilmoisture.calculateTimeDerivatives(InitialValues, SoilVariables, HeatMatrices, Delt_t, P_gg, GroundwaterSettings); + [RHS, HeatMatrices, boundaryFluxArray] = soilmoisture.calculateTimeDerivatives(InitialValues, SoilVariables, HeatMatrices, Delt_t, P_gg, ModelSettings, GroundwaterSettings); % When BoundaryCondition.NBCh == 3, otherwise Rn_SOIL, Evap, EVAP, Trap, % r_a_SOIL, Srt will be empty arrays! see issue 98, item 2 if BoundaryCondition.NBCh == 3 - [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, RWUs, RWUg] = soilmoisture.calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt, GroundwaterSettings); + [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, RWUs, RWUg] = soilmoisture.calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt, ModelSettings, GroundwaterSettings); else Rn_SOIL = InitialValues.Rn_SOIL; Evap = InitialValues.Evap; @@ -26,16 +26,15 @@ end [AVAIL0, RHS, HeatMatrices, Precip, ForcingData] = soilmoisture.calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ... - TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap, GroundwaterSettings); + TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap, ModelSettings, GroundwaterSettings); + + [CHK, hh, C4] = soilmoisture.solveTridiagonalMatrixEquations(HeatMatrices.C4, SoilVariables.hh, HeatMatrices.C4_a, RHS, ModelSettings, GroundwaterSettings); - [CHK, hh, C4] = soilmoisture.solveTridiagonalMatrixEquations(HeatMatrices.C4, SoilVariables.hh, HeatMatrices.C4_a, RHS, GroundwaterSettings); % update structures SoilVariables.hh = hh; HeatMatrices.C4 = C4; % fix hh values - ModelSettings = io.getModelSettings(); - if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top else % Groundwater Coupling is activated @@ -53,6 +52,6 @@ end % calculate boundary flux - HBoundaryFlux = soilmoisture.calculatesSoilWaterFluxes(boundaryFluxArray, SoilVariables.hh, GroundwaterSettings); + HBoundaryFlux = soilmoisture.calculatesSoilWaterFluxes(boundaryFluxArray, SoilVariables.hh, ModelSettings, GroundwaterSettings); end diff --git a/src/+soilmoisture/solveTridiagonalMatrixEquations.m b/src/+soilmoisture/solveTridiagonalMatrixEquations.m index 52cb54ea..c655566c 100644 --- a/src/+soilmoisture/solveTridiagonalMatrixEquations.m +++ b/src/+soilmoisture/solveTridiagonalMatrixEquations.m @@ -1,9 +1,8 @@ -function [CHK, hh, C4] = solveTridiagonalMatrixEquations(C4, hh, C4_a, RHS, GroundwaterSettings) +function [CHK, hh, C4] = solveTridiagonalMatrixEquations(C4, hh, C4_a, RHS, ModelSettings, GroundwaterSettings) %{ Solve the tridiagonal matrix equations using Thomas algorithm, which is in the form of Equation 4.25 (STEMMUS Technical Notes, page 41). %} - ModelSettings = io.getModelSettings(); if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling, added by Mostafa indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top diff --git a/src/Max_Rootdepth.m b/src/Max_Rootdepth.m index fcad0cc8..ee6ff9db 100644 --- a/src/Max_Rootdepth.m +++ b/src/Max_Rootdepth.m @@ -1,13 +1,10 @@ -function [bbx] = Max_Rootdepth(bbx) +function [bbx] = Max_Rootdepth(bbx, ModelSettings) %{ The function identify if there is root in certain layer, bbx=0 indicates no roots while bbx=1 indicates roots existing %} - % get model settings - ModelSettings = io.getModelSettings(); - % BR = 10:1:650; %% [gC /m^2 PFT] % rroot = 0.5*1e-3 ; % 3.3*1e-4 ;%% [0.5-6 *10^-3] [m] root radius diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index f53b232e..351cbd24 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -46,16 +46,33 @@ % Load model settings: replacing "run Constants" ModelSettings = io.getModelSettings(); - NN = ModelSettings.NN; + + % Load soil layers settings + soilLayersFile = [InputPath, 'input_soilLayThick.csv']; + if isfile(soilLayersFile) + SoilLayerSettings = io.readSoilLayerSettings(soilLayersFile); + % Override the default settings + fields = fieldnames(SoilLayerSettings); + for i = 1:numel(fields) + ModelSettings.(fields{i}) = SoilLayerSettings.(fields{i}); + end + end + + % Calculate other ModelSettings + ModelSettings.NN = ModelSettings.NL + 1; % Number of nodes; + ModelSettings.mN = ModelSettings.NN + 1; + ModelSettings.mL = ModelSettings.NL + 1; % Number of elements. Preventing the exceeds of arrays size; + ModelSettings.nD = 2; % Load groundwater settings GroundwaterSettings = groundwater.initializeGroundwaterSettings(); + GroundwaterSettings.soilThick = io.calculateSoilLayerThickness(ModelSettings); % load forcing data ForcingData = io.loadForcingData(InputPath, TimeProperties, SoilProperties, ModelSettings.Tot_Depth, GroundwaterSettings); % Get initial values - InitialValues = init.defineInitialValues(TimeProperties.Dur_tot); + InitialValues = init.defineInitialValues(TimeProperties.Dur_tot, ModelSettings); T_CRIT = InitialValues.T_CRIT; TT_CRIT = InitialValues.TT_CRIT; hOLD = InitialValues.hOLD; @@ -234,14 +251,14 @@ SoilVariables = init.defineSoilVariables(InitialValues, SoilProperties, VanGenuchten); % Add initial soil moisture and soil temperature - [SoilInitialValues, BtmX, BtmT, Tss] = io.loadSoilInitialValues(InputPath, TimeProperties, SoilProperties, ForcingData); + [SoilInitialValues, BtmX, BtmT, Tss] = io.loadSoilInitialValues(InputPath, TimeProperties, SoilProperties, ForcingData, ModelSettings); SoilVariables.InitialValues = SoilInitialValues; SoilVariables.BtmX = BtmX; SoilVariables.BtmT = BtmT; SoilVariables.Tss = Tss; % Run StartInit - [SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten); + [SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten, ModelSettings); % Set SoilVariables that are used in the loop T = SoilVariables.T; @@ -255,7 +272,7 @@ SoilConstants = io.getSoilConstants(); %% The boundary condition information settings - BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, SiteProperties.landcoverClass(1)); + BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, SiteProperties.landcoverClass(1), ModelSettings); DSTOR = BoundaryCondition.DSTOR; DSTOR0 = BoundaryCondition.DSTOR0; RS = BoundaryCondition.RS; @@ -266,10 +283,10 @@ %% 14. Initialize matrices calculate = 1; TIMEOLD = 0; - SAVEhh_frez = zeros(NN + 1, 1); - FCHK = zeros(1, NN); - KCHK = zeros(1, NN); - hCHK = zeros(1, NN); + SAVEhh_frez = zeros(ModelSettings.NN + 1, 1); + FCHK = zeros(1, ModelSettings.NN); + KCHK = zeros(1, ModelSettings.NN); + hCHK = zeros(1, ModelSettings.NN); TIMELAST = 0; % the start of simulation period is from 0mins, while the input data start from 30mins. @@ -282,9 +299,6 @@ Delt_t0 = Delt_t; % Duration of last time step TOLD_CRIT = []; - % 15. Calculate soil layer thickness - GroundwaterSettings.soilThick = groundwater.calculateSoilLayerThickness(); - % for soil moisture and temperature outputs monitorDepthTemperature = ModelSettings.NL:-1:1; monitorDepthSoilMoisture = ModelSettings.NL:-1:1; @@ -343,8 +357,7 @@ [ForcingData.R_Dunn, ForcingData.Precip_msr] = groundwater.updateDunnianRunoff(ForcingData.Precip_msr, GroundwaterSettings.gw_Dep); % Calculate the index of the bottom layer level - [GroundwaterSettings.indxBotmLayer, GroundwaterSettings.indxBotmLayer_R] = groundwater.calculateIndexBottomLayer(GroundwaterSettings.soilThick, GroundwaterSettings.gw_Dep); - [depToGWT_strt, indxGWLay_strt] = groundwater.findPhreaticSurface(SoilVariables.hh, KT, GroundwaterSettings.soilThick, GroundwaterSettings.indxBotmLayer_R); + [GroundwaterSettings.indxBotmLayer, GroundwaterSettings.indxBotmLayer_R] = groundwater.calculateIndexBottomLayer(GroundwaterSettings.soilThick, GroundwaterSettings.gw_Dep, ModelSettings); % Check the position of the groundwater table if any(GroundwaterSettings.gw_Dep > ModelSettings.Tot_Depth) || any(isempty(GroundwaterSettings.indxBotmLayer)) @@ -376,10 +389,10 @@ end %%%%% Updating the state variables. %%%%%%%%%%%%%%%%%%%%%%%%%%%% L_f = 0; % ignore Freeze/Thaw, see issue 139 - TT_CRIT(NN) = ModelSettings.T0; % unit K + TT_CRIT(ModelSettings.NN) = ModelSettings.T0; % unit K hOLD_frez = []; if IRPT1 == 0 && IRPT2 == 0 && SoilVariables.ISFT == 0 - for MN = 1:NN + for MN = 1:ModelSettings.NN hOLD_frez(MN) = h_frez(MN); h_frez(MN) = hh_frez(MN); TOLD_CRIT(MN) = T_CRIT(MN); @@ -485,7 +498,7 @@ [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential] ... = ebal(iter, options, spectral, rad, gap, ... leafopt, angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ... - Rl, SoilVariables, VanGenuchten, InitialValues, GroundwaterSettings); + Rl, SoilVariables, VanGenuchten, InitialValues, ModelSettings, GroundwaterSettings); if options.calc_fluor if options.calc_vert_profiles [rad, profiles] = RTMf(spectral, rad, soil, leafopt, canopy, gap, angles, profiles); @@ -558,7 +571,7 @@ BoundaryCondition.DSTOR0 = DSTOR0; if KT > 1 - SoilVariables.XWRE = updateWettingHistory(SoilVariables, VanGenuchten); + SoilVariables.XWRE = updateWettingHistory(SoilVariables, VanGenuchten, ModelSettings); end for i = 1:ModelSettings.NL @@ -566,37 +579,37 @@ end end if Delt_t ~= Delt_t0 - for MN = 1:NN + for MN = 1:ModelSettings.NN hh(MN) = h(MN) + (h(MN) - hOLD(MN)) * Delt_t / Delt_t0; TT(MN) = T(MN) + (T(MN) - TOLD(MN)) * Delt_t / Delt_t0; end end - hSAVE = hh(NN); - TSAVE = TT(NN); + hSAVE = hh(ModelSettings.NN); + TSAVE = TT(ModelSettings.NN); % set "hN" empty when the "if statement" below does not happen, see issue 98, % item 5 hN = []; if BoundaryCondition.NBCh == 1 hN = BoundaryCondition.BCh; - hh(NN) = hN; + hh(ModelSettings.NN) = hN; hSAVE = hN; elseif BoundaryCondition.NBCh == 2 if BoundaryCondition.NBChh ~= 2 if BoundaryCondition.BCh < 0 hN = DSTOR0; - hh(NN) = hN; + hh(ModelSettings.NN) = hN; hSAVE = hN; else hN = -1e6; - hh(NN) = hN; + hh(ModelSettings.NN) = hN; hSAVE = hN; end end else if BoundaryCondition.NBChh ~= 2 hN = DSTOR0; - hh(NN) = hN; + hh(ModelSettings.NN) = hN; hSAVE = hN; end end @@ -610,40 +623,40 @@ SoilVariables.Tss(KT) = Tss; for KIT = 1:ModelSettings.NIT - [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, NN, SoilConstants.hd, ForcingData.Tmin); + [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, ModelSettings.NN, SoilConstants.hd, ForcingData.Tmin); % update inputs for UpdateSoilWaterContent SoilVariables.TT_CRIT = TT_CRIT; SoilVariables.hh_frez = hh_frez; - SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten); + SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten, ModelSettings); % Reset KL_T here. CondL_T script is replaced by this line % see issue 181, item 4 KL_T = InitialValues.KL_T; - [RHOV, DRHOVh, DRHOVT] = Density_V(SoilVariables.TT, SoilVariables.hh, Constants.g, Constants.Rv, NN); + [RHOV, DRHOVh, DRHOVT] = Density_V(SoilVariables.TT, SoilVariables.hh, Constants.g, Constants.Rv, ModelSettings.NN); - TransportCoefficient = conductivity.calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t); + TransportCoefficient = conductivity.calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t, ModelSettings); W = TransportCoefficient.W; D_Ta = TransportCoefficient.D_Ta; - [L] = Latent(SoilVariables.TT, NN); + [L] = Latent(SoilVariables.TT, ModelSettings.NN); % DRHODAt unused! - [Xaa, XaT, Xah, DRHODAt, DRHODAz, RHODA] = Density_DA(SoilVariables.T, Constants.RDA, P_g, Constants.Rv, ModelSettings.DeltZ, SoilVariables.h, SoilVariables.hh, SoilVariables.TT, P_gg, Delt_t, ModelSettings.NL, NN, DRHOVT, DRHOVh, RHOV); + [Xaa, XaT, Xah, DRHODAt, DRHODAz, RHODA] = Density_DA(SoilVariables.T, Constants.RDA, P_g, Constants.Rv, ModelSettings.DeltZ, SoilVariables.h, SoilVariables.hh, SoilVariables.TT, P_gg, Delt_t, ModelSettings.NL, ModelSettings.NN, DRHOVT, DRHOVh, RHOV); - ThermalConductivityCapacity = conductivity.calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV); + ThermalConductivityCapacity = conductivity.calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, ModelSettings, DRHOVT, L, RHOV); - k_g = conductivity.calculateGasConductivity(InitialValues, TransportCoefficient, VanGenuchten, SoilVariables); + k_g = conductivity.calculateGasConductivity(InitialValues, TransportCoefficient, VanGenuchten, SoilVariables, ModelSettings); - VaporVariables = conductivity.calculateVaporVariables(InitialValues, SoilVariables, VanGenuchten, ThermalConductivityCapacity, SoilVariables.TT); + VaporVariables = conductivity.calculateVaporVariables(InitialValues, SoilVariables, VanGenuchten, ModelSettings, ThermalConductivityCapacity, SoilVariables.TT); - GasDispersivity = conductivity.calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g); + GasDispersivity = conductivity.calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g, ModelSettings); % Srt is both input and output [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = ... soilmoisture.solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, ... TimeProperties, SoilProperties, BoundaryCondition, Delt_t, RHOV, DRHOVh, ... - DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, GroundwaterSettings); + DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings); if BoundaryCondition.NBCh == 1 DSTOR = 0; @@ -670,7 +683,7 @@ if ModelSettings.Soilairefc == 1 [AirVariabes, RHS, SAVE, P_gg] = dryair.solveDryAirEquations(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, VaporVariables, ... - BoundaryCondition, ForcingData, P_gg, P_g, Xah, XaT, Xaa, RHODA, KT, Delt_t, GroundwaterSettings); + BoundaryCondition, ForcingData, P_gg, P_g, Xah, XaT, Xaa, RHODA, KT, Delt_t, ModelSettings, GroundwaterSettings); else AirVariabes.KLhBAR = InitialValues.KLhBAR; AirVariabes.KLTBAR = InitialValues.KLTBAR; @@ -688,25 +701,25 @@ AirVariabes, VaporVariables, GasDispersivity, ThermalConductivityCapacity, ... HBoundaryFlux, BoundaryCondition, ForcingData, DRHOVh, DRHOVT, KL_T, ... Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, Delt_t, P_g, P_gg, ... - TOLD, Precip, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, GroundwaterSettings); + TOLD, Precip, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings); end if max(CHK) < 0.1 break end - hSAVE = SoilVariables.hh(NN); - TSAVE = SoilVariables.TT(NN); + hSAVE = SoilVariables.hh(ModelSettings.NN); + TSAVE = SoilVariables.TT(ModelSettings.NN); end TIMEOLD = KT; KIT; KIT = 0; - [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, NN, SoilConstants.hd, ForcingData.Tmin); + [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, ModelSettings.NN, SoilConstants.hd, ForcingData.Tmin); % updates inputs for UpdateSoilWaterContent SoilVariables.TT_CRIT = TT_CRIT; SoilVariables.hh_frez = hh_frez; - SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten); + SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten, ModelSettings); if IRPT1 == 0 && IRPT2 == 0 if KT % In case last time step is not convergent and needs to be repeated. @@ -742,7 +755,7 @@ Sim_qtot(KT, 1:length(monitorDepthSoilMoisture)) = qtot(monitorDepthSoilMoisture, KT); end if (TEND - TIME) < 1E-3 - for i = 1:NN + for i = 1:ModelSettings.NN hOLD(i) = SoilVariables.h(i); SoilVariables.h(i) = SoilVariables.hh(i); if ModelSettings.Thmrlefc == 1 @@ -764,8 +777,7 @@ % Recharge calculations, added by Mostafa if GroundwaterSettings.GroundwaterCoupling == 1 % Groundwater coupling is enabled - % update depToGWT_strt, indxGWLay_strt for next time step - [depToGWT_strt, indxGWLay_strt, gwfluxes] = groundwater.calculateGroundwaterRecharge(EnergyVariables, SoilVariables, depToGWT_strt, indxGWLay_strt, KT, GroundwaterSettings); + gwfluxes = groundwater.calculateGroundwaterRecharge(EnergyVariables, SoilVariables, KT, ModelSettings, GroundwaterSettings); if GroundwaterSettings.gw_Dep <= 1 % soil is fully saturated gwfluxes.recharge = 0; end diff --git a/src/StartInit.m b/src/StartInit.m index b1a6eaad..ee3fa77f 100644 --- a/src/StartInit.m +++ b/src/StartInit.m @@ -1,4 +1,4 @@ -function [SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten) +function [SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten, ModelSettings) Ksh = repelem(18 / (3600 * 24), 6); BtmKsh = Ksh(6); @@ -10,22 +10,22 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Considering soil hetero effect modify date: 20170103 %%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - [SoilVariables, VanGenuchten] = init.applySoilHeteroEffect(SoilProperties, SoilVariables, VanGenuchten); + [SoilVariables, VanGenuchten] = init.applySoilHeteroEffect(SoilProperties, SoilVariables, VanGenuchten, ModelSettings); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Considering soil hetero effect modify date: 20170103 %%%%%%%%%%%% %%%%%% Perform initial freezing temperature for each soil type.%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - SoilVariables = init.applySoilHeteroWithInitialFreezing(LatentHeatOfFreezing, SoilVariables); + SoilVariables = init.applySoilHeteroWithInitialFreezing(LatentHeatOfFreezing, SoilVariables, ModelSettings); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Perform initial thermal calculations for each soil type.%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - ThermalConductivity = init.calculateInitialThermal(SoilVariables, VanGenuchten); + ThermalConductivity = init.calculateInitialThermal(SoilVariables, VanGenuchten, ModelSettings); % SoilVariables will be updated in UpdateSoilWaterContent - SoilVariables = UpdateSoilWaterContent(KIT, LatentHeatOfFreezing, SoilVariables, VanGenuchten); + SoilVariables = UpdateSoilWaterContent(KIT, LatentHeatOfFreezing, SoilVariables, VanGenuchten, ModelSettings); Theta_L = SoilVariables.Theta_L; Theta_I = SoilVariables.Theta_I; @@ -35,9 +35,6 @@ Theta_UU = SoilVariables.Theta_UU; Theta_II = SoilVariables.Theta_II; - % get model settings - ModelSettings = io.getModelSettings(); - for i = 1:ModelSettings.NL Theta_L(i, 1) = Theta_LL(i, 1); Theta_L(i, 2) = Theta_LL(i, 2); diff --git a/src/UpdateSoilWaterContent.m b/src/UpdateSoilWaterContent.m index 813a6907..6acd3218 100644 --- a/src/UpdateSoilWaterContent.m +++ b/src/UpdateSoilWaterContent.m @@ -1,11 +1,8 @@ -function SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten) +function SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten, ModelSettings) %{ Update SoilWaterContent i.e. Theta_LL. %} - % get model settings - ModelSettings = io.getModelSettings(); - Theta_s = VanGenuchten.Theta_s; COR = []; @@ -31,7 +28,7 @@ SoilVariables.hh(MN) = hhU(MN); end - SoilVariables = conductivity.calculateHydraulicConductivity(SoilVariables, VanGenuchten, KIT, L_f); + SoilVariables = conductivity.calculateHydraulicConductivity(SoilVariables, VanGenuchten, KIT, L_f, ModelSettings); Theta_LL = SoilVariables.Theta_LL; DTheta_LLh = SoilVariables.DTheta_LLh; Theta_UU = SoilVariables.Theta_UU; diff --git a/src/ebal.m b/src/ebal.m index d7af30c8..cbb27dce 100644 --- a/src/ebal.m +++ b/src/ebal.m @@ -1,7 +1,7 @@ function [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential] ... = ebal(iter, options, spectral, rad, gap, leafopt, ... angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ... - Rl, SoilVariables, VanGenuchten, InitialValues, GroundwaterSettings) + Rl, SoilVariables, VanGenuchten, InitialValues, ModelSettings, GroundwaterSettings) %{ function ebal.m calculates the energy balance of a vegetated surface @@ -93,8 +93,6 @@ resistance of leaves (or biochemical_MD12: alternative) %} %% 1. initialisations and other preparations for the iteration loop - ModelSettings = io.getModelSettings(); - counter = 0; % Iteration counter of ebal maxit = iter.maxit; maxEBer = iter.maxEBer; @@ -173,7 +171,7 @@ resistance of leaves (or biochemical_MD12: alternative) LAI = canopy.LAI; PSI = 0; - [bbx] = Max_Rootdepth(InitialValues.bbx); + [bbx] = Max_Rootdepth(InitialValues.bbx, ModelSettings); [PSIs, rsss, rrr, rxx] = calc_rsoil(Rl, ModelSettings, SoilVariables, VanGenuchten, bbx, GroundwaterSettings); [sfactor] = calc_sfactor(Rl, VanGenuchten.Theta_s, VanGenuchten.Theta_r, SoilVariables.Theta_LL, bbx, Ta, VanGenuchten.Theta_f); PSIss = PSIs(ModelSettings.NL, 1); diff --git a/src/updateWettingHistory.m b/src/updateWettingHistory.m index f80db6c8..8797460f 100644 --- a/src/updateWettingHistory.m +++ b/src/updateWettingHistory.m @@ -1,4 +1,4 @@ -function XWRE = updateWettingHistory(SoilVariables, VanGenuchten) +function XWRE = updateWettingHistory(SoilVariables, VanGenuchten, ModelSettings) %{ This subroutine is caled after one time step to update the wetting history. If the change in average moisture content of the element during the @@ -7,8 +7,6 @@ scanning curves are used, subject to the constraint that matric head and moisture content be continuous in time. %} - % get model settings - ModelSettings = io.getModelSettings(); XOLD = SoilVariables.XOLD; Theta_L = SoilVariables.Theta_L;