#### 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){ 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){ 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) }