WebSVN - public iLand - Diff - Rev 1164 and 1165


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/saplings.cpp':
#include "global.h"
#include "saplings.h"

#include "globalsettings.h"
#include "model.h"
#include "resourceunit.h"
#include "resourceunitspecies.h"
#include "establishment.h"
#include "species.h"
#include "seeddispersal.h"

double Saplings::mRecruitmentVariation = 0.1; // +/- 10%
double Saplings::mBrowsingPressure = 0.;


Saplings::Saplings()
{

}

void Saplings::setup()
{
    //mGrid.setup(GlobalSettings::instance()->model()->grid()->metricRect(), GlobalSettings::instance()->model()->grid()->cellsize());
    FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();
    // mask out out-of-project areas
    HeightGrid *hg = GlobalSettings::instance()->model()->heightGrid();
    for (int i=0; i<lif_grid->count(); ++i) {
        SaplingCell *s = cell(lif_grid->indexOf(i), false); // false: retrieve also invalid cells
        if (s) {
            if (!hg->valueAtIndex(lif_grid->index5(i)).isValid())
                s->state = SaplingCell::CellInvalid;
            else
                s->state = SaplingCell::CellFree;
        }

    }

}

void Saplings::establishment(const ResourceUnit *ru)
{
    HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid();
    FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();

    QPoint imap = ru->cornerPointOffset(); // offset on LIF/saplings grid
    QPoint iseedmap = QPoint(imap.x()/10, imap.y()/10); // seed-map has 20m resolution, LIF 2m -> factor 10

    for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i)
        (*i)->saplingStat().clearStatistics();

    double lif_corr[cPxPerHectare];
    for (int i=0;i<cPxPerHectare;++i)
        lif_corr[i]=-1.;

    int species_idx;
    QVector<int>::const_iterator sbegin, send;
    ru->speciesSet()->randomSpeciesOrder(sbegin, send);
    for (QVector<int>::const_iterator s_idx=sbegin; s_idx!=send;++s_idx) {
 

        // start from a random species (and cycle through the available species)
        species_idx = *s_idx;

        ResourceUnitSpecies *rus = ru->ruSpecies()[species_idx];
        // check if there are seeds of the given species on the resource unit
        float seeds = 0.f;
        Grid<float> &seedmap =  const_cast<Grid<float>& >(rus->species()->seedDispersal()->seedMap());
        for (int iy=0;iy<5;++iy) {
            float *p = seedmap.ptr(iseedmap.x(), iseedmap.y());
            for (int ix=0;ix<5;++ix)
                seeds += *p++;
        }
        // if there are no seeds: no need to do more
        if (seeds==0.f)
            continue;

        // calculate the abiotic environment (TACA)
        rus->establishment().calculateAbioticEnvironment();
        double abiotic_env = rus->establishment().abioticEnvironment();
        if (abiotic_env==0.)
            continue;

        // loop over all 2m cells on this resource unit
        SaplingCell *sap_cells = ru->saplingCellArray();
        SaplingCell *s;
        int isc = 0; // index on 2m cell
        for (int iy=0; iy<cPxPerRU; ++iy) {
            //s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row
            s = &sap_cells[iy*cPxPerRU]; // pointer to a row
            isc = lif_grid->index(imap.x(), imap.y()+iy);

            for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) {
                if (s->state == SaplingCell::CellFree) {
                    // is a sapling of the current species already on the pixel?
                    // * test for sapling height already in cell state
                    // * test for grass-cover already in cell state
                    SaplingTree *stree=0;
                    SaplingTree *slot=s->saplings;
                    for (int i=0;i<NSAPCELLS;++i, ++slot) {
                        if (!stree && !slot->is_occupied())
                            stree=slot;
                        if (slot->species_index == species_idx) {
                            stree=0;
                            break;
                        }
                    }

                    if (stree) {
                        // grass cover?
                        float seed_map_value = seedmap[lif_grid->index10(isc)];
                        if (seed_map_value==0.f)
                            continue;
                        const HeightGridValue &hgv = (*height_grid)[lif_grid->index5(isc)];
                        float lif_value = (*lif_grid)[isc];

                        double &lif_corrected = lif_corr[iy*cPxPerRU+ix];
                        // calculate the LIFcorrected only once per pixel
                        if (lif_corrected<0.)
                            lif_corrected = rus->species()->speciesSet()->LRIcorrection(lif_value, 4. / hgv.height);

                        // check for the combination of seed availability and light on the forest floor
                        if (drandom() < seed_map_value*lif_corrected*abiotic_env ) {
                            // ok, lets add a sapling at the given position (age is incremented later)
                            stree->setSapling(0.05f, 0, species_idx);
                            s->checkState();
                            rus->saplingStat().mAdded++;

                        }

                    }

                }
            }
        }

    }

}

void Saplings::saplingGrowth(const ResourceUnit *ru)
{
    HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid();
    FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();

    QPoint imap = ru->cornerPointOffset();
    bool need_check=false;
    SaplingCell *sap_cells = ru->saplingCellArray();
    for (int iy=0; iy<cPxPerRU; ++iy) {
        //SaplingCell *s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row
        SaplingCell *s = &sap_cells[iy*cPxPerRU]; // ptr to row
        int isc = lif_grid->index(imap.x(), imap.y()+iy);

        for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) {
            if (s->state != SaplingCell::CellInvalid) {
                need_check=false;
                for (int i=0;i<NSAPCELLS;++i) {
                    if (s->saplings[i].is_occupied()) {
                        // growth of this sapling tree
                        const HeightGridValue &hgv = (*height_grid)[height_grid->index5(isc)];
                        float lif_value = (*lif_grid)[isc];

                        need_check |= growSapling(ru, *s, s->saplings[i], isc, hgv.height, lif_value);
                    }
                }
                if (need_check)
                    s->checkState();
            }
        }
    }




    // store statistics on saplings/regeneration
    for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i) {
        (*i)->saplingStat().calculate((*i)->species(), const_cast<ResourceUnit*>(ru));
        (*i)->statistics().add(&((*i)->saplingStat()));
    }
}

SaplingCell *Saplings::cell(QPoint lif_coords, bool only_valid, ResourceUnit **rRUPtr)
{
    FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();

    // in this case, getting the actual cell is quite cumbersome: first, retrieve the resource unit, then the
    // cell based on the offset of the given coordiantes relative to the corner of the resource unit.
    ResourceUnit *ru = GlobalSettings::instance()->model()->ru(lif_grid->cellCenterPoint(lif_coords));
    if (rRUPtr)
        *rRUPtr = ru;

    if (ru) {
        QPoint local_coords = lif_coords - ru->cornerPointOffset();
        int idx = local_coords.y() * cPxPerRU + local_coords.x();
        DBGMODE( if (idx<0 || idx>=cPxPerHectare)
                 qDebug("invalid coords in Saplings::cell");
                    );
        SaplingCell *s=&ru->saplingCellArray()[idx];
        if (s && (!only_valid || s->state!=SaplingCell::CellInvalid))
            return s;
    }
    return 0;
}

