Rev 440 | Rev 611 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
211 | werner | 2 | #include "global.h" |
3 | #include "phenology.h" |
||
4 | |||
5 | #include "climate.h" |
||
214 | werner | 6 | #include "floatingaverage.h" |
7 | |||
8 | /// "ramp" function: value < min: 0, value>max: 1: in between linear interpolation |
||
9 | double ramp(const double &value, const double minValue, const double maxValue) |
||
10 | { |
||
11 | Q_ASSERT(minValue!=maxValue); |
||
12 | if (value<minValue) return 0.; |
||
13 | if (value>maxValue) return 1.; |
||
14 | return (value-minValue) / (maxValue - minValue); |
||
15 | } |
||
16 | |||
17 | /** calculates the phenology according to Jolly et al. 2005. |
||
226 | werner | 18 | The calculation is performed for a given "group" and a present "climate". |
214 | werner | 19 | */ |
20 | void Phenology::calculate() |
||
21 | { |
||
436 | werner | 22 | if (id()==0) { |
23 | // for needles: just calculate the chilling requirement for the establishment |
||
440 | werner | 24 | // i.e.: use the "bottom line" of 10.5 hrs daylength for the end of the vegetation season |
25 | calculateChillDays(mClimate->sun().dayShorter10_5hrs()); |
||
214 | werner | 26 | return; |
436 | werner | 27 | } |
214 | werner | 28 | const ClimateDay *end = mClimate->end(); |
29 | double vpd,temp,daylength; |
||
30 | double gsi; // combined factor of effect of vpd, temperature and day length |
||
31 | int iday = 0; |
||
32 | bool inside_period = !mClimate->sun().northernHemishere(); // on northern hemisphere 1.1. is in winter |
||
33 | int day_start=-1, day_stop=-1; |
||
34 | int day_wait_for=-1; |
||
35 | |||
36 | FloatingAverage floater(21); // a three-week floating average |
||
37 | for (const ClimateDay *day = mClimate->begin(); day!=end; ++day, ++iday) { |
||
38 | |||
39 | if (day_wait_for>=0 && iday<day_wait_for) |
||
40 | continue; |
||
41 | vpd = 1. - ramp(day->vpd, mMinVpd, mMaxVpd); // high value for low vpd |
||
441 | werner | 42 | temp = ramp(day->min_temperature, mMinTemp, mMaxTemp); |
214 | werner | 43 | daylength = ramp( mClimate->sun().daylength(iday), mMinDayLength, mMaxDayLength); |
44 | gsi = vpd * temp * daylength; |
||
45 | gsi = floater.add(gsi); |
||
46 | if (!inside_period && gsi>0.5) { |
||
47 | // switch from winter -> summer |
||
48 | inside_period = true; |
||
49 | day_start = iday; |
||
50 | if (day_stop!=-1) |
||
51 | break; |
||
52 | day_wait_for = mClimate->sun().longestDay(); |
||
53 | } else if (inside_period && gsi<0.5) { |
||
54 | // switch from summer to winter |
||
55 | day_stop = iday; |
||
56 | if (day_start!=-1) |
||
57 | break; // finished |
||
58 | day_wait_for = mClimate->sun().longestDay(); |
||
59 | inside_period = false; |
||
60 | } |
||
61 | } |
||
62 | day_start -= 10; // three-week-floating average: subtract 10 days |
||
63 | day_stop -= 10; |
||
64 | 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())); |
||
66 | qDebug() << "Jolly-phenology. start" << mClimate->dayOfYear(day_start)->toString() << "stop" << mClimate->dayOfYear(day_stop)->toString(); |
||
67 | iday = 0; |
||
226 | werner | 68 | mDayStart = day_start; |
69 | mDayEnd = day_stop; |
||
214 | werner | 70 | int bDay, bMon, eDay, eMon; |
71 | // convert yeardays to dates |
||
72 | mClimate->toDate(day_start, &bDay, &bMon); |
||
73 | mClimate->toDate(day_stop, &eDay, &eMon); |
||
74 | for (int i=0;i<12;i++) { |
||
75 | if (i<bMon || i>eMon) { |
||
76 | mPhenoFraction[i] = 0; // out of season |
||
77 | } else if (i>bMon && i<eMon) { |
||
78 | mPhenoFraction[i] = 1.; // full inside of season |
||
79 | } else { |
||
80 | // fractions of month |
||
81 | mPhenoFraction[i] = 1.; |
||
82 | if (i==bMon) |
||
83 | mPhenoFraction[i] -= (bDay+1) / double(mClimate->days(bMon)); |
||
84 | if (i==eMon) |
||
85 | mPhenoFraction[i] -= (mClimate->days(eMon) - (bDay+1)) / double(mClimate->days(eMon)); |
||
86 | } |
||
87 | } |
||
434 | werner | 88 | |
89 | calculateChillDays(); |
||
214 | werner | 90 | } |
434 | werner | 91 | |
92 | |||
93 | // *********** Chill-day calculations ******** |
||
436 | werner | 94 | void Phenology::calculateChillDays(const int end_of_season) |
434 | werner | 95 | { |
96 | int iday = 0; |
||
97 | mChillDaysBefore = 0; |
||
98 | int days_after = 0; |
||
436 | werner | 99 | int last_day = end_of_season>0?end_of_season:mDayEnd; |
434 | werner | 100 | for (const ClimateDay *day = mClimate->begin(); day!=mClimate->end(); ++day, ++iday) { |
101 | if (day->temperature>=-5 && day->temperature<5) { |
||
102 | if (iday<mDayStart) |
||
103 | mChillDaysBefore++; |
||
436 | werner | 104 | if (iday>last_day) |
434 | werner | 105 | days_after++; |
106 | } |
||
107 | } |
||
108 | if (GlobalSettings::instance()->currentYear()==1) { |
||
109 | // for the first simulation year, use the value of this autumn for the last years autumn |
||
110 | mChillDaysAfterLastYear = days_after; |
||
111 | } else { |
||
112 | mChillDaysAfterLastYear = mChillDaysAfter; |
||
113 | } |
||
114 | mChillDaysAfter = days_after; |
||
115 | } |