Subversion Repositories public iLand

Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/snag.cpp':

Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/snag.cpp':

#include "snag.h"

#include "snag.h"

#include "tree.h"

#include "tree.h"

#include "species.h"

#include "species.h"

#include "globalsettings.h"

#include "globalsettings.h"

#include "expression.h"

#include "expression.h"

// for calculation of climate decomposition

// for calculation of climate decomposition

#include "resourceunit.h"

#include "resourceunit.h"

#include "watercycle.h"

#include "watercycle.h"

#include "climate.h"

#include "climate.h"

/** @class Snag

/** @class Snag

  Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.

  Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.

  Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags

  Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags

  is subsequently forwarded to the soil sub model.

  is subsequently forwarded to the soil sub model.

  Carbon is stored in three classes (depending on the size)

  Carbon is stored in three classes (depending on the size)

  */

  */

// static variables

// static variables

double Snag::mDBHLower = 0.;

double Snag::mDBHLower = 0.;

double Snag::mDBHHigher = 0.;

double Snag::mDBHHigher = 0.;

double Snag::mCarbonThreshold[] = {0., 0., 0.};

double Snag::mCarbonThreshold[] = {0., 0., 0.};

double CNPool::biomassCFraction = biomassCFraction; // get global from globalsettings.h

double CNPool::biomassCFraction = biomassCFraction; // get global from globalsettings.h

void Snag::setupThresholds(const double lower, const double upper)

void Snag::setupThresholds(const double lower, const double upper)

{

{

    if (mDBHLower == lower)

    if (mDBHLower == lower)

        return;

        return;

    mDBHLower = lower;

    mDBHLower = lower;

    mDBHHigher = upper;

    mDBHHigher = upper;

    mCarbonThreshold[0] = lower / 2.;

    mCarbonThreshold[0] = lower / 2.;

    mCarbonThreshold[1] = lower + (upper - lower)/2.;

    mCarbonThreshold[1] = lower + (upper - lower)/2.;

    mCarbonThreshold[2] = upper + (upper - lower)/2.;

    mCarbonThreshold[2] = upper + (upper - lower)/2.;

    //# threshold levels for emptying out the dbh-snag-classes

    //# threshold levels for emptying out the dbh-snag-classes

    //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%

    //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%

    //# values in kg!

    //# values in kg!

    for (int i=0;i<3;i++)

    for (int i=0;i<3;i++)

        mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1;

        mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1;

}

}

// species specific soil paramters

// species specific soil paramters

struct SoilParameters

struct SoilParameters

{

{

    SoilParameters(): kyl(0.15), kyr(0.0807), ksw(0.015), cnFoliage(75.), cnFineroot(40.), cnWood(300.), halfLife(15.) {}

    SoilParameters(): kyl(0.15), kyr(0.0807), ksw(0.015), cnFoliage(75.), cnFineroot(40.), cnWood(300.), halfLife(15.) {}

    double kyl; // litter decomposition rate

    double kyl; // litter decomposition rate

    double kyr; // downed woody debris (dwd) decomposition rate

    double kyr; // downed woody debris (dwd) decomposition rate

    double ksw; // standing woody debris (swd) decomposition rate

    double ksw; // standing woody debris (swd) decomposition rate

    double cnFoliage; //  C/N foliage litter

    double cnFoliage; //  C/N foliage litter

    double cnFineroot; // C/N ratio fine root

    double cnFineroot; // C/N ratio fine root

    double cnWood; // C/N Wood: used for brances, stem and coarse root

    double cnWood; // C/N Wood: used for brances, stem and coarse root

    double halfLife; // half-life period of the given species (for calculation of SWD->DWD transition)

    double halfLife; // half-life period of the given species (for calculation of SWD->DWD transition)

} soilparams;

} soilparams;

Snag::Snag()

Snag::Snag()

{

{

    mRU = 0;

    mRU = 0;

    CNPool::setCFraction(biomassCFraction);

    CNPool::setCFraction(biomassCFraction);

}

}

void Snag::setup( const ResourceUnit *ru)

void Snag::setup( const ResourceUnit *ru)

