Add a new csv to import soil layers settings

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 = calcuulateDTheta_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;

src/+dryair/solveDryAirEquations.m

@@ -1,20 +1,20 @@
 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
         exampled as for the soil moisture equation, which derived in 'STEMMUS
         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

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

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

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

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

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

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;