WebSVN - public iLand - Diff - Rev 697 and 720


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." ).arg(id()));
    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 697	Rev 720
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." ).arg(id()));	92	throw IException(QString("Phenology::calculation(): was not able to determine the length of the vegetation period for group %1." ).arg(id()));
93	if (logLevelDebug())	93	if (logLevelDebug())
94	qDebug() << "Jolly-phenology. start" << mClimate->dayOfYear(day_start)->toString() << "stop" << mClimate->dayOfYear(day_stop)->toString();	94	qDebug() << "Jolly-phenology. start" << mClimate->dayOfYear(day_start)->toString() << "stop" << mClimate->dayOfYear(day_stop)->toString();
95	iday = 0;	95	iday = 0;
96	mDayStart = day_start;	96	mDayStart = day_start;
97	mDayEnd = day_stop;	97	mDayEnd = day_stop;
98	int bDay, bMon, eDay, eMon;	98	int bDay, bMon, eDay, eMon;
99	// convert yeardays to dates	99	// convert yeardays to dates
100	mClimate->toDate(day_start, &bDay, &bMon);	100	mClimate->toDate(day_start, &bDay, &bMon);
101	mClimate->toDate(day_stop, &eDay, &eMon);	101	mClimate->toDate(day_stop, &eDay, &eMon);
102	for (int i=0;i<12;i++) {	102	for (int i=0;i<12;i++) {
103	if (i<bMon \|\| i>eMon) {	103	if (i<bMon \|\| i>eMon) {
104	mPhenoFraction[i] = 0; // out of season	104	mPhenoFraction[i] = 0; // out of season
105	} else if (i>bMon && i<eMon) {	105	} else if (i>bMon && i<eMon) {
106	mPhenoFraction[i] = 1.; // full inside of season	106	mPhenoFraction[i] = 1.; // full inside of season
107	} else {	107	} else {
108	// fractions of month	108	// fractions of month
109	mPhenoFraction[i] = 1.;	109	mPhenoFraction[i] = 1.;
110	if (i==bMon)	110	if (i==bMon)
111	mPhenoFraction[i] -= (bDay+1) / double(mClimate->days(bMon));	111	mPhenoFraction[i] -= (bDay+1) / double(mClimate->days(bMon));
112	if (i==eMon)	112	if (i==eMon)
113	mPhenoFraction[i] -= (mClimate->days(eMon) - (bDay+1)) / double(mClimate->days(eMon));	113	mPhenoFraction[i] -= (mClimate->days(eMon) - (bDay+1)) / double(mClimate->days(eMon));
114	}	114	}
115	}	115	}
116		116
117	calculateChillDays();	117	calculateChillDays();
118	}	118	}
119		119
120		120
121	// ********* Chill-day calculations ******	121	// ********* Chill-day calculations ******
122	void Phenology::calculateChillDays(const int end_of_season)	122	void Phenology::calculateChillDays(const int end_of_season)
123	{	123	{
124	int iday = 0;	124	int iday = 0;
125	mChillDaysBefore = 0;	125	mChillDaysBefore = 0;
126	int days_after = 0;	126	int days_after = 0;
127	int last_day = end_of_season>0?end_of_season:mDayEnd;	127	int last_day = end_of_season>0?end_of_season:mDayEnd;
128	for (const ClimateDay *day = mClimate->begin(); day!=mClimate->end(); ++day, ++iday) {	128	for (const ClimateDay *day = mClimate->begin(); day!=mClimate->end(); ++day, ++iday) {
129	if (day->temperature>=-5 && day->temperature<5) {	129	if (day->temperature>=-5. && day->temperature<5.) {
130	if (iday<mDayStart)	130	if (iday<mDayStart)
131	mChillDaysBefore++;	131	mChillDaysBefore++;
132	if (iday>last_day)	132	if (iday>last_day)
133	days_after++;	133	days_after++;
134	}	134	}
135	}	135	}
136	if (GlobalSettings::instance()->currentYear()==1) {	136	if (GlobalSettings::instance()->currentYear()==1) {
137	// for the first simulation year, use the value of this autumn for the last years autumn	137	// for the first simulation year, use the value of this autumn for the last years autumn
138	mChillDaysAfterLastYear = days_after;	138	mChillDaysAfterLastYear = days_after;
139	} else {	139	} else {
140	mChillDaysAfterLastYear = mChillDaysAfter;	140	mChillDaysAfterLastYear = mChillDaysAfter;
141	}	141	}
142	mChillDaysAfter = days_after;	142	mChillDaysAfter = days_after;
143	}	143	}
144		144

Subversion Repositories public iLand

(root)/src/core/phenology.cpp @ 1222 - Rev 697 → 720