WebSVN - public iLand - Diff - Rev 697 and 720


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;

    mSnag = 0;
    mSoil = 0;
}

ResourceUnit::ResourceUnit(const int index)
{
    qDeleteAll(mRUSpecies);
    mSpeciesSet = 0;
    mClimate = 0;
    mPixelCount=0;
    mStockedArea = 0;
    mStockedPixelCount = 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, the calculate a starting value for available n.
    if (mSoil && Model::settings().useDynamicAvailableNitrogen && Model::settings().carbonCycleEnabled) {
        mSoil->setClimateFactor(1.);
        mSoil->calculateYear();
        mUnitVariables.nitrogenAvailable = mSoil->availableNitrogen();
    }
    mHasDeadTrees = false;
    mAverageAging = 0.;

}
void ResourceUnit::setBoundingBox(const QRectF &bb)
{
    mBoundingBox = bb;
    mCornerCoord = 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()];
}

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()
{
    // start simple: just append to the vector...
    mTrees.append(Tree());
    return mTrees.count()-1;
}

/// 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) {
        // nothing to do...
        return;
    }

    // the pixel counters are filled during the height-grid-calculations
    mStockedArea = 100. * mStockedPixelCount; // m2 (1 height grid pixel = 10x10m)

    // 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.)
            qDebug() << "LAI factor > 1";
        );
        (*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) {
        (*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

}

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()
{
    for (int i=0;i<mRUSpecies.count();i++) {
        mRUSpecies[i]->statistics().clear();
    }
    foreach(const Tree &t, mTrees) {
        resourceUnitSpecies(t.species()).statistics().add(&t, 0);
    }

}

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

}

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

Subversion Repositories public iLand

(root)/src/core/resourceunit.cpp @ 1222 - Rev 697 → 720