{

{

    mRU = ru;

    mRU = ru;

    mClimateFactor = 0.;

    mClimateFactor = 0.;

    // branches

    // branches

    mBranchCounter=0;

    mBranchCounter=0;

    for (int i=0;i<3;i++) {

    for (int i=0;i<3;i++) {

        mTimeSinceDeath[i] = 0.;

        mTimeSinceDeath[i] = 0.;

        mNumberOfSnags[i] = 0.;

        mNumberOfSnags[i] = 0.;

        mAvgDbh[i] = 0.;

        mAvgDbh[i] = 0.;

        mAvgHeight[i] = 0.;

        mAvgHeight[i] = 0.;

        mAvgVolume[i] = 0.;

        mAvgVolume[i] = 0.;

        mKSW[i] = 0.;

        mKSW[i] = 0.;

        mCurrentKSW[i] = 0.;

        mCurrentKSW[i] = 0.;

    }

    }

    mTotalSnagCarbon = 0.;

    mTotalSnagCarbon = 0.;

}

}

// debug outputs

// debug outputs

QList<QVariant> Snag::debugList()

QList<QVariant> Snag::debugList()

{

{

    // list columns

    // list columns

    // for three pools

    // for three pools

    QList<QVariant> list;

    QList<QVariant> list;

    // fluxes to labile soil pool and to refractory soil pool

    // fluxes to labile soil pool and to refractory soil pool

    list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N << mSWDtoSoil.C << mSWDtoSoil.N;

    list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N << mSWDtoSoil.C << mSWDtoSoil.N;

    for (int i=0;i<3;i++) {

    for (int i=0;i<3;i++) {

        // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"

        // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"

        list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N;

        list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N;

    }

    }

    // branch pools (5 yrs)

    // branch pools (5 yrs)

    list << mBranches[mBranchCounter].C << mBranches[mBranchCounter].N

    list << mBranches[mBranchCounter].C << mBranches[mBranchCounter].N

            << mBranches[(mBranchCounter+1)%5].C << mBranches[(mBranchCounter+1)%5].N

            << mBranches[(mBranchCounter+1)%5].C << mBranches[(mBranchCounter+1)%5].N

            << mBranches[(mBranchCounter+2)%5].C << mBranches[(mBranchCounter+2)%5].N

            << mBranches[(mBranchCounter+2)%5].C << mBranches[(mBranchCounter+2)%5].N

            << mBranches[(mBranchCounter+3)%5].C << mBranches[(mBranchCounter+3)%5].N

            << mBranches[(mBranchCounter+3)%5].C << mBranches[(mBranchCounter+3)%5].N

            << mBranches[(mBranchCounter+4)%5].C << mBranches[(mBranchCounter+4)%5].N;

            << mBranches[(mBranchCounter+4)%5].C << mBranches[(mBranchCounter+4)%5].N;

    return list;

    return list;

}

}

void Snag::newYear()

void Snag::newYear()

{

{

    for (int i=0;i<3;i++) {

    for (int i=0;i<3;i++) {

        mToSWD[i].clear(); // clear transfer pools to standing-woody-debris

        mToSWD[i].clear(); // clear transfer pools to standing-woody-debris

        mCurrentKSW[i] = 0.;

        mCurrentKSW[i] = 0.;

    }

    }

    mLabileFlux.clear();

    mLabileFlux.clear();

    mRefractoryFlux.clear();

    mRefractoryFlux.clear();

    mTotalToAtm.clear();

    mTotalToAtm.clear();

    mTotalToExtern.clear();

    mTotalToExtern.clear();

    mTotalIn.clear();

    mTotalIn.clear();

    mSWDtoSoil.clear();

    mSWDtoSoil.clear();

}

}

/// calculate the dynamic climate modifier for decomposition 're'

/// calculate the dynamic climate modifier for decomposition 're'

/// calculation is done on the level of ResourceUnit because the water content per day is needed.

/// calculation is done on the level of ResourceUnit because the water content per day is needed.

double Snag::calculateClimateFactors()

double Snag::calculateClimateFactors()

