Merge pull request #180 from EcoExtreML/fix_issue_96

Fix issue 96
EcoExtreML · Jun 20, 2023 · 1918938 · 1918938
2 parents 750addd + e8a8c77
commit 1918938
Show file tree

Hide file tree

Showing 24 changed files with 482 additions and 462 deletions.
diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE
diff --git a/src/+init/applySoilHeteroEffect.m b/src/+init/applySoilHeteroEffect.m
@@ -1,7 +1,7 @@
-function [SoilVariables, VanGenuchten] = applySoilHeteroEffect(SoilProperties, SoilConstants, SoilData, SoilVariables, VanGenuchten)
+function [SoilVariables, VanGenuchten] = applySoilHeteroEffect(SoilProperties, SoilVariables, VanGenuchten)
 
-    initX = SoilData.InitialValues.initX;
-    initND = SoilData.InitialValues.initND;
+    initX = SoilVariables.InitialValues.initX;
+    initND = SoilVariables.InitialValues.initND;
 
     SoilVariables.Phi_s = []; % see issue 139, item 5
     SoilVariables.Lamda = []; % see issue 139, item 5
@@ -11,6 +11,9 @@
     ModelSettings = io.getModelSettings();
     Eqlspace = ModelSettings.Eqlspace;
 
+    % Get soil constants for StartInit
+    SoilConstants = io.getSoilConstants();
+
     for i = 1:6
         if ModelSettings.SWCC == 0
             if SoilConstants.CHST == 0
@@ -23,23 +26,23 @@
     if ~Eqlspace
         j = ModelSettings.J;
         for i = 1:length(initX)
-            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));
+            SoilVariables.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: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, SoilData, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i);
-                    SoilData.InitialValues.initH = initH;
+                    [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, 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, SoilConstants, SoilData, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i);
-                SoilData.InitialValues.initH = initH;
+                [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i);
+                SoilVariables.InitialValues.initH = initH;
             end
         end
     else

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

diff --git a/src/+init/calculateInitialThermal.m b/src/+init/calculateInitialThermal.m
@@ -1,4 +1,4 @@
-function ThermalConductivity = calculateInitialThermal(SoilConstants, SoilVariables, VanGenuchten)
+function ThermalConductivity = calculateInitialThermal(SoilVariables, VanGenuchten)
     FEHCAP = []; % see issue 139
 
     HCAP(1) = 0.998 * 4.182;

diff --git a/src/+init/defineSoilVariables.m b/src/+init/defineSoilVariables.m
@@ -0,0 +1,38 @@
+function SoilVariables = defineSoilVariables(InitialValues, SoilProperties, VanGenuchten)
+    %{
+        Create a SoilVariables structure and define some variables.
+        SoilVariables is input for initializing Temperature, Matric potential
+        and soil air pressure.
+    %}
+
+    % Create a SoilVariables structure
+    SoilVariables = struct();
+
+    % Add initial values to SoilVariables
+    soil_fields = {
+                   'P_g', 'P_gg', 'h', 'T', 'TT', 'h_frez', 'Theta_L', ...
+                   'Theta_LL', 'Theta_V', 'Theta_g', 'Se', 'KL_h', ...
+                   'DTheta_LLh', 'KfL_T', 'Theta_II', 'Theta_I', ...
+                   'Theta_UU', 'Theta_U', 'TT_CRIT', 'KfL_h', 'DTheta_UUh'
+                  };
+    for field = soil_fields
+        SoilVariables.(field{1}) = InitialValues.(field{1});
+    end
+
+    % Calculate some of SoilVariables
+    SoilVariables.XK = VanGenuchten.Theta_r + 0.02; % 0.11 This is for silt loam; For sand XK=0.025
+    SoilVariables.XWILT = equations.van_genuchten(VanGenuchten.Theta_s, VanGenuchten.Theta_r, VanGenuchten.Alpha, -1.5e4, VanGenuchten.n, VanGenuchten.m);
+    SoilVariables.XCAP = equations.van_genuchten(VanGenuchten.Theta_s, VanGenuchten.Theta_r, VanGenuchten.Alpha, -336, VanGenuchten.n, VanGenuchten.m);
+
+    SoilVariables.VPERSOC = equations.calc_msoc_fraction(SoilProperties.MSOC); % fraction of soil organic matter
+    SoilVariables.FOSL = 1 - SoilProperties.FOC - SoilProperties.FOS - SoilVariables.VPERSOC;
+
+    % Add variables to SoilVariables, see issue 141
+    SoilVariables.POR = SoilProperties.porosity;
+    SoilVariables.VPERS = SoilProperties.FOS' .* (1 - SoilVariables.POR);
+    SoilVariables.VPERSL = SoilVariables.FOSL' .* (1 - SoilVariables.POR);
+    SoilVariables.VPERC = SoilProperties.FOC' .* (1 - SoilVariables.POR);
+    SoilVariables.Ks = SoilProperties.SaturatedK;
+    SoilVariables.Vol_qtz = SoilProperties.FOS;
+
+end
diff --git a/src/+init/setBoundaryCondition.m b/src/+init/setBoundaryCondition.m
@@ -1,4 +1,4 @@
-function BoundaryCondition = setBoundaryCondition(SoilVariables, SoilConstants, initialLandcoverClass)
+function BoundaryCondition = setBoundaryCondition(SoilVariables, ForcingData, initialLandcoverClass)
     % Variables information for boundary condition settings
     % NBCh   Indicator for type of surface boundary condition on mass euqation to be applied;
     %        "1"--Specified matric head;