void Saplings::clearSaplings(const QRectF &rectangle, const bool remove_biomass)
{
    GridRunner<float> runner(GlobalSettings::instance()->model()->grid(), rectangle);
    ResourceUnit *ru;
    while (runner.next()) {
        SaplingCell *s = cell(runner.currentIndex(), true, &ru);
        if (s) {
            clearSaplings(s, ru, remove_biomass);
        }

 
 
 
 
 
 
 
 
 
    }
}

void Saplings::clearSaplings(SaplingCell *s, ResourceUnit *ru, const bool remove_biomass)
{
    if (s) {
        for (int i=0;i<NSAPCELLS;++i)
            if (s->saplings[i].is_occupied()) {
                if (!remove_biomass) {
                    ResourceUnitSpecies *rus = ru->resourceUnitSpecies(s->saplings[i].species_index);
                    if (!rus && !rus->species()) {
                        qDebug() << "Saplings::clearSaplings(): invalid resource unit!!!";
                        return;
                    }
 
                    rus->saplingStat().addCarbonOfDeadSapling( s->saplings[i].height / rus->species()->saplingGrowthParameters().hdSapling * 100.f );
                }
 
                s->saplings[i].clear();
            }
        s->checkState();

    }

}

int Saplings::addSprout(const Tree *t)
{
    if (t->species()->saplingGrowthParameters().sproutGrowth==0.)
        return 0;
    SaplingCell *sc = cell(t->positionIndex());
    if (!sc)
        return 0;
    clearSaplings(sc, const_cast<ResourceUnit*>(t->ru()), false );
    SaplingTree *st=sc->addSapling(0.05f, 0, t->species()->index());
    if (st)
        st->set_sprout(true);

    // neighboring cells
    double crown_area = t->crownRadius()*t->crownRadius() * M_PI; //m2
    // calculate how many cells on the ground are covered by the crown (this is a rather rough estimate)
    // n_cells: in addition to the original cell
    int n_cells = static_cast<int>(round( crown_area / static_cast<double>(cPxSize*cPxSize) - 1.));
    if (n_cells>0) {
        ResourceUnit *ru;
        static const int offsets_x[8] = {1,1,0,-1,-1,-1,0,1};
        static const int offsets_y[8] = {0,1,1,1,0,-1,-1,-1};
        int s=irandom(0,8);
        while(n_cells) {
            sc = cell(t->positionIndex()+QPoint(offsets_x[s], offsets_y[s]),true,&ru);
            if (sc) {
                clearSaplings(sc, ru, false );
                SaplingTree *st=sc->addSapling(0.05f, 0, t->species()->index());
                if (st)
                    st->set_sprout(true);
            }

            s = (s+1)%8; --n_cells;
        }
    }
    return 1;
}

void Saplings::updateBrowsingPressure()
{
    if (GlobalSettings::instance()->settings().valueBool("model.settings.browsing.enabled"))
        Saplings::mBrowsingPressure = GlobalSettings::instance()->settings().valueDouble("model.settings.browsing.browsingPressure");
    else
        Saplings::mBrowsingPressure = 0.;
}

bool Saplings::growSapling(const ResourceUnit *ru, SaplingCell &scell, SaplingTree &tree, int isc, float dom_height, float lif_value)
{
    ResourceUnitSpecies *rus = const_cast<ResourceUnitSpecies*>(ru->ruSpecies()[tree.species_index]);
    const Species *species = rus->species();

    // (1) calculate height growth potential for the tree (uses linerization of expressions...)
    double h_pot = species->saplingGrowthParameters().heightGrowthPotential.calculate(tree.height);
    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.
    if (dom_height==0.f)
        throw IException(QString("growSapling: height grid at %1/%2 has value 0").arg(isc));

    double rel_height = tree.height / dom_height;

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

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

    rus->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 = rus->prod3PG().fEnvYear();

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

    // sprouts grow faster. Sprouts therefore are less prone to stress (threshold), and can grow higher than the growth potential.
    if (tree.is_sprout())
        delta_h_factor = delta_h_factor *species->saplingGrowthParameters().sproutGrowth;

    // check browsing
    if (mBrowsingPressure>0. && tree.height<=2.f) {
        double p = rus->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.stress_years++;
        if (tree.stress_years > species->saplingGrowthParameters().maxStressYears) {
            // sapling dies...
            rus->saplingStat().addCarbonOfDeadSapling( tree.height / species->saplingGrowthParameters().hdSapling * 100.f );
            tree.clear();
            return true; // need cleanup
        }
    } else {
        tree.stress_years=0; // reset stress counter
    }
 
    DBG_IF(delta_h_pot*delta_h_factor < 0.f || (!tree.is_sprout() && delta_h_pot*delta_h_factor > 2.), "Sapling::growSapling", "inplausible height growth.");

    // grow
    tree.height += delta_h_pot * delta_h_factor;
    tree.age++; // increase age of sapling by 1

    // recruitment?
    if (tree.height > 4.f) {
        rus->saplingStat().mRecruited++;

        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 = static_cast<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 = const_cast<ResourceUnit*>(ru)->newTree();

            bigtree.setPosition(GlobalSettings::instance()->model()->grid()->indexOf(isc));
            // add variation: add +/-10% to dbh and *independently* to height.
            bigtree.setDbh(static_cast<float>(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)));
            bigtree.setHeight(static_cast<float>(tree.height * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)));
            bigtree.setSpecies( const_cast<Species*>(species) );
            bigtree.setAge(tree.age,tree.height);
            bigtree.setRU(const_cast<ResourceUnit*>(ru));
            bigtree.setup();
            const Tree *t = &bigtree;
            const_cast<ResourceUnitSpecies*>(rus)->statistics().add(t, 0); // count the newly created trees already in the stats
        }
        // clear all regeneration from this pixel (including this tree)
        tree.clear(); // clear this tree (no carbon flow to the ground)
        for (int i=0;i<NSAPCELLS;++i) {
            if (scell.saplings[i].is_occupied()) {
                // add carbon to the ground
                rus->saplingStat().addCarbonOfDeadSapling( scell.saplings[i].height / species->saplingGrowthParameters().hdSapling * 100.f );
                scell.saplings[i].clear();
            }
        }
        return true; // need cleanup
    }
    // book keeping (only for survivors) for the sapling of the resource unit / species
    SaplingStat &ss = rus->saplingStat();
    ss.mLiving++;
    ss.mAvgHeight+=tree.height;
    ss.mAvgAge+=tree.age;
    ss.mAvgDeltaHPot+=delta_h_pot;
    ss.mAvgHRealized += delta_h_pot * delta_h_factor;
    return false;
}

void SaplingStat::clearStatistics()
{
    mRecruited=mDied=mLiving=0;
    mSumDbhDied=0.;
    mAvgHeight=0.;
    mAvgAge=0.;
    mAvgDeltaHPot=mAvgHRealized=0.;
    mAdded=0;

}