{

{

    double rel_wc;

    double rel_wc;

    double ft, fw;

    double ft, fw;

    double f_sum = 0.;

    double f_sum = 0.;

    for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day)

    for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day)

    {

    {

        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)

        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)

        fw = pow(1.-exp(-0.2*rel_wc),5.); //  #  see Standcarb for the 'stable soil' pool

        fw = pow(1.-exp(-0.2*rel_wc),5.); //  #  see Standcarb for the 'stable soil' pool

        f_sum += ft*fw;

        f_sum += ft*fw;

    }

    }

    // the climate factor is defined as the arithmentic annual mean value

    // the climate factor is defined as the arithmentic annual mean value

    mClimateFactor = f_sum / double(mRU->climate()->daysOfYear());

    mClimateFactor = f_sum / double(mRU->climate()->daysOfYear());

    return mClimateFactor;

    return mClimateFactor;

}

}

/// do the yearly calculation

/// do the yearly calculation

/// see http://iland.boku.ac.at/snag+dynamics

/// see http://iland.boku.ac.at/snag+dynamics

void Snag::processYear()

void Snag::processYear()

{

{

    mSWDtoSoil.clear();

    mSWDtoSoil.clear();

    if (isEmpty()) // nothing to do

    if (isEmpty()) // nothing to do

        return;

        return;

    // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool

    // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool

    mRefractoryFlux+=mBranches[mBranchCounter];

    mRefractoryFlux+=mBranches[mBranchCounter];

    mBranches[mBranchCounter].clear();

    mBranches[mBranchCounter].clear();

    mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0.

    mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0.

    // process standing snags.

    // process standing snags.

    // the input of the current year is in the mToSWD-Pools

    // the input of the current year is in the mToSWD-Pools

    const double climate_factor_re = calculateClimateFactors();

    const double climate_factor_re = calculateClimateFactors();

    for (int i=0;i<3;i++) {

    for (int i=0;i<3;i++) {

        // update the swd-pool with this years' input

        // update the swd-pool with this years' input

        if (!mToSWD[i].isEmpty()) {

        if (!mToSWD[i].isEmpty()) {

            // update decay rate (apply average yearly input to the state parameters)

            // update decay rate (apply average yearly input to the state parameters)

            mKSW[i] = mKSW[i]*(mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i]*(mToSWD[i].C/(mSWD[i].C+mToSWD[i].C));

            mKSW[i] = mKSW[i]*(mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i]*(mToSWD[i].C/(mSWD[i].C+mToSWD[i].C));

            //move content to the SWD pool

            //move content to the SWD pool

            mSWD[i] += mToSWD[i];

            mSWD[i] += mToSWD[i];

        }

        }

        if (mSWD[i].C > 0) {

        if (mSWD[i].C > 0) {

            // reduce the Carbon (note: the N stays, thus the CN ratio changes)

            // reduce the Carbon (note: the N stays, thus the CN ratio changes)

            // use the decay rate that is derived as a weighted average of all standing woody debris

            // use the decay rate that is derived as a weighted average of all standing woody debris

            // transition to downed woody debris

            // transition to downed woody debris

            // update: use negative exponential decay, species parameter: half-life

            // update: use negative exponential decay, species parameter: half-life

            // modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa

            // modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa

            // this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),

            // this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),

            // where lag times for cohort transitions are linearly modified with re although here individual good or bad years have

            // where lag times for cohort transitions are linearly modified with re although here individual good or bad years have

            // an immediate effect, the average climatic influence should come through (and it is inherently transient)

            // an immediate effect, the average climatic influence should come through (and it is inherently transient)

            // note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)

            // note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)

            // needs to be calculated, followed by a weighted update of the previous swd.hl.

            // needs to be calculated, followed by a weighted update of the previous swd.hl.

            // As weights here we use stem number, as the processes here pertain individual snags

            // As weights here we use stem number, as the processes here pertain individual snags

            // calculate the transition probability of SWD to downed dead wood

            // calculate the transition probability of SWD to downed dead wood

            double half_life = mHalfLife[i] / climate_factor_re;

            double half_life = mHalfLife[i] / climate_factor_re;

            double rate = -M_LN2 / half_life; // M_LN2: math. constant

            double rate = -M_LN2 / half_life; // M_LN2: math. constant

            // higher decay rate for the class with smallest diameters

            // higher decay rate for the class with smallest diameters

            if (i==0)

            if (i==0)

                rate*=2.;

                rate*=2.;

            // calculate flow to soil pool...

            // calculate flow to soil pool...

            mSWDtoSoil += mSWD[i] * transfer;

            mSWDtoSoil += mSWD[i] * transfer;

            mRefractoryFlux += mSWD[i] * transfer;

            mRefractoryFlux += mSWD[i] * transfer;

            mSWD[i] *= (1.-transfer); // reduce pool

            mSWD[i] *= (1.-transfer); // reduce pool

            // calculate the stem number of remaining snags

            // calculate the stem number of remaining snags

            mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer);

            mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer);

            // if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats

            // if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats

            // also, if the Carbon of an average snag is less than 10% of the original average tree

            // also, if the Carbon of an average snag is less than 10% of the original average tree

            // (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD

            // (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD

            if (mNumberOfSnags[i] < 0.5 || mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) {

            if (mNumberOfSnags[i] < 0.5 || mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) {

                // clear the pool: add the rest to the soil, clear statistics of the pool

                // clear the pool: add the rest to the soil, clear statistics of the pool

                mRefractoryFlux += mSWD[i];

                mRefractoryFlux += mSWD[i];

                mSWDtoSoil += mSWD[i];

                mSWDtoSoil += mSWD[i];

                mSWD[i].clear();

                mSWD[i].clear();

                mAvgDbh[i] = 0.;

                mAvgDbh[i] = 0.;

                mAvgHeight[i] = 0.;

                mAvgHeight[i] = 0.;

                mAvgVolume[i] = 0.;

                mAvgVolume[i] = 0.;

                mKSW[i] = 0.;

                mKSW[i] = 0.;

                mCurrentKSW[i] = 0.;

                mCurrentKSW[i] = 0.;

                mHalfLife[i] = 0.;

                mHalfLife[i] = 0.;

                mTimeSinceDeath[i] = 0.;

                mTimeSinceDeath[i] = 0.;

            }

            }

        }

        }

    }

    }

    // total carbon in the snag-container on the RU *after* processing is the content of the

    // total carbon in the snag-container on the RU *after* processing is the content of the

    // standing woody debris pools + the branches

    // standing woody debris pools + the branches

    mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C +

    mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C +

                       mBranches[0].C + mBranches[1].C + mBranches[2].C + mBranches[3].C + mBranches[4].C;

                       mBranches[0].C + mBranches[1].C + mBranches[2].C + mBranches[3].C + mBranches[4].C;

}

}

