WebSVN - public iLand - Diff - Rev 1221 and 1222


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/phenology.cpp':
/********************************************************************************************
**    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 "global.h"
#include "phenology.h"

#include "climate.h"
#include "floatingaverage.h"

/// "ramp" function: value < min: 0, value>max: 1: in between linear interpolation
double ramp(const double &value, const double minValue, const double maxValue)
{
    Q_ASSERT(minValue!=maxValue);
    if (value<minValue) return 0.;
    if (value>maxValue) return 1.;
    return (value-minValue) / (maxValue - minValue);
}

/** @class Phenology phenology submodule.
  @ingroup core
  The Phenology submodule calculates the length of the growing season according to the model of Jolly et al (2005). The calculation
  is performed for species-groups (i.e.: species are lumped together to groups) and a given climate (i.e. worst case: for each ResourceUnit).

  See http://iland.boku.ac.at/phenology for details.
  */

/** calculates the phenology according to Jolly et al. 2005.
  The calculation is performed for a given "group" and a present "climate".
*/
void Phenology::calculate()
{
    if (id()==0) {
        // for needles: just calculate the chilling requirement for the establishment
        // i.e.: use the "bottom line" of 10.5 hrs daylength for the end of the vegetation season
        calculateChillDays(mClimate->sun().dayShorter10_5hrs());
        return;
    }
    const ClimateDay  *end = mClimate->end();
    double vpd,temp,daylength;
    double gsi; // combined factor of effect of vpd, temperature and day length
    int iday = 0;
    bool inside_period = !mClimate->sun().northernHemishere(); // on northern hemisphere 1.1. is in winter
    int day_start=-1, day_stop=-1;
    int day_wait_for=-1;

    FloatingAverage floater(21); // a three-week floating average
    for (const ClimateDay *day = mClimate->begin(); day!=end; ++day, ++iday) {

        if (day_wait_for>=0 && iday<day_wait_for)
            continue;
        vpd = 1. - ramp(day->vpd, mMinVpd, mMaxVpd); // high value for low vpd
        temp = ramp(day->min_temperature, mMinTemp, mMaxTemp);
        daylength = ramp( mClimate->sun().daylength(iday), mMinDayLength, mMaxDayLength);
        gsi = vpd * temp * daylength;
        gsi = floater.add(gsi);
        if (!inside_period && gsi>0.5) {
            // switch from winter -> summer
            inside_period = true;
            day_start = iday;
            if (day_stop!=-1)
                break;
            day_wait_for = mClimate->sun().longestDay();
        } else if (inside_period && gsi<0.5) {
            // switch from summer to winter
            day_stop = iday;
            if (day_start!=-1)
                break; // finished
            day_wait_for = mClimate->sun().longestDay();
            inside_period = false;
        }
    }
    day_start -= 10; // three-week-floating average: subtract 10 days
    day_stop -= 10;
    if (day_start < -1 || day_stop<-1) {
        //throw IException(QString("Phenology::calculation(): was not able to determine the length of the vegetation period for group %1. climate table: '%2'." ).arg(id()).arg(mClimate->name()));
        qDebug() << "Phenology::calculation(): vegetation period is 0 for group" << id() << ", climate table: " << mClimate->name();
        day_start = mClimate->daysOfYear()-1; // last day of the year, never reached
        day_stop = day_start; // never reached
    }
    if (logLevelDebug())
        qDebug() << "Jolly-phenology. start" << mClimate->dayOfYear(day_start)->toString() << "stop" << mClimate->dayOfYear(day_stop)->toString();
    iday = 0;
    mDayStart = day_start;
    mDayEnd = day_stop;
    int bDay, bMon, eDay, eMon;
    // convert yeardays to dates
    mClimate->toDate(day_start, &bDay, &bMon);
    mClimate->toDate(day_stop, &eDay, &eMon);
    for (int i=0;i<12;i++) {
        if (i<bMon || i>eMon) {
            mPhenoFraction[i] = 0; // out of season
        } else if (i>bMon && i<eMon) {
            mPhenoFraction[i] = 1.; // full inside of season
        } else {
            // fractions of month
            mPhenoFraction[i] = 1.;
            if (i==bMon)
                mPhenoFraction[i] -=  (bDay+1) / double(mClimate->days(bMon));
            if (i==eMon)
                mPhenoFraction[i] -=  (mClimate->days(eMon) - (bDay+1)) / double(mClimate->days(eMon));
        }
    }

    calculateChillDays();
}


