Subversion Repositories public iLand

#include "global.h"

#include "tree.h"

#include "grid.h"

#include "stamp.h"

#include "species.h"

#include "resourceunit.h"

#include "model.h"

#include "snag.h"

// static varaibles

FloatGrid *Tree::mGrid = 0;

HeightGrid *Tree::mHeightGrid = 0;

int Tree::m_statPrint=0;

int Tree::m_statAboveZ=0;

int Tree::m_statCreated=0;

int Tree::m_nextId=0;

/** get distance and direction between two points.

  returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */

float dist_and_direction(const QPointF &PStart, const QPointF &PEnd, float &rAngle)

{

    float dx = PEnd.x() - PStart.x();

    float dy = PEnd.y() - PStart.y();

    float d = sqrt(dx*dx + dy*dy);

    // direction:

    rAngle = atan2(dx, dy);

    return d;

}

// lifecycle

Tree::Tree()

{

    mDbh = mHeight = 0;

    mRU = 0; mSpecies = 0;

    mFlags = mAge = 0;

    mOpacity=mFoliageMass=mWoodyMass=mCoarseRootMass=mFineRootMass=mLeafArea=0.;

    mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.;

    mLightResponse = 0.;

    mId = m_nextId++;

    m_statCreated++;

}

float Tree::crownRadius() const

{

    Q_ASSERT(mStamp!=0);

    return mStamp->crownRadius();

}

void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid)

{

    mGrid = gridToStamp; mHeightGrid = dominanceGrid;

}

/// dumps some core variables of a tree to a string.

QString Tree::dump()

{

    QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8")

                            .arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight)

                            .arg(position().x()).arg(position().y())

                            .arg(mRU->index()).arg(mLRI);

    return result;

}

void Tree::dumpList(DebugList &rTargetList)

{

    rTargetList << mId << species()->id() << mDbh << mHeight  << position().x() << position().y()   << mRU->index() << mLRI

                << mWoodyMass << mCoarseRootMass << mFoliageMass << mLeafArea;

}

void Tree::setup()

{

    if (mDbh<=0 || mHeight<=0)

        return;

    // check stamp

    Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL");

    mStamp = species()->stamp(mDbh, mHeight);

    mFoliageMass = species()->biomassFoliage(mDbh);

    mCoarseRootMass = species()->biomassRoot(mDbh); // coarse root (allometry)

    mFineRootMass = mFoliageMass * species()->finerootFoliageRatio(); //  fine root (size defined  by finerootFoliageRatio)

    mWoodyMass = species()->biomassWoody(mDbh);

    // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg]

    mLeafArea = mFoliageMass * species()->specificLeafArea();

    mOpacity = 1. - exp(- Model::settings().lightExtinctionCoefficientOpacity * mLeafArea / mStamp->crownArea());

    mNPPReserve = (1+species()->finerootFoliageRatio())*mFoliageMass; // initial value

    mDbhDelta = 0.1; // initial value: used in growth() to estimate diameter increment

    // initial value for tree aging...

    mRU->addTreeAging(mLeafArea,mSpecies->aging(mHeight, mAge));

}

void Tree::setAge(const int age, const float treeheight)

{

    mAge = age;

    if (age==0) {

        // estimate age using the tree height

        mAge = mSpecies->estimateAge(treeheight);

    }

}

//////////////////////////////////////////////////

////  Light functions (Pattern-stuff)

//////////////////////////////////////////////////

#define NOFULLDBG

//#define NOFULLOPT

void Tree::applyLIP()