/// foliage and fineroot litter is transferred during tree growth.

/// foliage and fineroot litter is transferred during tree growth.

void Snag::addTurnoverLitter(const Tree *tree, const double litter_foliage, const double litter_fineroot)

void Snag::addTurnoverLitter(const Tree *tree, const double litter_foliage, const double litter_fineroot)

{

{

    const SoilParameters &soil_params = soilparams; // to be updated

    const SoilParameters &soil_params = soilparams; // to be updated

    mLabileFlux.addBiomass(litter_foliage, soil_params.cnFoliage);

    mLabileFlux.addBiomass(litter_foliage, soil_params.cnFoliage);

    mLabileFlux.addBiomass(litter_fineroot, soil_params.cnFineroot);

    mLabileFlux.addBiomass(litter_fineroot, soil_params.cnFineroot);

}

}

/// after the death of the tree the five biomass compartments are processed.

/// after the death of the tree the five biomass compartments are processed.

void Snag::addMortality(const Tree *tree)

void Snag::addMortality(const Tree *tree)

{

{

    const SoilParameters &soil_params = soilparams; // to be updated

    const SoilParameters &soil_params = soilparams; // to be updated

    // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool

    // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool

    // 100% of coarse root biomass: that goes to the refractory pool

    // 100% of coarse root biomass: that goes to the refractory pool

    mLabileFlux.addBiomass(tree->biomassFoliage(), soil_params.cnFoliage);

    mLabileFlux.addBiomass(tree->biomassFoliage(), soil_params.cnFoliage);

    mLabileFlux.addBiomass(tree->biomassFineRoot(), soil_params.cnFineroot);

    mLabileFlux.addBiomass(tree->biomassFineRoot(), soil_params.cnFineroot);

    mRefractoryFlux.addBiomass(tree->biomassCoarseRoot(), soil_params.cnWood);

    mRefractoryFlux.addBiomass(tree->biomassCoarseRoot(), soil_params.cnWood);

    // branches are equally distributed over five years:

    // branches are equally distributed over five years:

    double biomass_branch = tree->biomassBranch() * 0.2;

    double biomass_branch = tree->biomassBranch() * 0.2;

    for (int i=0;i<5; i++)

    for (int i=0;i<5; i++)

        mBranches[i].addBiomass(biomass_branch, soil_params.cnWood);

        mBranches[i].addBiomass(biomass_branch, soil_params.cnWood);

    // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool

    // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool

    int pi = poolIndex(tree->dbh()); // get right transfer pool

    int pi = poolIndex(tree->dbh()); // get right transfer pool

    // update statistics - stemnumber-weighted averages

    // update statistics - stemnumber-weighted averages

    // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)

    // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)

    double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers)

    double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers)

    double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1).

    double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1).

    mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new;

    mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new;

    mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new;

    mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new;

    mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new;

    mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new;

    mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new;

    mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new;

    mHalfLife[pi] = mHalfLife[pi]*p_old + soil_params.halfLife * p_new;

    mHalfLife[pi] = mHalfLife[pi]*p_old + soil_params.halfLife * p_new;

    // average the decay rate (ksw); this is done based on the carbon content

    // average the decay rate (ksw); this is done based on the carbon content

    // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)

    // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)

    if (tree->biomassStem()==0)

    if (tree->biomassStem()==0)

        throw IException("Snag::addMortality: tree without stem biomass!!");

        throw IException("Snag::addMortality: tree without stem biomass!!");

    p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);

    p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);

    p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);

    p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);

    mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + soil_params.ksw * p_new;

    mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + soil_params.ksw * p_new;

    mNumberOfSnags[pi]++;

    mNumberOfSnags[pi]++;

}

}

