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 "establishment.h"

#include "climate.h"

#include "species.h"

#include "resourceunit.h"

#include "resourceunitspecies.h"

#include "seeddispersal.h"

#include "model.h"

#include "grasscover.h"

#include "helper.h"

#include "debugtimer.h"

#include "watercycle.h"

/** @class Establishment

    Establishment deals with the establishment process of saplings.

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

    Prerequisites for establishment are:

    the availability of seeds: derived from the seed-maps per Species (@sa SeedDispersal)

    the quality of the abiotic environment (TACA-model): calculations are performend here, based on climate and species responses

    the quality of the biotic environment, mainly light: based on the LIF-values

  */

Establishment::Establishment()

{

    mPAbiotic = 0.;

}

Establishment::Establishment(const Climate *climate, const ResourceUnitSpecies *rus)

{

    setup(climate, rus);

}

void Establishment::setup(const Climate *climate, const ResourceUnitSpecies *rus)

{

    mClimate = climate;

    mRUS = rus;

    mPAbiotic = 0.;

    mPxDensity = 0.;

    mNumberEstablished = 0;

    if (climate==0)

        throw IException("Establishment::setup: no valid climate for a resource unit.");

    if (rus==0 || rus->species()==0 || rus->ru()==0)

        throw IException("Establishment::setup: important variable is null (are the species properly set up?).");

}

void Establishment::clear()

{

    mPAbiotic = 0.;

    mNumberEstablished = 0;

    mPxDensity = 0.;

    mTACA_min_temp=mTACA_chill=mTACA_gdd=mTACA_frostfree=false;

    mTACA_frostAfterBuds=0;

    mSumLIFvalue = 0.;

    mLIFcount = 0;

    mWaterLimitation = 0.;

}

double Establishment::calculateWaterLimitation(const int veg_period_start, const int veg_period_end)

{

    // return 1 if effect is disabled

    if (mRUS->species()->establishmentParameters().psi_min >= 0.)

        return 1.;

    double psi_min = mRUS->species()->establishmentParameters().psi_min;

    const WaterCycle *water = mRUS->ru()->waterCycle();

    int days = mRUS->ru()->climate()->daysOfYear();

    // two week (14 days) running average of actual psi-values on the resource unit

    static const int nwindow = 14;

    double psi_buffer[nwindow];

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

        psi_buffer[i] = 0.;

    double current_sum = 0.;

    int i_buffer = 0;

    double min_average = 9999999.;

    double current_avg = 0.;

    for (int day=0;day<days;++day) {

        // running average: remove oldest item, add new item in a ringbuffer

        current_sum -= psi_buffer[i_buffer];

        psi_buffer[i_buffer] = water->psi_kPa(day);

        current_sum += psi_buffer[i_buffer];

        if (day>=veg_period_start && day<=veg_period_end) {

            current_avg = day>0? current_sum / std::min(day, nwindow) : current_sum;

            min_average = std::min(min_average, current_avg);

        }

        // move to next value in the buffer

        i_buffer = ++i_buffer % nwindow;

    }

    if (min_average > 1000.)

        return 1.; // invalid vegetation period?

    // calculate the response of the species to this value of psi (see also Species::soilwaterResponse())

    const double psi_mpa = min_average / 1000.; // convert to MPa

    double result = limit( (psi_mpa - psi_min) / (-0.015 -  psi_min) , 0., 1.);

    return result;

}

/** Calculate the abiotic environemnt for seedling for a given species and a given resource unit.

 The model is closely based on the TACA approach of Nitschke and Innes (2008), Ecol. Model 210, 263-277

 more details: http://iland.boku.ac.at/establishment#abiotic_environment

 a model mockup in R: script_establishment.r

 */

void Establishment::calculateAbioticEnvironment()

