WebSVN - public iLand - Diff - Rev 639 and 662


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/resourceunit.cpp':
/** @class ResourceUnit
  ResourceUnit is the spatial unit that encapsulates a forest stand and links to several environmental components
  (Climate, Soil, Water, ...).

  */
#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();
    }
    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()
{
    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));
                qDebug() << "reduce tree storage of RU" << index() << " from " << mTrees.capacity() << "to" << mTrees.count();
                mTrees.squeeze();
            }
        }
    }
}

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) {
         (*i)->setLAIfactor((*i)->statistics().leafAreaIndex() / ru_lai);
    }

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

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

Subversion Repositories public iLand

(root)/src/core/resourceunit.cpp @ 1222 - Rev 639 → 662