{

    if (!mStamp)

        return;

    Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0);

    QPoint pos = mPositionIndex;

    int offset = mStamp->offset();

    pos-=QPoint(offset, offset);

    float local_dom; // height of Z* on the current position

    int x,y;

    float value, z, z_zstar;

    int gr_stamp = mStamp->size();

    int grid_x, grid_y;

    float *grid_value_ptr;

    if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) {

        // this should not happen because of the buffer

        return;

    }

    grid_y = pos.y();

    for (y=0;y<gr_stamp; ++y) {

        grid_value_ptr = mGrid->ptr(pos.x(), grid_y);

        grid_x = pos.x();

        for (x=0;x<gr_stamp;++x) {

            // suppose there is no stamping outside

            local_dom = mHeightGrid->valueAtIndex(grid_x/cPxPerHeight, grid_y/cPxPerHeight).height;

            z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45° line)

            z_zstar = (z>=local_dom)?1.f:z/local_dom;

            value = (*mStamp)(x,y); // stampvalue

            value = 1. - value*mOpacity * z_zstar; // calculated value

            // old: value = 1. - value*mOpacity / local_dom; // calculated value

            value = qMax(value, 0.02f); // limit value

            *grid_value_ptr++ *= value;

            grid_x++;

        }

        grid_y++;

    }

    m_statPrint++; // count # of stamp applications...

}

/// helper function for gluing the edges together

/// index: index at grid

/// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50)

/// buffer: size of buffer around simulation area (in pixels)

inline int torusIndex(int index, int count, int buffer, int ru_index)

{

    return buffer + ru_index + (index-buffer+count)%count;

}

/** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area.

  */

void Tree::applyLIP_torus()

{

    if (!mStamp)

        return;

    Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0);

    int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer

    QPoint pos = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU  + bufferOffset,

                        (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha

    QPoint ru_offset = QPoint(mPositionIndex.x() - pos.x(), mPositionIndex.y() - pos.y()); // offset of the corner of the resource index

    int offset = mStamp->offset();

    pos-=QPoint(offset, offset);

    float local_dom; // height of Z* on the current position

    int x,y;

    float value;

    int gr_stamp = mStamp->size();

    int grid_x, grid_y;

    float *grid_value;

    if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) {

        // todo: in this case we should use another algorithm!!! necessary????

        return;

    }

    float z, z_zstar;

    int xt, yt; // wraparound coordinates

    for (y=0;y<gr_stamp; ++y) {

        grid_y = pos.y() + y;

        yt = torusIndex(grid_y, cPxPerRU,bufferOffset, ru_offset.y()); // 50 cells per 100m

        for (x=0;x<gr_stamp;++x) {

            // suppose there is no stamping outside

            grid_x = pos.x() + x;

            xt = torusIndex(grid_x,cPxPerRU,bufferOffset, ru_offset.x());

            local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight,yt/cPxPerHeight).height;

            z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45° line)

            z_zstar = (z>=local_dom)?1.f:z/local_dom;

            value = (*mStamp)(x,y); // stampvalue

            value = 1. - value*mOpacity * z_zstar; // calculated value

            // old: value = 1. - value*mOpacity / local_dom; // calculated value

            value = qMax(value, 0.02f); // limit value

            grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates

            *grid_value *= value;

        }

    }

    m_statPrint++; // count # of stamp applications...

}

/** heightGrid()

  This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid.

*/

void Tree::heightGrid()

{

    QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid

    // count trees that are on height-grid cells (used for stockable area)

    mHeightGrid->valueAtIndex(p).increaseCount();

    if (mHeight > mHeightGrid->valueAtIndex(p).height)

        mHeightGrid->valueAtIndex(p).height=mHeight;

    int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5

    int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge

    int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge

    if (index_eastwest - r < 0) { // east

        mHeightGrid->valueAtIndex(p.x()-1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()-1, p.y()).height,mHeight);

    }

    if (index_eastwest + r >= cPxPerHeight) {  // west

        mHeightGrid->valueAtIndex(p.x()+1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()+1, p.y()).height,mHeight);

    }

    if (index_northsouth - r < 0) {  // south

        mHeightGrid->valueAtIndex(p.x(), p.y()-1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()-1).height,mHeight);

    }

    if (index_northsouth + r >= cPxPerHeight) {  // north

        mHeightGrid->valueAtIndex(p.x(), p.y()+1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()+1).height,mHeight);

    }

    // without spread of the height grid

//    // height of Z*

//    const float cellsize = mHeightGrid->cellsize();

//

//    int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge

//    int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge

//    int dist[9];

//    dist[3] = index_northsouth * 2 + 1; // south

//    dist[1] = index_eastwest * 2 + 1; // west

