WebSVN - public iLand - Diff - Rev 552 and 553


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/snag.cpp':
#include "snag.h"
#include "tree.h"
#include "species.h"
#include "globalsettings.h"
#include "expression.h"
// for calculation of climate decomposition
#include "resourceunit.h"
#include "watercycle.h"
#include "climate.h"
#include "model.h"

/** @class Snag
  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
  is subsequently forwarded to the soil sub model.
  Carbon is stored in three classes (depending on the size)
  The Snag dynamics class uses the following species parameter:
  cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW

  */
// static variables
double Snag::mDBHLower = -1.;
double Snag::mDBHHigher = 0.;
double Snag::mCarbonThreshold[] = {0., 0., 0.};

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

/// add biomass and weigh the parameter_value with the current C-content of the pool
void CNPool::addBiomass(const double biomass, const double CNratio, const double parameter_value)
{
    if (biomass==0.)
        return;
    double new_c = biomass*biomassCFraction;
    double p_old = C / (new_c + C);
    mParameter = mParameter*p_old + parameter_value*(1.-p_old);
    CNPair::addBiomass(biomass, CNratio);
}

// increase pool (and weigh the value)
void CNPool::operator+=(const CNPool &s)
{
    if (s.C==0.)
        return;
    mParameter = parameter(s); // calculate weighted parameter
    C+=s.C;
    N+=s.N;
}

double CNPool::parameter(const CNPool &s) const
{
    if (s.C==0.)
        return parameter();
    double p_old = C / (s.C + C);
    double result =  mParameter*p_old + s.parameter()*(1.-p_old);
    return result;
}


void Snag::setupThresholds(const double lower, const double upper)
{
    if (mDBHLower == lower)
        return;
    mDBHLower = lower;
    mDBHHigher = upper;
    mCarbonThreshold[0] = lower / 2.;
    mCarbonThreshold[1] = lower + (upper - lower)/2.;
    mCarbonThreshold[2] = upper + (upper - lower)/2.;
    //# threshold levels for emptying out the dbh-snag-classes
    //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%
    //# values in kg!
    for (int i=0;i<3;i++)
        mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1;
}


Snag::Snag()
{
    mRU = 0;
    CNPair::setCFraction(biomassCFraction);
}

void Snag::setup( const ResourceUnit *ru)
{
    mRU = ru;
    mClimateFactor = 0.;
    // branches
    mBranchCounter=0;
    for (int i=0;i<3;i++) {
        mTimeSinceDeath[i] = 0.;
        mNumberOfSnags[i] = 0.;
        mAvgDbh[i] = 0.;
        mAvgHeight[i] = 0.;
        mAvgVolume[i] = 0.;
        mKSW[i] = 0.;
        mCurrentKSW[i] = 0.;
    }
    mTotalSnagCarbon = 0.;
    if (mDBHLower<=0)
        throw IException("Snag::setupThresholds() not called or called with invalid parameters.");
}

// debug outputs
QList<QVariant> Snag::debugList()
{
    // list columns
    // for three pools
    QList<QVariant> list;

    // totals
    list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N;
    // fluxes to labile soil pool and to refractory soil pool
    list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N;

    for (int i=0;i<3;i++) {
        // 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 << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i];
    }

    // branch/coarse wood pools (5 yrs)
    for (int i=0;i<5;i++) {
        list << mOtherWood[i].C << mOtherWood[i].N;
    }
//    list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N
//            << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N
//            << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N
//            << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N
//            << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;
    return list;
}

void Snag::newYear()
{
    for (int i=0;i<3;i++) {
        mToSWD[i].clear(); // clear transfer pools to standing-woody-debris
        mCurrentKSW[i] = 0.;
    }
    mLabileFlux.clear();
    mRefractoryFlux.clear();
    mTotalToAtm.clear();
    mTotalToExtern.clear();
    mTotalIn.clear();
    mSWDtoSoil.clear();
}