{

    //DebugTimer t("est_abiotic"); t.setSilent();

    // make sure that required calculations (e.g. watercycle are already performed)

    const_cast<ResourceUnitSpecies*>(mRUS)->calculate(true); // calculate the 3pg module and run the water cycle (this is done only if that did not happen up to now); true: call comes from regeneration

    const EstablishmentParameters &p = mRUS->species()->establishmentParameters();

    const Phenology &pheno = mClimate->phenology(mRUS->species()->phenologyClass());

    mTACA_min_temp = true; // minimum temperature threshold

    mTACA_chill = false;  // (total) chilling requirement

    mTACA_gdd = false;   // gdd-thresholds

    mTACA_frostfree = false; // frost free days in vegetation period

    mTACA_frostAfterBuds = 0; // frost days after bud birst

    const ClimateDay *day = mClimate->begin();

    int doy = 0;

    double GDD=0.;

    double GDD_BudBirst = 0.;

    int chill_days = pheno.chillingDaysLastYear(); // chilling days of the last autumn

    int frost_free = 0;

    mTACA_frostAfterBuds = 0;

    bool chill_ok = false;

    bool buds_are_birst = false;

    int veg_period_end = pheno.vegetationPeriodEnd();

    if (veg_period_end >= 365)

        veg_period_end = mClimate->sun().dayShorter10_5hrs();

    for (; day!=mClimate->end(); ++day, ++doy) {

        // minimum temperature: if temp too low -> set prob. to zero

        if (day->min_temperature < p.min_temp)

            mTACA_min_temp = false;

        // count frost free days

        if (day->min_temperature > 0.)

            frost_free++;

        // chilling requirement, GDD, bud birst

        if (day->temperature>=-5. && day->temperature<5.)

            chill_days++;

        if (chill_days>p.chill_requirement)

            chill_ok=true;

        // GDDs above the base temperature are counted if beginning from the day where the chilling requirements are met

        // up to a fixed day ending the veg period

        if (doy<=veg_period_end) {

            // accumulate growing degree days

            if (chill_ok && day->temperature > p.GDD_baseTemperature) {

                GDD += day->temperature - p.GDD_baseTemperature;

                GDD_BudBirst += day->temperature - p.GDD_baseTemperature;

            }

            // if day-frost occurs, the GDD counter for bud birst is reset

            if (day->temperature <= 0.)

                GDD_BudBirst = 0.;

            if (GDD_BudBirst > p.bud_birst)

                buds_are_birst = true;

            if (doy<veg_period_end && buds_are_birst && day->min_temperature <= 0.)

                mTACA_frostAfterBuds++;

        }

    }

    // chilling requirement

    if (chill_ok)

        mTACA_chill = true;

    // GDD requirements

    if (GDD>p.GDD_min && GDD<p.GDD_max)

        mTACA_gdd = true;

    // frost free days in the vegetation period

    if (frost_free > p.frost_free)

        mTACA_frostfree = true;

    // if all requirements are met:

    if (mTACA_chill && mTACA_min_temp && mTACA_gdd && mTACA_frostfree) {

        // negative effect of frost events after bud birst

        double frost_effect = 1.;

        if (mTACA_frostAfterBuds>0)

            frost_effect = pow(p.frost_tolerance, sqrt(double(mTACA_frostAfterBuds)));

        // negative effect due to water limitation on establishment [1: no effect]

        mWaterLimitation = calculateWaterLimitation(pheno.vegetationPeriodStart(), pheno.vegetationPeriodLength());

        // combine drought and frost effect multiplicatively

        mPAbiotic = frost_effect * mWaterLimitation;

    } else {

        mPAbiotic = 0.; // if any of the requirements is not met

    }

}

void Establishment::writeDebugOutputs()

{

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

        DebugList &out = GlobalSettings::instance()->debugList(mRUS->ru()->index(), GlobalSettings::dEstablishment);

        // establishment details

        out << mRUS->species()->id() << mRUS->ru()->index() << mRUS->ru()->id();

        out << avgSeedDensity();

        out << TACAminTemp() << TACAchill() << TACAfrostFree() << TACgdd();

        out << TACAfrostDaysAfterBudBirst() << waterLimitation() << abioticEnvironment();

        out << mRUS->prod3PG().fEnvYear() << mRUS->constSaplingStat().newSaplings();

        //out << mSaplingStat.livingSaplings() << mSaplingStat.averageHeight() << mSaplingStat.averageAge() << mSaplingStat.averageDeltaHPot() << mSaplingStat.averageDeltaHRealized();

        //out << mSaplingStat.newSaplings() << mSaplingStat.diedSaplings() << mSaplingStat.recruitedSaplings() << mSpecies->saplingGrowthParameters().referenceRatio;

    }

    //); // DBGMODE()

    if ( logLevelDebug() )

        qDebug() << "establishment of RU" << mRUS->ru()->index() << "species" << mRUS->species()->id()

                 << "seeds density:" << avgSeedDensity()

                 << "abiotic environment:" << abioticEnvironment()

                 << "f_env,yr:" << mRUS->prod3PG().fEnvYear()

                 << "N(established):" << numberEstablished();

}