WebSVN - public iLand - Diff - Rev 1202 and 1203


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':
/********************************************************************************************
**    iLand - an individual based forest landscape and disturbance model
**    http://iland.boku.ac.at
**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl
**
**    This program is free software: you can redistribute it and/or modify
**    it under the terms of the GNU General Public License as published by
**    the Free Software Foundation, either version 3 of the License, or
**    (at your option) any later version.
**
**    This program is distributed in the hope that it will be useful,
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
**    GNU General Public License for more details.
**
**    You should have received a copy of the GNU General Public License
**    along with this program.  If not, see <http://www.gnu.org/licenses/>.
********************************************************************************************/

/** @class ResourceUnit
  ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
  (Climate, Soil, Water, ...).
  @ingroup core
  A resource unit has a size of (currently) 100x100m. Many processes in iLand operate on the level of a ResourceUnit.
  Each resource unit has the same Climate and other properties (e.g. available nitrogen).
  Proceses on this level are, inter alia, NPP Production (see Production3PG), water calculations (WaterCycle), the modeling
  of dead trees (Snag) and soil processes (Soil).

  */
#include <QtCore>
#include "global.h"

#include "resourceunit.h"
#include "resourceunitspecies.h"
#include "speciesset.h"
#include "species.h"
#include "production3pg.h"
#include "model.h"
#include "climate.h"
#include "watercycle.h"
#include "snag.h"
#include "soil.h"
#include "helper.h"

ResourceUnit::~ResourceUnit()
{
    if (mWater)
        delete mWater;
    mWater = 0;
    if (mSnag)
        delete mSnag;
    if (mSoil)
        delete mSoil;

    qDeleteAll(mRUSpecies);

    if (mSaplings)
        delete[] mSaplings;

    mSnag = 0;
    mSoil = 0;
    mSaplings = 0;
}

ResourceUnit::ResourceUnit(const int index)
{
    qDeleteAll(mRUSpecies);
    mSpeciesSet = 0;
    mClimate = 0;
    mPixelCount=0;
    mStockedArea = 0;
    mStockedPixelCount = 0;
    mStockableArea = 0;
    mAggregatedWLA = 0.;
    mAggregatedLA = 0.;
    mAggregatedLR = 0.;
    mEffectiveArea = 0.;
    mLRI_modification = 0.;
    mIndex = index;
    mSaplingHeightMap = 0;
    mEffectiveArea_perWLA = 0.;
    mWater = new WaterCycle();
    mSnag = 0;
    mSoil = 0;
    mSaplings = 0;
    mID = 0;
}

void ResourceUnit::setup()
{
    mWater->setup(this);

    if (mSnag)
        delete mSnag;
    mSnag=0;
    if (mSoil)
        delete mSoil;
    mSoil=0;
    if (Model::settings().carbonCycleEnabled) {
        mSoil = new Soil(this);
        mSnag = new Snag;
        mSnag->setup(this);
        const XmlHelper &xml=GlobalSettings::instance()->settings();

        // setup contents of the soil of the RU; use values for C and N (kg/ha)
        mSoil->setInitialState(CNPool(xml.valueDouble("model.site.youngLabileC", -1),
                                      xml.valueDouble("model.site.youngLabileN", -1),
                                      xml.valueDouble("model.site.youngLabileDecompRate", -1)),
                               CNPool(xml.valueDouble("model.site.youngRefractoryC", -1),
                                      xml.valueDouble("model.site.youngRefractoryN", -1),
                                      xml.valueDouble("model.site.youngRefractoryDecompRate", -1)),
                               CNPair(xml.valueDouble("model.site.somC", -1), xml.valueDouble("model.site.somN", -1)));
    }

    if (mSaplings)
        delete mSaplings;
    if (Model::settings().regenerationEnabled) {
        mSaplings = new SaplingCell[cPxPerHectare];
    }

    // setup variables
    mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40);

    // if dynamic coupling of soil nitrogen is enabled, a starting value for available N is calculated
    if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) {
        mSoil->setClimateFactor(1.);
        mSoil->calculateYear();
        mUnitVariables.nitrogenAvailable = soil()->availableNitrogen();
    }
    mHasDeadTrees = false;
    mAverageAging = 0.;

}
void ResourceUnit::setBoundingBox(const QRectF &bb)
{
    mBoundingBox = bb;
    mCornerOffset = GlobalSettings::instance()->model()->grid()->indexAt(bb.topLeft());
}

/// return the sapling cell at given LIF-coordinates
SaplingCell *ResourceUnit::saplingCell(const QPoint &lifCoords) const
{
    // LIF-Coordinates are global, we here need (RU-)local coordinates
    int ix = lifCoords.x() % cPxPerRU;
    int iy = lifCoords.y() % cPxPerRU;
    int i = iy*cPxPerRU+ix;
    Q_ASSERT(i>=0 && i<cPxPerHectare);
    return &mSaplings[i];
}

