Skip to content

Commit

Permalink
Merge branch 'main' into fix_issue_99
Browse files Browse the repository at this point in the history
  • Loading branch information
@SarahAlidoost
SarahAlidoost committed Sep 12, 2023
2 parents d878304 + 5decaae commit 64d60b5
Show file tree
Hide file tree
Showing 25 changed files with 577 additions and 821 deletions.
3 changes: 3 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
Binary file modified run_model_on_snellius/exe/STEMMUS_SCOPE
View file
Binary file not shown.
3 changes: 3 additions & 0 deletions src/+init/defineInitialValues.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -236,4 +236,7 @@
InitialValues.(field{1}) = structure.(field{1});
end
end

% Arraies for calculating boundary flux;
InitialValues.SAVE = zeros(3, 3, 3);
end
4 changes: 4 additions & 0 deletions src/+io/define_constants.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
3 changes: 0 additions & 3 deletions src/+io/getModelSettings.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down
72 changes: 72 additions & 0 deletions src/+soilmoisture/assembleCoefficientMatrices.m
View file
Original file line number Diff line number Diff line change
@@ -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
11 changes: 11 additions & 0 deletions src/+soilmoisture/calculateAerodynamicResistance.m
View file
Original file line number Diff line number Diff line change
@@ -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
79 changes: 79 additions & 0 deletions src/+soilmoisture/calculateBoundaryConditions.m
View file
Original file line number Diff line number Diff line change
@@ -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
37 changes: 37 additions & 0 deletions src/+soilmoisture/calculateEvapotranspiration.m
View file
Original file line number Diff line number Diff line change
@@ -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
98 changes: 98 additions & 0 deletions src/+soilmoisture/calculateMatricCoefficients.m
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 64d60b5

Please sign in to comment.