/// calculate the dynamic climate modifier for decomposition 're'
/// calculation is done on the level of ResourceUnit because the water content per day is needed.
double Snag::calculateClimateFactors()
{
 
    double ft, fw;
 
    double f_sum = 0.;
    int iday=0;
    // calculate the water-factor for each month (see Adair et al 2008)
    double fw_month[12];
    double ratio;
    for (int m=0;m<12;m++) {
        if (mRU->waterCycle()->potentialEvapotranspiration()[m]>0.)
            ratio = mRU->climate()->precipitationMonth()[m] /  mRU->waterCycle()->potentialEvapotranspiration()[m];
        else
            ratio = 0;
        fw_month[m] = 1. / (1. + 30.*exp(-8.5*ratio));
        qDebug() <<"month"<< m << "PET" << mRU->waterCycle()->potentialEvapotranspiration()[m] << "prec" <<mRU->climate()->precipitationMonth()[m];
    }

    for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday)
    {
 
 
        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 = fw_month[day->month-1];
 
 

        f_sum += ft*fw;
    }
    // the climate factor is defined as the arithmentic annual mean value
    mClimateFactor = f_sum / double(mRU->climate()->daysOfYear());
    return mClimateFactor;
}

/// do the yearly calculation
/// see http://iland.boku.ac.at/snag+dynamics
void Snag::calculateYear()
{
    mSWDtoSoil.clear();
    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)
    if (isEmpty()) // nothing to do
        return;

    // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
    mRefractoryFlux+=mOtherWood[mBranchCounter];

    mOtherWood[mBranchCounter].clear();
    mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0.
    // decay of branches/coarse roots
    for (int i=0;i<5;i++) {
        if (mOtherWood[i].C>0.) {
            double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value...
            mOtherWood[i].C *= survive_rate;
        }
    }

    // process standing snags.
    // the input of the current year is in the mToSWD-Pools
    for (int i=0;i<3;i++) {

        // update the swd-pool with this years' input
        if (!mToSWD[i].isEmpty()) {
            // 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));
            //move content to the SWD pool
            mSWD[i] += mToSWD[i];
        }

        if (mSWD[i].C > 0) {
            // 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
            double survive_rate = exp(-mKSW[i] *climate_factor_re * 1. ); // 1: timestep
            mTotalToAtm.C += mSWD[i].C * (1. - survive_rate);
            mSWD[i].C *= survive_rate;

            // transition to downed woody debris
            // 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
            // 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
            // 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)
            // 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
            // calculate the transition probability of SWD to downed dead wood

            double half_life = mHalfLife[i] / climate_factor_re;
            double rate = -M_LN2 / half_life; // M_LN2: math. constant

            // higher decay rate for the class with smallest diameters
            if (i==0)
                rate*=2.;

            double transfer = 1. - exp(rate);

            // calculate flow to soil pool...
            mSWDtoSoil += mSWD[i] * transfer;
            mRefractoryFlux += mSWD[i] * transfer;
            mSWD[i] *= (1.-transfer); // reduce pool
            // calculate the stem number of remaining snags
            mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer);

            mTimeSinceDeath[i] += 1.;
            // 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
            // (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]) {
                // clear the pool: add the rest to the soil, clear statistics of the pool
                mRefractoryFlux += mSWD[i];
                mSWDtoSoil += mSWD[i];
                mSWD[i].clear();
                mAvgDbh[i] = 0.;
                mAvgHeight[i] = 0.;
                mAvgVolume[i] = 0.;
                mKSW[i] = 0.;
                mCurrentKSW[i] = 0.;
                mHalfLife[i] = 0.;
                mTimeSinceDeath[i] = 0.;
            }

        }

    }
    // total carbon in the snag-container on the RU *after* processing is the content of the
    // standing woody debris pools + the branches
    mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C +
                       mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C;
}

/// foliage and fineroot litter is transferred during tree growth.
void Snag::addTurnoverLitter(const Tree *tree, const double litter_foliage, const double litter_fineroot)
{
    mLabileFlux.addBiomass(litter_foliage, tree->species()->cnFoliage(), tree->species()->snagKyl());
    mLabileFlux.addBiomass(litter_fineroot, tree->species()->cnFineroot(), tree->species()->snagKyl());
}

/// after the death of the tree the five biomass compartments are processed.
void Snag::addMortality(const Tree *tree)
{
    const Species *species = tree->species();

    // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool
    mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl());
    mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());

    // branches and coarse roots are equally distributed over five years:
    double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2;
    for (int i=0;i<5; i++)
        mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());

    // just for book-keeping: keep track of all inputs into branches / roots / swd
    mTotalIn.addBiomass(tree->biomassBranch() + tree->biomassCoarseRoot() + tree->biomassStem(), species->cnWood());
    // 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

    // update statistics - stemnumber-weighted averages
    // 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_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;
    mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new;
    mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new;
    mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new;
    mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new;

    // 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)
    if (tree->biomassStem()==0)
        throw IException("Snag::addMortality: tree without stem biomass!!");
    p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
    p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
    mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new;
    mNumberOfSnags[pi]++;

    // finally add the biomass
    CNPool &to_swd = mToSWD[pi];
    to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr());
}