/// set species and setup the species-per-RU-data
void ResourceUnit::setSpeciesSet(SpeciesSet *set)
{
    mSpeciesSet = set;
    qDeleteAll(mRUSpecies);

    //mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
    for (int i=0;i<set->count();i++) {
        Species *s = const_cast<Species*>(mSpeciesSet->species(i));
        if (!s)
            throw IException("ResourceUnit::setSpeciesSet: invalid index!");

        ResourceUnitSpecies *rus = new ResourceUnitSpecies();
        mRUSpecies.push_back(rus);
        rus->setup(s, this);
        /* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
           If the container memory is relocated (QVector), the pointer gets invalid!!!
           Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
        //mRUSpecies[i].setup(s,this); // setup this element

    }
}

ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)
{
    return *mRUSpecies[species->index()];
}

const ResourceUnitSpecies *ResourceUnit::constResourceUnitSpecies(const Species *species) const
{
    return mRUSpecies[species->index()];
}

Tree &ResourceUnit::newTree()
{
    // start simple: just append to the vector...
    if (mTrees.isEmpty())
        mTrees.reserve(100); // reserve a junk of memory for trees

    mTrees.append(Tree());
    return mTrees.back();
}
int ResourceUnit::newTreeIndex()
{
    newTree();
    return mTrees.count()-1; // return index of the last tree
}

/// remove dead trees from tree list
/// reduce size of vector if lots of space is free
/// tests showed that this way of cleanup is very fast,
/// because no memory allocations are performed (simple memmove())
/// when trees are moved.
void ResourceUnit::cleanTreeList()
{
    if (!mHasDeadTrees)
        return;

    QVector<Tree>::iterator last=mTrees.end()-1;
    QVector<Tree>::iterator current = mTrees.begin();
    while (last>=current && (*last).isDead())
        --last;

    while (current<last) {
        if ((*current).isDead()) {
            *current = *last; // copy data!
            --last; //
            while (last>=current && (*last).isDead())
                --last;
        }
        ++current;
    }
    ++last; // last points now to the first dead tree

    // free ressources
    if (last!=mTrees.end()) {
        mTrees.erase(last, mTrees.end());
        if (mTrees.capacity()>100) {
            if (mTrees.count() / double(mTrees.capacity()) < 0.2) {
                //int target_size = mTrees.count()*2;
                //qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
                //mTrees.reserve(qMax(target_size, 100));
                if (logLevelDebug())
                    qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count();
                mTrees.squeeze();
            }
        }
    }
    mHasDeadTrees = false; // reset flag
}

void ResourceUnit::newYear()
{
    mAggregatedWLA = 0.;
    mAggregatedLA = 0.;
    mAggregatedLR = 0.;
    mEffectiveArea = 0.;
    mPixelCount = mStockedPixelCount = 0;
    snagNewYear();
    if (mSoil)
        mSoil->newYear();
    // clear statistics global and per species...
    QList<ResourceUnitSpecies*>::const_iterator i;
    QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();
    mStatistics.clear();
    for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
        (*i)->statisticsDead().clear();
        (*i)->statisticsMgmt().clear();
    }

}

/** production() is the "stand-level" part of the biomass production (3PG).
    - The amount of radiation intercepted by the stand is calculated
    - the water cycle is calculated
    - statistics for each species are cleared
    - The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production)
    see also: http://iland.boku.ac.at/individual+tree+light+availability */
void ResourceUnit::production()
{

    if (mAggregatedWLA==0. || mPixelCount==0) {
        // clear statistics of resourceunitspecies
        for ( QList<ResourceUnitSpecies*>::const_iterator i=mRUSpecies.constBegin(); i!=mRUSpecies.constEnd(); ++i)
            (*i)->statistics().clear();
        mEffectiveArea = 0.;
        mStockedArea = 0.;
        return;
    }

    // the pixel counters are filled during the height-grid-calculations
    mStockedArea = cHeightPerRU*cHeightPerRU * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)
    if (leafAreaIndex()<3.) {
        // estimate stocked area based on crown projections
        double crown_area = 0.;
        for (int i=0;i<mTrees.count();++i)
            crown_area += mTrees.at(i).isDead() ? 0. : mTrees.at(i).stamp()->reader()->crownArea();

        if (logLevelDebug())
            qDebug() << "crown area: lai" << leafAreaIndex() << "stocked area (pixels)" << mStockedArea << " area (crown)" << crown_area;
        if (leafAreaIndex()<1.) {
            mStockedArea = std::min(crown_area, mStockedArea);
        } else {
            // for LAI between 1 and 3:
            // interpolate between sum of crown area of trees (at LAI=1) and the pixel-based value (at LAI=3 and above)
            double px_frac = (leafAreaIndex()-1.)/2.; // 0 at LAI=1, 1 at LAI=3
            mStockedArea = mStockedArea * px_frac + std::min(crown_area, mStockedArea) * (1. - px_frac);
        }
        if (mStockedArea==0.)
            return;
    }

    // calculate the leaf area index (LAI)
    double LAI = mAggregatedLA / mStockedArea;
    // calculate the intercepted radiation fraction using the law of Beer Lambert
    const double k = Model::settings().lightExtinctionCoefficient;
    double interception_fraction = 1. - exp(-k * LAI);
    mEffectiveArea = mStockedArea * interception_fraction; // m2

    // calculate the total weighted leaf area on this RU:
    mLRI_modification = interception_fraction *  mStockedArea / mAggregatedWLA; // p_WLA
    if (mLRI_modification == 0.)
        qDebug() << "lri modifaction==0!";

    if (logLevelDebug()) {
    DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")
            .arg(LAI)
            .arg(interception_fraction)
            .arg(mLRI_modification)
            .arg(mStockedArea);
    );
    }

    // calculate LAI fractions
    QList<ResourceUnitSpecies*>::const_iterator i;
    QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();
    double ru_lai = leafAreaIndex();
    if (ru_lai < 1.)
        ru_lai = 1.;
    // note: LAIFactors are only 1 if sum of LAI is > 1. (see WaterCycle)
    for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
        double lai_factor = (*i)->statistics().leafAreaIndex() / ru_lai;

        //DBGMODE(
        if (lai_factor > 1.) {
                        const ResourceUnitSpecies* rus=*i;
                        qDebug() << "LAI factor > 1: species ru-index:" << rus->species()->name() << rus->ru()->index();
                    }
        //);
        (*i)->setLAIfactor( lai_factor );
    }

    // soil water model - this determines soil water contents needed for response calculations
    {
    mWater->run();
    }

    // invoke species specific calculation (3PG)
    for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
        //DBGMODE(
        if ((*i)->LAIfactor() > 1.) {
                    const ResourceUnitSpecies* rus=*i;
                    qDebug() << "LAI factor > 1: species ru-index value:" << rus->species()->name() << rus->ru()->index() << rus->LAIfactor();
                    }
        //);
        (*i)->calculate(); // CALCULATE 3PG

        // debug output related to production
        if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dStandGPP) && (*i)->LAIfactor()>0.) {
            DebugList &out = GlobalSettings::instance()->debugList(index(), GlobalSettings::dStandGPP);
            out << (*i)->species()->id() << index() << id();
            out << (*i)->LAIfactor() << (*i)->prod3PG().GPPperArea() << productiveArea()*(*i)->LAIfactor()*(*i)->prod3PG().GPPperArea() << averageAging() << (*i)->prod3PG().fEnvYear() ;

        }
    }
}

void ResourceUnit::calculateInterceptedArea()
{
    if (mAggregatedLR==0) {
        mEffectiveArea_perWLA = 0.;
        return;
    }
    Q_ASSERT(mAggregatedLR>0.);
    mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;
    if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR  << "eff.area./wla:" << mEffectiveArea_perWLA;
}

// function is called immediately before the growth of individuals
void ResourceUnit::beforeGrow()
{
    mAverageAging = 0.;
}

// function is called after finishing the indivdual growth / mortality.
void ResourceUnit::afterGrow()
{
    mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees)
    if (mAverageAging>0. && mAverageAging<0.00001)
        qDebug() << "ru" << mIndex << "aging <0.00001";
    if (mAverageAging<0. || mAverageAging>1.)
        qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();
}

void ResourceUnit::yearEnd()
{
    // calculate statistics for all tree species of the ressource unit
    int c = mRUSpecies.count();
    for (int i=0;i<c; i++) {
        mRUSpecies[i]->statisticsDead().calculate(); // calculate the dead trees
        mRUSpecies[i]->statisticsMgmt().calculate(); // stats of removed trees
        mRUSpecies[i]->updateGWL(); // get sum of dead trees (died + removed)
        mRUSpecies[i]->statistics().calculate(); // calculate the living (and add removed volume to gwl)
        mStatistics.add(mRUSpecies[i]->statistics());
    }
    mStatistics.calculate(); // aggreagte on stand level

    // update carbon flows
    if (soil() && GlobalSettings::instance()->model()->settings().carbonCycleEnabled) {
        double area_factor = stockableArea() / cRUArea; //conversion factor
        mUnitVariables.carbonUptake = statistics().npp() * biomassCFraction;
        mUnitVariables.carbonUptake += statistics().nppSaplings() * biomassCFraction;

        double to_atm = snag()->fluxToAtmosphere().C / area_factor; // from snags, kgC/ha
        to_atm += soil()->fluxToAtmosphere().C *cRUArea/10.; // soil: t/ha -> t/m2 -> kg/ha
        mUnitVariables.carbonToAtm = to_atm;

        double to_dist = snag()->fluxToDisturbance().C / area_factor;
        to_dist += soil()->fluxToDisturbance().C * cRUArea/10.;
        double to_harvest = snag()->fluxToExtern().C / area_factor;

        mUnitVariables.NEP = mUnitVariables.carbonUptake - to_atm - to_dist - to_harvest; // kgC/ha

        // incremental values....
        mUnitVariables.cumCarbonUptake += mUnitVariables.carbonUptake;
        mUnitVariables.cumCarbonToAtm += mUnitVariables.carbonToAtm;
        mUnitVariables.cumNEP += mUnitVariables.NEP;

    }

}

