Subversion Repositories public iLand

Rev

Rev 697 | Rev 770 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1
 
671 werner 2
/********************************************************************************************
3
**    iLand - an individual based forest landscape and disturbance model
4
**    http://iland.boku.ac.at
5
**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl
6
**
7
**    This program is free software: you can redistribute it and/or modify
8
**    it under the terms of the GNU General Public License as published by
9
**    the Free Software Foundation, either version 3 of the License, or
10
**    (at your option) any later version.
11
**
12
**    This program is distributed in the hope that it will be useful,
13
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
14
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15
**    GNU General Public License for more details.
16
**
17
**    You should have received a copy of the GNU General Public License
18
**    along with this program.  If not, see <http://www.gnu.org/licenses/>.
19
********************************************************************************************/
20
 
468 werner 21
#include "snag.h"
22
#include "tree.h"
23
#include "species.h"
24
#include "globalsettings.h"
25
#include "expression.h"
490 werner 26
// for calculation of climate decomposition
27
#include "resourceunit.h"
28
#include "watercycle.h"
29
#include "climate.h"
541 werner 30
#include "model.h"
468 werner 31
 
32
/** @class Snag
697 werner 33
  @ingroup core
468 werner 34
  Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
490 werner 35
  Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
468 werner 36
  is subsequently forwarded to the soil sub model.
522 werner 37
  Carbon is stored in three classes (depending on the size)
528 werner 38
  The Snag dynamics class uses the following species parameter:
39
  cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW
468 werner 40
 
41
  */
42
// static variables
528 werner 43
double Snag::mDBHLower = -1.;
522 werner 44
double Snag::mDBHHigher = 0.;
45
double Snag::mCarbonThreshold[] = {0., 0., 0.};
46
 
534 werner 47
double CNPair::biomassCFraction = biomassCFraction; // get global from globalsettings.h
468 werner 48
 
534 werner 49
/// add biomass and weigh the parameter_value with the current C-content of the pool
50
void CNPool::addBiomass(const double biomass, const double CNratio, const double parameter_value)
51
{
52
    if (biomass==0.)
53
        return;
54
    double new_c = biomass*biomassCFraction;
55
    double p_old = C / (new_c + C);
56
    mParameter = mParameter*p_old + parameter_value*(1.-p_old);
57
    CNPair::addBiomass(biomass, CNratio);
58
}
59
 
60
// increase pool (and weigh the value)
61
void CNPool::operator+=(const CNPool &s)
62
{
63
    if (s.C==0.)
64
        return;
65
    mParameter = parameter(s); // calculate weighted parameter
66
    C+=s.C;
67
    N+=s.N;
68
}
69
 
70
double CNPool::parameter(const CNPool &s) const
71
{
72
    if (s.C==0.)
73
        return parameter();
74
    double p_old = C / (s.C + C);
75
    double result =  mParameter*p_old + s.parameter()*(1.-p_old);
76
    return result;
77
}
78
 
79
 
522 werner 80
void Snag::setupThresholds(const double lower, const double upper)
81
{
82
    if (mDBHLower == lower)
83
        return;
84
    mDBHLower = lower;
85
    mDBHHigher = upper;
86
    mCarbonThreshold[0] = lower / 2.;
87
    mCarbonThreshold[1] = lower + (upper - lower)/2.;
88
    mCarbonThreshold[2] = upper + (upper - lower)/2.;
89
    //# threshold levels for emptying out the dbh-snag-classes
90
    //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%
91
    //# values in kg!
92
    for (int i=0;i<3;i++)
93
        mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1;
94
}
95
 
96
 
468 werner 97
Snag::Snag()
98
{
490 werner 99
    mRU = 0;
534 werner 100
    CNPair::setCFraction(biomassCFraction);
468 werner 101
}
102
 
