#### Calculate monthly growing degree-days using the sine curve method of Baskerville & Emin (1969)
#### This code written by Daniel Gavin, University of Oregon.
#### Use daily data from the historical climatology network (HCN).
#### http://cdiac.ornl.gov/epubs/ndp/ushcn/daily.html
#### The program will interpolate to fill in any missing values using mean values in the bracketing
#### 30-year period for that date.
#### Reference:
#### Baskerville, G. L., Emin, P., 1969. Rapid estimation of heat accumulation from maximum and
#### minimum temperatures. Ecology 50 (3), 514-517.
### Code updated Nov 18th 2009 to handle occurrence of months with no or too little data.
### As is below: months with less than 15 days of data will not be used. For months with
### more than 15 days of data, any missing days will be based on 30-year
### means for the +/- 15 years around the year with the missing value.
### COPY AND PASTE THE FUNCTION BELOW INTO R. SPECIFY MONTHLY LAPSE RATES AS SHOWN IN EXAMPLE BELOW.
### THEN, run as: gdd_results <- gdd(station.elev=90, site.elev=90, lapse=lr, gdd.min=5, gdd.max=35, method=1)
### REQUIRED INPUTS:
### 1. A dialog box will appear, where you can choose an HCN-formatted (and comma-delimited)
### input file (see example).
### 2. The program adjusts for elevation between the study site of interest and the elevation
### of the climate station.
### In the call of the program, you need to specify the values for:
### station.elev and site.elev (elevations of the climate data station and study site in meters)
### lapse = array of 12 monthly temperature lapse rates in deg C/km.
### gdd.min = base temperature for GDD summation
### gdd.max = maximum temperature for GDD summation (subtract area > gdd.max)
### method = 1 (use base, do not use gdd.max)
2 (use base and gdd.max)
3 (use base and maxT, but subtract twice the area exceeding gdd.max)
lr <- c(-2.33,-3.49,-5.06,-6.63,-5.65,-6.45,-5.82,-5.4,-5.75,-5.24,-5.87,-3.82) ## degC/1000m
gdd <- function(station.elev=90, site.elev=90, lapse=lr, gdd.min=5, gdd.max=35, method=1){
### vars to use: Day Month Year TMIN TMAX PCP
hcn<-read.table(file=file.choose(), header=TRUE, sep=",", na.strings=".")
## FILL DATA ARRAYS WITH VALUES
## [year,month,day]
year.1 <- hcn$Year[1] # first year in the data file
n.years <- hcn$Year[length(hcn$Year)]-year.1+1
tmin <- array(data=NA, dim=c(n.years, 12, 31))
tmax <- array(data=NA, dim=c(n.years, 12, 31))
pcp <- array(data=NA, dim=c(n.years, 12, 31))
### CONVERT TO DEG C AND TO MM PRECIP
for(i in 1:length(hcn$TMAX)){
ye <- hcn$Year[i]-year.1+1
mo <- hcn$Month[i]
da <- hcn$Day[i]
tmin[ye,mo,da] <- (hcn$TMIN[i]-32)*5/9
tmax[ye,mo,da] <- (hcn$TMAX[i]-32)*5/9
pcp[ye,mo,da] <- (hcn$PCP[i])*25.4
}
## len.month[year,month] = number of days in month
len.month <- matrix(NA,n.years,12)
for(ye in 1:n.years){
for(mo in 1:12){
len.month[ye,mo]<-length(subset(hcn$TMIN,as.integer(hcn$Year-year.1+1)==as.integer(ye) & as.integer(hcn$Month)==as.integer(mo)))
}
}
### GDD[year,month]
gdd <- matrix(NA,n.years,13)
colnames(gdd) <- c("Year","J","F","M","A","M","J","J","A","S","O","N","D")
###lapse-adjust temperatures
lapse.adjust <- lapse*(site.elev-station.elev)/1000
for(ye in 1:n.years){
gdd[ye,1]<-ye+year.1-1
for(mo in 1:12){
gdd.sum<-0
if(len.month[ye,mo] > 0){ ## if month exists
tmintmax<-tmin[ye,mo,]+tmax[ye,mo,]
if(length(tmintmax[!is.na(tmintmax)]) > 15){ ###only use months with more than 15 days of observations.
for(da in 1:len.month[ye,mo]){ ### assumes that if month exists, there are all days listed in input file.
tmin.site <- tmin[ye,mo,da]+lapse.adjust[mo]
tmax.site <- tmax[ye,mo,da]+lapse.adjust[mo]
# fill in missing values
# if value is missing, give it the bracketing mean over the 30yr normal for that day and lapse adjust this instead.
if(is.na(tmin.site)){
tmin.site <- mean(tmin[max(1,ye-15):min(ye+15,n.years),mo,da],na.rm=TRUE)+lapse.adjust[mo]
}
if(is.na(tmax.site)){
tmax.site <- mean(tmax[max(1,ye-15):min(ye+15,n.years),mo,da],na.rm=TRUE)+lapse.adjust[mo]
}
tmean.site <- (tmin.site+tmax.site)/2
## B.day calculation (A.day and B.day are parameters used in GDD calculation)
if( tmax.site>gdd.min && tmin.site<gdd.min && (tmean.site!=gdd.min) ){ ##
num1<-(gdd.min-tmean.site)/((tmax.site-tmin.site)/2)
B.day <- atan(num1/sqrt((-num1*num1)+1))
} else {
B.day <- 0
}
## A.day calculation
if(tmax.site<gdd.min && tmin.site<gdd.min){
A.day<-0
} else {
if(tmin.site>=gdd.min){
A.day <- 0
} else {
A.day <- B.day
}
}
## GDD calculation
if(tmax.site<=gdd.min){
gdd.day<-0
} else {
if(tmean.site>gdd.min && tmin.site>=gdd.min){
gdd.day <- tmean.site-gdd.min
} else {
gdd.day <- ( (((tmax.site-tmin.site)/2)*cos(A.day)) - (gdd.min-tmean.site) * ((3.14/2)-A.day) )/3.14
}
}
### If max cutoff is to be used temperature exceeds cutoff point.
if(method>1){
if(tmax.site>gdd.max){
if((tmax.site-gdd.max)>(gdd.max-gdd.min)){
tmax.site<-gdd.max+(gdd.max-gdd.min)
}
## B.day calculation (A.day and B.day are parameters used in GDD calculation)
if( tmax.site>gdd.max && tmin.site<gdd.max && (tmean.site!=gdd.max) ){ ##
num1<-(gdd.max-tmean.site)/((tmax.site-tmin.site)/2)
B.day <- atan(num1/sqrt((-num1*num1)+1))
} else {
B.day <- 0
}
## A.day calculation
if(tmax.site<gdd.max && tmin.site<gdd.max){
A.day <- 0
} else {
if(tmin.site>=gdd.max){
A.day <- 0
} else {
A.day <- B.day
}
}
## GDD.cutoff calculation
if(tmean.site>gdd.max && tmin.site >= gdd.max){
gdd.cutoff <- gdd.day
} else {
gdd.cutoff <- ( (((tmax.site-tmin.site)/2)*cos(A.day)) - (gdd.max-tmean.site) * ((3.14/2)-A.day) )/3.14
}
if(method==2){
gdd.day <- max(gdd.day-(1*gdd.cutoff),0) ## can't have negative DD for this version
} else {
gdd.day <- max(gdd.day-(2*gdd.cutoff),0) ## can't have negative DD for this version
}
}
}
if(!is.na(gdd.day)){
gdd.sum<-gdd.sum+gdd.day
}
} # next day
gdd[ye,(mo+1)]<-gdd.sum
} #end if (test if month has any values)
} #end if (test if month has more than 15 days of observations)
} # next month
} # next year
return(gdd)
}