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

#include "resourceunit.h"

#include "climate.h"

/** @class BBGenerations

    @ingroup beetlemodule

    BBGenerations calculates potential bark beetle generations based on climate data (including bark temperature).

  */

BBGenerations::BBGenerations()

{

}

/**

 * @brief BBGenerations::calculateGenerations

 * @param ru

 * @return the number of filial generation (i.e, main generations) + 0.5 if a sister brood develops also for the last generation

 */

double BBGenerations::calculateGenerations(const ResourceUnit *ru)

{

    calculateBarkTemperature(ru);

    // start at the 1. of April, and wait for 140.3 degree days (with a threhsold of 8.3 degrees)

    const ClimateDay *clim = ru->climate()->day(4-1,1-1); // 0-based indices

    const ClimateDay *last_day = ru->climate()->day(10-1,31-1); // 0-based indices -> Oct 31

    const ClimateDay * day_too_short = ru->climate()->dayOfYear(ru->climate()->sun().dayShorter14_5hrs()); // the first doy where the day is shorter than 14.5 hours

    double dd=0.;

    while (dd<140.3 && clim<last_day) {

        dd+=std::max(clim->max_temperature-8.3, 0.);

        ++clim;

    }

    // now wait for a decent warm day with tmax > 16.5 degrees

    while (clim->max_temperature<=16.5 && clim<last_day)

        ++clim;

    mGenerations.clear();

    // start with the first generation....

    BBGeneration bbgen;

    bbgen.start_day = ru->climate()->whichDayOfYear(clim);

    bbgen.is_sister_brood = false;

    bbgen.gen = 1; // the first generation

    mGenerations.append(bbgen);

    // process the generations

    int i=0;

    while (i<mGenerations.count()) {

        BBGeneration bb = mGenerations[i]; // deliberately make a copy and not a ref

        const ClimateDay *c = ru->climate()->dayOfYear(bb.start_day);

        int doy = bb.start_day;

        double base_temp = mEffectiveBarkTemp[doy];

        double t_sum=0.;

        bool added_sister_brood = false;

        while (c < last_day) {

            t_sum = (mEffectiveBarkTemp[doy]-base_temp) / 557.;

            if (t_sum>=1.) {

                if (c<day_too_short) {

                    // start a new parental generation (a full cycle), only if the current

                    // generation is also a filial generation.

                    if (bb.is_sister_brood)

                        mGenerations.append(BBGeneration(doy, true, bb.gen)); // sisterbrood: (true), keep counter of filial generation

                    else

                        mGenerations.append(BBGeneration(doy, false, bb.gen+1)); // new filial brood (false), start a new gen

                }

                break;

            } else if (t_sum>0.5 && !added_sister_brood) {

                // start a sister brood, *if* the maximum air temperature is high enough, and if the

                // length of the day > 14.5 hours

                if (c->max_temperature>16.5 && c<day_too_short) {

                    mGenerations.append(BBGeneration(doy, true, bb.gen)); // add a sister brood generation (true), keep gen. of originating brood

                    added_sister_brood = true;

                }

            }

            ++c; ++doy;

        }

        mGenerations[i].value = std::min(t_sum, 1.);

        ++i; // now process the next generation

    }

    mNSisterBroods = 0; mNFilialBroods = 0;

    for (i=0;i<mGenerations.count();++i) {

        if (mGenerations[i].value>0.6 && mGenerations[i].is_sister_brood==true)

            mNSisterBroods = std::max( mGenerations[i].gen, mNSisterBroods );

        if (mGenerations[i].value>0.6 && mGenerations[i].is_sister_brood==false)

            mNFilialBroods = mGenerations[i].gen;

    }

    // return the number of filial broods, and increase by 0.5 if a sister brood of the last filial generation has successfully developed.

    return mNFilialBroods + ( hasSisterBrood() ? 0.5 : 0.);

//    // now accumulate the generations

//    double filial_generations = 0.;

//    double sister_generations = 0.;

//    // accumulate possible number of offspring

//    const double offspring_factor = 25.; // assuming sex-ratio of 1:1 and 50 offspring per female (see p. 59)

//    int n=0;

//    int total_offspring = 0;

//    for (i=0;i<mGenerations.count();++i) {

//        if (mGenerations[i].value>0.6)  {

//            ++n;

//            total_offspring+= pow(offspring_factor, mGenerations[i].gen-1)*2.*offspring_factor;

//        }

//        if (mGenerations[i].value>0.6 && mGenerations[i].is_sister_brood==true)

//            sister_generations+=mGenerations[i].value;

//        if (mGenerations[i].value>0.6 && mGenerations[i].is_sister_brood==false)

//            filial_generations+=mGenerations[i].value;

//    }

//    if (logLevelDebug())

//        qDebug() << "rid" <<ru->id() << "filial/sister:" << filial_generations << sister_generations << "offspring:" << total_offspring << "started generations:" << n;

//    mNSisterBroods = sister_generations;

//    mNFilialBroods = filial_generations;

//    return filial_generations + (sister_generations>0?0.5:0);

}

/**

 * Calculate the bark temperatures for this year and a given resource unit.

 * Input: climate data (tmax (C), tmean (C), radiation (MJ/m2))

 * the LAI to estimate the radiation on the ground (Wh/m2)

 * Output: calculates for each day of the year the "effective" bark-temperature and saves a cumulative sum

 * Source: Schopf et al 2004: Risikoabschaetzung von Borkenkaefermassenkalamitaeten im Nationalpark Kalkalpen

 */

void BBGenerations::calculateBarkTemperature(const ResourceUnit *ru)

{

    // estimate the fraction of light on the ground (multiplier from 0..1)

    const double k = 0.5; // constant for the beer lambert function

    double ground_light_fraction = exp(-k * ru->leafAreaIndex() );

    mFrostDaysEarly=0;

    mFrostDaysLate=0;

    for (int i=0;i<ru->climate()->daysOfYear();++i) {

        const ClimateDay *clim = ru->climate()->dayOfYear(i);

        // radiation: MJ/m2/day -> the regression uses Wh/m2/day -> conversion-factor: 1/0.0036

        double rad_wh = clim->radiation*ground_light_fraction/0.0036;

        // calc. maximum bark temperature

        double bt_max=1.656 + 0.002955*rad_wh + 0.534*clim->max_temperature + 0.01884 * clim->max_temperature*clim->max_temperature;

        double diff_bt=0.;

        if (bt_max>=30.4)

            diff_bt=std::max(-310.667 + 9.603 * bt_max, 0.); // degree * hours

        // mean bark temperature

        double bt_mean=-0.173+0.0008518*rad_wh+1.054*clim->mean_temp();

        // effective:

        double bt_sum = std::max(bt_mean - 8.3, 0.)*24.; // degree * hours

        // corrected for days with very high bark temperature (>30.4 degrees)

        double bt_sum_eff = (bt_sum - diff_bt) / 24.; // degree * days

        mEffectiveBarkTemp[i] = (i>0? mEffectiveBarkTemp[i-1] : 0.) + bt_sum_eff;

        // frost days:

        if (clim->min_temperature < -15.) {

            if (i < ru->climate()->sun().longestDay())

                mFrostDaysEarly++;

            else

                mFrostDaysLate++;

        }

    }

}