Merge branch 'main' into add_csv_for_soil_layers_settings

EcoExtreML · Aug 14, 2024 · 18c39d3 · 18c39d3
2 parents b1069d9 + d4ddc41
commit 18c39d3
Show file tree

Hide file tree

Showing 8 changed files with 19 additions and 16 deletions.
diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE
diff --git a/src/+energy/assembleCoefficientMatrices.m b/src/+energy/assembleCoefficientMatrices.m
@@ -1,4 +1,4 @@
-function [RHS, EnergyMatrices, SAVE] = assembleCoefficientMatrices(EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, ModelSettings, 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
@@ -29,9 +29,10 @@
 
     % Alias of SoilVariables
     SV = SoilVariables;
+    RHS = InitialValues.RHS;
 
     if ModelSettings.Soilairefc && ModelSettings.Thmrlefc
-        RHS(indxBotm) = -C7(indxBotm) + (C2(indxBotm, 1) * SV.T(indxBotm) + C2(indxBotm, 2) * SV.T(2)) / Delt_t ...
+        RHS(indxBotm) = -C7(indxBotm) + (C2(indxBotm, 1) * SV.T(indxBotm) + C2(indxBotm, 2) * SV.T(indxBotm + 1)) / Delt_t ...
             - (C1(indxBotm, 1) / Delt_t + C4(indxBotm, 1)) * SV.hh(indxBotm) - (C1(indxBotm, 2) / Delt_t + C4(indxBotm, 2)) * SV.hh(indxBotm + 1) ...
             - (C3(indxBotm, 1) / Delt_t + C6(indxBotm, 1)) * P_gg(indxBotm) - (C3(indxBotm, 2) / Delt_t + C6(indxBotm, 2)) * P_gg(indxBotm + 1) ...
             + (C3(indxBotm, 1) / Delt_t) * P_g(indxBotm) + (C3(indxBotm, 2) / Delt_t) * P_g(indxBotm + 1) ...
@@ -78,7 +79,7 @@
             + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
 
     else
-        RHS(indxBotm) = -C7(indxBotm) + (C2(indxBotm, 1) * SV.T(indxBotm) + C2(1, 2) * SV.T(indxBotm + 1)) / Delt_t;
+        RHS(indxBotm) = -C7(indxBotm) + (C2(indxBotm, 1) * SV.T(indxBotm) + C2(indxBotm, 2) * SV.T(indxBotm + 1)) / Delt_t;
         for i = indxBotm + 1:ModelSettings.NL
             RHS(i) = -C7(i) + (C2(i - 1, 2) * SV.T(i - 1) + C2(i, 1) * SV.T(i) + C2(i, 2) * SV.T(i + 1)) / Delt_t;
         end

diff --git a/src/+energy/calculateBoundaryConditions.m b/src/+energy/calculateBoundaryConditions.m
@@ -15,10 +15,6 @@
         % index of bottom layer after neglecting saturated layers (from bottom to top)
         indxBotm = GroundwaterSettings.indxBotmLayer;
         tempBotm = GroundwaterSettings.tempBotm; % groundwater temperature
-
-        if isnan(GroundwaterSettings.tempBotm) % no available groundwater temperature data
-            tempBotm = SoilVariables.TT(indxBotm);
-        end
     end
 
     % Apply the bottom boundary condition called for by BoundaryCondition.NBCTB

diff --git a/src/+energy/solveEnergyBalanceEquations.m b/src/+energy/solveEnergyBalanceEquations.m
@@ -16,7 +16,7 @@
 
     EnergyMatrices = energy.calculateMatricCoefficients(EnergyVariables, InitialValues, ModelSettings, GroundwaterSettings);
 
-    [RHS, EnergyMatrices, SAVE] = energy.assembleCoefficientMatrices(EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, ModelSettings, 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, ModelSettings, GroundwaterSettings);

diff --git a/src/+energy/solveTridiagonalMatrixEquations.m b/src/+energy/solveTridiagonalMatrixEquations.m
@@ -26,8 +26,6 @@
         for i = indxBotm:-1:2
             RHS(i - 1) = RHS(i); % assign the groundwater temperature to the saturated layers
         end
-        % correct soil temp above groundwater table (avoid wrong temp behaviour in case RHS(indxBotm) > RHS(indxBotm + 1))
-        % RHS(indxBotm + 1) = (RHS(indxBotm) + RHS(indxBotm + 2)) / 2;
     end
 
     for i = 1:ModelSettings.NN

diff --git a/src/+groundwater/calculateGroundWaterDepth.m b/src/+groundwater/calculateGroundWaterDepth.m
@@ -6,8 +6,5 @@
     if groundWaterDepth <= 0
         warning('The soil is fully saturated up to the land surface level!');
         groundWaterDepth = 1.0; % to avoid model crashing, assign minimum groundWaterDepth value of 1 cm
-    elseif groundWaterDepth > Tot_Depth
-        warning('Groundwater table is below the end of the soil column!');
     end
-
 end
diff --git a/src/+soilmoisture/calculateTimeDerivatives.m b/src/+soilmoisture/calculateTimeDerivatives.m
@@ -60,7 +60,7 @@
     elseif ModelSettings.Thmrlefc && ModelSettings.Soilairefc
         RHS(indxBotm) = -C9(indxBotm) - C7(indxBotm) + (C1(indxBotm, 1) * h(indxBotm) + C1(indxBotm, 2) * h(indxBotm + 1)) / Delt_t ...
             - (C2(indxBotm, 1) / Delt_t + C5(indxBotm, 1)) * TT(indxBotm) - (C2(indxBotm, 2) / Delt_t + C5(indxBotm, 2)) * TT(indxBotm + 1) ...
-            - C6(indxBotm, 1) * P_gg(indxBotm) - C6(indxBotm, 2) * P_gg(2) ...
+            - C6(indxBotm, 1) * P_gg(indxBotm) - C6(indxBotm, 2) * P_gg(indxBotm + 1) ...
             + (C2(indxBotm, 1) / Delt_t) * T(indxBotm) + (C2(indxBotm, 2) / Delt_t) * T(indxBotm + 1);
         for i = indxBotm + 1:ModelSettings.NL
             ARG1 = C2(i - 1, 2) / Delt_t;

diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m
@@ -352,14 +352,25 @@
     if GroundwaterSettings.GroundwaterCoupling == 1  % Groundwater coupling is enabled
         BoundaryCondition.NBChB = 1;
 
-        % update GroundwaterSettings.headBotmLayer and GroundwaterSettings.tempBotm, from MODFLOW through BMI
+        % Update GroundwaterSettings.headBotmLayer and GroundwaterSettings.tempBotm, from MODFLOW through BMI
         GroundwaterSettings.gw_Dep = groundwater.calculateGroundWaterDepth(GroundwaterSettings.topLevel, GroundwaterSettings.headBotmLayer, ModelSettings.Tot_Depth);
 
-        % update Dunnian runoff and ForcingData.Precip_msr
+        % Update Dunnian runoff and ForcingData.Precip_msr
         [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, ModelSettings);
+
+        % Check the position of the groundwater table
+        if any(GroundwaterSettings.gw_Dep > ModelSettings.Tot_Depth) || any(isempty(GroundwaterSettings.indxBotmLayer))
+            warning('Groundwater table is below the end of the soil column. Please enlarge the total soil thickness!');
+            return
+        end
+
+        % Check if no available groundwater temperature data
+        if isnan(GroundwaterSettings.tempBotm)
+            GroundwaterSettings.tempBotm = SoilVariables.TT(GroundwaterSettings.indxBotmLayer);
+        end
     end
 
     % Will do one timestep in "update mode", and run until the end if in "full run" mode.