Subversion Repositories public iLand

/********************************************************************************************

**    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/>.

********************************************************************************************/

#include "sapling.h"

#include "model.h"

#include "species.h"

#include "resourceunit.h"

#include "resourceunitspecies.h"

#include "tree.h"

/** @class Sapling

  @ingroup core

    Sapling stores saplings per species and resource unit and computes sapling growth (before recruitment).

    http://iland.boku.ac.at/sapling+growth+and+competition

    Saplings are established in a separate step (@sa Regeneration). If sapling reach a height of 4m, they are recruited and become "real" iLand-trees.

    Within the regeneration layer, a cohort-approach is applied.

  */

double Sapling::mRecruitmentVariation = 0.1; // +/- 10%

double Sapling::mBrowsingPressure = 0.;

Sapling::Sapling()

{

    mRUS = 0;

    clearStatistics();

    mAdded = 0;

}

// reset statistics, called at newYear

void Sapling::clearStatistics()

{

    // mAdded: removed

    mRecruited=mDied=mLiving=0;

    mSumDbhDied=0.;

    mAvgHeight=0.;

    mAvgAge=0.;

    mAvgDeltaHPot=mAvgHRealized=0.;

}

void Sapling::updateBrowsingPressure()

{

    if (GlobalSettings::instance()->settings().valueBool("model.settings.browsing.enabled"))

        Sapling::mBrowsingPressure = GlobalSettings::instance()->settings().valueDouble("model.settings.browsing.browsingPressure");

    else

        Sapling::mBrowsingPressure = 0.;

}

/// get the *represented* (Reineke's Law) number of trees (N/ha)

double Sapling::livingStemNumber(double &rAvgDbh, double &rAvgHeight, double &rAvgAge) const

{

    double total = 0.;

    double dbh_sum = 0.;

    double h_sum = 0.;

    double age_sum = 0.;

    const SaplingGrowthParameters &p = mRUS->species()->saplingGrowthParameters();

    for (QVector<SaplingTreeOld>::const_iterator it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) {

        float dbh = it->height / p.hdSapling * 100.f;

        if (dbh<1.) // minimum size: 1cm

            continue;

        double n = p.representedStemNumber(dbh); // one cohort on the pixel represents that number of trees

        dbh_sum += n*dbh;

        h_sum += n*it->height;

        age_sum += n*it->age.age;

        total += n;

    }

    if (total>0.) {

        dbh_sum /= total;

        h_sum /= total;

        age_sum /= total;

    }

    rAvgDbh = dbh_sum;

    rAvgHeight = h_sum;

    rAvgAge = age_sum;

    return total;

}

double Sapling::representedStemNumber(float height) const

{

    const SaplingGrowthParameters &p = mRUS->species()->saplingGrowthParameters();

    float dbh = height / p.hdSapling * 100.f;

    double n = p.representedStemNumber(dbh);

    return n;

}

/// maintenance function to clear dead/recruited saplings from storage

void Sapling::cleanupStorage()

{

    QVector<SaplingTreeOld>::iterator forw=mSaplingTrees.begin();

    QVector<SaplingTreeOld>::iterator back;

    // seek last valid

    for (back=mSaplingTrees.end()-1; back>=mSaplingTrees.begin(); --back)

        if ((*back).isValid())

            break;

    if (back<mSaplingTrees.begin()) {

        mSaplingTrees.clear(); // no valid trees available

        return;

    }

    while (forw < back) {

        if (!(*forw).isValid()) {

            *forw = *back; // copy (fill gap)

            while (back>forw) // seek next valid

                if ((*--back).isValid())

                    break;

        }

        ++forw;

    }

    if (back != mSaplingTrees.end()-1) {

        // free resources...

        mSaplingTrees.erase(back+1, mSaplingTrees.end());

    }

}

// not a very good way of checking if sapling is present

// maybe better: use also a (local) maximum sapling height grid

// maybe better: use a bitset:

// position: index of pixel on LIF (absolute index)

bool Sapling::hasSapling(const QPoint &position) const

{

    const QPoint &offset = mRUS->ru()->cornerPointOffset();

    int index = (position.x()- offset.x())*cPxPerRU + (position.y() - offset.y());

    if (index<0)

        qDebug() << "Sapling error";

    return mSapBitset[index];

    /*

    float *target = GlobalSettings::instance()->model()->grid()->ptr(position.x(), position.y());

    QVector<SaplingTree>::const_iterator it;

    for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) {

        if (it->pixel==target)

            return true;

    }

    return false;

    */

}

/// retrieve the height of the sapling at the location 'position' (given in LIF-coordinates)

/// this is quite expensive and only done for initialization

double Sapling::heightAt(const QPoint &position) const

