#Read lake-level reconstruction
LL = read.csv("Deep_LLrecon_012414.csv",header=TRUE)

##Input variables
Area.W = 84677
#Watershed area in square meters

Area.L = 8925
#Lake area in square meters

##Calculate change in lake volume
Vol = LL$LakeElev.Anom/100 * Area.L

Vol.High = LL$UpperBound/100 * Area.L

Vol.Low = LL$LowerBound/100 * Area.L

##Calculate change in volume across the watershed in mm

Total = matrix(Vol*1000/Area.W)
Total.High = matrix(Vol.High*1000/Area.W)
Total.Low = matrix(Vol.Low*1000/Area.W)

##Incorporate equilibrium time (in years)

Eq.time.min = 0.1
Eq.time.max = 5
Eq.time.step = 0.1

Eq.times=seq(from=Eq.time.min, to=Eq.time.max,by=Eq.time.step)

mb = function(Eq.times,x) x/Eq.times

MB = apply(Total,1,mb,Eq.times=Eq.times)
MB.High = apply(Total.High,1,mb,Eq.times=Eq.times)
MB.Low = apply(Total.Low,1,mb,Eq.times=Eq.times)

##Calc mean P-E

MB.all=rbind(MB,MB.High,MB.Low)

P.E = apply(-MB.all,2,median)
Offset=P.E[1]
P.E=P.E-Offset

#1SD
P.E.682 = apply(-MB.all,2,quantile,prob=0.682,na.rm=T)
P.E.682=P.E.682-Offset
P.E.318 = apply(-MB.all,2,quantile,prob=0.318,na.rm=T)
P.E.318=P.E.318-Offset

#25% & 75%
P.E.75 = apply(-MB.all,2,quantile,prob=0.75,na.rm=T)
P.E.75=P.E.75-Offset
P.E.25 = apply(-MB.all,2,quantile,prob=0.25,na.rm=T)
P.E.25=P.E.25-Offset

Recon=cbind(LL$Age,P.E,P.E.75,P.E.682,P.E.318,P.E.25)
colnames(Recon)=c("Age","P-E.mm","75%","Plus1SD","Minus1SD","25%")
Recon=data.frame(Recon)

#Write output file
write.csv(Recon, "Deep_PErecon_012414.csv", row.names = FALSE)

###Plot

quartz(width=5,height=5)

par(mai=c(1,0.75,0.5,0.25))

Length=max(LL$Age)

plot(LL$Age,P.E,type="l",xlim=c(Length,0),xaxp=c(Length,0,Length/3000),ylim=c(min(P.E,na.rm=T),0),xlab="Cal yr before AD 1950",ylab="P-E (mm/yr)",frame.plot=F,col="blue",lwd=0.25)

P.E.682[1]=min(P.E.318,na.rm=T)
P.E.682[length(P.E.682[P.E.682<Inf])]=min(P.E.318,na.rm=T)
P.E.318[1]=min(P.E.318,na.rm=T)
P.E.318[length(P.E.318[P.E.318<Inf])]=min(P.E.318,na.rm=T)

lines(LL$Age,P.E.75,col="light blue")
polygon(LL$Age,P.E.682,col="light blue",border=F)
polygon(LL$Age,P.E.318,col="white",border=F)
lines(LL$Age,P.E.25,col="light blue")
lines(LL$Age,P.E,col="dark blue",lwd=2)