Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix issue 100 #197

Merged
merged 17 commits into from
Oct 10, 2023
Merged

Fix issue 100 #197

Show file tree
Hide file tree
Changes from 16 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
ab0c70b
rename all Engry scripts and move them to energy folder
SarahAlidoost Oct 2, 2023
8c8e314
rename functions and add doc
SarahAlidoost Oct 3, 2023
4a79748
Merge branch 'main' into fix_issue_100
SarahAlidoost Oct 3, 2023
f2fa459
refactor EngryParm
SarahAlidoost Oct 3, 2023
ea04820
refcator Engry_MAT
SarahAlidoost Oct 3, 2023
694462e
refactor Engry_EQ and Engry_BC
SarahAlidoost Oct 3, 2023
7d7b5c8
refactor Enrgy_sub and its functions
SarahAlidoost Oct 3, 2023
bcabd97
fix the inputs
SarahAlidoost Oct 3, 2023
29ca7f6
remove global vars and unused vars from main script
SarahAlidoost Oct 3, 2023
ddcc80c
fix linter errors
SarahAlidoost Oct 3, 2023
1d3b3f8
address todos
SarahAlidoost Oct 3, 2023
a0120de
debug the code
SarahAlidoost Oct 3, 2023
7085981
debug the code, polish main script
SarahAlidoost Oct 6, 2023
771497e
remove unused loop from previous refactoring
SarahAlidoost Oct 6, 2023
de756a5
add changes to changelog
SarahAlidoost Oct 6, 2023
5faac0e
regenerate exe file
SarahAlidoost Oct 6, 2023
ee139b1
add reference to docstrings
SarahAlidoost Oct 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,13 +8,15 @@
[#193](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/193)
- All `Air_*` functions are refcatored and moved to `+dryair` folder in
[194](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/194)
- All `Engry_*` functions are refcatored and moved to `+energy` folder in
[197](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/197)

**Fixed:**

- issue [#181](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/181)
- issue [#98](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/98)
- issue [#99](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/99)

- issue [#100](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/100)

<a name="1.3.0"></a>
# [1.3.0](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.3.0) - 22 Jun 2023
Expand Down
Binary file modified run_model_on_snellius/exe/STEMMUS_SCOPE
View file
Open in desktop
Binary file not shown.
90 changes: 90 additions & 0 deletions src/+energy/assembleCoefficientMatrices.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
function [RHS, EnergyMatrices, SAVE] = assembleCoefficientMatrices(EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg)
%{
assembles the coefficient matrices of Equation 4.32, STEMMUS Technical
Notes, page 44, the example was only shown for the soil moisture
equation, but here it is for the energy equation.
%}
C1 = EnergyMatrices.C1;
C2 = EnergyMatrices.C2;
C3 = EnergyMatrices.C3;
C4 = EnergyMatrices.C4;
C4_a = EnergyMatrices.C4_a;
C5 = EnergyMatrices.C5;
C5_a = EnergyMatrices.C5_a;
C6 = EnergyMatrices.C6;
C7 = EnergyMatrices.C7;

ModelSettings = io.getModelSettings();
n = ModelSettings.NN;

% Alias of SoilVariables
SV = SoilVariables;

if ModelSettings.Soilairefc && ModelSettings.Thmrlefc
RHS(1) = -C7(1) + (C2(1, 1) * SV.T(1) + C2(1, 2) * SV.T(2)) / Delt_t ...
- (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
- (C3(1, 1) / Delt_t + C6(1, 1)) * P_gg(1) - (C3(1, 2) / Delt_t + C6(1, 2)) * P_gg(2) ...
+ (C3(1, 1) / Delt_t) * P_g(1) + (C3(1, 2) / Delt_t) * P_g(2) ...
+ (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);

for i = 2:ModelSettings.NL
ARG1 = C3(i - 1, 2) / Delt_t;
ARG2 = C3(i, 1) / Delt_t;
ARG3 = C3(i, 2) / Delt_t;

ARG4 = C1(i - 1, 2) / Delt_t;
ARG5 = C1(i, 1) / Delt_t;
ARG6 = C1(i, 2) / Delt_t;

RHS(i) = -C7(i) + (C2(i - 1, 2) * SV.T(i - 1) + C2(i, 1) * SV.T(i) + C2(i, 2) * SV.T(i + 1)) / Delt_t ...
- (ARG1 + C6_a(i - 1)) * P_gg(i - 1) - (ARG2 + C6(i, 1)) * P_gg(i) - (ARG3 + C6(i, 2)) * P_gg(i + 1) ...
- (ARG4 + C4_a(i - 1)) * SV.hh(i - 1) - (ARG5 + C4(i, 1)) * SV.hh(i) - (ARG6 + C4(i, 2)) * SV.hh(i + 1) ...
+ ARG1 * P_g(i - 1) + ARG2 * P_g(i) + ARG3 * P_g(i + 1) ...
+ ARG4 * SV.h(i - 1) + ARG5 * SV.h(i) + ARG6 * SV.h(i + 1);
end

RHS(n) = -C7(n) + (C2(n - 1, 2) * SV.T(n - 1) + C2(n, 1) * SV.T(n)) / Delt_t ...
- (C3(n - 1, 2) / Delt_t + C6_a(n - 1)) * P_gg(n - 1) - (C3(n, 1) / Delt_t + C6(n, 1)) * P_gg(n) ...
- (C1(n - 1, 2) / Delt_t + C4_a(n - 1)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
+ (C3(n - 1, 2) / Delt_t) * P_g(n - 1) + (C3(n, 1) / Delt_t) * P_g(n) ...
+ (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
elseif ~ModelSettings.Soilairefc && ModelSettings.Thmrlefc
RHS(1) = -C7(1) + (C2(1, 1) * SV.T(1) + C2(1, 2) * SV.T(2)) / Delt_t ...
- (C1(1, 1) / Delt_t + C4(1, 1)) * SV.hh(1) - (C1(1, 2) / Delt_t + C4(1, 2)) * SV.hh(2) ...
+ (C1(1, 1) / Delt_t) * SV.h(1) + (C1(1, 2) / Delt_t) * SV.h(2);
for i = 2:ModelSettings.NL
ARG4 = C1(i - 1, 2) / Delt_t;
ARG5 = C1(i, 1) / Delt_t;
ARG6 = C1(i, 2) / Delt_t;

RHS(i) = -C7(i) + (C2(i - 1, 2) * SV.T(i - 1) + C2(i, 1) * SV.T(i) + C2(i, 2) * SV.T(i + 1)) / Delt_t ...
- (ARG4 + C4(i - 1, 2)) * SV.hh(i - 1) - (ARG5 + C4(i, 1)) * SV.hh(i) - (ARG6 + C4(i, 2)) * SV.hh(i + 1) ...
+ ARG4 * SV.h(i - 1) + ARG5 * SV.h(i) + ARG6 * SV.h(i + 1);
end

RHS(n) = -C7(n) + (C2(n - 1, 2) * SV.T(n - 1) + C2(n, 1) * SV.T(n)) / Delt_t ...
- (C1(n - 1, 2) / Delt_t + C4(n - 1, 2)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
+ (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
else
RHS(1) = -C7(1) + (C2(1, 1) * SV.T(1) + C2(1, 2) * SV.T(2)) / Delt_t;
for i = 2:ModelSettings.NL
RHS(i) = -C7(i) + (C2(i - 1, 2) * SV.T(i - 1) + C2(i, 1) * SV.T(i) + C2(i, 2) * SV.T(i + 1)) / Delt_t;
end
RHS(n) = -C7(n) + (C2(n - 1, 2) * SV.T(n - 1) + C2(n, 1) * SV.T(n)) / Delt_t;
end

for i = 1:ModelSettings.NN
for j = 1:ModelSettings.nD
C5(i, j) = C2(i, j) / Delt_t + C5(i, j);
end
end

EnergyMatrices.C5 = C5;

SAVE(1, 1, 2) = RHS(1);
SAVE(1, 2, 2) = C5(1, 1);
SAVE(1, 3, 2) = C5(1, 2);
SAVE(2, 1, 2) = RHS(n);
SAVE(2, 2, 2) = C5(n - 1, 2);
SAVE(2, 3, 2) = C5(n, 1);
end
43 changes: 43 additions & 0 deletions src/+energy/calculateBoundaryConditions.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
function [RHS, EnergyMatrices] = calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, ...
SoilVariables, Precip, EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT)
%{
Determine the boundary condition for solving the energy equation.
SarahAlidoost marked this conversation as resolved.
Show resolved Hide resolved
%}

ModelSettings = io.getModelSettings();
n = ModelSettings.NN;

Constants = io.define_constants();

% Apply the bottom boundary condition called for by BoundaryCondition.NBCTB
if BoundaryCondition.NBCTB == 1
RHS(1) = BoundaryCondition.BCTB;
EnergyMatrices.C5(1, 1) = 1;
RHS(2) = RHS(2) - EnergyMatrices.C5(1, 2) * RHS(1);
EnergyMatrices.C5(1, 2) = 0;
EnergyMatrices.C5_a(1) = 0;
elseif BoundaryCondition.NBCTB == 2
RHS(1) = RHS(1) + BoundaryCondition.BCTB;
else
EnergyMatrices.C5(1, 1) = EnergyMatrices.C5(1, 1) - Constants.c_L * Constants.RHOL * HBoundaryFlux.QMB;
end

% Apply the surface boundary condition called by BoundaryCondition.NBCT
if BoundaryCondition.NBCT == 1
if isreal(SoilVariables.Tss(KT))
RHS(n) = SoilVariables.Tss(KT);
else
RHS(n) = ForcingData.Ta_msr(KT);
end
EnergyMatrices.C5(n, 1) = 1;
RHS(n - 1) = RHS(n - 1) - EnergyMatrices.C5(n - 1, 2) * RHS(n);
EnergyMatrices.C5(n - 1, 2) = 0;
EnergyMatrices.C5_a(n - 1) = 0;
elseif BoundaryCondition.NBCT == 2
RHS(n) = RHS(n) - BoundaryCondition.BCT;
else
L_ts = L(n);
SH = 0.1200 * Constants.c_a * (SoilVariables.T(n) - ForcingData.Ta_msr(KT)) / r_a_SOIL(KT);
RHS(n) = RHS(n) + 100 * Rn_SOIL(KT) / 1800 - Constants.RHOL * L_ts * EVAP(KT) - SH + Constants.RHOL * Constants.c_L * (ForcingData.Ta_msr(KT) * Precip(KT) + BoundaryCondition.DSTOR0 * SoilVariables.T(n) / Delt_t); % J cm-2 s-1
end
end
10 changes: 10 additions & 0 deletions src/+energy/calculateEnergyFluxes.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
function [QET, QEB] = calculateEnergyFluxes(SAVE, TT)
%{
Calculate the energy fluxes on the boundary nodes.
SarahAlidoost marked this conversation as resolved.
Show resolved Hide resolved
%}

ModelSettings = io.getModelSettings();
n = ModelSettings.NN;
QET = SAVE(2, 1, 2) - SAVE(2, 2, 2) * TT(n - 1) - SAVE(2, 3, 2) * TT(n);
QEB = -SAVE(1, 1, 2) + SAVE(1, 2, 2) * TT(1) + SAVE(1, 3, 2) * TT(2);
end
154 changes: 154 additions & 0 deletions src/+energy/calculateEnergyParameters.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
function EnergyVariables = calculateEnergyParameters(InitialValues, SoilVariables, HeatVariables, TransportCoefficient, AirVariabes, ...
VaporVariables, GasDispersivity, ThermalConductivityCapacity, ...
DRHOVh, DRHOVT, KL_T, Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L)
%{
Calculate all the parameters related to energy balance equation e.Constants.g.,
Equation 3.65-3.73, STEMMUS Technical Notes, page 29-32.
%}

ModelSettings = io.getModelSettings();
Constants = io.define_constants();

% input
Kcva = InitialValues.Kcva;
Kcah = InitialValues.Kcah;
KcaT = InitialValues.KcaT;
Kcaa = InitialValues.Kcaa;
Ccah = InitialValues.Ccah;
CcaT = InitialValues.CcaT;
Ccaa = InitialValues.Ccaa;
CTT_PH = InitialValues.CTT_PH;
CTT_LT = InitialValues.CTT_LT;
CTT_Lg = InitialValues.CTT_Lg;
CTT_g = InitialValues.CTT_g;
KLhBAR = AirVariabes.KLhBAR;
KLTBAR = AirVariabes.KLTBAR;
DDhDZ = AirVariabes.DDhDZ;
DhDZ = AirVariabes.DhDZ;
DTDZ = AirVariabes.DTDZ;
Kaa = AirVariabes.Kaa;
Vaa = AirVariabes.Vaa;
QL = AirVariabes.QL;

% output
EnergyVariables.CTh = InitialValues.CTh;
EnergyVariables.CTa = InitialValues.CTa;
EnergyVariables.KTh = InitialValues.KTh;
EnergyVariables.KTT = InitialValues.KTT;
EnergyVariables.KTa = InitialValues.KTa;
EnergyVariables.VTT = InitialValues.VTT;
EnergyVariables.VTh = InitialValues.VTh;
EnergyVariables.VTa = InitialValues.VTa;
EnergyVariables.CTg = InitialValues.CTg;
EnergyVariables.CTT = InitialValues.CTT;

for i = 1:ModelSettings.NL
if ~ModelSettings.Soilairefc
KLhBAR(i) = (SoilVariables.KfL_h(i, 1) + SoilVariables.KfL_h(i, 2)) / 2;
KLTBAR(i) = (KL_T(i, 1) + KL_T(i, 2)) / 2;
DDhDZ(i) = (SoilVariables.hh(i + 1) - SoilVariables.hh(i)) / ModelSettings.DeltZ(i);
DhDZ(i) = (SoilVariables.hh(i + 1) + SoilVariables.hh_frez(i + 1) - SoilVariables.hh(i) - SoilVariables.hh_frez(i)) / ModelSettings.DeltZ(i);
DTDZ(i) = (SoilVariables.TT(i + 1) - SoilVariables.TT(i)) / ModelSettings.DeltZ(i);
end
DTDBAR(i) = (TransportCoefficient.D_Ta(i, 1) + TransportCoefficient.D_Ta(i, 2)) / 2;
DEhBAR = (VaporVariables.D_V(i, 1) + VaporVariables.D_V(i, 2)) / 2;
DRHOVhDz(i) = (DRHOVh(i + 1) + DRHOVh(i)) / 2;
DRHOVTDz(i) = (DRHOVT(i + 1) + DRHOVT(i)) / 2;
RHOVBAR = (RHOV(i + 1) + RHOV(i)) / 2;
EtaBAR = (VaporVariables.Eta(i, 1) + VaporVariables.Eta(i, 2)) / 2;

% The soil air gas in soil-pore is considered with Xah and XaT
% terms.(0.0003,volumetric heat capacity)
if ~ModelSettings.Soilairefc
QL(i) = -(KLhBAR(i) * DhDZ(i) + (KLTBAR(i) + DTDBAR(i)) * DTDZ(i) + KLhBAR(i));
Qa = 0;
else
Qa = -((DEhBAR + GasDispersivity.D_Vg(i)) * DRHODAz(i) - RHODA(i) * (GasDispersivity.V_A(i) + Constants.Hc * QL(i) / Constants.RHOL));
end

if SoilVariables.DVa_Switch == 1
QV = -(DEhBAR + GasDispersivity.D_Vg(i)) * DRHOVhDz(i) * DDhDZ(i) - (DEhBAR * EtaBAR + GasDispersivity.D_Vg(i)) * DRHOVTDz(i) * DTDZ(i) + RHOVBAR * GasDispersivity.V_A(i);
else
QV = -(DEhBAR + GasDispersivity.D_Vg(i)) * DRHOVhDz(i) * DDhDZ(i) - (DEhBAR * EtaBAR + GasDispersivity.D_Vg(i)) * DRHOVTDz(i) * DTDZ(i);
end

% These are unused vars, but I comment them for future reference,
% See issue 100, item 1
% DVH(i) = (DEhBAR) * DRHOVhDz(i);
% DVT(i) = (DEhBAR * EtaBAR) * DRHOVTDz(i);
% QVH(i) = -(DEhBAR + GasDispersivity.D_Vg(i)) * DRHOVhDz(i) * DDhDZ(i);
% QVT(i) = -(DEhBAR * EtaBAR + GasDispersivity.D_Vg(i)) * DRHOVTDz(i) * DTDZ(i);
for j = 1:ModelSettings.nD
MN = i + j - 1;
if ModelSettings.Soilairefc == 1
Kcah(i, j) = Constants.c_a * SoilVariables.TT(MN) * ((VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) * Xah(MN) + Constants.Hc * RHODA(MN) * SoilVariables.KfL_h(i, j));
KcaT(i, j) = Constants.c_a * SoilVariables.TT(MN) * ((VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) * XaT(MN) + Constants.Hc * RHODA(MN) * (KL_T(i, j) + TransportCoefficient.D_Ta(i, j))); %
Kcaa(i, j) = Constants.c_a * SoilVariables.TT(MN) * Kaa(i, j);
if SoilVariables.DVa_Switch == 1
Kcva(i, j) = L(MN) * RHOV(MN) * GasDispersivity.Beta_g(i, j);
else
Kcva(i, j) = 0;
end
Ccah(i, j) = Constants.c_a * SoilVariables.TT(MN) * (-GasDispersivity.V_A(i) - Constants.Hc * QL(i) / Constants.RHOL) * Xah(MN);
CcaT(i, j) = Constants.c_a * SoilVariables.TT(MN) * (-GasDispersivity.V_A(i) - Constants.Hc * QL(i) / Constants.RHOL) * XaT(MN);
Ccaa(i, j) = Constants.c_a * SoilVariables.TT(MN) * Vaa(i, j);
end

if abs(SoilVariables.SAVEDTheta_LLh(i, j) - SoilVariables.SAVEDTheta_UUh(i, j)) ~= 0
CTT_PH(i, j) = (10 * L_f^2 * Constants.RHOI / (Constants.g * (ModelSettings.T0 + SoilVariables.TT(MN)))) * SoilVariables.DTheta_UUh(i, j);
CTT_Lg(i, j) = (Constants.c_L * SoilVariables.TT(MN) + L(MN)) * SoilVariables.Theta_g(i, j) * DRHOVT(MN);
CTT_g(i, j) = Constants.c_a * SoilVariables.TT(MN) * SoilVariables.Theta_g(i, j) * XaT(MN);

CTT_LT(i, j) = (((Constants.c_L * SoilVariables.TT(MN) - TransportCoefficient.WW(i, j)) * Constants.RHOL - ((Constants.c_L * SoilVariables.TT(MN) + L(MN)) * RHOV(MN) + Constants.c_a * RHODA(MN) * SoilVariables.TT(MN))) * (1 - Constants.RHOI / Constants.RHOL) - Constants.RHOI * Constants.c_i * SoilVariables.TT(MN)) * 1e4 * L_f / (Constants.g * (ModelSettings.T0 + SoilVariables.TT(MN))) * SoilVariables.DTheta_UUh(i, j);
if CTT_PH(i, j) < 0
CTT_PH(i, j) = 0;
end
EnergyVariables.CTT(i, j) = ThermalConductivityCapacity.c_unsat(i, j) + CTT_Lg(i, j) + CTT_g(i, j) + CTT_LT(i, j) + CTT_PH(i, j);
EnergyVariables.CTh(i, j) = (Constants.c_L * SoilVariables.TT(MN) + L(MN)) * SoilVariables.Theta_g(i, j) * DRHOVh(MN) + Constants.c_a * SoilVariables.TT(MN) * SoilVariables.Theta_g(i, j) * Xah(MN);
EnergyVariables.CTa(i, j) = SoilVariables.TT(MN) * SoilVariables.Theta_V(i, j) * Constants.c_a * Xaa(MN); % This term isnot in Milly's work.

else
% Main coefficients for energy transport is here
CTT_Lg(i, j) = 0;
CTT_g(i, j) = 0;
CTT_LT(i, j) = 0;
CTT_PH(i, j) = 0;
EnergyVariables.CTh(i, j) = ((Constants.c_L * SoilVariables.TT(MN) - TransportCoefficient.WW(i, j)) * Constants.RHOL - (Constants.c_L * SoilVariables.TT(MN) + L(MN)) * RHOV(MN) - Constants.c_a * RHODA(MN) * SoilVariables.TT(MN)) * SoilVariables.DTheta_LLh(i, j);
+(Constants.c_L * SoilVariables.TT(MN) + L(MN)) * SoilVariables.Theta_g(i, j) * DRHOVh(MN) + Constants.c_a * SoilVariables.TT(MN) * SoilVariables.Theta_g(i, j) * Xah(MN);
EnergyVariables.CTT(i, j) = ThermalConductivityCapacity.c_unsat(i, j) + (Constants.c_L * SoilVariables.TT(MN) + L(MN)) * SoilVariables.Theta_g(i, j) * DRHOVT(MN) + Constants.c_a * SoilVariables.TT(MN) * SoilVariables.Theta_g(i, j) * XaT(MN) ...
+ ((Constants.c_L * SoilVariables.TT(MN) - TransportCoefficient.WW(i, j)) * Constants.RHOL - (Constants.c_L * SoilVariables.TT(MN) + L(MN)) * RHOV(MN) - Constants.c_a * RHODA(MN) * SoilVariables.TT(MN)) * SoilVariables.DTheta_LLT(i, j);
EnergyVariables.CTa(i, j) = SoilVariables.TT(MN) * SoilVariables.Theta_V(i, j) * Constants.c_a * Xaa(MN); % This term isnot in Milly's work.
end
if ModelSettings.SFCC == 0 % ice calculation use Sin function
if SoilVariables.TT(MN) + 273.15 > Tf1
CTT_PH(i, j) = 0;
elseif SoilVariables.TT(MN) + 273.15 >= Tf2
CTT_PH(i, j) = L_f * 10^(-3) * 0.5 * cos(pi() * (SoilVariables.TT(MN) + 273.15 - 0.5 * Tf1 - 0.5 * Tf2) / (Tf1 - Tf2)) * pi() / (Tf1 - Tf2);
else
CTT_PH(i, j) = 0;
end
CTT_Lg(i, j) = (Constants.c_L * SoilVariables.TT(MN) + L(MN)) * SoilVariables.Theta_g(i, j) * DRHOVT(MN);
CTT_g(i, j) = Constants.c_a * SoilVariables.TT(MN) * SoilVariables.Theta_g(i, j) * XaT(MN);
CTT_LT(i, j) = ((Constants.c_L * SoilVariables.TT(MN) - Constants.c_i * SoilVariables.TT(MN) - TransportCoefficient.WW(i, j)) * Constants.RHOL + ((Constants.c_L * SoilVariables.TT(MN) + L(MN)) * RHOV(MN) + Constants.c_a * RHODA(MN) * SoilVariables.TT(MN)) * (Constants.RHOL / Constants.RHOI - 1)) * 1e4 * L_f / (Constants.g * (ModelSettings.T0 + SoilVariables.TT(MN))) * SoilVariables.DTheta_UUh(i, j);

EnergyVariables.CTT(i, j) = ThermalConductivityCapacity.c_unsat(i, j) + CTT_Lg(i, j) + CTT_g(i, j) + CTT_LT(i, j) + CTT_PH(i, j);
EnergyVariables.CTh(i, j) = (Constants.c_L * SoilVariables.TT(MN) + L(MN)) * SoilVariables.Theta_g(i, j) * DRHOVh(MN) + Constants.c_a * SoilVariables.TT(MN) * SoilVariables.Theta_g(i, j) * Xah(MN);
EnergyVariables.CTa(i, j) = SoilVariables.TT(MN) * SoilVariables.Theta_V(i, j) * Constants.c_a * Xaa(MN); % This term isnot in Milly's work.
end
EnergyVariables.KTh(i, j) = L(MN) * (VaporVariables.D_V(i, j) + GasDispersivity.D_Vg(i)) * DRHOVh(MN) + Constants.c_L * SoilVariables.TT(MN) * Constants.RHOL * HeatVariables.Khh(i, j) + Kcah(i, j);
EnergyVariables.KTT(i, j) = ThermalConductivityCapacity.Lambda_eff(i, j) + Constants.c_L * SoilVariables.TT(MN) * Constants.RHOL * HeatVariables.KhT(i, j) + KcaT(i, j) + L(MN) * (VaporVariables.D_V(i, j) * VaporVariables.Eta(i, j) + GasDispersivity.D_Vg(i)) * DRHOVT(MN);
EnergyVariables.KTa(i, j) = Kcva(i, j) + Kcaa(i, j) + Constants.c_L * SoilVariables.TT(MN) * Constants.RHOL * HeatVariables.Kha(i, j); % This term isnot in Milly's work.

if SoilVariables.DVa_Switch == 1
EnergyVariables.VTh(i, j) = Constants.c_L * SoilVariables.TT(MN) * Constants.RHOL * HeatVariables.Vvh(i, j) + Ccah(i, j) - L(MN) * GasDispersivity.V_A(i) * DRHOVh(MN);
EnergyVariables.VTT(i, j) = Constants.c_L * SoilVariables.TT(MN) * Constants.RHOL * HeatVariables.VvT(i, j) + CcaT(i, j) - L(MN) * GasDispersivity.V_A(i) * DRHOVT(MN) - (Constants.c_L * (QL(i) + QV) + Constants.c_a * Qa - 2.369 * QV);
else
EnergyVariables.VTh(i, j) = Constants.c_L * SoilVariables.TT(MN) * Constants.RHOL * HeatVariables.Vvh(i, j) + Ccah(i, j);
EnergyVariables.VTT(i, j) = Constants.c_L * SoilVariables.TT(MN) * Constants.RHOL * HeatVariables.VvT(i, j) + CcaT(i, j) - (Constants.c_L * (QL(i) + QV) + Constants.c_a * Qa - 2.369 * QV);
end

EnergyVariables.VTa(i, j) = Ccaa(i, j);
EnergyVariables.CTg(i, j) = (Constants.c_L * Constants.RHOL + Constants.c_a * Constants.Hc * RHODA(MN)) * SoilVariables.KfL_h(i, j) * SoilVariables.TT(MN) - Constants.c_L * Srt(i, j) * SoilVariables.TT(MN);
end
end
end
Loading