-
Notifications
You must be signed in to change notification settings - Fork 6
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #197 from EcoExtreML/fix_issue_100
Fix issue 100
- Loading branch information
Showing
18 changed files
with
484 additions
and
522 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,44 @@ | ||
| 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, see | ||
| STEMMUS Technical Notes. | ||
| %} | ||
|
|
||
| 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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,11 @@ | ||
| function [QET, QEB] = calculateEnergyFluxes(SAVE, TT) | ||
| %{ | ||
| Calculate the energy fluxes on the boundary nodes, see STEMMUS Technical | ||
| Notes. | ||
| %} | ||
|
|
||
| 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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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 |
Oops, something went wrong.