Rev 126 | Rev 130 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
117 | Werner | 2 | #include "global.h" |
3 | #include "tree.h" |
||
3 | Werner | 4 | |
83 | Werner | 5 | #include "grid.h" |
3 | Werner | 6 | |
83 | Werner | 7 | #include "stamp.h" |
90 | Werner | 8 | #include "species.h" |
107 | Werner | 9 | #include "ressourceunit.h" |
38 | Werner | 10 | |
110 | Werner | 11 | // static varaibles |
106 | Werner | 12 | FloatGrid *Tree::mGrid = 0; |
13 | FloatGrid *Tree::mHeightGrid = 0; |
||
40 | Werner | 14 | int Tree::m_statPrint=0; |
48 | Werner | 15 | int Tree::m_statAboveZ=0; |
105 | Werner | 16 | int Tree::m_statCreated=0; |
40 | Werner | 17 | int Tree::m_nextId=0; |
53 | Werner | 18 | float Tree::lafactor = 1.; |
106 | Werner | 19 | int Tree::mDebugid = -1; |
40 | Werner | 20 | |
110 | Werner | 21 | // lifecycle |
3 | Werner | 22 | Tree::Tree() |
23 | { |
||
106 | Werner | 24 | mDbh = 0; |
25 | mHeight = 0; |
||
26 | mSpecies = 0; |
||
107 | Werner | 27 | mRU = 0; |
106 | Werner | 28 | mId = m_nextId++; |
105 | Werner | 29 | m_statCreated++; |
3 | Werner | 30 | } |
38 | Werner | 31 | |
15 | Werner | 32 | /** get distance and direction between two points. |
38 | Werner | 33 | returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */ |
3 | Werner | 34 | float dist_and_direction(const QPointF &PStart, const QPointF &PEnd, float &rAngle) |
35 | { |
||
36 | float dx = PEnd.x() - PStart.x(); |
||
37 | float dy = PEnd.y() - PStart.y(); |
||
38 | float d = sqrt(dx*dx + dy*dy); |
||
39 | // direction: |
||
40 | rAngle = atan2(dx, dy); |
||
41 | return d; |
||
42 | } |
||
43 | |||
125 | Werner | 44 | /// dumps some core variables of a tree to a string. |
45 | QString Tree::dump() |
||
46 | { |
||
47 | QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8") |
||
48 | .arg(mId).arg(mSpecies->id()).arg(mDbh).arg(mHeight) |
||
49 | .arg(mPosition.x()).arg(mPosition.y()) |
||
50 | .arg(mRU->index()).arg(mLRI); |
||
51 | return result; |
||
52 | } |
||
3 | Werner | 53 | |
129 | Werner | 54 | void Tree::dumpList(DebugList &rTargetList) |
55 | { |
||
56 | rTargetList << mId << mSpecies->id() << mDbh << mHeight << mPosition.x() << mPosition.y() << mRU->index() << mLRI |
||
57 | << mStemMass << mRootMass << mLeafMass << mLeafArea; |
||
58 | } |
||
59 | |||
38 | Werner | 60 | void Tree::setup() |
61 | { |
||
106 | Werner | 62 | if (mDbh<=0 || mHeight<=0) |
38 | Werner | 63 | return; |
64 | // check stamp |
||
106 | Werner | 65 | Q_ASSERT_X(mSpecies!=0, "Tree::setup()", "species is NULL"); |
66 | mStamp = mSpecies->stamp(mDbh, mHeight); |
||
110 | Werner | 67 | |
68 | calcBiomassCompartments(); |
||
115 | Werner | 69 | mNPPReserve = mLeafMass; // initial value |
125 | Werner | 70 | mDbhDelta = 0.1; // initial value: used in growth() to estimate diameter increment |
110 | Werner | 71 | |
38 | Werner | 72 | } |
39 | Werner | 73 | |
110 | Werner | 74 | ////////////////////////////////////////////////// |
75 | //// Light functions (Pattern-stuff) |
||
76 | ////////////////////////////////////////////////// |
||
77 | |||
70 | Werner | 78 | #define NOFULLDBG |
77 | Werner | 79 | //#define NOFULLOPT |
39 | Werner | 80 | |
40 | Werner | 81 | |
77 | Werner | 82 | void Tree::applyStamp() |
83 | { |
||
106 | Werner | 84 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
77 | Werner | 85 | |
106 | Werner | 86 | QPoint pos = mGrid->indexAt(mPosition); |
87 | int offset = mStamp->offset(); |
||
77 | Werner | 88 | pos-=QPoint(offset, offset); |
89 | |||
90 | float local_dom; // height of Z* on the current position |
||
91 | int x,y; |
||
92 | float value; |
||
106 | Werner | 93 | int gr_stamp = mStamp->size(); |
77 | Werner | 94 | int grid_x, grid_y; |
95 | float *grid_value; |
||
106 | Werner | 96 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
77 | Werner | 97 | // todo: in this case we should use another algorithm!!! |
98 | return; |
||
99 | } |
||
100 | |||
101 | for (y=0;y<gr_stamp; ++y) { |
||
102 | grid_y = pos.y() + y; |
||
106 | Werner | 103 | grid_value = mGrid->ptr(pos.x(), grid_y); |
77 | Werner | 104 | for (x=0;x<gr_stamp;++x) { |
105 | // suppose there is no stamping outside |
||
106 | grid_x = pos.x() + x; |
||
107 | |||
106 | Werner | 108 | local_dom = mHeightGrid->valueAtIndex(grid_x/5, grid_y/5); |
109 | value = (*mStamp)(x,y); // stampvalue |
||
77 | Werner | 110 | value = 1. - value*lafactor / local_dom; // calculated value |
111 | value = qMax(value, 0.02f); // limit value |
||
112 | |||
113 | *grid_value++ *= value; |
||
114 | } |
||
115 | } |
||
116 | |||
117 | m_statPrint++; // count # of stamp applications... |
||
118 | } |
||
119 | |||
62 | Werner | 120 | |
74 | Werner | 121 | /** heightGrid() |
122 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid. |
||
123 | |||
124 | */ |
||
125 | void Tree::heightGrid() |
||
126 | { |
||
127 | // height of Z* |
||
106 | Werner | 128 | const float cellsize = mHeightGrid->cellsize(); |
74 | Werner | 129 | |
106 | Werner | 130 | QPoint p = mHeightGrid->indexAt(mPosition); // pos of tree on height grid |
131 | QPoint competition_grid = mGrid->indexAt(mPosition); |
||
74 | Werner | 132 | |
133 | int index_eastwest = competition_grid.x() % 5; // 4: very west, 0 east edge |
||
134 | int index_northsouth = competition_grid.y() % 5; // 4: northern edge, 0: southern edge |
||
135 | int dist[9]; |
||
136 | dist[3] = index_northsouth * 2 + 1; // south |
||
137 | dist[1] = index_eastwest * 2 + 1; // west |
||
138 | dist[5] = 10 - dist[3]; // north |
||
139 | dist[7] = 10 - dist[1]; // east |
||
140 | dist[8] = qMax(dist[5], dist[7]); // north-east |
||
141 | dist[6] = qMax(dist[3], dist[7]); // south-east |
||
142 | dist[0] = qMax(dist[3], dist[1]); // south-west |
||
143 | dist[2] = qMax(dist[5], dist[1]); // north-west |
||
75 | Werner | 144 | dist[4] = 0; // center cell |
76 | Werner | 145 | /* 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: |
146 | index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above. |
||
147 | e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2 |
||
148 | */ |
||
74 | Werner | 149 | |
150 | |||
106 | Werner | 151 | int ringcount = int(floor(mHeight / cellsize)) + 1; |
74 | Werner | 152 | int ix, iy; |
153 | int ring; |
||
154 | QPoint pos; |
||
155 | float hdom; |
||
156 | |||
157 | for (ix=-ringcount;ix<=ringcount;ix++) |
||
158 | for (iy=-ringcount; iy<=+ringcount; iy++) { |
||
159 | ring = qMax(abs(ix), abs(iy)); |
||
160 | QPoint pos(ix+p.x(), iy+p.y()); |
||
106 | Werner | 161 | if (mHeightGrid->isIndexValid(pos)) { |
162 | float &rHGrid = mHeightGrid->valueAtIndex(pos); |
||
163 | if (rHGrid > mHeight) // skip calculation if grid is higher than tree |
||
74 | Werner | 164 | continue; |
165 | int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y) |
||
106 | Werner | 166 | hdom = mHeight - dist[direction]; |
74 | Werner | 167 | if (ring>1) |
168 | hdom -= (ring-1)*10; |
||
169 | |||
170 | rHGrid = qMax(rHGrid, hdom); // write value |
||
171 | } // is valid |
||
172 | } // for (y) |
||
39 | Werner | 173 | } |
40 | Werner | 174 | |
175 | double Tree::readStamp() |
||
176 | { |
||
106 | Werner | 177 | if (!mStamp) |
51 | Werner | 178 | return 0.; |
106 | Werner | 179 | const Stamp *stamp = mStamp->reader(); |
40 | Werner | 180 | if (!stamp) |
181 | return 0.; |
||
106 | Werner | 182 | QPoint pos = mGrid->indexAt(mPosition); |
40 | Werner | 183 | int offset = stamp->offset(); |
184 | pos-=QPoint(offset, offset); |
||
185 | QPoint p; |
||
186 | |||
187 | int x,y; |
||
188 | double sum=0.; |
||
189 | for (x=0;x<stamp->size();++x) { |
||
190 | for (y=0;y<stamp->size(); ++y) { |
||
191 | p = pos + QPoint(x,y); |
||
106 | Werner | 192 | if (mGrid->isIndexValid(p)) |
193 | sum += mGrid->valueAtIndex(p) * (*stamp)(x,y); |
||
40 | Werner | 194 | } |
195 | } |
||
106 | Werner | 196 | float eigenvalue = mStamp->readSum() * lafactor; |
197 | mLRI = sum - eigenvalue;// additive |
||
198 | float dom_height = (*mHeightGrid)[mPosition]; |
||
53 | Werner | 199 | if (dom_height>0.) |
106 | Werner | 200 | mLRI = mLRI / dom_height; |
53 | Werner | 201 | |
202 | //mImpact = sum + eigenvalue;// multiplicative |
||
48 | Werner | 203 | // read dominance field... |
53 | Werner | 204 | |
106 | Werner | 205 | if (dom_height < mHeight) { |
48 | Werner | 206 | // if tree is higher than Z*, the dominance height |
207 | // a part of the crown is in "full light". |
||
208 | // total value = zstar/treeheight*value + 1-zstar/height |
||
209 | // reformulated to: |
||
106 | Werner | 210 | mLRI = mLRI * dom_height/mHeight ; |
48 | Werner | 211 | m_statAboveZ++; |
212 | } |
||
106 | Werner | 213 | if (fabs(mLRI < 0.000001)) |
214 | mLRI = 0.f; |
||
215 | qDebug() << "Tree #"<< id() << "value" << sum << "eigenvalue" << eigenvalue << "Impact" << mLRI; |
||
216 | return mLRI; |
||
40 | Werner | 217 | } |
218 | |||
58 | Werner | 219 | |
107 | Werner | 220 | void Tree::readStampMul() |
78 | Werner | 221 | { |
106 | Werner | 222 | if (!mStamp) |
107 | Werner | 223 | return; |
106 | Werner | 224 | const Stamp *reader = mStamp->reader(); |
78 | Werner | 225 | if (!reader) |
107 | Werner | 226 | return; |
106 | Werner | 227 | QPoint pos_reader = mGrid->indexAt(mPosition); |
78 | Werner | 228 | |
229 | int offset_reader = reader->offset(); |
||
106 | Werner | 230 | int offset_writer = mStamp->offset(); |
78 | Werner | 231 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
232 | |||
233 | QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer); |
||
234 | pos_reader-=QPoint(offset_reader, offset_reader); |
||
235 | QPoint p; |
||
236 | |||
237 | //float dom_height = (*m_dominanceGrid)[m_Position]; |
||
238 | float local_dom; |
||
239 | |||
240 | int x,y; |
||
241 | double sum=0.; |
||
242 | double value, own_value; |
||
243 | float *grid_value; |
||
244 | int reader_size = reader->size(); |
||
245 | int rx = pos_reader.x(); |
||
246 | int ry = pos_reader.y(); |
||
247 | for (y=0;y<reader_size; ++y, ++ry) { |
||
106 | Werner | 248 | grid_value = mGrid->ptr(rx, ry); |
78 | Werner | 249 | for (x=0;x<reader_size;++x) { |
250 | |||
251 | //p = pos_reader + QPoint(x,y); |
||
252 | //if (m_grid->isIndexValid(p)) { |
||
125 | Werner | 253 | local_dom = mHeightGrid->valueAtIndex((rx+x)/5, ry/5); // ry: gets ++ed in outer loop, rx not |
78 | Werner | 254 | //local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) ); |
125 | Werner | 255 | |
256 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*lafactor / local_dom; // old: dom_height; |
||
78 | Werner | 257 | own_value = qMax(own_value, 0.02); |
258 | value = *grid_value++ / own_value; // remove impact of focal tree |
||
259 | //if (value>0.) |
||
260 | sum += value * (*reader)(x,y); |
||
261 | |||
262 | //} // isIndexValid |
||
263 | } |
||
264 | } |
||
106 | Werner | 265 | mLRI = sum; |
78 | Werner | 266 | // read dominance field... |
267 | // this applies only if some trees are potentially *higher* than the dominant height grid |
||
268 | //if (dom_height < m_Height) { |
||
269 | // if tree is higher than Z*, the dominance height |
||
270 | // a part of the crown is in "full light". |
||
271 | // m_statAboveZ++; |
||
272 | // mImpact = 1. - (1. - mImpact)*dom_height/m_Height; |
||
273 | //} |
||
106 | Werner | 274 | if (mLRI > 1) |
275 | mLRI = 1.f; |
||
78 | Werner | 276 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact; |
110 | Werner | 277 | mRU->addWLA(mLRI * mLeafArea, mLeafArea); |
58 | Werner | 278 | } |
279 | |||
40 | Werner | 280 | void Tree::resetStatistics() |
281 | { |
||
282 | m_statPrint=0; |
||
105 | Werner | 283 | m_statCreated=0; |
48 | Werner | 284 | m_statAboveZ=0; |
40 | Werner | 285 | m_nextId=1; |
286 | } |
||
107 | Werner | 287 | |
110 | Werner | 288 | ////////////////////////////////////////////////// |
289 | //// Growth Functions |
||
290 | ////////////////////////////////////////////////// |
||
107 | Werner | 291 | |
110 | Werner | 292 | /// evaluate allometries and calculate LeafArea |
293 | void Tree::calcBiomassCompartments() |
||
294 | { |
||
295 | mLeafMass = mSpecies->biomassLeaf(mDbh); |
||
296 | mRootMass = mSpecies->biomassRoot(mDbh); |
||
297 | mStemMass = mSpecies->biomassStem(mDbh); |
||
112 | Werner | 298 | // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg] |
110 | Werner | 299 | mLeafArea = mLeafMass * mSpecies->specificLeafArea(); |
300 | } |
||
301 | |||
302 | |||
107 | Werner | 303 | void Tree::grow() |
304 | { |
||
113 | Werner | 305 | // step 1: get radiation from ressource unit: radiaton (MJ/tree/year) total intercepted radiation for this tree per year! |
107 | Werner | 306 | double radiation = mRU->interceptedRadiation(mLRI * mLeafArea); |
113 | Werner | 307 | // step 2: get fraction of PARutilized, i.e. fraction of intercepted rad that is utiliziable (per year) |
107 | Werner | 308 | |
115 | Werner | 309 | double raw_gpp_per_rad = mRU->ressourceUnitSpecies(mSpecies).prod3PG().GPPperRad(); |
113 | Werner | 310 | // GPP (without aging-effect) [gC] / year |
311 | double raw_gpp = raw_gpp_per_rad * radiation; |
||
312 | /* |
||
313 | if (mRU->index()==3) { |
||
314 | qDebug() << "tree production: radiation: " << radiation << "gpp/rad:" << raw_gpp_per_rad << "gpp" << raw_gpp << "LRI:" << mLRI << "LeafArea:" << mLeafArea; |
||
315 | }*/ |
||
115 | Werner | 316 | // apply aging |
317 | double gpp = raw_gpp * 0.6; // aging |
||
318 | double npp = gpp * 0.47; // respiration loss |
||
113 | Werner | 319 | |
115 | Werner | 320 | partitioning(npp); |
321 | |||
125 | Werner | 322 | mStamp = mSpecies->stamp(mDbh, mHeight); // get new stamp for updated dimensions |
110 | Werner | 323 | |
107 | Werner | 324 | } |
325 | |||
117 | Werner | 326 | |
327 | // just used to test the DBG_IF_x macros... |
||
328 | QString test_cntr() |
||
329 | { |
||
330 | static int cnt = 0; |
||
331 | cnt++; |
||
332 | return QString::number(cnt); |
||
333 | } |
||
334 | |||
118 | Werner | 335 | #if !defined(DBGMODE) |
336 | # ifndef QT_NO_DEBUG |
||
337 | # define DBGMODE(stmts) { stmts } |
||
338 | # else |
||
339 | # define DBGMODE(stmts) qt_noop() |
||
340 | # endif |
||
341 | #endif |
||
342 | |||
115 | Werner | 343 | void Tree::partitioning(double npp) |
344 | { |
||
119 | Werner | 345 | DBGMODE( |
346 | if (mId==1) |
||
347 | test_cntr(); |
||
348 | ); |
||
115 | Werner | 349 | double harshness = mRU->ressourceUnitSpecies(mSpecies).prod3PG().harshness(); |
350 | // add content of reserve pool |
||
116 | Werner | 351 | npp += mNPPReserve; |
119 | Werner | 352 | const double reserve_size = mLeafMass; |
353 | |||
117 | Werner | 354 | double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) |
355 | // turnover rates |
||
356 | const double &to_fol = mSpecies->turnoverLeaf(); |
||
357 | const double &to_wood = mSpecies->turnoverStem(); |
||
358 | const double &to_root = mSpecies->turnoverRoot(); |
||
116 | Werner | 359 | |
117 | Werner | 360 | apct_root = harshness; |
361 | |||
362 | double b_wf = 1.32; // ratio of allometric exponents... now fixed |
||
363 | |||
364 | // Duursma 2007, Eq. (20) |
||
365 | apct_wood = (mLeafMass * to_wood / npp + b_wf*(1.-apct_root) - b_wf * to_fol/npp) / ( mStemMass / mStemMass + b_wf ); |
||
366 | apct_foliage = 1. - apct_root - apct_wood; |
||
367 | |||
119 | Werner | 368 | // senescence demands of the compartemnts |
369 | double senescence_foliage = mLeafMass * to_fol; |
||
370 | double senescence_stem = mStemMass * to_wood; |
||
371 | double senescence_root = mRootMass * to_root; |
||
372 | |||
373 | // brutto compartment increments |
||
374 | double gross_foliage = npp * apct_foliage; |
||
375 | double gross_stem = npp * apct_wood; |
||
376 | double gross_root = npp * apct_root; |
||
377 | |||
378 | // netto increments |
||
379 | double net_foliage = gross_foliage - senescence_foliage; |
||
380 | double net_root = gross_root - senescence_root; |
||
381 | double net_stem = gross_stem - senescence_stem; |
||
382 | |||
383 | // flow back to reserve pool: |
||
384 | double to_reserve = qMin(reserve_size, net_stem); |
||
385 | net_stem -= to_reserve; |
||
386 | |||
129 | Werner | 387 | /*DBG_IF_X(mId == 1 , "Tree::partitioning", "dump", dump() |
125 | Werner | 388 | + 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") |
389 | .arg(npp).arg(senescence_foliage).arg(senescence_stem).arg(senescence_root) |
||
129 | Werner | 390 | .arg(net_foliage).arg(net_stem).arg(net_root).arg(to_reserve).arg(mNPPReserve) );*/ |
125 | Werner | 391 | |
129 | Werner | 392 | DBGMODE( |
393 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreePartition) && mId<300) { |
||
394 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreePartition); |
||
395 | dumpList(out); // add tree headers |
||
396 | out << npp << senescence_foliage << senescence_stem << senescence_root << net_foliage << net_stem << net_root << to_reserve << mNPPReserve; |
||
397 | } |
||
398 | ); // DBGMODE() |
||
399 | |||
119 | Werner | 400 | // update of compartments |
401 | mLeafMass += net_foliage; |
||
125 | Werner | 402 | mLeafMass = qMax(mLeafMass, 0.f); |
403 | mLeafArea = mLeafMass * mSpecies->specificLeafArea(); |
||
404 | |||
119 | Werner | 405 | mRootMass += net_root; |
125 | Werner | 406 | mRootMass = qMax(mRootMass, 0.f); |
119 | Werner | 407 | |
125 | Werner | 408 | // calculate the dimensions of growth (diameter, height) |
119 | Werner | 409 | grow_diameter(net_stem); |
410 | |||
125 | Werner | 411 | // calculate stem biomass using the allometric equation |
412 | mStemMass = mSpecies->biomassStem(mDbh); |
||
121 | Werner | 413 | |
115 | Werner | 414 | } |
415 | |||
125 | Werner | 416 | |
417 | /** Determination of diamter and height growth based on increment of the stem mass (@net_stem_npp). |
||
418 | Refer to XXX for equations and variables. |
||
419 | This function updates the dbh and height of the tree. |
||
420 | */ |
||
119 | Werner | 421 | inline void Tree::grow_diameter(const double &net_stem_npp) |
422 | { |
||
423 | // determine dh-ratio of increment |
||
424 | // height increment is a function of light competition: |
||
125 | Werner | 425 | double hd_growth = relative_height_growth(); // hd of height growth |
426 | //DBG_IF_X(rel_height_growth<0 || rel_height_growth>1., "Tree::grow_dimater", "rel_height_growth out of bound.", dump()); |
||
115 | Werner | 427 | |
125 | Werner | 428 | const double volume_factor = mSpecies->volumeFactor(); |
123 | Werner | 429 | |
125 | Werner | 430 | double factor_diameter = 1. / ( volume_factor * (mDbh + mDbhDelta)*(mDbh + mDbhDelta) * ( 2. * mHeight/mDbh + hd_growth) ); |
431 | double delta_d_estimate = factor_diameter * net_stem_npp; // estimated dbh-inc using last years increment |
||
432 | |||
433 | // using that dbh-increment we estimate a stem-mass-increment and the residual (Eq. 9) |
||
434 | double stem_estimate = volume_factor * (mDbh + delta_d_estimate)*(mDbh + delta_d_estimate)*(mHeight + delta_d_estimate*hd_growth); |
||
435 | double stem_residual = stem_estimate - (mStemMass + net_stem_npp); |
||
436 | |||
437 | // the final increment is then: |
||
438 | double d_increment = factor_diameter * (net_stem_npp - stem_residual); // Eq. (11) |
||
439 | DBG_IF_X(mId == 1 || d_increment<0., "Tree::grow_dimater", "increment < 0.", dump() |
||
440 | + QString("\nhdz %1 factor_diameter %2 stem_residual %3 delta_d_estimate %4 d_increment %5 final residual(kg) %6") |
||
441 | .arg(hd_growth).arg(factor_diameter).arg(stem_residual).arg(delta_d_estimate).arg(d_increment) |
||
442 | .arg( volume_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((mStemMass + net_stem_npp)) )); |
||
443 | |||
444 | DBGMODE( |
||
445 | double res_final = volume_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((mStemMass + net_stem_npp)); |
||
446 | DBG_IF_X(res_final > 1, "Tree::grow_diameter", "final residual stem estimate > 1kg", dump()); |
||
447 | DBG_IF_X(d_increment > 10. || d_increment*hd_growth/100. >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()); |
||
448 | //dbgstruct["sen_demand"]=sen_demand; |
||
129 | Werner | 449 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeGrowth) && mId<300) { |
126 | Werner | 450 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeGrowth); |
129 | Werner | 451 | dumpList(out); // add tree headers |
452 | out << net_stem_npp << hd_growth << factor_diameter << delta_d_estimate << d_increment; |
||
126 | Werner | 453 | } |
125 | Werner | 454 | |
129 | Werner | 455 | ); // DBGMODE() |
125 | Werner | 456 | d_increment = qMax(d_increment, 0.); |
457 | |||
458 | // update state variables |
||
459 | mDbh += d_increment; |
||
460 | mDbhDelta = d_increment; // save for next year's growth |
||
461 | mHeight += d_increment * hd_growth * 0.01; |
||
119 | Werner | 462 | } |
463 | |||
125 | Werner | 464 | |
465 | /// return the HD ratio of this year's increment based on the light status. |
||
119 | Werner | 466 | inline double Tree::relative_height_growth() |
467 | { |
||
468 | double hd_low, hd_high; |
||
469 | mSpecies->hdRange(mDbh, hd_low, hd_high); |
||
470 | |||
125 | Werner | 471 | DBG_IF_X(hd_low>hd_high, "Tree::relative_height_growth", "hd low higher dann hd_high for ", dump()); |
472 | 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)); |
||
473 | |||
474 | // scale according to LRI: if receiving much light (LRI=1), the result is hd_low (for open grown trees) |
||
475 | double hd_ratio = hd_high - (hd_high-hd_low)*mLRI; |
||
476 | return hd_ratio; |
||
119 | Werner | 477 | } |