Rev 227 | Rev 251 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
117 | Werner | 2 | #include "global.h" |
3 | #include "tree.h" |
||
3 | Werner | 4 | |
83 | Werner | 5 | #include "grid.h" |
3 | Werner | 6 | |
83 | Werner | 7 | #include "stamp.h" |
90 | Werner | 8 | #include "species.h" |
189 | iland | 9 | #include "resourceunit.h" |
151 | iland | 10 | #include "model.h" |
38 | Werner | 11 | |
110 | Werner | 12 | // static varaibles |
106 | Werner | 13 | FloatGrid *Tree::mGrid = 0; |
151 | iland | 14 | HeightGrid *Tree::mHeightGrid = 0; |
40 | Werner | 15 | int Tree::m_statPrint=0; |
48 | Werner | 16 | int Tree::m_statAboveZ=0; |
105 | Werner | 17 | int Tree::m_statCreated=0; |
40 | Werner | 18 | int Tree::m_nextId=0; |
19 | |||
159 | werner | 20 | /// internal data structure which is passed between function |
21 | struct TreeGrowthData |
||
22 | { |
||
23 | double NPP; ///< total NPP |
||
24 | double NPP_stem; ///< NPP used for growth of stem (dbh,h) |
||
25 | double stress_index; ///< stress index used for mortality calculation |
||
26 | }; |
||
158 | werner | 27 | |
28 | /** get distance and direction between two points. |
||
29 | returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */ |
||
30 | float dist_and_direction(const QPointF &PStart, const QPointF &PEnd, float &rAngle) |
||
151 | iland | 31 | { |
158 | werner | 32 | float dx = PEnd.x() - PStart.x(); |
33 | float dy = PEnd.y() - PStart.y(); |
||
34 | float d = sqrt(dx*dx + dy*dy); |
||
35 | // direction: |
||
36 | rAngle = atan2(dx, dy); |
||
37 | return d; |
||
151 | iland | 38 | } |
39 | |||
158 | werner | 40 | |
110 | Werner | 41 | // lifecycle |
3 | Werner | 42 | Tree::Tree() |
43 | { |
||
149 | werner | 44 | mDbh = mHeight = 0; |
45 | mRU = 0; mSpecies = 0; |
||
169 | werner | 46 | mFlags = mAge = 0; |
149 | werner | 47 | mOpacity=mFoliageMass=mWoodyMass=mRootMass=mLeafArea=0.; |
159 | werner | 48 | mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.; |
106 | Werner | 49 | mId = m_nextId++; |
105 | Werner | 50 | m_statCreated++; |
3 | Werner | 51 | } |
38 | Werner | 52 | |
158 | werner | 53 | void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid) |
3 | Werner | 54 | { |
158 | werner | 55 | mGrid = gridToStamp; mHeightGrid = dominanceGrid; |
3 | Werner | 56 | } |
57 | |||
125 | Werner | 58 | /// dumps some core variables of a tree to a string. |
59 | QString Tree::dump() |
||
60 | { |
||
61 | QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8") |
||
159 | werner | 62 | .arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight) |
156 | werner | 63 | .arg(position().x()).arg(position().y()) |
125 | Werner | 64 | .arg(mRU->index()).arg(mLRI); |
65 | return result; |
||
66 | } |
||
3 | Werner | 67 | |
129 | Werner | 68 | void Tree::dumpList(DebugList &rTargetList) |
69 | { |
||
159 | werner | 70 | rTargetList << mId << species()->id() << mDbh << mHeight << position().x() << position().y() << mRU->index() << mLRI |
136 | Werner | 71 | << mWoodyMass << mRootMass << mFoliageMass << mLeafArea; |
129 | Werner | 72 | } |
73 | |||
38 | Werner | 74 | void Tree::setup() |
75 | { |
||
106 | Werner | 76 | if (mDbh<=0 || mHeight<=0) |
38 | Werner | 77 | return; |
78 | // check stamp |
||
159 | werner | 79 | Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL"); |
80 | mStamp = species()->stamp(mDbh, mHeight); |
||
110 | Werner | 81 | |
159 | werner | 82 | mFoliageMass = species()->biomassFoliage(mDbh); |
83 | mRootMass = species()->biomassRoot(mDbh) + mFoliageMass; // coarse root (allometry) + fine root (estimated size: foliage) |
||
84 | mWoodyMass = species()->biomassWoody(mDbh); |
||
110 | Werner | 85 | |
137 | Werner | 86 | // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg] |
159 | werner | 87 | mLeafArea = mFoliageMass * species()->specificLeafArea(); |
149 | werner | 88 | mOpacity = 1. - exp(-0.5 * mLeafArea / mStamp->crownArea()); |
137 | Werner | 89 | mNPPReserve = 2*mFoliageMass; // initial value |
90 | mDbhDelta = 0.1; // initial value: used in growth() to estimate diameter increment |
||
38 | Werner | 91 | } |
39 | Werner | 92 | |
110 | Werner | 93 | ////////////////////////////////////////////////// |
94 | //// Light functions (Pattern-stuff) |
||
95 | ////////////////////////////////////////////////// |
||
96 | |||
70 | Werner | 97 | #define NOFULLDBG |
77 | Werner | 98 | //#define NOFULLOPT |
39 | Werner | 99 | |
40 | Werner | 100 | |
158 | werner | 101 | void Tree::applyLIP() |
77 | Werner | 102 | { |
144 | Werner | 103 | if (!mStamp) |
104 | return; |
||
106 | Werner | 105 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
156 | werner | 106 | QPoint pos = mPositionIndex; |
106 | Werner | 107 | int offset = mStamp->offset(); |
77 | Werner | 108 | pos-=QPoint(offset, offset); |
109 | |||
110 | float local_dom; // height of Z* on the current position |
||
111 | int x,y; |
||
112 | float value; |
||
106 | Werner | 113 | int gr_stamp = mStamp->size(); |
77 | Werner | 114 | int grid_x, grid_y; |
115 | float *grid_value; |
||
106 | Werner | 116 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
77 | Werner | 117 | // todo: in this case we should use another algorithm!!! |
118 | return; |
||
119 | } |
||
120 | |||
121 | for (y=0;y<gr_stamp; ++y) { |
||
122 | grid_y = pos.y() + y; |
||
106 | Werner | 123 | grid_value = mGrid->ptr(pos.x(), grid_y); |
77 | Werner | 124 | for (x=0;x<gr_stamp;++x) { |
125 | // suppose there is no stamping outside |
||
126 | grid_x = pos.x() + x; |
||
127 | |||
151 | iland | 128 | local_dom = mHeightGrid->valueAtIndex(grid_x/5, grid_y/5).height; |
106 | Werner | 129 | value = (*mStamp)(x,y); // stampvalue |
149 | werner | 130 | value = 1. - value*mOpacity / local_dom; // calculated value |
77 | Werner | 131 | value = qMax(value, 0.02f); // limit value |
132 | |||
133 | *grid_value++ *= value; |
||
134 | } |
||
135 | } |
||
136 | |||
137 | m_statPrint++; // count # of stamp applications... |
||
138 | } |
||
139 | |||
155 | werner | 140 | /// helper function for gluing the edges together |
141 | /// index: index at grid |
||
142 | /// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50) |
||
143 | /// buffer: size of buffer around simulation area (in pixels) |
||
144 | int torusIndex(int index, int count, int buffer) |
||
145 | { |
||
146 | return buffer + (index-buffer+count)%count; |
||
147 | } |
||
62 | Werner | 148 | |
155 | werner | 149 | |
150 | /** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area. |
||
151 | */ |
||
158 | werner | 152 | void Tree::applyLIP_torus() |
155 | werner | 153 | { |
154 | if (!mStamp) |
||
155 | return; |
||
156 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
||
157 | |||
156 | werner | 158 | QPoint pos = mPositionIndex; |
155 | werner | 159 | int offset = mStamp->offset(); |
160 | pos-=QPoint(offset, offset); |
||
161 | |||
162 | float local_dom; // height of Z* on the current position |
||
163 | int x,y; |
||
164 | float value; |
||
165 | int gr_stamp = mStamp->size(); |
||
166 | int grid_x, grid_y; |
||
167 | float *grid_value; |
||
168 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
||
169 | // todo: in this case we should use another algorithm!!! necessary???? |
||
170 | return; |
||
171 | } |
||
172 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
||
173 | int xt, yt; // wraparound coordinates |
||
174 | for (y=0;y<gr_stamp; ++y) { |
||
175 | grid_y = pos.y() + y; |
||
176 | yt = torusIndex(grid_y, 50,bufferOffset); // 50 cells per 100m |
||
177 | for (x=0;x<gr_stamp;++x) { |
||
178 | // suppose there is no stamping outside |
||
179 | grid_x = pos.x() + x; |
||
180 | xt = torusIndex(grid_x,50,bufferOffset); |
||
181 | |||
182 | local_dom = mHeightGrid->valueAtIndex(xt/5,yt/5).height; |
||
183 | value = (*mStamp)(x,y); // stampvalue |
||
184 | value = 1. - value*mOpacity / local_dom; // calculated value |
||
185 | value = qMax(value, 0.02f); // limit value |
||
186 | |||
187 | grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates |
||
188 | *grid_value *= value; |
||
189 | } |
||
190 | } |
||
191 | |||
192 | m_statPrint++; // count # of stamp applications... |
||
193 | } |
||
194 | |||
74 | Werner | 195 | /** heightGrid() |
196 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid. |
||
197 | */ |
||
198 | void Tree::heightGrid() |
||
199 | { |
||
200 | // height of Z* |
||
106 | Werner | 201 | const float cellsize = mHeightGrid->cellsize(); |
74 | Werner | 202 | |
156 | werner | 203 | QPoint p = QPoint(mPositionIndex.x()/5, mPositionIndex.y()/5); // pos of tree on height grid |
74 | Werner | 204 | |
151 | iland | 205 | // count trees that are on height-grid cells (used for stockable area) |
206 | mHeightGrid->valueAtIndex(p).count++; |
||
207 | |||
156 | werner | 208 | int index_eastwest = mPositionIndex.x() % 5; // 4: very west, 0 east edge |
209 | int index_northsouth = mPositionIndex.y() % 5; // 4: northern edge, 0: southern edge |
||
74 | Werner | 210 | int dist[9]; |
211 | dist[3] = index_northsouth * 2 + 1; // south |
||
212 | dist[1] = index_eastwest * 2 + 1; // west |
||
213 | dist[5] = 10 - dist[3]; // north |
||
214 | dist[7] = 10 - dist[1]; // east |
||
215 | dist[8] = qMax(dist[5], dist[7]); // north-east |
||
216 | dist[6] = qMax(dist[3], dist[7]); // south-east |
||
217 | dist[0] = qMax(dist[3], dist[1]); // south-west |
||
218 | dist[2] = qMax(dist[5], dist[1]); // north-west |
||
75 | Werner | 219 | dist[4] = 0; // center cell |
76 | Werner | 220 | /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then: |
221 | index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above. |
||
222 | e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2 |
||
223 | */ |
||
74 | Werner | 224 | |
225 | |||
106 | Werner | 226 | int ringcount = int(floor(mHeight / cellsize)) + 1; |
74 | Werner | 227 | int ix, iy; |
228 | int ring; |
||
229 | QPoint pos; |
||
230 | float hdom; |
||
231 | |||
232 | for (ix=-ringcount;ix<=ringcount;ix++) |
||
233 | for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
234 | ring = qMax(abs(ix), abs(iy)); |
||
235 | QPoint pos(ix+p.x(), iy+p.y()); |
||
106 | Werner | 236 | if (mHeightGrid->isIndexValid(pos)) { |
151 | iland | 237 | float &rHGrid = mHeightGrid->valueAtIndex(pos).height; |
106 | Werner | 238 | if (rHGrid > mHeight) // skip calculation if grid is higher than tree |
74 | Werner | 239 | continue; |
240 | int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y) |
||
106 | Werner | 241 | hdom = mHeight - dist[direction]; |
74 | Werner | 242 | if (ring>1) |
243 | hdom -= (ring-1)*10; |
||
244 | |||
245 | rHGrid = qMax(rHGrid, hdom); // write value |
||
246 | } // is valid |
||
247 | } // for (y) |
||
39 | Werner | 248 | } |
40 | Werner | 249 | |
155 | werner | 250 | |
251 | |||
158 | werner | 252 | void Tree::readLIF() |
40 | Werner | 253 | { |
106 | Werner | 254 | if (!mStamp) |
155 | werner | 255 | return; |
256 | const Stamp *reader = mStamp->reader(); |
||
257 | if (!reader) |
||
258 | return; |
||
156 | werner | 259 | QPoint pos_reader = mPositionIndex; |
155 | werner | 260 | |
261 | int offset_reader = reader->offset(); |
||
262 | int offset_writer = mStamp->offset(); |
||
263 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
||
264 | |||
265 | QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer); |
||
266 | pos_reader-=QPoint(offset_reader, offset_reader); |
||
40 | Werner | 267 | QPoint p; |
268 | |||
155 | werner | 269 | //float dom_height = (*m_dominanceGrid)[m_Position]; |
270 | float local_dom; |
||
271 | |||
40 | Werner | 272 | int x,y; |
273 | double sum=0.; |
||
155 | werner | 274 | double value, own_value; |
275 | float *grid_value; |
||
276 | int reader_size = reader->size(); |
||
277 | int rx = pos_reader.x(); |
||
278 | int ry = pos_reader.y(); |
||
279 | for (y=0;y<reader_size; ++y, ++ry) { |
||
280 | grid_value = mGrid->ptr(rx, ry); |
||
281 | for (x=0;x<reader_size;++x) { |
||
282 | |||
283 | //p = pos_reader + QPoint(x,y); |
||
284 | //if (m_grid->isIndexValid(p)) { |
||
285 | local_dom = mHeightGrid->valueAtIndex((rx+x)/5, ry/5).height; // ry: gets ++ed in outer loop, rx not |
||
286 | //local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) ); |
||
287 | |||
288 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height; |
||
289 | own_value = qMax(own_value, 0.02); |
||
290 | value = *grid_value++ / own_value; // remove impact of focal tree |
||
291 | //if (value>0.) |
||
292 | sum += value * (*reader)(x,y); |
||
293 | |||
294 | //} // isIndexValid |
||
40 | Werner | 295 | } |
296 | } |
||
155 | werner | 297 | mLRI = sum; |
48 | Werner | 298 | // read dominance field... |
155 | werner | 299 | // this applies only if some trees are potentially *higher* than the dominant height grid |
300 | //if (dom_height < m_Height) { |
||
48 | Werner | 301 | // if tree is higher than Z*, the dominance height |
302 | // a part of the crown is in "full light". |
||
155 | werner | 303 | // m_statAboveZ++; |
304 | // mImpact = 1. - (1. - mImpact)*dom_height/m_Height; |
||
305 | //} |
||
306 | if (mLRI > 1.) |
||
307 | mLRI = 1.; |
||
206 | werner | 308 | |
212 | werner | 309 | mLightResponse = Model::settings().lightResponse->calculate(mLRI); |
206 | werner | 310 | // Finally, add LRI of this Tree to the ResourceUnit! |
212 | werner | 311 | mRU->addWLA(mLeafArea, mLightResponse); |
206 | werner | 312 | |
212 | werner | 313 | |
155 | werner | 314 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact; |
206 | werner | 315 | //mRU->addWLA(mLRI*mLeafArea, mLeafArea); |
40 | Werner | 316 | } |
317 | |||
158 | werner | 318 | void Tree::heightGrid_torus() |
155 | werner | 319 | { |
320 | // height of Z* |
||
321 | const float cellsize = mHeightGrid->cellsize(); |
||
58 | Werner | 322 | |
156 | werner | 323 | QPoint p = QPoint(mPositionIndex.x()/5, mPositionIndex.y()/5); // pos of tree on height grid |
155 | werner | 324 | |
325 | // count trees that are on height-grid cells (used for stockable area) |
||
326 | mHeightGrid->valueAtIndex(p).count++; |
||
327 | |||
156 | werner | 328 | int index_eastwest = mPositionIndex.x() % 5; // 4: very west, 0 east edge |
329 | int index_northsouth = mPositionIndex.y() % 5; // 4: northern edge, 0: southern edge |
||
155 | werner | 330 | int dist[9]; |
331 | dist[3] = index_northsouth * 2 + 1; // south |
||
332 | dist[1] = index_eastwest * 2 + 1; // west |
||
333 | dist[5] = 10 - dist[3]; // north |
||
334 | dist[7] = 10 - dist[1]; // east |
||
335 | dist[8] = qMax(dist[5], dist[7]); // north-east |
||
336 | dist[6] = qMax(dist[3], dist[7]); // south-east |
||
337 | dist[0] = qMax(dist[3], dist[1]); // south-west |
||
338 | dist[2] = qMax(dist[5], dist[1]); // north-west |
||
339 | dist[4] = 0; // center cell |
||
340 | /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then: |
||
341 | index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above. |
||
342 | e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2 |
||
343 | */ |
||
344 | |||
345 | |||
346 | int ringcount = int(floor(mHeight / cellsize)) + 1; |
||
347 | int ix, iy; |
||
348 | int ring; |
||
349 | QPoint pos; |
||
350 | float hdom; |
||
351 | int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
||
352 | for (ix=-ringcount;ix<=ringcount;ix++) |
||
353 | for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
354 | ring = qMax(abs(ix), abs(iy)); |
||
355 | QPoint pos(ix+p.x(), iy+p.y()); |
||
356 | if (mHeightGrid->isIndexValid(pos)) { |
||
357 | float &rHGrid = mHeightGrid->valueAtIndex(torusIndex(pos.x(),10,bufferOffset), torusIndex(pos.y(),10,bufferOffset)).height; |
||
358 | if (rHGrid > mHeight) // skip calculation if grid is higher than tree |
||
359 | continue; |
||
360 | int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y) |
||
361 | hdom = mHeight - dist[direction]; |
||
362 | if (ring>1) |
||
363 | hdom -= (ring-1)*10; |
||
364 | |||
365 | rHGrid = qMax(rHGrid, hdom); // write value |
||
366 | } // is valid |
||
367 | } // for (y) |
||
368 | } |
||
369 | |||
370 | /// Torus version of read stamp (glued edges) |
||
158 | werner | 371 | void Tree::readLIF_torus() |
78 | Werner | 372 | { |
106 | Werner | 373 | if (!mStamp) |
107 | Werner | 374 | return; |
106 | Werner | 375 | const Stamp *reader = mStamp->reader(); |
78 | Werner | 376 | if (!reader) |
107 | Werner | 377 | return; |
156 | werner | 378 | QPoint pos_reader = mPositionIndex; |
78 | Werner | 379 | |
380 | int offset_reader = reader->offset(); |
||
106 | Werner | 381 | int offset_writer = mStamp->offset(); |
78 | Werner | 382 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
383 | |||
384 | QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer); |
||
385 | pos_reader-=QPoint(offset_reader, offset_reader); |
||
386 | QPoint p; |
||
387 | |||
388 | //float dom_height = (*m_dominanceGrid)[m_Position]; |
||
389 | float local_dom; |
||
390 | |||
391 | int x,y; |
||
392 | double sum=0.; |
||
393 | double value, own_value; |
||
394 | float *grid_value; |
||
395 | int reader_size = reader->size(); |
||
396 | int rx = pos_reader.x(); |
||
397 | int ry = pos_reader.y(); |
||
155 | werner | 398 | int xt, yt; // wrapped coords |
399 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
||
400 | |||
78 | Werner | 401 | for (y=0;y<reader_size; ++y, ++ry) { |
106 | Werner | 402 | grid_value = mGrid->ptr(rx, ry); |
78 | Werner | 403 | for (x=0;x<reader_size;++x) { |
155 | werner | 404 | xt = torusIndex(rx+x,50, bufferOffset); |
405 | yt = torusIndex(ry+y,50, bufferOffset); |
||
406 | grid_value = mGrid->ptr(xt,yt); |
||
78 | Werner | 407 | //p = pos_reader + QPoint(x,y); |
408 | //if (m_grid->isIndexValid(p)) { |
||
155 | werner | 409 | local_dom = mHeightGrid->valueAtIndex(xt/5, yt/5).height; // ry: gets ++ed in outer loop, rx not |
78 | Werner | 410 | //local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) ); |
125 | Werner | 411 | |
149 | werner | 412 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height; |
78 | Werner | 413 | own_value = qMax(own_value, 0.02); |
155 | werner | 414 | value = *grid_value / own_value; // remove impact of focal tree |
78 | Werner | 415 | //if (value>0.) |
416 | sum += value * (*reader)(x,y); |
||
417 | |||
418 | //} // isIndexValid |
||
419 | } |
||
420 | } |
||
106 | Werner | 421 | mLRI = sum; |
78 | Werner | 422 | // read dominance field... |
423 | // this applies only if some trees are potentially *higher* than the dominant height grid |
||
424 | //if (dom_height < m_Height) { |
||
425 | // if tree is higher than Z*, the dominance height |
||
426 | // a part of the crown is in "full light". |
||
427 | // m_statAboveZ++; |
||
428 | // mImpact = 1. - (1. - mImpact)*dom_height/m_Height; |
||
429 | //} |
||
148 | iland | 430 | if (mLRI > 1.) |
431 | mLRI = 1.; |
||
78 | Werner | 432 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact; |
205 | werner | 433 | |
434 | // calculate a light response from lri: |
||
212 | werner | 435 | mLightResponse = Model::settings().lightResponse->calculate(mLRI); |
205 | werner | 436 | // Finally, add LRI of this Tree to the ResourceUnit! |
212 | werner | 437 | mRU->addWLA(mLeafArea, mLightResponse); |
58 | Werner | 438 | } |
439 | |||
155 | werner | 440 | |
40 | Werner | 441 | void Tree::resetStatistics() |
442 | { |
||
443 | m_statPrint=0; |
||
105 | Werner | 444 | m_statCreated=0; |
48 | Werner | 445 | m_statAboveZ=0; |
40 | Werner | 446 | m_nextId=1; |
447 | } |
||
107 | Werner | 448 | |
110 | Werner | 449 | ////////////////////////////////////////////////// |
450 | //// Growth Functions |
||
451 | ////////////////////////////////////////////////// |
||
107 | Werner | 452 | |
227 | werner | 453 | /** grow() is the main function of the yearly tree growth. |
454 | The main steps are: |
||
455 | - Production of GPP/NPP @sa http://iland.boku.ac.at/primary+production |
||
456 | - Partitioning of NPP to biomass compartments of the tree @sa http://iland.boku.ac.at/allocation (???) |
||
457 | - Growth of the stem http://iland.boku.ac.at/stem+growth (???) |
||
458 | Additionally, the age of the tree is increased and the mortality sub routine is executed.*/ |
||
107 | Werner | 459 | void Tree::grow() |
460 | { |
||
159 | werner | 461 | TreeGrowthData d; |
169 | werner | 462 | mAge++; // increase age |
230 | werner | 463 | // step 1: get "interception area" of the tree individual [m2] |
464 | // the sum of all area of all trees of a unit equal the total stocked area * interception_factor(Beer-Lambert) |
||
465 | double effective_area = mRU->interceptedArea(mLeafArea, mLightResponse); |
||
107 | Werner | 466 | |
230 | werner | 467 | // step 2: calculate GPP of the tree based |
468 | // (1) get the amount of GPP for a "unit area" of the tree species |
||
469 | double raw_gpp_per_area = mRU->resourceUnitSpecies(species()).prod3PG().GPPperArea(); |
||
470 | // (2) GPP (without aging-effect) in kg Biomass / year |
||
471 | double raw_gpp = raw_gpp_per_area * effective_area; |
||
161 | werner | 472 | |
227 | werner | 473 | // apply aging according to the state of the individuum |
169 | werner | 474 | const double aging_factor = mSpecies->aging(mHeight, mAge); |
227 | werner | 475 | double gpp = raw_gpp * aging_factor; // |
476 | d.NPP = gpp * 0.47; // respiration loss, cf. Waring et al 1998. |
||
113 | Werner | 477 | |
133 | Werner | 478 | DBGMODE( |
137 | Werner | 479 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) { |
133 | Werner | 480 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP); |
481 | dumpList(out); // add tree headers |
||
230 | werner | 482 | out << effective_area << raw_gpp << gpp << d.NPP << aging_factor; |
133 | Werner | 483 | } |
484 | ); // DBGMODE() |
||
217 | werner | 485 | if (d.NPP>0.) |
486 | partitioning(d); // split npp to compartments and grow (diameter, height) |
||
133 | Werner | 487 | |
200 | werner | 488 | if (Model::settings().mortalityEnabled) |
489 | mortality(d); |
||
110 | Werner | 490 | |
159 | werner | 491 | mStressIndex = d.stress_index; |
180 | werner | 492 | |
493 | if (!isDead()) |
||
200 | werner | 494 | mRU->resourceUnitSpecies(species()).statistics().add(this); |
107 | Werner | 495 | } |
496 | |||
227 | werner | 497 | /** partitioning of this years assimilates (NPP) to biomass compartments. |
498 | Conceptionally, the algorithm is based on Duursma, 2007. */ |
||
159 | werner | 499 | inline void Tree::partitioning(TreeGrowthData &d) |
115 | Werner | 500 | { |
164 | werner | 501 | if (isDebugging()) |
502 | enableDebugging(true); |
||
159 | werner | 503 | double npp = d.NPP; |
115 | Werner | 504 | // add content of reserve pool |
116 | Werner | 505 | npp += mNPPReserve; |
159 | werner | 506 | const double foliage_mass_allo = species()->biomassFoliage(mDbh); |
136 | Werner | 507 | const double reserve_size = 2 * foliage_mass_allo; |
163 | werner | 508 | double refill_reserve = qMin(reserve_size, 2.*mFoliageMass); // not always try to refill reserve 100% |
119 | Werner | 509 | |
136 | Werner | 510 | double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta) |
117 | Werner | 511 | // turnover rates |
159 | werner | 512 | const double &to_fol = species()->turnoverLeaf(); |
513 | const double &to_root = species()->turnoverRoot(); |
||
136 | Werner | 514 | // the turnover rate of wood depends on the size of the reserve pool: |
116 | Werner | 515 | |
136 | Werner | 516 | |
163 | werner | 517 | double to_wood = refill_reserve / (mWoodyMass + refill_reserve); |
518 | |||
227 | werner | 519 | apct_root = mRU->resourceUnitSpecies(species()).prod3PG().rootFraction(); |
159 | werner | 520 | double b_wf = species()->allometricRatio_wf(); // ratio of allometric exponents... now fixed |
117 | Werner | 521 | |
522 | // Duursma 2007, Eq. (20) |
||
167 | werner | 523 | apct_wood = (foliage_mass_allo*to_wood/npp + b_wf*(1.-apct_root) - b_wf*foliage_mass_allo*to_fol/npp) / ( foliage_mass_allo/mWoodyMass + b_wf ); |
163 | werner | 524 | if (apct_wood<0) |
525 | apct_wood = 0.; |
||
117 | Werner | 526 | apct_foliage = 1. - apct_root - apct_wood; |
527 | |||
163 | werner | 528 | |
529 | //DBGMODE( |
||
530 | if (apct_foliage<0 || apct_wood<0) |
||
531 | qDebug() << "transfer to foliage or wood < 0"; |
||
532 | if (npp<0) |
||
533 | qDebug() << "NPP < 0"; |
||
534 | // ); |
||
535 | |||
136 | Werner | 536 | // Change of biomass compartments |
159 | werner | 537 | double sen_root = mRootMass*to_root; |
538 | double sen_foliage = mFoliageMass*to_fol; |
||
136 | Werner | 539 | // Roots |
159 | werner | 540 | double delta_root = apct_root * npp - sen_root; |
137 | Werner | 541 | mRootMass += delta_root; |
119 | Werner | 542 | |
136 | Werner | 543 | // Foliage |
159 | werner | 544 | double delta_foliage = apct_foliage * npp - sen_foliage; |
137 | Werner | 545 | mFoliageMass += delta_foliage; |
217 | werner | 546 | if (_isnan(mFoliageMass)) |
547 | qDebug() << "foliage mass invalid!"; |
||
163 | werner | 548 | if (mFoliageMass<0.) mFoliageMass=0.; // limit to zero |
549 | |||
159 | werner | 550 | mLeafArea = mFoliageMass * species()->specificLeafArea(); // update leaf area |
119 | Werner | 551 | |
198 | werner | 552 | // stress index: different varaints at denominatior: to_fol*foliage_mass = leafmass to rebuild, |
553 | // foliage_mass_allo: simply higher chance for stress |
||
159 | werner | 554 | |
207 | werner | 555 | d.stress_index =qMax(1. - (npp) / ( to_fol*foliage_mass_allo + reserve_size), 0.); |
198 | werner | 556 | |
136 | Werner | 557 | // Woody compartments |
558 | // (1) transfer to reserve pool |
||
559 | double gross_woody = apct_wood * npp; |
||
560 | double to_reserve = qMin(reserve_size, gross_woody); |
||
561 | mNPPReserve = to_reserve; |
||
562 | double net_woody = gross_woody - to_reserve; |
||
137 | Werner | 563 | double net_stem = 0.; |
164 | werner | 564 | mDbhDelta = 0.; |
165 | werner | 565 | |
566 | |||
136 | Werner | 567 | if (net_woody > 0.) { |
568 | // (2) calculate part of increment that is dedicated to the stem (which is a function of diameter) |
||
159 | werner | 569 | net_stem = net_woody * species()->allometricFractionStem(mDbh); |
570 | d.NPP_stem = net_stem; |
||
137 | Werner | 571 | mWoodyMass += net_woody; |
136 | Werner | 572 | // (3) growth of diameter and height baseed on net stem increment |
159 | werner | 573 | grow_diameter(d); |
136 | Werner | 574 | } |
119 | Werner | 575 | |
129 | Werner | 576 | DBGMODE( |
137 | Werner | 577 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreePartition) |
578 | && isDebugging() ) { |
||
129 | Werner | 579 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreePartition); |
580 | dumpList(out); // add tree headers |
||
136 | Werner | 581 | out << npp << apct_foliage << apct_wood << apct_root |
137 | Werner | 582 | << delta_foliage << net_woody << delta_root << mNPPReserve << net_stem; |
583 | } |
||
144 | Werner | 584 | |
129 | Werner | 585 | ); // DBGMODE() |
144 | Werner | 586 | //DBGMODE( |
587 | if (mWoodyMass<0. || mWoodyMass>10000 || mFoliageMass<0. || mFoliageMass>1000. || mRootMass<0. || mRootMass>10000 |
||
588 | || mNPPReserve>2000.) { |
||
589 | qDebug() << "Tree:partitioning: invalid pools."; |
||
590 | qDebug() << GlobalSettings::instance()->debugListCaptions(GlobalSettings::DebugOutputs(0)); |
||
591 | DebugList dbg; dumpList(dbg); |
||
592 | qDebug() << dbg; |
||
593 | } //); |
||
594 | |||
136 | Werner | 595 | /*DBG_IF_X(mId == 1 , "Tree::partitioning", "dump", dump() |
596 | + QString("npp %1 npp_reserve %9 sen_fol %2 sen_stem %3 sen_root %4 net_fol %5 net_stem %6 net_root %7 to_reserve %8") |
||
597 | .arg(npp).arg(senescence_foliage).arg(senescence_stem).arg(senescence_root) |
||
598 | .arg(net_foliage).arg(net_stem).arg(net_root).arg(to_reserve).arg(mNPPReserve) );*/ |
||
129 | Werner | 599 | |
115 | Werner | 600 | } |
601 | |||
125 | Werner | 602 | |
134 | Werner | 603 | /** Determination of diamter and height growth based on increment of the stem mass (@p net_stem_npp). |
125 | Werner | 604 | Refer to XXX for equations and variables. |
605 | This function updates the dbh and height of the tree. |
||
227 | werner | 606 | The equations are based on dbh in meters! */ |
159 | werner | 607 | inline void Tree::grow_diameter(TreeGrowthData &d) |
119 | Werner | 608 | { |
609 | // determine dh-ratio of increment |
||
610 | // height increment is a function of light competition: |
||
125 | Werner | 611 | double hd_growth = relative_height_growth(); // hd of height growth |
153 | werner | 612 | double d_m = mDbh / 100.; // current diameter in [m] |
159 | werner | 613 | double net_stem_npp = d.NPP_stem; |
614 | |||
153 | werner | 615 | const double d_delta_m = mDbhDelta / 100.; // increment of last year in [m] |
115 | Werner | 616 | |
159 | werner | 617 | const double mass_factor = species()->volumeFactor() * species()->density(); |
153 | werner | 618 | double stem_mass = mass_factor * d_m*d_m * mHeight; // result: kg, dbh[cm], h[meter] |
123 | Werner | 619 | |
153 | werner | 620 | // factor is in diameter increment per NPP [m/kg] |
621 | double factor_diameter = 1. / ( mass_factor * (d_m + d_delta_m)*(d_m + d_delta_m) * ( 2. * mHeight/d_m + hd_growth) ); |
||
125 | Werner | 622 | double delta_d_estimate = factor_diameter * net_stem_npp; // estimated dbh-inc using last years increment |
623 | |||
624 | // using that dbh-increment we estimate a stem-mass-increment and the residual (Eq. 9) |
||
153 | werner | 625 | double stem_estimate = mass_factor * (d_m + delta_d_estimate)*(d_m + delta_d_estimate)*(mHeight + delta_d_estimate*hd_growth); |
137 | Werner | 626 | double stem_residual = stem_estimate - (stem_mass + net_stem_npp); |
125 | Werner | 627 | |
628 | // the final increment is then: |
||
629 | double d_increment = factor_diameter * (net_stem_npp - stem_residual); // Eq. (11) |
||
144 | Werner | 630 | DBG_IF_X(d_increment<0. || d_increment>0.1, "Tree::grow_dimater", "increment out of range.", dump() |
125 | Werner | 631 | + QString("\nhdz %1 factor_diameter %2 stem_residual %3 delta_d_estimate %4 d_increment %5 final residual(kg) %6") |
632 | .arg(hd_growth).arg(factor_diameter).arg(stem_residual).arg(delta_d_estimate).arg(d_increment) |
||
142 | Werner | 633 | .arg( mass_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)) )); |
125 | Werner | 634 | |
635 | DBGMODE( |
||
153 | werner | 636 | double res_final = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)); |
125 | Werner | 637 | DBG_IF_X(res_final > 1, "Tree::grow_diameter", "final residual stem estimate > 1kg", dump()); |
153 | werner | 638 | DBG_IF_X(d_increment > 10. || d_increment*hd_growth >10., "Tree::grow_diameter", "growth out of bound:",QString("d-increment %1 h-increment %2 ").arg(d_increment).arg(d_increment*hd_growth/100.) + dump()); |
158 | werner | 639 | |
137 | Werner | 640 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeGrowth) && isDebugging() ) { |
126 | Werner | 641 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeGrowth); |
129 | Werner | 642 | dumpList(out); // add tree headers |
143 | Werner | 643 | out << net_stem_npp << stem_mass << hd_growth << factor_diameter << delta_d_estimate*100 << d_increment*100; |
126 | Werner | 644 | } |
153 | werner | 645 | |
142 | Werner | 646 | ); // DBGMODE() |
125 | Werner | 647 | |
648 | d_increment = qMax(d_increment, 0.); |
||
649 | |||
650 | // update state variables |
||
153 | werner | 651 | mDbh += d_increment*100; // convert from [m] to [cm] |
652 | mDbhDelta = d_increment*100; // save for next year's growth |
||
653 | mHeight += d_increment * hd_growth; |
||
158 | werner | 654 | |
655 | // update state of LIP stamp and opacity |
||
159 | werner | 656 | mStamp = species()->stamp(mDbh, mHeight); // get new stamp for updated dimensions |
158 | werner | 657 | // calculate the CrownFactor which reflects the opacity of the crown |
200 | werner | 658 | const double k=Model::settings().lightExtinctionCoefficientOpacity; |
659 | mOpacity = 1. - exp(-k * mLeafArea / mStamp->crownArea()); |
||
158 | werner | 660 | |
119 | Werner | 661 | } |
662 | |||
125 | Werner | 663 | |
664 | /// return the HD ratio of this year's increment based on the light status. |
||
119 | Werner | 665 | inline double Tree::relative_height_growth() |
666 | { |
||
667 | double hd_low, hd_high; |
||
668 | mSpecies->hdRange(mDbh, hd_low, hd_high); |
||
669 | |||
125 | Werner | 670 | DBG_IF_X(hd_low>hd_high, "Tree::relative_height_growth", "hd low higher dann hd_high for ", dump()); |
671 | DBG_IF_X(hd_low < 20 || hd_high>250, "Tree::relative_height_growth", "hd out of range ", dump() + QString(" hd-low: %1 hd-high: %2").arg(hd_low).arg(hd_high)); |
||
672 | |||
673 | // scale according to LRI: if receiving much light (LRI=1), the result is hd_low (for open grown trees) |
||
674 | double hd_ratio = hd_high - (hd_high-hd_low)*mLRI; |
||
675 | return hd_ratio; |
||
119 | Werner | 676 | } |
141 | Werner | 677 | |
159 | werner | 678 | void Tree::mortality(TreeGrowthData &d) |
679 | { |
||
163 | werner | 680 | // death if leaf area is 0 |
681 | if (mFoliageMass<0.00001) |
||
682 | die(); |
||
683 | |||
159 | werner | 684 | double p_death, p_stress; |
170 | werner | 685 | p_stress = d.stress_index * species()->deathProb_stress(); |
686 | //if (d.stress_index>0) |
||
687 | // p_stress = species()->deathProb_stress(); |
||
159 | werner | 688 | p_death = species()->deathProb_intrinsic() + p_stress; |
190 | werner | 689 | double p = drandom(); //0..1 |
159 | werner | 690 | if (p<p_death) { |
691 | // die... |
||
692 | die(); |
||
693 | } |
||
694 | } |
||
141 | Werner | 695 | |
696 | ////////////////////////////////////////////////// |
||
697 | //// value functions |
||
698 | ////////////////////////////////////////////////// |
||
699 | |||
145 | Werner | 700 | double Tree::volume() const |
141 | Werner | 701 | { |
702 | /// @see Species::volumeFactor() for details |
||
159 | werner | 703 | const double volume_factor = species()->volumeFactor(); |
157 | werner | 704 | const double volume = volume_factor * mDbh*mDbh*mHeight * 0.0001; // dbh in cm: cm/100 * cm/100 = cm*cm * 0.0001 = m2 |
141 | Werner | 705 | return volume; |
706 | } |
||
180 | werner | 707 | |
708 | double Tree::basalArea() const |
||
709 | { |
||
710 | // A = r^2 * pi = d/2*pi; from cm->m: d/200 |
||
711 | const double b = (mDbh/200.)*(mDbh/200.)*M_PI; |
||
712 | return b; |
||
713 | } |