diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1f25150b..d10dbb1a 100644
--- a/CHANGELOG.md
+++ b/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)
diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE
index af863234..0c297e70 100755
Binary files a/run_model_on_snellius/exe/STEMMUS_SCOPE and b/run_model_on_snellius/exe/STEMMUS_SCOPE differ
diff --git a/src/+init/defineInitialValues.m b/src/+init/defineInitialValues.m
index 85057a28..af769214 100644
--- a/src/+init/defineInitialValues.m
+++ b/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
diff --git a/src/+io/define_constants.m b/src/+io/define_constants.m
index a248b84f..6e308e96 100644
--- a/src/+io/define_constants.m
+++ b/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
diff --git a/src/+io/getModelSettings.m b/src/+io/getModelSettings.m
index 951458cc..b39cc6a6 100644
--- a/src/+io/getModelSettings.m
+++ b/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;
diff --git a/src/+soilmoisture/assembleCoefficientMatrices.m b/src/+soilmoisture/assembleCoefficientMatrices.m
new file mode 100644
index 00000000..077458c3
--- /dev/null
+++ b/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
diff --git a/src/+soilmoisture/calculateAerodynamicResistance.m b/src/+soilmoisture/calculateAerodynamicResistance.m
new file mode 100644
index 00000000..0f4ec3ce
--- /dev/null
+++ b/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
diff --git a/src/+soilmoisture/calculateBoundaryConditions.m b/src/+soilmoisture/calculateBoundaryConditions.m
new file mode 100644
index 00000000..a8f1513c
--- /dev/null
+++ b/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
diff --git a/src/+soilmoisture/calculateEvapotranspiration.m b/src/+soilmoisture/calculateEvapotranspiration.m
new file mode 100644
index 00000000..7855d745
--- /dev/null
+++ b/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
diff --git a/src/+soilmoisture/calculateMatricCoefficients.m b/src/+soilmoisture/calculateMatricCoefficients.m
new file mode 100644
index 00000000..eb8afc62
--- /dev/null
+++ b/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
diff --git a/src/+soilmoisture/calculateTimeDerivatives.m b/src/+soilmoisture/calculateTimeDerivatives.m
new file mode 100644
index 00000000..eac31918
--- /dev/null
+++ b/src/+soilmoisture/calculateTimeDerivatives.m
@@ -0,0 +1,94 @@
+function [RHS, HeatMatrices, boundaryFluxArray] = calculateTimeDerivatives(InitialValues, SoilVariables, HeatMatrices, Delt_t)
+ %{
+ Perform the finite difference of the time derivatives of the matrix
+ equation as Equation 4.32 shows in 'STEMMUS Technical Notes', section 4.
+ %}
+ % Extract variables that are nedded in this fucntion
+ C1 = HeatMatrices.C1;
+ C2 = HeatMatrices.C2;
+ C4 = HeatMatrices.C4;
+ C5 = HeatMatrices.C5;
+ C6 = HeatMatrices.C6;
+ C7 = HeatMatrices.C7;
+ C5_a = HeatMatrices.C5_a;
+ C9 = HeatMatrices.C9;
+ P_gg = InitialValues.P_gg;
+ h = SoilVariables.h;
+ T = SoilVariables.T;
+ TT = SoilVariables.TT;
+
+ % These are output
+ boundaryFluxArray = InitialValues.SAVE;
+ RHS = InitialValues.RHS;
+
+ ModelSettings = io.getModelSettings();
+
+ % n is the index of n_th item
+ n = ModelSettings.NN;
+ if ModelSettings.Thmrlefc && ~ModelSettings.Soilairefc
+ RHS(1) = -C9(1) - C7(1) + (C1(1, 1) * h(1) + C1(1, 2) * h(2)) / Delt_t ...
+ - (C2(1, 1) / Delt_t + C5(1, 1)) * TT(1) - (C2(1, 2) / Delt_t + C5(1, 2)) * TT(2) ...
+ + (C2(1, 1) / Delt_t) * T(1) + (C2(1, 2) / Delt_t) * T(2);
+ for i = 2:ModelSettings.NL
+ ARG1 = C2(i - 1, 2) / Delt_t;
+ ARG2 = C2(i, 1) / Delt_t;
+ ARG3 = C2(i, 2) / Delt_t;
+
+ RHS(i) = -C9(i) - C7(i) + (C1(i - 1, 2) * h(i - 1) + C1(i, 1) * h(i) + C1(i, 2) * h(i + 1)) / Delt_t ...
+ - (ARG1 + C5(i - 1, 2)) * TT(i - 1) - (ARG2 + C5(i, 1)) * TT(i) - (ARG3 + C5(i, 2)) * TT(i + 1) ...
+ + ARG1 * T(i - 1) + ARG2 * T(i) + ARG3 * T(i + 1);
+ end
+ RHS(n) = -C9(n) - C7(n) + (C1(n - 1, 2) * h(n - 1) + C1(n, 1) * h(n)) / Delt_t ...
+ - (C2(n - 1, 2) / Delt_t + C5(n - 1, 2)) * TT(n - 1) - (C2(n, 1) / Delt_t + C5(n, 1)) * TT(n) ...
+ + (C2(n - 1, 2) / Delt_t) * T(n - 1) + (C2(n, 1) / Delt_t) * T(n);
+ elseif ~ModelSettings.Thmrlefc && ModelSettings.Soilairefc
+ RHS(1) = -C9(1) - C7(1) + (C1(1, 1) * h(1) + C1(1, 2) * h(2)) / Delt_t ...
+ - C6(1, 1) * P_gg(1) - C6(1, 2) * P_gg(2);
+ for i = 2:ModelSettings.NL
+ RHS(i) = -C9(i) - C7(i) + (C1(i - 1, 2) * h(i - 1) + C1(i, 1) * h(i) + C1(i, 2) * h(i + 1)) / Delt_t ...
+ - C6(i - 1, 2) * P_gg(i - 1) - C6(i, 1) * P_gg(i) - C6(i, 2) * P_gg(i + 1);
+ end
+ RHS(n) = -C9(n) - C7(n) + (C1(n - 1, 2) * h(n - 1) + C1(n, 1) * h(n)) / Delt_t ...
+ - C6(n - 1, 2) * P_gg(n - 1) - C6(n, 1) * P_gg(n);
+ elseif ModelSettings.Thmrlefc && ModelSettings.Soilairefc
+ RHS(1) = -C9(1) - C7(1) + (C1(1, 1) * h(1) + C1(1, 2) * h(2)) / Delt_t ...
+ - (C2(1, 1) / Delt_t + C5(1, 1)) * TT(1) - (C2(1, 2) / Delt_t + C5(1, 2)) * TT(2) ...
+ - C6(1, 1) * P_gg(1) - C6(1, 2) * P_gg(2) ...
+ + (C2(1, 1) / Delt_t) * T(1) + (C2(1, 2) / Delt_t) * T(2);
+ for i = 2:ModelSettings.NL
+ ARG1 = C2(i - 1, 2) / Delt_t;
+ ARG2 = C2(i, 1) / Delt_t;
+ ARG3 = C2(i, 2) / Delt_t;
+
+ RHS(i) = -C9(i) - C7(i) + (C1(i - 1, 2) * h(i - 1) + C1(i, 1) * h(i) + C1(i, 2) * h(i + 1)) / Delt_t ...
+ - (ARG1 + C5_a(i - 1)) * TT(i - 1) - (ARG2 + C5(i, 1)) * TT(i) - (ARG3 + C5(i, 2)) * TT(i + 1) ...
+ - C6(i - 1, 2) * P_gg(i - 1) - C6(i, 1) * P_gg(i) - C6(i, 2) * P_gg(i + 1) ...
+ + ARG1 * T(i - 1) + ARG2 * T(i) + ARG3 * T(i + 1);
+ end
+ RHS(n) = -C9(n) - C7(n) + (C1(n - 1, 2) * h(n - 1) + C1(n, 1) * h(n)) / Delt_t ...
+ - (C2(n - 1, 2) / Delt_t + C5_a(n - 1)) * TT(n - 1) - (C2(n, 1) / Delt_t + C5(n, 1)) * TT(n) ...
+ - C6(n - 1, 2) * P_gg(n - 1) - C6(n, 1) * P_gg(n) ...
+ + (C2(n - 1, 2) / Delt_t) * T(n - 1) + (C2(n, 1) / Delt_t) * T(n);
+ else
+ RHS(1) = -C9(1) - C7(1) + (C1(1, 1) * h(1) + C1(1, 2) * h(2)) / Delt_t;
+ for i = 2:ModelSettings.NL
+ RHS(i) = -C9(i) - C7(i) + (C1(i - 1, 2) * h(i - 1) + C1(i, 1) * h(i) + C1(i, 2) * h(i + 1)) / Delt_t;
+ end
+ RHS(n) = -C9(n) - C7(n) + (C1(n - 1, 2) * h(n - 1) + C1(n, 1) * h(n)) / Delt_t;
+ end
+
+ for i = 1:ModelSettings.NN
+ for j = 1:ModelSettings.nD
+ C4(i, j) = C1(i, j) / Delt_t + C4(i, j);
+ end
+ end
+
+ boundaryFluxArray(1, 1, 1) = RHS(1);
+ boundaryFluxArray(1, 2, 1) = C4(1, 1);
+ boundaryFluxArray(1, 3, 1) = C4(1, 2);
+ boundaryFluxArray(2, 1, 1) = RHS(n);
+ boundaryFluxArray(2, 2, 1) = C4(n - 1, 2);
+ boundaryFluxArray(2, 3, 1) = C4(n, 1);
+
+ HeatMatrices.C4 = C4;
+end
diff --git a/src/+soilmoisture/calculatesSoilWaterFluxes.m b/src/+soilmoisture/calculatesSoilWaterFluxes.m
new file mode 100644
index 00000000..d6ca4310
--- /dev/null
+++ b/src/+soilmoisture/calculatesSoilWaterFluxes.m
@@ -0,0 +1,9 @@
+function HBoundaryFlux = calculatesSoilWaterFluxes(SAVE, hh)
+ %{
+ Calculate the soil water fluxes on the boundary node.
+ %}
+ ModelSettings = io.getModelSettings();
+
+ HBoundaryFlux.QMT = SAVE(2, 1, 1) - SAVE(2, 2, 1) * hh(ModelSettings.NN - 1) - SAVE(2, 3, 1) * hh(ModelSettings.NN);
+ HBoundaryFlux.QMB = -SAVE(1, 1, 1) + SAVE(1, 2, 1) * hh(1) + SAVE(1, 3, 1) * hh(2);
+end
diff --git a/src/+soilmoisture/solveSoilMoistureBalance.m b/src/+soilmoisture/solveSoilMoistureBalance.m
new file mode 100644
index 00000000..06253c3b
--- /dev/null
+++ b/src/+soilmoisture/solveSoilMoistureBalance.m
@@ -0,0 +1,50 @@
+function [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip] = solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, ...
+ BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt)
+ %{
+ Solve the soil moisture balance equation with the Thomas algorithm to
+ update the soil matric potential 'hh', the finite difference
+ time-stepping scheme of this soil moisture equation is derived in
+ 'STEMMUS Technical Notes' section 4, Equation 4.32.
+ %}
+ [HeatVariables, SoilVariables] = soilmoisture.calculateMatricCoefficients(SoilVariables, VaporVariables, GasDispersivity, InitialValues, ...
+ RHOV, DRHOVh, DRHOVT, D_Ta);
+
+ HeatMatrices = soilmoisture.assembleCoefficientMatrices(HeatVariables, InitialValues, Srt);
+
+ [RHS, HeatMatrices, boundaryFluxArray] = soilmoisture.calculateTimeDerivatives(InitialValues, SoilVariables, HeatMatrices, Delt_t);
+
+ % When BoundaryCondition.NBCh == 3, otherwise Rn_SOIL, Evap, EVAP, Trap,
+ % r_a_SOIL, Srt will be empty arrays! see issue 98, item 2
+ if BoundaryCondition.NBCh == 3
+ [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt] = soilmoisture.calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt);
+ else
+ Rn_SOIL = InitialValues.Rn_SOIL;
+ Evap = InitialValues.Evap;
+ EVAP = InitialValues.EVAP;
+ Trap = [];
+ r_a_SOIL = InitialValues.r_a_SOIL;
+ end
+
+ [AVAIL0, RHS, HeatMatrices, Precip] = soilmoisture.calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ...
+ TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap);
+
+ [CHK, hh, C4] = soilmoisture.solveTridiagonalMatrixEquations(HeatMatrices.C4, SoilVariables.hh, HeatMatrices.C4_a, RHS);
+ % update structures
+ SoilVariables.hh = hh;
+ HeatMatrices.C4 = C4;
+
+ % fix hh values
+ ModelSettings = io.getModelSettings();
+ for i = 1:ModelSettings.NN
+ if isnan(SoilVariables.hh(i)) || SoilVariables.hh(i) <= -1E12
+ SoilVariables.hh(i) = hOLD(i);
+ end
+ if SoilVariables.hh(i) > -1e-6
+ SoilVariables.hh(i) = -1e-6;
+ end
+ end
+
+ % calculate boundary flux
+ HBoundaryFlux = soilmoisture.calculatesSoilWaterFluxes(boundaryFluxArray, SoilVariables.hh);
+
+end
diff --git a/src/+soilmoisture/solveTridiagonalMatrixEquations.m b/src/+soilmoisture/solveTridiagonalMatrixEquations.m
new file mode 100644
index 00000000..4a86287f
--- /dev/null
+++ b/src/+soilmoisture/solveTridiagonalMatrixEquations.m
@@ -0,0 +1,27 @@
+function [CHK, hh, C4] = solveTridiagonalMatrixEquations(C4, hh, C4_a, RHS)
+ %{
+ Solve the tridiagonal matrix equations using Thomas algorithm, which is
+ in the form of Equation 4.25 (STEMMUS Technical Notes, page 41).
+ %}
+ ModelSettings = io.getModelSettings();
+
+ RHS(1) = RHS(1) / C4(1, 1);
+
+ for i = 2:ModelSettings.NN
+ C4(i, 1) = C4(i, 1) - C4_a(i - 1) * C4(i - 1, 2) / C4(i - 1, 1);
+ RHS(i) = (RHS(i) - C4_a(i - 1) * RHS(i - 1)) / C4(i, 1);
+ end
+
+ for i = ModelSettings.NL:-1:1
+ RHS(i) = RHS(i) - C4(i, 2) * RHS(i + 1) / C4(i, 1);
+ end
+ for i = 1:ModelSettings.NN
+ CHK(i) = abs(RHS(i) - hh(i));
+ hh(i) = RHS(i);
+ SAVEhh(i) = hh(i);
+ end
+
+ if isnan(SAVEhh) == 1 || ~isreal(SAVEhh)
+ warning('\n SAVEhh == 1 or not real\r');
+ end
+end
diff --git a/src/Enrgy_BC.m b/src/Enrgy_BC.m
index 1652fa88..6accefb5 100644
--- a/src/Enrgy_BC.m
+++ b/src/Enrgy_BC.m
@@ -12,7 +12,7 @@
elseif NBCTB == 2
RHS(1) = RHS(1) + BCTB;
else
- C5(1, 1) = C5(1, 1) - c_L * RHOL * QMB(KT);
+ C5(1, 1) = C5(1, 1) - c_L * RHOL * QMB;
end
%%%%%%%%%% Apply the surface boundary condition called by NBCT %%%%%%%%%%%%
diff --git a/src/Evap_Cal.m b/src/Evap_Cal.m
deleted file mode 100644
index 79acf577..00000000
--- a/src/Evap_Cal.m
+++ /dev/null
@@ -1,332 +0,0 @@
-function [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Resis_a, Srt] = Evap_Cal(DeltZ, TIME, RHOV, Ta, HR_a, U, Theta_LL, Ts, Rv, g, NL, NN, KT, hh, rwuef, Theta_UU, Rn, T, TT, Gvc, Rns, Srt)
- global LAI rl_min RWU lEstot lEctot Tss
-
- %%%%%%% LAI and light extinction coefficient calculation %%%%%%%%%%%%%%%%%%
- %%%%%%% LAI and light extinction coefficient calculation %%%%%%%%%%%%%%%%%%
- Tao = 0.56; % light attenual coefficient
-
- LAI(KT) = Gvc(KT); % 0.9;
-
- if LAI(KT) <= 0.01
- LAI(KT) = 0.01;
- end
- % LAI(KT)=1.2;
- if LAI(KT) <= 2
- LAI_act(KT) = LAI(KT);
- elseif LAI(KT) <= 4
- LAI_act(KT) = 2;
- else
- LAI_act(KT) = 0.5 * LAI(KT);
- end
- % LAI(KT)=0;
- Gvc(KT) = 1 - exp(-1 * (Tao * LAI(KT)));
-
- % Set constants
- sigma = 4.903e-9; % Stefan Boltzmann constant MJ.m-2.day-1 FAO56 pag 74
- lambdav = 2.45; % latent heat of evaporation [MJ.kg-1] FAO56 pag 31
- % Gieske 2003 pag 74 Eq33/DKTgman 2002
- % lambda=2.501-2.361E-3*t, with t temperature evaporative surface (?C)
- % see script Lambda_function_t.py
- Gsc = 0.082; % solar constant [MJ.m-2.mKT-1] FAO56 pag 47 Eq28
- eps = 0.622; % ratio molecular weigth of vapour/dry air FAO56 p26 BOX6
- R = 0.287; % specific gas [kJ.kg-1.K-1] FAO56 p26 box6
- Cp = 1.013E-3; % specific heat at cte pressure [MJ.kg-1.?C-1] FAO56 p26 box6
- k = 0.41; % karman's cte [] FAO 56 Eq4
- Z = 3421; % altitute of the location(m)
- as = 0.25; % regression constant, expressKTg the fraction of extraterrestrial radiation FAO56 pag50
- bs = 0.5;
- alfa = 0.23; % albeo of vegetation set as 0.23
- z_m = 10; % observation height of wKTd speed; 10m
- Lz = 240 * pi() / 180; % latitude of Beijing time zone west of Greenwich
- Lm = 252 * pi() / 180; % latitude of Local time, west of Greenwich
- % albedo of soil calculation;
- Theta_LL_sur(KT) = Theta_LL(NL, 2);
- Theta_UU_sur(KT) = Theta_UU(NL, 2);
-
- if Theta_LL_sur(KT) < 0.1
- alfa_s(KT) = 0.25;
- elseif Theta_LL_sur(KT) < 0.25
- alfa_s(KT) = 0.35 - Theta_LL_sur(KT);
- else
- alfa_s(KT) = 0.1;
- end
- %% AIR PARAMETERS CALCULATION
- % compute DELTA - SLOPE OF SATURATION VAPOUR PRESSURE CURVE
- % [kPa.?C-1]
- % FAO56 pag 37 Eq13
- DELTA(KT) = 4098 * (0.6108 * exp(17.27 * Ta(KT) / (Ta(KT) + 237.3))) / (Ta(KT) + 237.3)^2;
- % ro_a - MEAN AIR DENSITY AT CTE PRESSURE
- % [kg.m-3]
- % FAO56 pag26 box6
- Pa = 101.3 * ((293 - 0.0065 * Z) / 293)^5.26;
- ro_a(KT) = Pa / (R * 1.01 * (Ta(KT) + 273.16));
- % compute e0_Ta - saturation vapour pressure at actual air temperature
- % [kPa]
- % FAO56 pag36 Eq11
- e0_Ta(KT) = 0.6108 * exp(17.27 * Ta(KT) / (Ta(KT) + 237.3));
- e0_Ts(KT) = 0.6108 * exp(17.27 * Ts(KT) / (Ts(KT) + 237.3));
- % compute e_a - ACTUAL VAPOUR PRESSURE
- % [kPa]
- % FAO56 pag74 Eq54
- e_a(KT) = e0_Ta(KT) * HR_a(KT);
- e_a_Ts(KT) = e0_Ts(KT) * HR_a(KT);
-
- % gama - PSYCHROMETRIC CONSTANT
- % [kPa.?C-1]
- % FAO56 pag31 eq8
- gama = 0.664742 * 1e-3 * Pa;
-
- hh_v(KT) = 0.12;
- rl_min(KT) = 100;
- % Rn(KT) =Rn_msr(KT);
- if TIME <= 1800 * 3600
- Rn_SOIL(KT) = Rn(KT) * 0.68; % exp(-1*(Tao*LAI(KT)))net radiation for soil
- else
- Rn_SOIL(KT) = Rn(KT) * 0.68;
- end
- %% SURFACE RESISTANCE PARAMETERS CALCULATION
- R_a = 0.81;
- R_b = 0.004 * 24 * 11.6;
- R_c = 0.05;
- rl(KT) = rl_min(KT) / ((R_b * Rn(KT) + R_c) / (R_a * (R_b * Rn(KT) + 1)));
-
- % r_s - SURFACE RESISTANCE
- % [s.m-1]
- % VEG: Dingman pag 208 (canopy conductance) (equivalent to FAO56 pag21 Eq5)
- r_s_VEG(KT) = rl(KT) / LAI_act(KT);
-
- % SOIL: equation 20 of van de Griend and Owe, 1994
-
- if KT <= 1047
- r_s_SOIL(KT) = 10.0 * exp(0.3563 * 100.0 * (0.2050 - Theta_LL_sur(KT))); % 0.25 set as minmum soil moisture for potential evaporation
- elseif KT <= 4300
- r_s_SOIL(KT) = 10.0 * exp(0.3563 * 100.0 * (0.205 - Theta_LL_sur(KT))); % 0.25 set as minmum soil moisture for potential evaporation
- else
- r_s_SOIL(KT) = 10.0 * exp(0.3563 * 100.0 * (0.205 - Theta_LL_sur(KT))); % 0.25 set as minmum soil moisture for potential evaporation
- end
- % r_a - AERODYNAMIC RESISTANCE
- % [s.m-1]
- % 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
- r_a_VEG(KT) = log((2 - (2 * hh_v(KT) / 3)) / (0.123 * hh_v(KT))) * log((2 - (2 * hh_v(KT) / 3)) / (0.0123 * hh_v(KT))) / ((k^2) * U(KT)) * 100; % s m-1
- % r_a of SOIL
- % Liu www.hydrol-earth-syst-sci.net/11/769/2007/
- % equation for neutral conditions (eq. 9)
- % only function of ws, it is assumed that roughness are the same for any type of soil
-
- RHOV_sur(KT) = RHOV(NN);
- % Theta_LL_sur(KT)=Theta_LL(NL,2);
- %
- % Theta_UU_sur(KT)=Theta_UU(NL,2);
- P_Va(KT) = 0.611 * exp(17.27 * Ta(KT) / (Ta(KT) + 237.3)) * HR_a(KT); % The atmospheric vapor pressure (KPa) (1000Pa=1000.1000.g.100^-1.cm^-1.s^-2)
-
- RHOV_A(KT) = P_Va(KT) * 1e4 / (Rv * (Ta(KT) + 273.15)); % g.cm^-3; Rv-cm^2.s^-2.Cels^-1
-
- z_ref = 200; % cm The reference height of tempearature measurement (usually 2 m)
- d0_disp = 0; % cm The zero-plane displacement (=0 m)
- z_srT = 0.1; % cm The surface roughness for the heat flux (=0.001m) 0.01m
- VK_Const = 0.41; % The von Karman constant (=0.41)
- z_srm = 0.1; % cm The surface roughness for momentum flux (=0.001m) 0.01m
- U_wind = 198.4597; % The mean wind speed at reference height.(cm.s^-1)
-
- MO(KT) = (Ta(KT) * U(KT)^2) / (g * (Ta(KT) - T(NN)) * log(z_ref / z_srm)); % Wind speed should be in cm.s^-1, MO-cm;
-
- Zeta_MO(KT) = z_ref / MO(KT);
-
- if abs(Ta(KT) - T(NN)) <= 0.01
- Stab_m(KT) = 0;
- Stab_T(KT) = 0;
- elseif Zeta_MO(KT) < 0 % Ta(KT) 1
- Stab_T(KT) = 5;
- Stab_m(KT) = 5;
- else
- Stab_T(KT) = 5 * Zeta_MO(KT);
- Stab_m(KT) = 5 * Zeta_MO(KT);
- end
- end
-
- Velo_fric(KT) = U(KT) * VK_Const / (log((z_ref - d0_disp + z_srm) / z_srm) + Stab_m(KT));
-
- Resis_a(KT) = (log((z_ref - d0_disp + z_srT) / z_srT) + Stab_T(KT)) / (VK_Const * Velo_fric(KT)); % (s.cm^-1)
-
- Resis_s(KT) = 10 * exp(35.63 * (0.205 - Theta_LL_sur(KT))) / 100; % (-805+4140*(Theta_s(J)-Theta_LL_sur(KT)))/100; % s.m^-1----->0.001s.cm^-1
-
- r_a_SOIL(KT) = log((2.0) / 0.001) * log(2.0 / 0.001) / ((k^2) * U(KT)) * 100; % (s.m^-1) ORIGINAL Zom=0.0058
- Evapo(KT) = (RHOV_sur(KT) - RHOV_A(KT)) / (Resis_s(KT) + r_a_SOIL(KT) / 100);
-
- % PT/PE - Penman-Montheith
- % mm.day-1
- % FAO56 pag19 eq3
- % VEG
- PT_PM_VEG(KT) = (DELTA(KT) * (Rn(KT)) + 3600 * ro_a(KT) * Cp * (e0_Ta(KT) - e_a(KT)) / r_a_VEG(KT)) / (lambdav * (DELTA(KT) + gama * (1 + r_s_VEG(KT) / r_a_VEG(KT)))) / 3600;
-
- if LAI(KT) == 0 || hh_v(KT) == 0
- PT_PM_VEG(KT) = 0;
- end
- PE_PM_SOIL(KT) = (DELTA(KT) * (Rn_SOIL(KT)) + 3600 * ro_a(KT) * Cp * (e0_Ta(KT) - e_a(KT)) / r_a_SOIL(KT)) / (lambdav * (DELTA(KT) + gama * (1 + r_s_SOIL(KT) / r_a_SOIL(KT)))) / 3600;
- Evap(KT) = 0.1 * PE_PM_SOIL(KT); % transfer to second value-G_SOIL(KT)
- % Evap(KT)=LE_msr(KT);
- EVAP(KT, 1) = Evap(KT);
- Tp_t(KT) = 0.1 * PT_PM_VEG(KT); % transfer to second value
- TP_t(KT, 1) = Tp_t(KT);
-
- if rwuef == 1
- if KT <= 3288 + 1103
- H1 = 0;
- H2 = -31;
- H4 = -8000;
- H3L = -600;
- H3H = -300;
- if Tp_t(KT) < 0.02 / 3600
- H3 = H3L;
- elseif Tp_t(KT) > 0.05 / 3600
- H3 = H3H;
- else
- H3 = H3H + (H3L - H3H) / (0.03 / 3600) * (0.05 / 3600 - Tp_t(KT));
- end
- else
- H1 = -1;
- H2 = -5;
- H4 = -16000;
- H3L = -600;
- H3H = -300;
- if Tp_t(KT) < 0.02 / 3600
- H3 = H3L;
- elseif Tp_t(KT) > 0.05 / 3600
- H3 = H3H;
- else
- H3 = H3H + (H3L - H3H) / (0.03 / 3600) * (0.05 / 3600 - Tp_t(KT));
- end
- end
- % piecewise linear reduction function
- MN = 0;
- for ML = 1:NL
- for ND = 1:2
- MN = ML + ND - 1;
- if hh(MN) >= H1
- alpha_h(ML, ND) = 0;
- elseif hh(MN) >= H2
- alpha_h(ML, ND) = (H1 - hh(MN)) / (H1 - H2);
- elseif hh(MN) >= H3
- alpha_h(ML, ND) = 1;
- elseif hh(MN) >= H4
- alpha_h(ML, ND) = (hh(MN) - H4) / (H3 - H4);
- else
- alpha_h(ML, ND) = 0;
- end
- end
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%% Root lenth distribution %%%%%%%%%%%%%%%%%%%%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- Elmn_Lnth = 0;
- RL = sum(DeltZ);
- Elmn_Lnth(1) = 0;
- RB = 0.9;
- LR = log(0.01) / log(RB);
- RY = 1 - RB^(LR);
-
- if LR <= 1
- for ML = 1:NL - 1 % ignore the surface root water uptake 1cm
- for ND = 1:2
- MN = ML + ND - 1;
- bx(ML, ND) = 0;
- end
- end
- else
- FRY = zeros(NL, 1);
- bx = zeros(NL, 2);
- for ML = 2:NL
- Elmn_Lnth(ML) = Elmn_Lnth(ML - 1) + DeltZ(ML - 1);
- if Elmn_Lnth < RL - LR % (KT)
- % bx(ML)=0;
- FRY(ML) = 1;
- % elseif Elmn_Lnth>=RL-LR(KT) && Elmn_Lnth -300 && lEstot > -300 && any(TT > 5)
- lambda1 = (2.501 - 0.002361 * Ta(KT)) * 1E6;
- lambda2 = (2.501 - 0.002361 * Tss) * 1E6;
- Evap(KT) = lEstot / lambda2 * 0.1; % transfer to second value unit: cm s-1
- EVAP(KT, 1) = Evap(KT);
- Tp_t(KT) = lEctot / lambda1 * 0.1; % transfer to second value
- Srt1 = RWU * 100 ./ DeltZ';
- else
- Evap(KT) = 0; % transfer to second value unit: cm s-1
- EVAP(KT, 1) = Evap(KT);
- Tp_t(KT) = 0; % transfer to second value
- Srt1 = 0 ./ DeltZ';
- end
- for ML = 1:NL
- for ND = 1:2
- MN = ML + ND - 1;
- Srt(ML, ND) = Srt1(ML, 1);
- end
- end
- Trap(KT) = 0;
- Trap(KT) = Tp_t(KT); % root water uptake integration by DeltZ;
-
-end
diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m
index 6c6709f7..3ed21a07 100644
--- a/src/STEMMUS_SCOPE.m
+++ b/src/STEMMUS_SCOPE.m
@@ -32,39 +32,33 @@
[InputPath, OutputPath, InitialConditionPath] = io.read_config(CFG);
% Prepare forcing and soil data
-global DELT SaturatedMC ResidualMC fieldMC theta_s0 Ks0
+global SaturatedMC ResidualMC fieldMC theta_s0
[SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath);
landcoverClass = SiteProperties.landcoverClass;
SaturatedMC = SoilProperties.SaturatedMC; % used in calc_rssrbs
ResidualMC = SoilProperties.ResidualMC; % used in calc_rssrbs
fieldMC = SoilProperties.fieldMC; % used in calc_rssrbs
-theta_s0 = SoilProperties.theta_s0; % used in h_BC
-Ks0 = SoilProperties.Ks0; % used in h_BC
-DELT = TimeProperties.DELT; % used in h_BC
+theta_s0 = SoilProperties.theta_s0; % used in select_input
% Load model settings: replacing "run Constants"
ModelSettings = io.getModelSettings();
-global J rwuef SWCC Thmrlefc Soilairefc hThmrl KT TIME Delt_t NN ML nD
-global fc T0 rroot SAVE NL DeltZ
+global J Thmrlefc Soilairefc KT TIME Delt_t NN ML rwuef
+global T0 rroot SAVE NL DeltZ
NL = ModelSettings.NL;
DeltZ = ModelSettings.DeltZ;
DeltZ_R = ModelSettings.DeltZ_R;
J = ModelSettings.J;
-SWCC = ModelSettings.SWCC;
Thmrlefc = ModelSettings.Thmrlefc;
Soilairefc = ModelSettings.Soilairefc;
-hThmrl = ModelSettings.hThmrl;
-fc = ModelSettings.fc;
T0 = ModelSettings.T0;
-rwuef = ModelSettings.rwuef;
rroot = ModelSettings.rroot;
-SAVE = ModelSettings.SAVE;
NIT = ModelSettings.NIT;
KT = ModelSettings.KT;
NN = ModelSettings.NN;
ML = ModelSettings.ML;
-nD = ModelSettings.nD;
+rwuef = ModelSettings.rwuef;
+
% defined as global, and used in other scripts
TIME = 0; % Time of simulation released;
Delt_t = TimeProperties.DELT; % Duration of time step [Unit of second]
@@ -72,43 +66,32 @@
% load forcing data
ForcingData = io.loadForcingData(InputPath, TimeProperties, SoilProperties.fmax, ModelSettings.Tot_Depth);
-global Tmin LAI_msr G_msr Precip_msr
+global Tmin LAI_msr
LAI_msr = ForcingData.LAI_msr; % used in Root_properties
-Precip_msr = ForcingData.Precip_msr; % used in h_BC and h_sub
Tmin = ForcingData.Tmin; % used in Enrgy_sub
-global MN ND hOLD TOLD h hh T TT P_g P_gg Evap QMT hN Trap RWU EVAP theta_s0 Ks0
-global Precip frac SUMTIME TTT Theta_LLL CHK Theta_LL Theta_L Theta_UUU Theta_UU
-global Theta_U Theta_III Theta_II AVAIL0 SRT alpha_h bx Srt CTT_PH
-global CTT_LT CTT_g CTT_Lg c_unsat DhDZ DTDZ DRHOVZ QL QL_h QL_T QV Qa KL_h Chh ChT
-global Khh KhT Resis_a KfL_h KfL_T TT_CRIT h_frez L_f CTT EPCT DTheta_LLh DTheta_LLT
+global MN ND TOLD h hh T TT P_g P_gg RWU EVAP QMB
+global Precip frac TTT Theta_LLL CHK Theta_LL Theta_UUU
+global Theta_III Theta_II SRT Srt CTT_PH
+global CTT_LT CTT_g CTT_Lg c_unsat DhDZ DTDZ DRHOVZ QL QL_h QL_T QV Qa KL_h
+global Khh KhT Resis_a KfL_h TT_CRIT h_frez L_f CTT EPCT DTheta_LLh DTheta_LLT
global DTheta_UUh Lambda_eff DDhDZ DEhBAR DRHOVhDz EtaBAR D_Vg
global DRHOVTDz KLhBAR KLTBAR DTDBAR SAVEDTheta_LLh SAVEDTheta_UUh QVT QVH HR QVa
global QLH QLT DVH DVT Se QL_a DPgDZ V_A Theta_V W WW D_Ta thermal Xaa
global XaT Xah KL_T DRHOVT DRHOVh DRHODAt DRHODAz Theta_g Beta_g D_V Eta
-global Ks RHODA RHOV L Evapo Gvc
-global sfactortot sfactor fluxes lEstot lEctot Tss
+global Ks RHODA RHOV L
+global sfactor fluxes Tss
% Get initial values
InitialValues = init.defineInitialValues(TimeProperties.Dur_tot);
-alpha_h = InitialValues.alpha_h;
-bx = InitialValues.bx;
-Srt = InitialValues.Srt;
-SAVEDTheta_UUh = InitialValues.SAVEDTheta_UUh;
-SAVEDTheta_LLh = InitialValues.SAVEDTheta_LLh;
-Lambda_eff = InitialValues.Lambda_eff;
D_V = InitialValues.D_V;
Eta = InitialValues.Eta;
-Chh = InitialValues.Chh;
-ChT = InitialValues.ChT;
Khh = InitialValues.Khh;
KhT = InitialValues.KhT;
QL = InitialValues.QL;
QL_h = InitialValues.QL_h;
QL_T = InitialValues.QL_T;
V_A = InitialValues.V_A;
-Beta_g = InitialValues.Beta_g;
-c_unsat = InitialValues.c_unsat;
CTT_PH = InitialValues.CTT_PH;
CTT_Lg = InitialValues.CTT_Lg;
CTT_g = InitialValues.CTT_g;
@@ -137,17 +120,12 @@
DPgDZ = InitialValues.DPgDZ;
QL_a = InitialValues.QL_a;
frac = InitialValues.frac;
-Precip = InitialValues.Precip;
-h_SUR = InitialValues.h_SUR;
-Evap = InitialValues.Evap;
-sfactortot = InitialValues.sfactortot;
-EVAP = InitialValues.EVAP;
P_g = InitialValues.P_g;
P_gg = InitialValues.P_gg;
T_CRIT = InitialValues.T_CRIT;
TT_CRIT = InitialValues.TT_CRIT;
EPCT = InitialValues.EPCT;
-HR = InitialValues.HR;
+HR = InitialValues.HR; % used in Density_V
RHOV = InitialValues.RHOV;
DRHOVh = InitialValues.DRHOVh;
DRHOVT = InitialValues.DRHOVT;
@@ -160,20 +138,11 @@
L = InitialValues.L;
hOLD = InitialValues.hOLD;
TOLD = InitialValues.TOLD;
+SAVE = InitialValues.SAVE;
-global Kha Vvh VvT Chg C1 C2 C3 C4 C5 C6 Cah CaT Caa Kah KaT Kaa Vah VaT Vaa Cag CTh CTa KTh KTT KTa
-global VTT VTh VTa CTg Kcva Kcah KcaT Kcaa Ccah CcaT Ccaa SMC bbx Ta Ts U HR_a Rns Rn
-global RHOV_s DRHOV_sT Tbtm r_a_SOIL Rn_SOIL SH MO Zeta_MO TopPg Tp_t RHS C7 C9
-Kha = InitialValues.Kha;
-Vvh = InitialValues.Vvh;
-VvT = InitialValues.VvT;
-Chg = InitialValues.Chg;
-C1 = InitialValues.C1;
-C2 = InitialValues.C2;
-C3 = InitialValues.C3;
-C4 = InitialValues.C4;
-C5 = InitialValues.C5;
-C6 = InitialValues.C6;
+global Kha Vvh VvT C1 C2 C3 C4 C5 C5_a C6 Cah CaT Caa Kah KaT Kaa Vah VaT Vaa Cag CTh CTa KTh KTT KTa
+global VTT VTh VTa CTg Kcva Kcah KcaT Kcaa Ccah CcaT Ccaa SMC bbx Ta Ts
+global RHOV_s DRHOV_sT r_a_SOIL Rn_SOIL SH TopPg RHS C7
Cah = InitialValues.Cah;
CaT = InitialValues.CaT;
Caa = InitialValues.Caa;
@@ -204,28 +173,14 @@
bbx = InitialValues.bbx;
Ta = InitialValues.Ta;
Ts = InitialValues.Ts;
-U = InitialValues.U;
-HR_a = InitialValues.HR_a;
-Rns = InitialValues.Rns;
-Rn = InitialValues.Rn;
-SH = InitialValues.SH;
-MO = InitialValues.MO;
-Zeta_MO = InitialValues.Zeta_MO;
-TopPg = InitialValues.TopPg;
-Tp_t = InitialValues.Tp_t;
RHS = InitialValues.RHS;
-C7 = InitialValues.C7;
-C9 = InitialValues.C9;
RHOV_s = InitialValues.RHOV_s;
DRHOV_sT = InitialValues.DRHOV_sT;
P_gOLD = InitialValues.P_gOLD;
-Tbtm = InitialValues.Tbtm;
-r_a_SOIL = InitialValues.r_a_SOIL;
-Rn_SOIL = InitialValues.Rn_SOIL;
%% 1. define Constants
Constants = io.define_constants();
-global g RHOL RHOI Rv RDA c_a c_V c_L Hc c_i Gamma0 Gamma_w Rl
+global g RHOL RHOI Rv RDA c_a c_V c_L Hc c_i Gamma_w Rl
g = Constants.g;
RHOL = Constants.RHOL;
RHOI = Constants.RHOI;
@@ -237,10 +192,8 @@
Hc = Constants.Hc;
% used in other scripts not here!
-Gamma0 = Constants.Gamma0; % used in other scripts!
Gamma_w = Constants.Gamma_w; % used in other scripts!
c_i = Constants.c_i; % used in EnrgyPARM!
-RHO_bulk = Constants.RHO_bulk;
RTB = 1000; % initial root total biomass (g m-2)
% Rl used in ebal
@@ -416,8 +369,8 @@
[SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten);
%% get variables that are defined global and are used by other scripts
-global hd hh_frez POR KaT_Switch XSOC
-global XCAP SAVEhh COR CORh m n Alpha
+global hd hh_frez POR
+global XCAP m n Alpha
global Theta_s Theta_r Theta_f
% get soil constants
@@ -426,10 +379,7 @@
POR = SoilVariables.POR;
XK = SoilVariables.XK;
-KaT_Switch = SoilVariables.KaT_Switch;
-XSOC = SoilVariables.XSOC;
XCAP = SoilVariables.XCAP;
-SAVEhh = SoilVariables.SAVEhh;
Theta_s = VanGenuchten.Theta_s;
Theta_r = VanGenuchten.Theta_r;
Theta_f = VanGenuchten.Theta_f;
@@ -452,21 +402,16 @@
BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, landcoverClass(1));
%% get global vars
-global NBCh NBCT NBChB NBCTB BCh DSTOR0 NBChh NBCP BChB BCTB BCPB BCT BCP BtmPg
-NBCh = BoundaryCondition.NBCh;
+global NBCT NBCTB DSTOR0 NBCP BCTB BCPB BCT BCP BtmPg
NBCT = BoundaryCondition.NBCT;
-NBChB = BoundaryCondition.NBChB;
NBCTB = BoundaryCondition.NBCTB;
-BCh = BoundaryCondition.BCh;
DSTOR = BoundaryCondition.DSTOR;
DSTOR0 = BoundaryCondition.DSTOR0;
RS = BoundaryCondition.RS;
-NBChh = BoundaryCondition.NBChh;
DSTMAX = BoundaryCondition.DSTMAX;
IRPT1 = BoundaryCondition.IRPT1;
IRPT2 = BoundaryCondition.IRPT2;
NBCP = BoundaryCondition.NBCP;
-BChB = BoundaryCondition.BChB;
BCTB = BoundaryCondition.BCTB;
BCPB = BoundaryCondition.BCPB;
BCT = BoundaryCondition.BCT;
@@ -495,8 +440,6 @@
% Convert unit to Centimeter-Gram-Second system
% see issue 188 to refactor these lines
-HR_a = 0.01 .* (ForcingData.RH_msr);
-U = 100 .* (ForcingData.WS_msr);
TopPg = 100 .* (ForcingData.Pg_msr);
% the start of simulation period is from 0mins, while the input data start from 30mins.
@@ -505,6 +448,10 @@
TEND = TIME + TimeProperties.DELT * TimeProperties.Dur_tot; % Time to be reached at the end of simulation period
Delt_t0 = Delt_t; % Duration of last time step
TOLD_CRIT = [];
+
+% Srt, root water uptake;
+Srt = InitialValues.Srt; % will be updated!
+
for i = 1:1:TimeProperties.Dur_tot
KT = KT + 1; % Counting Number of timesteps
if KT > 1 && Delt_t > (TEND - TIME)
@@ -514,7 +461,7 @@
TimeStep(KT, 1) = Delt_t;
SUMTIME(KT, 1) = TIME;
Processing = TIME / TEND;
- NoTime(KT) = fix(SUMTIME(KT) / DELT);
+ NoTime(KT) = fix(SUMTIME(KT) / TimeProperties.DELT);
if NoTime(KT) == 0
k = 1;
else
@@ -543,7 +490,7 @@
P_gOLD(MN) = P_g(MN);
P_g(MN) = P_gg(MN);
end
- if rwuef == 1
+ if ModelSettings.rwuef == 1
SRT(MN, KT) = Srt(MN, 1);
end
end
@@ -679,31 +626,27 @@
if KT == 1
if isreal(fluxes.Actot) && isreal(thermal.Tsave) && isreal(fluxes.lEstot) && isreal(fluxes.lEctot)
Acc = fluxes.Actot;
- lEstot = fluxes.lEstot;
- lEctot = fluxes.lEctot;
Tss = thermal.Tsave;
else
Acc = 0;
- lEstot = 0;
- lEctot = 0;
+ fluxes.lEstot = 0;
+ fluxes.lEctot = 0;
Tss = ForcingData.Ta_msr(KT);
end
elseif NoTime(KT) > NoTime(KT - 1)
if isreal(fluxes.Actot) && isreal(thermal.Tsave) && isreal(fluxes.lEstot) && isreal(fluxes.lEctot)
Acc = fluxes.Actot;
- lEstot = fluxes.lEstot;
- lEctot = fluxes.lEctot;
Tss = thermal.Tsave;
else
Acc = 0;
- lEstot = 0;
- lEctot = 0;
+ fluxes.lEstot = 0;
+ fluxes.lEctot = 0;
Tss = ForcingData.Ta_msr(KT);
end
end
- sfactortot(KT) = sfactor;
DSTOR0 = DSTOR;
+ BoundaryCondition.DSTOR0 = DSTOR0;
if KT > 1
SoilVariables.XWRE = updateWettingHistory(SoilVariables, VanGenuchten);
@@ -721,13 +664,17 @@
end
hSAVE = hh(NN);
TSAVE = TT(NN);
- if NBCh == 1
- hN = BCh;
+
+ % set "hN" empty when the "if statement" below does not happen, see issue 98,
+ % item 5
+ hN = [];
+ if BoundaryCondition.NBCh == 1
+ hN = BoundaryCondition.BCh;
hh(NN) = hN;
hSAVE = hN;
- elseif NBCh == 2
- if NBChh ~= 2
- if BCh < 0
+ elseif BoundaryCondition.NBCh == 2
+ if BoundaryCondition.NBChh ~= 2
+ if BoundaryCondition.BCh < 0
hN = DSTOR0;
hh(NN) = hN;
hSAVE = hN;
@@ -738,7 +685,7 @@
end
end
else
- if NBChh ~= 2
+ if BoundaryCondition.NBChh ~= 2
hN = DSTOR0;
hh(NN) = hN;
hSAVE = hN;
@@ -747,10 +694,10 @@
Ts(KT) = Tss; % Tss is calculated above
Ta(KT) = ForcingData.Ta_msr(KT); % it is reset here because Ta is a gloval var
- Gvc(KT) = ForcingData.LAI_msr(KT); % it is reset here because Gvc is a gloval var
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for KIT = 1:NIT % Start the iteration procedure in a time step.
+
[TT_CRIT, hh_frez] = HT_frez(hh, T0, g, L_f, TT, NN, hd, Tmin);
% update inputs for UpdateSoilWaterContent
@@ -761,13 +708,14 @@
SoilVariables.KL_h = KL_h;
SoilVariables.KfL_h = KfL_h;
SoilVariables.TT = TT;
+ SoilVariables.T = T;
SoilVariables.h_frez = h_frez;
+
SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten);
+
% these can be removed after refactoring functions below
h = SoilVariables.h;
hh = SoilVariables.hh;
- COR = SoilVariables.COR;
- CORh = SoilVariables.CORh;
Theta_V = SoilVariables.Theta_V;
Theta_g = SoilVariables.Theta_g;
Theta_LL = SoilVariables.Theta_LL;
@@ -775,9 +723,7 @@
KL_h = SoilVariables.KL_h;
DTheta_LLh = SoilVariables.DTheta_LLh;
KfL_h = SoilVariables.KfL_h;
- KfL_T = SoilVariables.KfL_T;
hh_frez = SoilVariables.hh_frez;
- Theta_UU = SoilVariables.Theta_UU;
DTheta_UUh = SoilVariables.DTheta_UUh;
Theta_II = SoilVariables.Theta_II;
@@ -787,17 +733,15 @@
[RHOV, DRHOVh, DRHOVT] = Density_V(TT, hh, g, Rv, NN);
- % update inputs
- SoilVariables.Theta_L = Theta_L;
TransportCoefficient = conductivity.calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t);
W = TransportCoefficient.W;
WW = TransportCoefficient.WW;
D_Ta = TransportCoefficient.D_Ta;
[L] = Latent(TT, NN);
+ % DRHODAt unused!
[Xaa, XaT, Xah, DRHODAt, DRHODAz, RHODA] = Density_DA(T, RDA, P_g, Rv, DeltZ, h, hh, TT, P_gg, Delt_t, NL, NN, DRHOVT, DRHOVh, RHOV);
- Theta_LL = SoilVariables.Theta_LL;
ThermalConductivityCapacity = conductivity.calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV);
c_unsat = ThermalConductivityCapacity.c_unsat;
Lambda_eff = ThermalConductivityCapacity.Lambda_eff;
@@ -814,13 +758,36 @@
Beta_g = GasDispersivity.Beta_g;
DPgDZ = GasDispersivity.DPgDZ;
- run h_sub;
- if NBCh == 1
+ SoilVariables.Tss(KT) = Tss;
+ % After refactoring Enrgy_sub, the input/output of this function can be
+ % polished
+ % Srt is both input and output
+ % Replace run h_sub;
+ [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip] = soilmoisture.solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, ...
+ BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt);
+
+ DTheta_LLh = SoilVariables.DTheta_LLh;
+ DTheta_LLT = SoilVariables.DTheta_LLT;
+ DTheta_UUh = SoilVariables.DTheta_UUh;
+ SAVEDTheta_UUh = SoilVariables.SAVEDTheta_UUh;
+ SAVEDTheta_LLh = SoilVariables.SAVEDTheta_LLh;
+ QMB = HBoundaryFlux.QMB; % used in Enrgy_BC.m
+ QMT = HBoundaryFlux.QMT;
+ C5_a = HeatMatrices.C5_a;
+ C4 = HeatMatrices.C4;
+ hh = SoilVariables.hh;
+ Kha = HeatVariables.Kha;
+ KhT = HeatVariables.KhT;
+ Khh = HeatVariables.Khh;
+ Vvh = HeatVariables.Vvh;
+ VvT = HeatVariables.VvT;
+
+ if BoundaryCondition.NBCh == 1
DSTOR = 0;
RS = 0;
- elseif NBCh == 2
- AVAIL = -BCh;
- EXCESS = (AVAIL + QMT(KT)) * Delt_t;
+ elseif BoundaryCondition.NBCh == 2
+ AVAIL = -BoundaryCondition.BCh;
+ EXCESS = (AVAIL + QMT) * Delt_t;
if abs(EXCESS / Delt_t) <= 1e-10
EXCESS = 0;
end
@@ -828,7 +795,7 @@
RS = (EXCESS - DSTOR) / Delt_t;
else
AVAIL = AVAIL0 - Evap(KT);
- EXCESS = (AVAIL + QMT(KT)) * Delt_t;
+ EXCESS = (AVAIL + QMT) * Delt_t;
if abs(EXCESS / Delt_t) <= 1e-10
EXCESS = 0;
end
@@ -850,6 +817,7 @@
hSAVE = hh(NN);
TSAVE = TT(NN);
end
+
TIMEOLD = KT;
KIT;
KIT = 0;
@@ -861,14 +829,13 @@
SoilVariables.h = h;
SoilVariables.hh = hh;
SoilVariables.TT = TT;
+ SoilVariables.T = T;
SoilVariables.h_frez = h_frez;
SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten);
% these can be removed after refactoring codes below
h = SoilVariables.h;
hh = SoilVariables.hh;
- COR = SoilVariables.COR;
- CORh = SoilVariables.CORh;
Theta_V = SoilVariables.Theta_V;
Theta_g = SoilVariables.Theta_g;
Se = SoilVariables.Se;
@@ -876,7 +843,6 @@
Theta_LL = SoilVariables.Theta_LL;
DTheta_LLh = SoilVariables.DTheta_LLh;
KfL_h = SoilVariables.KfL_h;
- KfL_T = SoilVariables.KfL_T;
hh_frez = SoilVariables.hh_frez;
Theta_UU = SoilVariables.Theta_UU;
DTheta_UUh = SoilVariables.DTheta_UUh;
@@ -896,6 +862,11 @@
Theta_I(ML, ND) = Theta_II(ML, ND);
end
end
+ % update SoilVariables
+ SoilVariables.Theta_L = Theta_L;
+ SoilVariables.Theta_U = Theta_U;
+ SoilVariables.Theta_I = Theta_I;
+
run ObservationPoints;
end
if (TEND - TIME) < 1E-3
diff --git a/src/StartInit.m b/src/StartInit.m
index 7c6ec780..b1a6eaad 100644
--- a/src/StartInit.m
+++ b/src/StartInit.m
@@ -55,12 +55,12 @@
SoilVariables.KLT_Switch = 1;
SoilVariables.DVT_Switch = 1;
SoilVariables.KaT_Switch = [];
+ SoilVariables.DVa_Switch = [];
+ SoilVariables.KLa_Switch = [];
if ModelSettings.Soilairefc
SoilVariables.KaT_Switch = 1;
- % these vars are not used in the main script!
- Kaa_Switch = 1;
- DVa_Switch = 1;
- KLa_Switch = 1;
+ SoilVariables.DVa_Switch = 1;
+ SoilVariables.KLa_Switch = 1;
end
SoilVariables.Theta_I = Theta_I;
diff --git a/src/hPARM.m b/src/hPARM.m
deleted file mode 100644
index b60b60e0..00000000
--- a/src/hPARM.m
+++ /dev/null
@@ -1,95 +0,0 @@
-function [Chh, ChT, Khh, KhT, Kha, Vvh, VvT, Chg, DTheta_LLh, DTheta_LLT, DTheta_UUh, SAVEDTheta_UUh, SAVEDTheta_LLh] = hPARM(NL, hh, ...
- h, TT, T, Theta_LL, Theta_L, DTheta_LLh, DTheta_LLT, RHOV, RHOL, Theta_V, V_A, Eta, DRHOVh, ...
- DRHOVT, KfL_h, D_Ta, KL_T, D_V, D_Vg, COR, Beta_g, Gamma0, Gamma_w, KLa_Switch, DVa_Switch, ...
- hThmrl, Thmrlefc, nD, TT_CRIT, T0, h_frez, hh_frez, Theta_UU, Theta_U, CORh, DTheta_UUh, ...
- Chh, ChT, Khh, KhT, Kha, Vvh, VvT, Chg)
-
- % piecewise linear reduction function parameters of h;(Wesseling
- % 1991,Veenhof and McBride 1994)
- MN = 0;
- for ML = 1:NL
- for ND = 1:2
- MN = ML + ND - 1;
- if hThmrl
- CORhh = -1 * CORh(MN); % -0.0068;
- DhUU = COR(MN) * (hh(MN) + hh_frez(MN) - h(MN) - h_frez(MN) + (hh(MN) + hh_frez(MN)) * CORhh * (TT(MN) - T(MN)));
- DhU = COR(MN) * (hh(MN) - h(MN) + hh(MN) * CORhh * (TT(MN) - T(MN)));
- if DhU ~= 0 && DhUU ~= 0 && abs(Theta_LL(ML, ND) - Theta_L(ML, ND)) > 1e-6 && DTheta_UUh(ML, ND) ~= 0 % && abs(TT(MN)-T(MN))<1e-4
- DTheta_UUh(ML, ND) = (Theta_LL(ML, ND) - Theta_L(ML, ND)) * COR(MN) / DhUU;
- end
- if DhU ~= 0 && DhUU ~= 0 && abs(Theta_LL(ML, ND) - Theta_L(ML, ND)) > 1e-6 && DTheta_LLh(ML, ND) ~= 0 % && abs(TT(MN)-T(MN))<1e-4
- DTheta_LLh(ML, ND) = (Theta_UU(ML, ND) - Theta_U(ML, ND)) * COR(MN) / DhU;
- end
- %%
- DTheta_LLT(ML, ND) = DTheta_LLh(ML, ND) * (hh(MN) * CORhh);
- SAVEDTheta_LLT(ML, ND) = DTheta_LLT(ML, ND);
- SAVEDTheta_LLh(ML, ND) = DTheta_LLh(ML, ND);
- SAVEDTheta_UUh(ML, ND) = DTheta_UUh(ML, ND);
- else
- if abs(TT(MN) - T(MN)) > 1e-6
- DTheta_LLT(ML, ND) = DTheta_LLh(ML, ND) * (hh(MN) / Gamma0) * (0.1425 + 4.76 * 10^(-4) * TT(MN));
- else
- DTheta_LLT(ML, ND) = (Theta_LL(ML, ND) - Theta_L(ML, ND)) / (TT(MN) - T(MN));
- end
- end
- end
- end
-
- MN = 0;
- for ML = 1:NL
- for ND = 1:nD
- MN = ML + ND - 1;
- Chh(ML, ND) = (1 - RHOV(MN) / RHOL) * DTheta_LLh(ML, ND) + Theta_V(ML, ND) * DRHOVh(MN) / RHOL; % DTheta_LLh==>DTheta_UUh
- Khh(ML, ND) = (D_V(ML, ND) + D_Vg(ML)) * DRHOVh(MN) / RHOL + KfL_h(ML, ND); % KL_h==>KfL_h
- Chg(ML, ND) = KfL_h(ML, ND); % KL_h==>KfL_h
-
- if Thmrlefc == 1
- ChT(ML, ND) = (1 - RHOV(MN) / RHOL) * DTheta_LLT(ML, ND) + Theta_V(ML, ND) * DRHOVT(MN) / RHOL; % -RHOI/RHOL*KfL_T(ML,ND)*DTheta_LLh(ML,ND)/1e4DTheta_LLT==>DTheta_UUT
- KhT(ML, ND) = (D_V(ML, ND) * Eta(ML, ND) + D_Vg(ML)) * DRHOVT(MN) / RHOL + KL_T(ML, ND) + D_Ta(ML, ND); % -Khh(ML,ND)*1e4*L_f/(g*(T0+TT(MN)))*max(helpers.heaviside1(TT_CRIT(MN)-(TT(MN)+T0)))considering ice freezing pressure*KfL_T(ML,ND);%;%();% KfL_T considering the ice
- end
-
- if KLa_Switch == 1
- Kha(ML, ND) = RHOV(MN) * Beta_g(ML, ND) / RHOL + KfL_h(ML, ND) / Gamma_w; % KL_h==>KfL_h
- else
- Kha(ML, ND) = 0;
- end
-
- if DVa_Switch == 1
- Vvh(ML, ND) = -V_A(ML) * DRHOVh(MN) / RHOL;
- VvT(ML, ND) = -V_A(ML) * DRHOVT(MN) / RHOL;
- else
- Vvh(ML, ND) = 0;
- VvT(ML, ND) = 0;
- end
- if isnan(Chh(ML, ND))
- warning('\n Warning: FIX warning message \r');
- end
- if isnan(Khh(ML, ND)) == 1
- warning('\n Warning: FIX warning message \r');
- end
- if isnan(Chg(ML, ND))
- warning('\n Warning: FIX warning message\r');
- end
- if isnan(ChT(ML, ND))
- warning('\n Warning: FIX warning message \r');
- end
- if isnan(KhT(ML, ND))
- warning('\n Warning: FIX warning message \r');
- end
- if ~isreal(Chh(ML, ND))
- warning('\n Warning: FIX warning message \r');
- end
- if ~isreal(Khh(ML, ND)) == 1
- warning('\n Warning: FIX warning message \r');
- end
- if ~isreal(Chg(ML, ND))
- warning('\n Warning: FIX warning message \r');
- end
- if ~isreal(ChT(ML, ND))
- warning('\n Warning: FIX warning message \r');
- end
- if ~isreal(KhT(ML, ND))
- warning('\n Warning: FIX warning message \r');
- end
- end
- end
diff --git a/src/h_BC.m b/src/h_BC.m
deleted file mode 100644
index 284ced98..00000000
--- a/src/h_BC.m
+++ /dev/null
@@ -1,66 +0,0 @@
-function [AVAIL0, RHS, C4, C4_a, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Resis_a, Srt, Precip] = h_BC(RHS, NBCh, NBChB, BCh, BChB, hN, KT, Delt_t, DSTOR0, NBChh, TIME, h_SUR, C4, KL_h, NN, C4_a, DeltZ, RHOV, Ta, HR_a, U, Theta_LL, Ts, Rv, g, NL, hh, rwuef, Theta_UU, Rn, T, TT, Gvc, Rns, Srt, Precip_msr, SUMTIME) % [AVAIL0,RHS,C4,C4_a,Evap,EVAP,Trap,Precip,bx,r_a_SOIL,Resis_a]=h_BC(DeltZ,bx,RHS,NBCh,NBChB,BCh,BChB,hN,KT,Delt_t,DSTOR0,NBChh,TIME,h_SUR,C4,KL_h,Precip,NN,AVAIL0,C4_a,Evap,RHOV,Ta,HR_a,U,Ts,Theta_LL,Rv,g,NL,hh,rwuef,Theta_UU,Rn,T,TT,Gvc,Rns)
- global DELT Ks0 theta_s0
-
- Precipp = 0;
- %%%%%%%%%% Apply the bottom boundary condition called for by NBChB %%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if NBChB == 1 % -----> Specify matric head at bottom to be ---BChB;
- RHS(1) = BChB;
- C4(1, 1) = 1;
- RHS(2) = RHS(2) - C4(1, 2) * RHS(1);
- C4(1, 2) = 0;
- C4_a(1) = 0;
- elseif NBChB == 2 % -----> Specify flux at bottom to be ---BChB (Positive upwards);
- RHS(1) = RHS(1) + BChB;
- elseif NBChB == 3 % -----> NBChB=3,Gravity drainage at bottom--specify flux= hydraulic conductivity;
- RHS(1) = RHS(1) - KL_h(1, 1);
- end
-
- %%%%%%%%%% Apply the surface boundary condition called for by NBCh %%%%%%%
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if NBCh == 1 % -----> Specified matric head at surface---equal to hN;
- RHS(NN) = h_SUR(KT);
- C4(NN, 1) = 1;
- RHS(NN - 1) = RHS(NN - 1) - C4(NN - 1, 2) * RHS(NN);
- C4(NN - 1, 2) = 0;
- C4_a(NN - 1) = 0;
- elseif NBCh == 2
- if NBChh == 1
- RHS(NN) = hN;
- C4(NN, 1) = 1;
- RHS(NN - 1) = RHS(NN - 1) - C4(NN - 1, 2) * RHS(NN);
- C4(NN - 1, 2) = 0;
- else
- RHS(NN) = RHS(NN) - BCh; % > and a specified matric head (saturation or dryness)was applied;
- end
- else
- [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Resis_a, Srt] = Evap_Cal(DeltZ, TIME, RHOV, Ta, HR_a, U, Theta_LL, Ts, Rv, g, NL, NN, KT, hh, rwuef, Theta_UU, Rn, T, TT, Gvc, Rns, Srt);
-
- Precip_msr(KT) = min(Precip_msr(KT), Ks0 / (3600 * 24) * DELT * 10);
- Precip_msr(KT) = min(Precip_msr(KT), theta_s0 * 50 - DeltZ(51:54) * Theta_UU(51:54, 1) * 10);
-
- if Ts(KT) > 0
- Precip(KT) = Precip_msr(KT) * 0.1 / DELT;
- else
- Precip(KT) = Precip_msr(KT) * 0.1 / DELT;
- Precipp = Precipp + Precip(KT);
- Precip(KT) = 0;
- end
-
- if Ts(KT) > 0
- AVAIL0 = Precip(KT) + Precipp + DSTOR0 / Delt_t;
- Precipp = 0;
- else
- AVAIL0 = Precip(KT) + DSTOR0 / Delt_t;
- end
-
- if NBChh == 1
- RHS(NN) = hN;
- C4(NN, 1) = 1;
- RHS(NN - 1) = RHS(NN - 1) - C4(NN - 1, 2) * RHS(NN);
- C4(NN - 1, 2) = 0;
- C4_a(NN - 1) = 0;
- else
- RHS(NN) = RHS(NN) + AVAIL0 - Evap(KT);
- end
- end
diff --git a/src/h_Bndry_Flux.m b/src/h_Bndry_Flux.m
deleted file mode 100644
index 781d53e5..00000000
--- a/src/h_Bndry_Flux.m
+++ /dev/null
@@ -1,4 +0,0 @@
-function [QMT, QMB] = h_Bndry_Flux(SAVE, hh, NN, KT)
-
- QMT(KT) = SAVE(2, 1, 1) - SAVE(2, 2, 1) * hh(NN - 1) - SAVE(2, 3, 1) * hh(NN);
- QMB(KT) = -SAVE(1, 1, 1) + SAVE(1, 2, 1) * hh(1) + SAVE(1, 3, 1) * hh(2);
diff --git a/src/h_EQ.m b/src/h_EQ.m
deleted file mode 100644
index 81923047..00000000
--- a/src/h_EQ.m
+++ /dev/null
@@ -1,67 +0,0 @@
-function [RHS, C4, SAVE] = h_EQ(C1, C2, C4, C5, C6, C7, C5_a, C9, NL, NN, Delt_t, T, TT, h, P_gg, Thmrlefc, Soilairefc)
-
- if Thmrlefc && ~Soilairefc
- RHS(1) = -C9(1) - C7(1) + (C1(1, 1) * h(1) + C1(1, 2) * h(2)) / Delt_t ...
- - (C2(1, 1) / Delt_t + C5(1, 1)) * TT(1) - (C2(1, 2) / Delt_t + C5(1, 2)) * TT(2) ...
- + (C2(1, 1) / Delt_t) * T(1) + (C2(1, 2) / Delt_t) * T(2);
- for ML = 2:NL
- ARG1 = C2(ML - 1, 2) / Delt_t;
- ARG2 = C2(ML, 1) / Delt_t;
- ARG3 = C2(ML, 2) / Delt_t;
-
- RHS(ML) = -C9(ML) - C7(ML) + (C1(ML - 1, 2) * h(ML - 1) + C1(ML, 1) * h(ML) + C1(ML, 2) * h(ML + 1)) / Delt_t ...
- - (ARG1 + C5(ML - 1, 2)) * TT(ML - 1) - (ARG2 + C5(ML, 1)) * TT(ML) - (ARG3 + C5(ML, 2)) * TT(ML + 1) ...
- + ARG1 * T(ML - 1) + ARG2 * T(ML) + ARG3 * T(ML + 1);
- end
- RHS(NN) = -C9(NN) - C7(NN) + (C1(NN - 1, 2) * h(NN - 1) + C1(NN, 1) * h(NN)) / Delt_t ...
- - (C2(NN - 1, 2) / Delt_t + C5(NN - 1, 2)) * TT(NN - 1) - (C2(NN, 1) / Delt_t + C5(NN, 1)) * TT(NN) ...
- + (C2(NN - 1, 2) / Delt_t) * T(NN - 1) + (C2(NN, 1) / Delt_t) * T(NN);
- elseif ~Thmrlefc && Soilairefc
- RHS(1) = -C9(1) - C7(1) + (C1(1, 1) * h(1) + C1(1, 2) * h(2)) / Delt_t ...
- - C6(1, 1) * P_gg(1) - C6(1, 2) * P_gg(2);
- for ML = 2:NL
- RHS(ML) = -C9(ML) - C7(ML) + (C1(ML - 1, 2) * h(ML - 1) + C1(ML, 1) * h(ML) + C1(ML, 2) * h(ML + 1)) / Delt_t ...
- - C6(ML - 1, 2) * P_gg(ML - 1) - C6(ML, 1) * P_gg(ML) - C6(ML, 2) * P_gg(ML + 1);
- end
- RHS(NN) = -C9(NN) - C7(NN) + (C1(NN - 1, 2) * h(NN - 1) + C1(NN, 1) * h(NN)) / Delt_t ...
- - C6(NN - 1, 2) * P_gg(NN - 1) - C6(NN, 1) * P_gg(NN);
- elseif Thmrlefc && Soilairefc
- RHS(1) = -C9(1) - C7(1) + (C1(1, 1) * h(1) + C1(1, 2) * h(2)) / Delt_t ...
- - (C2(1, 1) / Delt_t + C5(1, 1)) * TT(1) - (C2(1, 2) / Delt_t + C5(1, 2)) * TT(2) ...
- - C6(1, 1) * P_gg(1) - C6(1, 2) * P_gg(2) ...
- + (C2(1, 1) / Delt_t) * T(1) + (C2(1, 2) / Delt_t) * T(2);
- for ML = 2:NL
- ARG1 = C2(ML - 1, 2) / Delt_t;
- ARG2 = C2(ML, 1) / Delt_t;
- ARG3 = C2(ML, 2) / Delt_t;
-
- RHS(ML) = -C9(ML) - C7(ML) + (C1(ML - 1, 2) * h(ML - 1) + C1(ML, 1) * h(ML) + C1(ML, 2) * h(ML + 1)) / Delt_t ...
- - (ARG1 + C5_a(ML - 1)) * TT(ML - 1) - (ARG2 + C5(ML, 1)) * TT(ML) - (ARG3 + C5(ML, 2)) * TT(ML + 1) ...
- - C6(ML - 1, 2) * P_gg(ML - 1) - C6(ML, 1) * P_gg(ML) - C6(ML, 2) * P_gg(ML + 1) ...
- + ARG1 * T(ML - 1) + ARG2 * T(ML) + ARG3 * T(ML + 1);
- end
-
- RHS(NN) = -C9(NN) - C7(NN) + (C1(NN - 1, 2) * h(NN - 1) + C1(NN, 1) * h(NN)) / Delt_t ...
- - (C2(NN - 1, 2) / Delt_t + C5_a(NN - 1)) * TT(NN - 1) - (C2(NN, 1) / Delt_t + C5(NN, 1)) * TT(NN) ...
- - C6(NN - 1, 2) * P_gg(NN - 1) - C6(NN, 1) * P_gg(NN) ...
- + (C2(NN - 1, 2) / Delt_t) * T(NN - 1) + (C2(NN, 1) / Delt_t) * T(NN);
- else
- RHS(1) = -C9(1) - C7(1) + (C1(1, 1) * h(1) + C1(1, 2) * h(2)) / Delt_t;
- for ML = 2:NL
- RHS(ML) = -C9(ML) - C7(ML) + (C1(ML - 1, 2) * h(ML - 1) + C1(ML, 1) * h(ML) + C1(ML, 2) * h(ML + 1)) / Delt_t;
- end
- RHS(NN) = -C9(NN) - C7(NN) + (C1(NN - 1, 2) * h(NN - 1) + C1(NN, 1) * h(NN)) / Delt_t;
- end
-
- for MN = 1:NN
- for ND = 1:2
- C4(MN, ND) = C1(MN, ND) / Delt_t + C4(MN, ND);
- end
- end
-
- SAVE(1, 1, 1) = RHS(1);
- SAVE(1, 2, 1) = C4(1, 1);
- SAVE(1, 3, 1) = C4(1, 2);
- SAVE(2, 1, 1) = RHS(NN);
- SAVE(2, 2, 1) = C4(NN - 1, 2);
- SAVE(2, 3, 1) = C4(NN, 1);
diff --git a/src/h_MAT.m b/src/h_MAT.m
deleted file mode 100644
index 121a8ba9..00000000
--- a/src/h_MAT.m
+++ /dev/null
@@ -1,54 +0,0 @@
-function [C1, C2, C4, C3, C4_a, C5, C6, C7, C5_a, C9] = h_MAT(Chh, ChT, Khh, KhT, Kha, Vvh, VvT, Chg, DeltZ, NL, NN, Srt)
-
- for MN = 1:NN % Clean the space in C1-7 every iteration,otherwise, in *.PARM files,
- for ND = 1:2 % C1-7 will be mixed up with pre-storaged data, which will cause extremly crazy for computation, which exactly results in NAN.
- C1(MN, ND) = 0;
- C7(MN) = 0;
- C9(MN) = 0; % C9 is the matrix coefficient of root water uptake;
- C4(MN, ND) = 0;
- C4_a(MN) = 0;
- C5_a(MN) = 0;
- C2(MN, ND) = 0;
- C3(MN, ND) = 0;
- C5(MN, ND) = 0;
- C6(MN, ND) = 0;
- end
- end
- for ML = 1:NL
- C1(ML, 1) = C1(ML, 1) + Chh(ML, 1) * DeltZ(ML) / 2;
- C1(ML + 1, 1) = C1(ML + 1, 1) + Chh(ML, 2) * DeltZ(ML) / 2; %
-
- C2(ML, 1) = C2(ML, 1) + ChT(ML, 1) * DeltZ(ML) / 2;
- C2(ML + 1, 1) = C2(ML + 1, 1) + ChT(ML, 2) * DeltZ(ML) / 2; %
-
- C4ARG1 = (Khh(ML, 1) + Khh(ML, 2)) / (2 * DeltZ(ML)); % sqrt(Khh(ML,1)*Khh(ML,2))/(DeltZ(ML));%
- C4ARG2_1 = Vvh(ML, 1) / 3 + Vvh(ML, 2) / 6;
- C4ARG2_2 = Vvh(ML, 1) / 6 + Vvh(ML, 2) / 3;
- C4(ML, 1) = C4(ML, 1) + C4ARG1 - C4ARG2_1;
- C4(ML, 2) = C4(ML, 2) - C4ARG1 - C4ARG2_2;
- C4(ML + 1, 1) = C4(ML + 1, 1) + C4ARG1 + C4ARG2_2;
- C4_a(ML) = -C4ARG1 + C4ARG2_1;
-
- C5ARG1 = (KhT(ML, 1) + KhT(ML, 2)) / (2 * DeltZ(ML)); % sqrt(KhT(ML,1)*KhT(ML,2))/(DeltZ(ML));%
- C5ARG2_1 = VvT(ML, 1) / 3 + VvT(ML, 2) / 6;
- C5ARG2_2 = VvT(ML, 1) / 6 + VvT(ML, 2) / 3;
- C5(ML, 1) = C5(ML, 1) + C5ARG1 - C5ARG2_1;
- C5(ML, 2) = C5(ML, 2) - C5ARG1 - C5ARG2_2;
- C5(ML + 1, 1) = C5(ML + 1, 1) + C5ARG1 + C5ARG2_2;
- C5_a(ML) = -C5ARG1 + C5ARG2_1;
-
- C6ARG = (Kha(ML, 1) + Kha(ML, 2)) / (2 * DeltZ(ML)); % sqrt(Kha(ML,1)*Kha(ML,2))/(DeltZ(ML));%
- C6(ML, 1) = C6(ML, 1) + C6ARG;
- C6(ML, 2) = C6(ML, 2) - C6ARG;
- C6(ML + 1, 1) = C6(ML + 1, 1) + C6ARG;
-
- C7ARG = (Chg(ML, 1) + Chg(ML, 2)) / 2; % sqrt(Chg(ML,1)*Chg(ML,2));%
- C7(ML) = C7(ML) - C7ARG;
- C7(ML + 1) = C7(ML + 1) + C7ARG;
-
- % Srt, root water uptake;
- C9ARG1 = (2 * Srt(ML, 1) + Srt(ML, 2)) * DeltZ(ML) / 6; % sqrt(Chg(ML,1)*Chg(ML,2));%
- C9ARG2 = (Srt(ML, 1) + 2 * Srt(ML, 2)) * DeltZ(ML) / 6;
- C9(ML) = C9(ML) + C9ARG1;
- C9(ML + 1) = C9(ML + 1) + C9ARG2;
- end
diff --git a/src/h_sub.m b/src/h_sub.m
deleted file mode 100644
index 56f474ca..00000000
--- a/src/h_sub.m
+++ /dev/null
@@ -1,49 +0,0 @@
-% function [Chh,ChT,Khh,KhT,Kha,Vvh,VvT,Chg,DTheta_LLh,DTheta_LLT,DTheta_UUh,SAVEDTheta_UUh,SAVEDTheta_LLh,C1,C2,C4,C3,C4_a,C5,C6,C7,C5_a,C9,RHS,SAVE,AVAIL0,...
-% Rn_SOIL,Evap,EVAP,Trap,r_a_SOIL,Resis_a,Srt,Precip,...
-% CHK,hh,SAVEhh,QMT,QMB]= h_sub_SPAT(NL,...
-% h,TT,T,Theta_LL,Theta_L,DTheta_LLh,DTheta_LLT,RHOV,RHOL,Theta_V,V_A,Eta,DRHOVh,...
-% DRHOVT,KfL_h,D_Ta,KL_T,D_V,D_Vg,COR,Beta_g,Gamma0,Gamma_w,KLa_Switch,DVa_Switch,hThmrl,Thmrlefc,nD,TT_CRIT,T0,h_frez,hh_frez,Theta_UU,Theta_U,CORh,DTheta_UUh,...
-% Chh,ChT,Khh,KhT,DeltZ,Delt_t,P_gg,Soilairefc,...
-% NBCh,NBChB,BCh,BChB,hN,KT,DSTOR0,NBChh,TIME,h_SUR,C4,KL_h,Ta,HR_a,U,Ts,Rv,g,hh,rwuef,Rn,Gvc,Rns,Srt,Precip_msr,SUMTIME)
-
-function h_sub
- global hh MN NN
- global C1 C2 C4 C3 C4_a C5 C6 C7
- global Chh ChT Khh KhT Kha Vvh VvT Chg DeltZ C5_a
- global NL nD bx
- global Delt_t RHS T TT h P_gg SAVE Thmrlefc Soilairefc
- global RHOL Gamma_w DTheta_LLh DTheta_LLT
- global Theta_L Theta_LL Theta_V Eta V_A
- global RHOV DRHOVh DRHOVT KL_h D_Ta KL_T D_V D_Vg
- global COR hThmrl Beta_g Gamma0 KLa_Switch DVa_Switch
- global Precip Evap CHK EVAP
- global NBCh NBChB BCh BChB hN KT DSTOR0 NBChh TIME h_SUR AVAIL0
- global QMT QMB QMTT QMBB Evapo trap RnSOIL PrecipO
- global Ta U Ts Rv g HR_a Srt C9 % U_wind is the mean wind speed at height z_ref (m/s^-1), U is the wind speed at each time step.
- global rwuef Trap Rn_SOIL hOLD
- global TT_CRIT T0 h_frez hh_frez Theta_UU Theta_U DTheta_UUh KfL_h Rn SAVEhh CORh r_a_SOIL Resis_a Precip_msr SUMTIME Gvc SAVEDTheta_UUh SAVEDTheta_LLh Rns
-
- [Chh, ChT, Khh, KhT, Kha, Vvh, VvT, Chg, DTheta_LLh, DTheta_LLT, DTheta_UUh, SAVEDTheta_UUh, SAVEDTheta_LLh] = hPARM(NL, hh, ...
- h, TT, T, Theta_LL, Theta_L, DTheta_LLh, DTheta_LLT, RHOV, RHOL, Theta_V, V_A, Eta, DRHOVh, ...
- DRHOVT, KfL_h, D_Ta, KL_T, D_V, D_Vg, COR, Beta_g, Gamma0, Gamma_w, KLa_Switch, DVa_Switch, hThmrl, Thmrlefc, nD, TT_CRIT, T0, h_frez, hh_frez, Theta_UU, Theta_U, CORh, DTheta_UUh, Chh, ChT, Khh, KhT, Kha, Vvh, VvT, Chg);
- [C1, C2, C4, C3, C4_a, C5, C6, C7, C5_a, C9] = h_MAT(Chh, ChT, Khh, KhT, Kha, Vvh, VvT, Chg, DeltZ, NL, NN, Srt);
- [RHS, C4, SAVE] = h_EQ(C1, C2, C4, C5, C6, C7, C5_a, C9, NL, NN, Delt_t, T, TT, h, P_gg, Thmrlefc, Soilairefc);
- [AVAIL0, RHS, C4, C4_a, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Resis_a, Srt, Precip] = h_BC(RHS, NBCh, NBChB, BCh, BChB, hN, KT, Delt_t, DSTOR0, NBChh, TIME, h_SUR, C4, KL_h, NN, C4_a, DeltZ, RHOV, Ta, HR_a, U, Theta_LL, Ts, Rv, g, NL, hh, rwuef, Theta_UU, Rn, T, TT, Gvc, Rns, Srt, Precip_msr, SUMTIME); % h_BC;
- [CHK, hh, C4, SAVEhh] = hh_Solve(C4, hh, NN, NL, C4_a, RHS);
- if any(isnan(hh)) || any(hh(1:NN) <= -1E12)
- for MN = 1:NN
- hh(MN) = hOLD(MN);
- end
- end
- for MN = 1:NN
- if hh(MN) >= -1e-6
- hh(MN) = -1e-6;
- end
- end
- [QMT, QMB] = h_Bndry_Flux(SAVE, hh, NN, KT);
- PrecipO(KT) = Precip(KT);
- QMTT(KT) = QMT(KT);
- QMBB(KT) = QMB(KT);
- Evapo(KT) = Evap(KT);
- trap(KT) = Trap(KT);
- RnSOIL(KT) = Rn_SOIL(KT);
diff --git a/src/hh_Solve.m b/src/hh_Solve.m
deleted file mode 100644
index 1430c693..00000000
--- a/src/hh_Solve.m
+++ /dev/null
@@ -1,32 +0,0 @@
-function [CHK, hh, C4, SAVEhh] = hh_Solve(C4, hh, NN, NL, C4_a, RHS)
-
- TTheta_UU = zeros(NL + 1, 2);
- hdry(1:NN) = -5e5;
- hwet(1:NN) = -1;
- DTTheta_UU = zeros(NL + 1, 2);
- TGama_hh = zeros(NL + 1, 1);
- TTheta_m = zeros(NL, 1);
- RHS(1) = RHS(1) / C4(1, 1);
- SAVDTheta_LLh = zeros(NL + 1, 2);
- INDICATOR = 0;
-
- for ML = 2:NN
- C4(ML, 1) = C4(ML, 1) - C4_a(ML - 1) * C4(ML - 1, 2) / C4(ML - 1, 1);
- RHS(ML) = (RHS(ML) - C4_a(ML - 1) * RHS(ML - 1)) / C4(ML, 1);
- end
-
- for ML = NL:-1:1
- RHS(ML) = RHS(ML) - C4(ML, 2) * RHS(ML + 1) / C4(ML, 1);
- end
- for MN = 1:NN
- CHK(MN) = abs(RHS(MN) - hh(MN));
- hh(MN) = RHS(MN);
- SAVEhh(MN) = hh(MN);
- end
-
- if isnan(SAVEhh) == 1
- warning('\n case "isnan(SAVEhh) == 1" happens, dont know what to do!\r');
- end
- if ~isreal(SAVEhh)
- warning('\n case "~isreal(SAVEhh)" happens, dont know what to do! \r');
- end