// *********** Chill-day calculations ********
void Phenology::calculateChillDays(const int end_of_season)
{
    int iday = 0;
    mChillDaysBefore = 0;
    int days_after = 0;
    int last_day = end_of_season>0?end_of_season:mDayEnd;
    for (const ClimateDay *day = mClimate->begin(); day!=mClimate->end(); ++day, ++iday) {
        if (day->temperature>=-5. && day->temperature<5.) {
            if (iday<mDayStart)
                mChillDaysBefore++;
            if (iday>last_day)
                days_after++;
        }
    }
    if (GlobalSettings::instance()->currentYear()==1) {
        // for the first simulation year, use the value of this autumn for the last years autumn
        mChillDaysAfterLastYear = days_after;
    } else {
        mChillDaysAfterLastYear  = mChillDaysAfter;
    }
    mChillDaysAfter = days_after;
}
 

Rev 1221	Rev 1222
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/phenology.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/phenology.cpp':
2	/********************************************************************************************	2	/********************************************************************************************
3	** iLand - an individual based forest landscape and disturbance model	3	** iLand - an individual based forest landscape and disturbance model
4	** http://iland.boku.ac.at	4	** http://iland.boku.ac.at
5	** Copyright (C) 2009- Werner Rammer, Rupert Seidl	5	** Copyright (C) 2009- Werner Rammer, Rupert Seidl
6	**	6	**
7	** This program is free software: you can redistribute it and/or modify	7	** This program is free software: you can redistribute it and/or modify
8	** it under the terms of the GNU General Public License as published by	8	** it under the terms of the GNU General Public License as published by
9	** the Free Software Foundation, either version 3 of the License, or	9	** the Free Software Foundation, either version 3 of the License, or
10	** (at your option) any later version.	10	** (at your option) any later version.
11	**	11	**
12	** This program is distributed in the hope that it will be useful,	12	** This program is distributed in the hope that it will be useful,
13	** but WITHOUT ANY WARRANTY; without even the implied warranty of	13	** but WITHOUT ANY WARRANTY; without even the implied warranty of
14	** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the	14	** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15	** GNU General Public License for more details.	15	** GNU General Public License for more details.
16	**	16	**
17	** You should have received a copy of the GNU General Public License	17	** You should have received a copy of the GNU General Public License
18	** along with this program. If not, see <http://www.gnu.org/licenses/>.	18	** along with this program. If not, see <http://www.gnu.org/licenses/>.
19	********************************************************************************************/	19	********************************************************************************************/
20		20
21	#include "global.h"	21	#include "global.h"
22	#include "phenology.h"	22	#include "phenology.h"
23		23
24	#include "climate.h"	24	#include "climate.h"
25	#include "floatingaverage.h"	25	#include "floatingaverage.h"
26		26
27	/// "ramp" function: value < min: 0, value>max: 1: in between linear interpolation	27	/// "ramp" function: value < min: 0, value>max: 1: in between linear interpolation
28	double ramp(const double &value, const double minValue, const double maxValue)	28	double ramp(const double &value, const double minValue, const double maxValue)
29	{	29	{
30	Q_ASSERT(minValue!=maxValue);	30	Q_ASSERT(minValue!=maxValue);
31	if (value<minValue) return 0.;	31	if (value<minValue) return 0.;
32	if (value>maxValue) return 1.;	32	if (value>maxValue) return 1.;
33	return (value-minValue) / (maxValue - minValue);	33	return (value-minValue) / (maxValue - minValue);
34	}	34	}
35		35
36	/** @class Phenology phenology submodule.	36	/** @class Phenology phenology submodule.
37	@ingroup core	37	@ingroup core
38	The Phenology submodule calculates the length of the growing season according to the model of Jolly et al (2005). The calculation	38	The Phenology submodule calculates the length of the growing season according to the model of Jolly et al (2005). The calculation
39	is performed for species-groups (i.e.: species are lumped together to groups) and a given climate (i.e. worst case: for each ResourceUnit).	39	is performed for species-groups (i.e.: species are lumped together to groups) and a given climate (i.e. worst case: for each ResourceUnit).
40		40
41	See http://iland.boku.ac.at/phenology for details.	41	See http://iland.boku.ac.at/phenology for details.
42	*/	42	*/
43		43
44	/** calculates the phenology according to Jolly et al. 2005.	44	/** calculates the phenology according to Jolly et al. 2005.
45	The calculation is performed for a given "group" and a present "climate".	45	The calculation is performed for a given "group" and a present "climate".
46	*/	46	*/
47	void Phenology::calculate()	47	void Phenology::calculate()
48	{	48	{
49	if (id()==0) {	49	if (id()==0) {
50	// for needles: just calculate the chilling requirement for the establishment	50	// for needles: just calculate the chilling requirement for the establishment
51	// i.e.: use the "bottom line" of 10.5 hrs daylength for the end of the vegetation season	51	// i.e.: use the "bottom line" of 10.5 hrs daylength for the end of the vegetation season
52	calculateChillDays(mClimate->sun().dayShorter10_5hrs());	52	calculateChillDays(mClimate->sun().dayShorter10_5hrs());
53	return;	53	return;
54	}	54	}
55	const ClimateDay *end = mClimate->end();	55	const ClimateDay *end = mClimate->end();
56	double vpd,temp,daylength;	56	double vpd,temp,daylength;
57	double gsi; // combined factor of effect of vpd, temperature and day length	57	double gsi; // combined factor of effect of vpd, temperature and day length
58	int iday = 0;	58	int iday = 0;
59	bool inside_period = !mClimate->sun().northernHemishere(); // on northern hemisphere 1.1. is in winter	59	bool inside_period = !mClimate->sun().northernHemishere(); // on northern hemisphere 1.1. is in winter
60	int day_start=-1, day_stop=-1;	60	int day_start=-1, day_stop=-1;
61	int day_wait_for=-1;	61	int day_wait_for=-1;
62		62
63	FloatingAverage floater(21); // a three-week floating average	63	FloatingAverage floater(21); // a three-week floating average
64	for (const ClimateDay *day = mClimate->begin(); day!=end; ++day, ++iday) {	64	for (const ClimateDay *day = mClimate->begin(); day!=end; ++day, ++iday) {
65		65
66	if (day_wait_for>=0 && iday<day_wait_for)	66	if (day_wait_for>=0 && iday<day_wait_for)
67	continue;	67	continue;
68	vpd = 1. - ramp(day->vpd, mMinVpd, mMaxVpd); // high value for low vpd	68	vpd = 1. - ramp(day->vpd, mMinVpd, mMaxVpd); // high value for low vpd
69	temp = ramp(day->min_temperature, mMinTemp, mMaxTemp);	69	temp = ramp(day->min_temperature, mMinTemp, mMaxTemp);
70	daylength = ramp( mClimate->sun().daylength(iday), mMinDayLength, mMaxDayLength);	70	daylength = ramp( mClimate->sun().daylength(iday), mMinDayLength, mMaxDayLength);
71	gsi = vpd * temp * daylength;	71	gsi = vpd * temp * daylength;
72	gsi = floater.add(gsi);	72	gsi = floater.add(gsi);
73	if (!inside_period && gsi>0.5) {	73	if (!inside_period && gsi>0.5) {
74	// switch from winter -> summer	74	// switch from winter -> summer
75	inside_period = true;	75	inside_period = true;
76	day_start = iday;	76	day_start = iday;
77	if (day_stop!=-1)	77	if (day_stop!=-1)
78	break;	78	break;
79	day_wait_for = mClimate->sun().longestDay();	79	day_wait_for = mClimate->sun().longestDay();
80	} else if (inside_period && gsi<0.5) {	80	} else if (inside_period && gsi<0.5) {
81	// switch from summer to winter	81	// switch from summer to winter
82	day_stop = iday;	82	day_stop = iday;
83	if (day_start!=-1)	83	if (day_start!=-1)
84	break; // finished	84	break; // finished
85	day_wait_for = mClimate->sun().longestDay();	85	day_wait_for = mClimate->sun().longestDay();
86	inside_period = false;	86	inside_period = false;
87	}	87	}
88	}	88	}
89	day_start -= 10; // three-week-floating average: subtract 10 days	89	day_start -= 10; // three-week-floating average: subtract 10 days
90	day_stop -= 10;	90	day_stop -= 10;
91	if (day_start < -1 \|\| day_stop<-1) {	91	if (day_start < -1 \|\| day_stop<-1) {
92	//throw IException(QString("Phenology::calculation(): was not able to determine the length of the vegetation period for group %1. climate table: '%2'." ).arg(id()).arg(mClimate->name()));	92	//throw IException(QString("Phenology::calculation(): was not able to determine the length of the vegetation period for group %1. climate table: '%2'." ).arg(id()).arg(mClimate->name()));
93	qDebug() << "Phenology::calculation(): vegetation period is 0 for group" << id() << ", climate table: " << mClimate->name();	93	qDebug() << "Phenology::calculation(): vegetation period is 0 for group" << id() << ", climate table: " << mClimate->name();
94	day_start = mClimate->daysOfYear()-1; // last day of the year, never reached	94	day_start = mClimate->daysOfYear()-1; // last day of the year, never reached
95	day_stop = day_start; // never reached	95	day_stop = day_start; // never reached
96	}	96	}
97	if (logLevelDebug())	97	if (logLevelDebug())
98	qDebug() << "Jolly-phenology. start" << mClimate->dayOfYear(day_start)->toString() << "stop" << mClimate->dayOfYear(day_stop)->toString();	98	qDebug() << "Jolly-phenology. start" << mClimate->dayOfYear(day_start)->toString() << "stop" << mClimate->dayOfYear(day_stop)->toString();
99	iday = 0;	99	iday = 0;
100	mDayStart = day_start;	100	mDayStart = day_start;
101	mDayEnd = day_stop;	101	mDayEnd = day_stop;
102	int bDay, bMon, eDay, eMon;	102	int bDay, bMon, eDay, eMon;
103	// convert yeardays to dates	103	// convert yeardays to dates
104	mClimate->toDate(day_start, &bDay, &bMon);	104	mClimate->toDate(day_start, &bDay, &bMon);
105	mClimate->toDate(day_stop, &eDay, &eMon);	105	mClimate->toDate(day_stop, &eDay, &eMon);
106	for (int i=0;i<12;i++) {	106	for (int i=0;i<12;i++) {
107	if (i<bMon \|\| i>eMon) {	107	if (i<bMon \|\| i>eMon) {
108	mPhenoFraction[i] = 0; // out of season	108	mPhenoFraction[i] = 0; // out of season
109	} else if (i>bMon && i<eMon) {	109	} else if (i>bMon && i<eMon) {
110	mPhenoFraction[i] = 1.; // full inside of season	110	mPhenoFraction[i] = 1.; // full inside of season
111	} else {	111	} else {
112	// fractions of month	112	// fractions of month
113	mPhenoFraction[i] = 1.;	113	mPhenoFraction[i] = 1.;
114	if (i==bMon)	114	if (i==bMon)
115	mPhenoFraction[i] -= (bDay+1) / double(mClimate->days(bMon));	115	mPhenoFraction[i] -= (bDay+1) / double(mClimate->days(bMon));
116	if (i==eMon)	116	if (i==eMon)
117	mPhenoFraction[i] -= (mClimate->days(eMon) - (bDay+1)) / double(mClimate->days(eMon));	117	mPhenoFraction[i] -= (mClimate->days(eMon) - (bDay+1)) / double(mClimate->days(eMon));
118	}	118	}
119	}	119	}
120		120
121	calculateChillDays();	121	calculateChillDays();
122	}	122	}
123		123
124		124
125	// ********* Chill-day calculations ******	125	// ********* Chill-day calculations ******
126	void Phenology::calculateChillDays(const int end_of_season)	126	void Phenology::calculateChillDays(const int end_of_season)
127	{	127	{
128	int iday = 0;	128	int iday = 0;
129	mChillDaysBefore = 0;	129	mChillDaysBefore = 0;
130	int days_after = 0;	130	int days_after = 0;
131	int last_day = end_of_season>0?end_of_season:mDayEnd;	131	int last_day = end_of_season>0?end_of_season:mDayEnd;
132	for (const ClimateDay *day = mClimate->begin(); day!=mClimate->end(); ++day, ++iday) {	132	for (const ClimateDay *day = mClimate->begin(); day!=mClimate->end(); ++day, ++iday) {
133	if (day->temperature>=-5. && day->temperature<5.) {	133	if (day->temperature>=-5. && day->temperature<5.) {
134	if (iday<mDayStart)	134	if (iday<mDayStart)
135	mChillDaysBefore++;	135	mChillDaysBefore++;
136	if (iday>last_day)	136	if (iday>last_day)
137	days_after++;	137	days_after++;
138	}	138	}
139	}	139	}
140	if (GlobalSettings::instance()->currentYear()==1) {	140	if (GlobalSettings::instance()->currentYear()==1) {
141	// for the first simulation year, use the value of this autumn for the last years autumn	141	// for the first simulation year, use the value of this autumn for the last years autumn
142	mChillDaysAfterLastYear = days_after;	142	mChillDaysAfterLastYear = days_after;
143	} else {	143	} else {
144	mChillDaysAfterLastYear = mChillDaysAfter;	144	mChillDaysAfterLastYear = mChillDaysAfter;
145	}	145	}
146	mChillDaysAfter = days_after;	146	mChillDaysAfter = days_after;
147	}	147	}
148		148

Subversion Repositories public iLand

(root)/src/core/phenology.cpp - Rev 1221 → 1222