//    dist[5] = 10 - dist[3]; // north

//    dist[7] = 10 - dist[1]; // east

//    dist[8] = qMax(dist[5], dist[7]); // north-east

//    dist[6] = qMax(dist[3], dist[7]); // south-east

//    dist[0] = qMax(dist[3], dist[1]); // south-west

//    dist[2] = qMax(dist[5], dist[1]); // north-west

//    dist[4] = 0; // center cell

//    /* the scheme of indices is as follows:  if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:

//       index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.

//        e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2

//    */

//

//

//    int ringcount = int(floor(mHeight / cellsize)) + 1;

//    int ix, iy;

//    int ring;

//    float hdom;

//

//    for (ix=-ringcount;ix<=ringcount;ix++)

//        for (iy=-ringcount; iy<=+ringcount; iy++) {

//        ring = qMax(abs(ix), abs(iy));

//        QPoint pos(ix+p.x(), iy+p.y());

//        if (mHeightGrid->isIndexValid(pos)) {

//            float &rHGrid = mHeightGrid->valueAtIndex(pos).height;

//            if (rHGrid > mHeight) // skip calculation if grid is higher than tree

//                continue;

//            int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)

//            hdom = mHeight - dist[direction];

//            if (ring>1)

//                hdom -= (ring-1)*10;

//

//            rHGrid = qMax(rHGrid, hdom); // write value

//        } // is valid

//    } // for (y)

}

void Tree::heightGrid_torus()

{

    // height of Z*

    QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid

    int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer (i.e.: size of buffer in height-pixels)

    p.setX((p.x()-bufferOffset)%10 + bufferOffset); // 10: 10 x 10m pixeln in 100m

    p.setY((p.y()-bufferOffset)%10 + bufferOffset);

    // torus coordinates: ru_offset = coords of lower left corner of 1ha patch

    QPoint ru_offset =QPoint(mPositionIndex.x()/cPxPerHeight - p.x(), mPositionIndex.y()/cPxPerHeight - p.y());

    // count trees that are on height-grid cells (used for stockable area)

    HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()),

                                                   torusIndex(p.y(),10,bufferOffset,ru_offset.y()));

    v.increaseCount();

    v.height = qMax(v.height, mHeight);

    int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5

    int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge

    int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge

    if (index_eastwest - r < 0) { // east

        HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()-1,10,bufferOffset,ru_offset.x()),

                                                       torusIndex(p.y(),10,bufferOffset,ru_offset.y()));

        v.height = qMax(v.height, mHeight);

    }

    if (index_eastwest + r >= cPxPerHeight) {  // west

        HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()+1,10,bufferOffset,ru_offset.x()),

                                                       torusIndex(p.y(),10,bufferOffset,ru_offset.y()));

        v.height = qMax(v.height, mHeight);

    }

    if (index_northsouth - r < 0) {  // south

        HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()),

                                                       torusIndex(p.y()-1,10,bufferOffset,ru_offset.y()));

        v.height = qMax(v.height, mHeight);

    }

    if (index_northsouth + r >= cPxPerHeight) {  // north

        HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()),

                                                       torusIndex(p.y()+1,10,bufferOffset,ru_offset.y()));

        v.height = qMax(v.height, mHeight);

    }

//    int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge

//    int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge

//    int dist[9];

//    dist[3] = index_northsouth * 2 + 1; // south

//    dist[1] = index_eastwest * 2 + 1; // west

//    dist[5] = 10 - dist[3]; // north

//    dist[7] = 10 - dist[1]; // east

//    dist[8] = qMax(dist[5], dist[7]); // north-east

//    dist[6] = qMax(dist[3], dist[7]); // south-east

//    dist[0] = qMax(dist[3], dist[1]); // south-west

//    dist[2] = qMax(dist[5], dist[1]); // north-west

//    dist[4] = 0; // center cell

//    /* the scheme of indices is as follows:  if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:

//       index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.

//        e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2

//    */

//

//

//    int ringcount = int(floor(mHeight / cellsize)) + 1;

//    int ix, iy;

//    int ring;

//    float hdom;

