WebSVN - public iLand - Diff - Rev 167 and 169


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

#include "grid.h"

#include "stamp.h"
#include "species.h"
#include "ressourceunit.h"
#include "model.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;

/// internal data structure which is passed between function
struct TreeGrowthData
{
    double NPP; ///< total NPP
    double NPP_stem;  ///< NPP used for growth of stem (dbh,h)
    double stress_index; ///< stress index used for mortality calculation
};

/** 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=mRootMass=mLeafArea=0.;
    mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.;
    mId = m_nextId++;
    m_statCreated++;
}

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 << mRootMass << 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);
    mRootMass = species()->biomassRoot(mDbh) + mFoliageMass; // coarse root (allometry) + fine root (estimated size: foliage)
    mWoodyMass = species()->biomassWoody(mDbh);

    // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg]
    mLeafArea = mFoliageMass * species()->specificLeafArea();
    mOpacity = 1. - exp(-0.5 * mLeafArea / mStamp->crownArea());
    mNPPReserve = 2*mFoliageMass; // initial value
    mDbhDelta = 0.1; // initial value: used in growth() to estimate diameter increment
}

//////////////////////////////////////////////////
////  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;
    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!!!
        return;
    }

    for (y=0;y<gr_stamp; ++y) {
        grid_y = pos.y() + y;
        grid_value = mGrid->ptr(pos.x(), grid_y);
        for (x=0;x<gr_stamp;++x) {
            // suppose there is no stamping outside
            grid_x = pos.x() + x;

            local_dom = mHeightGrid->valueAtIndex(grid_x/5, grid_y/5).height;
            value = (*mStamp)(x,y); // stampvalue
            value = 1. - value*mOpacity / local_dom; // calculated value
            value = qMax(value, 0.02f); // limit value

            *grid_value++ *= value;
        }
    }

    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)
int torusIndex(int index, int count, int buffer)
{
    return buffer + (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);

    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;
    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;
    }
    int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer
    int xt, yt; // wraparound coordinates
    for (y=0;y<gr_stamp; ++y) {
        grid_y = pos.y() + y;
        yt = torusIndex(grid_y, 50,bufferOffset); // 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,50,bufferOffset);

            local_dom = mHeightGrid->valueAtIndex(xt/5,yt/5).height;
            value = (*mStamp)(x,y); // stampvalue
            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()
{
    // height of Z*
    const float cellsize = mHeightGrid->cellsize();

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

    // count trees that are on height-grid cells (used for stockable area)
    mHeightGrid->valueAtIndex(p).count++;

    int index_eastwest = mPositionIndex.x() % 5; // 4: very west, 0 east edge
    int index_northsouth = mPositionIndex.y() % 5; // 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;
    QPoint pos;
    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::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

    QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer);
    pos_reader-=QPoint(offset_reader, offset_reader);
    QPoint p;

    //float dom_height = (*m_dominanceGrid)[m_Position];
    float local_dom;

    int x,y;
    double sum=0.;
    double value, own_value;
    float *grid_value;
    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) {

            //p = pos_reader + QPoint(x,y);
            //if (m_grid->isIndexValid(p)) {
            local_dom = mHeightGrid->valueAtIndex((rx+x)/5, ry/5).height; // ry: gets ++ed in outer loop, rx not
            //local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) );

            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
            //if (value>0.)
            sum += value * (*reader)(x,y);

            //} // isIndexValid
        }
    }
    mLRI = sum;
    // 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.;
    //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
    mRU->addWLA(mLRI*mLeafArea, mLeafArea);
}

void Tree::heightGrid_torus()
{
    // height of Z*
    const float cellsize = mHeightGrid->cellsize();

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

    // count trees that are on height-grid cells (used for stockable area)
    mHeightGrid->valueAtIndex(p).count++;

    int index_eastwest = mPositionIndex.x() % 5; // 4: very west, 0 east edge
    int index_northsouth = mPositionIndex.y() % 5; // 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;
    QPoint pos;
    float hdom;
    int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer
    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(torusIndex(pos.x(),10,bufferOffset), torusIndex(pos.y(),10,bufferOffset)).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)
}

/// Torus version of read stamp (glued edges)
void Tree::readLIF_torus()
{
    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

    QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer);
    pos_reader-=QPoint(offset_reader, offset_reader);
    QPoint p;

    //float dom_height = (*m_dominanceGrid)[m_Position];
    float local_dom;

    int x,y;
    double sum=0.;
    double value, own_value;
    float *grid_value;
    int reader_size = reader->size();
    int rx = pos_reader.x();
    int ry = pos_reader.y();
    int xt, yt; // wrapped coords
    int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer

    for (y=0;y<reader_size; ++y, ++ry) {
        grid_value = mGrid->ptr(rx, ry);
        for (x=0;x<reader_size;++x) {
            xt = torusIndex(rx+x,50, bufferOffset);
            yt = torusIndex(ry+y,50, bufferOffset);
            grid_value = mGrid->ptr(xt,yt);
            //p = pos_reader + QPoint(x,y);
            //if (m_grid->isIndexValid(p)) {
            local_dom = mHeightGrid->valueAtIndex(xt/5, yt/5).height; // ry: gets ++ed in outer loop, rx not
            //local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) );

            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
            //if (value>0.)
            sum += value * (*reader)(x,y);

            //} // isIndexValid
        }
    }
    mLRI = sum;
    // 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.;
    //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
    mRU->addWLA(mLRI*mLeafArea, mLeafArea);
}


void Tree::resetStatistics()
{
    m_statPrint=0;
    m_statCreated=0;
    m_statAboveZ=0;
    m_nextId=1;
}

//////////////////////////////////////////////////
////  Growth Functions
//////////////////////////////////////////////////


void Tree::grow()
{
    TreeGrowthData d;
    mAge++; // increase age
    // step 1: get radiation from ressource unit: radiation (MJ/tree/year) total intercepted radiation for this tree per year!
    double radiation = mRU->interceptedRadiation(mLeafArea, mLRI);
    // step 2: get fraction of PARutilized, i.e. fraction of intercepted rad that is utiliziable (per year)

    double raw_gpp_per_rad = mRU->ressourceUnitSpecies(species()).prod3PG().GPPperRad();
    // GPP (without aging-effect) [gC] / year -> kg Biomass /GPP (*0.001 *2)
    double raw_gpp = raw_gpp_per_rad * radiation * 0.001 * 2;

    // apply aging
 
    const double aging_factor = mSpecies->aging(mHeight, mAge);
    double gpp = raw_gpp * aging_factor;
    d.NPP = gpp * 0.47; // respiration loss, transformation kg

    DBGMODE(
        if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) {
            DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP);
            dumpList(out); // add tree headers
 
            out << radiation << raw_gpp << gpp << d.NPP << aging_factor;
        }
    ); // DBGMODE()


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

    mortality(d);

    mStressIndex = d.stress_index;
}


inline void Tree::partitioning(TreeGrowthData &d)
{
    if (isDebugging())
        enableDebugging(true);
    double npp = d.NPP;
    double harshness = mRU->ressourceUnitSpecies(species()).prod3PG().harshness();
    // add content of reserve pool
    npp += mNPPReserve;
    const double foliage_mass_allo = species()->biomassFoliage(mDbh);
    const double reserve_size = 2 * foliage_mass_allo;
    double refill_reserve = qMin(reserve_size, 2.*mFoliageMass); // not always try to refill reserve 100%

    double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta)
    // 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 = harshness;
    double b_wf = species()->allometricRatio_wf(); // ratio of allometric exponents... now fixed

    // 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 = mRootMass*to_root;
    double sen_foliage = mFoliageMass*to_fol;
    // Roots
    double delta_root = apct_root * npp - sen_root;
    mRootMass += delta_root;

    // Foliage
    double delta_foliage = apct_foliage * npp - sen_foliage;
    mFoliageMass += delta_foliage;
    if (mFoliageMass<0.) mFoliageMass=0.; // limit to zero

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

    // stress index
    d.stress_index =qMax(1. - (npp-reserve_size) / to_fol*foliage_mass_allo, 0.);

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

    ); // DBGMODE()
    //DBGMODE(
      if (mWoodyMass<0. || mWoodyMass>10000 || mFoliageMass<0. || mFoliageMass>1000. || mRootMass<0. || mRootMass>10000
         || mNPPReserve>2000.) {
         qDebug() << "Tree:partitioning: invalid 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)
    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(
        double res_final = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp));
        DBG_IF_X(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
    mOpacity = 1. - exp(-0.7 * 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)
    double hd_ratio = hd_high - (hd_high-hd_low)*mLRI;
    return hd_ratio;
}

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

    double p_death,  p_stress;
    //p_stress = d.stress_index * species()->deathProb_stress();
    if (d.stress_index>0)
        p_stress = species()->deathProb_stress();
    p_death = species()->deathProb_intrinsic() + p_stress;
    double p = random(); //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;
}
 

Rev 167	Rev 169
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/tree.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/tree.cpp':
2	#include "global.h"	2	#include "global.h"
3	#include "tree.h"	3	#include "tree.h"
4		4
5	#include "grid.h"	5	#include "grid.h"
6		6
7	#include "stamp.h"	7	#include "stamp.h"
8	#include "species.h"	8	#include "species.h"
9	#include "ressourceunit.h"	9	#include "ressourceunit.h"
10	#include "model.h"	10	#include "model.h"
11		11
12	// static varaibles	12	// static varaibles
13	FloatGrid *Tree::mGrid = 0;	13	FloatGrid *Tree::mGrid = 0;
14	HeightGrid *Tree::mHeightGrid = 0;	14	HeightGrid *Tree::mHeightGrid = 0;
15	int Tree::m_statPrint=0;	15	int Tree::m_statPrint=0;
16	int Tree::m_statAboveZ=0;	16	int Tree::m_statAboveZ=0;
17	int Tree::m_statCreated=0;	17	int Tree::m_statCreated=0;
18	int Tree::m_nextId=0;	18	int Tree::m_nextId=0;
19		19
20	/// internal data structure which is passed between function	20	/// internal data structure which is passed between function
21	struct TreeGrowthData	21	struct TreeGrowthData
22	{	22	{
23	double NPP; ///< total NPP	23	double NPP; ///< total NPP
24	double NPP_stem; ///< NPP used for growth of stem (dbh,h)	24	double NPP_stem; ///< NPP used for growth of stem (dbh,h)
25	double stress_index; ///< stress index used for mortality calculation	25	double stress_index; ///< stress index used for mortality calculation
26	};	26	};
27		27
28	/** get distance and direction between two points.	28	/** get distance and direction between two points.
29	returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */	29	returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */
30	float dist_and_direction(const QPointF &PStart, const QPointF &PEnd, float &rAngle)	30	float dist_and_direction(const QPointF &PStart, const QPointF &PEnd, float &rAngle)
31	{	31	{
32	float dx = PEnd.x() - PStart.x();	32	float dx = PEnd.x() - PStart.x();
33	float dy = PEnd.y() - PStart.y();	33	float dy = PEnd.y() - PStart.y();
34	float d = sqrt(dxdx + dydy);	34	float d = sqrt(dxdx + dydy);
35	// direction:	35	// direction:
36	rAngle = atan2(dx, dy);	36	rAngle = atan2(dx, dy);
37	return d;	37	return d;
38	}	38	}
39		39
40		40
41	// lifecycle	41	// lifecycle
42	Tree::Tree()	42	Tree::Tree()
43	{	43	{
44	mDbh = mHeight = 0;	44	mDbh = mHeight = 0;
45	mRU = 0; mSpecies = 0;	45	mRU = 0; mSpecies = 0;
46	mFlags = 0;	-
-		46	mFlags = mAge = 0;
47	mOpacity=mFoliageMass=mWoodyMass=mRootMass=mLeafArea=0.;	47	mOpacity=mFoliageMass=mWoodyMass=mRootMass=mLeafArea=0.;
48	mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.;	48	mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.;
49	mId = m_nextId++;	49	mId = m_nextId++;
50	m_statCreated++;	50	m_statCreated++;
51	}	51	}
52		52
53	void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid)	53	void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid)
54	{	54	{
55	mGrid = gridToStamp; mHeightGrid = dominanceGrid;	55	mGrid = gridToStamp; mHeightGrid = dominanceGrid;
56	}	56	}
57		57
58	/// dumps some core variables of a tree to a string.	58	/// dumps some core variables of a tree to a string.
59	QString Tree::dump()	59	QString Tree::dump()
60	{	60	{
61	QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8")	61	QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8")
62	.arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight)	62	.arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight)
63	.arg(position().x()).arg(position().y())	63	.arg(position().x()).arg(position().y())
64	.arg(mRU->index()).arg(mLRI);	64	.arg(mRU->index()).arg(mLRI);
65	return result;	65	return result;
66	}	66	}
67		67
68	void Tree::dumpList(DebugList &rTargetList)	68	void Tree::dumpList(DebugList &rTargetList)
69	{	69	{
70	rTargetList << mId << species()->id() << mDbh << mHeight << position().x() << position().y() << mRU->index() << mLRI	70	rTargetList << mId << species()->id() << mDbh << mHeight << position().x() << position().y() << mRU->index() << mLRI
71	<< mWoodyMass << mRootMass << mFoliageMass << mLeafArea;	71	<< mWoodyMass << mRootMass << mFoliageMass << mLeafArea;
72	}	72	}
73		73
74	void Tree::setup()	74	void Tree::setup()
75	{	75	{
76	if (mDbh<=0 \|\| mHeight<=0)	76	if (mDbh<=0 \|\| mHeight<=0)
77	return;	77	return;
78	// check stamp	78	// check stamp
79	Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL");	79	Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL");
80	mStamp = species()->stamp(mDbh, mHeight);	80	mStamp = species()->stamp(mDbh, mHeight);
81		81
82	mFoliageMass = species()->biomassFoliage(mDbh);	82	mFoliageMass = species()->biomassFoliage(mDbh);
83	mRootMass = species()->biomassRoot(mDbh) + mFoliageMass; // coarse root (allometry) + fine root (estimated size: foliage)	83	mRootMass = species()->biomassRoot(mDbh) + mFoliageMass; // coarse root (allometry) + fine root (estimated size: foliage)
84	mWoodyMass = species()->biomassWoody(mDbh);	84	mWoodyMass = species()->biomassWoody(mDbh);
85		85
86	// LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg]	86	// LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg]
87	mLeafArea = mFoliageMass * species()->specificLeafArea();	87	mLeafArea = mFoliageMass * species()->specificLeafArea();
88	mOpacity = 1. - exp(-0.5 * mLeafArea / mStamp->crownArea());	88	mOpacity = 1. - exp(-0.5 * mLeafArea / mStamp->crownArea());
89	mNPPReserve = 2*mFoliageMass; // initial value	89	mNPPReserve = 2*mFoliageMass; // initial value
90	mDbhDelta = 0.1; // initial value: used in growth() to estimate diameter increment	90	mDbhDelta = 0.1; // initial value: used in growth() to estimate diameter increment
91	}	91	}
92		92
93	//////////////////////////////////////////////////	93	//////////////////////////////////////////////////
94	//// Light functions (Pattern-stuff)	94	//// Light functions (Pattern-stuff)
95	//////////////////////////////////////////////////	95	//////////////////////////////////////////////////
96		96
97	#define NOFULLDBG	97	#define NOFULLDBG
98	//#define NOFULLOPT	98	//#define NOFULLOPT
99		99
100		100
101	void Tree::applyLIP()	101	void Tree::applyLIP()
102	{	102	{
103	if (!mStamp)	103	if (!mStamp)
104	return;	104	return;
105	Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0);	105	Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0);
106	QPoint pos = mPositionIndex;	106	QPoint pos = mPositionIndex;
107	int offset = mStamp->offset();	107	int offset = mStamp->offset();
108	pos-=QPoint(offset, offset);	108	pos-=QPoint(offset, offset);
109		109
110	float local_dom; // height of Z* on the current position	110	float local_dom; // height of Z* on the current position
111	int x,y;	111	int x,y;
112	float value;	112	float value;
113	int gr_stamp = mStamp->size();	113	int gr_stamp = mStamp->size();
114	int grid_x, grid_y;	114	int grid_x, grid_y;
115	float *grid_value;	115	float *grid_value;
116	if (!mGrid->isIndexValid(pos) \|\| !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) {	116	if (!mGrid->isIndexValid(pos) \|\| !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) {
117	// todo: in this case we should use another algorithm!!!	117	// todo: in this case we should use another algorithm!!!
118	return;	118	return;
119	}	119	}
120		120
121	for (y=0;y<gr_stamp; ++y) {	121	for (y=0;y<gr_stamp; ++y) {
122	grid_y = pos.y() + y;	122	grid_y = pos.y() + y;
123	grid_value = mGrid->ptr(pos.x(), grid_y);	123	grid_value = mGrid->ptr(pos.x(), grid_y);
124	for (x=0;x<gr_stamp;++x) {	124	for (x=0;x<gr_stamp;++x) {
125	// suppose there is no stamping outside	125	// suppose there is no stamping outside
126	grid_x = pos.x() + x;	126	grid_x = pos.x() + x;
127		127
128	local_dom = mHeightGrid->valueAtIndex(grid_x/5, grid_y/5).height;	128	local_dom = mHeightGrid->valueAtIndex(grid_x/5, grid_y/5).height;
129	value = (*mStamp)(x,y); // stampvalue	129	value = (*mStamp)(x,y); // stampvalue
130	value = 1. - value*mOpacity / local_dom; // calculated value	130	value = 1. - value*mOpacity / local_dom; // calculated value
131	value = qMax(value, 0.02f); // limit value	131	value = qMax(value, 0.02f); // limit value
132		132
133	grid_value++ = value;	133	grid_value++ = value;
134	}	134	}
135	}	135	}
136		136
137	m_statPrint++; // count # of stamp applications...	137	m_statPrint++; // count # of stamp applications...
138	}	138	}
139		139
140	/// helper function for gluing the edges together	140	/// helper function for gluing the edges together
141	/// index: index at grid	141	/// index: index at grid
142	/// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50)	142	/// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50)
143	/// buffer: size of buffer around simulation area (in pixels)	143	/// buffer: size of buffer around simulation area (in pixels)
144	int torusIndex(int index, int count, int buffer)	144	int torusIndex(int index, int count, int buffer)
145	{	145	{
146	return buffer + (index-buffer+count)%count;	146	return buffer + (index-buffer+count)%count;
147	}	147	}
148		148
149		149
150	/** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area.	150	/** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area.
151	*/	151	*/
152	void Tree::applyLIP_torus()	152	void Tree::applyLIP_torus()
153	{	153	{
154	if (!mStamp)	154	if (!mStamp)
155	return;	155	return;
156	Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0);	156	Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0);
157		157
158	QPoint pos = mPositionIndex;	158	QPoint pos = mPositionIndex;
159	int offset = mStamp->offset();	159	int offset = mStamp->offset();
160	pos-=QPoint(offset, offset);	160	pos-=QPoint(offset, offset);
161		161
162	float local_dom; // height of Z* on the current position	162	float local_dom; // height of Z* on the current position
163	int x,y;	163	int x,y;
164	float value;	164	float value;
165	int gr_stamp = mStamp->size();	165	int gr_stamp = mStamp->size();
166	int grid_x, grid_y;	166	int grid_x, grid_y;
167	float *grid_value;	167	float *grid_value;
168	if (!mGrid->isIndexValid(pos) \|\| !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) {	168	if (!mGrid->isIndexValid(pos) \|\| !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) {
169	// todo: in this case we should use another algorithm!!! necessary????	169	// todo: in this case we should use another algorithm!!! necessary????
170	return;	170	return;
171	}	171	}
172	int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer	172	int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer
173	int xt, yt; // wraparound coordinates	173	int xt, yt; // wraparound coordinates
174	for (y=0;y<gr_stamp; ++y) {	174	for (y=0;y<gr_stamp; ++y) {
175	grid_y = pos.y() + y;	175	grid_y = pos.y() + y;
176	yt = torusIndex(grid_y, 50,bufferOffset); // 50 cells per 100m	176	yt = torusIndex(grid_y, 50,bufferOffset); // 50 cells per 100m
177	for (x=0;x<gr_stamp;++x) {	177	for (x=0;x<gr_stamp;++x) {
178	// suppose there is no stamping outside	178	// suppose there is no stamping outside
179	grid_x = pos.x() + x;	179	grid_x = pos.x() + x;
180	xt = torusIndex(grid_x,50,bufferOffset);	180	xt = torusIndex(grid_x,50,bufferOffset);
181		181
182	local_dom = mHeightGrid->valueAtIndex(xt/5,yt/5).height;	182	local_dom = mHeightGrid->valueAtIndex(xt/5,yt/5).height;
183	value = (*mStamp)(x,y); // stampvalue	183	value = (*mStamp)(x,y); // stampvalue
184	value = 1. - value*mOpacity / local_dom; // calculated value	184	value = 1. - value*mOpacity / local_dom; // calculated value
185	value = qMax(value, 0.02f); // limit value	185	value = qMax(value, 0.02f); // limit value
186		186
187	grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates	187	grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates
188	grid_value = value;	188	grid_value = value;
189	}	189	}
190	}	190	}
191		191
192	m_statPrint++; // count # of stamp applications...	192	m_statPrint++; // count # of stamp applications...
193	}	193	}
194		194
195	/** heightGrid()	195	/** heightGrid()
196	This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid.	196	This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid.
197	*/	197	*/
198	void Tree::heightGrid()	198	void Tree::heightGrid()
199	{	199	{
200	// height of Z*	200	// height of Z*
201	const float cellsize = mHeightGrid->cellsize();	201	const float cellsize = mHeightGrid->cellsize();
202		202
203	QPoint p = QPoint(mPositionIndex.x()/5, mPositionIndex.y()/5); // pos of tree on height grid	203	QPoint p = QPoint(mPositionIndex.x()/5, mPositionIndex.y()/5); // pos of tree on height grid
204		204
205	// count trees that are on height-grid cells (used for stockable area)	205	// count trees that are on height-grid cells (used for stockable area)
206	mHeightGrid->valueAtIndex(p).count++;	206	mHeightGrid->valueAtIndex(p).count++;
207		207
208	int index_eastwest = mPositionIndex.x() % 5; // 4: very west, 0 east edge	208	int index_eastwest = mPositionIndex.x() % 5; // 4: very west, 0 east edge
209	int index_northsouth = mPositionIndex.y() % 5; // 4: northern edge, 0: southern edge	209	int index_northsouth = mPositionIndex.y() % 5; // 4: northern edge, 0: southern edge
210	int dist[9];	210	int dist[9];
211	dist[3] = index_northsouth * 2 + 1; // south	211	dist[3] = index_northsouth * 2 + 1; // south
212	dist[1] = index_eastwest * 2 + 1; // west	212	dist[1] = index_eastwest * 2 + 1; // west
213	dist[5] = 10 - dist[3]; // north	213	dist[5] = 10 - dist[3]; // north
214	dist[7] = 10 - dist[1]; // east	214	dist[7] = 10 - dist[1]; // east
215	dist[8] = qMax(dist[5], dist[7]); // north-east	215	dist[8] = qMax(dist[5], dist[7]); // north-east
216	dist[6] = qMax(dist[3], dist[7]); // south-east	216	dist[6] = qMax(dist[3], dist[7]); // south-east
217	dist[0] = qMax(dist[3], dist[1]); // south-west	217	dist[0] = qMax(dist[3], dist[1]); // south-west
218	dist[2] = qMax(dist[5], dist[1]); // north-west	218	dist[2] = qMax(dist[5], dist[1]); // north-west
219	dist[4] = 0; // center cell	219	dist[4] = 0; // center cell
220	/* 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:	220	/* 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:
221	index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.	221	index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
222	e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2	222	e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
223	*/	223	*/
224		224
225		225
226	int ringcount = int(floor(mHeight / cellsize)) + 1;	226	int ringcount = int(floor(mHeight / cellsize)) + 1;
227	int ix, iy;	227	int ix, iy;
228	int ring;	228	int ring;
229	QPoint pos;	229	QPoint pos;
230	float hdom;	230	float hdom;
231		231
232	for (ix=-ringcount;ix<=ringcount;ix++)	232	for (ix=-ringcount;ix<=ringcount;ix++)
233	for (iy=-ringcount; iy<=+ringcount; iy++) {	233	for (iy=-ringcount; iy<=+ringcount; iy++) {
234	ring = qMax(abs(ix), abs(iy));	234	ring = qMax(abs(ix), abs(iy));
235	QPoint pos(ix+p.x(), iy+p.y());	235	QPoint pos(ix+p.x(), iy+p.y());
236	if (mHeightGrid->isIndexValid(pos)) {	236	if (mHeightGrid->isIndexValid(pos)) {
237	float &rHGrid = mHeightGrid->valueAtIndex(pos).height;	237	float &rHGrid = mHeightGrid->valueAtIndex(pos).height;
238	if (rHGrid > mHeight) // skip calculation if grid is higher than tree	238	if (rHGrid > mHeight) // skip calculation if grid is higher than tree
239	continue;	239	continue;
240	int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)	240	int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
241	hdom = mHeight - dist[direction];	241	hdom = mHeight - dist[direction];
242	if (ring>1)	242	if (ring>1)
243	hdom -= (ring-1)*10;	243	hdom -= (ring-1)*10;
244		244
245	rHGrid = qMax(rHGrid, hdom); // write value	245	rHGrid = qMax(rHGrid, hdom); // write value
246	} // is valid	246	} // is valid
247	} // for (y)	247	} // for (y)
248	}	248	}
249		249
250		250
251		251
252	void Tree::readLIF()	252	void Tree::readLIF()
253	{	253	{
254	if (!mStamp)	254	if (!mStamp)
255	return;	255	return;
256	const Stamp *reader = mStamp->reader();	256	const Stamp *reader = mStamp->reader();
257	if (!reader)	257	if (!reader)
258	return;	258	return;
259	QPoint pos_reader = mPositionIndex;	259	QPoint pos_reader = mPositionIndex;
260		260
261	int offset_reader = reader->offset();	261	int offset_reader = reader->offset();
262	int offset_writer = mStamp->offset();	262	int offset_writer = mStamp->offset();
263	int d_offset = offset_writer - offset_reader; // offset on the stamp to the crown-cells	263	int d_offset = offset_writer - offset_reader; // offset on the stamp to the crown-cells
264		264
265	QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer);	265	QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer);
266	pos_reader-=QPoint(offset_reader, offset_reader);	266	pos_reader-=QPoint(offset_reader, offset_reader);
267	QPoint p;	267	QPoint p;
268		268
269	//float dom_height = (*m_dominanceGrid)[m_Position];	269	//float dom_height = (*m_dominanceGrid)[m_Position];
270	float local_dom;	270	float local_dom;
271		271
272	int x,y;	272	int x,y;
273	double sum=0.;	273	double sum=0.;
274	double value, own_value;	274	double value, own_value;
275	float *grid_value;	275	float *grid_value;
276	int reader_size = reader->size();	276	int reader_size = reader->size();
277	int rx = pos_reader.x();	277	int rx = pos_reader.x();
278	int ry = pos_reader.y();	278	int ry = pos_reader.y();
279	for (y=0;y<reader_size; ++y, ++ry) {	279	for (y=0;y<reader_size; ++y, ++ry) {
280	grid_value = mGrid->ptr(rx, ry);	280	grid_value = mGrid->ptr(rx, ry);
281	for (x=0;x<reader_size;++x) {	281	for (x=0;x<reader_size;++x) {
282		282
283	//p = pos_reader + QPoint(x,y);	283	//p = pos_reader + QPoint(x,y);
284	//if (m_grid->isIndexValid(p)) {	284	//if (m_grid->isIndexValid(p)) {
285	local_dom = mHeightGrid->valueAtIndex((rx+x)/5, ry/5).height; // ry: gets ++ed in outer loop, rx not	285	local_dom = mHeightGrid->valueAtIndex((rx+x)/5, ry/5).height; // ry: gets ++ed in outer loop, rx not
286	//local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) );	286	//local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) );
287		287
288	own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;	288	own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;
289	own_value = qMax(own_value, 0.02);	289	own_value = qMax(own_value, 0.02);
290	value = *grid_value++ / own_value; // remove impact of focal tree	290	value = *grid_value++ / own_value; // remove impact of focal tree
291	//if (value>0.)	291	//if (value>0.)
292	sum += value * (*reader)(x,y);	292	sum += value * (*reader)(x,y);
293		293
294	//} // isIndexValid	294	//} // isIndexValid
295	}	295	}
296	}	296	}
297	mLRI = sum;	297	mLRI = sum;
298	// read dominance field...	298	// read dominance field...
299	// this applies only if some trees are potentially higher than the dominant height grid	299	// this applies only if some trees are potentially higher than the dominant height grid
300	//if (dom_height < m_Height) {	300	//if (dom_height < m_Height) {
301	// if tree is higher than Z*, the dominance height	301	// if tree is higher than Z*, the dominance height
302	// a part of the crown is in "full light".	302	// a part of the crown is in "full light".
303	// m_statAboveZ++;	303	// m_statAboveZ++;
304	// mImpact = 1. - (1. - mImpact)*dom_height/m_Height;	304	// mImpact = 1. - (1. - mImpact)*dom_height/m_Height;
305	//}	305	//}
306	if (mLRI > 1.)	306	if (mLRI > 1.)
307	mLRI = 1.;	307	mLRI = 1.;
308	//qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;	308	//qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
309	mRU->addWLA(mLRI*mLeafArea, mLeafArea);	309	mRU->addWLA(mLRI*mLeafArea, mLeafArea);
310	}	310	}
311		311
312	void Tree::heightGrid_torus()	312	void Tree::heightGrid_torus()
313	{	313	{
314	// height of Z*	314	// height of Z*
315	const float cellsize = mHeightGrid->cellsize();	315	const float cellsize = mHeightGrid->cellsize();
316		316
317	QPoint p = QPoint(mPositionIndex.x()/5, mPositionIndex.y()/5); // pos of tree on height grid	317	QPoint p = QPoint(mPositionIndex.x()/5, mPositionIndex.y()/5); // pos of tree on height grid
318		318
319	// count trees that are on height-grid cells (used for stockable area)	319	// count trees that are on height-grid cells (used for stockable area)
320	mHeightGrid->valueAtIndex(p).count++;	320	mHeightGrid->valueAtIndex(p).count++;
321		321
322	int index_eastwest = mPositionIndex.x() % 5; // 4: very west, 0 east edge	322	int index_eastwest = mPositionIndex.x() % 5; // 4: very west, 0 east edge
323	int index_northsouth = mPositionIndex.y() % 5; // 4: northern edge, 0: southern edge	323	int index_northsouth = mPositionIndex.y() % 5; // 4: northern edge, 0: southern edge
324	int dist[9];	324	int dist[9];
325	dist[3] = index_northsouth * 2 + 1; // south	325	dist[3] = index_northsouth * 2 + 1; // south
326	dist[1] = index_eastwest * 2 + 1; // west	326	dist[1] = index_eastwest * 2 + 1; // west
327	dist[5] = 10 - dist[3]; // north	327	dist[5] = 10 - dist[3]; // north
328	dist[7] = 10 - dist[1]; // east	328	dist[7] = 10 - dist[1]; // east
329	dist[8] = qMax(dist[5], dist[7]); // north-east	329	dist[8] = qMax(dist[5], dist[7]); // north-east
330	dist[6] = qMax(dist[3], dist[7]); // south-east	330	dist[6] = qMax(dist[3], dist[7]); // south-east
331	dist[0] = qMax(dist[3], dist[1]); // south-west	331	dist[0] = qMax(dist[3], dist[1]); // south-west
332	dist[2] = qMax(dist[5], dist[1]); // north-west	332	dist[2] = qMax(dist[5], dist[1]); // north-west
333	dist[4] = 0; // center cell	333	dist[4] = 0; // center cell
334	/* 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:	334	/* 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:
335	index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.	335	index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
336	e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2	336	e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
337	*/	337	*/
338		338
339		339
340	int ringcount = int(floor(mHeight / cellsize)) + 1;	340	int ringcount = int(floor(mHeight / cellsize)) + 1;
341	int ix, iy;	341	int ix, iy;
342	int ring;	342	int ring;
343	QPoint pos;	343	QPoint pos;
344	float hdom;	344	float hdom;
345	int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer	345	int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer
346	for (ix=-ringcount;ix<=ringcount;ix++)	346	for (ix=-ringcount;ix<=ringcount;ix++)
347	for (iy=-ringcount; iy<=+ringcount; iy++) {	347	for (iy=-ringcount; iy<=+ringcount; iy++) {
348	ring = qMax(abs(ix), abs(iy));	348	ring = qMax(abs(ix), abs(iy));
349	QPoint pos(ix+p.x(), iy+p.y());	349	QPoint pos(ix+p.x(), iy+p.y());
350	if (mHeightGrid->isIndexValid(pos)) {	350	if (mHeightGrid->isIndexValid(pos)) {
351	float &rHGrid = mHeightGrid->valueAtIndex(torusIndex(pos.x(),10,bufferOffset), torusIndex(pos.y(),10,bufferOffset)).height;	351	float &rHGrid = mHeightGrid->valueAtIndex(torusIndex(pos.x(),10,bufferOffset), torusIndex(pos.y(),10,bufferOffset)).height;
352	if (rHGrid > mHeight) // skip calculation if grid is higher than tree	352	if (rHGrid > mHeight) // skip calculation if grid is higher than tree
353	continue;	353	continue;
354	int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)	354	int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
355	hdom = mHeight - dist[direction];	355	hdom = mHeight - dist[direction];
356	if (ring>1)	356	if (ring>1)
357	hdom -= (ring-1)*10;	357	hdom -= (ring-1)*10;
358		358
359	rHGrid = qMax(rHGrid, hdom); // write value	359	rHGrid = qMax(rHGrid, hdom); // write value
360	} // is valid	360	} // is valid
361	} // for (y)	361	} // for (y)
362	}	362	}
363		363
364	/// Torus version of read stamp (glued edges)	364	/// Torus version of read stamp (glued edges)
365	void Tree::readLIF_torus()	365	void Tree::readLIF_torus()
366	{	366	{
367	if (!mStamp)	367	if (!mStamp)
368	return;	368	return;
369	const Stamp *reader = mStamp->reader();	369	const Stamp *reader = mStamp->reader();
370	if (!reader)	370	if (!reader)
371	return;	371	return;
372	QPoint pos_reader = mPositionIndex;	372	QPoint pos_reader = mPositionIndex;
373		373
374	int offset_reader = reader->offset();	374	int offset_reader = reader->offset();
375	int offset_writer = mStamp->offset();	375	int offset_writer = mStamp->offset();
376	int d_offset = offset_writer - offset_reader; // offset on the stamp to the crown-cells	376	int d_offset = offset_writer - offset_reader; // offset on the stamp to the crown-cells
377		377
378	QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer);	378	QPoint pos_writer=pos_reader - QPoint(offset_writer, offset_writer);
379	pos_reader-=QPoint(offset_reader, offset_reader);	379	pos_reader-=QPoint(offset_reader, offset_reader);
380	QPoint p;	380	QPoint p;
381		381
382	//float dom_height = (*m_dominanceGrid)[m_Position];	382	//float dom_height = (*m_dominanceGrid)[m_Position];
383	float local_dom;	383	float local_dom;
384		384
385	int x,y;	385	int x,y;
386	double sum=0.;	386	double sum=0.;
387	double value, own_value;	387	double value, own_value;
388	float *grid_value;	388	float *grid_value;
389	int reader_size = reader->size();	389	int reader_size = reader->size();
390	int rx = pos_reader.x();	390	int rx = pos_reader.x();
391	int ry = pos_reader.y();	391	int ry = pos_reader.y();
392	int xt, yt; // wrapped coords	392	int xt, yt; // wrapped coords
393	int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer	393	int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer
394		394
395	for (y=0;y<reader_size; ++y, ++ry) {	395	for (y=0;y<reader_size; ++y, ++ry) {
396	grid_value = mGrid->ptr(rx, ry);	396	grid_value = mGrid->ptr(rx, ry);
397	for (x=0;x<reader_size;++x) {	397	for (x=0;x<reader_size;++x) {
398	xt = torusIndex(rx+x,50, bufferOffset);	398	xt = torusIndex(rx+x,50, bufferOffset);
399	yt = torusIndex(ry+y,50, bufferOffset);	399	yt = torusIndex(ry+y,50, bufferOffset);
400	grid_value = mGrid->ptr(xt,yt);	400	grid_value = mGrid->ptr(xt,yt);
401	//p = pos_reader + QPoint(x,y);	401	//p = pos_reader + QPoint(x,y);
402	//if (m_grid->isIndexValid(p)) {	402	//if (m_grid->isIndexValid(p)) {
403	local_dom = mHeightGrid->valueAtIndex(xt/5, yt/5).height; // ry: gets ++ed in outer loop, rx not	403	local_dom = mHeightGrid->valueAtIndex(xt/5, yt/5).height; // ry: gets ++ed in outer loop, rx not
404	//local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) );	404	//local_dom = m_dominanceGrid->valueAt( m_grid->cellCoordinates(p) );
405		405
406	own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;	406	own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;
407	own_value = qMax(own_value, 0.02);	407	own_value = qMax(own_value, 0.02);
408	value = *grid_value / own_value; // remove impact of focal tree	408	value = *grid_value / own_value; // remove impact of focal tree
409	//if (value>0.)	409	//if (value>0.)
410	sum += value * (*reader)(x,y);	410	sum += value * (*reader)(x,y);
411		411
412	//} // isIndexValid	412	//} // isIndexValid
413	}	413	}
414	}	414	}
415	mLRI = sum;	415	mLRI = sum;
416	// read dominance field...	416	// read dominance field...
417	// this applies only if some trees are potentially higher than the dominant height grid	417	// this applies only if some trees are potentially higher than the dominant height grid
418	//if (dom_height < m_Height) {	418	//if (dom_height < m_Height) {
419	// if tree is higher than Z*, the dominance height	419	// if tree is higher than Z*, the dominance height
420	// a part of the crown is in "full light".	420	// a part of the crown is in "full light".
421	// m_statAboveZ++;	421	// m_statAboveZ++;
422	// mImpact = 1. - (1. - mImpact)*dom_height/m_Height;	422	// mImpact = 1. - (1. - mImpact)*dom_height/m_Height;
423	//}	423	//}
424	if (mLRI > 1.)	424	if (mLRI > 1.)
425	mLRI = 1.;	425	mLRI = 1.;
426	//qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;	426	//qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
427	mRU->addWLA(mLRI*mLeafArea, mLeafArea);	427	mRU->addWLA(mLRI*mLeafArea, mLeafArea);
428	}	428	}
429		429
430		430
431	void Tree::resetStatistics()	431	void Tree::resetStatistics()
432	{	432	{
433	m_statPrint=0;	433	m_statPrint=0;
434	m_statCreated=0;	434	m_statCreated=0;
435	m_statAboveZ=0;	435	m_statAboveZ=0;
436	m_nextId=1;	436	m_nextId=1;
437	}	437	}
438		438
439	//////////////////////////////////////////////////	439	//////////////////////////////////////////////////
440	//// Growth Functions	440	//// Growth Functions
441	//////////////////////////////////////////////////	441	//////////////////////////////////////////////////
442		442
443		443
444	void Tree::grow()	444	void Tree::grow()
445	{	445	{
446	TreeGrowthData d;	446	TreeGrowthData d;
-		447	mAge++; // increase age
447	// step 1: get radiation from ressource unit: radiation (MJ/tree/year) total intercepted radiation for this tree per year!	448	// step 1: get radiation from ressource unit: radiation (MJ/tree/year) total intercepted radiation for this tree per year!
448	double radiation = mRU->interceptedRadiation(mLeafArea, mLRI);	449	double radiation = mRU->interceptedRadiation(mLeafArea, mLRI);
449	// step 2: get fraction of PARutilized, i.e. fraction of intercepted rad that is utiliziable (per year)	450	// step 2: get fraction of PARutilized, i.e. fraction of intercepted rad that is utiliziable (per year)
450		451
451	double raw_gpp_per_rad = mRU->ressourceUnitSpecies(species()).prod3PG().GPPperRad();	452	double raw_gpp_per_rad = mRU->ressourceUnitSpecies(species()).prod3PG().GPPperRad();
452	// GPP (without aging-effect) [gC] / year -> kg Biomass /GPP (0.001 2)	453	// GPP (without aging-effect) [gC] / year -> kg Biomass /GPP (0.001 2)
453	double raw_gpp = raw_gpp_per_rad * radiation * 0.001 * 2;	454	double raw_gpp = raw_gpp_per_rad * radiation * 0.001 * 2;
454		455
455	// apply aging	456	// apply aging
456	double gpp = raw_gpp * 0.6; // aging	-
-		457	const double aging_factor = mSpecies->aging(mHeight, mAge);
-		458	double gpp = raw_gpp * aging_factor;
457	d.NPP = gpp * 0.47; // respiration loss, transformation kg	459	d.NPP = gpp * 0.47; // respiration loss, transformation kg
458		460
459	DBGMODE(	461	DBGMODE(
460	if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) {	462	if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) {
461	DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP);	463	DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP);
462	dumpList(out); // add tree headers	464	dumpList(out); // add tree headers
463	out << radiation << raw_gpp << gpp << d.NPP;	-
-		465	out << radiation << raw_gpp << gpp << d.NPP << aging_factor;
464	}	466	}
465	); // DBGMODE()	467	); // DBGMODE()
466		468
467		469
468	partitioning(d); // split npp to compartments and grow (diameter, height)	470	partitioning(d); // split npp to compartments and grow (diameter, height)
469		471
470	mortality(d);	472	mortality(d);
471		473
472	mStressIndex = d.stress_index;	474	mStressIndex = d.stress_index;
473	}	475	}
474		476
475		477
476	inline void Tree::partitioning(TreeGrowthData &d)	478	inline void Tree::partitioning(TreeGrowthData &d)
477	{	479	{
478	if (isDebugging())	480	if (isDebugging())
479	enableDebugging(true);	481	enableDebugging(true);
480	double npp = d.NPP;	482	double npp = d.NPP;
481	double harshness = mRU->ressourceUnitSpecies(species()).prod3PG().harshness();	483	double harshness = mRU->ressourceUnitSpecies(species()).prod3PG().harshness();
482	// add content of reserve pool	484	// add content of reserve pool
483	npp += mNPPReserve;	485	npp += mNPPReserve;
484	const double foliage_mass_allo = species()->biomassFoliage(mDbh);	486	const double foliage_mass_allo = species()->biomassFoliage(mDbh);
485	const double reserve_size = 2 * foliage_mass_allo;	487	const double reserve_size = 2 * foliage_mass_allo;
486	double refill_reserve = qMin(reserve_size, 2.*mFoliageMass); // not always try to refill reserve 100%	488	double refill_reserve = qMin(reserve_size, 2.*mFoliageMass); // not always try to refill reserve 100%
487		489
488	double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta)	490	double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta)
489	// turnover rates	491	// turnover rates
490	const double &to_fol = species()->turnoverLeaf();	492	const double &to_fol = species()->turnoverLeaf();
491	const double &to_root = species()->turnoverRoot();	493	const double &to_root = species()->turnoverRoot();
492	// the turnover rate of wood depends on the size of the reserve pool:	494	// the turnover rate of wood depends on the size of the reserve pool:
493		495
494		496
495	double to_wood = refill_reserve / (mWoodyMass + refill_reserve);	497	double to_wood = refill_reserve / (mWoodyMass + refill_reserve);

Subversion Repositories public iLand

(root)/src/core/tree.cpp @ 1222 - Rev 167 → 169