Merge pull request #193 from EcoExtreML/fix_issue_98

Fix issue 98

CHANGELOG.md

@@ -4,10 +4,13 @@
 
 - All `cond*_` functions are refcatored and moved to `+conductivity` folder in
   [#189](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/189)
+- All `h_*` functions are refcatored and moved to `+soilmoisture` folder in
+  [#193](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/193)
 
 **Fixed:**
 
 - issue [#181](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/181)
+- issue [#98](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/98)
 
 
 <a name="1.3.0"></a>

src/+init/defineInitialValues.m

@@ -236,4 +236,7 @@
             InitialValues.(field{1}) = structure.(field{1});
         end
     end
+
+    % Arraies for calculating boundary flux;
+    InitialValues.SAVE = zeros(3, 3, 3);
 end

src/+io/define_constants.m

@@ -4,7 +4,9 @@
     const.h = 6.6262E-34; % [J s] Planck's constant
     const.c = 299792458; % [m s-1] Speed of light
     const.cp = 1004; % [J kg-1 K-1] Specific heat of dry air
+    const.cp_specific =  1.013E-3; % specific heat at cte pressure [MJ.kg-1.?C-1] FAO56 p26 box6
     const.R = 8.314; % [J mol-1K-1] Molar gas constant
+    const.R_specific = 0.287; % specific gas [kJ.kg-1.K-1]    FAO56 p26 box6
     const.rhoa = 1.2047; % [kg m-3] Specific mass of air
     const.kappa = 0.4; % [] Von Karman constant
     const.MH2O = 18; % [g mol-1] Molecular mass of water
@@ -37,4 +39,6 @@
     const.c_V = 1.870; % [J/g-1/Cels-1] Specific heat capacity of vapor
     const.c_a = 1.005; % [J/g-1/Cels-1] Specific heat capacity of dry air
     const.c_i = 2.0455; % [J/g-1/Cels-1] Specific heat capacity of ice
+    const.lambdav = 2.45; % latent heat of evaporation [MJ.kg-1] FAO56 pag 31
+    const.k = 0.41; % karman's cte   []  FAO 56 Eq4
 end

src/+io/getModelSettings.m

@@ -51,9 +51,6 @@
     ModelSettings.rroot = 1.5 * 1e-3;
     ModelSettings.SFCC = 1;
 
-    % Arraies for calculating boundary flux;
-    ModelSettings.SAVE = zeros(3, 3, 3);
-
     ModelSettings.Tot_Depth = 500; % Unit is cm. it should be usually bigger than 0.5m. Otherwise,
     ModelSettings.Eqlspace = 0; % Indicator for deciding is the space step equal or not; % the DeltZ would be reset in 50cm by hand;
 

src/+soilmoisture/assembleCoefficientMatrices.m

@@ -0,0 +1,72 @@
+function HeatMatrices = assembleCoefficientMatrices(HeatVariables, InitialValues, Srt)
+    %{
+        Assemble the coefficient matrices of Equation 4.32 (STEMMUS Technical
+        Notes, page 44).
+    %}
+    ModelSettings = io.getModelSettings();
+
+    % Define HeatMatrices structure
+    HeatMatrices.C1 = InitialValues.C1;
+    HeatMatrices.C2 = InitialValues.C2;
+    HeatMatrices.C3 = InitialValues.C3;
+    HeatMatrices.C4 = InitialValues.C4;
+    HeatMatrices.C5 = InitialValues.C5;
+    HeatMatrices.C6 = InitialValues.C6;
+    HeatMatrices.C7 = InitialValues.C7;
+    HeatMatrices.C9 = InitialValues.C9;
+    HeatMatrices.C4_a = [];
+    HeatMatrices.C5_a = [];
+
+    for i = 1:ModelSettings.NN
+        for j = 1:ModelSettings.nD
+            HeatMatrices.C1(i, j) = 0;
+            HeatMatrices.C7(i) = 0;
+            % C9 is the matrix coefficient of root water uptake
+            HeatMatrices.C9(i) = 0;
+            HeatMatrices.C4(i, j) = 0;
+            HeatMatrices.C4_a(i) = 0;
+            HeatMatrices.C5_a(i) = 0;
+            HeatMatrices.C2(i, j) = 0;
+            HeatMatrices.C3(i, j) = 0;
+            HeatMatrices.C5(i, j) = 0;
+            HeatMatrices.C6(i, j) = 0;
+        end
+    end
+    for i = 1:ModelSettings.NL
+        HeatMatrices.C1(i, 1) = HeatMatrices.C1(i, 1) + HeatVariables.Chh(i, 1) * ModelSettings.DeltZ(i) / 2;
+        HeatMatrices.C1(i + 1, 1) = HeatMatrices.C1(i + 1, 1) + HeatVariables.Chh(i, 2) * ModelSettings.DeltZ(i) / 2;
+
+        HeatMatrices.C2(i, 1) = HeatMatrices.C2(i, 1) + HeatVariables.ChT(i, 1) * ModelSettings.DeltZ(i) / 2;
+        HeatMatrices.C2(i + 1, 1) = HeatMatrices.C2(i + 1, 1) + HeatVariables.ChT(i, 2) * ModelSettings.DeltZ(i) / 2;
+
+        C4ARG1 = (HeatVariables.Khh(i, 1) + HeatVariables.Khh(i, 2)) / (2 * ModelSettings.DeltZ(i));
+        C4ARG2_1 = HeatVariables.Vvh(i, 1) / 3 + HeatVariables.Vvh(i, 2) / 6;
+        C4ARG2_2 = HeatVariables.Vvh(i, 1) / 6 + HeatVariables.Vvh(i, 2) / 3;
+        HeatMatrices.C4(i, 1) = HeatMatrices.C4(i, 1) + C4ARG1 - C4ARG2_1;
+        HeatMatrices.C4(i, 2) = HeatMatrices.C4(i, 2) - C4ARG1 - C4ARG2_2;
+        HeatMatrices.C4(i + 1, 1) = HeatMatrices.C4(i + 1, 1) + C4ARG1 + C4ARG2_2;
+        HeatMatrices.C4_a(i) = -C4ARG1 + C4ARG2_1;
+
+        C5ARG1 = (HeatVariables.KhT(i, 1) + HeatVariables.KhT(i, 2)) / (2 * ModelSettings.DeltZ(i));
+        C5ARG2_1 = HeatVariables.VvT(i, 1) / 3 + HeatVariables.VvT(i, 2) / 6;
+        C5ARG2_2 = HeatVariables.VvT(i, 1) / 6 + HeatVariables.VvT(i, 2) / 3;
+        HeatMatrices.C5(i, 1) = HeatMatrices.C5(i, 1) + C5ARG1 - C5ARG2_1;
+        HeatMatrices.C5(i, 2) = HeatMatrices.C5(i, 2) - C5ARG1 - C5ARG2_2;
+        HeatMatrices.C5(i + 1, 1) = HeatMatrices.C5(i + 1, 1) + C5ARG1 + C5ARG2_2;
+        HeatMatrices.C5_a(i) = -C5ARG1 + C5ARG2_1;
+
+        C6ARG = (HeatVariables.Kha(i, 1) + HeatVariables.Kha(i, 2)) / (2 * ModelSettings.DeltZ(i));
+        HeatMatrices.C6(i, 1) = HeatMatrices.C6(i, 1) + C6ARG;
+        HeatMatrices.C6(i, 2) = HeatMatrices.C6(i, 2) - C6ARG;
+        HeatMatrices.C6(i + 1, 1) = HeatMatrices.C6(i + 1, 1) + C6ARG;
+
+        C7ARG = (HeatVariables.Chg(i, 1) + HeatVariables.Chg(i, 2)) / 2;
+        HeatMatrices.C7(i) = HeatMatrices.C7(i) - C7ARG;
+        HeatMatrices.C7(i + 1) = HeatMatrices.C7(i + 1) + C7ARG;
+
+        C9ARG1 = (2 * Srt(i, 1) + Srt(i, 2)) * ModelSettings.DeltZ(i) / 6;
+        C9ARG2 = (Srt(i, 1) + 2 * Srt(i, 2)) * ModelSettings.DeltZ(i) / 6;
+        HeatMatrices.C9(i) = HeatMatrices.C9(i) + C9ARG1;
+        HeatMatrices.C9(i + 1) = HeatMatrices.C9(i + 1) + C9ARG2;
+    end
+end

src/+soilmoisture/calculateAerodynamicResistance.m

@@ -0,0 +1,11 @@
+function AerodynamicResistance = calculateAerodynamicResistance(U)
+    %{
+    %}
+    Constants = io.define_constants();
+    k = Constants.k;
+    hh_v = 0.12;
+    % FAO56 pag20 eq4- (d - zero displacement plane, z_0m - roughness length momentum transfer, z_0h - roughness length heat and vapour transfer, [m], FAO56 pag21 BOX4
+    AerodynamicResistance.vegetation = log((2 - (2 * hh_v / 3)) / (0.123 * hh_v)) * log((2 - (2 * hh_v / 3)) / (0.0123 * hh_v)) / ((k^2) * U) * 100;  % s m-1
+
+    AerodynamicResistance.soil  = log((2.0) / 0.001) * log(2.0 / 0.001) / ((k^2) * U) * 100;
+end

src/+soilmoisture/calculateBoundaryConditions.m

@@ -0,0 +1,79 @@
+function [AVAIL0, RHS, HeatMatrices, Precip] = calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ...
+                                                                           TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap)
+    %{
+        Determine the boundary condition for solving the soil moisture equation.
+    %}
+    ModelSettings = io.getModelSettings();
+    % n is the index of n_th item
+    n = ModelSettings.NN;
+
+    C4 = HeatMatrices.C4;
+    C4_a = HeatMatrices.C4_a;
+
+    Precip = InitialValues.Precip;
+    Precip_msr = ForcingData.Precip_msr;
+
+    Precipp = 0;
+    %  Apply the bottom boundary condition called for by BoundaryCondition.NBChB
+    if BoundaryCondition.NBChB == 1            %  Specify matric head at bottom to be ---BoundaryCondition.BChB;
+        RHS(1) = BoundaryCondition.BChB;
+        C4(1, 1) = 1;
+        RHS(2) = RHS(2) - C4(1, 2) * RHS(1);
+        C4(1, 2) = 0;
+        C4_a(1) = 0;
+    elseif BoundaryCondition.NBChB == 2        %  Specify flux at bottom to be ---BoundaryCondition.BChB (Positive upwards);
+        RHS(1) = RHS(1) + BoundaryCondition.BChB;
+    elseif BoundaryCondition.NBChB == 3        %  BoundaryCondition.NBChB=3,Gravity drainage at bottom--specify flux= hydraulic conductivity;
+        RHS(1) = RHS(1) - SoilVariables.KL_h(1, 1);
+    end
+
+    %  Apply the surface boundary condition called for by BoundaryCondition.NBCh
+    if BoundaryCondition.NBCh == 1             %  Specified matric head at surface---equal to hN;
+        % h_SUR: Observed matric potential at surface. This variable
+        % is not calculated anywhere! see issue 98, item 6
+        RHS(n) = InitialValues.h_SUR(KT);
+        C4(n, 1) = 1;
+        RHS(n - 1) = RHS(n - 1) - C4(n - 1, 2) * RHS(n);
+        C4(n - 1, 2) = 0;
+        C4_a(n - 1) = 0;
+    elseif BoundaryCondition.NBCh == 2
+        if BoundaryCondition.NBChh == 1
+            RHS(n) = hN;
+            C4(n, 1) = 1;
+            RHS(n - 1) = RHS(n - 1) - C4(n - 1, 2) * RHS(n);
+            C4(n - 1, 2) = 0;
+        else
+            RHS(n) = RHS(n) - BoundaryCondition.BCh;   % a specified matric head (saturation or dryness)was applied;
+        end
+    else
+        Precip_msr(KT) = min(Precip_msr(KT), SoilProperties.Ks0 / (3600 * 24) * TimeProperties.DELT * 10);
+        Precip_msr(KT) = min(Precip_msr(KT), SoilProperties.theta_s0 * 50 - ModelSettings.DeltZ(51:54) * SoilVariables.Theta_UU(51:54, 1) * 10);
+
+        if SoilVariables.Tss(KT) > 0
+            Precip(KT) = Precip_msr(KT) * 0.1 / TimeProperties.DELT;
+        else
+            Precip(KT) = Precip_msr(KT) * 0.1 / TimeProperties.DELT;
+            Precipp = Precipp + Precip(KT);
+            Precip(KT) = 0;
+        end
+
+        if SoilVariables.Tss(KT) > 0
+            AVAIL0 = Precip(KT) + Precipp + BoundaryCondition.DSTOR0 / Delt_t;
+            Precipp = 0;
+        else
+            AVAIL0 = Precip(KT) + BoundaryCondition.DSTOR0 / Delt_t;
+        end
+
+        if BoundaryCondition.NBChh == 1
+            RHS(n) = hN;
+            C4(n, 1) = 1;
+            RHS(n - 1) = RHS(n - 1) - C4(n - 1, 2) * RHS(n);
+            C4(n - 1, 2) = 0;
+            C4_a(n - 1) = 0;
+        else
+            RHS(n) = RHS(n) + AVAIL0 - Evap(KT);
+        end
+    end
+    HeatMatrices.C4 = C4;
+    HeatMatrices.C4_a = C4_a;
+end

src/+soilmoisture/calculateEvapotranspiration.m

@@ -0,0 +1,37 @@
+function [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt] = calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt)
+
+    ModelSettings = io.getModelSettings();
+
+    Rn = (ForcingData.Rn_msr(KT)) * 8.64 / 24 / 100 * 1;
+    Rn_SOIL = Rn * 0.68;
+
+    % calculate Aerodynamic Resistance
+    U = 100 .* (ForcingData.WS_msr(KT));
+    AR = soilmoisture.calculateAerodynamicResistance(U);
+    r_a_SOIL = AR.soil;
+
+    Ta = ForcingData.Ta_msr(KT);
+    Evap = InitialValues.Evap;
+    EVAP = InitialValues.EVAP;
+
+    if fluxes.lEctot < 1000 && fluxes.lEstot < 800 && fluxes.lEctot > -300 && fluxes.lEstot > -300 && any(SoilVariables.TT > 5)
+        lambda1      = (2.501 - 0.002361 * Ta) * 1E6;
+        lambda2      = (2.501 - 0.002361 * SoilVariables.Tss(KT)) * 1E6;
+        Evap(KT) = fluxes.lEstot / lambda2 * 0.1; % transfer to second value unit: cm s-1
+        EVAP(KT, 1) = Evap(KT);
+        Tp_t = fluxes.lEctot / lambda1 * 0.1; % transfer to second value
+        Srt1 = RWU * 100 ./ ModelSettings.DeltZ';
+    else
+        Evap(KT) = 0; % transfer to second value unit: cm s-1
+        EVAP(KT, 1) = Evap(KT);
+        Tp_t = 0; % transfer to second value
+        Srt1 = 0 ./ ModelSettings.DeltZ';
+    end
+
+    for i = 1:ModelSettings.NL
+        for j = 1:ModelSettings.nD
+            Srt(i, j) = Srt1(i, 1);
+        end
+    end
+    Trap(KT) = Tp_t;   % root water uptake integration by ModelSettings.DeltZ;
+end

src/+soilmoisture/calculateMatricCoefficients.m

@@ -0,0 +1,98 @@
+function [HeatVariables, SoilVariables] = calculateMatricCoefficients(SoilVariables, VaporVariables, GasDispersivity, InitialValues, ...
+                                                                      RHOV, DRHOVh, DRHOVT, D_Ta)
+    %{
+        Calculate all the parameters related to matric coefficients (e.g.,
+        c1-c8) as in Equation 4.32 (STEMMUS Technical Notes, page 44).
+    %}
+
+    ModelSettings = io.getModelSettings();
+    Constants = io.define_constants();
+
+    % Add a new variables to SoilVariables
+    SoilVariables.DTheta_LLT = [];
+    SoilVariables.SAVEDTheta_LLh = InitialValues.SAVEDTheta_LLh;
+    SoilVariables.SAVEDTheta_UUh = InitialValues.SAVEDTheta_UUh;
+
+    % Define new HeatVariables structure
+    HeatVariables.Chh = InitialValues.Chh;
+    HeatVariables.ChT = InitialValues.ChT;
+    HeatVariables.Khh = InitialValues.Khh;
+    HeatVariables.KhT = InitialValues.KhT;
+    HeatVariables.Kha = InitialValues.Kha;
+    HeatVariables.Vvh = InitialValues.Vvh;
+    HeatVariables.VvT = InitialValues.VvT;
+    HeatVariables.Chg = InitialValues.Chg;
+
+    % Make SV as an alias of SoilVariables to make codes shorter
+    SV = SoilVariables;
+
+    for i = 1:ModelSettings.NL
+        for j = 1:ModelSettings.nD
+            MN = i + j - 1;
+            if ModelSettings.hThmrl
+                CORhh = -1 * SV.CORh(MN);
+                DhUU = SV.COR(MN) * (SV.hh(MN) + SV.hh_frez(MN) - SV.h(MN) - SV.h_frez(MN) + (SV.hh(MN) + SV.hh_frez(MN)) * CORhh * (SV.TT(MN) - SV.T(MN)));
+                DhU = SV.COR(MN) * (SV.hh(MN) - SV.h(MN) + SV.hh(MN) * CORhh * (SV.TT(MN) - SV.T(MN)));
+
+                if DhU ~= 0 && DhUU ~= 0 && abs(SV.Theta_LL(i, j) - SV.Theta_L(i, j)) > 1e-6 && SV.DTheta_UUh(i, j) ~= 0
+                    SV.DTheta_UUh(i, j) = (SV.Theta_LL(i, j) - SV.Theta_L(i, j)) * SV.COR(MN) / DhUU;
+                end
+                if DhU ~= 0 && DhUU ~= 0 && abs(SV.Theta_LL(i, j) - SV.Theta_L(i, j)) > 1e-6 && SV.DTheta_LLh(i, j) ~= 0
+                    SV.DTheta_LLh(i, j) = (SV.Theta_UU(i, j) - SV.Theta_U(i, j)) * SV.COR(MN) / DhU;
+                end
+
+                SV.DTheta_LLT(i, j) = SV.DTheta_LLh(i, j) * (SV.hh(MN) * CORhh);
+
+                SV.SAVEDTheta_LLh(i, j) = SV.DTheta_LLh(i, j);
+                SV.SAVEDTheta_UUh(i, j) = SV.DTheta_UUh(i, j);
+            else
+                if abs(SV.TT(MN) - SV.T(MN)) > 1e-6
+                    SV.DTheta_LLT(i, j) = SV.DTheta_LLh(i, j) * (SV.hh(MN) / Constants.Gamma0) * (0.1425 + 4.76 * 10^(-4) * SV.TT(MN));
+                else
+                    SV.DTheta_LLT(i, j) = (SV.Theta_LL(i, j) - SV.Theta_L(i, j)) / (SV.TT(MN) - SV.T(MN));
+                end
+            end
+
+            HeatVariables.Chh(i, j) = (1 - RHOV(MN) / Constants.RHOL) * SV.DTheta_LLh(i, j) + SV.Theta_V(i, j) * DRHOVh(MN) / Constants.RHOL;
+            HeatVariables.Khh(i, j) = (VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) * DRHOVh(MN) / Constants.RHOL + SV.KfL_h(i, j);
+            HeatVariables.Chg(i, j) = SV.KfL_h(i, j);
+
+            if ModelSettings.Thmrlefc == 1
+                HeatVariables.ChT(i, j) = (1 - RHOV(MN) / Constants.RHOL) * SV.DTheta_LLT(i, j) + SV.Theta_V(i, j) * DRHOVT(MN) / Constants.RHOL;
+                HeatVariables.KhT(i, j) = (VaporVariables.D_V(i, j) * VaporVariables.Eta(i, j) + GasDispersivity.D_Vg(i)) * DRHOVT(MN) / Constants.RHOL + InitialValues.KL_T(i, j) + D_Ta(i, j);
+            end
+
+            if SV.KLa_Switch == 1
+                HeatVariables.Kha(i, j) = RHOV(MN) * GasDispersivity.Beta_g(i, j) / Constants.RHOL + SV.KfL_h(i, j) / Constants.Gamma_w;
+            else
+                HeatVariables.Kha(i, j) = 0;
+            end
+
+            if SV.DVa_Switch == 1
+                HeatVariables.Vvh(i, j) = -GasDispersivity.V_A(i) * DRHOVh(MN) / Constants.RHOL;
+                HeatVariables.VvT(i, j) = -GasDispersivity.V_A(i) * DRHOVT(MN) / Constants.RHOL;
+            else
+                HeatVariables.Vvh(i, j) = 0;
+                HeatVariables.VvT(i, j) = 0;
+            end
+            warningMsg = '%s is nan or not real';
+            if isnan(HeatVariables.Chh(i, j)) || ~isreal(HeatVariables.Chh(i, j))
+                warning(warningMsg, 'Chh(i, j)');
+            end
+            if isnan(HeatVariables.Khh(i, j)) || ~isreal(HeatVariables.Khh(i, j))
+                warning(warningMsg, 'Khh(i, j)');
+            end
+            if isnan(HeatVariables.Chg(i, j)) || ~isreal(HeatVariables.Chg(i, j))
+                warning(warningMsg, 'Chg(i, j)');
+            end
+            if isnan(HeatVariables.ChT(i, j)) || ~isreal(HeatVariables.ChT(i, j))
+                warning(warningMsg, 'ChT(i, j)');
+            end
+            if isnan(HeatVariables.KhT(i, j)) || ~isreal(HeatVariables.KhT(i, j))
+                warning(warningMsg, 'KhT(i, j)');
+            end
+        end
+    end
+    % Update SoilVariables using alias SV
+    SoilVariables = SV;
+end