/// 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)
/// records on harvested biomass is collected (mTotalToExtern-pool).
void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction )
{
    const Species *species = tree->species();

    // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool
    mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl());
    mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());

    // for branches, add all biomass that remains in the forest to the soil
    mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr());
    // the same treatment for stem residuals
    mRefractoryFlux.addBiomass(tree->biomassStem() * remove_stem_fraction, species->cnWood(), tree->species()->snagKyr());

    // split the corase wood biomass into parts (slower decay)
    double biomass_rest = (tree->biomassCoarseRoot()) * 0.2;
    for (int i=0;i<5; i++)
        mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());


    // for book-keeping...
    mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction +
                              tree->biomassBranch()*remove_branch_fraction +
                              tree->biomassStem()*remove_stem_fraction, species->cnWood());
}




 

Rev 552	Rev 553
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/snag.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/snag.cpp':
2	#include "snag.h"	2	#include "snag.h"
3	#include "tree.h"	3	#include "tree.h"
4	#include "species.h"	4	#include "species.h"
5	#include "globalsettings.h"	5	#include "globalsettings.h"
6	#include "expression.h"	6	#include "expression.h"
7	// for calculation of climate decomposition	7	// for calculation of climate decomposition
8	#include "resourceunit.h"	8	#include "resourceunit.h"
9	#include "watercycle.h"	9	#include "watercycle.h"
10	#include "climate.h"	10	#include "climate.h"
11	#include "model.h"	11	#include "model.h"
12		12
13	/** @class Snag	13	/** @class Snag
14	Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.	14	Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
15	Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags	15	Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
16	is subsequently forwarded to the soil sub model.	16	is subsequently forwarded to the soil sub model.
17	Carbon is stored in three classes (depending on the size)	17	Carbon is stored in three classes (depending on the size)
18	The Snag dynamics class uses the following species parameter:	18	The Snag dynamics class uses the following species parameter:
19	cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW	19	cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW
20		20
21	*/	21	*/
22	// static variables	22	// static variables
23	double Snag::mDBHLower = -1.;	23	double Snag::mDBHLower = -1.;
24	double Snag::mDBHHigher = 0.;	24	double Snag::mDBHHigher = 0.;
25	double Snag::mCarbonThreshold[] = {0., 0., 0.};	25	double Snag::mCarbonThreshold[] = {0., 0., 0.};
26		26
27	double CNPair::biomassCFraction = biomassCFraction; // get global from globalsettings.h	27	double CNPair::biomassCFraction = biomassCFraction; // get global from globalsettings.h
28		28
29	/// add biomass and weigh the parameter_value with the current C-content of the pool	29	/// add biomass and weigh the parameter_value with the current C-content of the pool
30	void CNPool::addBiomass(const double biomass, const double CNratio, const double parameter_value)	30	void CNPool::addBiomass(const double biomass, const double CNratio, const double parameter_value)
31	{	31	{
32	if (biomass==0.)	32	if (biomass==0.)
33	return;	33	return;
34	double new_c = biomass*biomassCFraction;	34	double new_c = biomass*biomassCFraction;
35	double p_old = C / (new_c + C);	35	double p_old = C / (new_c + C);
36	mParameter = mParameterp_old + parameter_value(1.-p_old);	36	mParameter = mParameterp_old + parameter_value(1.-p_old);
37	CNPair::addBiomass(biomass, CNratio);	37	CNPair::addBiomass(biomass, CNratio);
38	}	38	}
39		39
40	// increase pool (and weigh the value)	40	// increase pool (and weigh the value)
41	void CNPool::operator+=(const CNPool &s)	41	void CNPool::operator+=(const CNPool &s)
42	{	42	{
43	if (s.C==0.)	43	if (s.C==0.)
44	return;	44	return;
45	mParameter = parameter(s); // calculate weighted parameter	45	mParameter = parameter(s); // calculate weighted parameter
46	C+=s.C;	46	C+=s.C;
47	N+=s.N;	47	N+=s.N;
48	}	48	}
49		49
50	double CNPool::parameter(const CNPool &s) const	50	double CNPool::parameter(const CNPool &s) const
51	{	51	{
52	if (s.C==0.)	52	if (s.C==0.)
53	return parameter();	53	return parameter();
54	double p_old = C / (s.C + C);	54	double p_old = C / (s.C + C);
55	double result = mParameterp_old + s.parameter()(1.-p_old);	55	double result = mParameterp_old + s.parameter()(1.-p_old);
56	return result;	56	return result;
57	}	57	}
58		58
59		59
60	void Snag::setupThresholds(const double lower, const double upper)	60	void Snag::setupThresholds(const double lower, const double upper)
61	{	61	{
62	if (mDBHLower == lower)	62	if (mDBHLower == lower)
63	return;	63	return;
64	mDBHLower = lower;	64	mDBHLower = lower;
65	mDBHHigher = upper;	65	mDBHHigher = upper;
66	mCarbonThreshold[0] = lower / 2.;	66	mCarbonThreshold[0] = lower / 2.;
67	mCarbonThreshold[1] = lower + (upper - lower)/2.;	67	mCarbonThreshold[1] = lower + (upper - lower)/2.;
68	mCarbonThreshold[2] = upper + (upper - lower)/2.;	68	mCarbonThreshold[2] = upper + (upper - lower)/2.;
69	//# threshold levels for emptying out the dbh-snag-classes	69	//# threshold levels for emptying out the dbh-snag-classes
70	//# derived from Psme woody allometry, converted to C, with a threshold level set to 10%	70	//# derived from Psme woody allometry, converted to C, with a threshold level set to 10%
71	//# values in kg!	71	//# values in kg!
72	for (int i=0;i<3;i++)	72	for (int i=0;i<3;i++)
73	mCarbonThreshold[i] = 0.10568pow(mCarbonThreshold[i],2.4247)0.5*0.1;	73	mCarbonThreshold[i] = 0.10568pow(mCarbonThreshold[i],2.4247)0.5*0.1;
74	}	74	}
75		75
76		76
77	Snag::Snag()	77	Snag::Snag()
78	{	78	{
79	mRU = 0;	79	mRU = 0;
80	CNPair::setCFraction(biomassCFraction);	80	CNPair::setCFraction(biomassCFraction);
81	}	81	}
82		82
83	void Snag::setup( const ResourceUnit *ru)	83	void Snag::setup( const ResourceUnit *ru)
84	{	84	{
85	mRU = ru;	85	mRU = ru;
86	mClimateFactor = 0.;	86	mClimateFactor = 0.;
87	// branches	87	// branches
88	mBranchCounter=0;	88	mBranchCounter=0;
89	for (int i=0;i<3;i++) {	89	for (int i=0;i<3;i++) {
90	mTimeSinceDeath[i] = 0.;	90	mTimeSinceDeath[i] = 0.;
91	mNumberOfSnags[i] = 0.;	91	mNumberOfSnags[i] = 0.;
92	mAvgDbh[i] = 0.;	92	mAvgDbh[i] = 0.;
93	mAvgHeight[i] = 0.;	93	mAvgHeight[i] = 0.;
94	mAvgVolume[i] = 0.;	94	mAvgVolume[i] = 0.;
95	mKSW[i] = 0.;	95	mKSW[i] = 0.;
96	mCurrentKSW[i] = 0.;	96	mCurrentKSW[i] = 0.;
97	}	97	}
98	mTotalSnagCarbon = 0.;	98	mTotalSnagCarbon = 0.;
99	if (mDBHLower<=0)	99	if (mDBHLower<=0)
100	throw IException("Snag::setupThresholds() not called or called with invalid parameters.");	100	throw IException("Snag::setupThresholds() not called or called with invalid parameters.");
101	}	101	}
102		102
103	// debug outputs	103	// debug outputs
104	QList<QVariant> Snag::debugList()	104	QList<QVariant> Snag::debugList()
105	{	105	{
106	// list columns	106	// list columns
107	// for three pools	107	// for three pools
108	QList<QVariant> list;	108	QList<QVariant> list;
109		109
110	// totals	110	// totals
111	list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N;	111	list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N;
112	// fluxes to labile soil pool and to refractory soil pool	112	// fluxes to labile soil pool and to refractory soil pool
113	list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N;	113	list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N;
114		114
115	for (int i=0;i<3;i++) {	115	for (int i=0;i<3;i++) {
116	// pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"	116	// pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"
117	list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N;	117	list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N;
118	list << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i];	118	list << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i];
119	}	119	}
120		120
121	// branch/coarse wood pools (5 yrs)	121	// branch/coarse wood pools (5 yrs)
122	for (int i=0;i<5;i++) {	122	for (int i=0;i<5;i++) {
123	list << mOtherWood[i].C << mOtherWood[i].N;	123	list << mOtherWood[i].C << mOtherWood[i].N;
124	}	124	}
125	// list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N	125	// list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N
126	// << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N	126	// << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N
127	// << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N	127	// << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N
128	// << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N	128	// << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N
129	// << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;	129	// << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;
130	return list;	130	return list;
131	}	131	}
132		132
133	void Snag::newYear()	133	void Snag::newYear()
134	{	134	{
135	for (int i=0;i<3;i++) {	135	for (int i=0;i<3;i++) {
136	mToSWD[i].clear(); // clear transfer pools to standing-woody-debris	136	mToSWD[i].clear(); // clear transfer pools to standing-woody-debris
137	mCurrentKSW[i] = 0.;	137	mCurrentKSW[i] = 0.;
138	}	138	}
139	mLabileFlux.clear();	139	mLabileFlux.clear();
140	mRefractoryFlux.clear();	140	mRefractoryFlux.clear();
141	mTotalToAtm.clear();	141	mTotalToAtm.clear();
142	mTotalToExtern.clear();	142	mTotalToExtern.clear();
143	mTotalIn.clear();	143	mTotalIn.clear();
144	mSWDtoSoil.clear();	144	mSWDtoSoil.clear();
145	}	145	}
146		146
147	/// calculate the dynamic climate modifier for decomposition 're'	147	/// calculate the dynamic climate modifier for decomposition 're'
148	/// calculation is done on the level of ResourceUnit because the water content per day is needed.	148	/// calculation is done on the level of ResourceUnit because the water content per day is needed.
149	double Snag::calculateClimateFactors()	149	double Snag::calculateClimateFactors()
150	{	150	{
151	double deficit;	-
152	double ft, fw;	151	double ft, fw;
153	const double top_layer_content = mRU->waterCycle()->topLayerWaterContent();	-
154	double f_sum = 0.;	152	double f_sum = 0.;
155	int iday=0;	153	int iday=0;
-		154	// calculate the water-factor for each month (see Adair et al 2008)
-		155	double fw_month[12];
-		156	double ratio;
-		157	for (int m=0;m<12;m++) {
-		158	if (mRU->waterCycle()->potentialEvapotranspiration()[m]>0.)
-		159	ratio = mRU->climate()->precipitationMonth()[m] / mRU->waterCycle()->potentialEvapotranspiration()[m];
-		160	else
-		161	ratio = 0;
-		162	fw_month[m] = 1. / (1. + 30.exp(-8.5ratio));
-		163	qDebug() <<"month"<< m << "PET" << mRU->waterCycle()->potentialEvapotranspiration()[m] << "prec" <<mRU->climate()->precipitationMonth()[m];
-		164	}
-		165
156	for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday)	166	for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday)
157	{	167	{
158	deficit = mRU->waterCycle()->waterDeficit_mm(iday);	-
159		-
160	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)	168	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)
161	fw = 1. - ~~limit(deficit~~ / top_layer_content, 0., 1.);	169	fw = fw_month[day->month-1];
162	// the water effect: depends on the water deficit; if the deficit is higher than the parameterized	-
163	// content of the top layer (where most microbial activity is located), than then fw gets 0.	-
164		170
165	f_sum += ft*fw;	171	f_sum += ft*fw;
166	}	172	}
167	// the climate factor is defined as the arithmentic annual mean value	173	// the climate factor is defined as the arithmentic annual mean value
168	mClimateFactor = f_sum / double(mRU->climate()->daysOfYear());	174	mClimateFactor = f_sum / double(mRU->climate()->daysOfYear());
169	return mClimateFactor;	175	return mClimateFactor;
170	}	176	}
171		177
172	/// do the yearly calculation	178	/// do the yearly calculation
173	/// see http://iland.boku.ac.at/snag+dynamics	179	/// see http://iland.boku.ac.at/snag+dynamics
174	void Snag::calculateYear()	180	void Snag::calculateYear()
175	{	181	{
176	mSWDtoSoil.clear();	182	mSWDtoSoil.clear();
177	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)	183	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)
178	if (isEmpty()) // nothing to do	184	if (isEmpty()) // nothing to do
179	return;	185	return;
180		186
181	// process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool	187	// process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
182	mRefractoryFlux+=mOtherWood[mBranchCounter];	188	mRefractoryFlux+=mOtherWood[mBranchCounter];
183		189
184	mOtherWood[mBranchCounter].clear();	190	mOtherWood[mBranchCounter].clear();
185	mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0.	191	mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0.
186	// decay of branches/coarse roots	192	// decay of branches/coarse roots
187	for (int i=0;i<5;i++) {	193	for (int i=0;i<5;i++) {
188	if (mOtherWood[i].C>0.) {	194	if (mOtherWood[i].C>0.) {
189	double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value...	195	double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value...
190	mOtherWood[i].C *= survive_rate;	196	mOtherWood[i].C *= survive_rate;
191	}	197	}
192	}	198	}
193		199
194	// process standing snags.	200	// process standing snags.
195	// the input of the current year is in the mToSWD-Pools	201	// the input of the current year is in the mToSWD-Pools
196	for (int i=0;i<3;i++) {	202	for (int i=0;i<3;i++) {
197		203
198	// update the swd-pool with this years' input	204	// update the swd-pool with this years' input
199	if (!mToSWD[i].isEmpty()) {	205	if (!mToSWD[i].isEmpty()) {
200	// update decay rate (apply average yearly input to the state parameters)	206	// update decay rate (apply average yearly input to the state parameters)
201	mKSW[i] = mKSW[i](mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i](mToSWD[i].C/(mSWD[i].C+mToSWD[i].C));	207	mKSW[i] = mKSW[i](mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i](mToSWD[i].C/(mSWD[i].C+mToSWD[i].C));
202	//move content to the SWD pool	208	//move content to the SWD pool
203	mSWD[i] += mToSWD[i];	209	mSWD[i] += mToSWD[i];
204	}	210	}
205		211
206	if (mSWD[i].C > 0) {	212	if (mSWD[i].C > 0) {
207	// reduce the Carbon (note: the N stays, thus the CN ratio changes)	213	// reduce the Carbon (note: the N stays, thus the CN ratio changes)
208	// use the decay rate that is derived as a weighted average of all standing woody debris	214	// use the decay rate that is derived as a weighted average of all standing woody debris
209	double survive_rate = exp(-mKSW[i] climate_factor_re 1. ); // 1: timestep	215	double survive_rate = exp(-mKSW[i] climate_factor_re 1. ); // 1: timestep
210	mTotalToAtm.C += mSWD[i].C * (1. - survive_rate);	216	mTotalToAtm.C += mSWD[i].C * (1. - survive_rate);
211	mSWD[i].C *= survive_rate;	217	mSWD[i].C *= survive_rate;
212		218
213	// transition to downed woody debris	219	// transition to downed woody debris
214	// update: use negative exponential decay, species parameter: half-life	220	// update: use negative exponential decay, species parameter: half-life
215	// modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa	221	// modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa
216	// this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),	222	// this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),
217	// where lag times for cohort transitions are linearly modified with re although here individual good or bad years have	223	// where lag times for cohort transitions are linearly modified with re although here individual good or bad years have
218	// an immediate effect, the average climatic influence should come through (and it is inherently transient)	224	// an immediate effect, the average climatic influence should come through (and it is inherently transient)
219	// note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)	225	// note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)
220	// needs to be calculated, followed by a weighted update of the previous swd.hl.	226	// needs to be calculated, followed by a weighted update of the previous swd.hl.
221	// As weights here we use stem number, as the processes here pertain individual snags	227	// As weights here we use stem number, as the processes here pertain individual snags
222	// calculate the transition probability of SWD to downed dead wood	228	// calculate the transition probability of SWD to downed dead wood
223		229
224	double half_life = mHalfLife[i] / climate_factor_re;	230	double half_life = mHalfLife[i] / climate_factor_re;
225	double rate = -M_LN2 / half_life; // M_LN2: math. constant	231	double rate = -M_LN2 / half_life; // M_LN2: math. constant
226		232
227	// higher decay rate for the class with smallest diameters	233	// higher decay rate for the class with smallest diameters
228	if (i==0)	234	if (i==0)
229	rate*=2.;	235	rate*=2.;
230		236
231	double transfer = 1. - exp(rate);	237	double transfer = 1. - exp(rate);
232		238
233	// calculate flow to soil pool...	239	// calculate flow to soil pool...
234	mSWDtoSoil += mSWD[i] * transfer;	240	mSWDtoSoil += mSWD[i] * transfer;
235	mRefractoryFlux += mSWD[i] * transfer;	241	mRefractoryFlux += mSWD[i] * transfer;
236	mSWD[i] *= (1.-transfer); // reduce pool	242	mSWD[i] *= (1.-transfer); // reduce pool
237	// calculate the stem number of remaining snags	243	// calculate the stem number of remaining snags
238	mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer);	244	mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer);
239		245
240	mTimeSinceDeath[i] += 1.;	246	mTimeSinceDeath[i] += 1.;
241	// if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats	247	// if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats
242	// also, if the Carbon of an average snag is less than 10% of the original average tree	248	// also, if the Carbon of an average snag is less than 10% of the original average tree
243	// (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD	249	// (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD
244	if (mNumberOfSnags[i] < 0.5 \|\| mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) {	250	if (mNumberOfSnags[i] < 0.5 \|\| mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) {
245	// clear the pool: add the rest to the soil, clear statistics of the pool	251	// clear the pool: add the rest to the soil, clear statistics of the pool
246	mRefractoryFlux += mSWD[i];	252	mRefractoryFlux += mSWD[i];
247	mSWDtoSoil += mSWD[i];	253	mSWDtoSoil += mSWD[i];
248	mSWD[i].clear();	254	mSWD[i].clear();
249	mAvgDbh[i] = 0.;	255	mAvgDbh[i] = 0.;
250	mAvgHeight[i] = 0.;	256	mAvgHeight[i] = 0.;
251	mAvgVolume[i] = 0.;	257	mAvgVolume[i] = 0.;
252	mKSW[i] = 0.;	258	mKSW[i] = 0.;
253	mCurrentKSW[i] = 0.;	259	mCurrentKSW[i] = 0.;
254	mHalfLife[i] = 0.;	260	mHalfLife[i] = 0.;
255	mTimeSinceDeath[i] = 0.;	261	mTimeSinceDeath[i] = 0.;
256	}	262	}
257		263
258	}	264	}
259		265
260	}	266	}
261	// total carbon in the snag-container on the RU after processing is the content of the	267	// total carbon in the snag-container on the RU after processing is the content of the
262	// standing woody debris pools + the branches	268	// standing woody debris pools + the branches
263	mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C +	269	mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C +
264	mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C;	270	mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C;
265	}	271	}
266		272
267	/// foliage and fineroot litter is transferred during tree growth.	273	/// foliage and fineroot litter is transferred during tree growth.
268	void Snag::addTurnoverLitter(const Tree *tree, const double litter_foliage, const double litter_fineroot)	274	void Snag::addTurnoverLitter(const Tree *tree, const double litter_foliage, const double litter_fineroot)
269	{	275	{
270	mLabileFlux.addBiomass(litter_foliage, tree->species()->cnFoliage(), tree->species()->snagKyl());	276	mLabileFlux.addBiomass(litter_foliage, tree->species()->cnFoliage(), tree->species()->snagKyl());
271	mLabileFlux.addBiomass(litter_fineroot, tree->species()->cnFineroot(), tree->species()->snagKyl());	277	mLabileFlux.addBiomass(litter_fineroot, tree->species()->cnFineroot(), tree->species()->snagKyl());
272	}	278	}
273		279
274	/// after the death of the tree the five biomass compartments are processed.	280	/// after the death of the tree the five biomass compartments are processed.
275	void Snag::addMortality(const Tree *tree)	281	void Snag::addMortality(const Tree *tree)
276	{	282	{
277	const Species *species = tree->species();	283	const Species *species = tree->species();
278		284
279	// immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool	285	// immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool
280	mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl());	286	mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl());
281	mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());	287	mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
282		288
283	// branches and coarse roots are equally distributed over five years:	289	// branches and coarse roots are equally distributed over five years:
284	double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2;	290	double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2;
285	for (int i=0;i<5; i++)	291	for (int i=0;i<5; i++)
286	mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());	292	mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
287		293
288	// just for book-keeping: keep track of all inputs into branches / roots / swd	294	// just for book-keeping: keep track of all inputs into branches / roots / swd
289	mTotalIn.addBiomass(tree->biomassBranch() + tree->biomassCoarseRoot() + tree->biomassStem(), species->cnWood());	295	mTotalIn.addBiomass(tree->biomassBranch() + tree->biomassCoarseRoot() + tree->biomassStem(), species->cnWood());
290	// stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool	296	// stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
291	int pi = poolIndex(tree->dbh()); // get right transfer pool	297	int pi = poolIndex(tree->dbh()); // get right transfer pool
292		298
293	// update statistics - stemnumber-weighted averages	299	// update statistics - stemnumber-weighted averages
294	// note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)	300	// note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
295	double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers)	301	double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers)
296	double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1).	302	double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1).
297	mAvgDbh[pi] = mAvgDbh[pi]p_old + tree->dbh()p_new;	303	mAvgDbh[pi] = mAvgDbh[pi]p_old + tree->dbh()p_new;
298	mAvgHeight[pi] = mAvgHeight[pi]p_old + tree->height()p_new;	304	mAvgHeight[pi] = mAvgHeight[pi]p_old + tree->height()p_new;
299	mAvgVolume[pi] = mAvgVolume[pi]p_old + tree->volume()p_new;	305	mAvgVolume[pi] = mAvgVolume[pi]p_old + tree->volume()p_new;
300	mTimeSinceDeath[pi] = mTimeSinceDeath[pi]p_old + 1.p_new;	306	mTimeSinceDeath[pi] = mTimeSinceDeath[pi]p_old + 1.p_new;
301	mHalfLife[pi] = mHalfLife[pi]p_old + species->snagHalflife() p_new;	307	mHalfLife[pi] = mHalfLife[pi]p_old + species->snagHalflife() p_new;
302		308
303	// average the decay rate (ksw); this is done based on the carbon content	309	// average the decay rate (ksw); this is done based on the carbon content
304	// aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)	310	// aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
305	if (tree->biomassStem()==0)	311	if (tree->biomassStem()==0)
306	throw IException("Snag::addMortality: tree without stem biomass!!");	312	throw IException("Snag::addMortality: tree without stem biomass!!");
307	p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);	313	p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
308	p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);	314	p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
309	mCurrentKSW[pi] = mCurrentKSW[pi]p_old + species->snagKsw() p_new;	315	mCurrentKSW[pi] = mCurrentKSW[pi]p_old + species->snagKsw() p_new;
310	mNumberOfSnags[pi]++;	316	mNumberOfSnags[pi]++;
311		317
312	// finally add the biomass	318	// finally add the biomass
313	CNPool &to_swd = mToSWD[pi];	319	CNPool &to_swd = mToSWD[pi];
314	to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr());	320	to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr());
315	}	321	}
316		322
317	/// add residual biomass of 'tree' after harvesting.	323	/// add residual biomass of 'tree' after harvesting.
318	/// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is removed by the harvest operation (i.e.: not to stay in the system)	324	/// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is removed by the harvest operation (i.e.: not to stay in the system)
319	/// records on harvested biomass is collected (mTotalToExtern-pool).	325	/// records on harvested biomass is collected (mTotalToExtern-pool).
320	void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction )	326	void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction )
321	{	327	{
322	const Species *species = tree->species();	328	const Species *species = tree->species();
323		329
324	// immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool	330	// immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool
325	mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl());	331	mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl());
326	mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());	332	mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
327		333
328	// for branches, add all biomass that remains in the forest to the soil	334	// for branches, add all biomass that remains in the forest to the soil
329	mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr());	335	mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr());
330	// the same treatment for stem residuals	336	// the same treatment for stem residuals
331	mRefractoryFlux.addBiomass(tree->biomassStem() * remove_stem_fraction, species->cnWood(), tree->species()->snagKyr());	337	mRefractoryFlux.addBiomass(tree->biomassStem() * remove_stem_fraction, species->cnWood(), tree->species()->snagKyr());
332		338
333	// split the corase wood biomass into parts (slower decay)	339	// split the corase wood biomass into parts (slower decay)
334	double biomass_rest = (tree->biomassCoarseRoot()) * 0.2;	340	double biomass_rest = (tree->biomassCoarseRoot()) * 0.2;
335	for (int i=0;i<5; i++)	341	for (int i=0;i<5; i++)
336	mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());	342	mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
337		343
338		344
339	// for book-keeping...	345	// for book-keeping...
340	mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction +	346	mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction +
341	tree->biomassBranch()*remove_branch_fraction +	347	tree->biomassBranch()*remove_branch_fraction +
342	tree->biomassStem()*remove_stem_fraction, species->cnWood());	348	tree->biomassStem()*remove_stem_fraction, species->cnWood());
343	}	349	}
344		350
345		351
346		352
347		353
348		354

Subversion Repositories public iLand

(root)/src/core/snag.cpp @ 1222 - Rev 552 → 553