490 werner 103
void Snag::setup( const ResourceUnit *ru)
468 werner 104
{
490 werner 105
    mRU = ru;
106
    mClimateFactor = 0.;
468 werner 107
    // branches
108
    mBranchCounter=0;
109
    for (int i=0;i<3;i++) {
110
        mTimeSinceDeath[i] = 0.;
111
        mNumberOfSnags[i] = 0.;
522 werner 112
        mAvgDbh[i] = 0.;
113
        mAvgHeight[i] = 0.;
114
        mAvgVolume[i] = 0.;
115
        mKSW[i] = 0.;
116
        mCurrentKSW[i] = 0.;
557 werner 117
        mHalfLife[i] = 0.;
468 werner 118
    }
475 werner 119
    mTotalSnagCarbon = 0.;
528 werner 120
    if (mDBHLower<=0)
121
        throw IException("Snag::setupThresholds() not called or called with invalid parameters.");
557 werner 122
 
123
    // Inital values from XML file
124
    XmlHelper xml=GlobalSettings::instance()->settings();
574 werner 125
    double kyr = xml.valueDouble("model.site.youngRefractoryDecompRate", -1);
557 werner 126
    // put carbon of snags to the middle size class
127
    xml.setCurrentNode("model.initialization.snags");
128
    mSWD[1].C = xml.valueDouble(".swdC");
129
    mSWD[1].N = mSWD[1].C / xml.valueDouble(".swdCN", 50.);
130
    mSWD[1].setParameter(kyr);
131
    mKSW[1] = xml.valueDouble(".swdDecompRate");
132
    mNumberOfSnags[1] = xml.valueDouble(".swdCount");
133
    mHalfLife[1] = xml.valueDouble(".swdHalfLife");
134
    // and for the Branch/coarse root pools: split the init value into five chunks
135
    CNPool other(xml.valueDouble(".otherC"), xml.valueDouble(".otherC")/xml.valueDouble(".otherCN", 50.), kyr );
136
 
137
    mTotalSnagCarbon = other.C + mSWD[1].C;
138
 
139
    other *= 0.2;
140
    for (int i=0;i<5;i++)
141
        mOtherWood[i] = other;
468 werner 142
}
143
 
475 werner 144
// debug outputs
145
QList<QVariant> Snag::debugList()
146
{
147
    // list columns
148
    // for three pools
149
    QList<QVariant> list;
150
 
523 werner 151
    // totals
152
    list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N;
477 werner 153
    // fluxes to labile soil pool and to refractory soil pool
524 werner 154
    list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N;
475 werner 155
 
156
    for (int i=0;i<3;i++) {
157
        // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"
158
        list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N;
524 werner 159
        list << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i];
475 werner 160
    }
161
 
540 werner 162
    // branch/coarse wood pools (5 yrs)
163
    for (int i=0;i<5;i++) {
164
        list << mOtherWood[i].C << mOtherWood[i].N;
165
    }
166
//    list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N
167
//            << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N
168
//            << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N
169
//            << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N
170
//            << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;
475 werner 171
    return list;
172
}
173
 
713 werner 174
 
468 werner 175
void Snag::newYear()
176
{
177
    for (int i=0;i<3;i++) {
178
        mToSWD[i].clear(); // clear transfer pools to standing-woody-debris
522 werner 179
        mCurrentKSW[i] = 0.;
468 werner 180
    }
181
    mLabileFlux.clear();
182
    mRefractoryFlux.clear();
476 werner 183
    mTotalToAtm.clear();
184
    mTotalToExtern.clear();
609 werner 185
    mTotalToDisturbance.clear();
476 werner 186
    mTotalIn.clear();
477 werner 187
    mSWDtoSoil.clear();
468 werner 188
}
189
 
