Rev 1118 | Rev 1159 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
| Rev 1118 | Rev 1158 | ||
|---|---|---|---|
| Line 40... | Line 40... | ||
| 40 | HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid(); |
40 | HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid(); |
| 41 | FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid(); |
41 | FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid(); |
| 42 | 42 | ||
| 43 | QPoint imap = ru->cornerPointOffset(); // offset on LIF/saplings grid |
43 | QPoint imap = ru->cornerPointOffset(); // offset on LIF/saplings grid |
| 44 | QPoint iseedmap = QPoint(imap.x()/10, imap.y()/10); // seed-map has 20m resolution, LIF 2m -> factor 10 |
44 | QPoint iseedmap = QPoint(imap.x()/10, imap.y()/10); // seed-map has 20m resolution, LIF 2m -> factor 10 |
| - | 45 | ||
| - | 46 | for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i) |
|
| - | 47 | (*i)->saplingStat().clearStatistics(); |
|
| - | 48 | ||
| - | 49 | double lif_corr[cPxPerHectare]; |
|
| - | 50 | for (int i=0;i<cPxPerHectare;++i) |
|
| - | 51 | lif_corr[i]=-1.; |
|
| 45 | 52 | ||
| 46 | int species_idx = irandom(0, ru->ruSpecies().size()-1); |
53 | int species_idx = irandom(0, ru->ruSpecies().size()-1); |
| 47 | for (int s_idx = 0; s_idx<ru->ruSpecies().size(); ++s_idx) { |
54 | for (int s_idx = 0; s_idx<ru->ruSpecies().size(); ++s_idx) { |
| 48 | 55 | ||
| 49 | // start from a random species (and cycle through the available species)
|
56 | // start from a random species (and cycle through the available species)
|
| Line 73... | Line 80... | ||
| 73 | int isc = 0; // index on 2m cell |
80 | int isc = 0; // index on 2m cell |
| 74 | for (int iy=0; iy<cPxPerRU; ++iy) { |
81 | for (int iy=0; iy<cPxPerRU; ++iy) { |
| 75 | s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row |
82 | s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row |
| 76 | isc = mGrid.index(imap.x(), imap.y()+iy); |
83 | isc = mGrid.index(imap.x(), imap.y()+iy); |
| 77 | 84 | ||
| 78 | for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc, ++mTested) { |
- | |
| - | 85 | for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) { |
|
| 79 | if (s->state == SaplingCell::CellFree) { |
86 | if (s->state == SaplingCell::CellFree) { |
| 80 | bool viable = true; |
- | |
| 81 | // is a sapling of the current species already on the pixel?
|
87 | // is a sapling of the current species already on the pixel?
|
| 82 | // * test for sapling height already in cell state
|
88 | // * test for sapling height already in cell state
|
| 83 | // * test for grass-cover already in cell state
|
89 | // * test for grass-cover already in cell state
|
| 84 | int i_occupied = -1; |
- | |
| 85 | for (int i=0;i<NSAPCELLS;++i) { |
- | |
| 86 | if (!s->saplings[i].is_occupied() && i_occupied<0) |
- | |
| 87 | i_occupied=i; |
- | |
| 88 | if (s->saplings[i].species_index == species_idx) { |
- | |
| 89 | viable = false; |
- | |
| - | 90 | SaplingTree *stree=0; |
|
| - | 91 | SaplingTree *slot=s->saplings; |
|
| - | 92 | for (int i=0;i<NSAPCELLS;++i, ++slot) { |
|
| - | 93 | if (!stree && !slot->is_occupied()) |
|
| - | 94 | stree=slot; |
|
| - | 95 | if (slot->species_index == species_idx) { |
|
| - | 96 | stree=0; |
|
| - | 97 | break; |
|
| 90 | }
|
98 | }
|
| 91 | }
|
99 | }
|
| 92 | 100 | ||
| 93 | if (viable && i_occupied>=0) { |
- | |
| - | 101 | if (stree) { |
|
| 94 | // grass cover?
|
102 | // grass cover?
|
| 95 | DBG_IF(i_occupied<0, "establishment", "invalid value i_occupied<0"); |
- | |
| 96 | float seed_map_value = seedmap[mGrid.index10(isc)]; |
103 | float seed_map_value = seedmap[mGrid.index10(isc)]; |
| 97 | if (seed_map_value==0.f) |
104 | if (seed_map_value==0.f) |
| 98 | continue; |
105 | continue; |
| 99 | const HeightGridValue &hgv = (*height_grid)[mGrid.index5(isc)]; |
106 | const HeightGridValue &hgv = (*height_grid)[mGrid.index5(isc)]; |
| 100 | float lif_value = (*lif_grid)[isc]; |
107 | float lif_value = (*lif_grid)[isc]; |
| 101 | double lif_corrected = rus->species()->speciesSet()->LRIcorrection(lif_value, 4. / hgv.height); |
- | |
| - | 108 | ||
| - | 109 | double &lif_corrected = lif_corr[iy*cPxPerRU+ix]; |
|
| - | 110 | // calculate the LIFcorrected only once per pixel
|
|
| - | 111 | if (lif_corrected<0.) |
|
| - | 112 | lif_corrected = rus->species()->speciesSet()->LRIcorrection(lif_value, 4. / hgv.height); |
|
| - | 113 | ||
| 102 | // check for the combination of seed availability and light on the forest floor
|
114 | // check for the combination of seed availability and light on the forest floor
|
| 103 | if (drandom() < seed_map_value*lif_corrected*abiotic_env ) { |
- | |
| 104 | // ok, lets add a sapling at the given position
|
- | |
| 105 | s->saplings[i_occupied].setSapling(0.05f, 1, species_idx); |
- | |
| 106 | s->checkState(); |
- | |
| 107 | mAdded++; |
- | |
| - | 115 | if (drandom() < seed_map_value*lif_corrected*abiotic_env ) { |
|
| - | 116 | // ok, lets add a sapling at the given position (age is incremented later)
|
|
| - | 117 | stree->setSapling(0.05f, 0, species_idx); |
|
| - | 118 | s->checkState(); |
|
| - | 119 | rus->saplingStat().mAdded++; |
|
| 108 | 120 | ||
| 109 | }
|
- | |
| - | 121 | }
|
|
| 110 | 122 | ||
| 111 | }
|
123 | }
|
| 112 | 124 | ||
| 113 | }
|
125 | }
|
| 114 | }
|
126 | }
|
| Line 120... | Line 132... | ||
| 120 | 132 | ||
| 121 | void Saplings::saplingGrowth(const ResourceUnit *ru) |
133 | void Saplings::saplingGrowth(const ResourceUnit *ru) |
| 122 | {
|
134 | {
|
| 123 | HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid(); |
135 | HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid(); |
| 124 | FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid(); |
136 | FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid(); |
| 125 | - | ||
| 126 | for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i) |
- | |
| 127 | (*i)->saplingStat().clearStatistics(); |
- | |
| 128 | 137 | ||
| 129 | QPoint imap = mGrid.indexAt(ru->boundingBox().topLeft()); |
138 | QPoint imap = mGrid.indexAt(ru->boundingBox().topLeft()); |
| 130 | bool need_check=false; |
139 | bool need_check=false; |
| 131 | for (int iy=0; iy<cPxPerRU; ++iy) { |
140 | for (int iy=0; iy<cPxPerRU; ++iy) { |
| 132 | SaplingCell *s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row |
141 | SaplingCell *s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row |
| Line 148... | Line 157... | ||
| 148 | s->checkState(); |
157 | s->checkState(); |
| 149 | }
|
158 | }
|
| 150 | }
|
159 | }
|
| 151 | }
|
160 | }
|
| 152 | 161 | ||
| - | 162 | ||
| - | 163 | ||
| - | 164 | ||
| - | 165 | // store statistics on saplings/regeneration
|
|
| - | 166 | for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i) { |
|
| - | 167 | (*i)->saplingStat().calculate((*i)->species(), const_cast<ResourceUnit*>(ru)); |
|
| - | 168 | (*i)->statistics().add(&((*i)->saplingStat())); |
|
| - | 169 | }
|
|
| 153 | }
|
170 | }
|
| 154 | 171 | ||
| 155 | void Saplings::updateBrowsingPressure() |
172 | void Saplings::updateBrowsingPressure() |
| 156 | {
|
173 | {
|
| 157 | if (GlobalSettings::instance()->settings().valueBool("model.settings.browsing.enabled")) |
174 | if (GlobalSettings::instance()->settings().valueBool("model.settings.browsing.enabled")) |
| Line 203... | Line 220... | ||
| 203 | tree.stress_years++; |
220 | tree.stress_years++; |
| 204 | if (tree.stress_years > species->saplingGrowthParameters().maxStressYears) { |
221 | if (tree.stress_years > species->saplingGrowthParameters().maxStressYears) { |
| 205 | // sapling dies...
|
222 | // sapling dies...
|
| 206 | tree.clear(); |
223 | tree.clear(); |
| 207 | rus->saplingStat().addCarbonOfDeadSapling( tree.height / species->saplingGrowthParameters().hdSapling * 100.f ); |
224 | rus->saplingStat().addCarbonOfDeadSapling( tree.height / species->saplingGrowthParameters().hdSapling * 100.f ); |
| - | 225 | rus->saplingStat().mDied++; |
|
| 208 | return true; // need cleanup |
226 | return true; // need cleanup |
| 209 | }
|
227 | }
|
| 210 | } else { |
228 | } else { |
| 211 | tree.stress_years=0; // reset stress counter |
229 | tree.stress_years=0; // reset stress counter |
| 212 | }
|
230 | }
|
| Line 234... | Line 252... | ||
| 234 | for (int i=0;i<to_establish;i++) { |
252 | for (int i=0;i<to_establish;i++) { |
| 235 | Tree &bigtree = const_cast<ResourceUnit*>(ru)->newTree(); |
253 | Tree &bigtree = const_cast<ResourceUnit*>(ru)->newTree(); |
| 236 | 254 | ||
| 237 | bigtree.setPosition(mGrid.indexOf(isc)); |
255 | bigtree.setPosition(mGrid.indexOf(isc)); |
| 238 | // add variation: add +/-10% to dbh and *independently* to height.
|
256 | // add variation: add +/-10% to dbh and *independently* to height.
|
| 239 | bigtree.setDbh(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)); |
- | |
| 240 | bigtree.setHeight(tree.height * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)); |
- | |
| - | 257 | bigtree.setDbh(static_cast<float>(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation))); |
|
| - | 258 | bigtree.setHeight(static_cast<float>(tree.height * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation))); |
|
| 241 | bigtree.setSpecies( const_cast<Species*>(species) ); |
259 | bigtree.setSpecies( const_cast<Species*>(species) ); |
| 242 | bigtree.setAge(tree.age,tree.height); |
260 | bigtree.setAge(tree.age,tree.height); |
| 243 | bigtree.setRU(const_cast<ResourceUnit*>(ru)); |
261 | bigtree.setRU(const_cast<ResourceUnit*>(ru)); |
| 244 | bigtree.setup(); |
262 | bigtree.setup(); |
| 245 | const Tree *t = &bigtree; |
263 | const Tree *t = &bigtree; |
| Line 272... | Line 290... | ||
| 272 | mRecruited=mDied=mLiving=0; |
290 | mRecruited=mDied=mLiving=0; |
| 273 | mSumDbhDied=0.; |
291 | mSumDbhDied=0.; |
| 274 | mAvgHeight=0.; |
292 | mAvgHeight=0.; |
| 275 | mAvgAge=0.; |
293 | mAvgAge=0.; |
| 276 | mAvgDeltaHPot=mAvgHRealized=0.; |
294 | mAvgDeltaHPot=mAvgHRealized=0.; |
| - | 295 | mAdded=0; |
|
| - | 296 | ||
| - | 297 | }
|
|
| - | 298 | ||
| - | 299 | void SaplingStat::calculate(const Species *species, ResourceUnit *ru) |
|
| - | 300 | {
|
|
| - | 301 | if (mLiving) { |
|
| - | 302 | mAvgHeight /= double(mLiving); |
|
| - | 303 | mAvgAge /= double(mLiving); |
|
| - | 304 | mAvgDeltaHPot /= double(mLiving); |
|
| - | 305 | mAvgHRealized /= double(mLiving); |
|
| - | 306 | }
|
|
| - | 307 | ||
| - | 308 | // calculate carbon balance
|
|
| - | 309 | CNPair old_state = mCarbonLiving; |
|
| - | 310 | mCarbonLiving.clear(); |
|
| - | 311 | ||
| - | 312 | CNPair dead_wood, dead_fine; // pools for mortality |
|
| - | 313 | // average dbh
|
|
| - | 314 | if (mLiving>0) { |
|
| - | 315 | // calculate the avg dbh and number of stems
|
|
| - | 316 | double avg_dbh = mAvgHeight / species->saplingGrowthParameters().hdSapling * 100.; |
|
| - | 317 | double n = mLiving * species->saplingGrowthParameters().representedStemNumber( avg_dbh ); |
|
| - | 318 | // woody parts: stem, branchse and coarse roots
|
|
| - | 319 | double woody_bm = species->biomassWoody(avg_dbh) + species->biomassBranch(avg_dbh) + species->biomassRoot(avg_dbh); |
|
| - | 320 | double foliage = species->biomassFoliage(avg_dbh); |
|
| - | 321 | double fineroot = foliage*species->finerootFoliageRatio(); |
|
| - | 322 | ||
| - | 323 | mCarbonLiving.addBiomass( woody_bm*n, species->cnWood() ); |
|
| - | 324 | mCarbonLiving.addBiomass( foliage*n, species->cnFoliage() ); |
|
| - | 325 | mCarbonLiving.addBiomass( fineroot*n, species->cnFineroot() ); |
|
| - | 326 | ||
| - | 327 | // turnover
|
|
| - | 328 | if (ru->snag()) |
|
| - | 329 | ru->snag()->addTurnoverLitter(species, foliage*species->turnoverLeaf(), fineroot*species->turnoverRoot()); |
|
| - | 330 | ||
| - | 331 | // calculate the "mortality from competition", i.e. carbon that stems from reduction of stem numbers
|
|
| - | 332 | // from Reinekes formula.
|
|
| - | 333 | //
|
|
| - | 334 | if (avg_dbh>1.) { |
|
| - | 335 | double avg_dbh_before = (mAvgHeight - mAvgHRealized) / species->saplingGrowthParameters().hdSapling * 100.; |
|
| - | 336 | double n_before = mLiving * species->saplingGrowthParameters().representedStemNumber( qMax(1.,avg_dbh_before) ); |
|
| - | 337 | if (n<n_before) { |
|
| - | 338 | dead_wood.addBiomass( woody_bm * (n_before-n), species->cnWood() ); |
|
| - | 339 | dead_fine.addBiomass( foliage * (n_before-n), species->cnFoliage() ); |
|
| - | 340 | dead_fine.addBiomass( fineroot * (n_before-n), species->cnFineroot() ); |
|
| - | 341 | }
|
|
| - | 342 | }
|
|
| - | 343 | ||
| - | 344 | }
|
|
| - | 345 | if (mDied) { |
|
| - | 346 | double avg_dbh_dead = mSumDbhDied / double(mDied); |
|
| - | 347 | double n = mDied * species->saplingGrowthParameters().representedStemNumber( avg_dbh_dead ); |
|
| - | 348 | // woody parts: stem, branchse and coarse roots
|
|
| - | 349 | ||
| - | 350 | dead_wood.addBiomass( ( species->biomassWoody(avg_dbh_dead) + species->biomassBranch(avg_dbh_dead) + species->biomassRoot(avg_dbh_dead)) * n, species->cnWood() ); |
|
| - | 351 | double foliage = species->biomassFoliage(avg_dbh_dead)*n; |
|
| - | 352 | ||
| - | 353 | dead_fine.addBiomass( foliage, species->cnFoliage() ); |
|
| - | 354 | dead_fine.addBiomass( foliage*species->finerootFoliageRatio(), species->cnFineroot() ); |
|
| - | 355 | }
|
|
| - | 356 | if (!dead_wood.isEmpty() || !dead_fine.isEmpty()) |
|
| - | 357 | if (ru->snag()) |
|
| - | 358 | ru->snag()->addToSoil(species, dead_wood, dead_fine); |
|
| - | 359 | ||
| - | 360 | // calculate net growth:
|
|
| - | 361 | // delta of stocks
|
|
| - | 362 | mCarbonGain = mCarbonLiving + dead_fine + dead_wood - old_state; |
|
| - | 363 | if (mCarbonGain.C < 0) |
|
| - | 364 | mCarbonGain.clear(); |
|
| - | 365 | ||
| - | 366 | ||
| - | 367 | GlobalSettings::instance()->systemStatistics()->saplingCount+=mLiving; |
|
| - | 368 | GlobalSettings::instance()->systemStatistics()->newSaplings+=mAdded; |
|
| 277 | 369 | ||
| 278 | }
|
370 | }
|