Subversion Repositories public iLand

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

**    iLand - an individual based forest landscape and disturbance model

**    http://iland.boku.ac.at

**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl

**

**    This program is free software: you can redistribute it and/or modify

**    it under the terms of the GNU General Public License as published by

**    the Free Software Foundation, either version 3 of the License, or

**    (at your option) any later version.

**

**    This program is distributed in the hope that it will be useful,

**    but WITHOUT ANY WARRANTY; without even the implied warranty of

**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

**    GNU General Public License for more details.

**

**    You should have received a copy of the GNU General Public License

**    along with this program.  If not, see <http://www.gnu.org/licenses/>.

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

#include "carbonflowout.h"

#include "globalsettings.h"

#include "model.h"

#include "resourceunit.h"

#include "resourceunitspecies.h"

#include "production3pg.h"

#include "soil.h"

CarbonFlowOut::CarbonFlowOut()

{

    setName("Carbon fluxes per RU or landscape/yr", "carbonflow");

    setDescription("Carbon fluxes per resource unit and year and/or aggregated for the full landscape. All values are reported on a per hectare basis (use the area provided in carbon or stand outputs to scale to realized values on the respective resource unit). " \

                   "For results limited to the project area, the data values need to be scaled to the stockable area.\n" \

                   "For landsacpe level outputs, data is always given per ha of (stockable) project area (i.e. scaling with stockable area is already included).\n" \

                   "Furthermore, the following sign convention is used in iLand: fluxes "\

                   "from the atmosphere to the ecosystem are positive, while C leaving the ecosystem is reported as negative C flux.\n" \

                   "You can specify a 'condition' to limit output execution to specific years (variable 'year'). " \

                   "The 'conditionRU' can be used to suppress resource-unit-level details; eg. specifying 'in(year,100,200,300)' limits output on reosurce unit level to the years 100,200,300 " \

                   "(leaving 'conditionRU' blank enables details per default).");

    columns() << OutputColumn::year() << OutputColumn::ru() << OutputColumn::id()

              << OutputColumn("area_ha", "total stockable area of the resource unit (or landscape) (ha)", OutDouble)

              << OutputColumn("GPP", "actually realized gross primary production, kg C; ((primary production|GPP)) including " \

                                         "the effect of decreasing productivity with age; note that a rough estimate of "\

                                         "((sapling growth and competition|#sapling C and N dynamics|sapling GPP)) is added to the GPP of adult trees here.", OutDouble)

              << OutputColumn("NPP", "net primary production, kg C; calculated as NPP=GPP-Ra; Ra, the autotrophic respiration (kg C/ha) is calculated as"\

                                     " a fixed fraction of GPP in iLand (see ((primary production|here)) for details). ", OutDouble)

              << OutputColumn("Rh", "heterotrophic respiration, kg C; sum of C released to the atmosphere from detrital pools, i.e."\

                                    " ((snag dynamics|#Snag decomposition|snags)), ((soil C and N cycling|downed deadwood, litter, and mineral soil)).", OutDouble)

              << OutputColumn("dist_loss", "disturbance losses, kg C; C that leaves the ecosystem as a result of disturbances, e.g. fire consumption", OutDouble)

              << OutputColumn("mgmt_loss", "management losses, kg C; C that leaves the ecosystem as a result of management interventions, e.g. harvesting", OutDouble)

              << OutputColumn("NEP", "net ecosytem productivity kg C, NEP=NPP - Rh - disturbance losses - management losses. "\

                                     "Note that NEP is also equal to the total net changes over all ecosystem C pools, as reported in the "\

                                     "carbon output (cf. [http://www.jstor.org/stable/3061028|Randerson et al. 2002])", OutDouble)

              << OutputColumn("cumNPP", "cumulative NPP, kg C. This is a running sum of NPP (including tree NPP and sapling carbon gain).", OutDouble)

              << OutputColumn("cumRh", "cumulative flux to atmosphere (heterotrophic respiration), kg C. This is a running sum of Rh.", OutDouble)

              << OutputColumn("cumNEP", "cumulative NEP (net ecosystem productivity), kg C. This is a running sum of NEP (positive values: carbon gain, negative values: carbon loss).", OutDouble);

}

