Skip to content

Commit fe5fec3

Browse files
MostafaGomaa93MostafaGomaa93
authored
Merge pull request #282 from EcoExtreML/correct_precip_snow_backup
Correct Snow precipitation
2 parents 8d3ef75 + 47e2ed3 commit fe5fec3

16 files changed

Lines changed: 175 additions & 134 deletions

‎.gitignore‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1 +1,7 @@
11
.data/
2+
run_model_on_snellius/exe/includedSupportPackages.txt
3+
run_model_on_snellius/exe/mccExcludedFiles.log
4+
run_model_on_snellius/exe/readme.txt
5+
run_model_on_snellius/exe/requiredMCRProducts.txt
6+
run_model_on_snellius/exe/run_STEMMUS_SCOPE.sh
7+
run_model_on_snellius/exe/unresolvedSymbols.txt

‎CHANGELOG.md‎

Lines changed: 11 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,20 @@
11
# Unreleased
22

3-
**Added:**
3+
**Changed:**
4+
- Change array shape of precipitation and runoff fluxes from time-series array to 1D array (for BMI purposes) discussed in [#280](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/280) and added in [#282](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/282)
5+
6+
**Fixed:**
7+
- Correct snow precipitation calcuations discussed in [#279](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/279) and fixed in [#282](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/282)
8+
9+
<a name="1.6.2"></a>
10+
# [1.6.2](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.6.2) - 17 Dec 2024
411

12+
The 1.6.2 release comes with minor changes as:
13+
14+
**Added:**
515
- Documentation using mkdocs and a github action workflow to publish the
616
documentation in [#264](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/264)
717

8-
**Changed:**
9-
1018
<a name="1.6.1"></a>
1119
# [1.6.1](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.6.1) - 26 Sep 2024
1220

‎run_model_on_snellius/exe/STEMMUS_SCOPE‎

-671 Bytes
Binary file not shown.

‎src/+energy/calculateBoundaryConditions.m‎

Lines changed: 23 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,27 @@
11
function [RHS, EnergyMatrices] = calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, SoilVariables, ...
2-
Precip, Evap, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings)
2+
Evap, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings)
33
%{
4-
Determine the boundary condition for solving the energy equation, see
5-
STEMMUS Technical Notes.
4+
Determine the boundary condition for solving the energy equation, see
5+
STEMMUS Technical Notes.
6+
Inputs:
7+
BoundaryCondition : structure
8+
EnergyMatrices : structure
9+
HBoundaryFlux : structure
10+
ForcingData : structure
11+
SoilVariables : structure
12+
ModelSettings : structure
13+
GroundwaterSettings: structure
14+
15+
Evap : [g cm^-2 s^-1] Evaporation
16+
Delt_t : [s] Time step, here it is 1800 s
17+
r_a_SOIL: [s cm^-1] Soil surface aerodynamic resistance
18+
Rn_SOIL : [] Net radiation reaching the soil surface
19+
L : [] Latent heat of vaporization
20+
KT : [-] Number of time step
21+
Outputs:
22+
RHS : [degree C] Soil temperature
23+
24+
EnergyMatrices : structure
625
%}
726

827
n = ModelSettings.NN;
@@ -46,6 +65,6 @@
4665
else
4766
L_ts = L(n);
4867
SH = 0.1200 * Constants.c_a * (SoilVariables.T(n) - ForcingData.Ta_msr(KT)) / r_a_SOIL(KT);
49-
RHS(n) = RHS(n) + 100 * Rn_SOIL(KT) / 1800 - Constants.RHOL * L_ts * Evap - SH + Constants.RHOL * Constants.c_L * (ForcingData.Ta_msr(KT) * Precip(KT) + BoundaryCondition.DSTOR0 * SoilVariables.T(n) / Delt_t); % J cm-2 s-1
68+
RHS(n) = RHS(n) + 100 * Rn_SOIL(KT) / 1800 - Constants.RHOL * L_ts * Evap - SH + Constants.RHOL * Constants.c_L * (ForcingData.Ta_msr(KT) * ForcingData.effectivePrecip + BoundaryCondition.DSTOR0 * SoilVariables.T(n) / Delt_t); % J cm-2 s-1
5069
end
5170
end

‎src/+energy/calculateEnergyFluxes.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function [QET, QEB] = calculateEnergyFluxes(SAVE, TT, ModelSettings)
1+
function [QET, QEB] = calculateEnergyFluxes(SAVE, TT, ModelSettings, GroundwaterSettings)
22
%{
33
Calculate the energy fluxes on the boundary nodes, see STEMMUS Technical
44
Notes.

‎src/+energy/solveEnergyBalanceEquations.m‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
1-
function [RHS, SAVE, CHK, SoilVariables, EnergyVariables] = solveEnergyBalanceEquations(InitialValues, SoilVariables, HeatVariables, TransportCoefficient, ...
2-
AirVariabes, VaporVariables, GasDispersivity, ThermalConductivityCapacity, ...
3-
HBoundaryFlux, BoundaryCondition, ForcingData, DRHOVh, DRHOVT, KL_T, ...
4-
Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, Delt_t, P_g, P_gg, ...
5-
TOLD, Precip, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings)
1+
function [RHS, SAVE, CHK, SoilVariables, EnergyVariables, gwfluxes] = solveEnergyBalanceEquations(InitialValues, SoilVariables, HeatVariables, TransportCoefficient, ...
2+
AirVariabes, VaporVariables, GasDispersivity, ThermalConductivityCapacity, ...
3+
HBoundaryFlux, BoundaryCondition, ForcingData, DRHOVh, DRHOVT, KL_T, ...
4+
Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, Delt_t, P_g, P_gg, ...
5+
TOLD, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings)
66
%{
77
Solve the Energy balance equation with the Thomas algorithm to update
88
the soil temperature 'SoilVariables.TT', the finite difference
@@ -19,7 +19,7 @@
1919
[RHS, EnergyMatrices, SAVE] = energy.assembleCoefficientMatrices(InitialValues, EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, ModelSettings, GroundwaterSettings);
2020

2121
[RHS, EnergyMatrices] = energy.calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, SoilVariables, ...
22-
Precip, EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings);
22+
EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings);
2323

2424
[SoilVariables, CHK, RHS, EnergyMatrices] = energy.solveTridiagonalMatrixEquations(EnergyMatrices, SoilVariables, RHS, CHK, ModelSettings, GroundwaterSettings);
2525

@@ -42,7 +42,7 @@
4242
end
4343
end
4444

45-
% These are unused vars, but I comment them for future reference,
46-
% See issue 100, item 2
47-
% [QET, QEB] = energy.calculateEnergyFluxes(SAVE, TT)(SAVE, SoilVariables.TT);
45+
[QET, QEB] = energy.calculateEnergyFluxes(SAVE, SoilVariables.TT, ModelSettings, GroundwaterSettings);
46+
gwfluxes.energyTopflux = QET; % energy flux at the top boundary node
47+
gwfluxes.energyBotmflux = QEB; % energy flux at the bottom boundary node
4848
end

‎src/+groundwater/calculateGroundwaterRecharge.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function gwfluxes = calculateGroundwaterRecharge(EnergyVariables, SoilVariables, KT, ModelSettings, GroundwaterSettings)
1+
function gwfluxes = calculateGroundwaterRecharge(EnergyVariables, SoilVariables, KT, gwfluxes, ModelSettings, GroundwaterSettings)
22
%{
33
Added by Mostafa, modified after Lianyu
44
The concept followed to calculate groundwater recharge can be found in:

‎src/+groundwater/updateDunnianRunoff.m‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,13 @@
1-
function [dunnianRunoff, update_Precip_msr] = updateDunnianRunoff(Precip_msr, groundWaterDepth)
1+
function [runoffDunnian, update_Precip_msr] = updateDunnianRunoff(Precip_msr, groundWaterDepth)
22
% Dunnian runoff = Direct water input from precipitation + return flow
33
% (a) direct water input from precipitation when soil is fully saturated (depth to water table = 0)
44
% (b) Return flow (from groundwater exfiltration) calculated in MODFLOW and added to Dunnian runoff (through BMI)
55
% here approach (a) is implemented
6-
dunnianRunoff = zeros(size(Precip_msr));
6+
runoffDunnian = zeros(size(Precip_msr));
77
update_Precip_msr = Precip_msr;
88

99
if groundWaterDepth <= 1.0
10-
dunnianRunoff = Precip_msr;
10+
runoffDunnian = Precip_msr;
1111
update_Precip_msr = zeros(size(Precip_msr));
1212
end
1313

‎src/+init/defineInitialValues.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -171,7 +171,7 @@
171171

172172
Nmsrmn = Dur_tot * 10; % Here, it is made as big as possible, in case a long simulation period containing many time step is defined.
173173
fields = {
174-
'Precip', 'Ta', 'Ts', 'U', 'HR_a', 'Rns', 'Rnl', 'Rn', ...
174+
'Ta', 'Ts', 'U', 'HR_a', 'Rns', 'Rnl', 'Rn', ...
175175
'h_SUR', 'SH', 'MO', 'Zeta_MO', 'TopPg'
176176
};
177177
structures{4} = helpers.createStructure(zeros(Nmsrmn, 1), fields);

‎src/+io/loadForcingData.m‎

Lines changed: 5 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -19,25 +19,16 @@
1919
ForcingData.G_msr = Mdata{:, 7} * 0.15;
2020
Precip_msr = Mdata{:, 6}; % (cm/sec)
2121

22-
%%%%%%%%%% Adjust precipitation to get applied infiltration after removing: (1) saturation excess runoff %%%%%%%%%%
23-
%%%%%%%%%% (2) infiltration excess runoff %%%%%%%%%%
24-
% Note: Code changes are related to the issue: https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/232
25-
% Note: Adjusting the precipitation after the canopy interception is not implemented yet.
26-
27-
% (1) Saturation excess runoff (Dunnian runoff)
22+
% Calculate saturation excess runoff (Dunnian runoff)
2823
if ~GroundwaterSettings.GroundwaterCoupling % Groundwater Coupling is not activated
29-
% Concept is adopted from the CLM model (see issue 232 in GitHub for more explanation)
24+
% Concept is adopted from the CLM model (see issue 232, https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/232)
3025
% Check also the CLM documents (https://doi.org/10.5065/D6N877R0, https://doi.org/10.1029/2005JD006111)
31-
wat_Dep = Tot_Depth / 100; % (m)
26+
wat_Dep = Tot_Depth / 100; % (m), this assumption water depth = total soil depth is not fully correct (to be improved)
3227
fover = 0.5; % decay factor (fixed to 0.5 m-1)
3328
fmax = SoilProperties.fmax; % potential maximum value of fsat
3429
fsat = (fmax .* exp(-0.5 * fover * wat_Dep)); % fraction of saturated area (unitless)
35-
ForcingData.R_Dunn = Precip_msr .* fsat; % Dunnian runoff (saturation excess runoff, in c/sec)
36-
Precip_msr = Precip_msr .* (1 - fsat); % applied infiltration after removing Dunnian runoff
37-
end
38-
39-
% (2) Infiltration excess runoff (Hortonian runoff)
40-
ForcingData.R_Hort = zeros(size(Precip_msr)); % will be updated in +soilmoisture/calculateBoundaryConditions
30+
ForcingData.runoffDunnian = Precip_msr .* fsat; % Dunnian runoff (saturation excess runoff, in cm/sec)
31+
end % In case Groundwater Coupling is activated, Dunnian runoff is calculated in +groundwater/updateDunnianRunoff
4132

4233
% replace negative values
4334
for jj = 1:Dur_tot
@@ -49,7 +40,4 @@
4940
% Outputs to be used by other functions
5041
ForcingData.Tmin = min(ForcingData.Ta_msr);
5142
ForcingData.Precip_msr = Precip_msr;
52-
% Applied infiltration (= precipitation after removal of Dunnian runoff)
53-
ForcingData.applied_inf = Precip_msr; % later will be updated in the ....
54-
% +soilmoisture/calculateBoundaryConditions after removal of Hortonian runoff
5543
end

0 commit comments

Comments
 (0)