Merge pull request #158 from EcoExtreML/fix_issue_95

Refactor Constants.m
EcoExtreML · May 30, 2023 · a6cda26 · a6cda26
2 parents 2d64e58 + 5d0669a
commit a6cda26
Show file tree

Hide file tree

Showing 46 changed files with 1,072 additions and 1,048 deletions.
diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE
diff --git a/src/+equations/zo_and_d.m b/src/+equations/zo_and_d.m
@@ -33,9 +33,9 @@
     %   d           zero plane displacement (m)
     %
 
-    %% constants
-    global constants
-    kappa   = constants.kappa;
+    % load Constantns
+    Constants = io.define_constants();
+    kappa = Constants.kappa;
 
     %% parameters
     CR      = canopy.CR;

diff --git a/src/+helpers/createStructure.m b/src/+helpers/createStructure.m
@@ -0,0 +1,7 @@
+function structure = createStructure(values, fields)
+    %{
+        Create a structure with names in fields and values in values
+    %}
+    rep_values = repmat({values}, 1, numel(fields));
+    structure = cell2struct(rep_values, fields, 2);
+end
diff --git a/src/+init/applySoilHeteroEffect.m b/src/+init/applySoilHeteroEffect.m
@@ -1,15 +1,18 @@
-function [SoilVariables, VanGenuchten] = applySoilHeteroEffect(SoilProperties, SoilConstants, SoilVariables, VanGenuchten)
+function [SoilVariables, VanGenuchten] = applySoilHeteroEffect(SoilProperties, SoilConstants, SoilData, SoilVariables, VanGenuchten)
 
-    initX = SoilConstants.InitialValues.initX;
-    initND = SoilConstants.InitialValues.initND;
-    Eqlspace = SoilConstants.Eqlspace;
+    initX = SoilData.InitialValues.initX;
+    initND = SoilData.InitialValues.initND;
 
     SoilVariables.Phi_s = []; % see issue 139, item 5
     SoilVariables.Lamda = []; % see issue 139, item 5
     ImpedF = repelem(3, 6);
 
+    % get model settings
+    ModelSettings = io.getModelSettings();
+    Eqlspace = ModelSettings.Eqlspace;
+
     for i = 1:6
-        if SoilConstants.SWCC == 0
+        if ModelSettings.SWCC == 0
             if SoilConstants.CHST == 0
                 SoilVariables.Phi_s(i) = SoilConstants.Phi_S(i);
                 SoilVariables.Lamda(i) = SoilProperties.Coef_Lamda(i);
@@ -18,29 +21,29 @@
     end
 
     if ~Eqlspace
-        j = SoilConstants.J;
+        j = ModelSettings.J;
         for i = 1:length(initX)
-            SoilConstants.InitialValues.initH(i) = init.calcInitH(VanGenuchten.Theta_s(j), VanGenuchten.Theta_r(j), initX(i), VanGenuchten.n(j), VanGenuchten.m(j), VanGenuchten.Alpha(j));
+            SoilData.InitialValues.initH(i) = init.calcInitH(VanGenuchten.Theta_s(j), VanGenuchten.Theta_r(j), initX(i), VanGenuchten.n(j), VanGenuchten.m(j), VanGenuchten.Alpha(j));
         end
         Dmark = [];
-        for i = 1:SoilConstants.totalNumberOfElements % NL
-            SoilConstants.Elmn_Lnth = SoilConstants.Elmn_Lnth + SoilConstants.DeltZ(i);
-            InitLnth(i) = SoilConstants.Tot_Depth - SoilConstants.Elmn_Lnth;
+        for i = 1:ModelSettings.NL
+            SoilConstants.Elmn_Lnth = SoilConstants.Elmn_Lnth + ModelSettings.DeltZ(i);
+            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, SoilConstants, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i);
-                    SoilConstants.InitialValues.initH = initH;
+                    [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilConstants, SoilData, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i);
+                    SoilData.InitialValues.initH = initH;
                     Dmark = i + 2;
                 end
             end
             if abs(InitLnth(i)) < 1e-10
                 subRoutine = 0;
-                [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilConstants, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i);
-                SoilConstants.InitialValues.initH = initH;
+                [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilConstants, SoilData, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i);
+                SoilData.InitialValues.initH = initH;
             end
         end
     else
-        for i = 1:SoilConstants.numberOfNodes
+        for i = 1:ModelSettings.NN
             SoilVariables.h(i) = -95;
             SoilVariables.T(i) = 22;
             SoilVariables.TT(i) = SoilVariables.T(i);