490 werner 190
/// calculate the dynamic climate modifier for decomposition 're'
522 werner 191
/// calculation is done on the level of ResourceUnit because the water content per day is needed.
490 werner 192
double Snag::calculateClimateFactors()
193
{
194
    double ft, fw;
195
    double f_sum = 0.;
552 werner 196
    int iday=0;
553 werner 197
    // calculate the water-factor for each month (see Adair et al 2008)
198
    double fw_month[12];
199
    double ratio;
200
    for (int m=0;m<12;m++) {
562 werner 201
        if (mRU->waterCycle()->referenceEvapotranspiration()[m]>0.)
202
            ratio = mRU->climate()->precipitationMonth()[m] /  mRU->waterCycle()->referenceEvapotranspiration()[m];
553 werner 203
        else
204
            ratio = 0;
205
        fw_month[m] = 1. / (1. + 30.*exp(-8.5*ratio));
564 werner 206
        if (logLevelDebug()) qDebug() <<"month"<< m << "PET" << mRU->waterCycle()->referenceEvapotranspiration()[m] << "prec" <<mRU->climate()->precipitationMonth()[m];
553 werner 207
    }
208
 
552 werner 209
    for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday)
490 werner 210
    {
211
        ft = exp(308.56*(1./56.02-1./((273.+day->temperature)-227.13)));  // empirical variable Q10 model of Lloyd and Taylor (1994), see also Adair et al. (2008)
553 werner 212
        fw = fw_month[day->month-1];
540 werner 213
 
490 werner 214
        f_sum += ft*fw;
215
    }
216
    // the climate factor is defined as the arithmentic annual mean value
217
    mClimateFactor = f_sum / double(mRU->climate()->daysOfYear());
218
    return mClimateFactor;
219
}
220
 
522 werner 221
/// do the yearly calculation
222
/// see http://iland.boku.ac.at/snag+dynamics
526 werner 223
void Snag::calculateYear()
468 werner 224
{
522 werner 225
    mSWDtoSoil.clear();
532 werner 226
    const double climate_factor_re = calculateClimateFactors(); // calculate anyway, because also the soil module needs it (and currently one can have Snag and Soil only as a couple)
477 werner 227
    if (isEmpty()) // nothing to do
475 werner 228
        return;
229
 
468 werner 230
    // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
540 werner 231
    mRefractoryFlux+=mOtherWood[mBranchCounter];
232
 
233
    mOtherWood[mBranchCounter].clear();
468 werner 234
    mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0.
540 werner 235
    // decay of branches/coarse roots
236
    for (int i=0;i<5;i++) {
237
        if (mOtherWood[i].C>0.) {
238
            double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value...
239
            mOtherWood[i].C *= survive_rate;
240
        }
241
    }
468 werner 242
 
243
    // process standing snags.
244
    // the input of the current year is in the mToSWD-Pools
245
    for (int i=0;i<3;i++) {
246
 
522 werner 247
        // update the swd-pool with this years' input
248
        if (!mToSWD[i].isEmpty()) {
249
            // update decay rate (apply average yearly input to the state parameters)
250
            mKSW[i] = mKSW[i]*(mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i]*(mToSWD[i].C/(mSWD[i].C+mToSWD[i].C));
251
            //move content to the SWD pool
252
            mSWD[i] += mToSWD[i];
253
        }
475 werner 254
 
522 werner 255
        if (mSWD[i].C > 0) {
256
            // reduce the Carbon (note: the N stays, thus the CN ratio changes)
257
            // use the decay rate that is derived as a weighted average of all standing woody debris
523 werner 258
            double survive_rate = exp(-mKSW[i] *climate_factor_re * 1. ); // 1: timestep
259
            mTotalToAtm.C += mSWD[i].C * (1. - survive_rate);
260
            mSWD[i].C *= survive_rate;
468 werner 261
 
522 werner 262
            // transition to downed woody debris
263
            // update: use negative exponential decay, species parameter: half-life
264
            // modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa
265
            // this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),
266
            // where lag times for cohort transitions are linearly modified with re although here individual good or bad years have
267
            // an immediate effect, the average climatic influence should come through (and it is inherently transient)
268
            // note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)
269
            // needs to be calculated, followed by a weighted update of the previous swd.hl.
270
            // As weights here we use stem number, as the processes here pertain individual snags
271
            // calculate the transition probability of SWD to downed dead wood
468 werner 272
 
522 werner 273
            double half_life = mHalfLife[i] / climate_factor_re;
274
            double rate = -M_LN2 / half_life; // M_LN2: math. constant
275
 
276
            // higher decay rate for the class with smallest diameters
277
            if (i==0)
278
                rate*=2.;
279
 
523 werner 280
            double transfer = 1. - exp(rate);
522 werner 281
 
468 werner 282
            // calculate flow to soil pool...
522 werner 283
            mSWDtoSoil += mSWD[i] * transfer;
284
            mRefractoryFlux += mSWD[i] * transfer;
285
            mSWD[i] *= (1.-transfer); // reduce pool
468 werner 286
            // calculate the stem number of remaining snags
522 werner 287
            mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer);
523 werner 288
 
289
            mTimeSinceDeath[i] += 1.;
522 werner 290
            // if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats
291
            // also, if the Carbon of an average snag is less than 10% of the original average tree
292
            // (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD
293
            if (mNumberOfSnags[i] < 0.5 || mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) {
294
                // clear the pool: add the rest to the soil, clear statistics of the pool
468 werner 295
                mRefractoryFlux += mSWD[i];
522 werner 296
                mSWDtoSoil += mSWD[i];
468 werner 297
                mSWD[i].clear();
522 werner 298
                mAvgDbh[i] = 0.;
299
                mAvgHeight[i] = 0.;
300
                mAvgVolume[i] = 0.;
301
                mKSW[i] = 0.;
302
                mCurrentKSW[i] = 0.;
303
                mHalfLife[i] = 0.;
304
                mTimeSinceDeath[i] = 0.;
468 werner 305
            }
522 werner 306
 
468 werner 307
        }
522 werner 308
 
468 werner 309
    }
