WebSVN - public iLand - Diff - Rev 1118 and 1157


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);

    mSnag = 0;
    mSoil = 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;
    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)));
    }

    // 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());
}

/// 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();
        (*i)->changeSapling().newYear();
    }

}

/** 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 = 100. * 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 {

 
 
            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
        if (logLevelInfo() &&  (*i)->LAIfactor()>0)
            qDebug() << "ru" << mIndex << "species" << (*i)->species()->id() << "LAIfraction" << (*i)->LAIfactor() << "raw_gpp_m2"
                     << (*i)->prod3PG().GPPperArea() << "area:" << productiveArea() << "gpp:"
                     << productiveArea()*(*i)->prod3PG().GPPperArea()
                     << "aging(lastyear):" << averageAging() << "f_env,yr:" << (*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();
    }

    // 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);
    }
    // summarize statistics for the whole resource unit
    for (int i=0;i<mRUSpecies.count();i++) {
        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)
{
    for (int i=0;i<mRUSpecies.count();i++) {
        mRUSpecies[i]->statistics().clear();
    }
    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::setSaplingHeightMap(float *map_pointer)
{
    // set to zero
    mSaplingHeightMap=map_pointer;
    if (!mSaplingHeightMap)
        return;
    for (int i=0;i<cPxPerRU*cPxPerRU;i++)
        mSaplingHeightMap[i]=0.f;
    // flag all trees in the resource unit
    for (QVector<Tree>::const_iterator i=mTrees.cbegin(); i!= mTrees.cend(); ++i) {
        setMaxSaplingHeightAt(i->positionIndex(), 10.f);
    }
}

void ResourceUnit::setMaxSaplingHeightAt(const QPoint &position, const float height)
{
    Q_ASSERT(mSaplingHeightMap);
    int pixel_index = cPxPerRU*(position.x()-mCornerOffset.x())+(position.y()-mCornerOffset.y());
    if (pixel_index<0 || pixel_index>=cPxPerRU*cPxPerRU) {
        qDebug() << "setSaplingHeightAt-Error for position" << position << "for RU at" << boundingBox() << "with corner" << mCornerOffset;
    } else {
        if (mSaplingHeightMap[pixel_index]<height)
            mSaplingHeightMap[pixel_index]=height;
    }
}

/// clear all saplings of all species on a given position (after recruitment)
void ResourceUnit::clearSaplings(const QPoint &position)
{
    foreach(ResourceUnitSpecies* rus, mRUSpecies)
        rus->clearSaplings(position);

}

/// kill all saplings within a given rect
void ResourceUnit::clearSaplings(const QRectF pixel_rect, const bool remove_from_soil)
{
    foreach(ResourceUnitSpecies* rus, mRUSpecies) {
        rus->changeSapling().clearSaplings(pixel_rect, remove_from_soil);
    }

}

 
double ResourceUnit::saplingHeightForInit(const QPoint &position) const
{
    double maxh = 0.;
    foreach(ResourceUnitSpecies* rus, mRUSpecies)
        maxh = qMax(maxh, rus->sapling().heightAt(position));
    return maxh;
}

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 1118	Rev 1157
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	mSnag = 0;	58	mSnag = 0;
59	mSoil = 0;	59	mSoil = 0;
60	}	60	}
61		61
62	ResourceUnit::ResourceUnit(const int index)	62	ResourceUnit::ResourceUnit(const int index)
63	{	63	{
64	qDeleteAll(mRUSpecies);	64	qDeleteAll(mRUSpecies);
65	mSpeciesSet = 0;	65	mSpeciesSet = 0;
66	mClimate = 0;	66	mClimate = 0;
67	mPixelCount=0;	67	mPixelCount=0;
68	mStockedArea = 0;	68	mStockedArea = 0;
69	mStockedPixelCount = 0;	69	mStockedPixelCount = 0;
-		70	mStockableArea = 0;
70	mAggregatedWLA = 0.;	71	mAggregatedWLA = 0.;
71	mAggregatedLA = 0.;	72	mAggregatedLA = 0.;
72	mAggregatedLR = 0.;	73	mAggregatedLR = 0.;
73	mEffectiveArea = 0.;	74	mEffectiveArea = 0.;
74	mLRI_modification = 0.;	75	mLRI_modification = 0.;
75	mIndex = index;	76	mIndex = index;
76	mSaplingHeightMap = 0;	77	mSaplingHeightMap = 0;
77	mEffectiveArea_perWLA = 0.;	78	mEffectiveArea_perWLA = 0.;
78	mWater = new WaterCycle();	79	mWater = new WaterCycle();
79	mSnag = 0;	80	mSnag = 0;
80	mSoil = 0;	81	mSoil = 0;
81	mID = 0;	82	mID = 0;
82	}	83	}
83		84
84	void ResourceUnit::setup()	85	void ResourceUnit::setup()
85	{	86	{
86	mWater->setup(this);	87	mWater->setup(this);
87		88
88	if (mSnag)	89	if (mSnag)
89	delete mSnag;	90	delete mSnag;
90	mSnag=0;	91	mSnag=0;
91	if (mSoil)	92	if (mSoil)
92	delete mSoil;	93	delete mSoil;
93	mSoil=0;	94	mSoil=0;
94	if (Model::settings().carbonCycleEnabled) {	95	if (Model::settings().carbonCycleEnabled) {
95	mSoil = new Soil(this);	96	mSoil = new Soil(this);
96	mSnag = new Snag;	97	mSnag = new Snag;
97	mSnag->setup(this);	98	mSnag->setup(this);
98	const XmlHelper &xml=GlobalSettings::instance()->settings();	99	const XmlHelper &xml=GlobalSettings::instance()->settings();
99		100
100	// setup contents of the soil of the RU; use values for C and N (kg/ha)	101	// setup contents of the soil of the RU; use values for C and N (kg/ha)
101	mSoil->setInitialState(CNPool(xml.valueDouble("model.site.youngLabileC", -1),	102	mSoil->setInitialState(CNPool(xml.valueDouble("model.site.youngLabileC", -1),
102	xml.valueDouble("model.site.youngLabileN", -1),	103	xml.valueDouble("model.site.youngLabileN", -1),
103	xml.valueDouble("model.site.youngLabileDecompRate", -1)),	104	xml.valueDouble("model.site.youngLabileDecompRate", -1)),
104	CNPool(xml.valueDouble("model.site.youngRefractoryC", -1),	105	CNPool(xml.valueDouble("model.site.youngRefractoryC", -1),
105	xml.valueDouble("model.site.youngRefractoryN", -1),	106	xml.valueDouble("model.site.youngRefractoryN", -1),
106	xml.valueDouble("model.site.youngRefractoryDecompRate", -1)),	107	xml.valueDouble("model.site.youngRefractoryDecompRate", -1)),
107	CNPair(xml.valueDouble("model.site.somC", -1), xml.valueDouble("model.site.somN", -1)));	108	CNPair(xml.valueDouble("model.site.somC", -1), xml.valueDouble("model.site.somN", -1)));
108	}	109	}
109		110
110	// setup variables	111	// setup variables
111	mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40);	112	mUnitVariables.nitrogenAvailable = GlobalSettings::instance()->settings().valueDouble("model.site.availableNitrogen", 40);
112		113
113	// if dynamic coupling of soil nitrogen is enabled, a starting value for available N is calculated	114	// if dynamic coupling of soil nitrogen is enabled, a starting value for available N is calculated
114	if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) {	115	if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) {
115	mSoil->setClimateFactor(1.);	116	mSoil->setClimateFactor(1.);
116	mSoil->calculateYear();	117	mSoil->calculateYear();
117	mUnitVariables.nitrogenAvailable = soil()->availableNitrogen();	118	mUnitVariables.nitrogenAvailable = soil()->availableNitrogen();
118	}	119	}
119	mHasDeadTrees = false;	120	mHasDeadTrees = false;
120	mAverageAging = 0.;	121	mAverageAging = 0.;
121		122
122	}	123	}
123	void ResourceUnit::setBoundingBox(const QRectF &bb)	124	void ResourceUnit::setBoundingBox(const QRectF &bb)
124	{	125	{
125	mBoundingBox = bb;	126	mBoundingBox = bb;
126	mCornerOffset = GlobalSettings::instance()->model()->grid()->indexAt(bb.topLeft());	127	mCornerOffset = GlobalSettings::instance()->model()->grid()->indexAt(bb.topLeft());
127	}	128	}
128		129
129	/// set species and setup the species-per-RU-data	130	/// set species and setup the species-per-RU-data
130	void ResourceUnit::setSpeciesSet(SpeciesSet *set)	131	void ResourceUnit::setSpeciesSet(SpeciesSet *set)
131	{	132	{
132	mSpeciesSet = set;	133	mSpeciesSet = set;
133	qDeleteAll(mRUSpecies);	134	qDeleteAll(mRUSpecies);
134		135
135	//mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated	136	//mRUSpecies.resize(set->count()); // ensure that the vector space is not relocated
136	for (int i=0;i<set->count();i++) {	137	for (int i=0;i<set->count();i++) {
137	Species s = const_cast<Species>(mSpeciesSet->species(i));	138	Species s = const_cast<Species>(mSpeciesSet->species(i));
138	if (!s)	139	if (!s)
139	throw IException("ResourceUnit::setSpeciesSet: invalid index!");	140	throw IException("ResourceUnit::setSpeciesSet: invalid index!");
140		141
141	ResourceUnitSpecies *rus = new ResourceUnitSpecies();	142	ResourceUnitSpecies *rus = new ResourceUnitSpecies();
142	mRUSpecies.push_back(rus);	143	mRUSpecies.push_back(rus);
143	rus->setup(s, this);	144	rus->setup(s, this);
144	/* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.	145	/* be careful: setup() is called with a pointer somewhere to the content of the mRUSpecies container.
145	If the container memory is relocated (QVector), the pointer gets invalid!!!	146	If the container memory is relocated (QVector), the pointer gets invalid!!!
146	Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */	147	Therefore, a resize() is called before the loop (no resize()-operations during the loop)! */
147	//mRUSpecies[i].setup(s,this); // setup this element	148	//mRUSpecies[i].setup(s,this); // setup this element
148		149
149	}	150	}
150	}	151	}
151		152
152	ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)	153	ResourceUnitSpecies &ResourceUnit::resourceUnitSpecies(const Species *species)
153	{	154	{
154	return *mRUSpecies[species->index()];	155	return *mRUSpecies[species->index()];
155	}	156	}
156		157
157	const ResourceUnitSpecies ResourceUnit::constResourceUnitSpecies(const Species species) const	158	const ResourceUnitSpecies ResourceUnit::constResourceUnitSpecies(const Species species) const
158	{	159	{
159	return mRUSpecies[species->index()];	160	return mRUSpecies[species->index()];
160	}	161	}
161		162
162	Tree &ResourceUnit::newTree()	163	Tree &ResourceUnit::newTree()
163	{	164	{
164	// start simple: just append to the vector...	165	// start simple: just append to the vector...
165	if (mTrees.isEmpty())	166	if (mTrees.isEmpty())
166	mTrees.reserve(100); // reserve a junk of memory for trees	167	mTrees.reserve(100); // reserve a junk of memory for trees
167		168
168	mTrees.append(Tree());	169	mTrees.append(Tree());
169	return mTrees.back();	170	return mTrees.back();
170	}	171	}
171	int ResourceUnit::newTreeIndex()	172	int ResourceUnit::newTreeIndex()
172	{	173	{
173	newTree();	174	newTree();
174	return mTrees.count()-1; // return index of the last tree	175	return mTrees.count()-1; // return index of the last tree
175	}	176	}
176		177
177	/// remove dead trees from tree list	178	/// remove dead trees from tree list
178	/// reduce size of vector if lots of space is free	179	/// reduce size of vector if lots of space is free
179	/// tests showed that this way of cleanup is very fast,	180	/// tests showed that this way of cleanup is very fast,
180	/// because no memory allocations are performed (simple memmove())	181	/// because no memory allocations are performed (simple memmove())
181	/// when trees are moved.	182	/// when trees are moved.
182	void ResourceUnit::cleanTreeList()	183	void ResourceUnit::cleanTreeList()
183	{	184	{
184	if (!mHasDeadTrees)	185	if (!mHasDeadTrees)
185	return;	186	return;
186		187
187	QVector<Tree>::iterator last=mTrees.end()-1;	188	QVector<Tree>::iterator last=mTrees.end()-1;
188	QVector<Tree>::iterator current = mTrees.begin();	189	QVector<Tree>::iterator current = mTrees.begin();
189	while (last>=current && (*last).isDead())	190	while (last>=current && (*last).isDead())
190	--last;	191	--last;
191		192
192	while (current<last) {	193	while (current<last) {
193	if ((*current).isDead()) {	194	if ((*current).isDead()) {
194	current = last; // copy data!	195	current = last; // copy data!
195	--last; //	196	--last; //
196	while (last>=current && (*last).isDead())	197	while (last>=current && (*last).isDead())
197	--last;	198	--last;
198	}	199	}
199	++current;	200	++current;
200	}	201	}
201	++last; // last points now to the first dead tree	202	++last; // last points now to the first dead tree
202		203
203	// free ressources	204	// free ressources
204	if (last!=mTrees.end()) {	205	if (last!=mTrees.end()) {
205	mTrees.erase(last, mTrees.end());	206	mTrees.erase(last, mTrees.end());
206	if (mTrees.capacity()>100) {	207	if (mTrees.capacity()>100) {
207	if (mTrees.count() / double(mTrees.capacity()) < 0.2) {	208	if (mTrees.count() / double(mTrees.capacity()) < 0.2) {
208	//int target_size = mTrees.count()*2;	209	//int target_size = mTrees.count()*2;
209	//qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;	210	//qDebug() << "reduce size from "<<mTrees.capacity() << "to" << target_size;
210	//mTrees.reserve(qMax(target_size, 100));	211	//mTrees.reserve(qMax(target_size, 100));
211	if (logLevelDebug())	212	if (logLevelDebug())
212	qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count();	213	qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count();
213	mTrees.squeeze();	214	mTrees.squeeze();
214	}	215	}
215	}	216	}
216	}	217	}
217	mHasDeadTrees = false; // reset flag	218	mHasDeadTrees = false; // reset flag
218	}	219	}
219		220
220	void ResourceUnit::newYear()	221	void ResourceUnit::newYear()
221	{	222	{
222	mAggregatedWLA = 0.;	223	mAggregatedWLA = 0.;
223	mAggregatedLA = 0.;	224	mAggregatedLA = 0.;
224	mAggregatedLR = 0.;	225	mAggregatedLR = 0.;
225	mEffectiveArea = 0.;	226	mEffectiveArea = 0.;
226	mPixelCount = mStockedPixelCount = 0;	227	mPixelCount = mStockedPixelCount = 0;
227	snagNewYear();	228	snagNewYear();
228	if (mSoil)	229	if (mSoil)
229	mSoil->newYear();	230	mSoil->newYear();
230	// clear statistics global and per species...	231	// clear statistics global and per species...
231	QList<ResourceUnitSpecies*>::const_iterator i;	232	QList<ResourceUnitSpecies*>::const_iterator i;
232	QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();	233	QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();
233	mStatistics.clear();	234	mStatistics.clear();
234	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {	235	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
235	(*i)->statisticsDead().clear();	236	(*i)->statisticsDead().clear();
236	(*i)->statisticsMgmt().clear();	237	(*i)->statisticsMgmt().clear();
237	(*i)->changeSapling().newYear();	238	(*i)->changeSapling().newYear();
238	}	239	}
239		240
240	}	241	}
241		242
242	/** production() is the "stand-level" part of the biomass production (3PG).	243	/** production() is the "stand-level" part of the biomass production (3PG).
243	- The amount of radiation intercepted by the stand is calculated	244	- The amount of radiation intercepted by the stand is calculated
244	- the water cycle is calculated	245	- the water cycle is calculated
245	- statistics for each species are cleared	246	- statistics for each species are cleared
246	- The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production)	247	- The 3PG production for each species and ressource unit is called (calculates species-responses and NPP production)
247	see also: http://iland.boku.ac.at/individual+tree+light+availability */	248	see also: http://iland.boku.ac.at/individual+tree+light+availability */
248	void ResourceUnit::production()	249	void ResourceUnit::production()
249	{	250	{
250		251
251	if (mAggregatedWLA==0. \|\| mPixelCount==0) {	252	if (mAggregatedWLA==0. \|\| mPixelCount==0) {
252	// clear statistics of resourceunitspecies	253	// clear statistics of resourceunitspecies
253	for ( QList<ResourceUnitSpecies*>::const_iterator i=mRUSpecies.constBegin(); i!=mRUSpecies.constEnd(); ++i)	254	for ( QList<ResourceUnitSpecies*>::const_iterator i=mRUSpecies.constBegin(); i!=mRUSpecies.constEnd(); ++i)
254	(*i)->statistics().clear();	255	(*i)->statistics().clear();
255	mEffectiveArea = 0.;	256	mEffectiveArea = 0.;
256	mStockedArea = 0.;	257	mStockedArea = 0.;
257	return;	258	return;
258	}	259	}
259		260
260	// the pixel counters are filled during the height-grid-calculations	261	// the pixel counters are filled during the height-grid-calculations
261	mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)	262	mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)
262	if (leafAreaIndex()<3.) {	263	if (leafAreaIndex()<3.) {
263	// estimate stocked area based on crown projections	264	// estimate stocked area based on crown projections
264	double crown_area = 0.;	265	double crown_area = 0.;
265	for (int i=0;i<mTrees.count();++i)	266	for (int i=0;i<mTrees.count();++i)
266	crown_area += mTrees.at(i).isDead() ? 0. : mTrees.at(i).stamp()->reader()->crownArea();	267	crown_area += mTrees.at(i).isDead() ? 0. : mTrees.at(i).stamp()->reader()->crownArea();
267		268
268	qDebug() << "crown area: lai" << leafAreaIndex() << "stocked area (pixels)" << mStockedArea << " area (crown)" << crown_area;	-
269	if (leafAreaIndex()<2.) {	-
270	mStockedArea = crown_area;	-
-		269	if (logLevelDebug())
-		270	qDebug() << "crown area: lai" << leafAreaIndex() << "stocked area (pixels)" << mStockedArea << " area (crown)" << crown_area;
-		271	if (leafAreaIndex()<1.) {
-		272	mStockedArea = std::min(crown_area, mStockedArea);
271	} else {	273	} else {
272		274
273	double px_frac = leafAreaIndex()-2.; // 0 at LAI=2, 1 at LAI=3	-
274	mStockedArea = mStockedArea * px_frac + crown_area * (1. - px_frac);	-
-		275	double px_frac = (leafAreaIndex()-1.)/2.; // 0 at LAI=1, 1 at LAI=3
-		276	mStockedArea = mStockedArea * px_frac + std::min(crown_area, mStockedArea) * (1. - px_frac);
275	}	277	}
276	if (mStockedArea==0.)	278	if (mStockedArea==0.)
277	return;	279	return;
278	}	280	}
279		281
280	// calculate the leaf area index (LAI)	282	// calculate the leaf area index (LAI)
281	double LAI = mAggregatedLA / mStockedArea;	283	double LAI = mAggregatedLA / mStockedArea;
282	// calculate the intercepted radiation fraction using the law of Beer Lambert	284	// calculate the intercepted radiation fraction using the law of Beer Lambert
283	const double k = Model::settings().lightExtinctionCoefficient;	285	const double k = Model::settings().lightExtinctionCoefficient;
284	double interception_fraction = 1. - exp(-k * LAI);	286	double interception_fraction = 1. - exp(-k * LAI);
285	mEffectiveArea = mStockedArea * interception_fraction; // m2	287	mEffectiveArea = mStockedArea * interception_fraction; // m2
286		288
287	// calculate the total weighted leaf area on this RU:	289	// calculate the total weighted leaf area on this RU:
288	mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA; // p_WLA	290	mLRI_modification = interception_fraction * mStockedArea / mAggregatedWLA; // p_WLA
289	if (mLRI_modification == 0.)	291	if (mLRI_modification == 0.)
290	qDebug() << "lri modifaction==0!";	292	qDebug() << "lri modifaction==0!";
291		293
292	if (logLevelDebug()) {	294	if (logLevelDebug()) {
293	DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")	295	DBGMODE(qDebug() << QString("production: LAI: %1 (intercepted fraction: %2, stocked area: %4). LRI-Multiplier: %3")
294	.arg(LAI)	296	.arg(LAI)
295	.arg(interception_fraction)	297	.arg(interception_fraction)
296	.arg(mLRI_modification)	298	.arg(mLRI_modification)
297	.arg(mStockedArea);	299	.arg(mStockedArea);
298	);	300	);
299	}	301	}
300		302
301	// calculate LAI fractions	303	// calculate LAI fractions
302	QList<ResourceUnitSpecies*>::const_iterator i;	304	QList<ResourceUnitSpecies*>::const_iterator i;
303	QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();	305	QList<ResourceUnitSpecies*>::const_iterator iend = mRUSpecies.constEnd();
304	double ru_lai = leafAreaIndex();	306	double ru_lai = leafAreaIndex();
305	if (ru_lai < 1.)	307	if (ru_lai < 1.)
306	ru_lai = 1.;	308	ru_lai = 1.;
307	// note: LAIFactors are only 1 if sum of LAI is > 1. (see WaterCycle)	309	// note: LAIFactors are only 1 if sum of LAI is > 1. (see WaterCycle)
308	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {	310	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
309	double lai_factor = (*i)->statistics().leafAreaIndex() / ru_lai;	311	double lai_factor = (*i)->statistics().leafAreaIndex() / ru_lai;
310	DBGMODE(	-
311	if (lai_factor > 1.)	-
312	qDebug() << "LAI factor > 1";	-
313	);	-
-		312
-		313	//DBGMODE(
-		314	if (lai_factor > 1.) {
-		315	const ResourceUnitSpecies* rus=*i;
-		316	qDebug() << "LAI factor > 1: species ru-index:" << rus->species()->name() << rus->ru()->index();
-		317	}
-		318	//);
314	(*i)->setLAIfactor( lai_factor );	319	(*i)->setLAIfactor( lai_factor );
315	}	320	}
316		321
317	// soil water model - this determines soil water contents needed for response calculations	322	// soil water model - this determines soil water contents needed for response calculations
318	{	323	{
319	mWater->run();	324	mWater->run();
320	}	325	}
321		326
322	// invoke species specific calculation (3PG)	327	// invoke species specific calculation (3PG)
323	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {	328	for (i=mRUSpecies.constBegin(); i!=iend; ++i) {
-		329	//DBGMODE(
-		330	if ((*i)->LAIfactor() > 1.) {
-		331	const ResourceUnitSpecies* rus=*i;
-		332	qDebug() << "LAI factor > 1: species ru-index value:" << rus->species()->name() << rus->ru()->index() << rus->LAIfactor();
-		333	}
-		334	//);
324	(*i)->calculate(); // CALCULATE 3PG	335	(*i)->calculate(); // CALCULATE 3PG
325	if (logLevelInfo() && (*i)->LAIfactor()>0)	336	if (logLevelInfo() && (*i)->LAIfactor()>0)
326	qDebug() << "ru" << mIndex << "species" << (i)->species()->id() << "LAIfraction" << (i)->LAIfactor() << "raw_gpp_m2"	337	qDebug() << "ru" << mIndex << "species" << (i)->species()->id() << "LAIfraction" << (i)->LAIfactor() << "raw_gpp_m2"
327	<< (*i)->prod3PG().GPPperArea() << "area:" << productiveArea() << "gpp:"	338	<< (*i)->prod3PG().GPPperArea() << "area:" << productiveArea() << "gpp:"
328	<< productiveArea()(i)->prod3PG().GPPperArea()	339	<< productiveArea()(i)->prod3PG().GPPperArea()
329	<< "aging(lastyear):" << averageAging() << "f_env,yr:" << (*i)->prod3PG().fEnvYear();	340	<< "aging(lastyear):" << averageAging() << "f_env,yr:" << (*i)->prod3PG().fEnvYear();
330	}	341	}
331	}	342	}
332		343
333	void ResourceUnit::calculateInterceptedArea()	344	void ResourceUnit::calculateInterceptedArea()
334	{	345	{
335	if (mAggregatedLR==0) {	346	if (mAggregatedLR==0) {
336	mEffectiveArea_perWLA = 0.;	347	mEffectiveArea_perWLA = 0.;
337	return;	348	return;
338	}	349	}
339	Q_ASSERT(mAggregatedLR>0.);	350	Q_ASSERT(mAggregatedLR>0.);
340	mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;	351	mEffectiveArea_perWLA = mEffectiveArea / mAggregatedLR;
341	if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA;	352	if (logLevelDebug()) qDebug() << "RU: aggregated lightresponse:" << mAggregatedLR << "eff.area./wla:" << mEffectiveArea_perWLA;
342	}	353	}
343		354
344	// function is called immediately before the growth of individuals	355	// function is called immediately before the growth of individuals
345	void ResourceUnit::beforeGrow()	356	void ResourceUnit::beforeGrow()
346	{	357	{
347	mAverageAging = 0.;	358	mAverageAging = 0.;
348	}	359	}
349		360
350	// function is called after finishing the indivdual growth / mortality.	361	// function is called after finishing the indivdual growth / mortality.
351	void ResourceUnit::afterGrow()	362	void ResourceUnit::afterGrow()
352	{	363	{
353	mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees)	364	mAverageAging = leafArea()>0.?mAverageAging/leafArea():0; // calculate aging value (calls to addAverageAging() by individual trees)
354	if (mAverageAging>0. && mAverageAging<0.00001)	365	if (mAverageAging>0. && mAverageAging<0.00001)
355	qDebug() << "ru" << mIndex << "aging <0.00001";	366	qDebug() << "ru" << mIndex << "aging <0.00001";
356	if (mAverageAging<0. \|\| mAverageAging>1.)	367	if (mAverageAging<0. \|\| mAverageAging>1.)
357	qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();	368	qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();
358	}	369	}
359		370
360	void ResourceUnit::yearEnd()	371	void ResourceUnit::yearEnd()
361	{	372	{
362	// calculate statistics for all tree species of the ressource unit	373	// calculate statistics for all tree species of the ressource unit
363	int c = mRUSpecies.count();	374	int c = mRUSpecies.count();
364	for (int i=0;i<c; i++) {	375	for (int i=0;i<c; i++) {
365	mRUSpecies[i]->statisticsDead().calculate(); // calculate the dead trees	376	mRUSpecies[i]->statisticsDead().calculate(); // calculate the dead trees
366	mRUSpecies[i]->statisticsMgmt().calculate(); // stats of removed trees	377	mRUSpecies[i]->statisticsMgmt().calculate(); // stats of removed trees
367	mRUSpecies[i]->updateGWL(); // get sum of dead trees (died + removed)	378	mRUSpecies[i]->updateGWL(); // get sum of dead trees (died + removed)
368	mRUSpecies[i]->statistics().calculate(); // calculate the living (and add removed volume to gwl)	379	mRUSpecies[i]->statistics().calculate(); // calculate the living (and add removed volume to gwl)
369	mStatistics.add(mRUSpecies[i]->statistics());	380	mStatistics.add(mRUSpecies[i]->statistics());
370	}	381	}
371	mStatistics.calculate(); // aggreagte on stand level	382	mStatistics.calculate(); // aggreagte on stand level
-		383
-		384	// update carbon flows
-		385	if (soil() && GlobalSettings::instance()->model()->settings().carbonCycleEnabled) {
-		386	double area_factor = stockableArea() / cRUArea; //conversion factor
-		387	mUnitVariables.carbonUptake = statistics().npp() * biomassCFraction;
-		388	mUnitVariables.carbonUptake += statistics().nppSaplings() * biomassCFraction;
-		389
-		390	double to_atm = snag()->fluxToAtmosphere().C / area_factor; // from snags, kgC/ha
-		391	to_atm += soil()->fluxToAtmosphere().C *cRUArea/10.; // soil: t/ha -> t/m2 -> kg/ha
-		392	mUnitVariables.carbonToAtm = to_atm;
-		393
-		394	double to_dist = snag()->fluxToDisturbance().C / area_factor;
-		395	to_dist += soil()->fluxToDisturbance().C * cRUArea/10.;
-		396	double to_harvest = snag()->fluxToExtern().C / area_factor;
-		397
-		398	mUnitVariables.NEP = mUnitVariables.carbonUptake - to_atm - to_dist - to_harvest; // kgC/ha
-		399
-		400	// incremental values....
-		401	mUnitVariables.cumCarbonUptake += mUnitVariables.carbonUptake;
-		402	mUnitVariables.cumCarbonToAtm += mUnitVariables.carbonToAtm;
-		403	mUnitVariables.cumNEP += mUnitVariables.NEP;
-		404
-		405	}
372		406
373	}	407	}
374		408
375	void ResourceUnit::addTreeAgingForAllTrees()	409	void ResourceUnit::addTreeAgingForAllTrees()
376	{	410	{
377	mAverageAging = 0.;	411	mAverageAging = 0.;
378	foreach(const Tree &t, mTrees) {	412	foreach(const Tree &t, mTrees) {
379	addTreeAging(t.leafArea(), t.species()->aging(t.height(), t.age()));	413	addTreeAging(t.leafArea(), t.species()->aging(t.height(), t.age()));
380	}	414	}
381		415
382	}	416	}
383		417
384	/// refresh of tree based statistics.	418	/// refresh of tree based statistics.
385	/// WARNING: this function is only called once (during startup).	419	/// WARNING: this function is only called once (during startup).
386	/// see function "yearEnd()" above!!!	420	/// see function "yearEnd()" above!!!
387	void ResourceUnit::createStandStatistics()	421	void ResourceUnit::createStandStatistics()
388	{	422	{
389	// clear statistics (ru-level and ru-species level)	423	// clear statistics (ru-level and ru-species level)
390	mStatistics.clear();	424	mStatistics.clear();
391	for (int i=0;i<mRUSpecies.count();i++) {	425	for (int i=0;i<mRUSpecies.count();i++) {
392	mRUSpecies[i]->statistics().clear();	426	mRUSpecies[i]->statistics().clear();
393	mRUSpecies[i]->statisticsDead().clear();	427	mRUSpecies[i]->statisticsDead().clear();
394	mRUSpecies[i]->statisticsMgmt().clear();	428	mRUSpecies[i]->statisticsMgmt().clear();
395	}	429	}
396		430
397	// add all trees to the statistics objects of the species	431	// add all trees to the statistics objects of the species
398	foreach(const Tree &t, mTrees) {	432	foreach(const Tree &t, mTrees) {
399	if (!t.isDead())	433	if (!t.isDead())
400	resourceUnitSpecies(t.species()).statistics().add(&t, 0);	434	resourceUnitSpecies(t.species()).statistics().add(&t, 0);
401	}	435	}
402	// summarize statistics for the whole resource unit	436	// summarize statistics for the whole resource unit
403	for (int i=0;i<mRUSpecies.count();i++) {	437	for (int i=0;i<mRUSpecies.count();i++) {
404	mRUSpecies[i]->statistics().calculate();	438	mRUSpecies[i]->statistics().calculate();
405	mStatistics.add(mRUSpecies[i]->statistics());	439	mStatistics.add(mRUSpecies[i]->statistics());
406	}	440	}
407	mStatistics.calculate();	441	mStatistics.calculate();
408	mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*stockableArea()):0.;	442	mAverageAging = mStatistics.leafAreaIndex()>0.?mAverageAging / (mStatistics.leafAreaIndex()*stockableArea()):0.;
409	if (mAverageAging<0. \|\| mAverageAging>1.)	443	if (mAverageAging<0. \|\| mAverageAging>1.)
410	qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();	444	qDebug() << "Average aging invalid: (RU, LAI):" << index() << mStatistics.leafAreaIndex();
411	}	445	}
412		446
413	/** recreate statistics. This is necessary after events that changed the structure	447	/** recreate statistics. This is necessary after events that changed the structure
414	of the stand after the growth of trees (where stand statistics are updated).	448	of the stand after the growth of trees (where stand statistics are updated).
415	An example is after disturbances. */	449	An example is after disturbances. */
416	void ResourceUnit::recreateStandStatistics()	-
-		450	void ResourceUnit::recreateStandStatistics(bool recalculate_stats)
417	{	451	{
418	for (int i=0;i<mRUSpecies.count();i++) {	452	for (int i=0;i<mRUSpecies.count();i++) {
419	mRUSpecies[i]->statistics().clear();	453	mRUSpecies[i]->statistics().clear();
420	}	454	}
421	foreach(const Tree &t, mTrees) {	455	foreach(const Tree &t, mTrees) {
422	resourceUnitSpecies(t.species()).statistics().add(&t, 0);	456	resourceUnitSpecies(t.species()).statistics().add(&t, 0);
423	}	457	}
424	for (int i=0;i<mRUSpecies.count();i++) {	-
425	mRUSpecies[i]->statistics().calculate();	-
-		458
-		459	if (recalculate_stats) {
-		460	for (int i=0;i<mRUSpecies.count();i++) {
-		461	mRUSpecies[i]->statistics().calculate();
-		462	}
426	}	463	}
427	}	464	}
428		465
429	void ResourceUnit::setSaplingHeightMap(float *map_pointer)	466	void ResourceUnit::setSaplingHeightMap(float *map_pointer)
430	{	467	{
431	// set to zero	468	// set to zero
432	mSaplingHeightMap=map_pointer;	469	mSaplingHeightMap=map_pointer;
433	if (!mSaplingHeightMap)	470	if (!mSaplingHeightMap)
434	return;	471	return;
435	for (int i=0;i<cPxPerRU*cPxPerRU;i++)	472	for (int i=0;i<cPxPerRU*cPxPerRU;i++)
436	mSaplingHeightMap[i]=0.f;	473	mSaplingHeightMap[i]=0.f;
437	// flag all trees in the resource unit	474	// flag all trees in the resource unit
438	for (QVector<Tree>::const_iterator i=mTrees.cbegin(); i!= mTrees.cend(); ++i) {	475	for (QVector<Tree>::const_iterator i=mTrees.cbegin(); i!= mTrees.cend(); ++i) {
439	setMaxSaplingHeightAt(i->positionIndex(), 10.f);	476	setMaxSaplingHeightAt(i->positionIndex(), 10.f);
440	}	477	}
441	}	478	}
442		479
443	void ResourceUnit::setMaxSaplingHeightAt(const QPoint &position, const float height)	480	void ResourceUnit::setMaxSaplingHeightAt(const QPoint &position, const float height)
444	{	481	{
445	Q_ASSERT(mSaplingHeightMap);	482	Q_ASSERT(mSaplingHeightMap);
446	int pixel_index = cPxPerRU*(position.x()-mCornerOffset.x())+(position.y()-mCornerOffset.y());	483	int pixel_index = cPxPerRU*(position.x()-mCornerOffset.x())+(position.y()-mCornerOffset.y());
447	if (pixel_index<0 \|\| pixel_index>=cPxPerRU*cPxPerRU) {	484	if (pixel_index<0 \|\| pixel_index>=cPxPerRU*cPxPerRU) {
448	qDebug() << "setSaplingHeightAt-Error for position" << position << "for RU at" << boundingBox() << "with corner" << mCornerOffset;	485	qDebug() << "setSaplingHeightAt-Error for position" << position << "for RU at" << boundingBox() << "with corner" << mCornerOffset;
449	} else {	486	} else {
450	if (mSaplingHeightMap[pixel_index]<height)	487	if (mSaplingHeightMap[pixel_index]<height)
451	mSaplingHeightMap[pixel_index]=height;	488	mSaplingHeightMap[pixel_index]=height;

Subversion Repositories public iLand

(root)/src/core/resourceunit.cpp @ 1222 - Rev 1118 → 1157