Subversion Repositories public iLand

Rev

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, ++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;
90
                    SaplingTree *stree=0;
-
 
91
                    SaplingTree *slot=s->saplings;
85
                    for (int i=0;i<NSAPCELLS;++i) {
92
                    for (int i=0;i<NSAPCELLS;++i, ++slot) {
86
                        if (!s->saplings[i].is_occupied() && i_occupied<0)
93
                        if (!stree && !slot->is_occupied())
87
                            i_occupied=i;
94
                            stree=slot;
88
                        if (s->saplings[i].species_index == species_idx) {
95
                        if (slot->species_index == species_idx) {
-
 
96
                            stree=0;
89
                            viable = false;
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];
-
 
108
-
 
109
                        double &lif_corrected = lif_corr[iy*cPxPerRU+ix];
-
 
110
                        // calculate the LIFcorrected only once per pixel
-
 
111
                        if (lif_corrected<0.)
101
                        double lif_corrected = rus->species()->speciesSet()->LRIcorrection(lif_value, 4. / hgv.height);
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
                             // ok, lets add a sapling at the given position
116
                            // ok, lets add a sapling at the given position (age is incremented later)
105
                             s->saplings[i_occupied].setSapling(0.05f, 1, species_idx);
117
                            stree->setSapling(0.05f, 0, species_idx);
106
                             s->checkState();
118
                            s->checkState();
107
                             mAdded++;
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