522 werner 310
    // total carbon in the snag-container on the RU *after* processing is the content of the
475 werner 311
    // standing woody debris pools + the branches
312
    mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C +
540 werner 313
                       mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C;
587 werner 314
    mTotalSWD = mSWD[0] + mSWD[1] + mSWD[2];
315
    mTotalOther = mOtherWood[0] + mOtherWood[1] + mOtherWood[2] + mOtherWood[3] + mOtherWood[4];
468 werner 316
}
317
 
318
/// foliage and fineroot litter is transferred during tree growth.
588 werner 319
void Snag::addTurnoverLitter(const Species *species, const double litter_foliage, const double litter_fineroot)
468 werner 320
{
588 werner 321
    mLabileFlux.addBiomass(litter_foliage, species->cnFoliage(), species->snagKyl());
322
    mLabileFlux.addBiomass(litter_fineroot, species->cnFineroot(), species->snagKyl());
468 werner 323
}
324
 
595 werner 325
void Snag::addTurnoverWood(const Species *species, const double woody_biomass)
326
{
327
    mRefractoryFlux.addBiomass(woody_biomass, species->cnWood(), species->snagKyr());
328
}
329
 
713 werner 330
 
331
/** process the remnants of a single tree.
332
    The part of the stem / branch not covered by snag/soil fraction is removed from the system (e.g. harvest, fire)
333
  @param tree the tree to process
334
  @param stem_to_snag fraction (0..1) of the stem biomass that should be moved to a standing snag
335
  @param stem_to_soil fraction (0..1) of the stem biomass that should go directly to the soil
336
  @param branch_to_snag fraction (0..1) of the branch biomass that should be moved to a standing snag
337
  @param branch_to_soil fraction (0..1) of the branch biomass that should go directly to the soil
338
  @param foliage_to_soil fraction (0..1) of the foliage biomass that should go directly to the soil
339
 
340
*/
341
void Snag::addBiomassPools(const Tree *tree,
342
                           const double stem_to_snag, const double stem_to_soil,
343
                           const double branch_to_snag, const double branch_to_soil,
344
                           const double foliage_to_soil)
468 werner 345
{
528 werner 346
    const Species *species = tree->species();
468 werner 347
 
713 werner 348
    double branch_biomass = tree->biomassBranch();
349
    // fine roots go to the labile pool
350
    mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), species->snagKyl());
468 werner 351
 
713 werner 352
    // a part of the foliage goes to the soil
353
    mLabileFlux.addBiomass(tree->biomassFoliage() * foliage_to_soil, species->cnFoliage(), species->snagKyl());
