WebSVN - public iLand - Diff - Rev 639 and 671


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);
}

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

Subversion Repositories public iLand

(root)/src/core/phenology.cpp @ 1222 - Rev 639 → 671