Rev 1180 | Rev 1204 | 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 | ********************************************************************************************/ |
||
117 | Werner | 20 | #include "global.h" |
21 | #include "tree.h" |
||
3 | Werner | 22 | |
83 | Werner | 23 | #include "grid.h" |
3 | Werner | 24 | |
83 | Werner | 25 | #include "stamp.h" |
90 | Werner | 26 | #include "species.h" |
189 | iland | 27 | #include "resourceunit.h" |
151 | iland | 28 | #include "model.h" |
468 | werner | 29 | #include "snag.h" |
38 | Werner | 30 | |
904 | werner | 31 | #include "forestmanagementengine.h" |
1044 | werner | 32 | #include "modules.h" |
904 | werner | 33 | |
1102 | werner | 34 | #include "treeout.h" |
1114 | werner | 35 | #include "landscapeout.h" |
1071 | werner | 36 | |
110 | Werner | 37 | // static varaibles |
106 | Werner | 38 | FloatGrid *Tree::mGrid = 0; |
151 | iland | 39 | HeightGrid *Tree::mHeightGrid = 0; |
1102 | werner | 40 | TreeRemovedOut *Tree::mRemovalOutput = 0; |
1114 | werner | 41 | LandscapeRemovedOut *Tree::mLSRemovalOutput = 0; |
40 | Werner | 42 | int Tree::m_statPrint=0; |
48 | Werner | 43 | int Tree::m_statAboveZ=0; |
105 | Werner | 44 | int Tree::m_statCreated=0; |
40 | Werner | 45 | int Tree::m_nextId=0; |
46 | |||
1071 | werner | 47 | #ifdef ALT_TREE_MORTALITY |
1102 | werner | 48 | static double _stress_threshold = 0.05; |
49 | static int _stress_years = 5; |
||
50 | static double _stress_death_prob = 0.33; |
||
1071 | werner | 51 | #endif |
52 | |||
697 | werner | 53 | /** @class Tree |
54 | @ingroup core |
||
55 | A tree is the basic simulation entity of iLand and represents a single tree. |
||
56 | Trees in iLand are designed to be lightweight, thus the list of stored properties is limited. Basic properties |
||
57 | are dimensions (dbh, height), biomass pools (stem, leaves, roots), the reserve NPP pool. Additionally, the location and species are stored. |
||
58 | A Tree has a height of at least 4m; trees below this threshold are covered by the regeneration layer (see Sapling). |
||
59 | Trees are stored in lists managed at the resource unit level. |
||
158 | werner | 60 | |
697 | werner | 61 | */ |
257 | werner | 62 | |
158 | werner | 63 | /** get distance and direction between two points. |
64 | returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */ |
||
1102 | werner | 65 | double dist_and_direction(const QPointF &PStart, const QPointF &PEnd, double &rAngle) |
151 | iland | 66 | { |
1102 | werner | 67 | double dx = PEnd.x() - PStart.x(); |
68 | double dy = PEnd.y() - PStart.y(); |
||
69 | double d = sqrt(dx*dx + dy*dy); |
||
158 | werner | 70 | // direction: |
71 | rAngle = atan2(dx, dy); |
||
72 | return d; |
||
151 | iland | 73 | } |
74 | |||
158 | werner | 75 | |
110 | Werner | 76 | // lifecycle |
3 | Werner | 77 | Tree::Tree() |
78 | { |
||
149 | werner | 79 | mDbh = mHeight = 0; |
80 | mRU = 0; mSpecies = 0; |
||
169 | werner | 81 | mFlags = mAge = 0; |
276 | werner | 82 | mOpacity=mFoliageMass=mWoodyMass=mCoarseRootMass=mFineRootMass=mLeafArea=0.; |
159 | werner | 83 | mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.; |
264 | werner | 84 | mLightResponse = 0.; |
975 | werner | 85 | mStamp=0; |
106 | Werner | 86 | mId = m_nextId++; |
105 | Werner | 87 | m_statCreated++; |
3 | Werner | 88 | } |
38 | Werner | 89 | |
407 | werner | 90 | float Tree::crownRadius() const |
91 | { |
||
92 | Q_ASSERT(mStamp!=0); |
||
93 | return mStamp->crownRadius(); |
||
94 | } |
||
95 | |||
476 | werner | 96 | float Tree::biomassBranch() const |
97 | { |
||
1102 | werner | 98 | return static_cast<float>( mSpecies->biomassBranch(mDbh) ); |
476 | werner | 99 | } |
100 | |||
158 | werner | 101 | void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid) |
3 | Werner | 102 | { |
158 | werner | 103 | mGrid = gridToStamp; mHeightGrid = dominanceGrid; |
3 | Werner | 104 | } |
105 | |||
667 | werner | 106 | // calculate the thickness of the bark of the tree |
107 | double Tree::barkThickness() const |
||
108 | { |
||
109 | return mSpecies->barkThickness(mDbh); |
||
110 | } |
||
111 | |||
125 | Werner | 112 | /// dumps some core variables of a tree to a string. |
113 | QString Tree::dump() |
||
114 | { |
||
115 | QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8") |
||
159 | werner | 116 | .arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight) |
156 | werner | 117 | .arg(position().x()).arg(position().y()) |
125 | Werner | 118 | .arg(mRU->index()).arg(mLRI); |
119 | return result; |
||
120 | } |
||
3 | Werner | 121 | |
129 | Werner | 122 | void Tree::dumpList(DebugList &rTargetList) |
123 | { |
||
159 | werner | 124 | rTargetList << mId << species()->id() << mDbh << mHeight << position().x() << position().y() << mRU->index() << mLRI |
276 | werner | 125 | << mWoodyMass << mCoarseRootMass << mFoliageMass << mLeafArea; |
129 | Werner | 126 | } |
127 | |||
38 | Werner | 128 | void Tree::setup() |
129 | { |
||
975 | werner | 130 | if (mDbh<=0 || mHeight<=0) { |
131 | throw IException(QString("Error: trying to set up a tree with invalid dimensions: dbh: %1 height: %2 id: %3 RU-index: %4").arg(mDbh).arg(mHeight).arg(mId).arg(mRU->index())); |
||
132 | } |
||
38 | Werner | 133 | // check stamp |
159 | werner | 134 | Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL"); |
135 | mStamp = species()->stamp(mDbh, mHeight); |
||
505 | werner | 136 | if (!mStamp) { |
137 | throw IException("Tree::setup() with invalid stamp!"); |
||
138 | } |
||
110 | Werner | 139 | |
1102 | werner | 140 | mFoliageMass = static_cast<float>( species()->biomassFoliage(mDbh) ); |
141 | mCoarseRootMass = static_cast<float>( species()->biomassRoot(mDbh) ); // coarse root (allometry) |
||
142 | mFineRootMass = static_cast<float>( mFoliageMass * species()->finerootFoliageRatio() ); // fine root (size defined by finerootFoliageRatio) |
||
143 | mWoodyMass = static_cast<float>( species()->biomassWoody(mDbh) ); |
||
110 | Werner | 144 | |
137 | Werner | 145 | // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg] |
1102 | werner | 146 | mLeafArea = static_cast<float>( mFoliageMass * species()->specificLeafArea() ); |
147 | mOpacity = static_cast<float>( 1. - exp(- Model::settings().lightExtinctionCoefficientOpacity * mLeafArea / mStamp->crownArea()) ); |
||
148 | mNPPReserve = static_cast<float>( (1+species()->finerootFoliageRatio())*mFoliageMass ); // initial value |
||
780 | werner | 149 | mDbhDelta = 0.1f; // initial value: used in growth() to estimate diameter increment |
376 | werner | 150 | |
38 | Werner | 151 | } |
39 | Werner | 152 | |
388 | werner | 153 | void Tree::setAge(const int age, const float treeheight) |
154 | { |
||
155 | mAge = age; |
||
156 | if (age==0) { |
||
157 | // estimate age using the tree height |
||
158 | mAge = mSpecies->estimateAge(treeheight); |
||
159 | } |
||
160 | } |
||
161 | |||
110 | Werner | 162 | ////////////////////////////////////////////////// |
163 | //// Light functions (Pattern-stuff) |
||
164 | ////////////////////////////////////////////////// |
||
165 | |||
70 | Werner | 166 | #define NOFULLDBG |
77 | Werner | 167 | //#define NOFULLOPT |
39 | Werner | 168 | |
40 | Werner | 169 | |
158 | werner | 170 | void Tree::applyLIP() |
77 | Werner | 171 | { |
144 | Werner | 172 | if (!mStamp) |
173 | return; |
||
106 | Werner | 174 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
156 | werner | 175 | QPoint pos = mPositionIndex; |
106 | Werner | 176 | int offset = mStamp->offset(); |
77 | Werner | 177 | pos-=QPoint(offset, offset); |
178 | |||
179 | float local_dom; // height of Z* on the current position |
||
180 | int x,y; |
||
401 | werner | 181 | float value, z, z_zstar; |
106 | Werner | 182 | int gr_stamp = mStamp->size(); |
705 | werner | 183 | |
106 | Werner | 184 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
407 | werner | 185 | // this should not happen because of the buffer |
77 | Werner | 186 | return; |
187 | } |
||
705 | werner | 188 | int grid_y = pos.y(); |
77 | Werner | 189 | for (y=0;y<gr_stamp; ++y) { |
403 | werner | 190 | |
705 | werner | 191 | float *grid_value_ptr = mGrid->ptr(pos.x(), grid_y); |
192 | int grid_x = pos.x(); |
||
193 | for (x=0;x<gr_stamp;++x, ++grid_x, ++grid_value_ptr) { |
||
77 | Werner | 194 | // suppose there is no stamping outside |
106 | Werner | 195 | value = (*mStamp)(x,y); // stampvalue |
705 | werner | 196 | //if (value>0.f) { |
197 | local_dom = (*mHeightGrid)(grid_x/cPxPerHeight, grid_y/cPxPerHeight).height; |
||
884 | werner | 198 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
705 | werner | 199 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
1102 | werner | 200 | value = 1.f - value*mOpacity * z_zstar; // calculated value |
705 | werner | 201 | value = std::max(value, 0.02f); // limit value |
77 | Werner | 202 | |
705 | werner | 203 | *grid_value_ptr *= value; |
204 | //} |
||
403 | werner | 205 | |
77 | Werner | 206 | } |
403 | werner | 207 | grid_y++; |
77 | Werner | 208 | } |
209 | |||
210 | m_statPrint++; // count # of stamp applications... |
||
211 | } |
||
212 | |||
155 | werner | 213 | /// helper function for gluing the edges together |
214 | /// index: index at grid |
||
215 | /// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50) |
||
216 | /// buffer: size of buffer around simulation area (in pixels) |
||
295 | werner | 217 | inline int torusIndex(int index, int count, int buffer, int ru_index) |
155 | werner | 218 | { |
295 | werner | 219 | return buffer + ru_index + (index-buffer+count)%count; |
155 | werner | 220 | } |
62 | Werner | 221 | |
155 | werner | 222 | |
223 | /** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area. |
||
224 | */ |
||
158 | werner | 225 | void Tree::applyLIP_torus() |
155 | werner | 226 | { |
227 | if (!mStamp) |
||
228 | return; |
||
229 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
||
295 | werner | 230 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
387 | werner | 231 | QPoint pos = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
232 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
||
295 | werner | 233 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos.x(), mPositionIndex.y() - pos.y()); // offset of the corner of the resource index |
155 | werner | 234 | |
235 | int offset = mStamp->offset(); |
||
236 | pos-=QPoint(offset, offset); |
||
237 | |||
238 | float local_dom; // height of Z* on the current position |
||
239 | int x,y; |
||
240 | float value; |
||
241 | int gr_stamp = mStamp->size(); |
||
242 | int grid_x, grid_y; |
||
243 | float *grid_value; |
||
244 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
||
245 | // todo: in this case we should use another algorithm!!! necessary???? |
||
246 | return; |
||
247 | } |
||
407 | werner | 248 | float z, z_zstar; |
155 | werner | 249 | int xt, yt; // wraparound coordinates |
250 | for (y=0;y<gr_stamp; ++y) { |
||
251 | grid_y = pos.y() + y; |
||
387 | werner | 252 | yt = torusIndex(grid_y, cPxPerRU,bufferOffset, ru_offset.y()); // 50 cells per 100m |
155 | werner | 253 | for (x=0;x<gr_stamp;++x) { |
254 | // suppose there is no stamping outside |
||
255 | grid_x = pos.x() + x; |
||
387 | werner | 256 | xt = torusIndex(grid_x,cPxPerRU,bufferOffset, ru_offset.x()); |
155 | werner | 257 | |
387 | werner | 258 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight,yt/cPxPerHeight).height; |
407 | werner | 259 | |
884 | werner | 260 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
407 | werner | 261 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
155 | werner | 262 | value = (*mStamp)(x,y); // stampvalue |
1102 | werner | 263 | value = 1.f - value*mOpacity * z_zstar; // calculated value |
407 | werner | 264 | // old: value = 1. - value*mOpacity / local_dom; // calculated value |
155 | werner | 265 | value = qMax(value, 0.02f); // limit value |
266 | |||
267 | grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates |
||
268 | *grid_value *= value; |
||
269 | } |
||
270 | } |
||
271 | |||
272 | m_statPrint++; // count # of stamp applications... |
||
273 | } |
||
274 | |||
74 | Werner | 275 | /** heightGrid() |
276 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid. |
||
277 | */ |
||
278 | void Tree::heightGrid() |
||
279 | { |
||
280 | |||
387 | werner | 281 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
74 | Werner | 282 | |
151 | iland | 283 | // count trees that are on height-grid cells (used for stockable area) |
285 | werner | 284 | mHeightGrid->valueAtIndex(p).increaseCount(); |
401 | werner | 285 | if (mHeight > mHeightGrid->valueAtIndex(p).height) |
286 | mHeightGrid->valueAtIndex(p).height=mHeight; |
||
406 | werner | 287 | |
288 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
||
289 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
||
290 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
||
291 | if (index_eastwest - r < 0) { // east |
||
410 | werner | 292 | mHeightGrid->valueAtIndex(p.x()-1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()-1, p.y()).height,mHeight); |
406 | werner | 293 | } |
294 | if (index_eastwest + r >= cPxPerHeight) { // west |
||
410 | werner | 295 | mHeightGrid->valueAtIndex(p.x()+1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()+1, p.y()).height,mHeight); |
406 | werner | 296 | } |
297 | if (index_northsouth - r < 0) { // south |
||
410 | werner | 298 | mHeightGrid->valueAtIndex(p.x(), p.y()-1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()-1).height,mHeight); |
406 | werner | 299 | } |
300 | if (index_northsouth + r >= cPxPerHeight) { // north |
||
410 | werner | 301 | mHeightGrid->valueAtIndex(p.x(), p.y()+1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()+1).height,mHeight); |
406 | werner | 302 | } |
303 | |||
304 | |||
401 | werner | 305 | // without spread of the height grid |
151 | iland | 306 | |
401 | werner | 307 | // // height of Z* |
308 | // const float cellsize = mHeightGrid->cellsize(); |
||
309 | // |
||
310 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
||
311 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
||
312 | // int dist[9]; |
||
313 | // dist[3] = index_northsouth * 2 + 1; // south |
||
314 | // dist[1] = index_eastwest * 2 + 1; // west |
||
315 | // dist[5] = 10 - dist[3]; // north |
||
316 | // dist[7] = 10 - dist[1]; // east |
||
317 | // dist[8] = qMax(dist[5], dist[7]); // north-east |
||
318 | // dist[6] = qMax(dist[3], dist[7]); // south-east |
||
319 | // dist[0] = qMax(dist[3], dist[1]); // south-west |
||
320 | // dist[2] = qMax(dist[5], dist[1]); // north-west |
||
321 | // dist[4] = 0; // center cell |
||
322 | // /* 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: |
||
323 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above. |
||
324 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2 |
||
325 | // */ |
||
326 | // |
||
327 | // |
||
328 | // int ringcount = int(floor(mHeight / cellsize)) + 1; |
||
329 | // int ix, iy; |
||
330 | // int ring; |
||
331 | // float hdom; |
||
332 | // |
||
333 | // for (ix=-ringcount;ix<=ringcount;ix++) |
||
334 | // for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
335 | // ring = qMax(abs(ix), abs(iy)); |
||
336 | // QPoint pos(ix+p.x(), iy+p.y()); |
||
337 | // if (mHeightGrid->isIndexValid(pos)) { |
||
338 | // float &rHGrid = mHeightGrid->valueAtIndex(pos).height; |
||
339 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree |
||
340 | // continue; |
||
341 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y) |
||
342 | // hdom = mHeight - dist[direction]; |
||
343 | // if (ring>1) |
||
344 | // hdom -= (ring-1)*10; |
||
345 | // |
||
346 | // rHGrid = qMax(rHGrid, hdom); // write value |
||
347 | // } // is valid |
||
348 | // } // for (y) |
||
39 | Werner | 349 | } |
40 | Werner | 350 | |
407 | werner | 351 | void Tree::heightGrid_torus() |
352 | { |
||
353 | // height of Z* |
||
155 | werner | 354 | |
407 | werner | 355 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
356 | int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer (i.e.: size of buffer in height-pixels) |
||
357 | p.setX((p.x()-bufferOffset)%10 + bufferOffset); // 10: 10 x 10m pixeln in 100m |
||
358 | p.setY((p.y()-bufferOffset)%10 + bufferOffset); |
||
155 | werner | 359 | |
407 | werner | 360 | |
361 | // torus coordinates: ru_offset = coords of lower left corner of 1ha patch |
||
362 | QPoint ru_offset =QPoint(mPositionIndex.x()/cPxPerHeight - p.x(), mPositionIndex.y()/cPxPerHeight - p.y()); |
||
363 | |||
364 | // count trees that are on height-grid cells (used for stockable area) |
||
365 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
||
366 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
||
367 | v.increaseCount(); |
||
368 | v.height = qMax(v.height, mHeight); |
||
369 | |||
370 | |||
371 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
||
372 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
||
373 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
||
374 | if (index_eastwest - r < 0) { // east |
||
375 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()-1,10,bufferOffset,ru_offset.x()), |
||
376 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
||
410 | werner | 377 | v.height = qMax(v.height, mHeight); |
407 | werner | 378 | } |
379 | if (index_eastwest + r >= cPxPerHeight) { // west |
||
380 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()+1,10,bufferOffset,ru_offset.x()), |
||
381 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
||
410 | werner | 382 | v.height = qMax(v.height, mHeight); |
407 | werner | 383 | } |
384 | if (index_northsouth - r < 0) { // south |
||
385 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
||
386 | torusIndex(p.y()-1,10,bufferOffset,ru_offset.y())); |
||
410 | werner | 387 | v.height = qMax(v.height, mHeight); |
407 | werner | 388 | } |
389 | if (index_northsouth + r >= cPxPerHeight) { // north |
||
390 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
||
391 | torusIndex(p.y()+1,10,bufferOffset,ru_offset.y())); |
||
410 | werner | 392 | v.height = qMax(v.height, mHeight); |
407 | werner | 393 | } |
394 | |||
395 | |||
396 | |||
397 | |||
398 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
||
399 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
||
400 | // int dist[9]; |
||
401 | // dist[3] = index_northsouth * 2 + 1; // south |
||
402 | // dist[1] = index_eastwest * 2 + 1; // west |
||
403 | // dist[5] = 10 - dist[3]; // north |
||
404 | // dist[7] = 10 - dist[1]; // east |
||
405 | // dist[8] = qMax(dist[5], dist[7]); // north-east |
||
406 | // dist[6] = qMax(dist[3], dist[7]); // south-east |
||
407 | // dist[0] = qMax(dist[3], dist[1]); // south-west |
||
408 | // dist[2] = qMax(dist[5], dist[1]); // north-west |
||
409 | // dist[4] = 0; // center cell |
||
410 | // /* 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: |
||
411 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above. |
||
412 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2 |
||
413 | // */ |
||
414 | // |
||
415 | // |
||
416 | // int ringcount = int(floor(mHeight / cellsize)) + 1; |
||
417 | // int ix, iy; |
||
418 | // int ring; |
||
419 | // float hdom; |
||
420 | // for (ix=-ringcount;ix<=ringcount;ix++) |
||
421 | // for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
422 | // ring = qMax(abs(ix), abs(iy)); |
||
423 | // QPoint pos(ix+p.x(), iy+p.y()); |
||
424 | // QPoint p_torus(torusIndex(pos.x(),10,bufferOffset,ru_offset.x()), |
||
425 | // torusIndex(pos.y(),10,bufferOffset,ru_offset.y())); |
||
426 | // if (mHeightGrid->isIndexValid(p_torus)) { |
||
427 | // float &rHGrid = mHeightGrid->valueAtIndex(p_torus.x(),p_torus.y()).height; |
||
428 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree |
||
429 | // continue; |
||
430 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y) |
||
431 | // hdom = mHeight - dist[direction]; |
||
432 | // if (ring>1) |
||
433 | // hdom -= (ring-1)*10; |
||
434 | // |
||
435 | // rHGrid = qMax(rHGrid, hdom); // write value |
||
436 | // } // is valid |
||
437 | // } // for (y) |
||
438 | } |
||
439 | |||
440 | |||
718 | werner | 441 | /** reads the light influence field value for a tree. |
442 | The LIF field is scanned within the crown area of the focal tree, and the influence of |
||
443 | the focal tree is "subtracted" from the LIF values. |
||
444 | Finally, the "LRI correction" is applied. |
||
445 | see http://iland.boku.ac.at/competition+for+light for details. |
||
446 | */ |
||
158 | werner | 447 | void Tree::readLIF() |
40 | Werner | 448 | { |
106 | Werner | 449 | if (!mStamp) |
155 | werner | 450 | return; |
451 | const Stamp *reader = mStamp->reader(); |
||
452 | if (!reader) |
||
453 | return; |
||
156 | werner | 454 | QPoint pos_reader = mPositionIndex; |
780 | werner | 455 | const float outside_area_factor = 0.1f; // |
155 | werner | 456 | |
457 | int offset_reader = reader->offset(); |
||
458 | int offset_writer = mStamp->offset(); |
||
459 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
||
460 | |||
461 | pos_reader-=QPoint(offset_reader, offset_reader); |
||
40 | Werner | 462 | |
155 | werner | 463 | float local_dom; |
464 | |||
40 | Werner | 465 | int x,y; |
466 | double sum=0.; |
||
155 | werner | 467 | double value, own_value; |
468 | float *grid_value; |
||
403 | werner | 469 | float z, z_zstar; |
155 | werner | 470 | int reader_size = reader->size(); |
471 | int rx = pos_reader.x(); |
||
472 | int ry = pos_reader.y(); |
||
473 | for (y=0;y<reader_size; ++y, ++ry) { |
||
474 | grid_value = mGrid->ptr(rx, ry); |
||
475 | for (x=0;x<reader_size;++x) { |
||
476 | |||
718 | werner | 477 | const HeightGridValue &hgv = mHeightGrid->constValueAtIndex((rx+x)/cPxPerHeight, ry/cPxPerHeight); // the height grid value, ry: gets ++ed in outer loop, rx not |
478 | local_dom = hgv.height; |
||
884 | werner | 479 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
403 | werner | 480 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
155 | werner | 481 | |
403 | werner | 482 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
155 | werner | 483 | own_value = qMax(own_value, 0.02); |
484 | value = *grid_value++ / own_value; // remove impact of focal tree |
||
720 | werner | 485 | // additional punishment if pixel is outside: |
718 | werner | 486 | if (hgv.isForestOutside()) |
720 | werner | 487 | value *= outside_area_factor; |
403 | werner | 488 | |
489 | //qDebug() << x << y << local_dom << z << z_zstar << own_value << value << *(grid_value-1) << (*reader)(x,y) << mStamp->offsetValue(x,y,d_offset); |
||
155 | werner | 490 | //if (value>0.) |
491 | sum += value * (*reader)(x,y); |
||
492 | |||
40 | Werner | 493 | } |
494 | } |
||
1102 | werner | 495 | mLRI = static_cast<float>( sum ); |
426 | werner | 496 | // LRI correction... |
497 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
||
498 | if (hrel<1.) |
||
1102 | werner | 499 | mLRI = static_cast<float>( species()->speciesSet()->LRIcorrection(mLRI, hrel) ); |
426 | werner | 500 | |
403 | werner | 501 | |
155 | werner | 502 | if (mLRI > 1.) |
503 | mLRI = 1.; |
||
206 | werner | 504 | |
505 | // Finally, add LRI of this Tree to the ResourceUnit! |
||
251 | werner | 506 | mRU->addWLA(mLeafArea, mLRI); |
206 | werner | 507 | |
155 | werner | 508 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact; |
40 | Werner | 509 | } |
510 | |||
155 | werner | 511 | /// Torus version of read stamp (glued edges) |
158 | werner | 512 | void Tree::readLIF_torus() |
78 | Werner | 513 | { |
106 | Werner | 514 | if (!mStamp) |
107 | Werner | 515 | return; |
106 | Werner | 516 | const Stamp *reader = mStamp->reader(); |
78 | Werner | 517 | if (!reader) |
107 | Werner | 518 | return; |
295 | werner | 519 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
78 | Werner | 520 | |
387 | werner | 521 | QPoint pos_reader = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
522 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
||
295 | werner | 523 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos_reader.x(), mPositionIndex.y() - pos_reader.y()); // offset of the corner of the resource index |
524 | |||
78 | Werner | 525 | int offset_reader = reader->offset(); |
106 | Werner | 526 | int offset_writer = mStamp->offset(); |
78 | Werner | 527 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
528 | |||
529 | pos_reader-=QPoint(offset_reader, offset_reader); |
||
530 | |||
531 | float local_dom; |
||
532 | |||
533 | int x,y; |
||
534 | double sum=0.; |
||
535 | double value, own_value; |
||
536 | float *grid_value; |
||
407 | werner | 537 | float z, z_zstar; |
78 | Werner | 538 | int reader_size = reader->size(); |
539 | int rx = pos_reader.x(); |
||
540 | int ry = pos_reader.y(); |
||
155 | werner | 541 | int xt, yt; // wrapped coords |
542 | |||
397 | werner | 543 | for (y=0;y<reader_size; ++y) { |
544 | yt = torusIndex(ry+y,cPxPerRU, bufferOffset, ru_offset.y()); |
||
78 | Werner | 545 | for (x=0;x<reader_size;++x) { |
387 | werner | 546 | xt = torusIndex(rx+x,cPxPerRU, bufferOffset, ru_offset.x()); |
155 | werner | 547 | grid_value = mGrid->ptr(xt,yt); |
407 | werner | 548 | |
387 | werner | 549 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight, yt/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not |
884 | werner | 550 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
407 | werner | 551 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
125 | Werner | 552 | |
407 | werner | 553 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
554 | // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height; |
||
78 | Werner | 555 | own_value = qMax(own_value, 0.02); |
407 | werner | 556 | value = *grid_value++ / own_value; // remove impact of focal tree |
557 | |||
397 | werner | 558 | // debug for one tree in HJA |
559 | //if (id()==178020) |
||
560 | // qDebug() << x << y << xt << yt << *grid_value << local_dom << own_value << value << (*reader)(x,y); |
||
321 | werner | 561 | //if (_isnan(value)) |
562 | // qDebug() << "isnan" << id(); |
||
397 | werner | 563 | if (value * (*reader)(x,y)>1.) |
564 | qDebug() << "LIFTorus: value>1."; |
||
78 | Werner | 565 | sum += value * (*reader)(x,y); |
566 | |||
567 | //} // isIndexValid |
||
568 | } |
||
569 | } |
||
1102 | werner | 570 | mLRI = static_cast<float>( sum ); |
426 | werner | 571 | |
572 | // LRI correction... |
||
573 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
||
574 | if (hrel<1.) |
||
1102 | werner | 575 | mLRI = static_cast<float>( species()->speciesSet()->LRIcorrection(mLRI, hrel) ); |
426 | werner | 576 | |
615 | werner | 577 | if (isnan(mLRI)) { |
321 | werner | 578 | qDebug() << "LRI invalid (nan)!" << id(); |
579 | mLRI=0.; |
||
580 | //qDebug() << reader->dump(); |
||
581 | } |
||
148 | iland | 582 | if (mLRI > 1.) |
583 | mLRI = 1.; |
||
78 | Werner | 584 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact; |
205 | werner | 585 | |
586 | // Finally, add LRI of this Tree to the ResourceUnit! |
||
251 | werner | 587 | mRU->addWLA(mLeafArea, mLRI); |
58 | Werner | 588 | } |
589 | |||
155 | werner | 590 | |
40 | Werner | 591 | void Tree::resetStatistics() |
592 | { |
||
593 | m_statPrint=0; |
||
105 | Werner | 594 | m_statCreated=0; |
48 | Werner | 595 | m_statAboveZ=0; |
40 | Werner | 596 | m_nextId=1; |
597 | } |
||
107 | Werner | 598 | |
1072 | werner | 599 | #ifdef ALT_TREE_MORTALITY |
1071 | werner | 600 | void Tree::mortalityParams(double dbh_inc_threshold, int stress_years, double stress_mort_prob) |
601 | { |
||
602 | _stress_threshold = dbh_inc_threshold; |
||
603 | _stress_years = stress_years; |
||
604 | _stress_death_prob = stress_mort_prob; |
||
605 | qDebug() << "Alternative Mortality enabled: threshold" << dbh_inc_threshold << ", years:" << _stress_years << ", level:" << _stress_death_prob; |
||
606 | } |
||
1072 | werner | 607 | #endif |
1071 | werner | 608 | |
251 | werner | 609 | void Tree::calcLightResponse() |
610 | { |
||
611 | // calculate a light response from lri: |
||
298 | werner | 612 | // http://iland.boku.ac.at/individual+tree+light+availability |
470 | werner | 613 | double lri = limit(mLRI * mRU->LRImodifier(), 0., 1.); // Eq. (3) |
1102 | werner | 614 | mLightResponse = static_cast<float>( mSpecies->lightResponse(lri) ); // Eq. (4) |
251 | werner | 615 | mRU->addLR(mLeafArea, mLightResponse); |
616 | |||
617 | } |
||
618 | |||
110 | Werner | 619 | ////////////////////////////////////////////////// |
620 | //// Growth Functions |
||
621 | ////////////////////////////////////////////////// |
||
107 | Werner | 622 | |
227 | werner | 623 | /** grow() is the main function of the yearly tree growth. |
624 | The main steps are: |
||
298 | werner | 625 | - Production of GPP/NPP @sa http://iland.boku.ac.at/primary+production http://iland.boku.ac.at/individual+tree+light+availability |
626 | - Partitioning of NPP to biomass compartments of the tree @sa http://iland.boku.ac.at/allocation |
||
227 | werner | 627 | - Growth of the stem http://iland.boku.ac.at/stem+growth (???) |
387 | werner | 628 | Further activties: * the age of the tree is increased |
629 | * the mortality sub routine is executed |
||
630 | * seeds are produced */ |
||
107 | Werner | 631 | void Tree::grow() |
632 | { |
||
159 | werner | 633 | TreeGrowthData d; |
169 | werner | 634 | mAge++; // increase age |
230 | werner | 635 | // step 1: get "interception area" of the tree individual [m2] |
636 | // the sum of all area of all trees of a unit equal the total stocked area * interception_factor(Beer-Lambert) |
||
637 | double effective_area = mRU->interceptedArea(mLeafArea, mLightResponse); |
||
107 | Werner | 638 | |
230 | werner | 639 | // step 2: calculate GPP of the tree based |
640 | // (1) get the amount of GPP for a "unit area" of the tree species |
||
641 | double raw_gpp_per_area = mRU->resourceUnitSpecies(species()).prod3PG().GPPperArea(); |
||
642 | // (2) GPP (without aging-effect) in kg Biomass / year |
||
643 | double raw_gpp = raw_gpp_per_area * effective_area; |
||
161 | werner | 644 | |
227 | werner | 645 | // apply aging according to the state of the individuum |
388 | werner | 646 | const double aging_factor = mSpecies->aging(mHeight, mAge); |
376 | werner | 647 | mRU->addTreeAging(mLeafArea, aging_factor); |
227 | werner | 648 | double gpp = raw_gpp * aging_factor; // |
608 | werner | 649 | d.NPP = gpp * cAutotrophicRespiration; // respiration loss (0.47), cf. Waring et al 1998. |
113 | Werner | 650 | |
279 | werner | 651 | //DBGMODE( |
137 | Werner | 652 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) { |
133 | Werner | 653 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP); |
654 | dumpList(out); // add tree headers |
||
299 | werner | 655 | out << mLRI * mRU->LRImodifier() << mLightResponse << effective_area << raw_gpp << gpp << d.NPP << aging_factor; |
133 | Werner | 656 | } |
279 | werner | 657 | //); // DBGMODE() |
1196 | werner | 658 | if (Globals->model()->settings().growthEnabled) { |
659 | if (d.NPP>0.) |
||
660 | partitioning(d); // split npp to compartments and grow (diameter, height) |
||
661 | } |
||
133 | Werner | 662 | |
387 | werner | 663 | // mortality |
1071 | werner | 664 | #ifdef ALT_TREE_MORTALITY |
665 | // alternative variant of tree mortality (note: mStrssIndex used otherwise) |
||
666 | altMortality(d); |
||
667 | |||
668 | #else |
||
200 | werner | 669 | if (Model::settings().mortalityEnabled) |
670 | mortality(d); |
||
110 | Werner | 671 | |
159 | werner | 672 | mStressIndex = d.stress_index; |
1071 | werner | 673 | #endif |
180 | werner | 674 | |
675 | if (!isDead()) |
||
257 | werner | 676 | mRU->resourceUnitSpecies(species()).statistics().add(this, &d); |
277 | werner | 677 | |
387 | werner | 678 | // regeneration |
1180 | werner | 679 | mSpecies->seedProduction(this); |
387 | werner | 680 | |
107 | Werner | 681 | } |
682 | |||
227 | werner | 683 | /** partitioning of this years assimilates (NPP) to biomass compartments. |
298 | werner | 684 | Conceptionally, the algorithm is based on Duursma, 2007. |
685 | @sa http://iland.boku.ac.at/allocation */ |
||
159 | werner | 686 | inline void Tree::partitioning(TreeGrowthData &d) |
115 | Werner | 687 | { |
159 | werner | 688 | double npp = d.NPP; |
115 | Werner | 689 | // add content of reserve pool |
116 | Werner | 690 | npp += mNPPReserve; |
159 | werner | 691 | const double foliage_mass_allo = species()->biomassFoliage(mDbh); |
276 | werner | 692 | const double reserve_size = foliage_mass_allo * (1. + mSpecies->finerootFoliageRatio()); |
297 | werner | 693 | double refill_reserve = qMin(reserve_size, (1. + mSpecies->finerootFoliageRatio())*mFoliageMass); // not always try to refill reserve 100% |
119 | Werner | 694 | |
136 | Werner | 695 | double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta) |
468 | werner | 696 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
117 | Werner | 697 | // turnover rates |
159 | werner | 698 | const double &to_fol = species()->turnoverLeaf(); |
699 | const double &to_root = species()->turnoverRoot(); |
||
136 | Werner | 700 | // the turnover rate of wood depends on the size of the reserve pool: |
116 | Werner | 701 | |
136 | Werner | 702 | |
163 | werner | 703 | double to_wood = refill_reserve / (mWoodyMass + refill_reserve); |
704 | |||
468 | werner | 705 | apct_root = rus.prod3PG().rootFraction(); |
261 | werner | 706 | d.NPP_above = d.NPP * (1. - apct_root); // aboveground: total NPP - fraction to roots |
298 | werner | 707 | double b_wf = species()->allometricRatio_wf(); // ratio of allometric exponents (b_woody / b_foliage) |
117 | Werner | 708 | |
709 | // Duursma 2007, Eq. (20) |
||
167 | werner | 710 | 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 ); |
902 | werner | 711 | |
712 | apct_wood = limit(apct_wood, 0., 1.-apct_root); |
||
713 | |||
117 | Werner | 714 | apct_foliage = 1. - apct_root - apct_wood; |
715 | |||
163 | werner | 716 | |
902 | werner | 717 | DBGMODE( |
163 | werner | 718 | if (apct_foliage<0 || apct_wood<0) |
719 | qDebug() << "transfer to foliage or wood < 0"; |
||
720 | if (npp<0) |
||
721 | qDebug() << "NPP < 0"; |
||
902 | werner | 722 | ); |
163 | werner | 723 | |
136 | Werner | 724 | // Change of biomass compartments |
276 | werner | 725 | double sen_root = mFineRootMass * to_root; |
726 | double sen_foliage = mFoliageMass * to_fol; |
||
521 | werner | 727 | if (ru()->snag()) |
588 | werner | 728 | ru()->snag()->addTurnoverLitter(this->species(), sen_foliage, sen_root); |
298 | werner | 729 | |
136 | Werner | 730 | // Roots |
298 | werner | 731 | // http://iland.boku.ac.at/allocation#belowground_NPP |
276 | werner | 732 | mFineRootMass -= sen_root; // reduce only fine root pool |
733 | double delta_root = apct_root * npp; |
||
734 | // 1st, refill the fine root pool |
||
735 | double fineroot_miss = mFoliageMass * mSpecies->finerootFoliageRatio() - mFineRootMass; |
||
736 | if (fineroot_miss>0.){ |
||
737 | double delta_fineroot = qMin(fineroot_miss, delta_root); |
||
738 | mFineRootMass += delta_fineroot; |
||
739 | delta_root -= delta_fineroot; |
||
740 | } |
||
741 | // 2nd, the rest of NPP allocated to roots go to coarse roots |
||
595 | werner | 742 | double max_coarse_root = species()->biomassRoot(mDbh); |
276 | werner | 743 | mCoarseRootMass += delta_root; |
595 | werner | 744 | if (mCoarseRootMass > max_coarse_root) { |
745 | // if the coarse root pool exceeds the value given by the allometry, then the |
||
746 | // surplus is accounted as turnover |
||
747 | if (ru()->snag()) |
||
748 | ru()->snag()->addTurnoverWood(species(), mCoarseRootMass-max_coarse_root); |
||
119 | Werner | 749 | |
1102 | werner | 750 | mCoarseRootMass = static_cast<float>( max_coarse_root ); |
595 | werner | 751 | } |
752 | |||
136 | Werner | 753 | // Foliage |
159 | werner | 754 | double delta_foliage = apct_foliage * npp - sen_foliage; |
137 | Werner | 755 | mFoliageMass += delta_foliage; |
615 | werner | 756 | if (isnan(mFoliageMass)) |
217 | werner | 757 | qDebug() << "foliage mass invalid!"; |
163 | werner | 758 | if (mFoliageMass<0.) mFoliageMass=0.; // limit to zero |
759 | |||
1102 | werner | 760 | mLeafArea = static_cast<float>( mFoliageMass * species()->specificLeafArea() ); // update leaf area |
119 | Werner | 761 | |
271 | werner | 762 | // stress index: different varaints at denominator: to_fol*foliage_mass = leafmass to rebuild, |
198 | werner | 763 | // foliage_mass_allo: simply higher chance for stress |
271 | werner | 764 | // note: npp = NPP + reserve (see above) |
276 | werner | 765 | d.stress_index =qMax(1. - (npp) / ( to_fol*foliage_mass_allo + to_root*foliage_mass_allo*species()->finerootFoliageRatio() + reserve_size), 0.); |
198 | werner | 766 | |
136 | Werner | 767 | // Woody compartments |
298 | werner | 768 | // see also: http://iland.boku.ac.at/allocation#reserve_and_allocation_to_stem_growth |
136 | Werner | 769 | // (1) transfer to reserve pool |
770 | double gross_woody = apct_wood * npp; |
||
771 | double to_reserve = qMin(reserve_size, gross_woody); |
||
1102 | werner | 772 | mNPPReserve = static_cast<float>( to_reserve ); |
136 | Werner | 773 | double net_woody = gross_woody - to_reserve; |
137 | Werner | 774 | double net_stem = 0.; |
164 | werner | 775 | mDbhDelta = 0.; |
165 | werner | 776 | |
777 | |||
136 | Werner | 778 | if (net_woody > 0.) { |
779 | // (2) calculate part of increment that is dedicated to the stem (which is a function of diameter) |
||
159 | werner | 780 | net_stem = net_woody * species()->allometricFractionStem(mDbh); |
781 | d.NPP_stem = net_stem; |
||
137 | Werner | 782 | mWoodyMass += net_woody; |
136 | Werner | 783 | // (3) growth of diameter and height baseed on net stem increment |
159 | werner | 784 | grow_diameter(d); |
136 | Werner | 785 | } |
119 | Werner | 786 | |
279 | werner | 787 | //DBGMODE( |
137 | Werner | 788 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreePartition) |
789 | && isDebugging() ) { |
||
129 | Werner | 790 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreePartition); |
791 | dumpList(out); // add tree headers |
||
136 | Werner | 792 | out << npp << apct_foliage << apct_wood << apct_root |
276 | werner | 793 | << delta_foliage << net_woody << delta_root << mNPPReserve << net_stem << d.stress_index; |
137 | Werner | 794 | } |
144 | Werner | 795 | |
279 | werner | 796 | //); // DBGMODE() |
497 | werner | 797 | DBGMODE( |
428 | werner | 798 | if (mWoodyMass<0. || mWoodyMass>50000 || mFoliageMass<0. || mFoliageMass>2000. || mCoarseRootMass<0. || mCoarseRootMass>30000 |
393 | werner | 799 | || mNPPReserve>4000.) { |
389 | werner | 800 | qDebug() << "Tree:partitioning: invalid or unlikely pools."; |
144 | Werner | 801 | qDebug() << GlobalSettings::instance()->debugListCaptions(GlobalSettings::DebugOutputs(0)); |
802 | DebugList dbg; dumpList(dbg); |
||
803 | qDebug() << dbg; |
||
497 | werner | 804 | } ); |
144 | Werner | 805 | |
136 | Werner | 806 | /*DBG_IF_X(mId == 1 , "Tree::partitioning", "dump", dump() |
807 | + 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") |
||
808 | .arg(npp).arg(senescence_foliage).arg(senescence_stem).arg(senescence_root) |
||
809 | .arg(net_foliage).arg(net_stem).arg(net_root).arg(to_reserve).arg(mNPPReserve) );*/ |
||
129 | Werner | 810 | |
115 | Werner | 811 | } |
812 | |||
125 | Werner | 813 | |
134 | Werner | 814 | /** Determination of diamter and height growth based on increment of the stem mass (@p net_stem_npp). |
125 | Werner | 815 | Refer to XXX for equations and variables. |
816 | This function updates the dbh and height of the tree. |
||
227 | werner | 817 | The equations are based on dbh in meters! */ |
159 | werner | 818 | inline void Tree::grow_diameter(TreeGrowthData &d) |
119 | Werner | 819 | { |
820 | // determine dh-ratio of increment |
||
821 | // height increment is a function of light competition: |
||
125 | Werner | 822 | double hd_growth = relative_height_growth(); // hd of height growth |
153 | werner | 823 | double d_m = mDbh / 100.; // current diameter in [m] |
159 | werner | 824 | double net_stem_npp = d.NPP_stem; |
825 | |||
153 | werner | 826 | const double d_delta_m = mDbhDelta / 100.; // increment of last year in [m] |
115 | Werner | 827 | |
159 | werner | 828 | const double mass_factor = species()->volumeFactor() * species()->density(); |
153 | werner | 829 | double stem_mass = mass_factor * d_m*d_m * mHeight; // result: kg, dbh[cm], h[meter] |
123 | Werner | 830 | |
153 | werner | 831 | // factor is in diameter increment per NPP [m/kg] |
832 | double factor_diameter = 1. / ( mass_factor * (d_m + d_delta_m)*(d_m + d_delta_m) * ( 2. * mHeight/d_m + hd_growth) ); |
||
125 | Werner | 833 | double delta_d_estimate = factor_diameter * net_stem_npp; // estimated dbh-inc using last years increment |
834 | |||
835 | // using that dbh-increment we estimate a stem-mass-increment and the residual (Eq. 9) |
||
153 | werner | 836 | double stem_estimate = mass_factor * (d_m + delta_d_estimate)*(d_m + delta_d_estimate)*(mHeight + delta_d_estimate*hd_growth); |
137 | Werner | 837 | double stem_residual = stem_estimate - (stem_mass + net_stem_npp); |
125 | Werner | 838 | |
839 | // the final increment is then: |
||
840 | double d_increment = factor_diameter * (net_stem_npp - stem_residual); // Eq. (11) |
||
463 | werner | 841 | double res_final = 0.; |
1118 | werner | 842 | if (fabs(stem_residual) > std::min(1.,stem_mass)) { |
465 | werner | 843 | |
463 | werner | 844 | // calculate final residual in stem |
845 | res_final = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)); |
||
1118 | werner | 846 | if (fabs(res_final)>std::min(1.,stem_mass)) { |
847 | // for large errors in stem biomass due to errors in diameter increment (> 1kg or >stem mass), we solve the increment iteratively. |
||
463 | werner | 848 | // first, increase increment with constant step until we overestimate the first time |
849 | // then, |
||
850 | d_increment = 0.02; // start with 2cm increment |
||
851 | bool reached_error = false; |
||
852 | double step=0.01; // step-width 1cm |
||
853 | double est_stem; |
||
854 | do { |
||
855 | est_stem = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth); // estimate with current increment |
||
856 | stem_residual = est_stem - (stem_mass + net_stem_npp); |
||
857 | |||
858 | if (fabs(stem_residual) <1.) // finished, if stem residual below 1kg |
||
859 | break; |
||
860 | if (stem_residual > 0.) { |
||
861 | d_increment -= step; |
||
862 | reached_error=true; |
||
863 | } else { |
||
864 | d_increment += step; |
||
865 | } |
||
866 | if (reached_error) |
||
867 | step /= 2.; |
||
868 | } while (step>0.00001); // continue until diameter "accuracy" falls below 1/100mm |
||
869 | } |
||
870 | } |
||
871 | |||
872 | if (d_increment<0.f) |
||
873 | qDebug() << "Tree::grow_diameter: d_inc < 0."; |
||
144 | Werner | 874 | DBG_IF_X(d_increment<0. || d_increment>0.1, "Tree::grow_dimater", "increment out of range.", dump() |
125 | Werner | 875 | + QString("\nhdz %1 factor_diameter %2 stem_residual %3 delta_d_estimate %4 d_increment %5 final residual(kg) %6") |
876 | .arg(hd_growth).arg(factor_diameter).arg(stem_residual).arg(delta_d_estimate).arg(d_increment) |
||
142 | Werner | 877 | .arg( mass_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)) )); |
125 | Werner | 878 | |
303 | werner | 879 | //DBGMODE( |
463 | werner | 880 | // do not calculate res_final twice if already done |
881 | DBG_IF_X( (res_final==0.?fabs(mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp))):res_final) > 1, "Tree::grow_diameter", "final residual stem estimate > 1kg", dump()); |
||
153 | werner | 882 | 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 | 883 | |
137 | Werner | 884 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeGrowth) && isDebugging() ) { |
126 | Werner | 885 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeGrowth); |
129 | Werner | 886 | dumpList(out); // add tree headers |
143 | Werner | 887 | out << net_stem_npp << stem_mass << hd_growth << factor_diameter << delta_d_estimate*100 << d_increment*100; |
126 | Werner | 888 | } |
153 | werner | 889 | |
303 | werner | 890 | //); // DBGMODE() |
125 | Werner | 891 | |
892 | d_increment = qMax(d_increment, 0.); |
||
893 | |||
894 | // update state variables |
||
1102 | werner | 895 | mDbh += d_increment*100.f; // convert from [m] to [cm] |
896 | mDbhDelta = static_cast<float>( d_increment*100. ); // save for next year's growth |
||
153 | werner | 897 | mHeight += d_increment * hd_growth; |
158 | werner | 898 | |
899 | // update state of LIP stamp and opacity |
||
159 | werner | 900 | mStamp = species()->stamp(mDbh, mHeight); // get new stamp for updated dimensions |
158 | werner | 901 | // calculate the CrownFactor which reflects the opacity of the crown |
200 | werner | 902 | const double k=Model::settings().lightExtinctionCoefficientOpacity; |
1102 | werner | 903 | mOpacity = static_cast<float>( 1. - exp(-k * mLeafArea / mStamp->crownArea()) ); |
158 | werner | 904 | |
119 | Werner | 905 | } |
906 | |||
125 | Werner | 907 | |
908 | /// return the HD ratio of this year's increment based on the light status. |
||
119 | Werner | 909 | inline double Tree::relative_height_growth() |
910 | { |
||
911 | double hd_low, hd_high; |
||
912 | mSpecies->hdRange(mDbh, hd_low, hd_high); |
||
913 | |||
125 | Werner | 914 | DBG_IF_X(hd_low>hd_high, "Tree::relative_height_growth", "hd low higher dann hd_high for ", dump()); |
949 | werner | 915 | DBG_IF_X(hd_low < 10 || 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)); |
125 | Werner | 916 | |
917 | // scale according to LRI: if receiving much light (LRI=1), the result is hd_low (for open grown trees) |
||
326 | werner | 918 | // use the corrected LRI (see tracker#11) |
919 | double lri = limit(mLRI * mRU->LRImodifier(),0.,1.); |
||
920 | double hd_ratio = hd_high - (hd_high-hd_low)*lri; |
||
125 | Werner | 921 | return hd_ratio; |
119 | Werner | 922 | } |
141 | Werner | 923 | |
278 | werner | 924 | /** This function is called if a tree dies. |
925 | @sa ResourceUnit::cleanTreeList(), remove() */ |
||
277 | werner | 926 | void Tree::die(TreeGrowthData *d) |
927 | { |
||
928 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
||
664 | werner | 929 | mRU->treeDied(); |
468 | werner | 930 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
931 | rus.statisticsDead().add(this, d); // add tree to statistics |
||
1064 | werner | 932 | notifyTreeRemoved(TreeDeath); |
521 | werner | 933 | if (ru()->snag()) |
934 | ru()->snag()->addMortality(this); |
||
277 | werner | 935 | } |
936 | |||
903 | werner | 937 | /// remove a tree (most likely due to harvest) from the system. |
564 | werner | 938 | void Tree::remove(double removeFoliage, double removeBranch, double removeStem ) |
278 | werner | 939 | { |
940 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
||
1050 | werner | 941 | setIsHarvested(); |
664 | werner | 942 | mRU->treeDied(); |
468 | werner | 943 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
944 | rus.statisticsMgmt().add(this, 0); |
||
1157 | werner | 945 | if (isCutdown()) |
946 | notifyTreeRemoved(TreeCutDown); |
||
947 | else |
||
948 | notifyTreeRemoved(TreeHarvest); |
||
903 | werner | 949 | |
1165 | werner | 950 | if (GlobalSettings::instance()->model()->saplings()) |
951 | GlobalSettings::instance()->model()->saplings()->addSprout(this); |
||
952 | |||
521 | werner | 953 | if (ru()->snag()) |
564 | werner | 954 | ru()->snag()->addHarvest(this, removeStem, removeBranch, removeFoliage); |
278 | werner | 955 | } |
956 | |||
713 | werner | 957 | /// remove the tree due to an special event (disturbance) |
903 | werner | 958 | /// this is +- the same as die(). |
713 | werner | 959 | void Tree::removeDisturbance(const double stem_to_soil_fraction, const double stem_to_snag_fraction, const double branch_to_soil_fraction, const double branch_to_snag_fraction, const double foliage_to_soil_fraction) |
960 | { |
||
961 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
||
962 | mRU->treeDied(); |
||
963 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
||
714 | werner | 964 | rus.statisticsDead().add(this, 0); |
1064 | werner | 965 | notifyTreeRemoved(TreeDisturbance); |
903 | werner | 966 | |
1165 | werner | 967 | if (GlobalSettings::instance()->model()->saplings()) |
968 | GlobalSettings::instance()->model()->saplings()->addSprout(this); |
||
1088 | werner | 969 | |
1100 | werner | 970 | if (ru()->snag()) { |
1102 | werner | 971 | if (isHarvested()) { // if the tree is harvested, do the same as in normal tree harvest (but with default values) |
1088 | werner | 972 | ru()->snag()->addHarvest(this, 1., 0., 0.); |
1102 | werner | 973 | } else { |
1088 | werner | 974 | ru()->snag()->addDisturbance(this, stem_to_snag_fraction, stem_to_soil_fraction, branch_to_snag_fraction, branch_to_soil_fraction, foliage_to_soil_fraction); |
1102 | werner | 975 | } |
1100 | werner | 976 | } |
713 | werner | 977 | } |
978 | |||
903 | werner | 979 | /// remove a part of the biomass of the tree, e.g. due to fire. |
980 | void Tree::removeBiomassOfTree(const double removeFoliageFraction, const double removeBranchFraction, const double removeStemFraction) |
||
668 | werner | 981 | { |
982 | mFoliageMass *= 1. - removeFoliageFraction; |
||
983 | mWoodyMass *= (1. - removeStemFraction); |
||
984 | // we have a problem with the branches: this currently cannot be done properly! |
||
766 | werner | 985 | (void) removeBranchFraction; // silence warning |
668 | werner | 986 | } |
987 | |||
975 | werner | 988 | void Tree::setHeight(const float height) |
989 | { |
||
990 | if (height<=0. || height>150.) |
||
991 | qWarning() << "trying to set tree height to invalid value:" << height << " for tree on RU:" << (mRU?mRU->boundingBox():QRect()); |
||
992 | mHeight=height; |
||
993 | } |
||
994 | |||
159 | werner | 995 | void Tree::mortality(TreeGrowthData &d) |
996 | { |
||
163 | werner | 997 | // death if leaf area is 0 |
998 | if (mFoliageMass<0.00001) |
||
999 | die(); |
||
1000 | |||
308 | werner | 1001 | double p_death, p_stress, p_intrinsic; |
1002 | p_intrinsic = species()->deathProb_intrinsic(); |
||
1003 | p_stress = species()->deathProb_stress(d.stress_index); |
||
1004 | p_death = p_intrinsic + p_stress; |
||
707 | werner | 1005 | double p = drandom(); //0..1 |
159 | werner | 1006 | if (p<p_death) { |
1007 | // die... |
||
1008 | die(); |
||
1009 | } |
||
1010 | } |
||
141 | Werner | 1011 | |
1072 | werner | 1012 | #ifdef ALT_TREE_MORTALITY |
1071 | werner | 1013 | void Tree::altMortality(TreeGrowthData &d) |
1014 | { |
||
1015 | // death if leaf area is 0 |
||
1016 | if (mFoliageMass<0.00001) |
||
1017 | die(); |
||
1018 | |||
1019 | double p_intrinsic, p_stress=0.; |
||
1020 | p_intrinsic = species()->deathProb_intrinsic(); |
||
1021 | |||
1022 | if (mDbhDelta < _stress_threshold) { |
||
1023 | mStressIndex++; |
||
1024 | if (mStressIndex> _stress_years) |
||
1025 | p_stress = _stress_death_prob; |
||
1026 | } else |
||
1027 | mStressIndex = 0; |
||
1028 | |||
1029 | double p = drandom(); //0..1 |
||
1030 | if (p<p_intrinsic + p_stress) { |
||
1031 | // die... |
||
1032 | die(); |
||
1033 | } |
||
1034 | } |
||
1072 | werner | 1035 | #endif |
1071 | werner | 1036 | |
1064 | werner | 1037 | void Tree::notifyTreeRemoved(TreeRemovalType reason) |
903 | werner | 1038 | { |
1039 | // add the volume of the current tree to the height grid |
||
1077 | werner | 1040 | // this information is used to track the removed volume for stands based on grids (and for salvaging operations) |
909 | werner | 1041 | ABE::ForestManagementEngine *abe = GlobalSettings::instance()->model()->ABEngine(); |
1042 | if (abe) |
||
1102 | werner | 1043 | abe->notifyTreeRemoval(this, static_cast<int>(reason)); |
1044 | werner | 1044 | |
1045 | // tell disturbance modules that a tree died |
||
1102 | werner | 1046 | GlobalSettings::instance()->model()->modules()->treeDeath(this, static_cast<int>(reason) ); |
1047 | |||
1157 | werner | 1048 | // update reason, if ABE handled the tree |
1049 | if (reason==TreeDisturbance && isHarvested()) |
||
1050 | reason = TreeSalavaged; |
||
1051 | if (isCutdown()) |
||
1052 | reason = TreeCutDown; |
||
1102 | werner | 1053 | // create output for tree removals |
1054 | if (mRemovalOutput && mRemovalOutput->isEnabled()) |
||
1055 | mRemovalOutput->execRemovedTree(this, static_cast<int>(reason)); |
||
1114 | werner | 1056 | if (mLSRemovalOutput && mLSRemovalOutput->isEnabled()) |
1057 | mLSRemovalOutput->execRemovedTree(this, static_cast<int>(reason)); |
||
903 | werner | 1058 | } |
1059 | |||
141 | Werner | 1060 | ////////////////////////////////////////////////// |
1061 | //// value functions |
||
1062 | ////////////////////////////////////////////////// |
||
1063 | |||
145 | Werner | 1064 | double Tree::volume() const |
141 | Werner | 1065 | { |
1066 | /// @see Species::volumeFactor() for details |
||
159 | werner | 1067 | const double volume_factor = species()->volumeFactor(); |
157 | werner | 1068 | const double volume = volume_factor * mDbh*mDbh*mHeight * 0.0001; // dbh in cm: cm/100 * cm/100 = cm*cm * 0.0001 = m2 |
141 | Werner | 1069 | return volume; |
1070 | } |
||
180 | werner | 1071 | |
579 | werner | 1072 | /// return the basal area in m2 |
180 | werner | 1073 | double Tree::basalArea() const |
1074 | { |
||
1075 | // A = r^2 * pi = d/2*pi; from cm->m: d/200 |
||
1076 | const double b = (mDbh/200.)*(mDbh/200.)*M_PI; |
||
1077 | return b; |
||
1078 | } |
||
668 | werner | 1079 |