@@ -47,7 +47,7 @@
     IRPT1 = 0;
     IRPT2 = 0;
 
-    Ta_msr = SoilConstants.Ta_msr;
+    Ta_msr = ForcingData.Ta_msr;
 
     % get model settings
     ModelSettings = io.getModelSettings();

diff --git a/src/+init/setSoilConstants.m b/src/+init/setSoilConstants.m
diff --git a/src/+init/setSoilVariables.m b/src/+init/setSoilVariables.m
diff --git a/src/+init/soilHeteroSubroutine.m b/src/+init/soilHeteroSubroutine.m
@@ -1,11 +1,11 @@
-function [SoilVariables, VanGenuchten, initH] = soilHeteroSubroutine(subRoutine, SoilConstants, SoilData, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, ML)
+function [SoilVariables, VanGenuchten, initH] = soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, ML)
     %{
         considering soil heterogeneity effect
     %}
 
-    initX = SoilData.InitialValues.initX;
-    initT = SoilData.InitialValues.initT;
-    initH = SoilData.InitialValues.initH;
+    initX = SoilVariables.InitialValues.initX;
+    initT = SoilVariables.InitialValues.initT;
+    initH = SoilVariables.InitialValues.initH;
 
     % get model settings
     ModelSettings = io.getModelSettings();
@@ -56,14 +56,14 @@
         end
 
         J = SoilVariables.IS(i);
-        [SoilVariables, VanGenuchten] = init.updateSoilVariables(SoilVariables, VanGenuchten, SoilConstants, SoilProperties, j, J);
+        [SoilVariables, VanGenuchten] = init.updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, j, J);
         SoilVariables.Imped(i) = ImpedF(J);
 
-        initH(indexOfInit) = init.updateInitH(initX(indexOfInit), VanGenuchten, SoilConstants, SoilVariables, j);
+        initH(indexOfInit) = init.updateInitH(initX(indexOfInit), VanGenuchten, SoilVariables, j);
 
         if subRoutine == 5
-            Btmh = init.updateBtmh(VanGenuchten, SoilData, SoilVariables, i);
-            SoilVariables.T(i) = SoilData.BtmT + (i - 1) * (initT(indexOfInit) - SoilData.BtmT) / ML;
+            Btmh = init.updateBtmh(VanGenuchten, SoilVariables, i);
+            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
         else

diff --git a/src/+init/updateBtmh.m b/src/+init/updateBtmh.m
@@ -1,10 +1,10 @@
-function Btmh = updateBtmh(VanGenuchten, SoilData, SoilVariables, i)
+function Btmh = updateBtmh(VanGenuchten, SoilVariables, i)
     % get model settings
     ModelSettings = io.getModelSettings();
 
     if ModelSettings.SWCC == 1   % VG soil water retention model
-        Btmh = init.calcInitH(VanGenuchten.Theta_s(i), VanGenuchten.Theta_r(i), SoilData.BtmX, VanGenuchten.n(i), VanGenuchten.m(i), VanGenuchten.Alpha(i));
+        Btmh = init.calcInitH(VanGenuchten.Theta_s(i), VanGenuchten.Theta_r(i), SoilVariables.BtmX, VanGenuchten.n(i), VanGenuchten.m(i), VanGenuchten.Alpha(i));
     else
-        Btmh = SoilVariables.Phi_s(i) * (SoilData.BtmX / VanGenuchten.Theta_s(i))^(-1 / SoilVariables.Lamda(i));
+        Btmh = SoilVariables.Phi_s(i) * (SoilVariables.BtmX / VanGenuchten.Theta_s(i))^(-1 / SoilVariables.Lamda(i));
     end
 end