diff --git a/src/+init/applySoilHeteroWithInitialFreezing.m b/src/+init/applySoilHeteroWithInitialFreezing.m
@@ -2,7 +2,10 @@
 
     SoilVariables.ISFT = 0;
 
-    for i = 1:SoilConstants.numberOfNodes
+    % get model settings
+    ModelSettings = io.getModelSettings();
+
+    for i = 1:ModelSettings.NN
         SoilVariables.h_frez = init.updateHfreez(i, SoilVariables, SoilConstants);
         SoilVariables.hh_frez(i) = SoilVariables.h_frez(i);
         SoilVariables.h(i) = SoilVariables.h(i) - SoilVariables.h_frez(i);
@@ -12,14 +15,14 @@
 
         [SoilVariables.Gama_h, SoilVariables.Gama_hh] = init.updateGamaH(i, SoilConstants, SoilVariables);
 
-        if SoilConstants.Thmrlefc == 1
+        if ModelSettings.Thmrlefc == 1
             SoilVariables.TT(i) = SoilVariables.T(i);
         end
-        if SoilConstants.Soilairefc == 1
+        if ModelSettings.Soilairefc == 1
             SoilConstants.P_g(i) = 95197.850;
             SoilConstants.P_gg(i) = SoilConstants.P_g(i);
         end
-        if i < SoilConstants.numberOfNodes
+        if i < ModelSettings.NN
             SoilVariables.XWRE(i, 1) = 0;
             SoilVariables.XWRE(i, 2) = 0;
         end

diff --git a/src/+init/calculateInitialThermal.m b/src/+init/calculateInitialThermal.m
@@ -25,8 +25,11 @@
     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:SoilConstants.totalNumberOfElements % NL
+    for j = 1:ModelSettings.NL
         % and the sums in the dry thermal conductivity;
         S1(j) = SoilVariables.POR(j) * TCA;
         S2(j) = SoilVariables.POR(j);
@@ -56,8 +59,8 @@
         GB2(j) = (GA1 - GA2) / SoilVariables.XWILT(j) + GB1(j);
 
         %%%%%%%% Johansen thermal conductivity method %%%%%%%
-        RHo_bulk(j) = (1 - VanGenuchten.Theta_s(j)) * 2.7 * 1000;         % Unit g.cm^-3
-        TCON_dry(j) = (0.135 * RHo_bulk(j) + 64.7) / (2700 - 0.947 * RHo_bulk(j));   % Unit W m-1 K-1 ==> j cm^-1 s^-1 Cels^-1
+        RHO_bulk(j) = (1 - VanGenuchten.Theta_s(j)) * 2.7 * 1000;         % Unit g.cm^-3
+        TCON_dry(j) = (0.135 * RHO_bulk(j) + 64.7) / (2700 - 0.947 * RHO_bulk(j));   % Unit W m-1 K-1 ==> j cm^-1 s^-1 Cels^-1
 
         %%%%%%%% organic thermal conductivity method %%%%%%%
         TCON_Soc = 0.05;
@@ -95,7 +98,7 @@
         TPS2(j) = TPS2(j) + TTERM(j);
 
         % Farouki thermal parameter method
-        if SoilConstants.ThermCond == 4
+        if ModelSettings.ThermCond == 4
             FEHCAP(j) = (2.128 * Theta_sa(j) + 2.385 * Theta_cl(j)) / (Theta_sa(j) + Theta_cl(j)) * 1e6;  % j m-3 K-1
             FEHCAP(j) = FEHCAP(j) * (1 - SoilVariables.XSOC(j)) + SoilVariables.XSOC(j) * 2.5 * 1e6;  % organic effect j m-3 K-1
             TCON_s(j) = (8.8 * Theta_sa(j) + 2.92 * Theta_cl(j)) / (Theta_sa(j) + Theta_cl(j)); % W m-1 K-1
@@ -110,7 +113,7 @@
     ThermalConductivity.GA2 = GA2;
     ThermalConductivity.GB1 = GB1;
     ThermalConductivity.GB2 = GB2;
-    ThermalConductivity.RHo_bulk = RHo_bulk;
+    ThermalConductivity.RHO_bulk = RHO_bulk;
     ThermalConductivity.TCON_dry = TCON_dry;
     ThermalConductivity.S1 = S1;
     ThermalConductivity.S2 = S2;