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