void CarbonFlowOut::setup()

{

    // use a condition for to control execuation for the current year

    QString condition = settings().value(".condition", "");

    mCondition.setExpression(condition);

    condition = settings().value(".conditionRU", "");

    mConditionDetails.setExpression(condition);

}

void CarbonFlowOut::exec()

{

    Model *m = GlobalSettings::instance()->model();

    // global condition

    if (!mCondition.isEmpty() && mCondition.calculate(GlobalSettings::instance()->currentYear())==0.)

        return;

    bool ru_level = true;

    // switch off details if this is indicated in the conditionRU option

    if (!mConditionDetails.isEmpty() && mConditionDetails.calculate(GlobalSettings::instance()->currentYear())==0.)

        ru_level = false;

    double npp = 0.;

    int ru_count = 0;

    QVector<double> v(10, 0.); // 11 data values

    QVector<double>::iterator vit;

    foreach(ResourceUnit *ru, m->ruList()) {

        if (ru->id()==-1)

            continue; // do not include if out of project area

        if (!ru->soil() || !ru->snag()) {

            qDebug() << "CarbonFlowOut::exec: resource unit without soil or snags module - no output generated.";

            continue;

        }

        npp = 0.;

        double area_factor = ru->stockableArea() / cRUArea; //conversion factor

        npp += ru->statistics().npp() * biomassCFraction; // kg C/ha

        npp += ru->statistics().nppSaplings() * biomassCFraction; // kgC/ha

        // Snag pools are not scaled per ha (but refer to the stockable RU), soil pools and biomass statistics (NPP, ...) are scaled.

        double to_atm = ru->snag()->fluxToAtmosphere().C / area_factor; // from snags, kg/ha

        to_atm += ru->soil()->fluxToAtmosphere().C * cRUArea/10.; // soil: t/ha -> t/m2 -> kg/ha

        double to_dist = ru->snag()->fluxToDisturbance().C / area_factor; // convert to kgC/ha

        to_dist += ru->soil()->fluxToDisturbance().C * cRUArea/10.; // kgC/ha

        double to_harvest = ru->snag()->fluxToExtern().C / area_factor; // kgC/ha

        double nep = npp - to_atm - to_harvest - to_dist; // kgC/ha

        if (ru_level) {

            *this << currentYear() << ru->index() << ru->id() << area_factor; // keys

            *this <<  npp / cAutotrophicRespiration // GPP_act

                    << npp // NPP

                    << -to_atm // rh

                    << -to_dist // disturbance

                    << -to_harvest // management loss

                    << nep; // nep

            *this << ru->resouceUnitVariables().cumCarbonUptake << ru->resouceUnitVariables().cumCarbonToAtm << ru->resouceUnitVariables().cumNEP;

            writeRow();

        }

        // landscape level

        ++ru_count;

        vit = v.begin();

        *vit++ += area_factor; // total area in ha

        *vit++ += npp / cAutotrophicRespiration * area_factor; // GPP_act

        *vit++ += npp * area_factor; // NPP

        *vit++ += -to_atm * area_factor; // rh

        *vit++ += -to_dist * area_factor; // disturbance

        *vit++ += -to_harvest * area_factor; // management loss

        *vit++ += nep *area_factor; // net ecosystem productivity

        *vit++ += ru->resouceUnitVariables().cumCarbonUptake * area_factor; // cum. NPP

        *vit++ += ru->resouceUnitVariables().cumCarbonToAtm * area_factor; // cum. Rh

        *vit++ += ru->resouceUnitVariables().cumNEP * area_factor; // cum. NEP

    }

    // write landscape sums

    double total_stockable_area = v[0]; // total ha of stockable area

    if (ru_count==0. || total_stockable_area==0.)

        return;

    *this << currentYear() << -1 << -1; // codes -1/-1 for landscape level

    *this << v[0]; // stockable area [m2]

    for (int i=1;i<v.size();++i)

        *this << v[i] / total_stockable_area;

    writeRow();

}