Rev 734 | Rev 779 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 734 | Rev 738 | ||
---|---|---|---|
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp': |
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.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 | /** @class ResourceUnit
|
21 | /** @class ResourceUnit
|
22 | ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
|
22 | ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
|
23 | (Climate, Soil, Water, ...).
|
23 | (Climate, Soil, Water, ...).
|
24 | @ingroup core
|
24 | @ingroup core
|
25 | A resource unit has a size of (currently) 100x100m. Many processes in iLand operate on the level of a ResourceUnit.
|
25 | A resource unit has a size of (currently) 100x100m. Many processes in iLand operate on the level of a ResourceUnit.
|
26 | Each resource unit has the same Climate and other properties (e.g. available nitrogen).
|
26 | Each resource unit has the same Climate and other properties (e.g. available nitrogen).
|
27 | Proceses on this level are, inter alia, NPP Production (see Production3PG), water calculations (WaterCycle), the modeling
|
27 | Proceses on this level are, inter alia, NPP Production (see Production3PG), water calculations (WaterCycle), the modeling
|
28 | of dead trees (Snag) and soil processes (Soil).
|
28 | of dead trees (Snag) and soil processes (Soil).
|
29 | 29 | ||
30 | */
|
30 | */
|
31 | #include <QtCore>
|
31 | #include <QtCore>
|
32 | #include "global.h"
|
32 | #include "global.h"
|
33 | 33 | ||
34 | #include "resourceunit.h"
|
34 | #include "resourceunit.h"
|
35 | #include "resourceunitspecies.h"
|
35 | #include "resourceunitspecies.h"
|
36 | #include "speciesset.h"
|
36 | #include "speciesset.h"
|
37 | #include "species.h"
|
37 | #include "species.h"
|
38 | #include "production3pg.h"
|
38 | #include "production3pg.h"
|
39 | #include "model.h"
|
39 | #include "model.h"
|
40 | #include "climate.h"
|
40 | #include "climate.h"
|
41 | #include "watercycle.h"
|
41 | #include "watercycle.h"
|
42 | #include "snag.h"
|
42 | #include "snag.h"
|
43 | #include "soil.h"
|
43 | #include "soil.h"
|
44 | #include "helper.h"
|
44 | #include "helper.h"
|
45 | 45 | ||
46 | ResourceUnit::~ResourceUnit() |
46 | ResourceUnit::~ResourceUnit() |
47 | {
|
47 | {
|
48 | if (mWater) |
48 | if (mWater) |
49 | delete mWater; |
49 | delete mWater; |
50 | mWater = 0; |
50 | mWater = 0; |
51 | if (mSnag) |
51 | if (mSnag) |
52 | delete mSnag; |
52 | delete mSnag; |
53 | if (mSoil) |
53 | if (mSoil) |
54 | delete mSoil; |
54 | delete mSoil; |
- | 55 | ||
- | 56 | qDeleteAll(mRUSpecies); |
|
55 | 57 | ||
56 | mSnag = 0; |
58 | mSnag = 0; |
57 | mSoil = 0; |
59 | mSoil = 0; |
58 | }
|
60 | }
|
59 | 61 | ||
60 | ResourceUnit::ResourceUnit(const int index) |
62 | ResourceUnit::ResourceUnit(const int index) |
61 | {
|
63 | {
|
62 | qDeleteAll(mRUSpecies); |
64 | qDeleteAll(mRUSpecies); |
63 | mSpeciesSet = 0; |
65 | mSpeciesSet = 0; |
64 | mClimate = 0; |
66 | mClimate = 0; |
65 | mPixelCount=0; |
67 | mPixelCount=0; |
66 | mStockedArea = 0; |
68 | mStockedArea = 0; |
67 | mStockedPixelCount = 0; |
69 | mStockedPixelCount = 0; |
68 | mIndex = index; |
70 | mIndex = index; |
69 | mSaplingHeightMap = 0; |
71 | mSaplingHeightMap = 0; |
70 | mEffectiveArea_perWLA = 0.; |
72 | mEffectiveArea_perWLA = 0.; |
71 | mWater = new WaterCycle(); |
73 | mWater = new WaterCycle(); |
72 | mSnag = 0; |
74 | mSnag = 0; |
73 | mSoil = 0; |
75 | mSoil = 0; |
74 | mID = 0; |
76 | mID = 0; |
75 | }
|
77 | }
|
76 | 78 | ||
77 | void ResourceUnit::setup() |
79 | void ResourceUnit::setup() |
78 | {
|
80 | {
|
79 | mWater->setup(this); |
81 | mWater->setup(this); |
80 | 82 | ||
81 | if (mSnag) |
83 | if (mSnag) |
82 | delete mSnag; |
84 | delete mSnag; |
83 | mSnag=0; |
85 | mSnag=0; |
84 | if (mSoil) |
86 | if (mSoil) |
85 | delete mSoil; |
87 | delete mSoil; |
86 | mSoil=0; |
88 | mSoil=0; |
87 | if (Model::settings().carbonCycleEnabled) { |
89 | if (Model::settings().carbonCycleEnabled) { |
88 | mSoil = new Soil(this); |
90 | mSoil = new Soil(this); |
89 | mSnag = new Snag; |
91 | mSnag = new Snag; |
90 | mSnag->setup(this); |
92 | mSnag->setup(this); |
91 | const XmlHelper &xml=GlobalSettings::instance()->settings(); |
93 | const XmlHelper &xml=GlobalSettings::instance()->settings(); |
92 | 94 | ||
93 | // setup contents of the soil of the RU; use values for C and N (kg/ha)
|
95 | // setup contents of the soil of the RU; use values for C and N (kg/ha)
|
94 | mSoil->setInitialState(CNPool(xml.valueDouble("model.site.youngLabileC", -1), |
96 | mSoil->setInitialState(CNPool(xml.valueDouble("model.site.youngLabileC", -1), |
95 | xml.valueDouble("model.site.youngLabileN", -1), |
97 | xml.valueDouble("model.site.youngLabileN", -1), |
96 | xml.valueDouble("model.site.youngLabileDecompRate", -1)), |
98 | xml.valueDouble("model.site.youngLabileDecompRate", -1)), |
97 | CNPool(xml.valueDouble("model.site.youngRefractoryC", -1), |
99 | CNPool(xml.valueDouble("model.site.youngRefractoryC", -1), |
98 | xml.valueDouble("model.site.youngRefractoryN", -1), |
100 | xml.valueDouble("model.site.youngRefractoryN", -1), |
99 | xml.valueDouble("model.site.youngRefractoryDecompRate", -1)), |
101 | xml.valueDouble("model.site.youngRefractoryDecompRate", -1)), |
100 | CNPair(xml.valueDouble("model.site.somC", -1), xml.valueDouble("model.site.somN", -1))); |
102 | CNPair(xml.valueDouble("model.site.somC", -1), xml.valueDouble("model.site.somN", -1))); |
101 | }
|
103 | }
|
102 | 104 | ||
103 | // setup variables
|
105 | // setup variables
|
104 | mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40); |
106 | mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40); |
105 | 107 | ||
106 | // if dynamic coupling of soil nitrogen is enabled, the calculate a starting value for available n.
|
108 | // if dynamic coupling of soil nitrogen is enabled, the calculate a starting value for available n.
|
107 | if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) { |
109 | if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) { |
108 | mSoil->setClimateFactor(1.); |
110 | mSoil->setClimateFactor(1.); |
109 | mSoil->calculateYear(); |
111 | mSoil->calculateYear(); |
110 | mUnitVariables.nitrogenAvailable = mSoil->availableNitrogen(); |
112 | mUnitVariables.nitrogenAvailable = mSoil->availableNitrogen(); |
111 | }
|
113 | }
|
112 | mHasDeadTrees = false; |
114 | mHasDeadTrees = false; |
113 | mAverageAging = 0.; |
115 | mAverageAging = 0.; |
114 | 116 | ||
115 | }
|
117 | }
|
116 | void ResourceUnit::setBoundingBox(const QRectF &bb) |
118 | void ResourceUnit::setBoundingBox(const QRectF &bb) |
117 | {
|
119 | {
|
118 | mBoundingBox = bb; |
120 | mBoundingBox = bb; |
119 | mCornerCoord = GlobalSettings::instance()->model()->grid()->indexAt(bb.topLeft()); |
121 | mCornerCoord = GlobalSettings::instance()->model()->grid()->indexAt(bb.topLeft()); |
120 | }
|
122 | }
|
121 | 123 | ||
122 | /// set species and setup the species-per-RU-data
|
124 | /// set species and setup the species-per-RU-data
|
123 | void ResourceUnit::setSpeciesSet(SpeciesSet *set) |
125 | void ResourceUnit::setSpeciesSet(SpeciesSet *set) |
124 | {
|
126 | {
|
125 | mSpeciesSet = set; |
127 | mSpeciesSet = set; |
126 | qDeleteAll(mRUSpecies); |
128 | qDeleteAll(mRUSpecies); |
127 | 129 | ||
128 | //mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
|
130 | //mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
|
129 | for (int i=0;i<set->count();i++) { |
131 | for (int i=0;i<set->count();i++) { |
130 | Species *s = const_cast<Species*>(mSpeciesSet->species(i)); |
132 | Species *s = const_cast<Species*>(mSpeciesSet->species(i)); |
131 | if (!s) |
133 | if (!s) |
132 | throw IException("ResourceUnit::setSpeciesSet: invalid index!"); |
134 | throw IException("ResourceUnit::setSpeciesSet: invalid index!"); |
133 | 135 | ||
134 | ResourceUnitSpecies *rus = new ResourceUnitSpecies(); |
136 | ResourceUnitSpecies *rus = new ResourceUnitSpecies(); |
135 | mRUSpecies.push_back(rus); |
137 | mRUSpecies.push_back(rus); |
136 | rus->setup(s, this); |
138 | rus->setup(s, this); |
137 | /* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
|
139 | /* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
|
138 | If the container memory is relocated (QVector), the pointer gets invalid!!!
|
140 | If the container memory is relocated (QVector), the pointer gets invalid!!!
|
139 | Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
|
141 | Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
|
140 | //mRUSpecies[i].setup(s,this); // setup this element
|
142 | //mRUSpecies[i].setup(s,this); // setup this element
|
141 | 143 | ||
142 | }
|
144 | }
|
143 | }
|
145 | }
|
144 | 146 | ||
145 | ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species) |
147 | ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species) |
146 | {
|
148 | {
|
147 | return *mRUSpecies[species->index()]; |
149 | return *mRUSpecies[species->index()]; |
148 | }
|
150 | }
|
149 | 151 | ||
150 | Tree &ResourceUnit::newTree() |
152 | Tree &ResourceUnit::newTree() |
151 | {
|
153 | {
|
152 | // start simple: just append to the vector...
|
154 | // start simple: just append to the vector...
|
153 | if (mTrees.isEmpty()) |
155 | if (mTrees.isEmpty()) |
154 | mTrees.reserve(100); // reserve a junk of memory for trees |
156 | mTrees.reserve(100); // reserve a junk of memory for trees |
155 | 157 | ||
156 | mTrees.append(Tree()); |
158 | mTrees.append(Tree()); |
157 | return mTrees.back(); |
159 | return mTrees.back(); |
158 | }
|
160 | }
|
159 | int ResourceUnit::newTreeIndex() |
161 | int ResourceUnit::newTreeIndex() |
160 | {
|
162 | {
|
161 | newTree(); |
163 | newTree(); |
162 | return mTrees.count()-1; // return index of the last tree |
164 | return mTrees.count()-1; // return index of the last tree |
163 | }
|
165 | }
|
164 | 166 | ||
165 | /// remove dead trees from tree list
|
167 | /// remove dead trees from tree list
|
166 | /// reduce size of vector if lots of space is free
|
168 | /// reduce size of vector if lots of space is free
|
167 | /// tests showed that this way of cleanup is very fast,
|
169 | /// tests showed that this way of cleanup is very fast,
|
168 | /// because no memory allocations are performed (simple memmove())
|
170 | /// because no memory allocations are performed (simple memmove())
|
169 | /// when trees are moved.
|
171 | /// when trees are moved.
|
170 | void ResourceUnit::cleanTreeList() |
172 | void ResourceUnit::cleanTreeList() |
171 | {
|
173 | {
|
172 | if (!mHasDeadTrees) |
174 | if (!mHasDeadTrees) |
173 | return; |
175 | return; |
174 | 176 | ||
175 | QVector<Tree>::iterator last=mTrees.end()-1; |
177 | QVector<Tree>::iterator last=mTrees.end()-1; |
176 | QVector<Tree>::iterator current = mTrees.begin(); |
178 | QVector<Tree>::iterator current = mTrees.begin(); |
177 | while (last>=current && (*last).isDead()) |
179 | while (last>=current && (*last).isDead()) |
178 | --last; |
180 | --last; |
179 | 181 | ||
180 | while (current<last) { |
182 | while (current<last) { |
181 | if ((*current).isDead()) { |
183 | if ((*current).isDead()) { |
182 | *current = *last; // copy data! |
184 | *current = *last; // copy data! |
183 | --last; // |
185 | --last; // |
184 | while (last>=current && (*last).isDead()) |
186 | while (last>=current && (*last).isDead()) |
185 | --last; |
187 | --last; |
186 | }
|
188 | }
|
187 | ++current; |
189 | ++current; |
188 | }
|
190 | }
|
189 | ++last; // last points now to the first dead tree |
191 | ++last; // last points now to the first dead tree |
190 | 192 | ||
191 | // free ressources
|
193 | // free ressources
|
192 | if (last!=mTrees.end()) { |
194 | if (last!=mTrees.end()) { |
193 | mTrees.erase(last, mTrees.end()); |
195 | mTrees.erase(last, mTrees.end()); |
194 | if (mTrees.capacity()>100) { |
196 | if (mTrees.capacity()>100) { |
195 | if (mTrees.count() / double(mTrees.capacity()) < 0.2) { |
197 | if (mTrees.count() / double(mTrees.capacity()) < 0.2) { |
196 | //int target_size = mTrees.count()*2;
|
198 | //int target_size = mTrees.count()*2;
|
197 | //qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
|
199 | //qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
|
198 | //mTrees.reserve(qMax(target_size, 100));
|
200 | //mTrees.reserve(qMax(target_size, 100));
|
199 | if (logLevelDebug()) |
201 | if (logLevelDebug()) |
200 | qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count(); |
202 | qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count(); |
201 | mTrees.squeeze(); |
203 | mTrees.squeeze(); |
202 | }
|
204 | }
|
203 | }
|
205 | }
|
204 | }
|
206 | }
|
205 | mHasDeadTrees = false; // reset flag |
207 | mHasDeadTrees = false; // reset flag |
206 | }
|
208 | }
|
207 | 209 | ||
208 | void ResourceUnit::newYear() |
210 | void ResourceUnit::newYear() |
209 | {
|
211 | {
|
210 | mAggregatedWLA = 0.; |
212 | mAggregatedWLA = 0.; |
211 | mAggregatedLA = 0.; |
213 | mAggregatedLA = 0.; |
212 | mAggregatedLR = 0.; |
214 | mAggregatedLR = 0.; |
213 | mEffectiveArea = 0.; |
215 | mEffectiveArea = 0.; |
214 | mPixelCount = mStockedPixelCount = 0; |
216 | mPixelCount = mStockedPixelCount = 0; |
215 | snagNewYear(); |
217 | snagNewYear(); |
216 | if (mSoil) |
218 | if (mSoil) |
217 | mSoil->newYear(); |
219 | mSoil->newYear(); |
218 | // clear statistics global and per species...
|
220 | // clear statistics global and per species...
|
219 | QList<ResourceUnitSpecies*>::const_iterator i; |
221 | QList<ResourceUnitSpecies*>::const_iterator i; |
220 | QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd(); |
222 | QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd(); |
221 | mStatistics.clear(); |
223 | mStatistics.clear(); |
222 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
224 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
223 | (*i)->statisticsDead().clear(); |
225 | (*i)->statisticsDead().clear(); |
224 | (*i)->statisticsMgmt().clear(); |
226 | (*i)->statisticsMgmt().clear(); |
225 | (*i)->changeSapling().newYear(); |
227 | (*i)->changeSapling().newYear(); |
226 | }
|
228 | }
|
227 | 229 | ||
228 | }
|
230 | }
|
229 | 231 | ||
230 | /** production() is the "stand-level" part of the biomass production (3PG).
|
232 | /** production() is the "stand-level" part of the biomass production (3PG).
|
231 | - The amount of radiation intercepted by the stand is calculated
|
233 | - The amount of radiation intercepted by the stand is calculated
|
232 | - the water cycle is calculated
|
234 | - the water cycle is calculated
|
233 | - statistics for each species are cleared
|
235 | - statistics for each species are cleared
|
234 | - The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production)
|
236 | - The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production)
|
235 | see also: http://iland.boku.ac.at/individual+tree+light+availability */
|
237 | see also: http://iland.boku.ac.at/individual+tree+light+availability */
|
236 | void ResourceUnit::production() |
238 | void ResourceUnit::production() |
237 | {
|
239 | {
|
238 | 240 | ||
239 | if (mAggregatedWLA==0 || mPixelCount==0) { |
241 | if (mAggregatedWLA==0 || mPixelCount==0) { |
240 | // nothing to do...
|
242 | // nothing to do...
|
241 | return; |
243 | return; |
242 | }
|
244 | }
|
243 | 245 | ||
244 | // the pixel counters are filled during the height-grid-calculations
|
246 | // the pixel counters are filled during the height-grid-calculations
|
245 | mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m) |
247 | mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m) |
246 | 248 | ||
247 | // calculate the leaf area index (LAI)
|
249 | // calculate the leaf area index (LAI)
|
248 | double LAI = mAggregatedLA / mStockedArea; |
250 | double LAI = mAggregatedLA / mStockedArea; |
249 | // calculate the intercepted radiation fraction using the law of Beer Lambert
|
251 | // calculate the intercepted radiation fraction using the law of Beer Lambert
|
250 | const double k = Model::settings().lightExtinctionCoefficient; |
252 | const double k = Model::settings().lightExtinctionCoefficient; |
251 | double interception_fraction = 1. - exp(-k * LAI); |
253 | double interception_fraction = 1. - exp(-k * LAI); |
252 | mEffectiveArea = mStockedArea * interception_fraction; // m2 |
254 | mEffectiveArea = mStockedArea * interception_fraction; // m2 |
253 | 255 | ||
254 | // calculate the total weighted leaf area on this RU:
|
256 | // calculate the total weighted leaf area on this RU:
|
255 | mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA; // p_WLA |
257 | mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA; // p_WLA |
256 | if (mLRI_modification == 0.) |
258 | if (mLRI_modification == 0.) |
257 | qDebug() << "lri modifaction==0!"; |
259 | qDebug() << "lri modifaction==0!"; |
258 | 260 | ||
259 | if (logLevelDebug()) { |
261 | if (logLevelDebug()) { |
260 | DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3") |
262 | DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3") |
261 | .arg(LAI) |
263 | .arg(LAI) |
262 | .arg(interception_fraction) |
264 | .arg(interception_fraction) |
263 | .arg(mLRI_modification) |
265 | .arg(mLRI_modification) |
264 | .arg(mStockedArea); |
266 | .arg(mStockedArea); |
265 | ); |
267 | ); |
266 | }
|
268 | }
|
267 | 269 | ||
268 | // calculate LAI fractions
|
270 | // calculate LAI fractions
|
269 | QList<ResourceUnitSpecies*>::const_iterator i; |
271 | QList<ResourceUnitSpecies*>::const_iterator i; |
270 | QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd(); |
272 | QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd(); |
271 | double ru_lai = leafAreaIndex(); |
273 | double ru_lai = leafAreaIndex(); |
272 | if (ru_lai < 1.) |
274 | if (ru_lai < 1.) |
273 | ru_lai = 1.; |
275 | ru_lai = 1.; |
274 | // note: LAIFactors are only 1 if sum of LAI is > 1. (see WaterCycle)
|
276 | // note: LAIFactors are only 1 if sum of LAI is > 1. (see WaterCycle)
|
275 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
277 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
276 | double lai_factor = (*i)->statistics().leafAreaIndex() / ru_lai; |
278 | double lai_factor = (*i)->statistics().leafAreaIndex() / ru_lai; |
277 | DBGMODE(
|
279 | DBGMODE(
|
278 | if (lai_factor > 1.) |
280 | if (lai_factor > 1.) |
279 | qDebug() << "LAI factor > 1"; |
281 | qDebug() << "LAI factor > 1"; |
280 | ); |
282 | ); |
281 | (*i)->setLAIfactor( lai_factor ); |
283 | (*i)->setLAIfactor( lai_factor ); |
282 | }
|
284 | }
|
283 | 285 | ||
284 | // soil water model - this determines soil water contents needed for response calculations
|
286 | // soil water model - this determines soil water contents needed for response calculations
|
285 | {
|
287 | {
|
286 | mWater->run(); |
288 | mWater->run(); |
287 | }
|
289 | }
|
288 | 290 | ||
289 | // invoke species specific calculation (3PG)
|
291 | // invoke species specific calculation (3PG)
|
290 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
292 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
291 | (*i)->calculate(); // CALCULATE 3PG |
293 | (*i)->calculate(); // CALCULATE 3PG |
292 | if (logLevelInfo() && (*i)->LAIfactor()>0) |
294 | if (logLevelInfo() && (*i)->LAIfactor()>0) |
293 | qDebug() << "ru" << mIndex << "species" << (*i)->species()->id() << "LAIfraction" << (*i)->LAIfactor() << "raw_gpp_m2" |
295 | qDebug() << "ru" << mIndex << "species" << (*i)->species()->id() << "LAIfraction" << (*i)->LAIfactor() << "raw_gpp_m2" |
294 | << (*i)->prod3PG().GPPperArea() << "area:" << productiveArea() << "gpp:" |
296 | << (*i)->prod3PG().GPPperArea() << "area:" << productiveArea() << "gpp:" |
295 | << productiveArea()*(*i)->prod3PG().GPPperArea() |
297 | << productiveArea()*(*i)->prod3PG().GPPperArea() |
296 | << "aging(lastyear):" << averageAging() << "f_env,yr:" << (*i)->prod3PG().fEnvYear(); |
298 | << "aging(lastyear):" << averageAging() << "f_env,yr:" << (*i)->prod3PG().fEnvYear(); |
297 | }
|
299 | }
|
298 | }
|
300 | }
|
299 | 301 | ||
300 | void ResourceUnit::calculateInterceptedArea() |
302 | void ResourceUnit::calculateInterceptedArea() |
301 | {
|
303 | {
|
302 | if (mAggregatedLR==0) { |
304 | if (mAggregatedLR==0) { |
303 | mEffectiveArea_perWLA = 0.; |
305 | mEffectiveArea_perWLA = 0.; |
304 | return; |
306 | return; |
305 | }
|
307 | }
|
306 | Q_ASSERT(mAggregatedLR>0.); |
308 | Q_ASSERT(mAggregatedLR>0.); |
307 | mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR; |
309 | mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR; |
308 | if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA; |
310 | if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA; |
309 | }
|
311 | }
|
310 | 312 | ||
311 | // function is called immediately before the growth of individuals
|
313 | // function is called immediately before the growth of individuals
|
312 | void ResourceUnit::beforeGrow() |
314 | void ResourceUnit::beforeGrow() |
313 | {
|
315 | {
|
314 | mAverageAging = 0.; |
316 | mAverageAging = 0.; |
315 | }
|
317 | }
|
316 | 318 | ||
317 | // function is called after finishing the indivdual growth / mortality.
|
319 | // function is called after finishing the indivdual growth / mortality.
|
318 | void ResourceUnit::afterGrow() |
320 | void ResourceUnit::afterGrow() |
319 | {
|
321 | {
|
320 | mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees) |
322 | mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees) |
321 | if (mAverageAging>0. && mAverageAging<0.00001) |
323 | if (mAverageAging>0. && mAverageAging<0.00001) |
322 | qDebug() << "ru" << mIndex << "aging <0.00001"; |
324 | qDebug() << "ru" << mIndex << "aging <0.00001"; |
323 | if (mAverageAging<0. || mAverageAging>1.) |
325 | if (mAverageAging<0. || mAverageAging>1.) |
324 | qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex(); |
326 | qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex(); |
325 | }
|
327 | }
|
326 | 328 | ||
327 | void ResourceUnit::yearEnd() |
329 | void ResourceUnit::yearEnd() |
328 | {
|
330 | {
|
329 | // calculate statistics for all tree species of the ressource unit
|
331 | // calculate statistics for all tree species of the ressource unit
|
330 | int c = mRUSpecies.count(); |
332 | int c = mRUSpecies.count(); |
331 | for (int i=0;i<c; i++) { |
333 | for (int i=0;i<c; i++) { |
332 | mRUSpecies[i]->statisticsDead().calculate(); // calculate the dead trees |
334 | mRUSpecies[i]->statisticsDead().calculate(); // calculate the dead trees |
333 | mRUSpecies[i]->statisticsMgmt().calculate(); // stats of removed trees |
335 | mRUSpecies[i]->statisticsMgmt().calculate(); // stats of removed trees |
334 | mRUSpecies[i]->updateGWL(); // get sum of dead trees (died + removed) |
336 | mRUSpecies[i]->updateGWL(); // get sum of dead trees (died + removed) |
335 | mRUSpecies[i]->statistics().calculate(); // calculate the living (and add removed volume to gwl) |
337 | mRUSpecies[i]->statistics().calculate(); // calculate the living (and add removed volume to gwl) |
336 | mStatistics.add(mRUSpecies[i]->statistics()); |
338 | mStatistics.add(mRUSpecies[i]->statistics()); |
337 | }
|
339 | }
|
338 | mStatistics.calculate(); // aggreagte on stand level |
340 | mStatistics.calculate(); // aggreagte on stand level |
339 | 341 | ||
340 | }
|
342 | }
|
341 | 343 | ||
342 | void ResourceUnit::addTreeAgingForAllTrees() |
344 | void ResourceUnit::addTreeAgingForAllTrees() |
343 | {
|
345 | {
|
344 | mAverageAging = 0.; |
346 | mAverageAging = 0.; |
345 | foreach(const Tree &t, mTrees) { |
347 | foreach(const Tree &t, mTrees) { |
346 | addTreeAging(t.leafArea(), t.species()->aging(t.height(), t.age())); |
348 | addTreeAging(t.leafArea(), t.species()->aging(t.height(), t.age())); |
347 | }
|
349 | }
|
348 | 350 | ||
349 | }
|
351 | }
|
350 | 352 | ||
351 | /// refresh of tree based statistics.
|
353 | /// refresh of tree based statistics.
|
352 | /// WARNING: this function is only called once (during startup).
|
354 | /// WARNING: this function is only called once (during startup).
|
353 | /// see function "yearEnd()" above!!!
|
355 | /// see function "yearEnd()" above!!!
|
354 | void ResourceUnit::createStandStatistics() |
356 | void ResourceUnit::createStandStatistics() |
355 | {
|
357 | {
|
356 | // clear statistics (ru-level and ru-species level)
|
358 | // clear statistics (ru-level and ru-species level)
|
357 | mStatistics.clear(); |
359 | mStatistics.clear(); |
358 | for (int i=0;i<mRUSpecies.count();i++) { |
360 | for (int i=0;i<mRUSpecies.count();i++) { |
359 | mRUSpecies[i]->statistics().clear(); |
361 | mRUSpecies[i]->statistics().clear(); |
360 | mRUSpecies[i]->statisticsDead().clear(); |
362 | mRUSpecies[i]->statisticsDead().clear(); |
361 | mRUSpecies[i]->statisticsMgmt().clear(); |
363 | mRUSpecies[i]->statisticsMgmt().clear(); |
362 | }
|
364 | }
|
363 | 365 | ||
364 | // add all trees to the statistics objects of the species
|
366 | // add all trees to the statistics objects of the species
|
365 | foreach(const Tree &t, mTrees) { |
367 | foreach(const Tree &t, mTrees) { |
366 | if (!t.isDead()) |
368 | if (!t.isDead()) |
367 | resourceUnitSpecies(t.species()).statistics().add(&t, 0); |
369 | resourceUnitSpecies(t.species()).statistics().add(&t, 0); |
368 | }
|
370 | }
|
369 | // summarize statistics for the whole resource unit
|
371 | // summarize statistics for the whole resource unit
|
370 | for (int i=0;i<mRUSpecies.count();i++) { |
372 | for (int i=0;i<mRUSpecies.count();i++) { |
371 | mRUSpecies[i]->statistics().calculate(); |
373 | mRUSpecies[i]->statistics().calculate(); |
372 | mStatistics.add(mRUSpecies[i]->statistics()); |
374 | mStatistics.add(mRUSpecies[i]->statistics()); |
373 | }
|
375 | }
|
374 | mStatistics.calculate(); |
376 | mStatistics.calculate(); |
375 | mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*stockableArea()):0.; |
377 | mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*stockableArea()):0.; |
376 | if (mAverageAging<0. || mAverageAging>1.) |
378 | if (mAverageAging<0. || mAverageAging>1.) |
377 | qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex(); |
379 | qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex(); |
378 | }
|
380 | }
|
379 | 381 | ||
380 | /** recreate statistics. This is necessary after events that changed the structure
|
382 | /** recreate statistics. This is necessary after events that changed the structure
|
381 | of the stand *after* the growth of trees (where stand statistics are updated).
|
383 | of the stand *after* the growth of trees (where stand statistics are updated).
|
382 | An example is after disturbances. */
|
384 | An example is after disturbances. */
|
383 | void ResourceUnit::recreateStandStatistics() |
385 | void ResourceUnit::recreateStandStatistics() |
384 | {
|
386 | {
|
385 | for (int i=0;i<mRUSpecies.count();i++) { |
387 | for (int i=0;i<mRUSpecies.count();i++) { |
386 | mRUSpecies[i]->statistics().clear(); |
388 | mRUSpecies[i]->statistics().clear(); |
387 | }
|
389 | }
|
388 | foreach(const Tree &t, mTrees) { |
390 | foreach(const Tree &t, mTrees) { |
389 | resourceUnitSpecies(t.species()).statistics().add(&t, 0); |
391 | resourceUnitSpecies(t.species()).statistics().add(&t, 0); |
390 | }
|
392 | }
|
391 | 393 | ||
392 | }
|
394 | }
|
393 | 395 | ||
394 | void ResourceUnit::setMaxSaplingHeightAt(const QPoint &position, const float height) |
396 | void ResourceUnit::setMaxSaplingHeightAt(const QPoint &position, const float height) |
395 | {
|
397 | {
|
396 | Q_ASSERT(mSaplingHeightMap); |
398 | Q_ASSERT(mSaplingHeightMap); |
397 | int pixel_index = cPxPerRU*(position.x()-mCornerCoord.x())+(position.y()-mCornerCoord.y()); |
399 | int pixel_index = cPxPerRU*(position.x()-mCornerCoord.x())+(position.y()-mCornerCoord.y()); |
398 | if (pixel_index<0 || pixel_index>=cPxPerRU*cPxPerRU) { |
400 | if (pixel_index<0 || pixel_index>=cPxPerRU*cPxPerRU) { |
399 | qDebug() << "setSaplingHeightAt-Error for position" << position << "for RU at" << boundingBox() << "with corner" << mCornerCoord; |
401 | qDebug() << "setSaplingHeightAt-Error for position" << position << "for RU at" << boundingBox() << "with corner" << mCornerCoord; |
400 | } else { |
402 | } else { |
401 | if (mSaplingHeightMap[pixel_index]<height) |
403 | if (mSaplingHeightMap[pixel_index]<height) |
402 | mSaplingHeightMap[pixel_index]=height; |
404 | mSaplingHeightMap[pixel_index]=height; |
403 | }
|
405 | }
|
404 | }
|
406 | }
|
405 | 407 | ||
406 | /// clear all saplings of all species on a given position (after recruitment)
|
408 | /// clear all saplings of all species on a given position (after recruitment)
|
407 | void ResourceUnit::clearSaplings(const QPoint &position) |
409 | void ResourceUnit::clearSaplings(const QPoint &position) |
408 | {
|
410 | {
|
409 | foreach(ResourceUnitSpecies* rus, mRUSpecies) |
411 | foreach(ResourceUnitSpecies* rus, mRUSpecies) |
410 | rus->clearSaplings(position); |
412 | rus->clearSaplings(position); |
411 | 413 | ||
412 | }
|
414 | }
|
413 | 415 | ||
414 | /// kill all saplings within a given rect
|
416 | /// kill all saplings within a given rect
|
415 | void ResourceUnit::clearSaplings(const QRectF pixel_rect, const bool remove_from_soil) |
417 | void ResourceUnit::clearSaplings(const QRectF pixel_rect, const bool remove_from_soil) |
416 | {
|
418 | {
|
417 | foreach(ResourceUnitSpecies* rus, mRUSpecies) { |
419 | foreach(ResourceUnitSpecies* rus, mRUSpecies) { |
418 | rus->changeSapling().clearSaplings(pixel_rect, remove_from_soil); |
420 | rus->changeSapling().clearSaplings(pixel_rect, remove_from_soil); |
419 | }
|
421 | }
|
420 | 422 | ||
421 | }
|
423 | }
|
422 | 424 | ||
423 | float ResourceUnit::saplingHeightForInit(const QPoint &position) const |
425 | float ResourceUnit::saplingHeightForInit(const QPoint &position) const |
424 | {
|
426 | {
|
425 | double maxh = 0.; |
427 | double maxh = 0.; |
426 | foreach(ResourceUnitSpecies* rus, mRUSpecies) |
428 | foreach(ResourceUnitSpecies* rus, mRUSpecies) |
427 | maxh = qMax(maxh, rus->sapling().heightAt(position)); |
429 | maxh = qMax(maxh, rus->sapling().heightAt(position)); |
428 | return maxh; |
430 | return maxh; |
429 | }
|
431 | }
|
430 | 432 | ||
431 | void ResourceUnit::calculateCarbonCycle() |
433 | void ResourceUnit::calculateCarbonCycle() |
432 | {
|
434 | {
|
433 | if (!snag()) |
435 | if (!snag()) |
434 | return; |
436 | return; |
435 | 437 | ||
436 | // (1) calculate the snag dynamics
|
438 | // (1) calculate the snag dynamics
|
437 | // because all carbon/nitrogen-flows from trees to the soil are routed through the snag-layer,
|
439 | // because all carbon/nitrogen-flows from trees to the soil are routed through the snag-layer,
|
438 | // all soil inputs (litter + deadwood) are collected in the Snag-object.
|
440 | // all soil inputs (litter + deadwood) are collected in the Snag-object.
|
439 | snag()->calculateYear(); |
441 | snag()->calculateYear(); |
440 | soil()->setClimateFactor( snag()->climateFactor() ); // the climate factor is only calculated once |
442 | soil()->setClimateFactor( snag()->climateFactor() ); // the climate factor is only calculated once |
441 | soil()->setSoilInput( snag()->labileFlux(), snag()->refractoryFlux()); |
443 | soil()->setSoilInput( snag()->labileFlux(), snag()->refractoryFlux()); |
442 | soil()->calculateYear(); // update the ICBM/2N model |
444 | soil()->calculateYear(); // update the ICBM/2N model |
443 | // use available nitrogen?
|
445 | // use available nitrogen?
|
444 | if (Model::settings().useDynamicAvailableNitrogen) |
446 | if (Model::settings().useDynamicAvailableNitrogen) |
445 | mUnitVariables.nitrogenAvailable = soil()->availableNitrogen(); |
447 | mUnitVariables.nitrogenAvailable = soil()->availableNitrogen(); |
446 | 448 | ||
447 | // debug output
|
449 | // debug output
|
448 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dCarbonCycle) && !snag()->isEmpty()) { |
450 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dCarbonCycle) && !snag()->isEmpty()) { |
449 | DebugList &out = GlobalSettings::instance()->debugList(index(), GlobalSettings::dCarbonCycle); |
451 | DebugList &out = GlobalSettings::instance()->debugList(index(), GlobalSettings::dCarbonCycle); |
450 | out << index() << id(); // resource unit index and id |
452 | out << index() << id(); // resource unit index and id |
451 | out << snag()->debugList(); // snag debug outs |
453 | out << snag()->debugList(); // snag debug outs |
452 | out << soil()->debugList(); // ICBM/2N debug outs |
454 | out << soil()->debugList(); // ICBM/2N debug outs |
453 | }
|
455 | }
|
454 | 456 | ||
455 | }
|
457 | }
|
456 | 458 | ||
457 | 459 | ||
458 | 460 |