void ResourceUnit::addTreeAgingForAllTrees()
{
    mAverageAging = 0.;
    foreach(const Tree &t, mTrees) {
        addTreeAging(t.leafArea(), t.species()->aging(t.height(), t.age()));
    }

}

/// refresh of tree based statistics.
/// WARNING: this function is only called once (during startup).
/// see function "yearEnd()" above!!!
void ResourceUnit::createStandStatistics()
{
    // clear statistics (ru-level and ru-species level)
    mStatistics.clear();
    for (int i=0;i<mRUSpecies.count();i++) {
        mRUSpecies[i]->statistics().clear();
        mRUSpecies[i]->statisticsDead().clear();
        mRUSpecies[i]->statisticsMgmt().clear();
        mRUSpecies[i]->saplingStat().clearStatistics();
    }

    // add all trees to the statistics objects of the species
    foreach(const Tree &t, mTrees) {
        if (!t.isDead())
            resourceUnitSpecies(t.species()).statistics().add(&t, 0);
    }
    // summarise sapling stats
    GlobalSettings::instance()->model()->saplings()->calculateInitialStatistics(this);

    // summarize statistics for the whole resource unit
    for (int i=0;i<mRUSpecies.count();i++) {
        mRUSpecies[i]->saplingStat().calculate(mRUSpecies[i]->species(), this);
        mRUSpecies[i]->statistics().add(&mRUSpecies[i]->saplingStat());
        mRUSpecies[i]->statistics().calculate();
        mStatistics.add(mRUSpecies[i]->statistics());
    }
    mStatistics.calculate();
    mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*stockableArea()):0.;
    if (mAverageAging<0. || mAverageAging>1.)
        qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();

}

/** recreate statistics. This is necessary after events that changed the structure
    of the stand *after* the growth of trees (where stand statistics are updated).
    An example is after disturbances.  */
void ResourceUnit::recreateStandStatistics(bool recalculate_stats)
{
    // when called after disturbances (recalculate_stats=false), we
    // clear only the tree-specific variables in the stats (i.e. we keep NPP, and regen carbon),
    // and then re-add all trees (since TreeGrowthData is NULL no NPP is available).
    // The statistics are not summarised here, because this happens for all resource units
    // in the yearEnd function of RU.
    for (int i=0;i<mRUSpecies.count();i++) {
        if (recalculate_stats)
            mRUSpecies[i]->statistics().clear();
        else
            mRUSpecies[i]->statistics().clearOnlyTrees();
    }
    foreach(const Tree &t, mTrees) {
        resourceUnitSpecies(t.species()).statistics().add(&t, 0);
    }

    if (recalculate_stats) {
        for (int i=0;i<mRUSpecies.count();i++) {
            mRUSpecies[i]->statistics().calculate();
        }
    }
}




void ResourceUnit::calculateCarbonCycle()
{
    if (!snag())
        return;

    // (1) calculate the snag dynamics
    // because all carbon/nitrogen-flows from trees to the soil are routed through the snag-layer,
    // all soil inputs (litter + deadwood) are collected in the Snag-object.
    snag()->calculateYear();
    soil()->setClimateFactor( snag()->climateFactor() ); // the climate factor is only calculated once
    soil()->setSoilInput( snag()->labileFlux(), snag()->refractoryFlux());
    soil()->calculateYear(); // update the ICBM/2N model
    // use available nitrogen?
    if (Model::settings().useDynamicAvailableNitrogen)
        mUnitVariables.nitrogenAvailable = soil()->availableNitrogen();

    // debug output
    if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dCarbonCycle) && !snag()->isEmpty()) {
        DebugList &out = GlobalSettings::instance()->debugList(index(), GlobalSettings::dCarbonCycle);
        out << index() << id(); // resource unit index and id
        out << snag()->debugList(); // snag debug outs
        out << soil()->debugList(); // ICBM/2N debug outs
    }

}


 