{

    if (!hasSapling(position))

        return 0.;

    // ok, we'll have to search through all saplings

    QVector<SaplingTreeOld>::const_iterator it;

    float *lif_ptr = GlobalSettings::instance()->model()->grid()->ptr(position.x(), position.y());

    for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) {

        if (it->isValid() && it->pixel == lif_ptr)

            return it->height;

    }

    return 0.;

}

void Sapling::setBit(const QPoint &pos_index, bool value)

{

    int index = (pos_index.x() - mRUS->ru()->cornerPointOffset().x()) * cPxPerRU +(pos_index.y() - mRUS->ru()->cornerPointOffset().y());

    mSapBitset.set(index,value); // set bit: now there is a sapling there

}

/// add a sapling at the given position (index on the LIF grid, i.e. 2x2m)

int Sapling::addSapling(const QPoint &pos_lif, const float height, const int age)

{

    // adds a sapling...

    mSaplingTrees.push_back(SaplingTreeOld());

    SaplingTreeOld &t = mSaplingTrees.back();

    t.height = height; // default is 5cm height

    t.age.age = age;

    Grid<float> &lif_map = *GlobalSettings::instance()->model()->grid();

    t.pixel = lif_map.ptr(pos_lif.x(), pos_lif.y());

    setBit(pos_lif, true);

    mAdded++;

    return mSaplingTrees.count()-1; // index of the newly added tree.

}

/// clear  saplings on a given position (after recruitment)

void Sapling::clearSaplings(const QPoint &position)

{

    float *target = GlobalSettings::instance()->model()->grid()->ptr(position.x(), position.y());

    QVector<SaplingTreeOld>::const_iterator it;

    for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) {

        if (it->pixel==target) {

            // trick: use a const iterator to avoid a deep copy of the vector; then do an ugly const_cast to actually write the data

            //const SaplingTree &t = *it;

            //const_cast<SaplingTree&>(t).pixel=0;

            clearSapling(const_cast<SaplingTreeOld&>(*it), false); // kill sapling and move carbon to soil

        }

    }

    setBit(position, false); // clear bit: now there is no sapling on this position

    //int index = (position.x() - mRUS->ru()->cornerPointOffset().x()) * cPxPerRU +(position.y() - mRUS->ru()->cornerPointOffset().y());

    //mSapBitset.set(index,false); // clear bit: now there is no sapling on this position

}

/// clear  saplings within a given rectangle

void Sapling::clearSaplings(const QRectF &rectangle, const bool remove_biomass)

{

    QVector<SaplingTreeOld>::const_iterator it;

    FloatGrid *grid = GlobalSettings::instance()->model()->grid();

    for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) {

        if (rectangle.contains(grid->cellCenterPoint(it->coords()))) {

            clearSapling(const_cast<SaplingTreeOld&>(*it), remove_biomass);

        }

    }

}

void Sapling::clearSapling(SaplingTreeOld &tree, const bool remove)

{

    QPoint p=tree.coords();

    tree.pixel=0;

    setBit(p, false); // no tree left

    if (!remove) {

        // killing of saplings:

        // if remove=false, then remember dbh/number of trees (used later in calculateGrowth() to estimate carbon flow)

        mDied++;

        mSumDbhDied+=tree.height / mRUS->species()->saplingGrowthParameters().hdSapling * 100.;

    }

}

//

void Sapling::clearSapling(int index, const bool remove)

{

    Q_ASSERT(index < mSaplingTrees.count());

    clearSapling(mSaplingTrees[index], remove);

}

/// growth function for an indivudal sapling.

/// returns true, if sapling survives, false if sapling dies or is recruited to iLand.

/// see also http://iland.boku.ac.at/recruitment

bool Sapling::growSapling(SaplingTreeOld &tree, const double f_env_yr, Species* species)

