Rev 441 | Rev 639 | 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())); |
||
| 611 | werner | 66 | if (logLevelDebug()) |
| 67 | qDebug() << "Jolly-phenology. start" << mClimate->dayOfYear(day_start)->toString() << "stop" << mClimate->dayOfYear(day_stop)->toString(); |
||
| 214 | werner | 68 | iday = 0; |
| 226 | werner | 69 | mDayStart = day_start; |
| 70 | mDayEnd = day_stop; |
||
| 214 | werner | 71 | int bDay, bMon, eDay, eMon; |
| 72 | // convert yeardays to dates |
||
| 73 | mClimate->toDate(day_start, &bDay, &bMon); |
||
| 74 | mClimate->toDate(day_stop, &eDay, &eMon); |
||
| 75 | for (int i=0;i<12;i++) { |
||
| 76 | if (i<bMon || i>eMon) { |
||
| 77 | mPhenoFraction[i] = 0; // out of season |
||
| 78 | } else if (i>bMon && i<eMon) { |
||
| 79 | mPhenoFraction[i] = 1.; // full inside of season |
||
| 80 | } else { |
||
| 81 | // fractions of month |
||
| 82 | mPhenoFraction[i] = 1.; |
||
| 83 | if (i==bMon) |
||
| 84 | mPhenoFraction[i] -= (bDay+1) / double(mClimate->days(bMon)); |
||
| 85 | if (i==eMon) |
||
| 86 | mPhenoFraction[i] -= (mClimate->days(eMon) - (bDay+1)) / double(mClimate->days(eMon)); |
||
| 87 | } |
||
| 88 | } |
||
| 434 | werner | 89 | |
| 90 | calculateChillDays(); |
||
| 214 | werner | 91 | } |
| 434 | werner | 92 | |
| 93 | |||
| 94 | // *********** Chill-day calculations ******** |
||
| 436 | werner | 95 | void Phenology::calculateChillDays(const int end_of_season) |
| 434 | werner | 96 | { |
| 97 | int iday = 0; |
||
| 98 | mChillDaysBefore = 0; |
||
| 99 | int days_after = 0; |
||
| 436 | werner | 100 | int last_day = end_of_season>0?end_of_season:mDayEnd; |
| 434 | werner | 101 | for (const ClimateDay *day = mClimate->begin(); day!=mClimate->end(); ++day, ++iday) { |
| 102 | if (day->temperature>=-5 && day->temperature<5) { |
||
| 103 | if (iday<mDayStart) |
||
| 104 | mChillDaysBefore++; |
||
| 436 | werner | 105 | if (iday>last_day) |
| 434 | werner | 106 | days_after++; |
| 107 | } |
||
| 108 | } |
||
| 109 | if (GlobalSettings::instance()->currentYear()==1) { |
||
| 110 | // for the first simulation year, use the value of this autumn for the last years autumn |
||
| 111 | mChillDaysAfterLastYear = days_after; |
||
| 112 | } else { |
||
| 113 | mChillDaysAfterLastYear = mChillDaysAfter; |
||
| 114 | } |
||
| 115 | mChillDaysAfter = days_after; |
||
| 116 | } |