//    for (ix=-ringcount;ix<=ringcount;ix++)

//        for (iy=-ringcount; iy<=+ringcount; iy++) {

//        ring = qMax(abs(ix), abs(iy));

//        QPoint pos(ix+p.x(), iy+p.y());

//        QPoint p_torus(torusIndex(pos.x(),10,bufferOffset,ru_offset.x()),

//                       torusIndex(pos.y(),10,bufferOffset,ru_offset.y()));

//        if (mHeightGrid->isIndexValid(p_torus)) {

//            float &rHGrid = mHeightGrid->valueAtIndex(p_torus.x(),p_torus.y()).height;

//            if (rHGrid > mHeight) // skip calculation if grid is higher than tree

//                continue;

//            int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)

//            hdom = mHeight - dist[direction];

//            if (ring>1)

//                hdom -= (ring-1)*10;

//

//            rHGrid = qMax(rHGrid, hdom); // write value

//        } // is valid

//    } // for (y)

}

void Tree::readLIF()

{

    if (!mStamp)

        return;

    const Stamp *reader = mStamp->reader();

    if (!reader)

        return;

    QPoint pos_reader = mPositionIndex;

    int offset_reader = reader->offset();

    int offset_writer = mStamp->offset();

    int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells

    pos_reader-=QPoint(offset_reader, offset_reader);

    float local_dom;

    int x,y;

    double sum=0.;

    double value, own_value;

    float *grid_value;

    float z, z_zstar;

    int reader_size = reader->size();

    int rx = pos_reader.x();

    int ry = pos_reader.y();

    for (y=0;y<reader_size; ++y, ++ry) {

        grid_value = mGrid->ptr(rx, ry);

        for (x=0;x<reader_size;++x) {

            local_dom = mHeightGrid->valueAtIndex((rx+x)/cPxPerHeight, ry/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not

            z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45° line)

            z_zstar = (z>=local_dom)?1.f:z/local_dom;

            own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar;

            // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;

            own_value = qMax(own_value, 0.02);

            value =  *grid_value++ / own_value; // remove impact of focal tree

            //qDebug() << x << y << local_dom << z << z_zstar << own_value << value << *(grid_value-1) << (*reader)(x,y) << mStamp->offsetValue(x,y,d_offset);

            //if (value>0.)

            sum += value * (*reader)(x,y);

        }

    }

    mLRI = sum;

    // LRI correction...

    double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height;

    if (hrel<1.)

        mLRI = species()->speciesSet()->LRIcorrection(mLRI, hrel);

    // read dominance field...

    // this applies only if some trees are potentially *higher* than the dominant height grid

    //if (dom_height < m_Height) {

        // if tree is higher than Z*, the dominance height

        // a part of the crown is in "full light".

    //    m_statAboveZ++;

    //    mImpact = 1. - (1. - mImpact)*dom_height/m_Height;

    //}

    if (mLRI > 1.)

        mLRI = 1.;

    // Finally, add LRI of this Tree to the ResourceUnit!

    mRU->addWLA(mLeafArea, mLRI);

    //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;

    //mRU->addWLA(mLRI*mLeafArea, mLeafArea);

}

/// Torus version of read stamp (glued edges)

void Tree::readLIF_torus()