/// add residual biomass of 'tree' after harvesting.

/// add residual biomass of 'tree' after harvesting.

/// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is *removed* by the harvest operation (i.e.: not to stay in the system)

/// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is *removed* by the harvest operation (i.e.: not to stay in the system)

/// the harvested biomass is collected.

/// the harvested biomass is collected.

void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction )

void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction )

{

{

    const SoilParameters &soil_params = soilparams; // to be updated

    const SoilParameters &soil_params = soilparams; // to be updated

    // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool

    // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool

    // 100% of coarse root biomass: that goes to the refractory pool

    // 100% of coarse root biomass: that goes to the refractory pool

    mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), soil_params.cnFoliage);

    mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), soil_params.cnFoliage);

    mLabileFlux.addBiomass(tree->biomassFineRoot(), soil_params.cnFineroot);

    mLabileFlux.addBiomass(tree->biomassFineRoot(), soil_params.cnFineroot);

    mRefractoryFlux.addBiomass(tree->biomassCoarseRoot(), soil_params.cnWood);

    mRefractoryFlux.addBiomass(tree->biomassCoarseRoot(), soil_params.cnWood);

    // residual branches are equally distributed over five years:

    // residual branches are equally distributed over five years:

    for (int i=0;i<5; i++)

    for (int i=0;i<5; i++)

    // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool

    // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool

    // TODO: what to do with harvest and stems??? I think harvested stems (that are not removed) should

    // TODO: what to do with harvest and stems??? I think harvested stems (that are not removed) should

    // go directly to the DWD??

    // go directly to the DWD??

    CNPool &swd = mToSWD[poolIndex(tree->dbh())]; // get right transfer pool

    CNPool &swd = mToSWD[poolIndex(tree->dbh())]; // get right transfer pool

    swd.addBiomass(tree->biomassStem() * remove_stem_fraction, soil_params.cnWood);

    swd.addBiomass(tree->biomassStem() * remove_stem_fraction, soil_params.cnWood);

    if (remove_stem_fraction < 1.)

    if (remove_stem_fraction < 1.)

        mNumberOfSnags[poolIndex(tree->dbh())]++;

        mNumberOfSnags[poolIndex(tree->dbh())]++;

}

}