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/>.

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

/** @class SpeciesResponse

  @ingroup core

    Environmental responses relevant for production of a tree species on resource unit level.

    SpeciesResponse combines data from different sources and converts information about the environment

    into responses of a species. The spatial level is the "ResourceUnit", where homogenetiy of environmental factors

    is assumed. The temporal aggregation depends on the factor, but usually, the daily environmental data is

    aggregated to monthly response values (which subsequently are used during 3PG production).

 Sources are:

    - vapour pressure deficit (dryness of atmosphere): directly from climate data (daily)

    - soil water status (dryness of soil)(daily)

    - temperature: directly from climate data (daily)

    - phenology: @sa Phenology, combines several sources (quasi-monthly)

    - CO2: @sa SpeciesSet::co2Response() based on ambient CO2 level (climate data), nitrogen and soil water responses (yearly)

    - nitrogen: based on the amount of available nitrogen (yearly)

*/

#include "speciesresponse.h"

#include "resourceunit.h"

#include "species.h"

#include "resourceunitspecies.h"

#include "climate.h"

#include "model.h"

#include "watercycle.h"

#include "debugtimer.h"

SpeciesResponse::SpeciesResponse()

{

    mSpecies=0;

    mRu=0;

}

void SpeciesResponse::clear()

{

    for (int i=0;i<12;i++)

        mCO2Response[i]=mSoilWaterResponse[i]=mTempResponse[i]=mRadiation[i]=mUtilizableRadiation[i]=mVpdResponse[i]=0.;

    mNitrogenResponse = 0.;

    mTotalRadiation = 0.;

    mTotalUtilizeableRadiation = 0.;

}

void SpeciesResponse::setup(ResourceUnitSpecies *rus)

{

    mSpecies = rus->species();

    mRu = rus->ru();

    clear();

}

/// calculate responses for VPD and Soil Water. Return the minimum of those responses

/// @param psi_kPa psi of the soil in kPa

/// @param vpd vapor pressure deficit in kPa

/// @return minimum of soil water and vpd response

void SpeciesResponse::soilAtmosphereResponses(const double psi_kPa, const double vpd, double &rMinResponse) const

{

    double water_resp = mSpecies->soilwaterResponse(psi_kPa);

    double vpd_resp = mSpecies->vpdResponse( vpd );

    rMinResponse = qMin(water_resp, vpd_resp);

}

/// Main function that calculates monthly / annual species responses

void SpeciesResponse::calculate()

{

    DebugTimer tpg("SpeciesResponse::calculate");

    clear(); // reset values

    // calculate yearly responses

    const WaterCycle *water = mRu->waterCycle();

    const Phenology &pheno = mRu->climate()->phenology(mSpecies->phenologyClass());

    int veg_begin = pheno.vegetationPeriodStart();

    int veg_end = pheno.vegetationPeriodEnd();

    // yearly response

    const double nitrogen = mRu->resouceUnitVariables().nitrogenAvailable;

    // Nitrogen response: a yearly value based on available nitrogen

    mNitrogenResponse = mSpecies->nitrogenResponse( nitrogen );

    const double ambient_co2 = mRu->climate()->begin()->co2; // CO2 level of first day of year (co2 is static)

    double water_resp, vpd_resp, temp_resp, min_resp;

    double  utilizeable_radiation;

    int doy=0;

    int month;

    const ClimateDay *end = mRu->climate()->end();

    for (const ClimateDay *day=mRu->climate()->begin(); day!=end; ++day) {

        month = day->month - 1;

        // environmental responses

        water_resp = mSpecies->soilwaterResponse(water->psi_kPa(doy));

        vpd_resp = mSpecies->vpdResponse( day->vpd );

        temp_resp = mSpecies->temperatureResponse(day->temp_delayed);

        mSoilWaterResponse[month] += water_resp;

        mTempResponse[month] += temp_resp;

        mVpdResponse[month] += vpd_resp;

        mRadiation[month] += day->radiation;

        if (doy>=veg_begin && doy<=veg_end) {

            // environmental responses for the day

            // combine responses

            min_resp = qMin(qMin(vpd_resp, temp_resp), water_resp);

            // calculate utilizable radiation, Eq. 4, http://iland.boku.ac.at/primary+production

            utilizeable_radiation = day->radiation * min_resp;

        } else {

            utilizeable_radiation = 0.; // no utilizable radiation outside of vegetation period

            min_resp = 0.;

        }

        mUtilizableRadiation[month]+= utilizeable_radiation;

        doy++;

        //DBGMODE(

            if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dDailyResponses)) {

                DebugList &out = GlobalSettings::instance()->debugList(day->id(), GlobalSettings::dDailyResponses);

                // climatic variables

                out << mSpecies->id() << day->id() << mRu->index() << mRu->id(); // date << day->temperature << day->vpd << day->preciptitation << day->radiation;

                out << water_resp << temp_resp << vpd_resp << day->radiation << utilizeable_radiation;

            }

        //); // DBGMODE()

    }

    mTotalRadiation = mRu->climate()->totalRadiation();

    // monthly values

    for (int i=0;i<12;i++) {

        double days = mRu->climate()->days(i);

        mTotalUtilizeableRadiation += mUtilizableRadiation[i];

        mSoilWaterResponse[i]/=days;

        mTempResponse[i]/=days;

        mVpdResponse[i]/=days;

        mCO2Response[i] = mSpecies->speciesSet()->co2Response(ambient_co2,

                                                           mNitrogenResponse,

                                                           mSoilWaterResponse[i]);

    }

}