Update the code to count for groundwater (SSM v.0.3.2)

CHANGELOG.md

@@ -3,10 +3,23 @@
 **Changed:**
 
 - Added changes to support groundwater coupling via BMI in
-  [#221](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/221)
+  [#221](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/221) and
+  [#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)
 
+**Fixed:**
+
+- Calculations of surface runoff discussed in
+  [#232](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/232) and fixed in
+  [#234](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/234)
+- The bug in the QVT calculations discussed in
+  [#230](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/230) and fixed in
+  [#234](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/234)
+- The bug in  activating the dry air calculations discussed in
+  [#227](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/230) and fixed in
+  [#234](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/227)
+
 
 <a name="1.5.0"></a>
 # [1.5.0](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.5.0) - 3 Jan 2024

src/+conductivity/calculateVaporVariables.m

@@ -41,6 +41,11 @@
 
                 if SoilVariables.DVT_Switch == 1
                     Eta(i, j) = ThermalConductivityCapacity.ZETA(i, j) * EnhnLiqIsland / (f0 * Theta_g(i, j));
+                    % When Theta_g(i, j) = 0, both EnhnLiqIsland and f0 have too low values -> leading to extreme large value of Eta and QVT, ...
+                    % The following line (if statement) tries to solve the issue mentioned in https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/230
+                    if Theta_g(i, j) <= 1e-2 % added by Mostafa
+                        Eta(i, j) = 0;
+                    end
                 else
                     Eta(i, j) = 0;
                 end

src/+dryair/assembleCoefficientMatrices.m

@@ -1,4 +1,4 @@
-function [RHS, AirMatrices, SAVE] = assembleCoefficientMatrices(AirMatrices, SoilVariables, Delt_t, P_g)
+function [RHS, AirMatrices, SAVE] = assembleCoefficientMatrices(AirMatrices, SoilVariables, Delt_t, P_g, GroundwaterSettings)
     %{
         Assemble the coefficient matrices of Equation 4.32 STEMMUS Technical
         Notes, page 44, for dry air equation.
@@ -20,14 +20,21 @@
     % Alias of SoilVariables
     SV = SoilVariables;
 
+    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
+        % index of bottom layer after neglecting saturated layers (from bottom to top)
+        indxBotm = GroundwaterSettings.indxBotmLayer;
+    end
+
     if ModelSettings.Thmrlefc
-        RHS(1) = -C7(1) + (C3(1, 1) * P_g(1) + C3(1, 2) * P_g(2)) / Delt_t ...
-            - (C2(1, 1) / Delt_t + C5(1, 1)) * SV.TT(1) - (C2(1, 2) / Delt_t + C5(1, 2)) * SV.TT(2) ...
-            - (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
-            + (C2(1, 1) / Delt_t) * SV.T(1) + (C2(1, 2) / Delt_t) * SV.T(2) ...
-            + (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
+        RHS(indxBotm) = -C7(indxBotm) + (C3(indxBotm, 1) * P_g(indxBotm) + C3(indxBotm, 2) * P_g(indxBotm + 1)) / Delt_t ...
+            - (C2(indxBotm, 1) / Delt_t + C5(indxBotm, 1)) * SV.TT(indxBotm) - (C2(indxBotm, 2) / Delt_t + C5(indxBotm, 2)) * SV.TT(indxBotm + 1) ...
+            - (C1(indxBotm, 1) / Delt_t + C4(indxBotm, 1)) * SV.hh(indxBotm) - (C1(indxBotm, 2) / Delt_t + C4(indxBotm, 2)) * SV.hh(indxBotm + 1) ...
+            + (C2(indxBotm, 1) / Delt_t) * SV.T(indxBotm) + (C2(indxBotm, 2) / Delt_t) * SV.T(indxBotm + 1) ...
+            + (C1(indxBotm, 1) / Delt_t) * SV.h(indxBotm) + (C1(indxBotm, 2) / Delt_t) * SV.h(indxBotm + 1);
 
-        for i = 2:ModelSettings.NL
+        for i = indxBotm + 1:ModelSettings.NL
             ARG1 = C2(i - 1, 2) / Delt_t;
             ARG2 = C2(i, 1) / Delt_t;
             ARG3 = C2(i, 2) / Delt_t;
@@ -49,10 +56,10 @@
             + (C2(n - 1, 2) / Delt_t) * SV.T(n - 1) + (C2(n, 1) / Delt_t) * SV.T(n) ...
             + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
     else
-        RHS(1) = -C7(1) + (C3(1, 1) * P_g(1) + C3(1, 2) * P_g(2)) / Delt_t ...
-            - (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
-            + (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
-        for i = 2:ModelSettings.NL
+        RHS(indxBotm) = -C7(indxBotm) + (C3(indxBotm, 1) * P_g(indxBotm) + C3(indxBotm, 2) * P_g(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) ...
+            + (C1(indxBotm, 1) / Delt_t) * SV.h(indxBotm) + (C1(indxBotm, 2) / Delt_t) * SV.h(indxBotm + 1);
+        for i = indxBotm + 1:ModelSettings.NL
             ARG4 = C1(i - 1, 2) / Delt_t;
             ARG5 = C1(i, 1) / Delt_t;
             ARG6 = C1(i, 2) / Delt_t;
@@ -65,17 +72,17 @@
             + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
     end
 
-    for i = 1:ModelSettings.NN
+    for i = indxBotm:ModelSettings.NN
         for j = 1:ModelSettings.nD
             C6(i, j) = C3(i, j) / Delt_t + C6(i, j);
         end
     end
 
     AirMatrices.C6 = C6;
 
-    SAVE(1, 1, 3) = RHS(1);
-    SAVE(1, 2, 3) = C6(1, 1);
-    SAVE(1, 3, 3) = C6(1, 2);
+    SAVE(1, 1, 3) = RHS(indxBotm);
+    SAVE(1, 2, 3) = C6(indxBotm, 1);
+    SAVE(1, 3, 3) = C6(indxBotm, 2);
     SAVE(2, 1, 3) = RHS(n);
     SAVE(2, 2, 3) = C6(n - 1, 2);
     SAVE(2, 3, 3) = C6(n, 1);

src/+dryair/calculateBoundaryConditions.m

@@ -1,4 +1,4 @@
-function [RHS, AirMatrices] = calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT)
+function [RHS, AirMatrices] = calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT, P_gg, GroundwaterSettings)
     %{
         Determine the boundary condition for solving the dry air equation.
     %}
@@ -7,15 +7,24 @@
     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
+        BtmPg = BoundaryCondition.BtmPg;
+    else % Groundwater Coupling is activated
+        % index of bottom layer after neglecting saturated layers (from bottom to top)
+        indxBotm = GroundwaterSettings.indxBotmLayer;
+        BtmPg = P_gg(indxBotm);
+    end
+
     % Apply the bottom boundary condition called for by NBCPB
     if BoundaryCondition.NBCPB == 1  % Bounded bottom with the water table
-        RHS(1) = BtmPg;
-        AirMatrices.C6(1, 1) = 1;
-        RHS(2) = RHS(2) - AirMatrices.C6(1, 2) * RHS(1);
-        AirMatrices.C6(1, 2) = 0;
-        AirMatrices.C6_a(1) = 0;
+        RHS(indxBotm) = BtmPg;
+        AirMatrices.C6(indxBotm, 1) = 1;
+        RHS(indxBotm + 1) = RHS(indxBotm + 1) - AirMatrices.C6(indxBotm, 2) * RHS(indxBotm);
+        AirMatrices.C6(indxBotm, 2) = 0;
+        AirMatrices.C6_a(indxBotm) = 0;
     elseif BoundaryCondition.NBCPB == 2  % The soil air is allowed to escape from the bottom
-        RHS(1) = RHS(1) + BoundaryCondition.BCPB;
+        RHS(indxBotm) = RHS(indxBotm) + BoundaryCondition.BCPB;
     end
 
     % Apply the surface boundary condition called by NBCP

src/+dryair/calculateDryAirParameters.m

@@ -1,5 +1,5 @@
 function AirVariabes = calculateDryAirParameters(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, VaporVariables, ...
-                                                 P_gg, Xah, XaT, Xaa, RHODA)
+                                                 P_gg, Xah, XaT, Xaa, RHODA, GroundwaterSettings)
     %{
         Calculate all the parameters related to dry air equation e.g., Equation
         3.59-3.64, STEMMUS Technical Notes, page 27-28.
@@ -24,7 +24,14 @@
     AirVariabes.DTDZ = InitialValues.DTDZ;
     GasDispersivity.DPgDZ = InitialValues.DPgDZ;
 
-    for i = 1:ModelSettings.NL
+    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
+        % index of bottom layer after neglecting saturated layers (from bottom to top)
+        indxBotm = GroundwaterSettings.indxBotmLayer;
+    end
+
+    for i = indxBotm:ModelSettings.NL
         KLhBAR = (SoilVariables.KfL_h(i, 1) + SoilVariables.KfL_h(i, 2)) / 2;
         KLTBAR = (InitialValues.KL_T(i, 1) + InitialValues.KL_T(i, 2)) / 2;
         DDhDZ = (SoilVariables.hh(i + 1) - SoilVariables.hh(i)) / ModelSettings.DeltZ(i);
@@ -50,7 +57,11 @@
         AirVariabes.DDhDZ(i) = DDhDZ;
         AirVariabes.DhDZ(i) = DhDZ;
         AirVariabes.DTDZ(i) = DTDZ;
-        AirVariabes.QL(i) = QL;
+        AirVariabes.QL(i) = QL; % total liquid flux
+        % added by Mostafa
+        AirVariabes.QL_h(i) = QL_h(i); % liquid flux due to matric potential
+        AirVariabes.QL_T(i) = QL_T(i); % liquid flux due to temperature gradient
+        AirVariabes.QL_a(i) = QL_a(I); % liquid flux due to air pressure gradient
 
         for j = 1:ModelSettings.nD
             MN = i + j - 1;

src/+dryair/calculateMatricCoefficients.m

@@ -1,4 +1,4 @@
-function AirMatrices = calculateMatricCoefficients(AirVariabes, InitialValues)
+function AirMatrices = calculateMatricCoefficients(AirVariabes, InitialValues, 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
@@ -15,7 +15,14 @@
     AirMatrices.C6 = InitialValues.C6;
     AirMatrices.C7 = zeros(ModelSettings.NN);
 
-    for i = 1:ModelSettings.NL
+    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
+        % index of bottom layer after neglecting saturated layers (from bottom to top)
+        indxBotm = GroundwaterSettings.indxBotmLayer;
+    end
+
+    for i = indxBotm:ModelSettings.NL
         AirMatrices.C1(i, 1) = AirMatrices.C1(i, 1) + AirVariabes.Cah(i, 1) * ModelSettings.DeltZ(i) / 2;
         AirMatrices.C1(i + 1, 1) = AirMatrices.C1(i + 1, 1) + AirVariabes.Cah(i, 2) * ModelSettings.DeltZ(i) / 2;
 

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)
+                                                               BoundaryCondition, ForcingData, P_gg, P_g, Xah, XaT, Xaa, RHODA, KT, Delt_t, 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);
+                                                   P_gg, Xah, XaT, Xaa, RHODA, GroundwaterSettings);
 
-    AirMatrices = dryair.calculateMatricCoefficients(AirVariabes, InitialValues);
+    AirMatrices = dryair.calculateMatricCoefficients(AirVariabes, InitialValues, GroundwaterSettings);
 
-    [RHS, AirMatrices, SAVE] = dryair.assembleCoefficientMatrices(AirMatrices, SoilVariables, Delt_t, P_g);
+    [RHS, AirMatrices, SAVE] = dryair.assembleCoefficientMatrices(AirMatrices, SoilVariables, Delt_t, P_g, GroundwaterSettings);
 
-    [RHS, AirMatrices] = dryair.calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT);
+    [RHS, AirMatrices] = dryair.calculateBoundaryConditions(BoundaryCondition, AirMatrices, ForcingData, RHS, KT, P_gg, GroundwaterSettings);
 
-    [AirMatrices, P_gg, RHS] = dryair.solveTridiagonalMatrixEquations(RHS, AirMatrices);
+    [AirMatrices, P_gg, RHS] = dryair.solveTridiagonalMatrixEquations(RHS, AirMatrices, GroundwaterSettings);
 
 end

src/+dryair/solveTridiagonalMatrixEquations.m

@@ -1,22 +1,30 @@
-function [AirMatrices, P_gg, RHS] = solveTridiagonalMatrixEquations(RHS, AirMatrices)
+function [AirMatrices, P_gg, RHS] = solveTridiagonalMatrixEquations(RHS, AirMatrices, 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();
-    RHS(1) = RHS(1) / AirMatrices.C6(1, 1);
 
-    for i = 2: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
+        % index of bottom layer after neglecting saturated layers (from bottom to top)
+        indxBotm = GroundwaterSettings.indxBotmLayer;
+    end
+
+    RHS(indxBotm) = RHS(indxBotm) / AirMatrices.C6(indxBotm, 1);
+
+    for i = indxBotm + 1:ModelSettings.NN
         AirMatrices.C6(i, 1) = AirMatrices.C6(i, 1) - AirMatrices.C6_a(i - 1) * AirMatrices.C6(i - 1, 2) / AirMatrices.C6(i - 1, 1);
         RHS(i) = (RHS(i) - AirMatrices.C6_a(i - 1) * RHS(i - 1)) / AirMatrices.C6(i, 1);
     end
 
-    for i = ModelSettings.NL:-1:1
+    for i = ModelSettings.NL:-1:indxBotm
         RHS(i) = RHS(i) - AirMatrices.C6(i, 2) * RHS(i + 1) / AirMatrices.C6(i, 1);
     end
 
-    for i = 1:ModelSettings.NN
+    for i = indxBotm:ModelSettings.NN
         P_gg(i) = RHS(i);
     end
 end

src/+energy/assembleCoefficientMatrices.m

@@ -1,4 +1,4 @@
-function [RHS, EnergyMatrices, SAVE] = assembleCoefficientMatrices(EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg)
+function [RHS, EnergyMatrices, SAVE] = assembleCoefficientMatrices(EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, GroundwaterSettings)
     %{
         assembles the coefficient matrices of Equation 4.32, STEMMUS Technical
         Notes, page 44, the example was only shown for the soil moisture
@@ -17,17 +17,28 @@
     ModelSettings = io.getModelSettings();
     n = ModelSettings.NN;
 
+    if ModelSettings.Soilairefc == 1 % see https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/227
+        C6_a = EnergyMatrices.C6_a;
+    end
+
+    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
+        % index of bottom layer after neglecting saturated layers (from bottom to top)
+        indxBotm = GroundwaterSettings.indxBotmLayer;
+    end
+
     % Alias of SoilVariables
     SV = SoilVariables;
 
     if ModelSettings.Soilairefc && ModelSettings.Thmrlefc
-        RHS(1) = -C7(1) + (C2(1, 1) * SV.T(1) + C2(1, 2) * SV.T(2)) / Delt_t ...
-            - (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
-            - (C3(1, 1) / Delt_t + C6(1, 1)) * P_gg(1) - (C3(1, 2) / Delt_t + C6(1, 2)) * P_gg(2) ...
-            + (C3(1, 1) / Delt_t) * P_g(1) + (C3(1, 2) / Delt_t) * P_g(2) ...
-            + (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
+        RHS(indxBotm) = -C7(indxBotm) + (C2(indxBotm, 1) * SV.T(indxBotm) + C2(indxBotm, 2) * SV.T(2)) / 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) ...
+            + (C1(indxBotm, 1) / Delt_t) * SV.h(indxBotm) + (C1(indxBotm, 2) / Delt_t) * SV.h(indxBotm + 1);
 
-        for i = 2:ModelSettings.NL
+        for i = indxBotm + 1:ModelSettings.NL
             ARG1 = C3(i - 1, 2) / Delt_t;
             ARG2 = C3(i, 1) / Delt_t;
             ARG3 = C3(i, 2) / Delt_t;
@@ -49,9 +60,9 @@
             + (C3(n - 1, 2) / Delt_t) * P_g(n - 1) + (C3(n, 1) / Delt_t) * P_g(n) ...
             + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
     elseif ~ModelSettings.Soilairefc && ModelSettings.Thmrlefc
-        RHS(1) = -C7(1) + (C2(1, 1) * SV.T(1) + C2(1, 2) * SV.T(2)) / Delt_t ...
-            - (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
-            + (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
+        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) ...
+            + (C1(indxBotm, 1) / Delt_t) * SV.h(indxBotm) + (C1(indxBotm, 2) / Delt_t) * SV.h(indxBotm + 1);
         for i = 2:ModelSettings.NL
             ARG4 = C1(i - 1, 2) / Delt_t;
             ARG5 = C1(i, 1) / Delt_t;
@@ -66,24 +77,24 @@
             - (C1(n - 1, 2) / Delt_t + C4(n - 1, 2)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
             + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
     else
-        RHS(1) = -C7(1) + (C2(1, 1) * SV.T(1) + C2(1, 2) * SV.T(2)) / Delt_t;
-        for i = 2:ModelSettings.NL
+        RHS(indxBotm) = -C7(indxBotm) + (C2(indxBotm, 1) * SV.T(indxBotm) + C2(1, 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
         RHS(n) = -C7(n) + (C2(n - 1, 2) * SV.T(n - 1) + C2(n, 1) * SV.T(n)) / Delt_t;
     end
 
-    for i = 1:ModelSettings.NN
+    for i = indxBotm:ModelSettings.NN
         for j = 1:ModelSettings.nD
             C5(i, j) = C2(i, j) / Delt_t + C5(i, j);
         end
     end
 
     EnergyMatrices.C5 = C5;
 
-    SAVE(1, 1, 2) = RHS(1);
-    SAVE(1, 2, 2) = C5(1, 1);
-    SAVE(1, 3, 2) = C5(1, 2);
+    SAVE(1, 1, 2) = RHS(indxBotm);
+    SAVE(1, 2, 2) = C5(indxBotm, 1);
+    SAVE(1, 3, 2) = C5(indxBotm, 2);
     SAVE(2, 1, 2) = RHS(n);
     SAVE(2, 2, 2) = C5(n - 1, 2);
     SAVE(2, 3, 2) = C5(n, 1);