354
 
355
    //coarse roots and a part of branches are equally distributed over five years:
356
    double biomass_rest = (tree->biomassCoarseRoot() + branch_to_snag*branch_biomass) * 0.2;
468 werner 357
    for (int i=0;i<5; i++)
713 werner 358
        mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), species->snagKyr());
468 werner 359
 
713 werner 360
    // the other part of the branches goes directly to the soil
361
    mRefractoryFlux.addBiomass(branch_biomass*branch_to_soil, species->cnWood(), species->snagKyr() );
362
    // a part of the stem wood goes directly to the soil
363
    mRefractoryFlux.addBiomass(tree->biomassStem()*stem_to_soil, species->cnWood(), species->snagKyr() );
364
 
365
    // just for book-keeping: keep track of all inputs of branches / roots / swd into the "snag" pools
366
    mTotalIn.addBiomass(tree->biomassBranch()*branch_to_snag + tree->biomassCoarseRoot() + tree->biomassStem()*stem_to_snag, species->cnWood());
468 werner 367
    // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
522 werner 368
    int pi = poolIndex(tree->dbh()); // get right transfer pool
369
 
713 werner 370
    if (stem_to_snag>0.) {
371
        // update statistics - stemnumber-weighted averages
372
        // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
373
        double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers)
374
        double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1).
375
        mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new;
376
        mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new;
377
        mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new;
378
        mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new;
379
        mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new;
522 werner 380
 
713 werner 381
        // average the decay rate (ksw); this is done based on the carbon content
382
        // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
383
        p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
384
        p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
385
        mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new;
386
        mNumberOfSnags[pi]++;
387
    }
523 werner 388
 
713 werner 389
    // finally add the biomass of the stem to the standing snag pool
534 werner 390
    CNPool &to_swd = mToSWD[pi];
713 werner 391
    to_swd.addBiomass(tree->biomassStem()*stem_to_snag, species->cnWood(), species->snagKyr());
392
 
393
    // the biomass that is not routed to snags or directly to the soil
394
    // is removed from the system (to atmosphere or harvested)
395
    mTotalToExtern.addBiomass(tree->biomassFoliage()* (1. - foliage_to_soil) +
396
                              branch_biomass*(1. - branch_to_snag - branch_to_soil) +
397
                              tree->biomassStem()*(1. - stem_to_snag - stem_to_soil), species->cnWood());
398
 
468 werner 399
}
400
 
713 werner 401
 
402
/// after the death of the tree the five biomass compartments are processed.
403
void Snag::addMortality(const Tree *tree)
404
{
405
    addBiomassPools(tree, 1., 0.,  // all stem biomass goes to snag
406
                    1., 0.,        // all branch biomass to snag
407
                    1.);           // all foliage to soil
408
 
409
//    const Species *species = tree->species();
410
 
411
//    // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool
412
//    mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl());
413
//    mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
414
 
415
//    // branches and coarse roots are equally distributed over five years:
416
//    double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2;
417
//    for (int i=0;i<5; i++)
418
//        mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
419
 
420
//    // just for book-keeping: keep track of all inputs into branches / roots / swd
421
//    mTotalIn.addBiomass(tree->biomassBranch() + tree->biomassCoarseRoot() + tree->biomassStem(), species->cnWood());
422
//    // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
423
//    int pi = poolIndex(tree->dbh()); // get right transfer pool
424
 
425
//    // update statistics - stemnumber-weighted averages
426
//    // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
427
//    double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers)
428
//    double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1).
429
//    mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new;
430
//    mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new;
431
//    mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new;
432
//    mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new;
433
//    mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new;
434
 
435
//    // average the decay rate (ksw); this is done based on the carbon content
436
//    // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
437
//    if (tree->biomassStem()==0)
438
//        throw IException("Snag::addMortality: tree without stem biomass!!");
439
//    p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
440
//    p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
441
//    mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new;
442
//    mNumberOfSnags[pi]++;
443
 
444
//    // finally add the biomass
445
//    CNPool &to_swd = mToSWD[pi];
446
//    to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr());
447
}
448
 
