Rev 1202 | Rev 1217 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
671 | werner | 2 | /******************************************************************************************** |
3 | ** iLand - an individual based forest landscape and disturbance model |
||
4 | ** http://iland.boku.ac.at |
||
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl |
||
6 | ** |
||
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 |
||
9 | ** the Free Software Foundation, either version 3 of the License, or |
||
10 | ** (at your option) any later version. |
||
11 | ** |
||
12 | ** This program is distributed in the hope that it will be useful, |
||
13 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of |
||
14 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
||
15 | ** GNU General Public License for more details. |
||
16 | ** |
||
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/>. |
||
19 | ********************************************************************************************/ |
||
20 | |||
534 | werner | 21 | /** @class ResourceUnit |
22 | ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components |
||
23 | (Climate, Soil, Water, ...). |
||
697 | werner | 24 | @ingroup core |
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). |
||
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). |
||
534 | werner | 29 | |
30 | */ |
||
31 | #include <QtCore> |
||
32 | #include "global.h" |
||
33 | |||
34 | #include "resourceunit.h" |
||
35 | #include "resourceunitspecies.h" |
||
36 | #include "speciesset.h" |
||
37 | #include "species.h" |
||
38 | #include "production3pg.h" |
||
39 | #include "model.h" |
||
40 | #include "climate.h" |
||
41 | #include "watercycle.h" |
||
42 | #include "snag.h" |
||
43 | #include "soil.h" |
||
44 | #include "helper.h" |
||
45 | |||
46 | ResourceUnit::~ResourceUnit() |
||
47 | { |
||
48 | if (mWater) |
||
49 | delete mWater; |
||
50 | mWater = 0; |
||
51 | if (mSnag) |
||
52 | delete mSnag; |
||
53 | if (mSoil) |
||
54 | delete mSoil; |
||
55 | |||
738 | werner | 56 | qDeleteAll(mRUSpecies); |
57 | |||
1159 | werner | 58 | if (mSaplings) |
59 | delete[] mSaplings; |
||
60 | |||
534 | werner | 61 | mSnag = 0; |
62 | mSoil = 0; |
||
1159 | werner | 63 | mSaplings = 0; |
534 | werner | 64 | } |
65 | |||
66 | ResourceUnit::ResourceUnit(const int index) |
||
67 | { |
||
68 | qDeleteAll(mRUSpecies); |
||
69 | mSpeciesSet = 0; |
||
70 | mClimate = 0; |
||
71 | mPixelCount=0; |
||
72 | mStockedArea = 0; |
||
73 | mStockedPixelCount = 0; |
||
1157 | werner | 74 | mStockableArea = 0; |
1024 | werner | 75 | mAggregatedWLA = 0.; |
76 | mAggregatedLA = 0.; |
||
77 | mAggregatedLR = 0.; |
||
78 | mEffectiveArea = 0.; |
||
79 | mLRI_modification = 0.; |
||
534 | werner | 80 | mIndex = index; |
81 | mSaplingHeightMap = 0; |
||
82 | mEffectiveArea_perWLA = 0.; |
||
83 | mWater = new WaterCycle(); |
||
84 | mSnag = 0; |
||
85 | mSoil = 0; |
||
1159 | werner | 86 | mSaplings = 0; |
569 | werner | 87 | mID = 0; |
534 | werner | 88 | } |
89 | |||
90 | void ResourceUnit::setup() |
||
91 | { |
||
92 | mWater->setup(this); |
||
93 | |||
94 | if (mSnag) |
||
95 | delete mSnag; |
||
96 | mSnag=0; |
||
97 | if (mSoil) |
||
98 | delete mSoil; |
||
99 | mSoil=0; |
||
100 | if (Model::settings().carbonCycleEnabled) { |
||
591 | werner | 101 | mSoil = new Soil(this); |
534 | werner | 102 | mSnag = new Snag; |
103 | mSnag->setup(this); |
||
104 | const XmlHelper &xml=GlobalSettings::instance()->settings(); |
||
105 | |||
106 | // setup contents of the soil of the RU; use values for C and N (kg/ha) |
||
107 | mSoil->setInitialState(CNPool(xml.valueDouble("model.site.youngLabileC", -1), |
||
108 | xml.valueDouble("model.site.youngLabileN", -1), |
||
109 | xml.valueDouble("model.site.youngLabileDecompRate", -1)), |
||
110 | CNPool(xml.valueDouble("model.site.youngRefractoryC", -1), |
||
111 | xml.valueDouble("model.site.youngRefractoryN", -1), |
||
112 | xml.valueDouble("model.site.youngRefractoryDecompRate", -1)), |
||
113 | CNPair(xml.valueDouble("model.site.somC", -1), xml.valueDouble("model.site.somN", -1))); |
||
114 | } |
||
115 | |||
1159 | werner | 116 | if (mSaplings) |
117 | delete mSaplings; |
||
118 | if (Model::settings().regenerationEnabled) { |
||
119 | mSaplings = new SaplingCell[cPxPerHectare]; |
||
120 | } |
||
121 | |||
534 | werner | 122 | // setup variables |
123 | mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40); |
||
124 | |||
895 | werner | 125 | // if dynamic coupling of soil nitrogen is enabled, a starting value for available N is calculated |
534 | werner | 126 | if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) { |
127 | mSoil->setClimateFactor(1.); |
||
128 | mSoil->calculateYear(); |
||
895 | werner | 129 | mUnitVariables.nitrogenAvailable = soil()->availableNitrogen(); |
534 | werner | 130 | } |
664 | werner | 131 | mHasDeadTrees = false; |
534 | werner | 132 | mAverageAging = 0.; |
133 | |||
134 | } |
||
135 | void ResourceUnit::setBoundingBox(const QRectF &bb) |
||
136 | { |
||
137 | mBoundingBox = bb; |
||
1118 | werner | 138 | mCornerOffset = GlobalSettings::instance()->model()->grid()->indexAt(bb.topLeft()); |
534 | werner | 139 | } |
140 | |||
1203 | werner | 141 | /// return the sapling cell at given LIF-coordinates |
142 | SaplingCell *ResourceUnit::saplingCell(const QPoint &lifCoords) const |
||
143 | { |
||
144 | // LIF-Coordinates are global, we here need (RU-)local coordinates |
||
145 | int ix = lifCoords.x() % cPxPerRU; |
||
146 | int iy = lifCoords.y() % cPxPerRU; |
||
147 | int i = iy*cPxPerRU+ix; |
||
148 | Q_ASSERT(i>=0 && i<cPxPerHectare); |
||
149 | return &mSaplings[i]; |
||
150 | } |
||
151 | |||
534 | werner | 152 | /// set species and setup the species-per-RU-data |
153 | void ResourceUnit::setSpeciesSet(SpeciesSet *set) |
||
154 | { |
||
155 | mSpeciesSet = set; |
||
156 | qDeleteAll(mRUSpecies); |
||
157 | |||
158 | //mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated |
||
159 | for (int i=0;i<set->count();i++) { |
||
160 | Species *s = const_cast<Species*>(mSpeciesSet->species(i)); |
||
161 | if (!s) |
||
162 | throw IException("ResourceUnit::setSpeciesSet: invalid index!"); |
||
163 | |||
164 | ResourceUnitSpecies *rus = new ResourceUnitSpecies(); |
||
165 | mRUSpecies.push_back(rus); |
||
166 | rus->setup(s, this); |
||
167 | /* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container. |
||
168 | If the container memory is relocated (QVector), the pointer gets invalid!!! |
||
169 | Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */ |
||
170 | //mRUSpecies[i].setup(s,this); // setup this element |
||
171 | |||
172 | } |
||
173 | } |
||
174 | |||
175 | ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species) |
||
176 | { |
||
177 | return *mRUSpecies[species->index()]; |
||
178 | } |
||
179 | |||
1040 | werner | 180 | const ResourceUnitSpecies *ResourceUnit::constResourceUnitSpecies(const Species *species) const |
181 | { |
||
182 | return mRUSpecies[species->index()]; |
||
183 | } |
||
184 | |||
534 | werner | 185 | Tree &ResourceUnit::newTree() |
186 | { |
||
187 | // start simple: just append to the vector... |
||
188 | if (mTrees.isEmpty()) |
||
189 | mTrees.reserve(100); // reserve a junk of memory for trees |
||
190 | |||
191 | mTrees.append(Tree()); |
||
192 | return mTrees.back(); |
||
193 | } |
||
194 | int ResourceUnit::newTreeIndex() |
||
195 | { |
||
734 | werner | 196 | newTree(); |
197 | return mTrees.count()-1; // return index of the last tree |
||
534 | werner | 198 | } |
199 | |||
200 | /// remove dead trees from tree list |
||
201 | /// reduce size of vector if lots of space is free |
||
202 | /// tests showed that this way of cleanup is very fast, |
||
203 | /// because no memory allocations are performed (simple memmove()) |
||
204 | /// when trees are moved. |
||
205 | void ResourceUnit::cleanTreeList() |
||
206 | { |
||
664 | werner | 207 | if (!mHasDeadTrees) |
208 | return; |
||
209 | |||
534 | werner | 210 | QVector<Tree>::iterator last=mTrees.end()-1; |
211 | QVector<Tree>::iterator current = mTrees.begin(); |
||
212 | while (last>=current && (*last).isDead()) |
||
213 | --last; |
||
214 | |||
215 | while (current<last) { |
||
216 | if ((*current).isDead()) { |
||
217 | *current = *last; // copy data! |
||
218 | --last; // |
||
219 | while (last>=current && (*last).isDead()) |
||
220 | --last; |
||
221 | } |
||
222 | ++current; |
||
223 | } |
||
224 | ++last; // last points now to the first dead tree |
||
225 | |||
226 | // free ressources |
||
227 | if (last!=mTrees.end()) { |
||
228 | mTrees.erase(last, mTrees.end()); |
||
229 | if (mTrees.capacity()>100) { |
||
230 | if (mTrees.count() / double(mTrees.capacity()) < 0.2) { |
||
231 | //int target_size = mTrees.count()*2; |
||
232 | //qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size; |
||
233 | //mTrees.reserve(qMax(target_size, 100)); |
||
664 | werner | 234 | if (logLevelDebug()) |
235 | qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count(); |
||
534 | werner | 236 | mTrees.squeeze(); |
237 | } |
||
238 | } |
||
239 | } |
||
664 | werner | 240 | mHasDeadTrees = false; // reset flag |
534 | werner | 241 | } |
242 | |||
243 | void ResourceUnit::newYear() |
||
244 | { |
||
245 | mAggregatedWLA = 0.; |
||
246 | mAggregatedLA = 0.; |
||
247 | mAggregatedLR = 0.; |
||
248 | mEffectiveArea = 0.; |
||
249 | mPixelCount = mStockedPixelCount = 0; |
||
250 | snagNewYear(); |
||
609 | werner | 251 | if (mSoil) |
252 | mSoil->newYear(); |
||
534 | werner | 253 | // clear statistics global and per species... |
254 | QList<ResourceUnitSpecies*>::const_iterator i; |
||
255 | QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd(); |
||
256 | mStatistics.clear(); |
||
257 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
||
258 | (*i)->statisticsDead().clear(); |
||
259 | (*i)->statisticsMgmt().clear(); |
||
260 | } |
||
261 | |||
262 | } |
||
263 | |||
264 | /** production() is the "stand-level" part of the biomass production (3PG). |
||
265 | - The amount of radiation intercepted by the stand is calculated |
||
266 | - the water cycle is calculated |
||
267 | - statistics for each species are cleared |
||
268 | - The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production) |
||
269 | see also: http://iland.boku.ac.at/individual+tree+light+availability */ |
||
270 | void ResourceUnit::production() |
||
271 | { |
||
272 | |||
1107 | werner | 273 | if (mAggregatedWLA==0. || mPixelCount==0) { |
936 | werner | 274 | // clear statistics of resourceunitspecies |
275 | for ( QList<ResourceUnitSpecies*>::const_iterator i=mRUSpecies.constBegin(); i!=mRUSpecies.constEnd(); ++i) |
||
276 | (*i)->statistics().clear(); |
||
277 | mEffectiveArea = 0.; |
||
278 | mStockedArea = 0.; |
||
534 | werner | 279 | return; |
280 | } |
||
281 | |||
282 | // the pixel counters are filled during the height-grid-calculations |
||
1184 | werner | 283 | mStockedArea = cHeightPerRU*cHeightPerRU * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m) |
1107 | werner | 284 | if (leafAreaIndex()<3.) { |
285 | // estimate stocked area based on crown projections |
||
286 | double crown_area = 0.; |
||
287 | for (int i=0;i<mTrees.count();++i) |
||
288 | crown_area += mTrees.at(i).isDead() ? 0. : mTrees.at(i).stamp()->reader()->crownArea(); |
||
534 | werner | 289 | |
1157 | werner | 290 | if (logLevelDebug()) |
291 | qDebug() << "crown area: lai" << leafAreaIndex() << "stocked area (pixels)" << mStockedArea << " area (crown)" << crown_area; |
||
292 | if (leafAreaIndex()<1.) { |
||
293 | mStockedArea = std::min(crown_area, mStockedArea); |
||
1107 | werner | 294 | } else { |
1184 | werner | 295 | // for LAI between 1 and 3: |
296 | // interpolate between sum of crown area of trees (at LAI=1) and the pixel-based value (at LAI=3 and above) |
||
1157 | werner | 297 | double px_frac = (leafAreaIndex()-1.)/2.; // 0 at LAI=1, 1 at LAI=3 |
298 | mStockedArea = mStockedArea * px_frac + std::min(crown_area, mStockedArea) * (1. - px_frac); |
||
1107 | werner | 299 | } |
300 | if (mStockedArea==0.) |
||
301 | return; |
||
302 | } |
||
303 | |||
534 | werner | 304 | // calculate the leaf area index (LAI) |
305 | double LAI = mAggregatedLA / mStockedArea; |
||
306 | // calculate the intercepted radiation fraction using the law of Beer Lambert |
||
307 | const double k = Model::settings().lightExtinctionCoefficient; |
||
308 | double interception_fraction = 1. - exp(-k * LAI); |
||
309 | mEffectiveArea = mStockedArea * interception_fraction; // m2 |
||
310 | |||
311 | // calculate the total weighted leaf area on this RU: |
||
312 | mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA; // p_WLA |
||
313 | if (mLRI_modification == 0.) |
||
314 | qDebug() << "lri modifaction==0!"; |
||
315 | |||
611 | werner | 316 | if (logLevelDebug()) { |
534 | werner | 317 | DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3") |
318 | .arg(LAI) |
||
319 | .arg(interception_fraction) |
||
320 | .arg(mLRI_modification) |
||
321 | .arg(mStockedArea); |
||
322 | ); |
||
611 | werner | 323 | } |
534 | werner | 324 | |
325 | // calculate LAI fractions |
||
326 | QList<ResourceUnitSpecies*>::const_iterator i; |
||
327 | QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd(); |
||
328 | double ru_lai = leafAreaIndex(); |
||
329 | if (ru_lai < 1.) |
||
330 | ru_lai = 1.; |
||
331 | // note: LAIFactors are only 1 if sum of LAI is > 1. (see WaterCycle) |
||
332 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
||
720 | werner | 333 | double lai_factor = (*i)->statistics().leafAreaIndex() / ru_lai; |
1157 | werner | 334 | |
335 | //DBGMODE( |
||
336 | if (lai_factor > 1.) { |
||
337 | const ResourceUnitSpecies* rus=*i; |
||
338 | qDebug() << "LAI factor > 1: species ru-index:" << rus->species()->name() << rus->ru()->index(); |
||
339 | } |
||
340 | //); |
||
720 | werner | 341 | (*i)->setLAIfactor( lai_factor ); |
534 | werner | 342 | } |
343 | |||
344 | // soil water model - this determines soil water contents needed for response calculations |
||
345 | { |
||
346 | mWater->run(); |
||
347 | } |
||
348 | |||
349 | // invoke species specific calculation (3PG) |
||
350 | for (i=mRUSpecies.constBegin(); i!=iend; ++i) { |
||
1157 | werner | 351 | //DBGMODE( |
352 | if ((*i)->LAIfactor() > 1.) { |
||
353 | const ResourceUnitSpecies* rus=*i; |
||
354 | qDebug() << "LAI factor > 1: species ru-index value:" << rus->species()->name() << rus->ru()->index() << rus->LAIfactor(); |
||
355 | } |
||
356 | //); |
||
534 | werner | 357 | (*i)->calculate(); // CALCULATE 3PG |
1196 | werner | 358 | |
359 | // debug output related to production |
||
360 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dStandGPP) && (*i)->LAIfactor()>0.) { |
||
361 | DebugList &out = GlobalSettings::instance()->debugList(index(), GlobalSettings::dStandGPP); |
||
362 | out << (*i)->species()->id() << index() << id(); |
||
363 | out << (*i)->LAIfactor() << (*i)->prod3PG().GPPperArea() << productiveArea()*(*i)->LAIfactor()*(*i)->prod3PG().GPPperArea() << averageAging() << (*i)->prod3PG().fEnvYear() ; |
||
364 | |||
365 | } |
||
534 | werner | 366 | } |
367 | } |
||
368 | |||
369 | void ResourceUnit::calculateInterceptedArea() |
||
370 | { |
||
371 | if (mAggregatedLR==0) { |
||
372 | mEffectiveArea_perWLA = 0.; |
||
373 | return; |
||
374 | } |
||
375 | Q_ASSERT(mAggregatedLR>0.); |
||
376 | mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR; |
||
377 | if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA; |
||
378 | } |
||
379 | |||
380 | // function is called immediately before the growth of individuals |
||
381 | void ResourceUnit::beforeGrow() |
||
382 | { |
||
383 | mAverageAging = 0.; |
||
384 | } |
||
385 | |||
386 | // function is called after finishing the indivdual growth / mortality. |
||
387 | void ResourceUnit::afterGrow() |
||
388 | { |
||
389 | mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees) |
||
390 | if (mAverageAging>0. && mAverageAging<0.00001) |
||
391 | qDebug() << "ru" << mIndex << "aging <0.00001"; |
||
392 | if (mAverageAging<0. || mAverageAging>1.) |
||
393 | qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex(); |
||
394 | } |
||
395 | |||
396 | void ResourceUnit::yearEnd() |
||
397 | { |
||
398 | // calculate statistics for all tree species of the ressource unit |
||
399 | int c = mRUSpecies.count(); |
||
400 | for (int i=0;i<c; i++) { |
||
401 | mRUSpecies[i]->statisticsDead().calculate(); // calculate the dead trees |
||
402 | mRUSpecies[i]->statisticsMgmt().calculate(); // stats of removed trees |
||
403 | mRUSpecies[i]->updateGWL(); // get sum of dead trees (died + removed) |
||
404 | mRUSpecies[i]->statistics().calculate(); // calculate the living (and add removed volume to gwl) |
||
405 | mStatistics.add(mRUSpecies[i]->statistics()); |
||
406 | } |
||
407 | mStatistics.calculate(); // aggreagte on stand level |
||
408 | |||
1157 | werner | 409 | // update carbon flows |
410 | if (soil() && GlobalSettings::instance()->model()->settings().carbonCycleEnabled) { |
||
411 | double area_factor = stockableArea() / cRUArea; //conversion factor |
||
412 | mUnitVariables.carbonUptake = statistics().npp() * biomassCFraction; |
||
413 | mUnitVariables.carbonUptake += statistics().nppSaplings() * biomassCFraction; |
||
414 | |||
415 | double to_atm = snag()->fluxToAtmosphere().C / area_factor; // from snags, kgC/ha |
||
416 | to_atm += soil()->fluxToAtmosphere().C *cRUArea/10.; // soil: t/ha -> t/m2 -> kg/ha |
||
417 | mUnitVariables.carbonToAtm = to_atm; |
||
418 | |||
419 | double to_dist = snag()->fluxToDisturbance().C / area_factor; |
||
420 | to_dist += soil()->fluxToDisturbance().C * cRUArea/10.; |
||
421 | double to_harvest = snag()->fluxToExtern().C / area_factor; |
||
422 | |||
423 | mUnitVariables.NEP = mUnitVariables.carbonUptake - to_atm - to_dist - to_harvest; // kgC/ha |
||
424 | |||
425 | // incremental values.... |
||
426 | mUnitVariables.cumCarbonUptake += mUnitVariables.carbonUptake; |
||
427 | mUnitVariables.cumCarbonToAtm += mUnitVariables.carbonToAtm; |
||
428 | mUnitVariables.cumNEP += mUnitVariables.NEP; |
||
429 | |||
430 | } |
||
431 | |||
534 | werner | 432 | } |
433 | |||
434 | void ResourceUnit::addTreeAgingForAllTrees() |
||
435 | { |
||
436 | mAverageAging = 0.; |
||
437 | foreach(const Tree &t, mTrees) { |
||
438 | addTreeAging(t.leafArea(), t.species()->aging(t.height(), t.age())); |
||
439 | } |
||
440 | |||
441 | } |
||
442 | |||
443 | /// refresh of tree based statistics. |
||
444 | /// WARNING: this function is only called once (during startup). |
||
445 | /// see function "yearEnd()" above!!! |
||
446 | void ResourceUnit::createStandStatistics() |
||
447 | { |
||
448 | // clear statistics (ru-level and ru-species level) |
||
449 | mStatistics.clear(); |
||
450 | for (int i=0;i<mRUSpecies.count();i++) { |
||
451 | mRUSpecies[i]->statistics().clear(); |
||
452 | mRUSpecies[i]->statisticsDead().clear(); |
||
453 | mRUSpecies[i]->statisticsMgmt().clear(); |
||
1178 | werner | 454 | mRUSpecies[i]->saplingStat().clearStatistics(); |
534 | werner | 455 | } |
456 | |||
457 | // add all trees to the statistics objects of the species |
||
458 | foreach(const Tree &t, mTrees) { |
||
459 | if (!t.isDead()) |
||
460 | resourceUnitSpecies(t.species()).statistics().add(&t, 0); |
||
461 | } |
||
1178 | werner | 462 | // summarise sapling stats |
463 | GlobalSettings::instance()->model()->saplings()->calculateInitialStatistics(this); |
||
464 | |||
534 | werner | 465 | // summarize statistics for the whole resource unit |
466 | for (int i=0;i<mRUSpecies.count();i++) { |
||
1178 | werner | 467 | mRUSpecies[i]->saplingStat().calculate(mRUSpecies[i]->species(), this); |
468 | mRUSpecies[i]->statistics().add(&mRUSpecies[i]->saplingStat()); |
||
534 | werner | 469 | mRUSpecies[i]->statistics().calculate(); |
470 | mStatistics.add(mRUSpecies[i]->statistics()); |
||
471 | } |
||
472 | mStatistics.calculate(); |
||
575 | werner | 473 | mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*stockableArea()):0.; |
534 | werner | 474 | if (mAverageAging<0. || mAverageAging>1.) |
475 | qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex(); |
||
1178 | werner | 476 | |
534 | werner | 477 | } |
478 | |||
720 | werner | 479 | /** recreate statistics. This is necessary after events that changed the structure |
480 | of the stand *after* the growth of trees (where stand statistics are updated). |
||
481 | An example is after disturbances. */ |
||
1157 | werner | 482 | void ResourceUnit::recreateStandStatistics(bool recalculate_stats) |
720 | werner | 483 | { |
1202 | werner | 484 | // when called after disturbances (recalculate_stats=false), we |
485 | // clear only the tree-specific variables in the stats (i.e. we keep NPP, and regen carbon), |
||
486 | // and then re-add all trees (since TreeGrowthData is NULL no NPP is available). |
||
487 | // The statistics are not summarised here, because this happens for all resource units |
||
488 | // in the yearEnd function of RU. |
||
720 | werner | 489 | for (int i=0;i<mRUSpecies.count();i++) { |
1202 | werner | 490 | if (recalculate_stats) |
491 | mRUSpecies[i]->statistics().clear(); |
||
492 | else |
||
493 | mRUSpecies[i]->statistics().clearOnlyTrees(); |
||
720 | werner | 494 | } |
495 | foreach(const Tree &t, mTrees) { |
||
496 | resourceUnitSpecies(t.species()).statistics().add(&t, 0); |
||
497 | } |
||
1157 | werner | 498 | |
499 | if (recalculate_stats) { |
||
500 | for (int i=0;i<mRUSpecies.count();i++) { |
||
501 | mRUSpecies[i]->statistics().calculate(); |
||
502 | } |
||
937 | werner | 503 | } |
720 | werner | 504 | } |
505 | |||
824 | werner | 506 | |
534 | werner | 507 | |
508 | |||
509 | void ResourceUnit::calculateCarbonCycle() |
||
510 | { |
||
511 | if (!snag()) |
||
512 | return; |
||
513 | |||
514 | // (1) calculate the snag dynamics |
||
515 | // because all carbon/nitrogen-flows from trees to the soil are routed through the snag-layer, |
||
516 | // all soil inputs (litter + deadwood) are collected in the Snag-object. |
||
517 | snag()->calculateYear(); |
||
518 | soil()->setClimateFactor( snag()->climateFactor() ); // the climate factor is only calculated once |
||
519 | soil()->setSoilInput( snag()->labileFlux(), snag()->refractoryFlux()); |
||
520 | soil()->calculateYear(); // update the ICBM/2N model |
||
521 | // use available nitrogen? |
||
522 | if (Model::settings().useDynamicAvailableNitrogen) |
||
523 | mUnitVariables.nitrogenAvailable = soil()->availableNitrogen(); |
||
524 | |||
525 | // debug output |
||
526 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dCarbonCycle) && !snag()->isEmpty()) { |
||
527 | DebugList &out = GlobalSettings::instance()->debugList(index(), GlobalSettings::dCarbonCycle); |
||
605 | werner | 528 | out << index() << id(); // resource unit index and id |
534 | werner | 529 | out << snag()->debugList(); // snag debug outs |
530 | out << soil()->debugList(); // ICBM/2N debug outs |
||
531 | } |
||
532 | |||
533 | } |
||
600 | werner | 534 | |
535 |