Rev 680 | Rev 720 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 680 | Rev 697 | ||
---|---|---|---|
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/watercycle.cpp': |
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/watercycle.cpp': |
2 | /********************************************************************************************
|
2 | /********************************************************************************************
|
3 | ** iLand - an individual based forest landscape and disturbance model
|
3 | ** iLand - an individual based forest landscape and disturbance model
|
4 | ** http://iland.boku.ac.at
|
4 | ** http://iland.boku.ac.at
|
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl
|
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl
|
6 | **
|
6 | **
|
7 | ** This program is free software: you can redistribute it and/or modify
|
7 | ** This program is free software: you can redistribute it and/or modify
|
8 | ** it under the terms of the GNU General Public License as published by
|
8 | ** it under the terms of the GNU General Public License as published by
|
9 | ** the Free Software Foundation, either version 3 of the License, or
|
9 | ** the Free Software Foundation, either version 3 of the License, or
|
10 | ** (at your option) any later version.
|
10 | ** (at your option) any later version.
|
11 | **
|
11 | **
|
12 | ** This program is distributed in the hope that it will be useful,
|
12 | ** This program is distributed in the hope that it will be useful,
|
13 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of
|
13 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of
|
14 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
14 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
15 | ** GNU General Public License for more details.
|
15 | ** GNU General Public License for more details.
|
16 | **
|
16 | **
|
17 | ** You should have received a copy of the GNU General Public License
|
17 | ** You should have received a copy of the GNU General Public License
|
18 | ** along with this program. If not, see <http://www.gnu.org/licenses/>.
|
18 | ** along with this program. If not, see <http://www.gnu.org/licenses/>.
|
19 | ********************************************************************************************/
|
19 | ********************************************************************************************/
|
20 | 20 | ||
21 | #include "global.h"
|
21 | #include "global.h"
|
22 | #include "watercycle.h"
|
22 | #include "watercycle.h"
|
23 | #include "climate.h"
|
23 | #include "climate.h"
|
24 | #include "resourceunit.h"
|
24 | #include "resourceunit.h"
|
25 | #include "species.h"
|
25 | #include "species.h"
|
26 | #include "model.h"
|
26 | #include "model.h"
|
27 | #include "helper.h"
|
27 | #include "helper.h"
|
28 | #include "modules.h"
|
28 | #include "modules.h"
|
- | 29 | ||
- | 30 | /** @class WaterCycle
|
|
- | 31 | @ingroup core
|
|
- | 32 | simulates the water cycle on a ResourceUnit.
|
|
- | 33 | The WaterCycle is simulated with a daily time step on the spatial level of a ResourceUnit. Related are
|
|
- | 34 | the snow module (SnowPack), and Canopy module that simulates the interception (and evaporation) of precipitation and the
|
|
- | 35 | transpiration from the canopy.
|
|
- | 36 | The WaterCycle covers the "soil water bucket". Main entry function is run().
|
|
- | 37 | ||
- | 38 | See http://iland.boku.ac.at/water+cycle
|
|
- | 39 | */
|
|
29 | 40 | ||
30 | WaterCycle::WaterCycle() |
41 | WaterCycle::WaterCycle() |
31 | {
|
42 | {
|
32 | mSoilDepth = 0; |
43 | mSoilDepth = 0; |
33 | mLastYear = -1; |
44 | mLastYear = -1; |
34 | }
|
45 | }
|
35 | 46 | ||
36 | void WaterCycle::setup(const ResourceUnit *ru) |
47 | void WaterCycle::setup(const ResourceUnit *ru) |
37 | {
|
48 | {
|
38 | mRU = ru; |
49 | mRU = ru; |
39 | // get values...
|
50 | // get values...
|
40 | mFieldCapacity = 0.; // on top |
51 | mFieldCapacity = 0.; // on top |
41 | const XmlHelper &xml=GlobalSettings::instance()->settings(); |
52 | const XmlHelper &xml=GlobalSettings::instance()->settings(); |
42 | mSoilDepth = xml.valueDouble("model.site.soilDepth", 0.) * 10; // convert from cm to mm |
53 | mSoilDepth = xml.valueDouble("model.site.soilDepth", 0.) * 10; // convert from cm to mm |
43 | double pct_sand = xml.valueDouble("model.site.pctSand"); |
54 | double pct_sand = xml.valueDouble("model.site.pctSand"); |
44 | double pct_silt = xml.valueDouble("model.site.pctSilt"); |
55 | double pct_silt = xml.valueDouble("model.site.pctSilt"); |
45 | double pct_clay = xml.valueDouble("model.site.pctClay"); |
56 | double pct_clay = xml.valueDouble("model.site.pctClay"); |
46 | if (fabs(100. - (pct_sand + pct_silt + pct_clay)) > 0.01) |
57 | if (fabs(100. - (pct_sand + pct_silt + pct_clay)) > 0.01) |
47 | throw IException(QString("Setup Watercycle: soil composition percentages do not sum up to 100. Sand: %1, Silt: %2 Clay: %3").arg(pct_sand).arg(pct_silt).arg(pct_clay)); |
58 | throw IException(QString("Setup Watercycle: soil composition percentages do not sum up to 100. Sand: %1, Silt: %2 Clay: %3").arg(pct_sand).arg(pct_silt).arg(pct_clay)); |
48 | 59 | ||
49 | // calculate soil characteristics based on empirical functions (Schwalm & Ek, 2004)
|
60 | // calculate soil characteristics based on empirical functions (Schwalm & Ek, 2004)
|
50 | // note: the variables are percentages [0..100]
|
61 | // note: the variables are percentages [0..100]
|
51 | mPsi_ref = -exp((1.54 - 0.0095*pct_sand + 0.0063*pct_silt) * log(10)) * 0.000098; // Eq. 83 |
62 | mPsi_ref = -exp((1.54 - 0.0095*pct_sand + 0.0063*pct_silt) * log(10)) * 0.000098; // Eq. 83 |
52 | mPsi_koeff_b = -( 3.1 + 0.157*pct_clay - 0.003*pct_sand ); // Eq. 84 |
63 | mPsi_koeff_b = -( 3.1 + 0.157*pct_clay - 0.003*pct_sand ); // Eq. 84 |
53 | mRho_ref = 0.01 * (50.5 - 0.142*pct_sand - 0.037*pct_clay); // Eq. 78 |
64 | mRho_ref = 0.01 * (50.5 - 0.142*pct_sand - 0.037*pct_clay); // Eq. 78 |
54 | mCanopy.setup(); |
65 | mCanopy.setup(); |
55 | 66 | ||
56 | mPermanentWiltingPoint = heightFromPsi(-4000); // maximum psi is set to a constant of -4MPa |
67 | mPermanentWiltingPoint = heightFromPsi(-4000); // maximum psi is set to a constant of -4MPa |
57 | if (xml.valueBool("model.settings.waterUseSoilSaturation",false)==false) { |
68 | if (xml.valueBool("model.settings.waterUseSoilSaturation",false)==false) { |
58 | mFieldCapacity = heightFromPsi(-15); |
69 | mFieldCapacity = heightFromPsi(-15); |
59 | } else { |
70 | } else { |
60 | // =-EXP((1.54-0.0095* pctSand +0.0063* pctSilt)*LN(10))*0.000098
|
71 | // =-EXP((1.54-0.0095* pctSand +0.0063* pctSilt)*LN(10))*0.000098
|
61 | double psi_sat = -exp((1.54-0.0095 * pct_sand + 0.0063*pct_silt)*log(10.))*0.000098; |
72 | double psi_sat = -exp((1.54-0.0095 * pct_sand + 0.0063*pct_silt)*log(10.))*0.000098; |
62 | mFieldCapacity = heightFromPsi(psi_sat); |
73 | mFieldCapacity = heightFromPsi(psi_sat); |
63 | if (logLevelDebug()) qDebug() << "psi: saturation " << psi_sat << "field capacity:" << mFieldCapacity; |
74 | if (logLevelDebug()) qDebug() << "psi: saturation " << psi_sat << "field capacity:" << mFieldCapacity; |
64 | }
|
75 | }
|
65 | 76 | ||
66 | mContent = mFieldCapacity; // start with full water content (in the middle of winter) |
77 | mContent = mFieldCapacity; // start with full water content (in the middle of winter) |
67 | if (logLevelDebug()) qDebug() << "setup of water: Psi_ref (kPa)" << mPsi_ref << "Rho_ref" << mRho_ref << "coeff. b" << mPsi_koeff_b; |
78 | if (logLevelDebug()) qDebug() << "setup of water: Psi_ref (kPa)" << mPsi_ref << "Rho_ref" << mRho_ref << "coeff. b" << mPsi_koeff_b; |
68 | mCanopyConductance = 0.; |
79 | mCanopyConductance = 0.; |
69 | mLastYear = -1; |
80 | mLastYear = -1; |
70 | 81 | ||
71 | // canopy settings
|
82 | // canopy settings
|
72 | mCanopy.mNeedleFactor = xml.valueDouble("model.settings.interceptionStorageNeedle", 4.); |
83 | mCanopy.mNeedleFactor = xml.valueDouble("model.settings.interceptionStorageNeedle", 4.); |
73 | mCanopy.mDecidousFactor = xml.valueDouble("model.settings.interceptionStorageBroadleaf", 2.); |
84 | mCanopy.mDecidousFactor = xml.valueDouble("model.settings.interceptionStorageBroadleaf", 2.); |
74 | mSnowPack.mSnowTemperature = xml.valueDouble("model.settings.snowMeltTemperature", 0.); |
85 | mSnowPack.mSnowTemperature = xml.valueDouble("model.settings.snowMeltTemperature", 0.); |
75 | }
|
86 | }
|
76 | 87 | ||
77 | /** function to calculate the water pressure [saugspannung] for a given amount of water.
|
88 | /** function to calculate the water pressure [saugspannung] for a given amount of water.
|
78 | returns water potential in kPa.
|
89 | returns water potential in kPa.
|
79 | see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance */
|
90 | see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance */
|
80 | inline double WaterCycle::psiFromHeight(const double mm) const |
91 | inline double WaterCycle::psiFromHeight(const double mm) const |
81 | {
|
92 | {
|
82 | // psi_x = psi_ref * ( rho_x / rho_ref) ^ b
|
93 | // psi_x = psi_ref * ( rho_x / rho_ref) ^ b
|
83 | if (mm<0.001) |
94 | if (mm<0.001) |
84 | return -100000000; |
95 | return -100000000; |
85 | double psi_x = mPsi_ref * pow((mm / mSoilDepth / mRho_ref),mPsi_koeff_b); |
96 | double psi_x = mPsi_ref * pow((mm / mSoilDepth / mRho_ref),mPsi_koeff_b); |
86 | return psi_x; // pis |
97 | return psi_x; // pis |
87 | }
|
98 | }
|
88 | 99 | ||
89 | /// calculate the height of the water column for a given pressure
|
100 | /// calculate the height of the water column for a given pressure
|
90 | /// return water amount in mm
|
101 | /// return water amount in mm
|
91 | /// see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance
|
102 | /// see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance
|
92 | inline double WaterCycle::heightFromPsi(const double psi_kpa) const |
103 | inline double WaterCycle::heightFromPsi(const double psi_kpa) const |
93 | {
|
104 | {
|
94 | // rho_x = rho_ref * (psi_x / psi_ref)^(1/b)
|
105 | // rho_x = rho_ref * (psi_x / psi_ref)^(1/b)
|
95 | double h = mSoilDepth * mRho_ref * pow(psi_kpa / mPsi_ref, 1./mPsi_koeff_b); |
106 | double h = mSoilDepth * mRho_ref * pow(psi_kpa / mPsi_ref, 1./mPsi_koeff_b); |
96 | return h; |
107 | return h; |
97 | }
|
108 | }
|
98 | 109 | ||
99 | /// get canopy characteristics of the resource unit.
|
110 | /// get canopy characteristics of the resource unit.
|
100 | /// It is important, that species-statistics are valid when this function is called (LAI)!
|
111 | /// It is important, that species-statistics are valid when this function is called (LAI)!
|
101 | void WaterCycle::getStandValues() |
112 | void WaterCycle::getStandValues() |
102 | {
|
113 | {
|
103 | mLAINeedle=mLAIBroadleaved=0.; |
114 | mLAINeedle=mLAIBroadleaved=0.; |
104 | mCanopyConductance=0.; |
115 | mCanopyConductance=0.; |
105 | const double ground_vegetationCC = 0.02; |
116 | const double ground_vegetationCC = 0.02; |
106 | double lai; |
117 | double lai; |
107 | foreach(ResourceUnitSpecies *rus, mRU->ruSpecies()) { |
118 | foreach(ResourceUnitSpecies *rus, mRU->ruSpecies()) { |
108 | lai = rus->constStatistics().leafAreaIndex(); |
119 | lai = rus->constStatistics().leafAreaIndex(); |
109 | if (rus->species()->isConiferous()) |
120 | if (rus->species()->isConiferous()) |
110 | mLAINeedle+=lai; |
121 | mLAINeedle+=lai; |
111 | else
|
122 | else
|
112 | mLAIBroadleaved+=lai; |
123 | mLAIBroadleaved+=lai; |
113 | mCanopyConductance += rus->species()->canopyConductance() * lai; // weigh with LAI |
124 | mCanopyConductance += rus->species()->canopyConductance() * lai; // weigh with LAI |
114 | }
|
125 | }
|
115 | double total_lai = mLAIBroadleaved+mLAINeedle; |
126 | double total_lai = mLAIBroadleaved+mLAINeedle; |
116 | 127 | ||
117 | // handle cases with LAI < 1 (use generic "ground cover characteristics" instead)
|
128 | // handle cases with LAI < 1 (use generic "ground cover characteristics" instead)
|
118 | if (total_lai<1.) { |
129 | if (total_lai<1.) { |
119 | mCanopyConductance+=(ground_vegetationCC)*(1. - total_lai); |
130 | mCanopyConductance+=(ground_vegetationCC)*(1. - total_lai); |
120 | total_lai = 1.; |
131 | total_lai = 1.; |
121 | }
|
132 | }
|
122 | mCanopyConductance /= total_lai; |
133 | mCanopyConductance /= total_lai; |
123 | 134 | ||
124 | if (total_lai < Model::settings().laiThresholdForClosedStands) { |
135 | if (total_lai < Model::settings().laiThresholdForClosedStands) { |
125 | // following Landsberg and Waring: when LAI is < 3 (default for laiThresholdForClosedStands), a linear "ramp" from 0 to 3 is assumed
|
136 | // following Landsberg and Waring: when LAI is < 3 (default for laiThresholdForClosedStands), a linear "ramp" from 0 to 3 is assumed
|
126 | // http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance
|
137 | // http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance
|
127 | mCanopyConductance *= total_lai / Model::settings().laiThresholdForClosedStands; |
138 | mCanopyConductance *= total_lai / Model::settings().laiThresholdForClosedStands; |
128 | }
|
139 | }
|
129 | if (logLevelInfo()) qDebug() << "WaterCycle:getStandValues: LAI needle" << mLAINeedle << "LAI Broadl:"<< mLAIBroadleaved << "weighted avg. Conductance (m/2):" << mCanopyConductance; |
140 | if (logLevelInfo()) qDebug() << "WaterCycle:getStandValues: LAI needle" << mLAINeedle << "LAI Broadl:"<< mLAIBroadleaved << "weighted avg. Conductance (m/2):" << mCanopyConductance; |
130 | }
|
141 | }
|
131 | 142 | ||
132 | /// calculate responses for ground vegetation, i.e. for "unstocked" areas.
|
143 | /// calculate responses for ground vegetation, i.e. for "unstocked" areas.
|
133 | /// this duplicates calculations done in Species.
|
144 | /// this duplicates calculations done in Species.
|
134 | /// @return Minimum of vpd and soilwater response for default
|
145 | /// @return Minimum of vpd and soilwater response for default
|
135 | inline double WaterCycle::calculateBaseSoilAtmosphereResponse(const double psi_kpa, const double vpd_kpa) |
146 | inline double WaterCycle::calculateBaseSoilAtmosphereResponse(const double psi_kpa, const double vpd_kpa) |
136 | {
|
147 | {
|
137 | // constant parameters used for ground vegetation:
|
148 | // constant parameters used for ground vegetation:
|
138 | const double mPsiMin = 1.5; // MPa |
149 | const double mPsiMin = 1.5; // MPa |
139 | const double mRespVpdExponent = -0.6; |
150 | const double mRespVpdExponent = -0.6; |
140 | // see SpeciesResponse::soilAtmosphereResponses()
|
151 | // see SpeciesResponse::soilAtmosphereResponses()
|
141 | double water_resp; |
152 | double water_resp; |
142 | // see Species::soilwaterResponse:
|
153 | // see Species::soilwaterResponse:
|
143 | const double psi_mpa = psi_kpa / 1000.; // convert to MPa |
154 | const double psi_mpa = psi_kpa / 1000.; // convert to MPa |
144 | water_resp = limit( 1. - psi_mpa / mPsiMin, 0., 1.); |
155 | water_resp = limit( 1. - psi_mpa / mPsiMin, 0., 1.); |
145 | // see species::vpdResponse
|
156 | // see species::vpdResponse
|
146 | 157 | ||
147 | double vpd_resp; |
158 | double vpd_resp; |
148 | vpd_resp = exp(mRespVpdExponent * vpd_kpa); |
159 | vpd_resp = exp(mRespVpdExponent * vpd_kpa); |
149 | return qMin(water_resp, vpd_resp); |
160 | return qMin(water_resp, vpd_resp); |
150 | }
|
161 | }
|
151 | 162 | ||
152 | /// calculate combined VPD and soilwaterresponse for all species
|
163 | /// calculate combined VPD and soilwaterresponse for all species
|
153 | /// on the RU. This is used for the calc. of the transpiration.
|
164 | /// on the RU. This is used for the calc. of the transpiration.
|
154 | inline double WaterCycle::calculateSoilAtmosphereResponse(const double psi_kpa, const double vpd_kpa) |
165 | inline double WaterCycle::calculateSoilAtmosphereResponse(const double psi_kpa, const double vpd_kpa) |
155 | {
|
166 | {
|
156 | double min_response; |
167 | double min_response; |
157 | double total_response = 0; // LAI weighted minimum response for all speices on the RU |
168 | double total_response = 0; // LAI weighted minimum response for all speices on the RU |
158 | double total_lai_factor = 0.; |
169 | double total_lai_factor = 0.; |
159 | foreach(ResourceUnitSpecies *rus, mRU->ruSpecies()) { |
170 | foreach(ResourceUnitSpecies *rus, mRU->ruSpecies()) { |
160 | if (rus->LAIfactor()>0) { |
171 | if (rus->LAIfactor()>0) { |
161 | // retrieve the minimum of VPD / soil water response for that species
|
172 | // retrieve the minimum of VPD / soil water response for that species
|
162 | rus->speciesResponse()->soilAtmosphereResponses(psi_kpa, vpd_kpa, min_response); |
173 | rus->speciesResponse()->soilAtmosphereResponses(psi_kpa, vpd_kpa, min_response); |
163 | total_response += min_response * rus->LAIfactor(); |
174 | total_response += min_response * rus->LAIfactor(); |
164 | total_lai_factor += rus->LAIfactor(); |
175 | total_lai_factor += rus->LAIfactor(); |
165 | }
|
176 | }
|
166 | }
|
177 | }
|
167 | 178 | ||
168 | if (total_lai_factor<1.) { |
179 | if (total_lai_factor<1.) { |
169 | // the LAI is below 1: the rest is considered as "ground vegetation"
|
180 | // the LAI is below 1: the rest is considered as "ground vegetation"
|
170 | total_response += calculateBaseSoilAtmosphereResponse(psi_kpa, vpd_kpa) * (1. - total_lai_factor); |
181 | total_response += calculateBaseSoilAtmosphereResponse(psi_kpa, vpd_kpa) * (1. - total_lai_factor); |
171 | }
|
182 | }
|
172 | 183 | ||
173 | // add an aging factor to the total response (averageAging: leaf area weighted mean aging value):
|
184 | // add an aging factor to the total response (averageAging: leaf area weighted mean aging value):
|
174 | // conceptually: response = min(vpd_response, water_response)*aging
|
185 | // conceptually: response = min(vpd_response, water_response)*aging
|
175 | if (total_lai_factor==1.) |
186 | if (total_lai_factor==1.) |
176 | total_response *= mRU->averageAging(); // no ground cover: use aging value for all LA |
187 | total_response *= mRU->averageAging(); // no ground cover: use aging value for all LA |
177 | else if (total_lai_factor>0. && mRU->averageAging()>0.) |
188 | else if (total_lai_factor>0. && mRU->averageAging()>0.) |
178 | total_response *= (1.-total_lai_factor)*1. + (total_lai_factor * mRU->averageAging()); // between 0..1: a part of the LAI is "ground cover" (aging=1) |
189 | total_response *= (1.-total_lai_factor)*1. + (total_lai_factor * mRU->averageAging()); // between 0..1: a part of the LAI is "ground cover" (aging=1) |
179 | 190 | ||
180 | DBGMODE( if (mRU->averageAging()>1. || mRU->averageAging()<0. || total_response<0 || total_response>1.) |
191 | DBGMODE( if (mRU->averageAging()>1. || mRU->averageAging()<0. || total_response<0 || total_response>1.) |
181 | qDebug() << "water cycle: average aging invalid. aging:" << mRU->averageAging() << "total response" << total_response; |
192 | qDebug() << "water cycle: average aging invalid. aging:" << mRU->averageAging() << "total response" << total_response; |
182 | ); |
193 | ); |
183 | 194 | ||
184 | //DBG_IF(mRU->averageAging()>1. || mRU->averageAging()<0.,"water cycle", "average aging invalid!" );
|
195 | //DBG_IF(mRU->averageAging()>1. || mRU->averageAging()<0.,"water cycle", "average aging invalid!" );
|
185 | return total_response; |
196 | return total_response; |
186 | }
|
197 | }
|
187 | 198 | ||
188 | 199 | ||
189 | /// Main Water Cycle function. This function triggers all water related tasks for
|
200 | /// Main Water Cycle function. This function triggers all water related tasks for
|
190 | /// one simulation year.
|
201 | /// one simulation year.
|
191 | /// @sa http://iland.boku.ac.at/water+cycle
|
202 | /// @sa http://iland.boku.ac.at/water+cycle
|
192 | void WaterCycle::run() |
203 | void WaterCycle::run() |
193 | {
|
204 | {
|
194 | // necessary?
|
205 | // necessary?
|
195 | if (GlobalSettings::instance()->currentYear() == mLastYear) |
206 | if (GlobalSettings::instance()->currentYear() == mLastYear) |
196 | return; |
207 | return; |
197 | DebugTimer tw("water:run"); |
208 | DebugTimer tw("water:run"); |
198 | WaterCycleData add_data;
|
209 | WaterCycleData add_data;
|
199 | 210 | ||
200 | // preparations (once a year)
|
211 | // preparations (once a year)
|
201 | getStandValues(); // fetch canopy characteristics from iLand (including weighted average for mCanopyConductance) |
212 | getStandValues(); // fetch canopy characteristics from iLand (including weighted average for mCanopyConductance) |
202 | mCanopy.setStandParameters(mLAINeedle, |
213 | mCanopy.setStandParameters(mLAINeedle, |
203 | mLAIBroadleaved, |
214 | mLAIBroadleaved, |
204 | mCanopyConductance); |
215 | mCanopyConductance); |
205 | 216 | ||
206 | // main loop over all days of the year
|
217 | // main loop over all days of the year
|
207 | double prec_mm, prec_after_interception, prec_to_soil, et, excess; |
218 | double prec_mm, prec_after_interception, prec_to_soil, et, excess; |
208 | const Climate *climate = mRU->climate(); |
219 | const Climate *climate = mRU->climate(); |
209 | const ClimateDay *day = climate->begin(); |
220 | const ClimateDay *day = climate->begin(); |
210 | const ClimateDay *end = climate->end(); |
221 | const ClimateDay *end = climate->end(); |
211 | int doy=0; |
222 | int doy=0; |
212 | double total_excess = 0.; |
223 | double total_excess = 0.; |
213 | for (; day<end; ++day, ++doy) { |
224 | for (; day<end; ++day, ++doy) { |
214 | // (1) precipitation of the day
|
225 | // (1) precipitation of the day
|
215 | prec_mm = day->preciptitation; |
226 | prec_mm = day->preciptitation; |
216 | // (2) interception by the crown
|
227 | // (2) interception by the crown
|
217 | prec_after_interception = mCanopy.flow(prec_mm, day->temperature); |
228 | prec_after_interception = mCanopy.flow(prec_mm, day->temperature); |
218 | // (3) storage in the snow pack
|
229 | // (3) storage in the snow pack
|
219 | prec_to_soil = mSnowPack.flow(prec_after_interception, day->temperature); |
230 | prec_to_soil = mSnowPack.flow(prec_after_interception, day->temperature); |
220 | // save extra data (used by e.g. fire module)
|
231 | // save extra data (used by e.g. fire module)
|
221 | add_data.water_to_ground[doy] = prec_to_soil; |
232 | add_data.water_to_ground[doy] = prec_to_soil; |
222 | add_data.snow_cover[doy] = mSnowPack.snowPack(); |
233 | add_data.snow_cover[doy] = mSnowPack.snowPack(); |
223 | // (4) add rest to soil
|
234 | // (4) add rest to soil
|
224 | mContent += prec_to_soil; |
235 | mContent += prec_to_soil; |
225 | 236 | ||
226 | excess = 0.; |
237 | excess = 0.; |
227 | if (mContent>mFieldCapacity) { |
238 | if (mContent>mFieldCapacity) { |
228 | // excess water runoff
|
239 | // excess water runoff
|
229 | excess = mContent - mFieldCapacity; |
240 | excess = mContent - mFieldCapacity; |
230 | total_excess += excess; |
241 | total_excess += excess; |
231 | mContent = mFieldCapacity; |
242 | mContent = mFieldCapacity; |
232 | }
|
243 | }
|
233 | 244 | ||
234 | double current_psi = psiFromHeight(mContent); |
245 | double current_psi = psiFromHeight(mContent); |
235 | mPsi[doy] = current_psi; |
246 | mPsi[doy] = current_psi; |
236 | 247 | ||
237 | // (5) transpiration of the vegetation (and of water intercepted in canopy)
|
248 | // (5) transpiration of the vegetation (and of water intercepted in canopy)
|
238 | // calculate the LAI-weighted response values for soil water and vpd:
|
249 | // calculate the LAI-weighted response values for soil water and vpd:
|
239 | double interception_before_transpiration = mCanopy.interception(); |
250 | double interception_before_transpiration = mCanopy.interception(); |
240 | double combined_response = calculateSoilAtmosphereResponse( current_psi, day->vpd); |
251 | double combined_response = calculateSoilAtmosphereResponse( current_psi, day->vpd); |
241 | et = mCanopy.evapotranspiration3PG(day, climate->daylength_h(doy), combined_response); |
252 | et = mCanopy.evapotranspiration3PG(day, climate->daylength_h(doy), combined_response); |
242 | // if there is some flow from intercepted water to the ground -> add to "water_to_the_ground"
|
253 | // if there is some flow from intercepted water to the ground -> add to "water_to_the_ground"
|
243 | if (mCanopy.interception() < interception_before_transpiration) |
254 | if (mCanopy.interception() < interception_before_transpiration) |
244 | add_data.water_to_ground[doy]+= interception_before_transpiration - mCanopy.interception(); |
255 | add_data.water_to_ground[doy]+= interception_before_transpiration - mCanopy.interception(); |
245 | 256 | ||
246 | mContent -= et; // reduce content (transpiration) |
257 | mContent -= et; // reduce content (transpiration) |
247 | // add intercepted water (that is *not* evaporated) again to the soil (or add to snow if temp too low -> call to snowpack)
|
258 | // add intercepted water (that is *not* evaporated) again to the soil (or add to snow if temp too low -> call to snowpack)
|
248 | mContent += mSnowPack.add(mCanopy.interception(),day->temperature); |
259 | mContent += mSnowPack.add(mCanopy.interception(),day->temperature); |
249 | 260 | ||
250 | // do not remove water below the PWP (fixed value)
|
261 | // do not remove water below the PWP (fixed value)
|
251 | if (mContent<mPermanentWiltingPoint) { |
262 | if (mContent<mPermanentWiltingPoint) { |
252 | et -= mPermanentWiltingPoint - mContent; // reduce et (for bookkeeping) |
263 | et -= mPermanentWiltingPoint - mContent; // reduce et (for bookkeeping) |
253 | mContent = mPermanentWiltingPoint; |
264 | mContent = mPermanentWiltingPoint; |
254 | }
|
265 | }
|
255 | 266 | ||
256 | 267 | ||
257 | //DBGMODE(
|
268 | //DBGMODE(
|
258 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dWaterCycle)) { |
269 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dWaterCycle)) { |
259 | DebugList &out = GlobalSettings::instance()->debugList(day->id(), GlobalSettings::dWaterCycle); |
270 | DebugList &out = GlobalSettings::instance()->debugList(day->id(), GlobalSettings::dWaterCycle); |
260 | // climatic variables
|
271 | // climatic variables
|
261 | out << day->id() << mRU->index() << mRU->id() << day->temperature << day->vpd << day->preciptitation << day->radiation; |
272 | out << day->id() << mRU->index() << mRU->id() << day->temperature << day->vpd << day->preciptitation << day->radiation; |
262 | out << combined_response; // combined response of all species on RU (min(water, vpd)) |
273 | out << combined_response; // combined response of all species on RU (min(water, vpd)) |
263 | // fluxes
|
274 | // fluxes
|
264 | out << prec_after_interception << prec_to_soil << et << mCanopy.evaporationCanopy() |
275 | out << prec_after_interception << prec_to_soil << et << mCanopy.evaporationCanopy() |
265 | << mContent << mPsi[doy] << excess; |
276 | << mContent << mPsi[doy] << excess; |
266 | // other states
|
277 | // other states
|
267 | out << mSnowPack.snowPack(); |
278 | out << mSnowPack.snowPack(); |
268 | //special sanity check:
|
279 | //special sanity check:
|
269 | if (prec_to_soil>0. && mCanopy.interception()>0.) |
280 | if (prec_to_soil>0. && mCanopy.interception()>0.) |
270 | if (mSnowPack.snowPack()==0. && day->preciptitation==0) |
281 | if (mSnowPack.snowPack()==0. && day->preciptitation==0) |
271 | qDebug() << "watercontent increase without precipititaion"; |
282 | qDebug() << "watercontent increase without precipititaion"; |
272 | 283 | ||
273 | }
|
284 | }
|
274 | //); // DBGMODE()
|
285 | //); // DBGMODE()
|
275 | 286 | ||
276 | }
|
287 | }
|
277 | // call external modules
|
288 | // call external modules
|
278 | GlobalSettings::instance()->model()->modules()->calculateWater(mRU, &add_data); |
289 | GlobalSettings::instance()->model()->modules()->calculateWater(mRU, &add_data); |
279 | mLastYear = GlobalSettings::instance()->currentYear(); |
290 | mLastYear = GlobalSettings::instance()->currentYear(); |
280 | 291 | ||
281 | }
|
292 | }
|
282 | 293 | ||
283 | 294 | ||
284 | namespace Water { |
295 | namespace Water { |
285 | 296 | ||
286 | /** calculates the input/output of water that is stored in the snow pack.
|
297 | /** calculates the input/output of water that is stored in the snow pack.
|
287 | The approach is similar to Picus 1.3 and ForestBGC (Running, 1988).
|
298 | The approach is similar to Picus 1.3 and ForestBGC (Running, 1988).
|
288 | Returns the amount of water that exits the snowpack (precipitation, snow melt) */
|
299 | Returns the amount of water that exits the snowpack (precipitation, snow melt) */
|
289 | double SnowPack::flow(const double &preciptitation_mm, const double &temperature) |
300 | double SnowPack::flow(const double &preciptitation_mm, const double &temperature) |
290 | {
|
301 | {
|
291 | if (temperature>mSnowTemperature) { |
302 | if (temperature>mSnowTemperature) { |
292 | if (mSnowPack==0.) |
303 | if (mSnowPack==0.) |
293 | return preciptitation_mm; // no change |
304 | return preciptitation_mm; // no change |
294 | else { |
305 | else { |
295 | // snow melts
|
306 | // snow melts
|
296 | const double melting_coefficient = 0.7; // mm/°C |
307 | const double melting_coefficient = 0.7; // mm/°C |
297 | double melt = qMin( (temperature-mSnowTemperature) * melting_coefficient, mSnowPack); |
308 | double melt = qMin( (temperature-mSnowTemperature) * melting_coefficient, mSnowPack); |
298 | mSnowPack -=melt; |
309 | mSnowPack -=melt; |
299 | return preciptitation_mm + melt; |
310 | return preciptitation_mm + melt; |
300 | }
|
311 | }
|
301 | } else { |
312 | } else { |
302 | // snow:
|
313 | // snow:
|
303 | mSnowPack += preciptitation_mm; |
314 | mSnowPack += preciptitation_mm; |
304 | return 0.; // no output. |
315 | return 0.; // no output. |
305 | }
|
316 | }
|
306 | 317 | ||
307 | }
|
318 | }
|
308 | 319 | ||
309 | 320 | ||
310 | inline double SnowPack::add(const double &preciptitation_mm, const double &temperature) |
321 | inline double SnowPack::add(const double &preciptitation_mm, const double &temperature) |
311 | {
|
322 | {
|
312 | // do nothing for temps > 0°
|
323 | // do nothing for temps > 0°
|
313 | if (temperature>mSnowTemperature) |
324 | if (temperature>mSnowTemperature) |
314 | return preciptitation_mm; |
325 | return preciptitation_mm; |
315 | 326 | ||
316 | // temps < 0°: add to snow
|
327 | // temps < 0°: add to snow
|
317 | mSnowPack += preciptitation_mm; |
328 | mSnowPack += preciptitation_mm; |
318 | return 0.; |
329 | return 0.; |
319 | }
|
330 | }
|
320 | 331 | ||
321 | /** Interception in crown canopy.
|
332 | /** Interception in crown canopy.
|
322 | Calculates the amount of preciptitation that does not reach the ground and
|
333 | Calculates the amount of preciptitation that does not reach the ground and
|
323 | is stored in the canopy. The approach is adopted from Picus 1.3.
|
334 | is stored in the canopy. The approach is adopted from Picus 1.3.
|
324 | Returns the amount of precipitation (mm) that surpasses the canopy layer.
|
335 | Returns the amount of precipitation (mm) that surpasses the canopy layer.
|
325 | @sa http://iland.boku.ac.at/water+cycle#precipitation_and_interception */
|
336 | @sa http://iland.boku.ac.at/water+cycle#precipitation_and_interception */
|
326 | double Canopy::flow(const double &preciptitation_mm, const double &temperature) |
337 | double Canopy::flow(const double &preciptitation_mm, const double &temperature) |
327 | {
|
338 | {
|
328 | // sanity checks
|
339 | // sanity checks
|
329 | mInterception = 0.; |
340 | mInterception = 0.; |
330 | mEvaporation = 0.; |
341 | mEvaporation = 0.; |
331 | if (!mLAI) |
342 | if (!mLAI) |
332 | return preciptitation_mm; |
343 | return preciptitation_mm; |
333 | if (!preciptitation_mm) |
344 | if (!preciptitation_mm) |
334 | return 0.; |
345 | return 0.; |
335 | double max_interception_mm=0.; // maximum interception based on the current foliage |
346 | double max_interception_mm=0.; // maximum interception based on the current foliage |
336 | double max_storage_mm=0.; // maximum storage in canopy |
347 | double max_storage_mm=0.; // maximum storage in canopy |
337 | 348 | ||
338 | if (mLAINeedle>0.) { |
349 | if (mLAINeedle>0.) { |
339 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation)
|
350 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation)
|
340 | double max_flow_needle = 0.9 * sqrt(1.03 - exp(-0.055*preciptitation_mm)); |
351 | double max_flow_needle = 0.9 * sqrt(1.03 - exp(-0.055*preciptitation_mm)); |
341 | max_interception_mm += preciptitation_mm * (1. - max_flow_needle * mLAINeedle/mLAI); |
352 | max_interception_mm += preciptitation_mm * (1. - max_flow_needle * mLAINeedle/mLAI); |
342 | // (2) calculate maximum storage potential based on the current LAI
|
353 | // (2) calculate maximum storage potential based on the current LAI
|
343 | double max_storage_needle = mNeedleFactor * (1. - exp(-0.55*mLAINeedle) ); |
354 | double max_storage_needle = mNeedleFactor * (1. - exp(-0.55*mLAINeedle) ); |
344 | max_storage_mm += max_storage_needle; |
355 | max_storage_mm += max_storage_needle; |
345 | }
|
356 | }
|
346 | 357 | ||
347 | if (mLAIBroadleaved>0.) { |
358 | if (mLAIBroadleaved>0.) { |
348 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation)
|
359 | // (1) calculate maximum fraction of thru-flow the crown (based on precipitation)
|
349 | double max_flow_broad = 0.9 * pow(1.22 - exp(-0.055*preciptitation_mm), 0.35); |
360 | double max_flow_broad = 0.9 * pow(1.22 - exp(-0.055*preciptitation_mm), 0.35); |
350 | max_interception_mm += preciptitation_mm * (1. - max_flow_broad * mLAIBroadleaved/mLAI); |
361 | max_interception_mm += preciptitation_mm * (1. - max_flow_broad * mLAIBroadleaved/mLAI); |
351 | // (2) calculate maximum storage potential based on the current LAI
|
362 | // (2) calculate maximum storage potential based on the current LAI
|
352 | double max_storage_broad = mDecidousFactor * (1. - exp(-0.5*mLAIBroadleaved) ); |
363 | double max_storage_broad = mDecidousFactor * (1. - exp(-0.5*mLAIBroadleaved) ); |
353 | max_storage_mm += max_storage_broad; |
364 | max_storage_mm += max_storage_broad; |
354 | }
|
365 | }
|
355 | 366 | ||
356 | // (3) calculate actual interception and store for evaporation calculation
|
367 | // (3) calculate actual interception and store for evaporation calculation
|
357 | mInterception = qMin( max_storage_mm, max_interception_mm ); |
368 | mInterception = qMin( max_storage_mm, max_interception_mm ); |
358 | 369 | ||
359 | // (4) limit interception with amount of precipitation
|
370 | // (4) limit interception with amount of precipitation
|
360 | mInterception = qMin( mInterception, preciptitation_mm); |
371 | mInterception = qMin( mInterception, preciptitation_mm); |
361 | 372 | ||
362 | // (5) reduce preciptitaion by the amount that is intercepted by the canopy
|
373 | // (5) reduce preciptitaion by the amount that is intercepted by the canopy
|
363 | return preciptitation_mm - mInterception; |
374 | return preciptitation_mm - mInterception; |
364 | 375 | ||
365 | }
|
376 | }
|
366 | 377 | ||
367 | /// sets up the canopy. fetch some global parameter values...
|
378 | /// sets up the canopy. fetch some global parameter values...
|
368 | void Canopy::setup() |
379 | void Canopy::setup() |
369 | {
|
380 | {
|
370 | mAirDensity = Model::settings().airDensity; // kg / m3 |
381 | mAirDensity = Model::settings().airDensity; // kg / m3 |
371 | }
|
382 | }
|
372 | 383 | ||
373 | void Canopy::setStandParameters(const double LAIneedle, const double LAIbroadleave, const double maxCanopyConductance) |
384 | void Canopy::setStandParameters(const double LAIneedle, const double LAIbroadleave, const double maxCanopyConductance) |
374 | {
|
385 | {
|
375 | mLAINeedle = LAIneedle; |
386 | mLAINeedle = LAIneedle; |
376 | mLAIBroadleaved=LAIbroadleave; |
387 | mLAIBroadleaved=LAIbroadleave; |
377 | mLAI=LAIneedle+LAIbroadleave; |
388 | mLAI=LAIneedle+LAIbroadleave; |
378 | mAvgMaxCanopyConductance = maxCanopyConductance; |
389 | mAvgMaxCanopyConductance = maxCanopyConductance; |
379 | 390 | ||
380 | // clear aggregation containers
|
391 | // clear aggregation containers
|
381 | for (int i=0;i<12;++i) mET0[i]=0.; |
392 | for (int i=0;i<12;++i) mET0[i]=0.; |
382 | 393 | ||
383 | }
|
394 | }
|
384 | 395 | ||
385 | 396 | ||
386 | 397 | ||
387 | /** calculate the daily evaporation/transpiration using the Penman-Monteith-Equation.
|
398 | /** calculate the daily evaporation/transpiration using the Penman-Monteith-Equation.
|
388 | This version is based on 3PG. See the Visual Basic Code in 3PGjs.xls.
|
399 | This version is based on 3PG. See the Visual Basic Code in 3PGjs.xls.
|
389 | Returns the total sum of evaporation+transpiration in mm of the day. */
|
400 | Returns the total sum of evaporation+transpiration in mm of the day. */
|
390 | double Canopy::evapotranspiration3PG(const ClimateDay *climate, const double daylength_h, const double combined_response) |
401 | double Canopy::evapotranspiration3PG(const ClimateDay *climate, const double daylength_h, const double combined_response) |
391 | {
|
402 | {
|
392 | double vpd_mbar = climate->vpd * 10.; // convert from kPa to mbar |
403 | double vpd_mbar = climate->vpd * 10.; // convert from kPa to mbar |
393 | double temperature = climate->temperature; // average temperature of the day (°C) |
404 | double temperature = climate->temperature; // average temperature of the day (°C) |
394 | double daylength = daylength_h * 3600.; // daylength in seconds (convert from length in hours) |
405 | double daylength = daylength_h * 3600.; // daylength in seconds (convert from length in hours) |
395 | double rad = climate->radiation / daylength * 1000000; //convert from MJ/m2 (day sum) to average radiation flow W/m2 [MJ=MWs -> /s * 1,000,000 |
406 | double rad = climate->radiation / daylength * 1000000; //convert from MJ/m2 (day sum) to average radiation flow W/m2 [MJ=MWs -> /s * 1,000,000 |
396 | 407 | ||
397 | // the radiation: based on linear empirical function
|
408 | // the radiation: based on linear empirical function
|
398 | const double qa = -90.; |
409 | const double qa = -90.; |
399 | const double qb = 0.8; |
410 | const double qb = 0.8; |
400 | double net_rad = qa + qb*rad; |
411 | double net_rad = qa + qb*rad; |
401 | 412 | ||
402 | //: Landsberg original: const double e20 = 2.2; //rate of change of saturated VP with T at 20C
|
413 | //: Landsberg original: const double e20 = 2.2; //rate of change of saturated VP with T at 20C
|
403 | const double VPDconv = 0.000622; //convert VPD to saturation deficit = 18/29/1000 |
414 | const double VPDconv = 0.000622; //convert VPD to saturation deficit = 18/29/1000 |
404 | const double latent_heat = 2460000.; // Latent heat of vaporization. Energy required per unit mass of water vaporized [J kg-1] |
415 | const double latent_heat = 2460000.; // Latent heat of vaporization. Energy required per unit mass of water vaporized [J kg-1] |
405 | 416 | ||
406 | double gBL = Model::settings().boundaryLayerConductance; // boundary layer conductance |
417 | double gBL = Model::settings().boundaryLayerConductance; // boundary layer conductance |
407 | 418 | ||
408 | // canopy conductance.
|
419 | // canopy conductance.
|
409 | // The species traits are weighted by LAI on the RU.
|
420 | // The species traits are weighted by LAI on the RU.
|
410 | // maximum canopy conductance: see getStandValues()
|
421 | // maximum canopy conductance: see getStandValues()
|
411 | // current response: see calculateSoilAtmosphereResponse(). This is basically a weighted average of min(water_response, vpd_response) for each species
|
422 | // current response: see calculateSoilAtmosphereResponse(). This is basically a weighted average of min(water_response, vpd_response) for each species
|
412 | double gC = mAvgMaxCanopyConductance * combined_response; |
423 | double gC = mAvgMaxCanopyConductance * combined_response; |
413 | 424 | ||
414 | 425 | ||
415 | double defTerm = mAirDensity * latent_heat * (vpd_mbar * VPDconv) * gBL; |
426 | double defTerm = mAirDensity * latent_heat * (vpd_mbar * VPDconv) * gBL; |
416 | 427 | ||
417 | // with temperature-dependent slope of vapor pressure saturation curve
|
428 | // with temperature-dependent slope of vapor pressure saturation curve
|
418 | // (following Allen et al. (1998), http://www.fao.org/docrep/x0490e/x0490e07.htm#atmospheric%20parameters)
|
429 | // (following Allen et al. (1998), http://www.fao.org/docrep/x0490e/x0490e07.htm#atmospheric%20parameters)
|
419 | // svp_slope in mbar.
|
430 | // svp_slope in mbar.
|
420 | //double svp_slope = 4098. * (6.1078 * exp(17.269 * temperature / (temperature + 237.3))) / ((237.3+temperature)*(237.3+temperature));
|
431 | //double svp_slope = 4098. * (6.1078 * exp(17.269 * temperature / (temperature + 237.3))) / ((237.3+temperature)*(237.3+temperature));
|
421 | 432 | ||
422 | // alternatively: very simple variant (following here the original 3PG code). This
|
433 | // alternatively: very simple variant (following here the original 3PG code). This
|
423 | // keeps yields +- same results for summer, but slightly lower values in winter (2011/03/16)
|
434 | // keeps yields +- same results for summer, but slightly lower values in winter (2011/03/16)
|
424 | double svp_slope = 2.2; |
435 | double svp_slope = 2.2; |
425 | 436 | ||
426 | double div = (1. + svp_slope + gBL / gC); |
437 | double div = (1. + svp_slope + gBL / gC); |
427 | double Etransp = (svp_slope * net_rad + defTerm) / div; |
438 | double Etransp = (svp_slope * net_rad + defTerm) / div; |
428 | double canopy_transpiration = Etransp / latent_heat * daylength; |
439 | double canopy_transpiration = Etransp / latent_heat * daylength; |
429 | 440 | ||
430 | // calculate reference evapotranspiration
|
441 | // calculate reference evapotranspiration
|
431 | // see Adair et al 2008
|
442 | // see Adair et al 2008
|
432 | const double psychrometric_const = 0.0672718682328237; // kPa/degC |
443 | const double psychrometric_const = 0.0672718682328237; // kPa/degC |
433 | const double windspeed = 2.; // m/s |
444 | const double windspeed = 2.; // m/s |
434 | double net_rad_mj_day = net_rad*daylength/1000000.; // convert W/m2 again to MJ/m2*day |
445 | double net_rad_mj_day = net_rad*daylength/1000000.; // convert W/m2 again to MJ/m2*day |
435 | double et0_day = 0.408*svp_slope*net_rad_mj_day + psychrometric_const*900./(temperature+273.)*windspeed*climate->vpd; |
446 | double et0_day = 0.408*svp_slope*net_rad_mj_day + psychrometric_const*900./(temperature+273.)*windspeed*climate->vpd; |
436 | double et0_div = svp_slope+psychrometric_const*(1.+0.34*windspeed); |
447 | double et0_div = svp_slope+psychrometric_const*(1.+0.34*windspeed); |
437 | et0_day = et0_day / et0_div; |
448 | et0_day = et0_day / et0_div; |
438 | mET0[climate->month-1] += et0_day; |
449 | mET0[climate->month-1] += et0_day; |
439 | 450 | ||
440 | if (mInterception>0.) { |
451 | if (mInterception>0.) { |
441 | // we assume that for evaporation from leaf surface gBL/gC -> 0
|
452 | // we assume that for evaporation from leaf surface gBL/gC -> 0
|
442 | double div_evap = 1. + svp_slope; |
453 | double div_evap = 1. + svp_slope; |
443 | double evap_canopy_potential = (svp_slope*net_rad + defTerm) / div_evap / latent_heat * daylength; |
454 | double evap_canopy_potential = (svp_slope*net_rad + defTerm) / div_evap / latent_heat * daylength; |
444 | // reduce the amount of transpiration on a wet day based on the approach of
|
455 | // reduce the amount of transpiration on a wet day based on the approach of
|
445 | // Wigmosta et al (1994). see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance
|
456 | // Wigmosta et al (1994). see http://iland.boku.ac.at/water+cycle#transpiration_and_canopy_conductance
|
446 | 457 | ||
447 | double ratio_T_E = canopy_transpiration / evap_canopy_potential; |
458 | double ratio_T_E = canopy_transpiration / evap_canopy_potential; |
448 | double evap_canopy = qMin(evap_canopy_potential, mInterception); |
459 | double evap_canopy = qMin(evap_canopy_potential, mInterception); |
449 | 460 | ||
450 | // for interception -> 0, the canopy transpiration is unchanged
|
461 | // for interception -> 0, the canopy transpiration is unchanged
|
451 | canopy_transpiration = (evap_canopy_potential - evap_canopy) * ratio_T_E; |
462 | canopy_transpiration = (evap_canopy_potential - evap_canopy) * ratio_T_E; |
452 | 463 | ||
453 | mInterception -= evap_canopy; // reduce interception |
464 | mInterception -= evap_canopy; // reduce interception |
454 | mEvaporation = evap_canopy; // evaporation from intercepted water |
465 | mEvaporation = evap_canopy; // evaporation from intercepted water |
455 | 466 | ||
456 | }
|
467 | }
|
457 | return canopy_transpiration; |
468 | return canopy_transpiration; |
458 | }
|
469 | }
|
459 | 470 | ||
460 | } // end namespace |
471 | } // end namespace |
461 | 472 |