void SaplingStat::calculate(const Species *species, ResourceUnit *ru)
{
    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>0) {
        // 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()  );

        DBGMODE(
        if (isnan(mCarbonLiving.C))
            qDebug("carbon NaN in SaplingStat::calculate (living trees).");
                );

        // turnover
        if (ru->snag())
            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()  );
                DBGMODE(
                if (isnan(dead_fine.C))
                    qDebug("carbon NaN in SaplingStat::calculate (self thinning).");
                        );

            }
        }

    }
    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()  );
        DBGMODE(
        if (isnan(dead_fine.C))
            qDebug("carbon NaN in SaplingStat::calculate (died trees).");
                );

    }
    if (!dead_wood.isEmpty() || !dead_fine.isEmpty())
        if (ru->snag())
            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();


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

}

double SaplingStat::livingStemNumber(const Species *species, double &rAvgDbh, double &rAvgHeight, double &rAvgAge) const
{
     rAvgHeight = averageHeight();
     rAvgDbh = rAvgHeight / species->saplingGrowthParameters().hdSapling * 100.f;
     rAvgAge = averageAge();
     double n= species->saplingGrowthParameters().representedStemNumber(rAvgDbh);
     return n;
// *** old code (sapling.cpp) ***
//    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;
}
 

Rev 1164	Rev 1165
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/saplings.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/saplings.cpp':
2	#include "global.h"	2	#include "global.h"
3	#include "saplings.h"	3	#include "saplings.h"
4		4
5	#include "globalsettings.h"	5	#include "globalsettings.h"
6	#include "model.h"	6	#include "model.h"
7	#include "resourceunit.h"	7	#include "resourceunit.h"
8	#include "resourceunitspecies.h"	8	#include "resourceunitspecies.h"
9	#include "establishment.h"	9	#include "establishment.h"
10	#include "species.h"	10	#include "species.h"
11	#include "seeddispersal.h"	11	#include "seeddispersal.h"
12		12
13	double Saplings::mRecruitmentVariation = 0.1; // +/- 10%	13	double Saplings::mRecruitmentVariation = 0.1; // +/- 10%
14	double Saplings::mBrowsingPressure = 0.;	14	double Saplings::mBrowsingPressure = 0.;
15		15
16		16
17	Saplings::Saplings()	17	Saplings::Saplings()
18	{	18	{
19		19
20	}	20	}
21		21
22	void Saplings::setup()	22	void Saplings::setup()
23	{	23	{
24	//mGrid.setup(GlobalSettings::instance()->model()->grid()->metricRect(), GlobalSettings::instance()->model()->grid()->cellsize());	24	//mGrid.setup(GlobalSettings::instance()->model()->grid()->metricRect(), GlobalSettings::instance()->model()->grid()->cellsize());
25	FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();	25	FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();
26	// mask out out-of-project areas	26	// mask out out-of-project areas
27	HeightGrid *hg = GlobalSettings::instance()->model()->heightGrid();	27	HeightGrid *hg = GlobalSettings::instance()->model()->heightGrid();
28	for (int i=0; i<lif_grid->count(); ++i) {	28	for (int i=0; i<lif_grid->count(); ++i) {
29	SaplingCell *s = cell(lif_grid->indexOf(i), false); // false: retrieve also invalid cells	29	SaplingCell *s = cell(lif_grid->indexOf(i), false); // false: retrieve also invalid cells
30	if (s) {	30	if (s) {
31	if (!hg->valueAtIndex(lif_grid->index5(i)).isValid())	31	if (!hg->valueAtIndex(lif_grid->index5(i)).isValid())
32	s->state = SaplingCell::CellInvalid;	32	s->state = SaplingCell::CellInvalid;
33	else	33	else
34	s->state = SaplingCell::CellFree;	34	s->state = SaplingCell::CellFree;
35	}	35	}
36		36
37	}	37	}
38		38
39	}	39	}
40		40
41	void Saplings::establishment(const ResourceUnit *ru)	41	void Saplings::establishment(const ResourceUnit *ru)
42	{	42	{
43	HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid();	43	HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid();
44	FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();	44	FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();
45		45
46	QPoint imap = ru->cornerPointOffset(); // offset on LIF/saplings grid	46	QPoint imap = ru->cornerPointOffset(); // offset on LIF/saplings grid
47	QPoint iseedmap = QPoint(imap.x()/10, imap.y()/10); // seed-map has 20m resolution, LIF 2m -> factor 10	47	QPoint iseedmap = QPoint(imap.x()/10, imap.y()/10); // seed-map has 20m resolution, LIF 2m -> factor 10
48		48
49	for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i)	49	for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i)
50	(*i)->saplingStat().clearStatistics();	50	(*i)->saplingStat().clearStatistics();
51		51
52	double lif_corr[cPxPerHectare];	52	double lif_corr[cPxPerHectare];
53	for (int i=0;i<cPxPerHectare;++i)	53	for (int i=0;i<cPxPerHectare;++i)
54	lif_corr[i]=-1.;	54	lif_corr[i]=-1.;
55		55
56	int species_idx;	56	int species_idx;
57	QVector<int>::const_iterator sbegin, send;	57	QVector<int>::const_iterator sbegin, send;
58	ru->speciesSet()->randomSpeciesOrder(sbegin, send);	58	ru->speciesSet()->randomSpeciesOrder(sbegin, send);
59	for (QVector<int>::const_iterator s_idx=sbegin; s_idx!=send;++s_idx) {	59	for (QVector<int>::const_iterator s_idx=sbegin; s_idx!=send;++s_idx) {
60	//for (int s_idx = 0; s_idx<ru->ruSpecies().size(); ++s_idx) {	-
61		60
62	// start from a random species (and cycle through the available species)	61	// start from a random species (and cycle through the available species)
63	species_idx = *s_idx;	62	species_idx = *s_idx;
64		63
65	ResourceUnitSpecies *rus = ru->ruSpecies()[species_idx];	64	ResourceUnitSpecies *rus = ru->ruSpecies()[species_idx];
66	// check if there are seeds of the given species on the resource unit	65	// check if there are seeds of the given species on the resource unit
67	float seeds = 0.f;	66	float seeds = 0.f;
68	Grid<float> &seedmap = const_cast<Grid<float>& >(rus->species()->seedDispersal()->seedMap());	67	Grid<float> &seedmap = const_cast<Grid<float>& >(rus->species()->seedDispersal()->seedMap());
69	for (int iy=0;iy<5;++iy) {	68	for (int iy=0;iy<5;++iy) {
70	float *p = seedmap.ptr(iseedmap.x(), iseedmap.y());	69	float *p = seedmap.ptr(iseedmap.x(), iseedmap.y());
71	for (int ix=0;ix<5;++ix)	70	for (int ix=0;ix<5;++ix)
72	seeds += *p++;	71	seeds += *p++;
73	}	72	}
74	// if there are no seeds: no need to do more	73	// if there are no seeds: no need to do more
75	if (seeds==0.f)	74	if (seeds==0.f)
76	continue;	75	continue;
77		76
78	// calculate the abiotic environment (TACA)	77	// calculate the abiotic environment (TACA)
79	rus->establishment().calculateAbioticEnvironment();	78	rus->establishment().calculateAbioticEnvironment();
80	double abiotic_env = rus->establishment().abioticEnvironment();	79	double abiotic_env = rus->establishment().abioticEnvironment();
81	if (abiotic_env==0.)	80	if (abiotic_env==0.)
82	continue;	81	continue;
83		82
84	// loop over all 2m cells on this resource unit	83	// loop over all 2m cells on this resource unit
85	SaplingCell *sap_cells = ru->saplingCellArray();	84	SaplingCell *sap_cells = ru->saplingCellArray();
86	SaplingCell *s;	85	SaplingCell *s;
87	int isc = 0; // index on 2m cell	86	int isc = 0; // index on 2m cell
88	for (int iy=0; iy<cPxPerRU; ++iy) {	87	for (int iy=0; iy<cPxPerRU; ++iy) {
89	//s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row	88	//s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row
90	s = &sap_cells[iy*cPxPerRU]; // pointer to a row	89	s = &sap_cells[iy*cPxPerRU]; // pointer to a row
91	isc = lif_grid->index(imap.x(), imap.y()+iy);	90	isc = lif_grid->index(imap.x(), imap.y()+iy);
92		91
93	for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) {	92	for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) {
94	if (s->state == SaplingCell::CellFree) {	93	if (s->state == SaplingCell::CellFree) {
95	// is a sapling of the current species already on the pixel?	94	// is a sapling of the current species already on the pixel?
96	// * test for sapling height already in cell state	95	// * test for sapling height already in cell state
97	// * test for grass-cover already in cell state	96	// * test for grass-cover already in cell state
98	SaplingTree *stree=0;	97	SaplingTree *stree=0;
99	SaplingTree *slot=s->saplings;	98	SaplingTree *slot=s->saplings;
100	for (int i=0;i<NSAPCELLS;++i, ++slot) {	99	for (int i=0;i<NSAPCELLS;++i, ++slot) {
101	if (!stree && !slot->is_occupied())	100	if (!stree && !slot->is_occupied())
102	stree=slot;	101	stree=slot;
103	if (slot->species_index == species_idx) {	102	if (slot->species_index == species_idx) {
104	stree=0;	103	stree=0;
105	break;	104	break;
106	}	105	}
107	}	106	}
108		107
109	if (stree) {	108	if (stree) {
110	// grass cover?	109	// grass cover?
111	float seed_map_value = seedmap[lif_grid->index10(isc)];	110	float seed_map_value = seedmap[lif_grid->index10(isc)];
112	if (seed_map_value==0.f)	111	if (seed_map_value==0.f)
113	continue;	112	continue;
114	const HeightGridValue &hgv = (*height_grid)[lif_grid->index5(isc)];	113	const HeightGridValue &hgv = (*height_grid)[lif_grid->index5(isc)];
115	float lif_value = (*lif_grid)[isc];	114	float lif_value = (*lif_grid)[isc];
116		115
117	double &lif_corrected = lif_corr[iy*cPxPerRU+ix];	116	double &lif_corrected = lif_corr[iy*cPxPerRU+ix];
118	// calculate the LIFcorrected only once per pixel	117	// calculate the LIFcorrected only once per pixel
119	if (lif_corrected<0.)	118	if (lif_corrected<0.)
120	lif_corrected = rus->species()->speciesSet()->LRIcorrection(lif_value, 4. / hgv.height);	119	lif_corrected = rus->species()->speciesSet()->LRIcorrection(lif_value, 4. / hgv.height);
121		120
122	// check for the combination of seed availability and light on the forest floor	121	// check for the combination of seed availability and light on the forest floor
123	if (drandom() < seed_map_valuelif_correctedabiotic_env ) {	122	if (drandom() < seed_map_valuelif_correctedabiotic_env ) {
124	// ok, lets add a sapling at the given position (age is incremented later)	123	// ok, lets add a sapling at the given position (age is incremented later)
125	stree->setSapling(0.05f, 0, species_idx);	124	stree->setSapling(0.05f, 0, species_idx);
126	s->checkState();	125	s->checkState();
127	rus->saplingStat().mAdded++;	126	rus->saplingStat().mAdded++;
128		127
129	}	128	}
130		129
131	}	130	}
132		131
133	}	132	}
134	}	133	}
135	}	134	}
136		135
137	}	136	}
138		137
139	}	138	}
140		139
141	void Saplings::saplingGrowth(const ResourceUnit *ru)	140	void Saplings::saplingGrowth(const ResourceUnit *ru)
142	{	141	{
143	HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid();	142	HeightGrid *height_grid = GlobalSettings::instance()->model()->heightGrid();
144	FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();	143	FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();
145		144
146	QPoint imap = ru->cornerPointOffset();	145	QPoint imap = ru->cornerPointOffset();
147	bool need_check=false;	146	bool need_check=false;
148	SaplingCell *sap_cells = ru->saplingCellArray();	147	SaplingCell *sap_cells = ru->saplingCellArray();
149	for (int iy=0; iy<cPxPerRU; ++iy) {	148	for (int iy=0; iy<cPxPerRU; ++iy) {
150	//SaplingCell *s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row	149	//SaplingCell *s = mGrid.ptr(imap.x(), imap.y()+iy); // ptr to the row
151	SaplingCell s = &sap_cells[iycPxPerRU]; // ptr to row	150	SaplingCell s = &sap_cells[iycPxPerRU]; // ptr to row
152	int isc = lif_grid->index(imap.x(), imap.y()+iy);	151	int isc = lif_grid->index(imap.x(), imap.y()+iy);
153		152
154	for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) {	153	for (int ix=0;ix<cPxPerRU; ++ix, ++s, ++isc) {
155	if (s->state != SaplingCell::CellInvalid) {	154	if (s->state != SaplingCell::CellInvalid) {
156	need_check=false;	155	need_check=false;
157	for (int i=0;i<NSAPCELLS;++i) {	156	for (int i=0;i<NSAPCELLS;++i) {
158	if (s->saplings[i].is_occupied()) {	157	if (s->saplings[i].is_occupied()) {
159	// growth of this sapling tree	158	// growth of this sapling tree
160	const HeightGridValue &hgv = (*height_grid)[height_grid->index5(isc)];	159	const HeightGridValue &hgv = (*height_grid)[height_grid->index5(isc)];
161	float lif_value = (*lif_grid)[isc];	160	float lif_value = (*lif_grid)[isc];
162		161
163	need_check \|= growSapling(ru, *s, s->saplings[i], isc, hgv.height, lif_value);	162	need_check \|= growSapling(ru, *s, s->saplings[i], isc, hgv.height, lif_value);
164	}	163	}
165	}	164	}
166	if (need_check)	165	if (need_check)
167	s->checkState();	166	s->checkState();
168	}	167	}
169	}	168	}
170	}	169	}
171		170
172		171
173		172
174		173
175	// store statistics on saplings/regeneration	174	// store statistics on saplings/regeneration
176	for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i) {	175	for (QList<ResourceUnitSpecies*>::const_iterator i=ru->ruSpecies().constBegin(); i!=ru->ruSpecies().constEnd(); ++i) {
177	(i)->saplingStat().calculate((i)->species(), const_cast<ResourceUnit*>(ru));	176	(i)->saplingStat().calculate((i)->species(), const_cast<ResourceUnit*>(ru));
178	(i)->statistics().add(&((i)->saplingStat()));	177	(i)->statistics().add(&((i)->saplingStat()));
179	}	178	}
180	}	179	}
181		180
182	SaplingCell Saplings::cell(QPoint lif_coords, bool only_valid, ResourceUnit *rRUPtr)	181	SaplingCell Saplings::cell(QPoint lif_coords, bool only_valid, ResourceUnit *rRUPtr)
183	{	182	{
184	FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();	183	FloatGrid *lif_grid = GlobalSettings::instance()->model()->grid();
185		184
186	// in this case, getting the actual cell is quite cumbersome: first, retrieve the resource unit, then the	185	// in this case, getting the actual cell is quite cumbersome: first, retrieve the resource unit, then the
187	// cell based on the offset of the given coordiantes relative to the corner of the resource unit.	186	// cell based on the offset of the given coordiantes relative to the corner of the resource unit.
188	ResourceUnit *ru = GlobalSettings::instance()->model()->ru(lif_grid->cellCenterPoint(lif_coords));	187	ResourceUnit *ru = GlobalSettings::instance()->model()->ru(lif_grid->cellCenterPoint(lif_coords));
189	if (rRUPtr)	188	if (rRUPtr)
190	*rRUPtr = ru;	189	*rRUPtr = ru;
191		190
192	if (ru) {	191	if (ru) {
193	QPoint local_coords = lif_coords - ru->cornerPointOffset();	192	QPoint local_coords = lif_coords - ru->cornerPointOffset();
194	int idx = local_coords.y() * cPxPerRU + local_coords.x();	193	int idx = local_coords.y() * cPxPerRU + local_coords.x();
195	DBGMODE( if (idx<0 \|\| idx>=cPxPerHectare)	194	DBGMODE( if (idx<0 \|\| idx>=cPxPerHectare)
196	qDebug("invalid coords in Saplings::cell");	195	qDebug("invalid coords in Saplings::cell");
197	);	196	);
198	SaplingCell *s=&ru->saplingCellArray()[idx];	197	SaplingCell *s=&ru->saplingCellArray()[idx];
199	if (s && (!only_valid \|\| s->state!=SaplingCell::CellInvalid))	198	if (s && (!only_valid \|\| s->state!=SaplingCell::CellInvalid))
200	return s;	199	return s;
201	}	200	}
202	return 0;	201	return 0;
203	}	202	}
204		203
205	void Saplings::clearSaplings(const QRectF &rectangle, const bool remove_biomass)	204	void Saplings::clearSaplings(const QRectF &rectangle, const bool remove_biomass)
206	{	205	{
207	GridRunner<float> runner(GlobalSettings::instance()->model()->grid(), rectangle);	206	GridRunner<float> runner(GlobalSettings::instance()->model()->grid(), rectangle);
208	ResourceUnit *ru;	207	ResourceUnit *ru;
209	while (runner.next()) {	208	while (runner.next()) {
210	SaplingCell *s = cell(runner.currentIndex(), true, &ru);	209	SaplingCell *s = cell(runner.currentIndex(), true, &ru);
211	if (s) {	210	if (s) {
-		211	clearSaplings(s, ru, remove_biomass);
-		212	}
212		213
213	for (int i=0;i<NSAPCELLS;++i)	-
214	if (s->saplings[i].is_occupied()) {	-
215	if (remove_biomass) {	-
216	ResourceUnitSpecies *rus = ru->resourceUnitSpecies(s->saplings[i].species_index);	-
217	if (!rus && !rus->species()) {	-
218	qDebug() << "Saplings::clearSaplings(): invalid resource unit!!!";	-
219	return;	-
220	}	-
221	rus->saplingStat().addCarbonOfDeadSapling( s->saplings[i].height / rus->species()->saplingGrowthParameters().hdSapling * 100.f );	-
-		214	}
-		215	}
-		216
-		217	void Saplings::clearSaplings(SaplingCell s, ResourceUnit ru, const bool remove_biomass)
-		218	{
-		219	if (s) {
-		220	for (int i=0;i<NSAPCELLS;++i)
-		221	if (s->saplings[i].is_occupied()) {
-		222	if (!remove_biomass) {
-		223	ResourceUnitSpecies *rus = ru->resourceUnitSpecies(s->saplings[i].species_index);
-		224	if (!rus && !rus->species()) {
-		225	qDebug() << "Saplings::clearSaplings(): invalid resource unit!!!";
-		226	return;
222	}	227	}
223	s->saplings[i].clear();	-
-		228	rus->saplingStat().addCarbonOfDeadSapling( s->saplings[i].height / rus->species()->saplingGrowthParameters().hdSapling * 100.f );
224	}	229	}
225	s->checkState();	-
-		230	s->saplings[i].clear();
-		231	}
-		232	s->checkState();
-		233
-		234	}
-		235
-		236	}
-		237
-		238	int Saplings::addSprout(const Tree *t)
-		239	{
-		240	if (t->species()->saplingGrowthParameters().sproutGrowth==0.)
-		241	return 0;
-		242	SaplingCell *sc = cell(t->positionIndex());
-		243	if (!sc)
-		244	return 0;
-		245	clearSaplings(sc, const_cast<ResourceUnit*>(t->ru()), false );
-		246	SaplingTree *st=sc->addSapling(0.05f, 0, t->species()->index());
-		247	if (st)
-		248	st->set_sprout(true);
-		249
-		250	// neighboring cells
-		251	double crown_area = t->crownRadius()t->crownRadius() M_PI; //m2
-		252	// calculate how many cells on the ground are covered by the crown (this is a rather rough estimate)
-		253	// n_cells: in addition to the original cell
-		254	int n_cells = static_cast<int>(round( crown_area / static_cast<double>(cPxSize*cPxSize) - 1.));
-		255	if (n_cells>0) {
-		256	ResourceUnit *ru;
-		257	static const int offsets_x[8] = {1,1,0,-1,-1,-1,0,1};
-		258	static const int offsets_y[8] = {0,1,1,1,0,-1,-1,-1};
-		259	int s=irandom(0,8);
-		260	while(n_cells) {
-		261	sc = cell(t->positionIndex()+QPoint(offsets_x[s], offsets_y[s]),true,&ru);
-		262	if (sc) {
-		263	clearSaplings(sc, ru, false );
-		264	SaplingTree *st=sc->addSapling(0.05f, 0, t->species()->index());
-		265	if (st)
-		266	st->set_sprout(true);
-		267	}
226		268
-		269	s = (s+1)%8; --n_cells;
227	}	270	}
228	}	271	}
-		272	return 1;
229	}	273	}
230		274
231	void Saplings::updateBrowsingPressure()	275	void Saplings::updateBrowsingPressure()
232	{	276	{
233	if (GlobalSettings::instance()->settings().valueBool("model.settings.browsing.enabled"))	277	if (GlobalSettings::instance()->settings().valueBool("model.settings.browsing.enabled"))
234	Saplings::mBrowsingPressure = GlobalSettings::instance()->settings().valueDouble("model.settings.browsing.browsingPressure");	278	Saplings::mBrowsingPressure = GlobalSettings::instance()->settings().valueDouble("model.settings.browsing.browsingPressure");
235	else	279	else
236	Saplings::mBrowsingPressure = 0.;	280	Saplings::mBrowsingPressure = 0.;
237	}	281	}
238		282
239	bool Saplings::growSapling(const ResourceUnit *ru, SaplingCell &scell, SaplingTree &tree, int isc, float dom_height, float lif_value)	283	bool Saplings::growSapling(const ResourceUnit *ru, SaplingCell &scell, SaplingTree &tree, int isc, float dom_height, float lif_value)
240	{	284	{
241	ResourceUnitSpecies rus = const_cast<ResourceUnitSpecies>(ru->ruSpecies()[tree.species_index]);	285	ResourceUnitSpecies rus = const_cast<ResourceUnitSpecies>(ru->ruSpecies()[tree.species_index]);
242	const Species *species = rus->species();	286	const Species *species = rus->species();
243		287
244	// (1) calculate height growth potential for the tree (uses linerization of expressions...)	288	// (1) calculate height growth potential for the tree (uses linerization of expressions...)
245	double h_pot = species->saplingGrowthParameters().heightGrowthPotential.calculate(tree.height);	289	double h_pot = species->saplingGrowthParameters().heightGrowthPotential.calculate(tree.height);
246	double delta_h_pot = h_pot - tree.height;	290	double delta_h_pot = h_pot - tree.height;
247		291
248	// (2) reduce height growth potential with species growth response f_env_yr and with light state (i.e. LIF-value) of home-pixel.	292	// (2) reduce height growth potential with species growth response f_env_yr and with light state (i.e. LIF-value) of home-pixel.
249	if (dom_height==0.f)	293	if (dom_height==0.f)
250	throw IException(QString("growSapling: height grid at %1/%2 has value 0").arg(isc));	294	throw IException(QString("growSapling: height grid at %1/%2 has value 0").arg(isc));
251		295
252	double rel_height = tree.height / dom_height;	296	double rel_height = tree.height / dom_height;
253		297
254	double lif_corrected = species->speciesSet()->LRIcorrection(lif_value, rel_height); // correction based on height	298	double lif_corrected = species->speciesSet()->LRIcorrection(lif_value, rel_height); // correction based on height
255		299
256	double lr = species->lightResponse(lif_corrected); // species specific light response (LUI, light utilization index)	300	double lr = species->lightResponse(lif_corrected); // species specific light response (LUI, light utilization index)
257		301
258	rus->calculate(true); // calculate the 3pg module (this is done only if that did not happen up to now); true: call comes from regeneration	302	rus->calculate(true); // calculate the 3pg module (this is done only if that did not happen up to now); true: call comes from regeneration
259	double f_env_yr = rus->prod3PG().fEnvYear();	303	double f_env_yr = rus->prod3PG().fEnvYear();
260		304
261	double delta_h_factor = f_env_yr * lr; // relative growth	305	double delta_h_factor = f_env_yr * lr; // relative growth
262		306
263	if (h_pot<0. \|\| delta_h_pot<0. \|\| lif_corrected<0. \|\| lif_corrected>1. \|\| delta_h_factor<0. \|\| delta_h_factor>1. )	307	if (h_pot<0. \|\| delta_h_pot<0. \|\| lif_corrected<0. \|\| lif_corrected>1. \|\| delta_h_factor<0. \|\| delta_h_factor>1. )
264	qDebug() << "invalid values in Sapling::growSapling";	308	qDebug() << "invalid values in Sapling::growSapling";
-		309
-		310	// sprouts grow faster. Sprouts therefore are less prone to stress (threshold), and can grow higher than the growth potential.
-		311	if (tree.is_sprout())
-		312	delta_h_factor = delta_h_factor *species->saplingGrowthParameters().sproutGrowth;
265		313
266	// check browsing	314	// check browsing
267	if (mBrowsingPressure>0. && tree.height<=2.f) {	315	if (mBrowsingPressure>0. && tree.height<=2.f) {
268	double p = rus->species()->saplingGrowthParameters().browsingProbability;	316	double p = rus->species()->saplingGrowthParameters().browsingProbability;
269	// calculate modifed annual browsing probability via odds-ratios	317	// calculate modifed annual browsing probability via odds-ratios
270	// odds = p/(1-p) -> odds_mod = odds * browsingPressure -> p_mod = odds_mod /( 1 + odds_mod) === ppressure/(1-p+ppressure)	318	// odds = p/(1-p) -> odds_mod = odds * browsingPressure -> p_mod = odds_mod /( 1 + odds_mod) === ppressure/(1-p+ppressure)
271	double p_browse = pmBrowsingPressure / (1. - p + pmBrowsingPressure);	319	double p_browse = pmBrowsingPressure / (1. - p + pmBrowsingPressure);
272	if (drandom() < p_browse) {	320	if (drandom() < p_browse) {
273	delta_h_factor = 0.;	321	delta_h_factor = 0.;
274	}	322	}
275	}	323	}
276		324
277	// check mortality of saplings	325	// check mortality of saplings
278	if (delta_h_factor < species->saplingGrowthParameters().stressThreshold) {	326	if (delta_h_factor < species->saplingGrowthParameters().stressThreshold) {
279	tree.stress_years++;	327	tree.stress_years++;
280	if (tree.stress_years > species->saplingGrowthParameters().maxStressYears) {	328	if (tree.stress_years > species->saplingGrowthParameters().maxStressYears) {
281	// sapling dies...	329	// sapling dies...
282	rus->saplingStat().addCarbonOfDeadSapling( tree.height / species->saplingGrowthParameters().hdSapling * 100.f );	330	rus->saplingStat().addCarbonOfDeadSapling( tree.height / species->saplingGrowthParameters().hdSapling * 100.f );
283	tree.clear();	331	tree.clear();
284	return true; // need cleanup	332	return true; // need cleanup
285	}	333	}
286	} else {	334	} else {
287	tree.stress_years=0; // reset stress counter	335	tree.stress_years=0; // reset stress counter
288	}	336	}
289	DBG_IF(delta_h_potdelta_h_factor < 0.f \|\| delta_h_potdelta_h_factor > 2., "Sapling::growSapling", "inplausible height growth.");	-
-		337	DBG_IF(delta_h_potdelta_h_factor < 0.f \|\| (!tree.is_sprout() && delta_h_potdelta_h_factor > 2.), "Sapling::growSapling", "inplausible height growth.");
290		338
291	// grow	339	// grow
292	tree.height += delta_h_pot * delta_h_factor;	340	tree.height += delta_h_pot * delta_h_factor;
293	tree.age++; // increase age of sapling by 1	341	tree.age++; // increase age of sapling by 1
294		342
295	// recruitment?	343	// recruitment?
296	if (tree.height > 4.f) {	344	if (tree.height > 4.f) {
297	rus->saplingStat().mRecruited++;	345	rus->saplingStat().mRecruited++;
298		346
299	float dbh = tree.height / species->saplingGrowthParameters().hdSapling * 100.f;	347	float dbh = tree.height / species->saplingGrowthParameters().hdSapling * 100.f;
300	// the number of trees to create (result is in trees per pixel)	348	// the number of trees to create (result is in trees per pixel)
301	double n_trees = species->saplingGrowthParameters().representedStemNumber(dbh);	349	double n_trees = species->saplingGrowthParameters().representedStemNumber(dbh);
302	int to_establish = static_cast<int>( n_trees );	350	int to_establish = static_cast<int>( n_trees );
303		351
304	// if n_trees is not an integer, choose randomly if we should add a tree.	352	// if n_trees is not an integer, choose randomly if we should add a tree.
305	// e.g.: n_trees = 2.3 -> add 2 trees with 70% probability, and add 3 trees with p=30%.	353	// e.g.: n_trees = 2.3 -> add 2 trees with 70% probability, and add 3 trees with p=30%.
306	if (drandom() < (n_trees-to_establish) \|\| to_establish==0)	354	if (drandom() < (n_trees-to_establish) \|\| to_establish==0)
307	to_establish++;	355	to_establish++;
308		356
309	// add a new tree	357	// add a new tree
310	for (int i=0;i<to_establish;i++) {	358	for (int i=0;i<to_establish;i++) {
311	Tree &bigtree = const_cast<ResourceUnit*>(ru)->newTree();	359	Tree &bigtree = const_cast<ResourceUnit*>(ru)->newTree();
312		360
313	bigtree.setPosition(GlobalSettings::instance()->model()->grid()->indexOf(isc));	361	bigtree.setPosition(GlobalSettings::instance()->model()->grid()->indexOf(isc));
314	// add variation: add +/-10% to dbh and independently to height.	362	// add variation: add +/-10% to dbh and independently to height.
315	bigtree.setDbh(static_cast<float>(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)));	363	bigtree.setDbh(static_cast<float>(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)));
316	bigtree.setHeight(static_cast<float>(tree.height * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)));	364	bigtree.setHeight(static_cast<float>(tree.height * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)));
317	bigtree.setSpecies( const_cast<Species*>(species) );	365	bigtree.setSpecies( const_cast<Species*>(species) );
318	bigtree.setAge(tree.age,tree.height);	366	bigtree.setAge(tree.age,tree.height);
319	bigtree.setRU(const_cast<ResourceUnit*>(ru));	367	bigtree.setRU(const_cast<ResourceUnit*>(ru));
320	bigtree.setup();	368	bigtree.setup();
321	const Tree *t = &bigtree;	369	const Tree *t = &bigtree;
322	const_cast<ResourceUnitSpecies*>(rus)->statistics().add(t, 0); // count the newly created trees already in the stats	370	const_cast<ResourceUnitSpecies*>(rus)->statistics().add(t, 0); // count the newly created trees already in the stats
323	}	371	}
324	// clear all regeneration from this pixel (including this tree)	372	// clear all regeneration from this pixel (including this tree)
325	tree.clear(); // clear this tree (no carbon flow to the ground)	373	tree.clear(); // clear this tree (no carbon flow to the ground)
326	for (int i=0;i<NSAPCELLS;++i) {	374	for (int i=0;i<NSAPCELLS;++i) {
327	if (scell.saplings[i].is_occupied()) {	375	if (scell.saplings[i].is_occupied()) {
328	// add carbon to the ground	376	// add carbon to the ground
329	rus->saplingStat().addCarbonOfDeadSapling( scell.saplings[i].height / species->saplingGrowthParameters().hdSapling * 100.f );	377	rus->saplingStat().addCarbonOfDeadSapling( scell.saplings[i].height / species->saplingGrowthParameters().hdSapling * 100.f );
330	scell.saplings[i].clear();	378	scell.saplings[i].clear();
331	}	379	}
332	}	380	}
333	return true; // need cleanup	381	return true; // need cleanup
334	}	382	}
335	// book keeping (only for survivors) for the sapling of the resource unit / species	383	// book keeping (only for survivors) for the sapling of the resource unit / species
336	SaplingStat &ss = rus->saplingStat();	384	SaplingStat &ss = rus->saplingStat();
337	ss.mLiving++;	385	ss.mLiving++;
338	ss.mAvgHeight+=tree.height;	386	ss.mAvgHeight+=tree.height;
339	ss.mAvgAge+=tree.age;	387	ss.mAvgAge+=tree.age;
340	ss.mAvgDeltaHPot+=delta_h_pot;	388	ss.mAvgDeltaHPot+=delta_h_pot;
341	ss.mAvgHRealized += delta_h_pot * delta_h_factor;	389	ss.mAvgHRealized += delta_h_pot * delta_h_factor;
342	return false;	390	return false;
343	}	391	}
344		392
345	void SaplingStat::clearStatistics()	393	void SaplingStat::clearStatistics()
346	{	394	{
347	mRecruited=mDied=mLiving=0;	395	mRecruited=mDied=mLiving=0;
348	mSumDbhDied=0.;	396	mSumDbhDied=0.;
349	mAvgHeight=0.;	397	mAvgHeight=0.;
350	mAvgAge=0.;	398	mAvgAge=0.;
351	mAvgDeltaHPot=mAvgHRealized=0.;	399	mAvgDeltaHPot=mAvgHRealized=0.;
352	mAdded=0;	400	mAdded=0;
353		401
354	}	402	}
355		403
356	void SaplingStat::calculate(const Species species, ResourceUnit ru)	404	void SaplingStat::calculate(const Species species, ResourceUnit ru)
357	{	405	{
358	if (mLiving) {	406	if (mLiving) {
359	mAvgHeight /= double(mLiving);	407	mAvgHeight /= double(mLiving);
360	mAvgAge /= double(mLiving);	408	mAvgAge /= double(mLiving);
361	mAvgDeltaHPot /= double(mLiving);	409	mAvgDeltaHPot /= double(mLiving);
362	mAvgHRealized /= double(mLiving);	410	mAvgHRealized /= double(mLiving);
363	}	411	}
364		412
365	// calculate carbon balance	413	// calculate carbon balance
366	CNPair old_state = mCarbonLiving;	414	CNPair old_state = mCarbonLiving;
367	mCarbonLiving.clear();	415	mCarbonLiving.clear();
368		416
369	CNPair dead_wood, dead_fine; // pools for mortality	417	CNPair dead_wood, dead_fine; // pools for mortality
370	// average dbh	418	// average dbh
371	if (mLiving>0) {	419	if (mLiving>0) {
372	// calculate the avg dbh and number of stems	420	// calculate the avg dbh and number of stems
373	double avg_dbh = mAvgHeight / species->saplingGrowthParameters().hdSapling * 100.;	421	double avg_dbh = mAvgHeight / species->saplingGrowthParameters().hdSapling * 100.;
374	double n = mLiving * species->saplingGrowthParameters().representedStemNumber( avg_dbh );	422	double n = mLiving * species->saplingGrowthParameters().representedStemNumber( avg_dbh );
375	// woody parts: stem, branchse and coarse roots	423	// woody parts: stem, branchse and coarse roots
376	double woody_bm = species->biomassWoody(avg_dbh) + species->biomassBranch(avg_dbh) + species->biomassRoot(avg_dbh);	424	double woody_bm = species->biomassWoody(avg_dbh) + species->biomassBranch(avg_dbh) + species->biomassRoot(avg_dbh);
377	double foliage = species->biomassFoliage(avg_dbh);	425	double foliage = species->biomassFoliage(avg_dbh);
378	double fineroot = foliage*species->finerootFoliageRatio();	426	double fineroot = foliage*species->finerootFoliageRatio();
379		427
380	mCarbonLiving.addBiomass( woody_bm*n, species->cnWood() );	428	mCarbonLiving.addBiomass( woody_bm*n, species->cnWood() );
381	mCarbonLiving.addBiomass( foliage*n, species->cnFoliage() );	429	mCarbonLiving.addBiomass( foliage*n, species->cnFoliage() );
382	mCarbonLiving.addBiomass( fineroot*n, species->cnFineroot() );	430	mCarbonLiving.addBiomass( fineroot*n, species->cnFineroot() );
383		431
384	DBGMODE(	432	DBGMODE(
385	if (isnan(mCarbonLiving.C))	433	if (isnan(mCarbonLiving.C))
386	qDebug("carbon NaN in SaplingStat::calculate (living trees).");	434	qDebug("carbon NaN in SaplingStat::calculate (living trees).");
387	);	435	);
388		436
389	// turnover	437	// turnover
390	if (ru->snag())	438	if (ru->snag())
391	ru->snag()->addTurnoverLitter(species, foliagespecies->turnoverLeaf(), finerootspecies->turnoverRoot());	439	ru->snag()->addTurnoverLitter(species, foliagespecies->turnoverLeaf(), finerootspecies->turnoverRoot());
392		440
393	// calculate the "mortality from competition", i.e. carbon that stems from reduction of stem numbers	441	// calculate the "mortality from competition", i.e. carbon that stems from reduction of stem numbers
394	// from Reinekes formula.	442	// from Reinekes formula.
395	//	443	//
396	if (avg_dbh>1.) {	444	if (avg_dbh>1.) {
397	double avg_dbh_before = (mAvgHeight - mAvgHRealized) / species->saplingGrowthParameters().hdSapling * 100.;	445	double avg_dbh_before = (mAvgHeight - mAvgHRealized) / species->saplingGrowthParameters().hdSapling * 100.;
398	double n_before = mLiving * species->saplingGrowthParameters().representedStemNumber( qMax(1.,avg_dbh_before) );	446	double n_before = mLiving * species->saplingGrowthParameters().representedStemNumber( qMax(1.,avg_dbh_before) );
399	if (n<n_before) {	447	if (n<n_before) {
400	dead_wood.addBiomass( woody_bm * (n_before-n), species->cnWood() );	448	dead_wood.addBiomass( woody_bm * (n_before-n), species->cnWood() );
401	dead_fine.addBiomass( foliage * (n_before-n), species->cnFoliage() );	449	dead_fine.addBiomass( foliage * (n_before-n), species->cnFoliage() );
402	dead_fine.addBiomass( fineroot * (n_before-n), species->cnFineroot() );	450	dead_fine.addBiomass( fineroot * (n_before-n), species->cnFineroot() );
403	DBGMODE(	451	DBGMODE(
404	if (isnan(dead_fine.C))	452	if (isnan(dead_fine.C))
405	qDebug("carbon NaN in SaplingStat::calculate (self thinning).");	453	qDebug("carbon NaN in SaplingStat::calculate (self thinning).");
406	);	454	);
407		455
408	}	456	}
409	}	457	}
410		458
411	}	459	}
412	if (mDied) {	460	if (mDied) {
413	double avg_dbh_dead = mSumDbhDied / double(mDied);	461	double avg_dbh_dead = mSumDbhDied / double(mDied);
414	double n = mDied * species->saplingGrowthParameters().representedStemNumber( avg_dbh_dead );	462	double n = mDied * species->saplingGrowthParameters().representedStemNumber( avg_dbh_dead );
415	// woody parts: stem, branchse and coarse roots	463	// woody parts: stem, branchse and coarse roots
416		464
417	dead_wood.addBiomass( ( species->biomassWoody(avg_dbh_dead) + species->biomassBranch(avg_dbh_dead) + species->biomassRoot(avg_dbh_dead)) * n, species->cnWood() );	465	dead_wood.addBiomass( ( species->biomassWoody(avg_dbh_dead) + species->biomassBranch(avg_dbh_dead) + species->biomassRoot(avg_dbh_dead)) * n, species->cnWood() );
418	double foliage = species->biomassFoliage(avg_dbh_dead)*n;	466	double foliage = species->biomassFoliage(avg_dbh_dead)*n;
419		467
420	dead_fine.addBiomass( foliage, species->cnFoliage() );	468	dead_fine.addBiomass( foliage, species->cnFoliage() );
421	dead_fine.addBiomass( foliage*species->finerootFoliageRatio(), species->cnFineroot() );	469	dead_fine.addBiomass( foliage*species->finerootFoliageRatio(), species->cnFineroot() );
422	DBGMODE(	470	DBGMODE(
423	if (isnan(dead_fine.C))	471	if (isnan(dead_fine.C))
424	qDebug("carbon NaN in SaplingStat::calculate (died trees).");	472	qDebug("carbon NaN in SaplingStat::calculate (died trees).");
425	);	473	);
426		474
427	}	475	}
428	if (!dead_wood.isEmpty() \|\| !dead_fine.isEmpty())	476	if (!dead_wood.isEmpty() \|\| !dead_fine.isEmpty())
429	if (ru->snag())	477	if (ru->snag())
430	ru->snag()->addToSoil(species, dead_wood, dead_fine);	478	ru->snag()->addToSoil(species, dead_wood, dead_fine);
431		479
432	// calculate net growth:	480	// calculate net growth:
433	// delta of stocks	481	// delta of stocks
434	mCarbonGain = mCarbonLiving + dead_fine + dead_wood - old_state;	482	mCarbonGain = mCarbonLiving + dead_fine + dead_wood - old_state;
435	if (mCarbonGain.C < 0)	483	if (mCarbonGain.C < 0)
436	mCarbonGain.clear();	484	mCarbonGain.clear();
437		485
438		486
439	GlobalSettings::instance()->systemStatistics()->saplingCount+=mLiving;	487	GlobalSettings::instance()->systemStatistics()->saplingCount+=mLiving;

Subversion Repositories public iLand

(root)/src/core/saplings.cpp @ 1222 - Rev 1164 → 1165