WebSVN - public iLand - Diff - Rev 639 and 664


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/tree.cpp':
#include "global.h"
#include "math.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();
}

float Tree::biomassBranch() const
{
    return mSpecies->biomassBranch(mDbh);
}

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);
    if (!mStamp) {
        throw IException("Tree::setup() with invalid stamp!");
    }

    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

}

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.); // Eq. (3)
    mLightResponse = mSpecies->lightResponse(lri); // Eq. (4)
    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 * cAutotrophicRespiration; // respiration loss (0.47), 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)
{
    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 (ru()->snag())
        ru()->snag()->addTurnoverLitter(this->species(), 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
    double max_coarse_root = species()->biomassRoot(mDbh);
    mCoarseRootMass += delta_root;
    if (mCoarseRootMass > max_coarse_root) {
        // if the coarse root pool exceeds the value given by the allometry, then the
        // surplus is accounted as turnover
        if (ru()->snag())
            ru()->snag()->addTurnoverWood(species(), mCoarseRootMass-max_coarse_root);

        mCoarseRootMass = max_coarse_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);
    mNPPReserve = to_reserve;
    double net_woody = gross_woody - to_reserve;
    double net_stem = 0.;
    mDbhDelta = 0.;


    if (net_woody > 0.) {
        // (2) calculate part of increment that is dedicated to the stem (which is a function of diameter)
        net_stem = net_woody * species()->allometricFractionStem(mDbh);
        d.NPP_stem = net_stem;
        mWoodyMass += net_woody;
        //  (3) growth of diameter and height baseed on net stem increment
        grow_diameter(d);
    }

    //DBGMODE(
     if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreePartition)
         && isDebugging() ) {
            DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreePartition);
            dumpList(out); // add tree headers
            out << npp << apct_foliage << apct_wood << apct_root
                    << delta_foliage << net_woody << delta_root << mNPPReserve << net_stem << d.stress_index;
     }

    //); // DBGMODE()
    DBGMODE(
      if (mWoodyMass<0. || mWoodyMass>50000 || mFoliageMass<0. || mFoliageMass>2000. || mCoarseRootMass<0. || mCoarseRootMass>30000
         || mNPPReserve>4000.) {
         qDebug() << "Tree:partitioning: invalid or unlikely pools.";
         qDebug() << GlobalSettings::instance()->debugListCaptions(GlobalSettings::DebugOutputs(0));
         DebugList dbg; dumpList(dbg);
         qDebug() << dbg;
     } );

    /*DBG_IF_X(mId == 1 , "Tree::partitioning", "dump", dump()
             + QString("npp %1 npp_reserve %9 sen_fol %2 sen_stem %3 sen_root %4 net_fol %5 net_stem %6 net_root %7 to_reserve %8")
               .arg(npp).arg(senescence_foliage).arg(senescence_stem).arg(senescence_root)
               .arg(net_foliage).arg(net_stem).arg(net_root).arg(to_reserve).arg(mNPPReserve) );*/

}


/** Determination of diamter and height growth based on increment of the stem mass (@p net_stem_npp).
    Refer to XXX for equations and variables.
    This function updates the dbh and height of the tree.
    The equations are based on dbh in meters! */