{

    if (!mStamp)

        return;

    const Stamp *reader = mStamp->reader();

    if (!reader)

        return;

    int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer

    QPoint pos_reader = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset,

                               (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha

    QPoint ru_offset = QPoint(mPositionIndex.x() - pos_reader.x(), mPositionIndex.y() - pos_reader.y()); // offset of the corner of the resource index

    int offset_reader = reader->offset();

    int offset_writer = mStamp->offset();

    int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells

    pos_reader-=QPoint(offset_reader, offset_reader);

    float local_dom;

    int x,y;

    double sum=0.;

    double value, own_value;

    float *grid_value;

    float z, z_zstar;

    int reader_size = reader->size();

    int rx = pos_reader.x();

    int ry = pos_reader.y();

    int xt, yt; // wrapped coords

    for (y=0;y<reader_size; ++y) {

        yt = torusIndex(ry+y,cPxPerRU, bufferOffset, ru_offset.y());

        for (x=0;x<reader_size;++x) {

            xt = torusIndex(rx+x,cPxPerRU, bufferOffset, ru_offset.x());

            grid_value = mGrid->ptr(xt,yt);

            local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight, yt/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not

            z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45° line)

            z_zstar = (z>=local_dom)?1.f:z/local_dom;

            own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar;

            // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;

            own_value = qMax(own_value, 0.02);

            value =  *grid_value++ / own_value; // remove impact of focal tree

            // debug for one tree in HJA

            //if (id()==178020)

            //    qDebug() << x << y << xt << yt << *grid_value << local_dom << own_value << value << (*reader)(x,y);

            //if (_isnan(value))

            //    qDebug() << "isnan" << id();

            if (value * (*reader)(x,y)>1.)

                qDebug() << "LIFTorus: value>1.";

            sum += value * (*reader)(x,y);

            //} // isIndexValid

        }

    }

    mLRI = sum;

    // LRI correction...

    double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height;

    if (hrel<1.)

        mLRI = species()->speciesSet()->LRIcorrection(mLRI, hrel);

    if (_isnan(mLRI)) {

        qDebug() << "LRI invalid (nan)!" << id();

        mLRI=0.;

        //qDebug() << reader->dump();

    }

    if (mLRI > 1.)

        mLRI = 1.;

    //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;

    // Finally, add LRI of this Tree to the ResourceUnit!

    mRU->addWLA(mLeafArea, mLRI);

}

void Tree::resetStatistics()

{

    m_statPrint=0;

    m_statCreated=0;

    m_statAboveZ=0;

    m_nextId=1;

}

void Tree::calcLightResponse()

{

    // calculate a light response from lri:

    // http://iland.boku.ac.at/individual+tree+light+availability

    double lri = limit(mLRI * mRU->LRImodifier(), 0., 1.);

    mLightResponse = mSpecies->lightResponse(lri);

    mRU->addLR(mLeafArea, mLightResponse);

}

//////////////////////////////////////////////////

////  Growth Functions

//////////////////////////////////////////////////

/** grow() is the main function of the yearly tree growth.

  The main steps are:

  - Production of GPP/NPP   @sa http://iland.boku.ac.at/primary+production http://iland.boku.ac.at/individual+tree+light+availability

  - Partitioning of NPP to biomass compartments of the tree @sa http://iland.boku.ac.at/allocation

  - Growth of the stem http://iland.boku.ac.at/stem+growth (???)

  Further activties: * the age of the tree is increased

                     * the mortality sub routine is executed

                     * seeds are produced */

void Tree::grow()

{

    TreeGrowthData d;

    mAge++; // increase age

    // step 1: get "interception area" of the tree individual [m2]

    // the sum of all area of all trees of a unit equal the total stocked area * interception_factor(Beer-Lambert)

    double effective_area = mRU->interceptedArea(mLeafArea, mLightResponse);

    // step 2: calculate GPP of the tree based

    // (1) get the amount of GPP for a "unit area" of the tree species

    double raw_gpp_per_area = mRU->resourceUnitSpecies(species()).prod3PG().GPPperArea();

    // (2) GPP (without aging-effect) in kg Biomass / year

    double raw_gpp = raw_gpp_per_area * effective_area;

    // apply aging according to the state of the individuum

    const double aging_factor = mSpecies->aging(mHeight, mAge);

    mRU->addTreeAging(mLeafArea, aging_factor);

    double gpp = raw_gpp * aging_factor; //

    d.NPP = gpp * 0.47; // respiration loss, cf. Waring et al 1998.

    //DBGMODE(

        if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) {

            DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP);

            dumpList(out); // add tree headers

            out << mLRI * mRU->LRImodifier() << mLightResponse << effective_area << raw_gpp << gpp << d.NPP << aging_factor;

        }

    //); // DBGMODE()

    if (d.NPP>0.)

        partitioning(d); // split npp to compartments and grow (diameter, height)

    // mortality

    if (Model::settings().mortalityEnabled)

        mortality(d);

    mStressIndex = d.stress_index;

    if (!isDead())

        mRU->resourceUnitSpecies(species()).statistics().add(this, &d);

    // regeneration

    mSpecies->seedProduction(mAge, mHeight, mPositionIndex);

}

