# Program for simulating the effects of climate change on masting in valley oaks (Quercus lobata)
# based on the relationships between temperature, phenological synchrony, and the acorn crop

# described in Koenig et al.

# Written by Walt Koenig, August 2013, in R2.15.1

### This version does 100 runs of 33 years each, matching the field data from

### Hastings Reservation, central coastal California

# set up a dataframe for the results
tres <- rep(NA,100*8*2)
dim(tres) <- c(800,2)
colnames(tres) <- c('tem','cv')
tres <- as.data.frame(tres)

# do the entire set of trials 100 times
for (trial in 1:100){
# lets one know something is happening
if ((round(trial/10)*10)==trial){
print(trial)
flush.console()
}

# set up a dataframe for the results of each trial
runres <- rep(NA,6)
dim(runres) <- c(1,6)
runres <- as.data.frame(runres)
colnames(runres) <- c('yrs','xtemp','xcvtemp','xcvphen','xcrop','cv')

# Trials consist of a series of 800 runs; each series of 100 uses a value of the mean maximum temperature (tmean) of 16 to 23oC 
# Thus, the simulation involves 8 sets of 100 trials of 33 years each

# Values of tmean are normally distributed (SD = 2.0)

# The actual mean maximum temperature from 1 March - 30 April value (1980 - 2012) was 17.9 +/- 2.0oC (SD)

# based on data from Hastings Reservation

# This sets the value for the mean maximum temperature

for (tmean in 16:23){
# set up a dataframe for each year of the 33 years of the trial
res <- rep(NA,33*6)
dim(res) <- c(33,6)
res <- as.data.frame(res)
colnames(res) <- c('yr','xtemp','cvtemp','cvphen','acorn1','acorn')

# pick a random starting value for the acorn crop based on the actual data, which is normally 
# distributed with a mean of 1.866 +/- 1.052 (SD)

p.crop <- rnorm(1,1.866,1.052)
# make sure the acorn crop doesn't go below zero
p.crop <- ifelse(p.crop<0,0,p.crop)

# simulate values for the year's acorn crop

for (j in 1:33){

# pick a value for mean max April temp

xt <- rnorm(1,tmean,2.0)
# Use xt to calculate the cv in microclimate (cvmicro), based on the regression of cvmicro on tmean: 
# cvmicro <- (-.347 +/- .104SE) * tmean + (14.228 +/- 1.839SE)
v1 <- xt * rnorm(1,-.347,0.104)

v2 <- rnorm(1,14.228,1.839)
cvmicro <- v1 + v2

# Limit values to a reasonable range (between 3 and 15; see Fig. 3a)
# Note: altering these values changes the absolute values of the results,
# but does not alter the main conclusion (that is, the decline in the CV
# of acorn production with increasing tmean)

cvmicro <- ifelse(cvmicro<3,3,cvmicro)
cvmicro <- ifelse(cvmicro>15,15,cvmicro)

# Use cvmicro to calculate the cv in phenology (cvphen)
# actual lm is cvphen <- (3.810 +/- 1.202SE) * cvmicro + (-16.964 +/- 9.815SE)

v1 <- cvmicro * rnorm(1,3.810,1.202)
v2 <- rnorm(1,-16.964,9.815)
cvphen <- (v1 + v2)

# Again, limit cvphen to values to a reasonable range (between 3 and 25; see Fig. 3b)
# Again, this does not change the conclusions of the simulation

cvphen <- ifelse(cvphen<3,3,cvphen)
cvphen <- ifelse(cvphen>25,25,cvphen)

# use cvphen (and prior year's acorn crop) to calculate the acorn crop for the year

ac <- 4.446 - 0.1162*cvphen - 0.40485*p.crop

# make sure the crop doesn't go below zero
ac <- ifelse(ac<0,0,ac)

# save results for the year in dataframe res

res[j,] <- c(j,xt,cvmicro,cvphen,p.crop,ac)

# save the year's acorn crop for use the following year
p.crop <- ac
}     # end of the loop for the 33 years

# res now has the results from each year

# standardize the overall size of the acorn crop to 2 (an arbitrary, but not unreasonable, value)

m <- mean(res$acorn)
res$acorn.st <- res$acorn
res$acorn.st <- res$acorn* (2/m)

# calculate the cv of the acorn crop over the 33 years
m <- mean(res$acorn.st)
s <- sd(res$acorn.st)
cv <- s*100/m

# save the summary data for the run
runres <- rbind(runres,c(33,mean(res$xtemp),mean(res$cvtemp),mean(res$cvphen),m,cv))

}   # loop for the different values of tmean 

# save the data for the entire sets of trials

results <- runres[-1,]
ct <- c(16:23)
results <- cbind(results,ct)
results
nn <- ((trial-1)*8)+1
nn8 <- nn+7
tres[c(nn:nn8),1] <- results[,7]
tres[c(nn:nn8),2] <- results[,6]

}    # loop for the 100 sets of trials

# calculate means and SD
res2 <- results[,1:4]
colnames(res2) <- c('tem','cv','sd','se')

for (i in 16:23){
ii <- i-15
res2[ii,1] <- i
x <- subset(tres,tem==i)
res2[ii,2] <- mean(x$cv)
res2[ii,3] <- sd(x$cv)
res2[ii,4] <- sd(x$cv)/sqrt(100)
}

# summary of the results
res2

### Plot the results (FIG.  5 in the paper)

par(mfrow=c(2,2),las=1,cex=1,cex.lab=1.3,mar=c(0,4,6,1),las=0)
plot(results$xtemp,results$cv,xlab='',ylab='',pch=16,cex=0,cex.lab=.8,ylim=c(45,62),xlim=c(15.5,23.5),col.axis='white',bty='L')
par(las=1)
axis(2,cex.axis=0.8,hadj=.5,las=1)
axis(1,cex.axis=0.8,padj=-1)
par(las=0)
mtext('CV of acorn crop (%)',2,font=1,cex=1.1,line=2)
xx <- expression(paste("Mean max spring temp (",degree,"C)",sep=""))
mtext(xx,1,cex=1.1,las=1,line=1.5)
points(res2$tem,res2$cv,pch=21,col='Black',bg='Red',cex=1.3)
tic <- .1
for (i in 1:8){
yy <- res2[i,1]
pp <- res2[i,2]+2*res2[i,4]
pm <- res2[i,2]-2*res2[i,4]
lines(c(yy,yy),c(pp,pm))
lines(c(yy+tic,yy-tic),c(pp,pp))
lines(c(yy+tic,yy-tic),c(pm,pm))
}

# Plot the observed value
points(17.9,56.4,pch='+',cex=1.5)

##########################################################################################
##########################################################################################
##########################################################################################