Add a new csv to import soil layers settings

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)
-
-
+
 <a name="1.5.0"></a>
 # [1.5.0](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.5.0) - 3 Jan 2024
 

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;"
-```
+```

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`.
+This creates a log file per each forcing file in the folder `slurm`.

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,

src/+conductivity/+hydraulicConductivity/calcuulateDTheta_LLh.m → 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

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

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();

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;

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();

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

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

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();
 

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

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

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

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();
 

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;

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

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

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;

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;