/** partitioning of this years assimilates (NPP) to biomass compartments.

  Conceptionally, the algorithm is based on Duursma, 2007.

  @sa http://iland.boku.ac.at/allocation */

inline void Tree::partitioning(TreeGrowthData &d)

{

    if (isDebugging())

        enableDebugging(true);

    double npp = d.NPP;

    // add content of reserve pool

    npp += mNPPReserve;

    const double foliage_mass_allo = species()->biomassFoliage(mDbh);

    const double reserve_size = foliage_mass_allo * (1. + mSpecies->finerootFoliageRatio());

    double refill_reserve = qMin(reserve_size, (1. + mSpecies->finerootFoliageRatio())*mFoliageMass); // not always try to refill reserve 100%

    double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta)

    ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species());

    // turnover rates

    const double &to_fol = species()->turnoverLeaf();

    const double &to_root = species()->turnoverRoot();

    // the turnover rate of wood depends on the size of the reserve pool:

    double to_wood = refill_reserve / (mWoodyMass + refill_reserve);

    apct_root = rus.prod3PG().rootFraction();

    d.NPP_above = d.NPP * (1. - apct_root); // aboveground: total NPP - fraction to roots

    double b_wf = species()->allometricRatio_wf(); // ratio of allometric exponents (b_woody / b_foliage)

    // Duursma 2007, Eq. (20)

    apct_wood = (foliage_mass_allo*to_wood/npp + b_wf*(1.-apct_root) - b_wf*foliage_mass_allo*to_fol/npp) / ( foliage_mass_allo/mWoodyMass + b_wf );

    if (apct_wood<0)

        apct_wood = 0.;

    apct_foliage = 1. - apct_root - apct_wood;

    //DBGMODE(

            if (apct_foliage<0 || apct_wood<0)

                qDebug() << "transfer to foliage or wood < 0";

             if (npp<0)

                 qDebug() << "NPP < 0";

         //   );

    // Change of biomass compartments

    double sen_root = mFineRootMass * to_root;

    double sen_foliage = mFoliageMass * to_fol;

    if (rus.snag())

        rus.snag()->addTurnoverLitter(this, sen_foliage, sen_root);

    // Roots

    // http://iland.boku.ac.at/allocation#belowground_NPP

    mFineRootMass -= sen_root; // reduce only fine root pool

    double delta_root = apct_root * npp;

    // 1st, refill the fine root pool

    double fineroot_miss = mFoliageMass * mSpecies->finerootFoliageRatio() - mFineRootMass;

    if (fineroot_miss>0.){

        double delta_fineroot = qMin(fineroot_miss, delta_root);

        mFineRootMass += delta_fineroot;

        delta_root -= delta_fineroot;

    }

    // 2nd, the rest of NPP allocated to roots go to coarse roots

    mCoarseRootMass += delta_root;

    // Foliage

    double delta_foliage = apct_foliage * npp - sen_foliage;

    mFoliageMass += delta_foliage;

    if (_isnan(mFoliageMass))

        qDebug() << "foliage mass invalid!";

    if (mFoliageMass<0.) mFoliageMass=0.; // limit to zero

    mLeafArea = mFoliageMass * species()->specificLeafArea(); // update leaf area

    // stress index: different varaints at denominator: to_fol*foliage_mass = leafmass to rebuild,

    // foliage_mass_allo: simply higher chance for stress

    // note: npp = NPP + reserve (see above)

    d.stress_index =qMax(1. - (npp) / ( to_fol*foliage_mass_allo + to_root*foliage_mass_allo*species()->finerootFoliageRatio() + reserve_size), 0.);

    // Woody compartments

    // see also: http://iland.boku.ac.at/allocation#reserve_and_allocation_to_stem_growth

    // (1) transfer to reserve pool

    double gross_woody = apct_wood * npp;

    double to_reserve = qMin(reserve_size, gross_woody);