468 werner 449
/// add residual biomass of 'tree' after harvesting.
522 werner 450
/// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is *removed* by the harvest operation (i.e.: not to stay in the system)
528 werner 451
/// records on harvested biomass is collected (mTotalToExtern-pool).
468 werner 452
void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction )
453
{
713 werner 454
    addBiomassPools(tree,
455
                    0., 1.-remove_stem_fraction, // "remove_stem_fraction" is removed -> the rest goes to soil
456
                    0., 1.-remove_branch_fraction, // "remove_branch_fraction" is removed -> the rest goes directly to the soil
457
                    1.-remove_foliage_fraction); // the rest of foliage is routed to the soil
458
//    const Species *species = tree->species();
468 werner 459
 
713 werner 460
//    // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool
461
//    mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl());
462
//    mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
540 werner 463
 
713 werner 464
//    // for branches, add all biomass that remains in the forest to the soil
465
//    mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr());
466
//    // the same treatment for stem residuals
467
//    mRefractoryFlux.addBiomass(tree->biomassStem() * (1. - remove_stem_fraction), species->cnWood(), tree->species()->snagKyr());
468 werner 468
 
713 werner 469
//    // split the corase wood biomass into parts (slower decay)
470
//    double biomass_rest = (tree->biomassCoarseRoot()) * 0.2;
471
//    for (int i=0;i<5; i++)
472
//        mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
540 werner 473
 
474
 
713 werner 475
//    // for book-keeping...
476
//    mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction +
477
//                              tree->biomassBranch()*remove_branch_fraction +
478
//                              tree->biomassStem()*remove_stem_fraction, species->cnWood());
468 werner 479
}
480
 
588 werner 481
// add flow from regeneration layer (dead trees) to soil
595 werner 482
void Snag::addToSoil(const Species *species, const CNPair &woody_pool, const CNPair &litter_pool)
588 werner 483
{
484
    mLabileFlux.add(litter_pool, species->snagKyl());
485
    mRefractoryFlux.add(woody_pool, species->snagKyr());
486
}
534 werner 487
 
607 werner 488
/// disturbance function: remove the fraction of 'factor' of biomass from the SWD pools; 0: remove nothing, 1: remove all
489
/// biomass removed by this function goes to the atmosphere
490
void Snag::removeCarbon(const double factor)
491
{
492
    // reduce pools of currently standing dead wood and also of pools that are added
493
    // during (previous) management operations of the current year
494
    for (int i=0;i<3;i++) {
609 werner 495
        mTotalToDisturbance += (mSWD[i] + mToSWD[i]) * factor;
607 werner 496
        mSWD[i] *= (1. - factor);
497
        mToSWD[i] *= (1. - factor);
498
    }
534 werner 499
 
607 werner 500
    for (int i=0;i<5;i++) {
609 werner 501
        mTotalToDisturbance += mOtherWood[i]*factor;
607 werner 502
        mOtherWood[i] *= (1. - factor);
503
    }
504
}
505
 
506
 
507
/// cut down swd (and branches) and move to soil pools
508
/// @param factor 0: cut 0%, 1: cut and slash 100% of the wood
509
void Snag::management(const double factor)
510
{
511
    if (factor<0. || factor>1.)
512
        throw IException(QString("Invalid factor in Snag::management: '%1'").arg(factor));
513
    // swd pools
514
    for (int i=0;i<3;i++) {
515
        mSWDtoSoil += mSWD[i] * factor;
516
        mSWD[i] *= (1. - factor);
517
        mSWDtoSoil += mToSWD[i] * factor;
518
        mToSWD[i] *= (1. - factor);
519
    }
520
    // what to do with the branches: now move also all wood to soil (note: this is note
521
    // very good w.r.t the coarse roots...
522
    for (int i=0;i<5;i++) {
523
        mRefractoryFlux+=mOtherWood[i]*factor;
524
        mOtherWood[i]*=(1. - factor);
525
    }
526
 
527
}
528
 
529