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