Rev 697 | Rev 779 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
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 |