Rev 376 | Rev 436 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
189 | iland | 2 | /** @class ResourceUnit |
3 | ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components |
||
92 | Werner | 4 | (Climate, Soil, Water, ...). |
5 | |||
6 | */ |
||
7 | #include <QtCore> |
||
8 | #include "global.h" |
||
9 | |||
189 | iland | 10 | #include "resourceunit.h" |
229 | werner | 11 | #include "resourceunitspecies.h" |
111 | Werner | 12 | #include "speciesset.h" |
13 | #include "species.h" |
||
113 | Werner | 14 | #include "production3pg.h" |
200 | werner | 15 | #include "model.h" |
208 | werner | 16 | #include "climate.h" |
241 | werner | 17 | #include "watercycle.h" |
18 | #include "helper.h" |
||
92 | Werner | 19 | |
241 | werner | 20 | ResourceUnit::~ResourceUnit() |
21 | { |
||
22 | if (mWater) |
||
23 | delete mWater; |
||
24 | } |
||
111 | Werner | 25 | |
189 | iland | 26 | ResourceUnit::ResourceUnit(const int index) |
92 | Werner | 27 | { |
94 | Werner | 28 | mSpeciesSet = 0; |
208 | werner | 29 | mClimate = 0; |
331 | werner | 30 | mPixelCount=0; |
113 | Werner | 31 | mIndex = index; |
241 | werner | 32 | mWater = new WaterCycle(); |
33 | |||
157 | werner | 34 | mTrees.reserve(100); // start with space for 100 trees. |
92 | Werner | 35 | } |
105 | Werner | 36 | |
241 | werner | 37 | void ResourceUnit::setup() |
38 | { |
||
39 | mWater->setup(this); |
||
281 | werner | 40 | // setup variables |
41 | mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40); |
||
376 | werner | 42 | mAverageAging = 0.; |
281 | werner | 43 | |
241 | werner | 44 | } |
45 | |||
111 | Werner | 46 | /// set species and setup the species-per-RU-data |
189 | iland | 47 | void ResourceUnit::setSpeciesSet(SpeciesSet *set) |
111 | Werner | 48 | { |
49 | mSpeciesSet = set; |
||
50 | mRUSpecies.clear(); |
||
229 | werner | 51 | mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated |
111 | Werner | 52 | for (int i=0;i<set->count();i++) { |
53 | Species *s = const_cast<Species*>(mSpeciesSet->species(i)); |
||
54 | if (!s) |
||
189 | iland | 55 | throw IException("ResourceUnit::setSpeciesSet: invalid index!"); |
229 | werner | 56 | |
57 | /* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container. |
||
58 | If the container memory is relocated (QVector), the pointer gets invalid!!! |
||
59 | Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */ |
||
60 | mRUSpecies[i].setup(s,this); // setup this element |
||
277 | werner | 61 | |
111 | Werner | 62 | } |
63 | } |
||
64 | |||
200 | werner | 65 | ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species) |
111 | Werner | 66 | { |
67 | return mRUSpecies[species->index()]; |
||
68 | } |
||
69 | |||
189 | iland | 70 | Tree &ResourceUnit::newTree() |
105 | Werner | 71 | { |
72 | // start simple: just append to the vector... |
||
73 | mTrees.append(Tree()); |
||
74 | return mTrees.back(); |
||
75 | } |
||
287 | werner | 76 | int ResourceUnit::newTreeIndex() |
77 | { |
||
78 | // start simple: just append to the vector... |
||
79 | mTrees.append(Tree()); |
||
80 | return mTrees.count()-1; |
||
81 | } |
||
107 | Werner | 82 | |
157 | werner | 83 | /// remove dead trees from tree list |
84 | /// reduce size of vector if lots of space is free |
||
85 | /// tests showed that this way of cleanup is very fast, |
||
86 | /// because no memory allocations are performed (simple memmove()) |
||
87 | /// when trees are moved. |
||
189 | iland | 88 | void ResourceUnit::cleanTreeList() |
157 | werner | 89 | { |
90 | QVector<Tree>::iterator last=mTrees.end()-1; |
||
91 | QVector<Tree>::iterator current = mTrees.begin(); |
||
158 | werner | 92 | while (last>=current && (*last).isDead()) |
157 | werner | 93 | --last; |
107 | Werner | 94 | |
157 | werner | 95 | while (current<last) { |
158 | werner | 96 | if ((*current).isDead()) { |
157 | werner | 97 | *current = *last; // copy data! |
98 | --last; // |
||
158 | werner | 99 | while (last>=current && (*last).isDead()) |
157 | werner | 100 | --last; |
101 | } |
||
102 | ++current; |
||
103 | } |
||
104 | ++last; // last points now to the first dead tree |
||
105 | |||
106 | // free ressources |
||
278 | werner | 107 | if (last!=mTrees.end()) { |
108 | mTrees.erase(last, mTrees.end()); |
||
109 | if (mTrees.capacity()>100) { |
||
110 | if (mTrees.count() / double(mTrees.capacity()) < 0.2) { |
||
111 | int target_size = mTrees.count()*2; |
||
112 | qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size; |
||
113 | mTrees.reserve(qMax(target_size, 100)); |
||
114 | } |
||
157 | werner | 115 | } |
116 | } |
||
117 | } |
||
118 | |||
189 | iland | 119 | void ResourceUnit::newYear() |
107 | Werner | 120 | { |
251 | werner | 121 | mAggregatedWLA = 0.; |
122 | mAggregatedLA = 0.; |
||
123 | mAggregatedLR = 0.; |
||
124 | mEffectiveArea = 0.; |
||
151 | iland | 125 | mPixelCount = mStockedPixelCount = 0; |
111 | Werner | 126 | // clear statistics global and per species... |
278 | werner | 127 | ResourceUnitSpecies *i; |
128 | QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end(); |
||
129 | mStatistics.clear(); |
||
130 | for (i=mRUSpecies.begin(); i!=iend; ++i) { |
||
131 | i->statisticsDead().clear(); |
||
132 | i->statisticsMgmt().clear(); |
||
133 | } |
||
134 | |||
107 | Werner | 135 | } |
110 | Werner | 136 | |
112 | Werner | 137 | /** production() is the "stand-level" part of the biomass production (3PG). |
138 | - The amount of radiation intercepted by the stand is calculated |
||
331 | werner | 139 | - the water cycle is calculated |
140 | - statistics for each species are cleared |
||
141 | - The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production) |
||
298 | werner | 142 | see also: http://iland.boku.ac.at/individual+tree+light+availability */ |
189 | iland | 143 | void ResourceUnit::production() |
110 | Werner | 144 | { |
241 | werner | 145 | |
151 | iland | 146 | if (mAggregatedWLA==0 || mPixelCount==0) { |
112 | Werner | 147 | // nothing to do... |
148 | return; |
||
149 | } |
||
151 | iland | 150 | |
151 | // the pixel counters are filled during the height-grid-calculations |
||
230 | werner | 152 | mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m) |
153 | |||
112 | Werner | 154 | // calculate the leaf area index (LAI) |
151 | iland | 155 | double LAI = mAggregatedLA / mStockedArea; |
112 | Werner | 156 | // calculate the intercepted radiation fraction using the law of Beer Lambert |
200 | werner | 157 | const double k = Model::settings().lightExtinctionCoefficient; |
112 | Werner | 158 | double interception_fraction = 1. - exp(-k * LAI); |
251 | werner | 159 | mEffectiveArea = mStockedArea * interception_fraction; // m2 |
112 | Werner | 160 | |
230 | werner | 161 | // calculate the total weighted leaf area on this RU: |
251 | werner | 162 | mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA; |
265 | werner | 163 | if (mLRI_modification == 0.) |
164 | qDebug() << "lri modifaction==0!"; |
||
205 | werner | 165 | |
251 | werner | 166 | |
167 | DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3") |
||
230 | werner | 168 | .arg(LAI) |
169 | .arg(interception_fraction) |
||
251 | werner | 170 | .arg(mLRI_modification) |
230 | werner | 171 | .arg(mStockedArea); |
172 | ); |
||
367 | werner | 173 | |
174 | // calculate LAI fractions |
||
175 | ResourceUnitSpecies *i; |
||
176 | QVector<ResourceUnitSpecies>::iterator iend = mRUSpecies.end(); |
||
177 | for (i=mRUSpecies.begin(); i!=iend; ++i) { |
||
178 | i->setLAIfactor(i->statistics().leafAreaIndex() / leafAreaIndex()); |
||
179 | } |
||
180 | |||
241 | werner | 181 | // soil water model - this determines soil water contents needed for response calculations |
182 | { |
||
183 | DebugTimer tw("water:run"); |
||
184 | mWater->run(); |
||
185 | } |
||
112 | Werner | 186 | |
187 | // invoke species specific calculation (3PG) |
||
188 | for (i=mRUSpecies.begin(); i!=iend; ++i) { |
||
331 | werner | 189 | i->statistics().clear(); |
229 | werner | 190 | i->calculate(); |
431 | werner | 191 | if (logLevelInfo() && i->LAIfactor()>0) |
376 | werner | 192 | qDebug() << "ru" << mIndex << "species" << (*i).species()->id() << "LAIfraction" << i->LAIfactor() << "raw_gpp_m2" |
193 | << i->prod3PG().GPPperArea() << "area:" << productiveArea() << "gpp:" |
||
194 | << productiveArea()*i->prod3PG().GPPperArea() |
||
195 | << "aging(lastyear):" << averageAging(); |
||
112 | Werner | 196 | } |
110 | Werner | 197 | } |
198 | |||
251 | werner | 199 | void ResourceUnit::calculateInterceptedArea() |
200 | { |
||
265 | werner | 201 | if (mAggregatedLR==0) { |
202 | mEffectiveArea_perWLA = 0.; |
||
203 | return; |
||
204 | } |
||
251 | werner | 205 | Q_ASSERT(mAggregatedLR>0.); |
206 | mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR; |
||
431 | werner | 207 | if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA; |
251 | werner | 208 | } |
209 | |||
376 | werner | 210 | // function is called immediately before the growth of individuals |
211 | void ResourceUnit::beforeGrow() |
||
212 | { |
||
213 | mAverageAging = 0.; |
||
214 | } |
||
215 | |||
216 | // function is called after finishing the indivdual growth / mortality. |
||
217 | void ResourceUnit::afterGrow() |
||
218 | { |
||
219 | mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees) |
||
220 | if (mAverageAging>0. && mAverageAging<0.00001) |
||
221 | qDebug() << "ru" << mIndex << "aging <0.00001"; |
||
222 | } |
||
223 | |||
189 | iland | 224 | void ResourceUnit::yearEnd() |
180 | werner | 225 | { |
226 | // calculate statistics for all tree species of the ressource unit |
||
227 | int c = mRUSpecies.count(); |
||
228 | for (int i=0;i<c; i++) { |
||
277 | werner | 229 | mRUSpecies[i].statisticsDead().calculate(); // calculate the dead trees |
278 | werner | 230 | mRUSpecies[i].statisticsMgmt().calculate(); // stats of removed trees |
231 | mRUSpecies[i].updateGWL(); // get sum of dead trees (died + removed) |
||
277 | werner | 232 | mRUSpecies[i].statistics().calculate(); // calculate the living (and add removed volume to gwl) |
180 | werner | 233 | mStatistics.add(mRUSpecies[i].statistics()); |
234 | } |
||
235 | mStatistics.calculate(); // aggreagte on stand level |
||
236 | } |
||
237 | |||
241 | werner | 238 | /// refresh of tree based statistics. |
240 | werner | 239 | void ResourceUnit::createStandStatistics() |
240 | { |
||
241 | werner | 241 | // clear statistics (ru-level and ru-species level) |
240 | werner | 242 | mStatistics.clear(); |
262 | werner | 243 | for (int i=0;i<mRUSpecies.count();i++) { |
240 | werner | 244 | mRUSpecies[i].statistics().clear(); |
262 | werner | 245 | mRUSpecies[i].statisticsDead().clear(); |
278 | werner | 246 | mRUSpecies[i].statisticsMgmt().clear(); |
262 | werner | 247 | } |
241 | werner | 248 | |
249 | // add all trees to the statistics objects of the species |
||
240 | werner | 250 | foreach(const Tree &t, mTrees) { |
251 | if (!t.isDead()) |
||
257 | werner | 252 | resourceUnitSpecies(t.species()).statistics().add(&t, 0); |
240 | werner | 253 | } |
241 | werner | 254 | // summarize statistics for the whole resource unit |
240 | werner | 255 | for (int i=0;i<mRUSpecies.count();i++) { |
256 | mRUSpecies[i].statistics().calculate(); |
||
257 | mStatistics.add(mRUSpecies[i].statistics()); |
||
258 | } |
||
331 | werner | 259 | mStatistics.calculate(); |
376 | werner | 260 | mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*area()):0.; |
240 | werner | 261 | } |