{

    QPoint p=GlobalSettings::instance()->model()->grid()->indexOf(tree.pixel);

    //GlobalSettings::instance()->model()->heightGrid()[Grid::index5(tree.pixel-GlobalSettings::instance()->model()->grid()->begin())];

    // (1) calculate height growth potential for the tree (uses linerization of expressions...)

    double h_pot = species->saplingGrowthParameters().heightGrowthPotential.calculate(tree.height); // TODO check if this can be source of crashes (race condition)

    double delta_h_pot = h_pot - tree.height;

    // (2) reduce height growth potential with species growth response f_env_yr and with light state (i.e. LIF-value) of home-pixel.

    double lif_value = *tree.pixel;

    double h_height_grid = GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(p.x()/cPxPerHeight, p.y()/cPxPerHeight).height;

    if (h_height_grid==0.)

        throw IException(QString("growSapling: height grid at %1/%2 has value 0").arg(p.x()).arg(p.y()));

    double rel_height = tree.height / h_height_grid;

    double lif_corrected = mRUS->species()->speciesSet()->LRIcorrection(lif_value, rel_height); // correction based on height

    double lr = mRUS->species()->lightResponse(lif_corrected); // species specific light response (LUI, light utilization index)

    double delta_h_factor = f_env_yr * lr; // relative growth

    if (h_pot<0. || delta_h_pot<0. || lif_corrected<0. || lif_corrected>1. || delta_h_factor<0. || delta_h_factor>1. )

        qDebug() << "invalid values in Sapling::growSapling";

    // check browsing

    if (mBrowsingPressure>0. && tree.height<=2.f) {

        double p = mRUS->species()->saplingGrowthParameters().browsingProbability;

        // calculate modifed annual browsing probability via odds-ratios

        // odds = p/(1-p) -> odds_mod = odds * browsingPressure -> p_mod = odds_mod /( 1 + odds_mod) === p*pressure/(1-p+p*pressure)

        double p_browse = p*mBrowsingPressure / (1. - p + p*mBrowsingPressure);

        if (drandom() < p_browse) {

            delta_h_factor = 0.;

        }

    }

    // check mortality of saplings

    if (delta_h_factor < species->saplingGrowthParameters().stressThreshold) {

        tree.age.stress_years++;

        if (tree.age.stress_years > species->saplingGrowthParameters().maxStressYears) {

            // sapling dies...

            clearSapling(tree, false); // false: put carbon to the soil

            return false;

        }

    } else {

        tree.age.stress_years=0; // reset stress counter

    }

    DBG_IF(delta_h_pot*delta_h_factor < 0.f || delta_h_pot*delta_h_factor > 2., "Sapling::growSapling", "inplausible height growth.");

    // grow

    tree.height += delta_h_pot * delta_h_factor;

    tree.age.age++; // increase age of sapling by 1

    // recruitment?

    if (tree.height > 4.f) {

        mRecruited++;

        ResourceUnit *ru = const_cast<ResourceUnit*> (mRUS->ru());

        float dbh = tree.height / species->saplingGrowthParameters().hdSapling * 100.f;

        // the number of trees to create (result is in trees per pixel)

        double n_trees = species->saplingGrowthParameters().representedStemNumber(dbh);

        int to_establish = (int) n_trees;

        // if n_trees is not an integer, choose randomly if we should add a tree.

        // e.g.: n_trees = 2.3 -> add 2 trees with 70% probability, and add 3 trees with p=30%.

        if (drandom() < (n_trees-to_establish) || to_establish==0)

            to_establish++;

        // add a new tree

        for (int i=0;i<to_establish;i++) {

            Tree &bigtree = ru->newTree();

            bigtree.setPosition(p);

            // add variation: add +/-10% to dbh and *independently* to height.

            bigtree.setDbh(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation));

            bigtree.setHeight(tree.height * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation));

            bigtree.setSpecies( species );

            bigtree.setAge(tree.age.age,tree.height);

            bigtree.setRU(ru);

            bigtree.setup();

            const Tree *t = &bigtree;

            mRUS->statistics().add(t, 0); // count the newly created trees already in the stats

        }

        // clear all regeneration from this pixel (including this tree)

        clearSapling(tree, true); // remove this tree (but do not move biomass to soil)

//        ru->clearSaplings(p); // remove all other saplings on the same pixel

        return false;

    }

    // book keeping (only for survivors)

    mLiving++;

    mAvgHeight+=tree.height;

    mAvgAge+=tree.age.age;

    mAvgDeltaHPot+=delta_h_pot;

    mAvgHRealized += delta_h_pot * delta_h_factor;

    return true;

}

/** main growth function for saplings.

    Statistics are cleared at the beginning of the year.

*/

void Sapling::calculateGrowth()