inline void Tree::grow_diameter(TreeGrowthData &d)
{
    // determine dh-ratio of increment
    // height increment is a function of light competition:
    double hd_growth = relative_height_growth(); // hd of height growth
    double d_m = mDbh / 100.; // current diameter in [m]
    double net_stem_npp = d.NPP_stem;

    const double d_delta_m = mDbhDelta / 100.; // increment of last year in [m]

    const double mass_factor = species()->volumeFactor() * species()->density();
    double stem_mass = mass_factor * d_m*d_m * mHeight; // result: kg, dbh[cm], h[meter]

    // factor is in diameter increment per NPP [m/kg]
    double factor_diameter = 1. / (  mass_factor * (d_m + d_delta_m)*(d_m + d_delta_m) * ( 2. * mHeight/d_m + hd_growth) );
    double delta_d_estimate = factor_diameter * net_stem_npp; // estimated dbh-inc using last years increment

    // using that dbh-increment we estimate a stem-mass-increment and the residual (Eq. 9)
    double stem_estimate = mass_factor * (d_m + delta_d_estimate)*(d_m + delta_d_estimate)*(mHeight + delta_d_estimate*hd_growth);
    double stem_residual = stem_estimate - (stem_mass + net_stem_npp);

    // the final increment is then:
    double d_increment = factor_diameter * (net_stem_npp - stem_residual); // Eq. (11)
    double res_final  = 0.;
    if (fabs(stem_residual) > 1.) {

        // calculate final residual in stem
        res_final = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp));
        if (fabs(res_final)>1.) {
            // for large errors in stem biomass due to errors in diameter increment (> 1kg), we solve the increment iteratively.
            // first, increase increment with constant step until we overestimate the first time
            // then,
            d_increment = 0.02; // start with 2cm increment
            bool reached_error = false;
            double step=0.01; // step-width 1cm
            double est_stem;
            do {
                est_stem = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth); // estimate with current increment
                stem_residual = est_stem - (stem_mass + net_stem_npp);

                if (fabs(stem_residual) <1.) // finished, if stem residual below 1kg
                    break;
                if (stem_residual > 0.) {
                    d_increment -= step;
                    reached_error=true;
                } else {
                    d_increment += step;
                }
                if (reached_error)
                    step /= 2.;
            } while (step>0.00001); // continue until diameter "accuracy" falls below 1/100mm
        }
    }

    if (d_increment<0.f)
        qDebug() << "Tree::grow_diameter: d_inc < 0.";
    DBG_IF_X(d_increment<0. || d_increment>0.1, "Tree::grow_dimater", "increment out of range.", dump()
             + QString("\nhdz %1 factor_diameter %2 stem_residual %3 delta_d_estimate %4 d_increment %5 final residual(kg) %6")
               .arg(hd_growth).arg(factor_diameter).arg(stem_residual).arg(delta_d_estimate).arg(d_increment)
               .arg( mass_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)) ));

    //DBGMODE(
        // do not calculate res_final twice if already done
        DBG_IF_X( (res_final==0.?fabs(mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp))):res_final) > 1, "Tree::grow_diameter", "final residual stem estimate > 1kg", dump());
        DBG_IF_X(d_increment > 10. || d_increment*hd_growth >10., "Tree::grow_diameter", "growth out of bound:",QString("d-increment %1 h-increment %2 ").arg(d_increment).arg(d_increment*hd_growth/100.) + dump());

        if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeGrowth) && isDebugging() ) {
            DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeGrowth);
            dumpList(out); // add tree headers
            out << net_stem_npp << stem_mass << hd_growth << factor_diameter << delta_d_estimate*100 << d_increment*100;
        }

    //); // DBGMODE()

    d_increment = qMax(d_increment, 0.);

    // update state variables
    mDbh += d_increment*100; // convert from [m] to [cm]
    mDbhDelta = d_increment*100; // save for next year's growth
    mHeight += d_increment * hd_growth;

    // update state of LIP stamp and opacity
    mStamp = species()->stamp(mDbh, mHeight); // get new stamp for updated dimensions
    // calculate the CrownFactor which reflects the opacity of the crown
    const double k=Model::settings().lightExtinctionCoefficientOpacity;
    mOpacity = 1. - exp(-k * mLeafArea / mStamp->crownArea());

}


/// return the HD ratio of this year's increment based on the light status.
inline double Tree::relative_height_growth()
{
    double hd_low, hd_high;
    mSpecies->hdRange(mDbh, hd_low, hd_high);

    DBG_IF_X(hd_low>hd_high, "Tree::relative_height_growth", "hd low higher dann hd_high for ", dump());
    DBG_IF_X(hd_low < 20 || hd_high>250, "Tree::relative_height_growth", "hd out of range ", dump() + QString(" hd-low: %1 hd-high: %2").arg(hd_low).arg(hd_high));

    // scale according to LRI: if receiving much light (LRI=1), the result is hd_low (for open grown trees)
    // use the corrected LRI (see tracker#11)
    double lri = limit(mLRI * mRU->LRImodifier(),0.,1.);
    double hd_ratio = hd_high - (hd_high-hd_low)*lri;
    return hd_ratio;
}

/** This function is called if a tree dies.
  @sa ResourceUnit::cleanTreeList(), remove() */
void Tree::die(TreeGrowthData *d)
{
    setFlag(Tree::TreeDead, true); // set flag that tree is dead
    mRU->treeDied();
    ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species());
    rus.statisticsDead().add(this, d); // add tree to statistics
    if (ru()->snag())
        ru()->snag()->addMortality(this);
}

void Tree::remove(double removeFoliage, double removeBranch, double removeStem )
{
    setFlag(Tree::TreeDead, true); // set flag that tree is dead
    mRU->treeDied();
    ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species());
    rus.statisticsMgmt().add(this, 0);
    if (ru()->snag())
        ru()->snag()->addHarvest(this, removeStem, removeBranch, removeFoliage);
}

void Tree::mortality(TreeGrowthData &d)
{
    // death if leaf area is 0
    if (mFoliageMass<0.00001)
        die();

    double p_death,  p_stress, p_intrinsic;
    p_intrinsic = species()->deathProb_intrinsic();
    p_stress = species()->deathProb_stress(d.stress_index);
    p_death =  p_intrinsic + p_stress;
    double p = drandom(ru()->randomGenerator()); //0..1
    if (p<p_death) {
        // die...
        die();
    }
}

//////////////////////////////////////////////////
////  value functions
//////////////////////////////////////////////////

double Tree::volume() const
{
    /// @see Species::volumeFactor() for details
    const double volume_factor = species()->volumeFactor();
    const double volume =  volume_factor * mDbh*mDbh*mHeight * 0.0001; // dbh in cm: cm/100 * cm/100 = cm*cm * 0.0001 = m2
    return volume;
}

/// return the basal area in m2
double Tree::basalArea() const
{
    // A = r^2 * pi = d/2*pi; from cm->m: d/200
    const double b = (mDbh/200.)*(mDbh/200.)*M_PI;
    return b;
}
 
Subversion Repositories public iLand

(root)/src/core/tree.cpp @ 1222 - Rev 639 → 664