Rev 1202	Rev 1203
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':
2	/********************************************************************************************	2	/********************************************************************************************
3	** iLand - an individual based forest landscape and disturbance model	3	** iLand - an individual based forest landscape and disturbance model
4	** http://iland.boku.ac.at	4	** http://iland.boku.ac.at
5	** Copyright (C) 2009- Werner Rammer, Rupert Seidl	5	** Copyright (C) 2009- Werner Rammer, Rupert Seidl
6	**	6	**
7	** This program is free software: you can redistribute it and/or modify	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	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	9	** the Free Software Foundation, either version 3 of the License, or
10	** (at your option) any later version.	10	** (at your option) any later version.
11	**	11	**
12	** This program is distributed in the hope that it will be useful,	12	** This program is distributed in the hope that it will be useful,
13	** but WITHOUT ANY WARRANTY; without even the implied warranty of	13	** but WITHOUT ANY WARRANTY; without even the implied warranty of
14	** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the	14	** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15	** GNU General Public License for more details.	15	** GNU General Public License for more details.
16	**	16	**
17	** You should have received a copy of the GNU General Public License	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/>.	18	** along with this program. If not, see <http://www.gnu.org/licenses/>.
19	********************************************************************************************/	19	********************************************************************************************/
20		20
21	/** @class ResourceUnit	21	/** @class ResourceUnit
22	ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components	22	ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
23	(Climate, Soil, Water, ...).	23	(Climate, Soil, Water, ...).
24	@ingroup core	24	@ingroup core
25	A resource unit has a size of (currently) 100x100m. Many processes in iLand operate on the level of a ResourceUnit.	25	A resource unit has a size of (currently) 100x100m. Many processes in iLand operate on the level of a ResourceUnit.
26	Each resource unit has the same Climate and other properties (e.g. available nitrogen).	26	Each resource unit has the same Climate and other properties (e.g. available nitrogen).
27	Proceses on this level are, inter alia, NPP Production (see Production3PG), water calculations (WaterCycle), the modeling	27	Proceses on this level are, inter alia, NPP Production (see Production3PG), water calculations (WaterCycle), the modeling
28	of dead trees (Snag) and soil processes (Soil).	28	of dead trees (Snag) and soil processes (Soil).
29		29
30	*/	30	*/
31	#include <QtCore>	31	#include <QtCore>
32	#include "global.h"	32	#include "global.h"
33		33
34	#include "resourceunit.h"	34	#include "resourceunit.h"
35	#include "resourceunitspecies.h"	35	#include "resourceunitspecies.h"
36	#include "speciesset.h"	36	#include "speciesset.h"
37	#include "species.h"	37	#include "species.h"
38	#include "production3pg.h"	38	#include "production3pg.h"
39	#include "model.h"	39	#include "model.h"
40	#include "climate.h"	40	#include "climate.h"
41	#include "watercycle.h"	41	#include "watercycle.h"
42	#include "snag.h"	42	#include "snag.h"
43	#include "soil.h"	43	#include "soil.h"
44	#include "helper.h"	44	#include "helper.h"
45		45
46	ResourceUnit::~ResourceUnit()	46	ResourceUnit::~ResourceUnit()
47	{	47	{
48	if (mWater)	48	if (mWater)
49	delete mWater;	49	delete mWater;
50	mWater = 0;	50	mWater = 0;
51	if (mSnag)	51	if (mSnag)
52	delete mSnag;	52	delete mSnag;
53	if (mSoil)	53	if (mSoil)
54	delete mSoil;	54	delete mSoil;
55		55
56	qDeleteAll(mRUSpecies);	56	qDeleteAll(mRUSpecies);
57		57
58	if (mSaplings)	58	if (mSaplings)
59	delete[] mSaplings;	59	delete[] mSaplings;
60		60
61	mSnag = 0;	61	mSnag = 0;
62	mSoil = 0;	62	mSoil = 0;
63	mSaplings = 0;	63	mSaplings = 0;
64	}	64	}
65		65
66	ResourceUnit::ResourceUnit(const int index)	66	ResourceUnit::ResourceUnit(const int index)
67	{	67	{
68	qDeleteAll(mRUSpecies);	68	qDeleteAll(mRUSpecies);
69	mSpeciesSet = 0;	69	mSpeciesSet = 0;
70	mClimate = 0;	70	mClimate = 0;
71	mPixelCount=0;	71	mPixelCount=0;
72	mStockedArea = 0;	72	mStockedArea = 0;
73	mStockedPixelCount = 0;	73	mStockedPixelCount = 0;
74	mStockableArea = 0;	74	mStockableArea = 0;
75	mAggregatedWLA = 0.;	75	mAggregatedWLA = 0.;
76	mAggregatedLA = 0.;	76	mAggregatedLA = 0.;
77	mAggregatedLR = 0.;	77	mAggregatedLR = 0.;
78	mEffectiveArea = 0.;	78	mEffectiveArea = 0.;
79	mLRI_modification = 0.;	79	mLRI_modification = 0.;
80	mIndex = index;	80	mIndex = index;
81	mSaplingHeightMap = 0;	81	mSaplingHeightMap = 0;
82	mEffectiveArea_perWLA = 0.;	82	mEffectiveArea_perWLA = 0.;
83	mWater = new WaterCycle();	83	mWater = new WaterCycle();
84	mSnag = 0;	84	mSnag = 0;
85	mSoil = 0;	85	mSoil = 0;
86	mSaplings = 0;	86	mSaplings = 0;
87	mID = 0;	87	mID = 0;
88	}	88	}
89		89
90	void ResourceUnit::setup()	90	void ResourceUnit::setup()
91	{	91	{
92	mWater->setup(this);	92	mWater->setup(this);
93		93
94	if (mSnag)	94	if (mSnag)
95	delete mSnag;	95	delete mSnag;
96	mSnag=0;	96	mSnag=0;
97	if (mSoil)	97	if (mSoil)
98	delete mSoil;	98	delete mSoil;
99	mSoil=0;	99	mSoil=0;
100	if (Model::settings().carbonCycleEnabled) {	100	if (Model::settings().carbonCycleEnabled) {
101	mSoil = new Soil(this);	101	mSoil = new Soil(this);
102	mSnag = new Snag;	102	mSnag = new Snag;
103	mSnag->setup(this);	103	mSnag->setup(this);
104	const XmlHelper &xml=GlobalSettings::instance()->settings();	104	const XmlHelper &xml=GlobalSettings::instance()->settings();
105		105
106	// setup contents of the soil of the RU; use values for C and N (kg/ha)	106	// setup contents of the soil of the RU; use values for C and N (kg/ha)
107	mSoil->setInitialState(CNPool(xml.valueDouble("model.site.youngLabileC", -1),	107	mSoil->setInitialState(CNPool(xml.valueDouble("model.site.youngLabileC", -1),
108	xml.valueDouble("model.site.youngLabileN", -1),	108	xml.valueDouble("model.site.youngLabileN", -1),
109	xml.valueDouble("model.site.youngLabileDecompRate", -1)),	109	xml.valueDouble("model.site.youngLabileDecompRate", -1)),
110	CNPool(xml.valueDouble("model.site.youngRefractoryC", -1),	110	CNPool(xml.valueDouble("model.site.youngRefractoryC", -1),
111	xml.valueDouble("model.site.youngRefractoryN", -1),	111	xml.valueDouble("model.site.youngRefractoryN", -1),
112	xml.valueDouble("model.site.youngRefractoryDecompRate", -1)),	112	xml.valueDouble("model.site.youngRefractoryDecompRate", -1)),
113	CNPair(xml.valueDouble("model.site.somC", -1), xml.valueDouble("model.site.somN", -1)));	113	CNPair(xml.valueDouble("model.site.somC", -1), xml.valueDouble("model.site.somN", -1)));
114	}	114	}
115		115
116	if (mSaplings)	116	if (mSaplings)
117	delete mSaplings;	117	delete mSaplings;
118	if (Model::settings().regenerationEnabled) {	118	if (Model::settings().regenerationEnabled) {
119	mSaplings = new SaplingCell[cPxPerHectare];	119	mSaplings = new SaplingCell[cPxPerHectare];
120	}	120	}
121		121
122	// setup variables	122	// setup variables
123	mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40);	123	mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40);
124		124
125	// if dynamic coupling of soil nitrogen is enabled, a starting value for available N is calculated	125	// if dynamic coupling of soil nitrogen is enabled, a starting value for available N is calculated
126	if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) {	126	if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) {
127	mSoil->setClimateFactor(1.);	127	mSoil->setClimateFactor(1.);
128	mSoil->calculateYear();	128	mSoil->calculateYear();
129	mUnitVariables.nitrogenAvailable = soil()->availableNitrogen();	129	mUnitVariables.nitrogenAvailable = soil()->availableNitrogen();
130	}	130	}
131	mHasDeadTrees = false;	131	mHasDeadTrees = false;
132	mAverageAging = 0.;	132	mAverageAging = 0.;
133		133
134	}	134	}
135	void ResourceUnit::setBoundingBox(const QRectF &bb)	135	void ResourceUnit::setBoundingBox(const QRectF &bb)
136	{	136	{
137	mBoundingBox = bb;	137	mBoundingBox = bb;
138	mCornerOffset = GlobalSettings::instance()->model()->grid()->indexAt(bb.topLeft());	138	mCornerOffset = GlobalSettings::instance()->model()->grid()->indexAt(bb.topLeft());
-		139	}
-		140
-		141	/// return the sapling cell at given LIF-coordinates
-		142	SaplingCell *ResourceUnit::saplingCell(const QPoint &lifCoords) const
-		143	{
-		144	// LIF-Coordinates are global, we here need (RU-)local coordinates
-		145	int ix = lifCoords.x() % cPxPerRU;
-		146	int iy = lifCoords.y() % cPxPerRU;
-		147	int i = iy*cPxPerRU+ix;
-		148	Q_ASSERT(i>=0 && i<cPxPerHectare);
-		149	return &mSaplings[i];
139	}	150	}
140		151
141	/// set species and setup the species-per-RU-data	152	/// set species and setup the species-per-RU-data
142	void ResourceUnit::setSpeciesSet(SpeciesSet *set)	153	void ResourceUnit::setSpeciesSet(SpeciesSet *set)
143	{	154	{
144	mSpeciesSet = set;	155	mSpeciesSet = set;
145	qDeleteAll(mRUSpecies);	156	qDeleteAll(mRUSpecies);
146		157
147	//mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated	158	//mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
148	for (int i=0;i<set->count();i++) {	159	for (int i=0;i<set->count();i++) {
149	Species s = const_cast<Species>(mSpeciesSet->species(i));	160	Species s = const_cast<Species>(mSpeciesSet->species(i));
150	if (!s)	161	if (!s)
151	throw IException("ResourceUnit::setSpeciesSet: invalid index!");	162	throw IException("ResourceUnit::setSpeciesSet: invalid index!");
152		163
153	ResourceUnitSpecies *rus = new ResourceUnitSpecies();	164	ResourceUnitSpecies *rus = new ResourceUnitSpecies();
154	mRUSpecies.push_back(rus);	165	mRUSpecies.push_back(rus);
155	rus->setup(s, this);	166	rus->setup(s, this);
156	/* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.	167	/* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
157	If the container memory is relocated (QVector), the pointer gets invalid!!!	168	If the container memory is relocated (QVector), the pointer gets invalid!!!
158	Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */	169	Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
159	//mRUSpecies[i].setup(s,this); // setup this element	170	//mRUSpecies[i].setup(s,this); // setup this element
160		171
161	}	172	}
162	}	173	}
163		174
164	ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)	175	ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)
165	{	176	{
166	return *mRUSpecies[species->index()];	177	return *mRUSpecies[species->index()];
167	}	178	}
168		179
169	const ResourceUnitSpecies ResourceUnit::constResourceUnitSpecies(const Species species) const	180	const ResourceUnitSpecies ResourceUnit::constResourceUnitSpecies(const Species species) const
170	{	181	{
171	return mRUSpecies[species->index()];	182	return mRUSpecies[species->index()];
172	}	183	}
173		184
174	Tree &ResourceUnit::newTree()	185	Tree &ResourceUnit::newTree()
175	{	186	{
176	// start simple: just append to the vector...	187	// start simple: just append to the vector...
177	if (mTrees.isEmpty())	188	if (mTrees.isEmpty())
178	mTrees.reserve(100); // reserve a junk of memory for trees	189	mTrees.reserve(100); // reserve a junk of memory for trees
179		190
180	mTrees.append(Tree());	191	mTrees.append(Tree());
181	return mTrees.back();	192	return mTrees.back();
182	}	193	}
183	int ResourceUnit::newTreeIndex()	194	int ResourceUnit::newTreeIndex()
184	{	195	{
185	newTree();	196	newTree();
186	return mTrees.count()-1; // return index of the last tree	197	return mTrees.count()-1; // return index of the last tree
187	}	198	}
188		199
189	/// remove dead trees from tree list	200	/// remove dead trees from tree list
190	/// reduce size of vector if lots of space is free	201	/// reduce size of vector if lots of space is free
191	/// tests showed that this way of cleanup is very fast,	202	/// tests showed that this way of cleanup is very fast,
192	/// because no memory allocations are performed (simple memmove())	203	/// because no memory allocations are performed (simple memmove())
193	/// when trees are moved.	204	/// when trees are moved.
194	void ResourceUnit::cleanTreeList()	205	void ResourceUnit::cleanTreeList()
195	{	206	{
196	if (!mHasDeadTrees)	207	if (!mHasDeadTrees)
197	return;	208	return;
198		209
199	QVector<Tree>::iterator last=mTrees.end()-1;	210	QVector<Tree>::iterator last=mTrees.end()-1;
200	QVector<Tree>::iterator current = mTrees.begin();	211	QVector<Tree>::iterator current = mTrees.begin();
201	while (last>=current && (*last).isDead())	212	while (last>=current && (*last).isDead())
202	--last;	213	--last;
203		214
204	while (current<last) {	215	while (current<last) {
205	if ((*current).isDead()) {	216	if ((*current).isDead()) {
206	current = last; // copy data!	217	current = last; // copy data!
207	--last; //	218	--last; //
208	while (last>=current && (*last).isDead())	219	while (last>=current && (*last).isDead())
209	--last;	220	--last;
210	}	221	}
211	++current;	222	++current;
212	}	223	}
213	++last; // last points now to the first dead tree	224	++last; // last points now to the first dead tree
214		225
215	// free ressources	226	// free ressources
216	if (last!=mTrees.end()) {	227	if (last!=mTrees.end()) {
217	mTrees.erase(last, mTrees.end());	228	mTrees.erase(last, mTrees.end());
218	if (mTrees.capacity()>100) {	229	if (mTrees.capacity()>100) {
219	if (mTrees.count() / double(mTrees.capacity()) < 0.2) {	230	if (mTrees.count() / double(mTrees.capacity()) < 0.2) {
220	//int target_size = mTrees.count()*2;	231	//int target_size = mTrees.count()*2;
221	//qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;	232	//qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
222	//mTrees.reserve(qMax(target_size, 100));	233	//mTrees.reserve(qMax(target_size, 100));
223	if (logLevelDebug())	234	if (logLevelDebug())
224	qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count();	235	qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count();
225	mTrees.squeeze();	236	mTrees.squeeze();
226	}	237	}
227	}	238	}
228	}	239	}
229	mHasDeadTrees = false; // reset flag	240	mHasDeadTrees = false; // reset flag
230	}	241	}
231		242
232	void ResourceUnit::newYear()	243	void ResourceUnit::newYear()
233	{	244	{
234	mAggregatedWLA = 0.;	245	mAggregatedWLA = 0.;
235	mAggregatedLA = 0.;	246	mAggregatedLA = 0.;
236	mAggregatedLR = 0.;	247	mAggregatedLR = 0.;
237	mEffectiveArea = 0.;	248	mEffectiveArea = 0.;
238	mPixelCount = mStockedPixelCount = 0;	249	mPixelCount = mStockedPixelCount = 0;
239	snagNewYear();	250	snagNewYear();
240	if (mSoil)	251	if (mSoil)
241	mSoil->newYear();	252	mSoil->newYear();
242	// clear statistics global and per species...	253	// clear statistics global and per species...
243	QList<ResourceUnitSpecies*>::const_iterator i;	254	QList<ResourceUnitSpecies*>::const_iterator i;
244	QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();	255	QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();
245	mStatistics.clear();	256	mStatistics.clear();
246	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {	257	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
247	(*i)->statisticsDead().clear();	258	(*i)->statisticsDead().clear();
248	(*i)->statisticsMgmt().clear();	259	(*i)->statisticsMgmt().clear();
249	}	260	}
250		261
251	}	262	}
252		263
253	/** production() is the "stand-level" part of the biomass production (3PG).	264	/** production() is the "stand-level" part of the biomass production (3PG).
254	- The amount of radiation intercepted by the stand is calculated	265	- The amount of radiation intercepted by the stand is calculated
255	- the water cycle is calculated	266	- the water cycle is calculated
256	- statistics for each species are cleared	267	- statistics for each species are cleared
257	- The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production)	268	- The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production)
258	see also: http://iland.boku.ac.at/individual+tree+light+availability */	269	see also: http://iland.boku.ac.at/individual+tree+light+availability */
259	void ResourceUnit::production()	270	void ResourceUnit::production()
260	{	271	{
261		272
262	if (mAggregatedWLA==0. \|\| mPixelCount==0) {	273	if (mAggregatedWLA==0. \|\| mPixelCount==0) {
263	// clear statistics of resourceunitspecies	274	// clear statistics of resourceunitspecies
264	for ( QList<ResourceUnitSpecies*>::const_iterator i=mRUSpecies.constBegin(); i!=mRUSpecies.constEnd(); ++i)	275	for ( QList<ResourceUnitSpecies*>::const_iterator i=mRUSpecies.constBegin(); i!=mRUSpecies.constEnd(); ++i)
265	(*i)->statistics().clear();	276	(*i)->statistics().clear();
266	mEffectiveArea = 0.;	277	mEffectiveArea = 0.;
267	mStockedArea = 0.;	278	mStockedArea = 0.;
268	return;	279	return;
269	}	280	}
270		281
271	// the pixel counters are filled during the height-grid-calculations	282	// the pixel counters are filled during the height-grid-calculations
272	mStockedArea = cHeightPerRUcHeightPerRU mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)	283	mStockedArea = cHeightPerRUcHeightPerRU mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)
273	if (leafAreaIndex()<3.) {	284	if (leafAreaIndex()<3.) {
274	// estimate stocked area based on crown projections	285	// estimate stocked area based on crown projections
275	double crown_area = 0.;	286	double crown_area = 0.;
276	for (int i=0;i<mTrees.count();++i)	287	for (int i=0;i<mTrees.count();++i)
277	crown_area += mTrees.at(i).isDead() ? 0. : mTrees.at(i).stamp()->reader()->crownArea();	288	crown_area += mTrees.at(i).isDead() ? 0. : mTrees.at(i).stamp()->reader()->crownArea();
278		289
279	if (logLevelDebug())	290	if (logLevelDebug())
280	qDebug() << "crown area: lai" << leafAreaIndex() << "stocked area (pixels)" << mStockedArea << " area (crown)" << crown_area;	291	qDebug() << "crown area: lai" << leafAreaIndex() << "stocked area (pixels)" << mStockedArea << " area (crown)" << crown_area;
281	if (leafAreaIndex()<1.) {	292	if (leafAreaIndex()<1.) {
282	mStockedArea = std::min(crown_area, mStockedArea);	293	mStockedArea = std::min(crown_area, mStockedArea);
283	} else {	294	} else {
284	// for LAI between 1 and 3:	295	// for LAI between 1 and 3:
285	// interpolate between sum of crown area of trees (at LAI=1) and the pixel-based value (at LAI=3 and above)	296	// interpolate between sum of crown area of trees (at LAI=1) and the pixel-based value (at LAI=3 and above)
286	double px_frac = (leafAreaIndex()-1.)/2.; // 0 at LAI=1, 1 at LAI=3	297	double px_frac = (leafAreaIndex()-1.)/2.; // 0 at LAI=1, 1 at LAI=3
287	mStockedArea = mStockedArea * px_frac + std::min(crown_area, mStockedArea) * (1. - px_frac);	298	mStockedArea = mStockedArea * px_frac + std::min(crown_area, mStockedArea) * (1. - px_frac);
288	}	299	}
289	if (mStockedArea==0.)	300	if (mStockedArea==0.)
290	return;	301	return;
291	}	302	}
292		303
293	// calculate the leaf area index (LAI)	304	// calculate the leaf area index (LAI)
294	double LAI = mAggregatedLA / mStockedArea;	305	double LAI = mAggregatedLA / mStockedArea;
295	// calculate the intercepted radiation fraction using the law of Beer Lambert	306	// calculate the intercepted radiation fraction using the law of Beer Lambert
296	const double k = Model::settings().lightExtinctionCoefficient;	307	const double k = Model::settings().lightExtinctionCoefficient;
297	double interception_fraction = 1. - exp(-k * LAI);	308	double interception_fraction = 1. - exp(-k * LAI);
298	mEffectiveArea = mStockedArea * interception_fraction; // m2	309	mEffectiveArea = mStockedArea * interception_fraction; // m2
299		310
300	// calculate the total weighted leaf area on this RU:	311	// calculate the total weighted leaf area on this RU:
301	mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA; // p_WLA	312	mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA; // p_WLA
302	if (mLRI_modification == 0.)	313	if (mLRI_modification == 0.)
303	qDebug() << "lri modifaction==0!";	314	qDebug() << "lri modifaction==0!";
304		315
305	if (logLevelDebug()) {	316	if (logLevelDebug()) {
306	DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")	317	DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")
307	.arg(LAI)	318	.arg(LAI)
308	.arg(interception_fraction)	319	.arg(interception_fraction)
309	.arg(mLRI_modification)	320	.arg(mLRI_modification)
310	.arg(mStockedArea);	321	.arg(mStockedArea);
311	);	322	);
312	}	323	}
313		324
314	// calculate LAI fractions	325	// calculate LAI fractions
315	QList<ResourceUnitSpecies*>::const_iterator i;	326	QList<ResourceUnitSpecies*>::const_iterator i;
316	QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();	327	QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();
317	double ru_lai = leafAreaIndex();	328	double ru_lai = leafAreaIndex();
318	if (ru_lai < 1.)	329	if (ru_lai < 1.)
319	ru_lai = 1.;	330	ru_lai = 1.;
320	// note: LAIFactors are only 1 if sum of LAI is > 1. (see WaterCycle)	331	// note: LAIFactors are only 1 if sum of LAI is > 1. (see WaterCycle)
321	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {	332	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
322	double lai_factor = (*i)->statistics().leafAreaIndex() / ru_lai;	333	double lai_factor = (*i)->statistics().leafAreaIndex() / ru_lai;
323		334
324	//DBGMODE(	335	//DBGMODE(
325	if (lai_factor > 1.) {	336	if (lai_factor > 1.) {
326	const ResourceUnitSpecies* rus=*i;	337	const ResourceUnitSpecies* rus=*i;
327	qDebug() << "LAI factor > 1: species ru-index:" << rus->species()->name() << rus->ru()->index();	338	qDebug() << "LAI factor > 1: species ru-index:" << rus->species()->name() << rus->ru()->index();
328	}	339	}
329	//);	340	//);
330	(*i)->setLAIfactor( lai_factor );	341	(*i)->setLAIfactor( lai_factor );
331	}	342	}
332		343
333	// soil water model - this determines soil water contents needed for response calculations	344	// soil water model - this determines soil water contents needed for response calculations
334	{	345	{
335	mWater->run();	346	mWater->run();
336	}	347	}
337		348
338	// invoke species specific calculation (3PG)	349	// invoke species specific calculation (3PG)
339	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {	350	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
340	//DBGMODE(	351	//DBGMODE(
341	if ((*i)->LAIfactor() > 1.) {	352	if ((*i)->LAIfactor() > 1.) {
342	const ResourceUnitSpecies* rus=*i;	353	const ResourceUnitSpecies* rus=*i;
343	qDebug() << "LAI factor > 1: species ru-index value:" << rus->species()->name() << rus->ru()->index() << rus->LAIfactor();	354	qDebug() << "LAI factor > 1: species ru-index value:" << rus->species()->name() << rus->ru()->index() << rus->LAIfactor();
344	}	355	}
345	//);	356	//);
346	(*i)->calculate(); // CALCULATE 3PG	357	(*i)->calculate(); // CALCULATE 3PG
347		358
348	// debug output related to production	359	// debug output related to production
349	if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dStandGPP) && (*i)->LAIfactor()>0.) {	360	if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dStandGPP) && (*i)->LAIfactor()>0.) {
350	DebugList &out = GlobalSettings::instance()->debugList(index(), GlobalSettings::dStandGPP);	361	DebugList &out = GlobalSettings::instance()->debugList(index(), GlobalSettings::dStandGPP);
351	out << (*i)->species()->id() << index() << id();	362	out << (*i)->species()->id() << index() << id();
352	out << (i)->LAIfactor() << (i)->prod3PG().GPPperArea() << productiveArea()(i)->LAIfactor()(i)->prod3PG().GPPperArea() << averageAging() << (*i)->prod3PG().fEnvYear() ;	363	out << (i)->LAIfactor() << (i)->prod3PG().GPPperArea() << productiveArea()(i)->LAIfactor()(i)->prod3PG().GPPperArea() << averageAging() << (*i)->prod3PG().fEnvYear() ;
353		364
354	}	365	}
355	}	366	}
356	}	367	}
357		368
358	void ResourceUnit::calculateInterceptedArea()	369	void ResourceUnit::calculateInterceptedArea()
359	{	370	{
360	if (mAggregatedLR==0) {	371	if (mAggregatedLR==0) {
361	mEffectiveArea_perWLA = 0.;	372	mEffectiveArea_perWLA = 0.;
362	return;	373	return;
363	}	374	}
364	Q_ASSERT(mAggregatedLR>0.);	375	Q_ASSERT(mAggregatedLR>0.);
365	mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;	376	mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;
366	if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA;	377	if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA;
367	}	378	}
368		379
369	// function is called immediately before the growth of individuals	380	// function is called immediately before the growth of individuals
370	void ResourceUnit::beforeGrow()	381	void ResourceUnit::beforeGrow()
371	{	382	{
372	mAverageAging = 0.;	383	mAverageAging = 0.;
373	}	384	}
374		385
375	// function is called after finishing the indivdual growth / mortality.	386	// function is called after finishing the indivdual growth / mortality.
376	void ResourceUnit::afterGrow()	387	void ResourceUnit::afterGrow()
377	{	388	{
378	mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees)	389	mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees)
379	if (mAverageAging>0. && mAverageAging<0.00001)	390	if (mAverageAging>0. && mAverageAging<0.00001)
380	qDebug() << "ru" << mIndex << "aging <0.00001";	391	qDebug() << "ru" << mIndex << "aging <0.00001";
381	if (mAverageAging<0. \|\| mAverageAging>1.)	392	if (mAverageAging<0. \|\| mAverageAging>1.)
382	qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();	393	qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();
383	}	394	}
384		395
385	void ResourceUnit::yearEnd()	396	void ResourceUnit::yearEnd()
386	{	397	{
387	// calculate statistics for all tree species of the ressource unit	398	// calculate statistics for all tree species of the ressource unit
388	int c = mRUSpecies.count();	399	int c = mRUSpecies.count();
389	for (int i=0;i<c; i++) {	400	for (int i=0;i<c; i++) {
390	mRUSpecies[i]->statisticsDead().calculate(); // calculate the dead trees	401	mRUSpecies[i]->statisticsDead().calculate(); // calculate the dead trees
391	mRUSpecies[i]->statisticsMgmt().calculate(); // stats of removed trees	402	mRUSpecies[i]->statisticsMgmt().calculate(); // stats of removed trees
392	mRUSpecies[i]->updateGWL(); // get sum of dead trees (died + removed)	403	mRUSpecies[i]->updateGWL(); // get sum of dead trees (died + removed)
393	mRUSpecies[i]->statistics().calculate(); // calculate the living (and add removed volume to gwl)	404	mRUSpecies[i]->statistics().calculate(); // calculate the living (and add removed volume to gwl)
394	mStatistics.add(mRUSpecies[i]->statistics());	405	mStatistics.add(mRUSpecies[i]->statistics());
395	}	406	}
396	mStatistics.calculate(); // aggreagte on stand level	407	mStatistics.calculate(); // aggreagte on stand level
397		408
398	// update carbon flows	409	// update carbon flows
399	if (soil() && GlobalSettings::instance()->model()->settings().carbonCycleEnabled) {	410	if (soil() && GlobalSettings::instance()->model()->settings().carbonCycleEnabled) {
400	double area_factor = stockableArea() / cRUArea; //conversion factor	411	double area_factor = stockableArea() / cRUArea; //conversion factor
401	mUnitVariables.carbonUptake = statistics().npp() * biomassCFraction;	412	mUnitVariables.carbonUptake = statistics().npp() * biomassCFraction;
402	mUnitVariables.carbonUptake += statistics().nppSaplings() * biomassCFraction;	413	mUnitVariables.carbonUptake += statistics().nppSaplings() * biomassCFraction;
403		414
404	double to_atm = snag()->fluxToAtmosphere().C / area_factor; // from snags, kgC/ha	415	double to_atm = snag()->fluxToAtmosphere().C / area_factor; // from snags, kgC/ha
405	to_atm += soil()->fluxToAtmosphere().C *cRUArea/10.; // soil: t/ha -> t/m2 -> kg/ha	416	to_atm += soil()->fluxToAtmosphere().C *cRUArea/10.; // soil: t/ha -> t/m2 -> kg/ha
406	mUnitVariables.carbonToAtm = to_atm;	417	mUnitVariables.carbonToAtm = to_atm;
407		418
408	double to_dist = snag()->fluxToDisturbance().C / area_factor;	419	double to_dist = snag()->fluxToDisturbance().C / area_factor;
409	to_dist += soil()->fluxToDisturbance().C * cRUArea/10.;	420	to_dist += soil()->fluxToDisturbance().C * cRUArea/10.;
410	double to_harvest = snag()->fluxToExtern().C / area_factor;	421	double to_harvest = snag()->fluxToExtern().C / area_factor;
411		422
412	mUnitVariables.NEP = mUnitVariables.carbonUptake - to_atm - to_dist - to_harvest; // kgC/ha	423	mUnitVariables.NEP = mUnitVariables.carbonUptake - to_atm - to_dist - to_harvest; // kgC/ha
413		424
414	// incremental values....	425	// incremental values....
415	mUnitVariables.cumCarbonUptake += mUnitVariables.carbonUptake;	426	mUnitVariables.cumCarbonUptake += mUnitVariables.carbonUptake;
416	mUnitVariables.cumCarbonToAtm += mUnitVariables.carbonToAtm;	427	mUnitVariables.cumCarbonToAtm += mUnitVariables.carbonToAtm;
417	mUnitVariables.cumNEP += mUnitVariables.NEP;	428	mUnitVariables.cumNEP += mUnitVariables.NEP;
418		429
419	}	430	}
420		431
421	}	432	}
422		433
423	void ResourceUnit::addTreeAgingForAllTrees()	434	void ResourceUnit::addTreeAgingForAllTrees()
424	{	435	{
425	mAverageAging = 0.;	436	mAverageAging = 0.;
426	foreach(const Tree &t, mTrees) {	437	foreach(const Tree &t, mTrees) {
427	addTreeAging(t.leafArea(), t.species()->aging(t.height(), t.age()));	438	addTreeAging(t.leafArea(), t.species()->aging(t.height(), t.age()));
428	}	439	}
429		440
430	}	441	}
431		442
432	/// refresh of tree based statistics.	443	/// refresh of tree based statistics.
433	/// WARNING: this function is only called once (during startup).	444	/// WARNING: this function is only called once (during startup).
434	/// see function "yearEnd()" above!!!	445	/// see function "yearEnd()" above!!!
435	void ResourceUnit::createStandStatistics()	446	void ResourceUnit::createStandStatistics()
436	{	447	{
437	// clear statistics (ru-level and ru-species level)	448	// clear statistics (ru-level and ru-species level)
438	mStatistics.clear();	449	mStatistics.clear();
439	for (int i=0;i<mRUSpecies.count();i++) {	450	for (int i=0;i<mRUSpecies.count();i++) {
440	mRUSpecies[i]->statistics().clear();	451	mRUSpecies[i]->statistics().clear();
441	mRUSpecies[i]->statisticsDead().clear();	452	mRUSpecies[i]->statisticsDead().clear();
442	mRUSpecies[i]->statisticsMgmt().clear();	453	mRUSpecies[i]->statisticsMgmt().clear();
443	mRUSpecies[i]->saplingStat().clearStatistics();	454	mRUSpecies[i]->saplingStat().clearStatistics();
444	}	455	}
445		456
446	// add all trees to the statistics objects of the species	457	// add all trees to the statistics objects of the species
447	foreach(const Tree &t, mTrees) {	458	foreach(const Tree &t, mTrees) {
448	if (!t.isDead())	459	if (!t.isDead())
449	resourceUnitSpecies(t.species()).statistics().add(&t, 0);	460	resourceUnitSpecies(t.species()).statistics().add(&t, 0);
450	}	461	}
451	// summarise sapling stats	462	// summarise sapling stats
452	GlobalSettings::instance()->model()->saplings()->calculateInitialStatistics(this);	463	GlobalSettings::instance()->model()->saplings()->calculateInitialStatistics(this);
453		464
454	// summarize statistics for the whole resource unit	465	// summarize statistics for the whole resource unit
455	for (int i=0;i<mRUSpecies.count();i++) {	466	for (int i=0;i<mRUSpecies.count();i++) {
456	mRUSpecies[i]->saplingStat().calculate(mRUSpecies[i]->species(), this);	467	mRUSpecies[i]->saplingStat().calculate(mRUSpecies[i]->species(), this);
457	mRUSpecies[i]->statistics().add(&mRUSpecies[i]->saplingStat());	468	mRUSpecies[i]->statistics().add(&mRUSpecies[i]->saplingStat());
458	mRUSpecies[i]->statistics().calculate();	469	mRUSpecies[i]->statistics().calculate();
459	mStatistics.add(mRUSpecies[i]->statistics());	470	mStatistics.add(mRUSpecies[i]->statistics());
460	}	471	}
461	mStatistics.calculate();	472	mStatistics.calculate();
462	mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*stockableArea()):0.;	473	mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*stockableArea()):0.;
463	if (mAverageAging<0. \|\| mAverageAging>1.)	474	if (mAverageAging<0. \|\| mAverageAging>1.)
464	qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();	475	qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();
465		476
466	}	477	}
467		478
468	/** recreate statistics. This is necessary after events that changed the structure	479	/** recreate statistics. This is necessary after events that changed the structure
469	of the stand after the growth of trees (where stand statistics are updated).	480	of the stand after the growth of trees (where stand statistics are updated).
470	An example is after disturbances. */	481	An example is after disturbances. */
471	void ResourceUnit::recreateStandStatistics(bool recalculate_stats)	482	void ResourceUnit::recreateStandStatistics(bool recalculate_stats)
472	{	483	{
473	// when called after disturbances (recalculate_stats=false), we	484	// when called after disturbances (recalculate_stats=false), we
474	// clear only the tree-specific variables in the stats (i.e. we keep NPP, and regen carbon),	485	// clear only the tree-specific variables in the stats (i.e. we keep NPP, and regen carbon),
475	// and then re-add all trees (since TreeGrowthData is NULL no NPP is available).	486	// and then re-add all trees (since TreeGrowthData is NULL no NPP is available).
476	// The statistics are not summarised here, because this happens for all resource units	487	// The statistics are not summarised here, because this happens for all resource units
477	// in the yearEnd function of RU.	488	// in the yearEnd function of RU.
478	for (int i=0;i<mRUSpecies.count();i++) {	489	for (int i=0;i<mRUSpecies.count();i++) {
479	if (recalculate_stats)	490	if (recalculate_stats)
480	mRUSpecies[i]->statistics().clear();	491	mRUSpecies[i]->statistics().clear();
481	else	492	else
482	mRUSpecies[i]->statistics().clearOnlyTrees();	493	mRUSpecies[i]->statistics().clearOnlyTrees();
483	}	494	}
484	foreach(const Tree &t, mTrees) {	495	foreach(const Tree &t, mTrees) {
485	resourceUnitSpecies(t.species()).statistics().add(&t, 0);	496	resourceUnitSpecies(t.species()).statistics().add(&t, 0);
486	}	497	}
487		498
488	if (recalculate_stats) {	499	if (recalculate_stats) {
489	for (int i=0;i<mRUSpecies.count();i++) {	500	for (int i=0;i<mRUSpecies.count();i++) {

Subversion Repositories public iLand

(root)/src/core/resourceunit.cpp @ 1222 - Rev 1202 → 1203