{

    Q_ASSERT(mRUS);

    if (mSaplingTrees.count()==0)

        return;

    ResourceUnit *ru = const_cast<ResourceUnit*> (mRUS->ru() );

    Species *species = const_cast<Species*>(mRUS->species());

    // calculate necessary growth modifier (this is done only once per year)

    mRUS->calculate(true); // calculate the 3pg module (this is done only if that did not happen up to now); true: call comes from regeneration

    double f_env_yr = mRUS->prod3PG().fEnvYear();

    mLiving=0;

    QVector<SaplingTreeOld>::const_iterator it;

    for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) {

        const SaplingTreeOld &tree = *it;

        if (tree.height<0)

            qDebug() << "Sapling::calculateGrowth(): h<0";

        // if sapling is still living check execute growth routine

        if (tree.isValid()) {

            // growing (increases mLiving if tree did not die, mDied otherwise)

            if (growSapling(const_cast<SaplingTreeOld&>(tree), f_env_yr, species)) {

                // set the sapling height to the maximum value on the current pixel

//                ru->setMaxSaplingHeightAt(tree.coords(),tree.height);

            }

        }

    }

    if (mLiving) {

        mAvgHeight /= double(mLiving);

        mAvgAge /= double(mLiving);

        mAvgDeltaHPot /= double(mLiving);

        mAvgHRealized /= double(mLiving);

    }

    // calculate carbon balance

    CNPair old_state = mCarbonLiving;

    mCarbonLiving.clear();

    CNPair dead_wood, dead_fine; // pools for mortality

    // average dbh

    if (mLiving) {

        // calculate the avg dbh and number of stems

        double avg_dbh = mAvgHeight / species->saplingGrowthParameters().hdSapling * 100.;

        double n = mLiving * species->saplingGrowthParameters().representedStemNumber( avg_dbh );

        // woody parts: stem, branchse and coarse roots

        double woody_bm = species->biomassWoody(avg_dbh) + species->biomassBranch(avg_dbh) + species->biomassRoot(avg_dbh);

        double foliage = species->biomassFoliage(avg_dbh);

        double fineroot = foliage*species->finerootFoliageRatio();

        mCarbonLiving.addBiomass( woody_bm*n, species->cnWood()  );

        mCarbonLiving.addBiomass( foliage*n, species->cnFoliage()  );

        mCarbonLiving.addBiomass( fineroot*n, species->cnFineroot()  );

        // turnover

        if (mRUS->ru()->snag())

            mRUS->ru()->snag()->addTurnoverLitter(species, foliage*species->turnoverLeaf(), fineroot*species->turnoverRoot());

        // calculate the "mortality from competition", i.e. carbon that stems from reduction of stem numbers

        // from Reinekes formula.

        //

        if (avg_dbh>1.) {

            double avg_dbh_before = (mAvgHeight - mAvgHRealized) / species->saplingGrowthParameters().hdSapling * 100.;

            double n_before = mLiving * species->saplingGrowthParameters().representedStemNumber( qMax(1.,avg_dbh_before) );

            if (n<n_before) {

                dead_wood.addBiomass( woody_bm * (n_before-n), species->cnWood() );

                dead_fine.addBiomass( foliage * (n_before-n), species->cnFoliage()  );

                dead_fine.addBiomass( fineroot * (n_before-n), species->cnFineroot()  );

            }

        }

    }

    if (mDied) {

        double avg_dbh_dead = mSumDbhDied / double(mDied);

        double n = mDied * species->saplingGrowthParameters().representedStemNumber( avg_dbh_dead );

        // woody parts: stem, branchse and coarse roots

        dead_wood.addBiomass( ( species->biomassWoody(avg_dbh_dead) + species->biomassBranch(avg_dbh_dead) + species->biomassRoot(avg_dbh_dead)) * n, species->cnWood()  );

        double foliage = species->biomassFoliage(avg_dbh_dead)*n;

        dead_fine.addBiomass( foliage, species->cnFoliage()  );

        dead_fine.addBiomass( foliage*species->finerootFoliageRatio(), species->cnFineroot()  );

    }

    if (!dead_wood.isEmpty() || !dead_fine.isEmpty())

        if (mRUS->ru()->snag())

            mRUS->ru()->snag()->addToSoil(species, dead_wood, dead_fine);

    // calculate net growth:

    // delta of stocks

    mCarbonGain = mCarbonLiving + dead_fine + dead_wood - old_state;

    if (mCarbonGain.C < 0)

        mCarbonGain.clear();

    if (mSaplingTrees.count() > mLiving*1.3)

        cleanupStorage();

//    mRUS->statistics().add(this);

    GlobalSettings::instance()->systemStatistics()->saplingCount+=mLiving;

    GlobalSettings::instance()->systemStatistics()->newSaplings+=mAdded;

    mAdded = 0; // reset

    //qDebug() << ru->index() << species->id()<< ": (living/avg.height):" <<  mLiving << mAvgHeight;

}

/// fill a grid with the maximum height of saplings per pixel (2x2m).

/// this function is used for visualization only

void Sapling::fillMaxHeightGrid(Grid<float> &grid) const

{

    QVector<SaplingTreeOld>::const_iterator it;

    for (it = mSaplingTrees.begin(); it!=mSaplingTrees.end(); ++it) {

        if (it->isValid()) {

             QPoint p=it->coords();

             if (grid.valueAtIndex(p)<it->height)

                 grid.valueAtIndex(p) = it->height;

        }

    }

}