diff --git a/src/+init/updateGamaH.m b/src/+init/updateGamaH.m
@@ -1,4 +1,7 @@
-function [Gama_h, Gama_hh] = updateGamaH(i, SoilConstants, SoilVariables)
+function [Gama_h, Gama_hh] = updateGamaH(i, SoilVariables)
+
+    % get soil constants
+    SoilConstants = io.getSoilConstants();
 
     hd = SoilConstants.hd;
     hm = SoilConstants.hm;

diff --git a/src/+init/updateHfreez.m b/src/+init/updateHfreez.m
@@ -1,10 +1,7 @@
-function h_frez = updateHfreez(i, SoilVariables, SoilConstants)
+function h_frez = updateHfreez(i, LatentHeatOfFreezing, SoilVariables)
     % get Constants
     Constants = io.define_constants();
 
-    L_f = 3.34 * 1e5; % latent heat of freezing fusion J Kg-1
-    T0 = 273.15; % unit K
-
     T = SoilVariables.T;
     h_frez = SoilVariables.h_frez;
     h = SoilVariables.h;
@@ -14,7 +11,7 @@
     ModelSettings = io.getModelSettings();
 
     if T(i) <= 0
-        h_frez(i) = L_f * 1e4 * (T(i)) / Constants.g / T0;
+        h_frez(i) = LatentHeatOfFreezing * 1e4 * (T(i)) / Constants.g / ModelSettings.T0;
     else
         h_frez(i) = 0;
     end

diff --git a/src/+init/updateInitH.m b/src/+init/updateInitH.m
@@ -1,4 +1,4 @@
-function initH = updateInitH(initX, VanGenuchten, SoilConstants, SoilVariables, i)
+function initH = updateInitH(initX, VanGenuchten, SoilVariables, i)
     % get model settings
     ModelSettings = io.getModelSettings();
 

diff --git a/src/+init/updateSoilVariables.m b/src/+init/updateSoilVariables.m
@@ -1,12 +1,15 @@
-function [SoilVariables, VanGenuchten] = updateSoilVariables(SoilVariables, VanGenuchten, SoilConstants, SoilProperties, i, j)
+function [SoilVariables, VanGenuchten] = updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, i, j)
+
+    % Get soil constants for StartInit
+    SoilConstants = io.getSoilConstants();
 
     SoilVariables.POR(i) = SoilProperties.porosity(j);
     SoilVariables.Ks(i) = SoilProperties.SaturatedK(j);
-    SoilVariables.Theta_qtz(i) = SoilConstants.Vol_qtz(j);
+    SoilVariables.Theta_qtz(i) = SoilVariables.Vol_qtz(j);
     SoilVariables.VPER(i, 1) = SoilVariables.VPERS(j);
     SoilVariables.VPER(i, 2) = SoilVariables.VPERSL(j);
     SoilVariables.VPER(i, 3) = SoilVariables.VPERC(j);
-    SoilVariables.XSOC(i) = SoilConstants.VPERSOC(j);
+    SoilVariables.XSOC(i) = SoilVariables.VPERSOC(j);
     SoilVariables.XK(i) = SoilConstants.XK;
 
     % get model settings

diff --git a/src/+io/getModelSettings.m b/src/+io/getModelSettings.m
@@ -59,7 +59,6 @@
     ModelSettings.NS = 1; % Number of soil types;
 
     % The time and domain information setting
-    ModelSettings.KIT = 0; % KIT is used to count the number of iteration in a time step;
     ModelSettings.NIT = 30; % Desirable number of iterations in a time step;
     ModelSettings.KT = 0; % Number of time steps;
 

diff --git a/src/+io/getSoilConstants.m b/src/+io/getSoilConstants.m
@@ -0,0 +1,16 @@
+function SoilConstants = getSoilConstants()
+
+    SoilConstants.hd = -1e7;
+    SoilConstants.hm = -9899;
+    SoilConstants.CHST = 0;
+    SoilConstants.Elmn_Lnth = 0;
+    SoilConstants.Dmark = 0;
+    SoilConstants.Phi_S = [-17.9 -17 -17 -19 -10 -10];
+    SoilConstants.Phi_soc = -0.0103;
+    SoilConstants.Lamda_soc = 2.7;
+    SoilConstants.Theta_soc = 0.6;
+    % XK=0.11 for silt loam; For sand XK=0.025
+    % SoilConstants.XK is used in updateSoilVariables
+    SoilConstants.XK = 0.11;
+
+end