Rev 1118 | Rev 1159 | Go to most recent revision | Only display areas with differences | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 1118 | Rev 1158 | ||
---|---|---|---|
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/saplings.cpp': |
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/saplings.cpp': |
2 | #include "global.h"
|
2 | #include "global.h"
|
3 | #include "saplings.h"
|
3 | #include "saplings.h"
|
4 | 4 | ||
5 | #include "globalsettings.h"
|
5 | #include "globalsettings.h"
|
6 | #include "model.h"
|
6 | #include "model.h"
|
7 | #include "resourceunit.h"
|
7 | #include "resourceunit.h"
|
8 | #include "resourceunitspecies.h"
|
8 | #include "resourceunitspecies.h"
|
9 | #include "establishment.h"
|
9 | #include "establishment.h"
|
10 | #include "species.h"
|
10 | #include "species.h"
|
11 | #include "seeddispersal.h"
|
11 | #include "seeddispersal.h"
|
12 | 12 | ||
13 | double Saplings::mRecruitmentVariation = 0.1; // +/- 10% |
13 | double Saplings::mRecruitmentVariation = 0.1; // +/- 10% |
14 | double Saplings::mBrowsingPressure = 0.; |
14 | double Saplings::mBrowsingPressure = 0.; |
15 | 15 | ||
16 | 16 | ||
17 | Saplings::Saplings() |
17 | Saplings::Saplings() |
18 | {
|
18 | {
|
19 | 19 | ||
20 | }
|
20 | }
|
21 | 21 | ||
22 | void Saplings::setup() |
22 | void Saplings::setup() |
23 | {
|
23 | {
|
24 | mGrid.setup(GlobalSettings::instance()->model()->grid()->metricRect(), GlobalSettings::instance()->model()->grid()->cellsize()); |
24 | mGrid.setup(GlobalSettings::instance()->model()->grid()->metricRect(), GlobalSettings::instance()->model()->grid()->cellsize()); |
25 | 25 | ||
26 | // mask out out-of-project areas
|
26 | // mask out out-of-project areas
|
27 | HeightGrid *hg = GlobalSettings::instance()->model()->heightGrid(); |
27 | HeightGrid *hg = GlobalSettings::instance()->model()->heightGrid(); |
28 | for (int i=0;i<mGrid.count();++i) { |
28 | for (int i=0;i<mGrid.count();++i) { |
29 | if (!hg->valueAtIndex(mGrid.index5(i)).isValid()) |
29 | if (!hg->valueAtIndex(mGrid.index5(i)).isValid()) |
30 | mGrid[i].state = SaplingCell::CellInvalid; |
30 | mGrid[i].state = SaplingCell::CellInvalid; |
31 | else
|
31 | else
|
32 | mGrid[i].state = SaplingCell::CellFree; |
32 | mGrid[i].state = SaplingCell::CellFree; |
33 | }
|
33 | }
|
34 | 34 | ||
35 | 35 | ||
36 | }
|
36 | }
|
37 | 37 | ||
38 | void Saplings::establishment(const ResourceUnit *ru) |
38 | void Saplings::establishment(const ResourceUnit *ru) |
39 | {
|
39 | {
|
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)
|
50 | species_idx = ++species_idx % ru->ruSpecies().size(); |
57 | species_idx = ++species_idx % ru->ruSpecies().size(); |
51 | 58 | ||
52 | ResourceUnitSpecies *rus = ru->ruSpecies()[species_idx]; |
59 | ResourceUnitSpecies *rus = ru->ruSpecies()[species_idx]; |
53 | // check if there are seeds of the given species on the resource unit
|
60 | // check if there are seeds of the given species on the resource unit
|
54 | float seeds = 0.f; |
61 | float seeds = 0.f; |
55 | Grid<float> &seedmap = const_cast<Grid<float>& >(rus->species()->seedDispersal()->seedMap()); |
62 | Grid<float> &seedmap = const_cast<Grid<float>& >(rus->species()->seedDispersal()->seedMap()); |
56 | for (int iy=0;iy<5;++iy) { |
63 | for (int iy=0;iy<5;++iy) { |
57 | float *p = seedmap.ptr(iseedmap.x(), iseedmap.y()); |
64 | float *p = seedmap.ptr(iseedmap.x(), iseedmap.y()); |
58 | for (int ix=0;ix<5;++ix) |
65 | for (int ix=0;ix<5;++ix) |
59 | seeds += *p++; |
66 | seeds += *p++; |
60 | }
|
67 | }
|
61 | // if there are no seeds: no need to do more
|
68 | // if there are no seeds: no need to do more
|
62 | if (seeds==0.f) |
69 | if (seeds==0.f) |
63 | continue; |
70 | continue; |
64 | 71 | ||
65 | // calculate the abiotic environment (TACA)
|
72 | // calculate the abiotic environment (TACA)
|
66 | rus->establishment().calculateAbioticEnvironment(); |
73 | rus->establishment().calculateAbioticEnvironment(); |
67 | double abiotic_env = rus->establishment().abioticEnvironment(); |
74 | double abiotic_env = rus->establishment().abioticEnvironment(); |
68 | if (abiotic_env==0.) |
75 | if (abiotic_env==0.) |
69 | continue; |
76 | continue; |
70 | 77 | ||
71 | // loop over all 2m cells on this resource unit
|
78 | // loop over all 2m cells on this resource unit
|
72 | SaplingCell *s; |
79 | SaplingCell *s; |
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, ++ |
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 |
|
90 | SaplingTree *stree=0; |
- | 91 | SaplingTree *slot=s->saplings; |
|
85 | for (int i=0;i<NSAPCELLS;++ |
92 | for (int i=0;i<NSAPCELLS;++i, ++slot) { |
86 | if (! |
93 | if (!stree && !slot->is_occupied()) |
87 |
|
94 | stree=slot; |
88 | if ( |
95 | if (slot->species_index == species_idx) { |
- | 96 | stree=0; |
|
89 |
|
97 | break; |
90 | }
|
98 | }
|
91 | }
|
99 | }
|
92 | 100 | ||
93 | if ( |
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]; |
- | 108 | ||
- | 109 | double &lif_corrected = lif_corr[iy*cPxPerRU+ix]; |
|
- | 110 | // calculate the LIFcorrected only once per pixel
|
|
- | 111 | if (lif_corrected<0.) |
|
101 |
|
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 ) { |
115 | if (drandom() < seed_map_value*lif_corrected*abiotic_env ) { |
104 |
|
116 | // ok, lets add a sapling at the given position (age is incremented later)
|
105 |
|
117 | stree->setSapling(0.05f, 0, species_idx); |
106 | s->checkState(); |
118 | s->checkState(); |
107 |
|
119 | rus->saplingStat().mAdded++; |
108 | 120 | ||
109 | }
|
121 | }
|
110 | 122 | ||
111 | }
|
123 | }
|
112 | 124 | ||
113 | }
|
125 | }
|
114 | }
|
126 | }
|
115 | }
|
127 | }
|
116 | 128 | ||
117 | }
|
129 | }
|
118 | 130 | ||
119 | }
|
131 | }
|
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 |
133 | int isc = mGrid.index(imap.x(), imap.y()+iy); |
142 | int isc = mGrid.index(imap.x(), imap.y()+iy); |
134 | 143 | ||
135 | for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) { |
144 | for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) { |
136 | if (s->state != SaplingCell::CellInvalid) { |
145 | if (s->state != SaplingCell::CellInvalid) { |
137 | need_check=false; |
146 | need_check=false; |
138 | for (int i=0;i<NSAPCELLS;++i) { |
147 | for (int i=0;i<NSAPCELLS;++i) { |
139 | if (s->saplings[i].is_occupied()) { |
148 | if (s->saplings[i].is_occupied()) { |
140 | // growth of this sapling tree
|
149 | // growth of this sapling tree
|
141 | const HeightGridValue &hgv = (*height_grid)[height_grid->index5(isc)]; |
150 | const HeightGridValue &hgv = (*height_grid)[height_grid->index5(isc)]; |
142 | float lif_value = (*lif_grid)[isc]; |
151 | float lif_value = (*lif_grid)[isc]; |
143 | 152 | ||
144 | need_check |= growSapling(ru, s->saplings[i], isc, hgv.height, lif_value); |
153 | need_check |= growSapling(ru, s->saplings[i], isc, hgv.height, lif_value); |
145 | }
|
154 | }
|
146 | }
|
155 | }
|
147 | if (need_check) |
156 | if (need_check) |
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")) |
158 | Saplings::mBrowsingPressure = GlobalSettings::instance()->settings().valueDouble("model.settings.browsing.browsingPressure"); |
175 | Saplings::mBrowsingPressure = GlobalSettings::instance()->settings().valueDouble("model.settings.browsing.browsingPressure"); |
159 | else
|
176 | else
|
160 | Saplings::mBrowsingPressure = 0.; |
177 | Saplings::mBrowsingPressure = 0.; |
161 | }
|
178 | }
|
162 | 179 | ||
163 | bool Saplings::growSapling(const ResourceUnit *ru, SaplingTree &tree, int isc, float dom_height, float lif_value) |
180 | bool Saplings::growSapling(const ResourceUnit *ru, SaplingTree &tree, int isc, float dom_height, float lif_value) |
164 | {
|
181 | {
|
165 | ResourceUnitSpecies *rus = const_cast<ResourceUnitSpecies*>(ru->ruSpecies()[tree.species_index]); |
182 | ResourceUnitSpecies *rus = const_cast<ResourceUnitSpecies*>(ru->ruSpecies()[tree.species_index]); |
166 | const Species *species = rus->species(); |
183 | const Species *species = rus->species(); |
167 | 184 | ||
168 | // (1) calculate height growth potential for the tree (uses linerization of expressions...)
|
185 | // (1) calculate height growth potential for the tree (uses linerization of expressions...)
|
169 | double h_pot = species->saplingGrowthParameters().heightGrowthPotential.calculate(tree.height); |
186 | double h_pot = species->saplingGrowthParameters().heightGrowthPotential.calculate(tree.height); |
170 | double delta_h_pot = h_pot - tree.height; |
187 | double delta_h_pot = h_pot - tree.height; |
171 | 188 | ||
172 | // (2) reduce height growth potential with species growth response f_env_yr and with light state (i.e. LIF-value) of home-pixel.
|
189 | // (2) reduce height growth potential with species growth response f_env_yr and with light state (i.e. LIF-value) of home-pixel.
|
173 | if (dom_height==0.f) |
190 | if (dom_height==0.f) |
174 | throw IException(QString("growSapling: height grid at %1/%2 has value 0").arg(isc)); |
191 | throw IException(QString("growSapling: height grid at %1/%2 has value 0").arg(isc)); |
175 | 192 | ||
176 | double rel_height = tree.height / dom_height; |
193 | double rel_height = tree.height / dom_height; |
177 | 194 | ||
178 | double lif_corrected = species->speciesSet()->LRIcorrection(lif_value, rel_height); // correction based on height |
195 | double lif_corrected = species->speciesSet()->LRIcorrection(lif_value, rel_height); // correction based on height |
179 | 196 | ||
180 | double lr = species->lightResponse(lif_corrected); // species specific light response (LUI, light utilization index) |
197 | double lr = species->lightResponse(lif_corrected); // species specific light response (LUI, light utilization index) |
181 | 198 | ||
182 | rus->calculate(true); // calculate the 3pg module (this is done only if that did not happen up to now); true: call comes from regeneration |
199 | rus->calculate(true); // calculate the 3pg module (this is done only if that did not happen up to now); true: call comes from regeneration |
183 | double f_env_yr = rus->prod3PG().fEnvYear(); |
200 | double f_env_yr = rus->prod3PG().fEnvYear(); |
184 | 201 | ||
185 | double delta_h_factor = f_env_yr * lr; // relative growth |
202 | double delta_h_factor = f_env_yr * lr; // relative growth |
186 | 203 | ||
187 | if (h_pot<0. || delta_h_pot<0. || lif_corrected<0. || lif_corrected>1. || delta_h_factor<0. || delta_h_factor>1. ) |
204 | if (h_pot<0. || delta_h_pot<0. || lif_corrected<0. || lif_corrected>1. || delta_h_factor<0. || delta_h_factor>1. ) |
188 | qDebug() << "invalid values in Sapling::growSapling"; |
205 | qDebug() << "invalid values in Sapling::growSapling"; |
189 | 206 | ||
190 | // check browsing
|
207 | // check browsing
|
191 | if (mBrowsingPressure>0. && tree.height<=2.f) { |
208 | if (mBrowsingPressure>0. && tree.height<=2.f) { |
192 | double p = rus->species()->saplingGrowthParameters().browsingProbability; |
209 | double p = rus->species()->saplingGrowthParameters().browsingProbability; |
193 | // calculate modifed annual browsing probability via odds-ratios
|
210 | // calculate modifed annual browsing probability via odds-ratios
|
194 | // odds = p/(1-p) -> odds_mod = odds * browsingPressure -> p_mod = odds_mod /( 1 + odds_mod) === p*pressure/(1-p+p*pressure)
|
211 | // odds = p/(1-p) -> odds_mod = odds * browsingPressure -> p_mod = odds_mod /( 1 + odds_mod) === p*pressure/(1-p+p*pressure)
|
195 | double p_browse = p*mBrowsingPressure / (1. - p + p*mBrowsingPressure); |
212 | double p_browse = p*mBrowsingPressure / (1. - p + p*mBrowsingPressure); |
196 | if (drandom() < p_browse) { |
213 | if (drandom() < p_browse) { |
197 | delta_h_factor = 0.; |
214 | delta_h_factor = 0.; |
198 | }
|
215 | }
|
199 | }
|
216 | }
|
200 | 217 | ||
201 | // check mortality of saplings
|
218 | // check mortality of saplings
|
202 | if (delta_h_factor < species->saplingGrowthParameters().stressThreshold) { |
219 | if (delta_h_factor < species->saplingGrowthParameters().stressThreshold) { |
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 | }
|
213 | DBG_IF(delta_h_pot*delta_h_factor < 0.f || delta_h_pot*delta_h_factor > 2., "Sapling::growSapling", "inplausible height growth."); |
231 | DBG_IF(delta_h_pot*delta_h_factor < 0.f || delta_h_pot*delta_h_factor > 2., "Sapling::growSapling", "inplausible height growth."); |
214 | 232 | ||
215 | // grow
|
233 | // grow
|
216 | tree.height += delta_h_pot * delta_h_factor; |
234 | tree.height += delta_h_pot * delta_h_factor; |
217 | tree.age++; // increase age of sapling by 1 |
235 | tree.age++; // increase age of sapling by 1 |
218 | 236 | ||
219 | // recruitment?
|
237 | // recruitment?
|
220 | if (tree.height > 4.f) { |
238 | if (tree.height > 4.f) { |
221 | rus->saplingStat().mRecruited++; |
239 | rus->saplingStat().mRecruited++; |
222 | 240 | ||
223 | float dbh = tree.height / species->saplingGrowthParameters().hdSapling * 100.f; |
241 | float dbh = tree.height / species->saplingGrowthParameters().hdSapling * 100.f; |
224 | // the number of trees to create (result is in trees per pixel)
|
242 | // the number of trees to create (result is in trees per pixel)
|
225 | double n_trees = species->saplingGrowthParameters().representedStemNumber(dbh); |
243 | double n_trees = species->saplingGrowthParameters().representedStemNumber(dbh); |
226 | int to_establish = static_cast<int>( n_trees ); |
244 | int to_establish = static_cast<int>( n_trees ); |
227 | 245 | ||
228 | // if n_trees is not an integer, choose randomly if we should add a tree.
|
246 | // if n_trees is not an integer, choose randomly if we should add a tree.
|
229 | // e.g.: n_trees = 2.3 -> add 2 trees with 70% probability, and add 3 trees with p=30%.
|
247 | // e.g.: n_trees = 2.3 -> add 2 trees with 70% probability, and add 3 trees with p=30%.
|
230 | if (drandom() < (n_trees-to_establish) || to_establish==0) |
248 | if (drandom() < (n_trees-to_establish) || to_establish==0) |
231 | to_establish++; |
249 | to_establish++; |
232 | 250 | ||
233 | // add a new tree
|
251 | // add a new tree
|
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)); |
257 | bigtree.setDbh(static_cast<float>(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation))); |
240 | bigtree.setHeight(tree.height * 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; |
246 | const_cast<ResourceUnitSpecies*>(rus)->statistics().add(t, 0); // count the newly created trees already in the stats |
264 | const_cast<ResourceUnitSpecies*>(rus)->statistics().add(t, 0); // count the newly created trees already in the stats |
247 | }
|
265 | }
|
248 | // clear all regeneration from this pixel (including this tree)
|
266 | // clear all regeneration from this pixel (including this tree)
|
249 | tree.clear(); // clear this tree (no carbon flow to the ground) |
267 | tree.clear(); // clear this tree (no carbon flow to the ground) |
250 | SaplingCell &s=mGrid[isc]; |
268 | SaplingCell &s=mGrid[isc]; |
251 | for (int i=0;i<NSAPCELLS;++i) { |
269 | for (int i=0;i<NSAPCELLS;++i) { |
252 | if (s.saplings[i].is_occupied()) { |
270 | if (s.saplings[i].is_occupied()) { |
253 | // add carbon to the ground
|
271 | // add carbon to the ground
|
254 | rus->saplingStat().addCarbonOfDeadSapling( s.saplings[i].height / species->saplingGrowthParameters().hdSapling * 100.f ); |
272 | rus->saplingStat().addCarbonOfDeadSapling( s.saplings[i].height / species->saplingGrowthParameters().hdSapling * 100.f ); |
255 | s.saplings[i].clear(); |
273 | s.saplings[i].clear(); |
256 | }
|
274 | }
|
257 | }
|
275 | }
|
258 | return true; // need cleanup |
276 | return true; // need cleanup |
259 | }
|
277 | }
|
260 | // book keeping (only for survivors) for the sapling of the resource unit / species
|
278 | // book keeping (only for survivors) for the sapling of the resource unit / species
|
261 | SaplingStat &ss = rus->saplingStat(); |
279 | SaplingStat &ss = rus->saplingStat(); |
262 | ss.mLiving++; |
280 | ss.mLiving++; |
263 | ss.mAvgHeight+=tree.height; |
281 | ss.mAvgHeight+=tree.height; |
264 | ss.mAvgAge+=tree.age; |
282 | ss.mAvgAge+=tree.age; |
265 | ss.mAvgDeltaHPot+=delta_h_pot; |
283 | ss.mAvgDeltaHPot+=delta_h_pot; |
266 | ss.mAvgHRealized += delta_h_pot * delta_h_factor; |
284 | ss.mAvgHRealized += delta_h_pot * delta_h_factor; |
267 | return false; |
285 | return false; |
268 | }
|
286 | }
|
269 | 287 | ||
270 | void SaplingStat::clearStatistics() |
288 | void SaplingStat::clearStatistics() |
271 | {
|
289 | {
|
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 | }
|
279 | 371 |