BigMallardBOOT <- function(metdata=BigMallard_prec_pet,lakedata=BigMallard_vol,metdata21=ESTmonprpet1921to2017)  {
  
  #   To make plots, BigMallardRun(), you can leave () blank and run the defaults
  #   metdata contains monthly precip and pet for water years 1991-2016
  #   lakedata contains annual lake level and volume for the end of water years 1991-2016
  #   no missing values allowed!
  
  # Get working directory
  WD <- getwd()
  
  # ************* append climate division data to the end of the Big Mallard Marsh metdata ************
  
  metdataxx <- as.matrix(metdata21[,c(1,2,5,6)])
  row.names(metdataxx) <- NULL
  namesxx <- names(metdata)
  metdata <- as.matrix(metdata)
  row.names(metdata) <- NULL
  metdataxx <- data.frame(rbind(metdata,metdataxx))
  names(metdataxx) <- namesxx
  
  mprec <- metdataxx$prec
  mpet <- metdataxx$pet
  
  # *************************************
  
  alev <- lakedata$elev
  avol <- lakedata$vol/1000000
  
  #  Determine elevation/volume relation
  
  eft <- c(1820.7, 1822.7, 1825.7, 1827.7, 1829.7, 1832.7, 1835.7)
  vmcf <- c(0,5, 97, 214, 350, 582, 838)
  evcoefs <- summary(nls(vmcf~a*(eft-1820.7)^b,start=list(a=2,b=2)))$coef
  xxtmp <- seq(1820.7, 1835.7,0.1); yytmp <- evcoefs[1]*(xxtmp-1820.7)^evcoefs[2]
  
  # replace volumes with computed values from equation
  
  wyr <- 1991:2016
  nyrs <- 26
  for (i in 1:nyrs) { avol[i] <- bigmalvol(alev[i])}
  
  # compute quarterly precip and pet from monthly data and plot it
  
  qprec <- apply(matrix(mprec,ncol=3,byrow=T),1,sum)
  qpet <- apply(matrix(mpet,ncol=3,byrow=T),1,sum)
  
  dyrqtr <- 1991+ seq(0.5,103.5,1)/4
  
  # Ready to fit the water balance model. The parameters are estimated by
  #  minimizing the rmse between the estimated and observed annual volume
  #  changes, subject to parameter constraints. Parameter defmax is the maximum
  #  deficit (in inches), which is also the holding capacity of the basin
  #  feeding the lake. Parameter barea is the basin area, expressed as a multiple
  #  of the reference lake area. The reference lake area for Big Mallard Marsh is 
  #  the lake area at the spill elevation (see code for bigmalwatbal)
  # used weights to block out number 6 for 1997, we have a very large volume change between 1996 and 1997 
 
  obsvolchng <- avol[2:nyrs]  
  obsvoldif <-  (avol[-1] - avol[-nyrs])  
  
  tmpout <- nls(obsvolchng~bigmalwatbal(qprec,qpet,avol,barea,3.5),algorithm="port",lower=c(5),
                upper=c(6),start=list(barea=5.5),weights=c(rep(1,5),1,rep(1,19)))
  
  barea <- summary(tmpout)$coef[1]
  maxdef <- 3.5
  estvolchng <- bigmalwatbal(qprec,qpet,avol,barea,3.5)
  estvoldif <-  estvolchng - c(avol[1],estvolchng[-(nyrs-1)]) 
  nyrs <- 25 ; pkfirst <- seq(1,nyrs*4,4)
  
  
  rho <- as.character(round(cor.test(estvoldif,obsvoldif)$estimate,3))
  bias <- as.character(round(mean(obsvoldif-estvoldif),2))
  rmse <- as.character(round(mean((obsvoldif-estvoldif)^2)^0.5,2))
  bch1 <- as.character(round(barea,2)); bch2 <- as.character(round(maxdef,2))

  ##  Plots estimated vs. observed lake volumes
  #Create JPEG
  jpeg(paste(getwd(),"/BigMallard_Est_Obs.jpg", sep = ""), height = 5, width = 5, units = "in", res = 300)
  #svg(paste(getwd(),"/BigMallard_Est_Obs.svg", sep = ""), height = 5, width = 5)
  #postscript(paste(getwd(),"/BigMallard_Est_Obs.eps", sep = ""), height = 5, width = 5)
  
  plot(watbalxx$dqtrxx,watbalxx$vest*22.9568753,type="l",cex=0.6,xlim=c(1990,2020),ylim=c(0,1200*22.9568753), col = "red",
       ylab="Lake volume (Acre-feet)", xlab = "Year", tck = 0.05, las = 1, cex.axis = 0.8, xaxs = "i", yaxs = "i")
  points(wyr+1,avol*22.9568753,cex=1,lwd=1.5) 
  mtext(side=3,line=-1,"Big Mallard Marsh")
  mtext(side = 3, line = -2, cex = 0.8, paste("rho =", rho))
  
  dev.off()
  
  ### long-term simulations 
  estvolxx <- bigmalwatbal1(qprec,qpet,avol,barea,3.5,nyrs=123)
  nyrs <- 122
  plot(watbalxx$dqtrxx,watbalxx$vest,type="l", col = "red", ylim = c(0,1200), 
       ylab = "Lake Volume (MCF)", xaxs = "i", xlim = c(1940,2017))
  estvolchng <- bigmalwatbal(qprec,qpet,avol,barea,3.5,prmult=1)
  points(watbalxx$dqtrxx,watbalxx$vest,type="l", col="black")
  mtext(side=3,line=1,"Big Mallard Long-term simulation using climate division data") 
  
  #Create PDF
  pdf(paste(getwd(),"/BigMallard_Long.pdf", sep = ""), width = 9, height = 7)
  #svg(paste(getwd(),"/BigMallard_Long.svg", sep = ""), width = 9, height = 7)
  #postscript(paste(getwd(),"/BigMallard_Long.eps", sep = ""), width = 9, height = 7)
  
  
  ### Start of bootstrapping plot
  
  nboot <- 1000  # number of bootstrap simulations
  
  estvolxx <- bigmalwatbal1(qprec,qpet,avol,barea,3.5,nyrs=123)
  nyrs <- 122
  pkr <- watbalxx$dqtrxx > 1940
  options(scipen = 10)
  plot(watbalxx$dqtrxx[pkr],watbalxx$vest[pkr]*22.9568753,type="l", col = "purple", ylim = c(0,1000*22.9568753), 
       ylab = "Lake volume, acre-feet", xlab = "Year", tck = 0.05, las = 1, cex.axis = 0.8,
       xlim = c(1940,2100), xaxs = "i", yaxs = "i", xaxt = "n")
  axis(1, at = c(1940, 1992, 2016, 2040, 2067), tck = 0.05,  cex.axis = 0.8)
  axis(1, at =c(2020, 2025, 2030, 2035), tck = 0.02, cex.axis = 0.6, padj = -3)
  #mtext(side=3,line=1,"Big Mallard Marsh Long-term Simulation") 
  mycol <- c("black", "blue", "cyan", "chartreuse2")
  
 ##### Use bigmalbootvols saved dataframe. Uncomment lines below if first time running script.
  
 ##  Bootstrapping
 ##   keep the same initial climate inputs (first 26 years) but bootstrap random segments onto the end
 ##   first, replace 30's drought (30-39) with repeat of 1921-29 to reduce weight
  
 #  qprec[(36*4+1):(46*4)] <- qprec[(26*4+1):(36*4)]
 #  qpet[(36*4+1):(46*4)] <- qpet[(26*4+1):(36*4)]
 # 
 #  nboot <- 1000
 #  simvols <- matrix(nrow=75*4,ncol=nboot)
 #  for (isim in 1:nboot) {
 #  qprecboot <- rep(NA,76*4)
 #  qpetboot <- rep(NA,76*4)
 #  qprecboot[1:(26*4)] <- qprec[1:(26*4)]
 #  qpetboot[1:(26*4)] <- qpet[1:(26*4)]
 #  randur <- floor(runif(1)*49)+1
 #  ranyr <- 19+floor(runif(1)*78)
 #  eyrs <- 97-ranyr
 #      if(eyrs<randur)  {
 #  qprecboot[(26*4+1):((26+eyrs)*4)] <- qprec[(ranyr*4+1):((ranyr+eyrs)*4)]
 #  qprecboot[((26+eyrs)*4+1):((26+randur)*4)] <- qprec[1:((randur-eyrs)*4)]
 #  qpetboot[(26*4+1):((26+eyrs)*4)] <- qpet[(ranyr*4+1):((ranyr+eyrs)*4)]
 #  qpetboot[((26+eyrs)*4+1):((26+randur)*4)] <- qpet[1:((randur-eyrs)*4)]
 #     } else {
 #  qprecboot[(26*4+1):((26+randur)*4)] <- qprec[(ranyr*4+1):((ranyr+randur)*4)]
 #  qpetboot[(26*4+1):((26+randur)*4)] <- qpet[(ranyr*4+1):((ranyr+randur)*4)]
 #                   }
 #  ranyr <- 19+floor(runif(1)*78)
 #  randur2 <- 50-randur
 #  eyrs <- 97-ranyr
 #      if(eyrs<randur2)  {
 #  qprecboot[((26+randur)*4+1):((26+randur+eyrs)*4)] <- qprec[(ranyr*4+1):((ranyr+eyrs)*4)]
 #  qprecboot[((26+randur+eyrs)*4+1):(76*4)] <- qprec[1:((randur2-eyrs)*4)]
 #  qpetboot[((26+randur)*4+1):((26+randur+eyrs)*4)] <- qpet[(ranyr*4+1):((ranyr+eyrs)*4)]
 #  qpetboot[((26+randur+eyrs)*4+1):(76*4)] <- qpet[1:((randur2-eyrs)*4)]
 #     } else {
 #  qprecboot[((26+randur)*4+1):(76*4)] <- qprec[(ranyr*4+1):((ranyr+randur2)*4)]
 #  qpetboot[((26+randur)*4+1):(76*4)] <- qpet[(ranyr*4+1):((ranyr+randur2)*4)]
 #                   }
 #  
 #  estvolxx <- bigmalwatbal1(qprecboot,qpetboot,avol,barea,3.5,nyrs=76)
 #  simvols[,isim] <- watbalxx$vest
 #                       }
 # 
 #  simvols <- round(simvols,1)
 #  quants <- matrix(nrow=300,ncol=3)
 #  for (j in 1:200) {quants[100+j,] <- quantile(simvols[100+j,],p=c(.5,.25,.1)) }
 #  quants <- round(quants,1)
 #  simvols <- cbind(97+watbalxx$dqtrxx,simvols,quants)
 #  
 #  # Use bigmalbootvols saved dataframe. Uncomment lines below if first time running script.
 #  bigmalbootvols <<- data.frame(simvols)
 #  bigmalbootvolsyear <- bigmalbootvols[,1]
 #  bigmalbootvolsnoyear <- bigmalbootvols[,-1]
    # Convert volumes from cubic feet to acre feet
 #  bigmalbootvolsnoyear <- bigmalbootvolsnoyear*22.9568753
 #  bigmalbootvols <- cbind(bigmalbootvolsyear, bigmalbootvolsnoyear)
 #  bigmalbootvols <<- data.frame(bigmalbootvols)
  
  ############################################################################
  #### plot the first 4 simulations ("spaghetti plot")
  # remove 1-96 to remove the 1992-2017 part of the forecast

  # bigmalbootvolsyear <- bigmalbootvols[,1]
  # bigmalbootvolsnoyear <- bigmalbootvols[,-1]
  # # Convert volumes from cubic feet to acre feet
  # bigmalbootvolsnoyear <- bigmalbootvolsnoyear*22.9568753
  # bigmalbootvols <- cbind(bigmalbootvolsyear, bigmalbootvolsnoyear)
  
  colors <- c("orange", "blue", "cyan", "chartreuse2")
  for (ii in 1:4) { lines(bigmalbootvols[96:1004,1], bigmalbootvols[96:1004,1+ii], col= colors[ii], lwd=0.5) }

  #places exceedance text at the end of the line
  xtext <- bigmalbootvols[300,1]
  ytext50 <- bigmalbootvols[300, 1+nboot+1]
  ytext75 <- bigmalbootvols[300, 1+nboot+2]
  ytext90 <- bigmalbootvols[300, 1+nboot+3]
  text(xtext+15,ytext50, label = "50 percent exceedance level", cex = 0.6)
  text(xtext+15,ytext75, label = "75 percent exceedance level", cex = 0.6)
  text(xtext+15,ytext90, label = "90 percent exceedance level", cex = 0.6)
  abline(v = 2016, cex = 0.5)
  abline(v = 1992, cex = 0.5)
  
  #### overlay the 50 percent, 75 percent, and 90 percent exceedance levels
  lines(supsmu(bigmalbootvols[,1],bigmalbootvols[,1+nboot+1], span = 0.1),col="red",lwd=2.5)
  lines(supsmu(bigmalbootvols[,1],bigmalbootvols[,1+nboot+2], span = 0.1),col="red",lwd=2.5)
  lines(supsmu(bigmalbootvols[,1],bigmalbootvols[,1+nboot+3], span = 0.1),col="red",lwd=2.5)
  
  # Lake is 50% full of WY2017
  vol50per <- bigmalbootvols[97, 1001]*0.5
  abline(h = vol50per, cex = 0.5, type = "l", lty = 2)
  text(1942, vol50per-(vol50per*.05), label = "50 percent lake volume", cex = 0.8, adj = 0)
  
  # Lake is 75% full of WY2017
  vol75per <- bigmalbootvols[97, 1001]*0.75
  abline(h = vol75per, cex = 0.5, type = "l", lty = 2)
  text(1942, vol75per-(vol75per*.035), label = "75 percent lake volume", cex = 0.8, adj = 0)
  
  # The two below line are not needed, but uncommented if a line is to be drawn at the
  #   elevation instead of volume
  
  # Lake is 12 feet deep
  #vol12ft <- 6.52*(12)^1.8  
  #abline(h = vol12ft, cex = 0.5, type = "l", lty = 2)
  #text(2085, vol12ft-(vol12ft*.05), label = "12 feet deep", cex = 0.6)
  
  # Lake is 14 feet deep
  #vol14ft <- 6.52*(14)^1.8  
  #abline(h = vol14ft, cex = 0.5, type = "l", lty = 2)
  #text(2085, vol14ft-(vol14ft*.05), label = "14 feet deep", cex = 0.6)
  
  dev.off()

  #------------------------------------------------------------------------------------------------#
  ##  when done, there will be a dataframe (bigmalbootvols) with the simulated volumes for all       ##
  ###  nboot simulations and another dataframe (bigmalbootqntann) with the annual exceedance levels ###
  ##   for the end of each water year                                                             ##
  #------------------------------------------------------------------------------------------------#
  
  #################################################################################
  # The below plots are not used in the report, but are useful when adapting each
  #   script to the individual lake. They can be exported if needed.
  #################################################################################

  
  # #  Creates capacity curve with adjusted a and b 
  # 
  # plot(xxtmp,yytmp,type="l",xlim=c(1820, 1840),ylim=c(0,1200),ylab="Lake volume (mcf)",xlab="Lake elevation (ft)")
  # points(eft,vmcf)
  # axx <- as.character(round(evcoefs[1],2)); bxx <- as.character(round(evcoefs[2],2))
  # mtext(side=3,line=0,paste("Capacity curve: a=",axx,", b=",bxx))
  # summary(tmpout)
  # 
  # #  Plots estimated vs. observed annual volume changes
  # 
  # xxlow <- floor(min(c(estvoldif,obsvoldif))/50)*50
  # xxup <- ceiling(max(c(estvoldif,obsvoldif))/50)*50
  # plot(estvoldif,obsvoldif,xlim=c(xxlow,xxup),ylim=c(xxlow,xxup))
  # abline(0,1)
  # mtext(side=3,line=0,cex=0.8,paste("barea = ",bch1,", maxdef = ",bch2,", rho = ",rho,", bias = ",bias,", rmse = ",rmse,sep=""))
  # mtext(side=3,line=1,"Estimated versus observed annual volume changes")
  # 
  # #  Plots difference in observed to estimated volume
  # 
  # plot(1992:2016,obsvoldif-estvoldif)
  # mtext(side=3,line=1,"Water balance model residuals")
  # summary(tmpout)
  # 
  # #  Plots various charts based on quarterly time series
  # 
  # wyr <- 1992:2016
  # nyrs <- 25
  # xtmp <- watbalxx$dqtrxx 
  # ytmp <- watbalxx$qin
  # plot(xtmp,ytmp,type="o",xlim=c(1990,2020),cex=0.6,pch=16,ylab="Inflow (MCF)")
  # points(xtmp[seq(1,nyrs*4,4)],ytmp[seq(1,nyrs*4,4)],pch=16,cex=.6,col=2)
  # mtext(side=3,line=1,"Estimated quarterly inflow")
  # ytmp <- watbalxx$qout
  # plot(xtmp,ytmp,type="o",xlim=c(1990,2020),cex=0.6,pch=16,ylab="Outflow (MCF)")
  # points(xtmp[seq(1,nyrs*4,4)],ytmp[seq(1,nyrs*4,4)],pch=16,cex=.6,col=2)
  # mtext(side=3,line=1,"Estimated quarterly outflow")
  # ytmp <- watbalxx$defout*12
  # ytmp <- watbalxx$evest
  # plot(xtmp,ytmp,type="o",xlim=c(1990,2020),cex=0.6,pch=16,ylab="Net lake evaporation (MCF)")
  # points(xtmp[seq(1,nyrs*4,4)],ytmp[seq(1,nyrs*4,4)],pch=16,cex=.6,col=2)
  # mtext(side=3,line=1,"Estimated quarterly net lake evaporation")
  # 
  # # Model sensitivity runs and plots
  # 
  # wyr <- 1991:2016
  # estvolchng <- bigmalwatbal(qprec,qpet,avol,barea,3.5,prmult=0.9)
  # plot(watbalxx$dqtrxx,watbalxx$vest,type="l",xlim=c(1990,2020),ylim=c(0,1200),ylab="Lake volume (MCF)")
  # points(wyr+1,avol,cex=0.8)
  # mtext(side=3,line=1,"Sensitivity to 10 pct reduction in precip",cex=0.7)
  # estvolchng <- bigmalwatbal(qprec,qpet,avol,barea,3.5,prmult=0.9,etmult=1.1)
  # plot(watbalxx$dqtrxx,watbalxx$vest,type="l",xlim=c(1990,2020),ylim=c(0,1200),ylab="Lake volume (MCF)")
  # points(wyr+1,avol,cex=0.8)
  # mtext(side=3,line=1,"Sensitivity to 10 pct reduction in precip and 10 pct increase in PET",cex=0.7)
  # 
  # summary(tmpout)
} 



###################################################################################################
# Two functions below (bigmalwatbal and bigmalwatbal1) used in the water balance model
###################################################################################################

bigmalwatbal <- function(qprec,qpet,avol,barea,defmax,cfdgq= 0,cgwleak=0,prmult=1.0,etmult=1.0,nyrs=26) {
  
  b2 <- barea ; b1 <- cgwleak ; sscap <- defmax/12 ; b3 <- cfdgq
  # nyrs <- 26
  
  #  Adjust precipitation upward to account for undercatch during snowy months
  
  qprec[seq(1,nyrs*4,4)] <- qprec[seq(1,nyrs*4,4)]*1.1
  qprec[seq(2,nyrs*4,4)] <- qprec[seq(2,nyrs*4,4)]*1.2
  qprec[seq(3,nyrs*4,4)] <- qprec[seq(3,nyrs*4,4)]*1.1
  qprec[seq(4,nyrs*4,4)] <- qprec[seq(4,nyrs*4,4)]*1.0
  
  # Apply adjustments to potential ET to approximate actual ET
  
  qpet[seq(1,nyrs*4,4)] <- qpet[seq(1,nyrs*4,4)]*0.4
  qpet[seq(2,nyrs*4,4)] <- qpet[seq(2,nyrs*4,4)]*0
  qpet[seq(3,nyrs*4,4)] <- qpet[seq(3,nyrs*4,4)]*0.3
  qpet[seq(4,nyrs*4,4)] <- qpet[seq(4,nyrs*4,4)]*.8
  
  # Multiply inputs by optional factors (used for sensitivity analysis)
  
  qprec <- qprec*prmult
  qpet <- qpet*etmult
  
  #  initialize storage variables 
  #    v0 is lake volume for the beginning of wyr 1992
  #    prec0 and pet0 are precipitation for qtr 4, wyr 1991
  
  v0 <- avol[1] ; v0save <- v0
  e0 <- bigmalelev(v0)
  a0 <- bigmalarea(e0)
  prec0 <- qprec[4] ; pet0 <- qpet[4]
  
  #  strip off first year of precip and pet data and convert to units of ft
  
  qprec <- qprec[-c(1:4)]/12
  qpet <- qpet[-c(1:4)]/12
  
  # set spill elevation, volume, and area
  
  espl <- 1834.5  #assuming elevation of 1836 strip of land separating from Barnes Lake
  vspl <- bigmalvol(espl)
  aspl <- bigmalarea(espl)
  emin <- 1820.7
  amax <- aspl
  vmax <- vspl
  
  
  # initialize quarterly water balance variables, including estimated
  #  lake volume (vest), surface inflow volume (qin),
  #  lake area (aest), lake elevation (elest), net lake evaporation (evest), groundwater
  #  loss (qgw), surface outflow volume (qout), and storage deficit (defout)
  
  ##### nyrs <- 25
  nyrs <- nyrs-1
  vest <- rep(NA,nyrs*4)
  qin <- rep(NA,nyrs*4)
  aest <- rep(NA,nyrs*4)
  elest <- rep(NA,nyrs*4)
  evest <- rep(NA,nyrs*4)
  qgw <- rep(NA,nyrs*4)
  qout <- rep(NA,nyrs*4)
  defout <- rep(NA,nyrs*4)
  
  #  Ready to compute water balance variables
  #   i is the year, j is the quarter, and k is the cumulative time index
  #   petlag is the left-over storage deficit carried from one quarter to the next
  #   preclag is the residual fall/winter (quarters 1 and 2) precipitation that
  #     remains in frozen storage until the following spring (qtr 3)
  
  #  Initialize petlag and preclag for the first quarter (OND of WY1992)
  
  petlag <- sscap ; preclag <- 0
  qinlag <- rep(NA,nyrs*4)
  storx <- 0
  maxstorx <- 400
  cgainx <-  0.1
  cleakx <- 0.0
  cgwloss <- 0.1
  vgwloss <- 750
  
  for (i in 1:nyrs)  {
    for (j in 1:4)  {
      
      k <- (i-1)*4 + j
      
      #  quarter 1 (OND)
      #  no surface inflow or outflow
      #  carry over excess precipitation, if any
      
      if (j==1)  {
        
        evest[k] <- (1.5*qpet[k]-qprec[k])*bigmalarea(e0+5)
        evest[k] <- min(evest[k],v0)
        defout[k] <- petlag
        if(qprec[k] >= qpet[k]+petlag) { 
          preclag <- qprec[k]-qpet[k]-petlag
          petlag <- 0
        } else if (qprec[k]>=qpet[k]) {
          preclag <- 0
          petlag <- petlag - (qprec[k]-qpet[k]) 
        } else {
          preclag <- 0
          petlag = min(petlag + 1.0*(qpet[k] - qprec[k]),sscap)
        }
        
        qgw[k] <- cgwloss*max(v0-vgwloss,0)
        
        qin[k] <- 0
        storx <- storx-cleakx*exp(-2*(1-storx/(maxstorx+1)))*storx+qin[k]
        if(storx>maxstorx) {
          qinlag[k] <- (storx-maxstorx) + cgainx*maxstorx
          storx <- storx-qinlag[k]
        } else {
          qinlag[k] <- cgainx*exp(-2*(1-storx/(maxstorx+1)))*storx
          storx <- storx-qinlag[k]
        }  
        qin[k] <-  qinlag[k]
        qout[k] <- 0
        vest[k] <- max(v0 + qin[k] - evest[k] - qgw[k] - qout[k],0) 
        qout[k] <- qout[k]+qgw[k]
        elest[k] <- bigmalelev(vest[k]) ; aest[k] <- bigmalarea(elest[k])
        v0 <- vest[k]; e0 <- bigmalelev(v0); a0 <- bigmalarea(e0)
      }
      
      #  quarter 2 (JFM)
      #  no surface inflow or outflow
      #  no ET, no change in deficit 
      #  carry over precipitation
      
      if (j==2)  {
        
        evest[k] <- (1.5*qpet[k]-qprec[k])*bigmalarea(e0+5)
        evest[k] <- min(evest[k],v0)
        defout[k] <- petlag
        preclag <- preclag+qprec[k]
        
        qgw[k] <- cgwloss*max(v0-vgwloss,0)
        qin[k] <- 0
        
        storx <- storx-cleakx*exp(-2*(1-storx/(maxstorx+1)))*storx+qin[k]
        if(storx>maxstorx) {
          qinlag[k] <-(storx-maxstorx) + cgainx*maxstorx
          storx <- storx-qinlag[k]
        } else {
          qinlag[k] <- cgainx*exp(-2*(1-storx/(maxstorx+1)))*storx
          storx <- storx-qinlag[k]
        }  
        
        qin[k] <-  qinlag[k]
        qout[k] <- 0
        vest[k] <- max(v0 + qin[k] - evest[k] - qgw[k] - qout[k],0) 
        qout[k] <- qout[k]+qgw[k]
        elest[k] <- bigmalelev(vest[k]) ; aest[k] <- bigmalarea(elest[k])
        v0 <- vest[k]; e0 <- bigmalelev(v0); a0 <- bigmalarea(e0)
        
      }
      
      #  quarter 3 (AMJ)
      #  compute surface inflow and outflow
      #  update storage deficit (petlag)
      #  Note that computed outflow us one-half of the volume above the spill elevation.
      #  The rest of the excess volume remains until the next time step (accounts for
      #    for time it takes to raise the lake above spill elevation and for
      #    the water to flow out of the outlet)
      #  
      
      if (j==3)  {
        
        evest[k] <- (1.5*qpet[k]-qprec[k])*bigmalarea(e0+5)
        evest[k] <- min(evest[k],v0)
        qgw[k] <- cgwloss*max(v0-vgwloss,0)
        
        cumprec <- qprec[k] + preclag
        defout[k] <- petlag
        if(cumprec >= qpet[k] + petlag) { 
          qin[k] <- max(b2*amax-bigmalarea(e0+5),0)*(cumprec-qpet[k]-petlag)
          petlag <- 0
        } else if (cumprec >= qpet[k]) {
          qin[k] <- 0
          petlag <- petlag - (cumprec-qpet[k]) 
        } else {
          qin[k] <- 0
          petlag = min(petlag + 1.0*(qpet[k] - cumprec),sscap)
        }
        
        storx <- storx-cleakx*exp(-2*(1-storx/(maxstorx+1)))*storx+qin[k]
        if(storx>maxstorx) {
          qinlag[k] <- (storx-maxstorx) + cgainx*maxstorx
          storx <- storx-qinlag[k]
        } else {
          qinlag[k] <- cgainx*exp(-2*(1-storx/(maxstorx+1)))*storx
          storx <- storx-qinlag[k]
        }  
        
        qin[k] <- qinlag[k]+0*qin[k]
        qout[k] <- 0.5*max(0,v0+qin[k]-evest[k]-qgw[k]-vspl)
        vest[k] <- max(v0 + qin[k] - evest[k] - qgw[k] - qout[k],0) 
        qout[k] <- qout[k]+qgw[k]
        elest[k] <- bigmalelev(vest[k]) ; aest[k] <- bigmalarea(elest[k])
        v0 <- vest[k]; e0 <- bigmalelev(v0); a0 <- bigmalarea(e0)
        
      }
      
      #  quarter 4 (JAS)
      #  Same as quarter 3, except all of the excess volume above
      #    the spill elevation is assumed to become outflow, leaving the
      #    lake volume equal to or less than the spill volume at the end
      #    of the water year.
      
      if (j==4)  {
        
        evest[k] <- (1.0*qpet[k]-qprec[k])*bigmalarea(e0+5)
        evest[k] <- min(evest[k],v0)
        qgw[k] <- cgwloss*max(v0-vgwloss,0)
        
        cumprec <- qprec[k] 
        defout[k] <- petlag
        if(cumprec >= qpet[k] + petlag) { 
          qin[k] <- max(b2*amax-bigmalarea(e0+5),0)*(cumprec-qpet[k]-petlag)
          petlag <- 0
        } else if (cumprec >= qpet[k]) {
          qin[k] <- 0
          petlag <- petlag - (cumprec-qpet[k]) 
        } else {
          qin[k] <- 0
          petlag = min(petlag + 1.0*(qpet[k] - cumprec),sscap)
        }
        
        storx <- storx-cleakx*exp(-2*(1-storx/(maxstorx+1)))*storx+qin[k]
        if(storx>maxstorx) {
          qinlag[k] <- (storx-maxstorx) + cgainx*maxstorx
          storx <- storx-qinlag[k]
        } else {
          qinlag[k] <- cgainx*exp(-2*(1-storx/(maxstorx+1)))*storx
          storx <- storx-qinlag[k]
        }  
        
        qin[k] <- qinlag[k] +0*qin[k]
        qout[k] <- 1.0*max(0,v0+qin[k]-evest[k]-qgw[k]-vspl)
        vest[k] <- max(v0 + qin[k] - evest[k] - qgw[k] - qout[k],0)
        qout[k] <- qout[k]+qgw[k]
        elest[k] <- bigmalelev(vest[k]) ; aest[k] <- bigmalarea(elest[k])
        v0 <- vest[k]; e0 <- bigmalelev(v0); a0 <- bigmalarea(e0)
        
      }
      
      
    }
  }
  
  #  For model diagnostic plots, save the quarterly water balance variables in a
  #  dataframe called watbalxx. Note the double assignment arrow, which means the
  #  dataframe is saved in the workspace when the function is done.
  
  dqtrxx <- 1992+seq(0,nyrs*4-1,1)/4
  yrxx <- floor(dqtrxx)
  qtrxx <- (dqtrxx-floor(dqtrxx))*4 +1
  watbalxx <<- data.frame(yrxx,qtrxx,dqtrxx,elest,aest,vest,qin,qout,evest,qgw,defout)
  
  # Function output consists of the fitted annual lake volume changes for each water year
  # print(qinaccum)
  # plot(dqtrxx,qinlag,type="l") ; abline(h=600) ; abline(h=1200)
  volannual <- c(v0save,vest[seq(4,nyrs*4,4)])
  volannual <- volannual[-1]  #  volannual[-1] - volannual[-(nyrs+1)]
  volannual
  
}

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

bigmalwatbal1 <- function(qprec,qpet,avol,barea,defmax,cfdgq= 0,cgwleak=0,prmult=1.0,etmult=1.0,nyrs=26) {
  
  b2 <- barea ; b1 <- cgwleak ; sscap <- defmax/12 ; b3 <- cfdgq
  # nyrs <- 26
  
  #  Adjust precipitation upward to account for undercatch during snowy months
  
  qprec[seq(1,nyrs*4,4)] <- qprec[seq(1,nyrs*4,4)]*1.1
  qprec[seq(2,nyrs*4,4)] <- qprec[seq(2,nyrs*4,4)]*1.2
  qprec[seq(3,nyrs*4,4)] <- qprec[seq(3,nyrs*4,4)]*1.1
  qprec[seq(4,nyrs*4,4)] <- qprec[seq(4,nyrs*4,4)]*1.0
  
  # Apply adjustments to potential ET to approximate actual ET
  
  qpet[seq(1,nyrs*4,4)] <- qpet[seq(1,nyrs*4,4)]*0.4
  qpet[seq(2,nyrs*4,4)] <- qpet[seq(2,nyrs*4,4)]*0
  qpet[seq(3,nyrs*4,4)] <- qpet[seq(3,nyrs*4,4)]*0.3
  qpet[seq(4,nyrs*4,4)] <- qpet[seq(4,nyrs*4,4)]*.8
  
  # Multiply inputs by optional factors (used for sensitivity analysis)
  
  qprec <- qprec*prmult
  qpet <- qpet*etmult
  
  #  initialize storage variables 
  #    v0 is lake volume for the beginning of wyr 1992
  #    prec0 and pet0 are precipitation for qtr 4, wyr 1991
  
  v0 <- avol[1] ; v0save <- v0
  e0 <- bigmalelev(v0)
  a0 <- bigmalarea(e0)
  prec0 <- qprec[4] ; pet0 <- qpet[4]
  
  #  strip off first year of precip and pet data and convert to units of ft
  
  qprec <- qprec[-c(1:4)]/12
  qpet <- qpet[-c(1:4)]/12
  
  # set spill elevation, volume, and area
  
  espl <- 1834.5  #assuming elevation of 1836 strip of land separating from Barnes Lake
  vspl <- bigmalvol(espl)
  aspl <- bigmalarea(espl)
  emin <- 1820.7
  amax <- aspl
  vmax <- vspl
  
  # initialize quarterly water balance variables, including estimated
  #  lake volume (vest), surface inflow volume (qin),
  #  lake area (aest), lake elevation (elest), net lake evaporation (evest), groundwater
  #  loss (qgw), surface outflow volume (qout), and storage deficit (defout)
  
  ##### nyrs <- 25
  nyrs <- nyrs-1
  vest <- rep(NA,nyrs*4)
  qin <- rep(NA,nyrs*4)
  aest <- rep(NA,nyrs*4)
  elest <- rep(NA,nyrs*4)
  evest <- rep(NA,nyrs*4)
  qgw <- rep(NA,nyrs*4)
  qout <- rep(NA,nyrs*4)
  defout <- rep(NA,nyrs*4)
  
  #  Ready to compute water balance variables
  #   i is the year, j is the quarter, and k is the cumulative time index
  #   petlag is the left-over storage deficit carried from one quarter to the next
  #   preclag is the residual fall/winter (quarters 1 and 2) precipitation that
  #     remains in frozen storage until the following spring (qtr 3)
  
  #  Initialize petlag and preclag for the first quarter (OND of WY1992)
  
  petlag <- sscap ; preclag <- 0
  qinlag <- rep(NA,nyrs*4)
  storx <- 0
  maxstorx <- 400
  cgainx <-  0.1
  cleakx <- 0.0
  cgwloss <- 0.1
  vgwloss <- 750
  
  for (i in 1:nyrs)  {
    for (j in 1:4)  {
      
      k <- (i-1)*4 + j
      
      #  quarter 1 (OND)
      #  no surface inflow or outflow
      #  carry over excess precipitation, if any
      
      if (j==1)  {
        
        evest[k] <- (1.5*qpet[k]-qprec[k])*bigmalarea(e0+5)
        evest[k] <- min(evest[k],v0)
        defout[k] <- petlag
        if(qprec[k] >= qpet[k]+petlag) { 
          preclag <- qprec[k]-qpet[k]-petlag
          petlag <- 0
        } else if (qprec[k]>=qpet[k]) {
          preclag <- 0
          petlag <- petlag - (qprec[k]-qpet[k]) 
        } else {
          preclag <- 0
          petlag = min(petlag + 1.0*(qpet[k] - qprec[k]),sscap)
        }
        
        qgw[k] <- cgwloss*max(v0-vgwloss,0)
        
        qin[k] <- 0
        storx <- storx-cleakx*exp(-2*(1-storx/(maxstorx+1)))*storx+qin[k]
        if(storx>maxstorx) {
          qinlag[k] <- (storx-maxstorx) + cgainx*maxstorx
          storx <- storx-qinlag[k]
        } else {
          qinlag[k] <- cgainx*exp(-2*(1-storx/(maxstorx+1)))*storx
          storx <- storx-qinlag[k]
        }  
        qin[k] <-  qinlag[k]
        qout[k] <- 0
        vest[k] <- max(v0 + qin[k] - evest[k] - qgw[k] - qout[k],0) 
        qout[k] <- qout[k]+qgw[k]
        elest[k] <- bigmalelev(vest[k]) ; aest[k] <- bigmalarea(elest[k])
        v0 <- vest[k]; e0 <- bigmalelev(v0); a0 <- bigmalarea(e0)
      }
      
      #  quarter 2 (JFM)
      #  no surface inflow or outflow
      #  no ET, no change in deficit 
      #  carry over precipitation
      
      if (j==2)  {
        
        evest[k] <- (1.5*qpet[k]-qprec[k])*bigmalarea(e0+5)
        evest[k] <- min(evest[k],v0)
        defout[k] <- petlag
        preclag <- preclag+qprec[k]
        
        qgw[k] <- cgwloss*max(v0-vgwloss,0)
        qin[k] <- 0
        
        storx <- storx-cleakx*exp(-2*(1-storx/(maxstorx+1)))*storx+qin[k]
        if(storx>maxstorx) {
          qinlag[k] <-(storx-maxstorx) + cgainx*maxstorx
          storx <- storx-qinlag[k]
        } else {
          qinlag[k] <- cgainx*exp(-2*(1-storx/(maxstorx+1)))*storx
          storx <- storx-qinlag[k]
        }  
        
        qin[k] <-  qinlag[k]
        qout[k] <- 0
        vest[k] <- max(v0 + qin[k] - evest[k] - qgw[k] - qout[k],0) 
        qout[k] <- qout[k]+qgw[k]
        elest[k] <- bigmalelev(vest[k]) ; aest[k] <- bigmalarea(elest[k])
        v0 <- vest[k]; e0 <- bigmalelev(v0); a0 <- bigmalarea(e0)
        
      }
      
      #  quarter 3 (AMJ)
      #  compute surface inflow and outflow
      #  update storage deficit (petlag)
      #  Note that computed outflow us one-half of the volume above the spill elevation.
      #  The rest of the excess volume remains until the next time step (accounts for
      #    for time it takes to raise the lake above spill elevation and for
      #    the water to flow out of the outlet)
      #  
      
      if (j==3)  {
        
        evest[k] <- (1.5*qpet[k]-qprec[k])*bigmalarea(e0+5)
        evest[k] <- min(evest[k],v0)
        qgw[k] <- cgwloss*max(v0-vgwloss,0)
        
        cumprec <- qprec[k] + preclag
        defout[k] <- petlag
        if(cumprec >= qpet[k] + petlag) { 
          qin[k] <- max(b2*amax-bigmalarea(e0+5),0)*(cumprec-qpet[k]-petlag)
          petlag <- 0
        } else if (cumprec >= qpet[k]) {
          qin[k] <- 0
          petlag <- petlag - (cumprec-qpet[k]) 
        } else {
          qin[k] <- 0
          petlag = min(petlag + 1.0*(qpet[k] - cumprec),sscap)
        }
        
        storx <- storx-cleakx*exp(-2*(1-storx/(maxstorx+1)))*storx+qin[k]
        if(storx>maxstorx) {
          qinlag[k] <- (storx-maxstorx) + cgainx*maxstorx
          storx <- storx-qinlag[k]
        } else {
          qinlag[k] <- cgainx*exp(-2*(1-storx/(maxstorx+1)))*storx
          storx <- storx-qinlag[k]
        }  
        
        qin[k] <- qinlag[k]+0*qin[k]
        qout[k] <- 0.5*max(0,v0+qin[k]-evest[k]-qgw[k]-vspl)
        vest[k] <- max(v0 + qin[k] - evest[k] - qgw[k] - qout[k],0) 
        qout[k] <- qout[k]+qgw[k]
        elest[k] <- bigmalelev(vest[k]) ; aest[k] <- bigmalarea(elest[k])
        v0 <- vest[k]; e0 <- bigmalelev(v0); a0 <- bigmalarea(e0)
        
      }
      
      #  quarter 4 (JAS)
      #  Same as quarter 3, except all of the excess volume above
      #    the spill elevation is assumed to become outflow, leaving the
      #    lake volume equal to or less than the spill volume at the end
      #    of the water year.
      
      if (j==4)  {
        
        evest[k] <- (1.0*qpet[k]-qprec[k])*bigmalarea(e0+5)
        evest[k] <- min(evest[k],v0)
        qgw[k] <- cgwloss*max(v0-vgwloss,0)
        
        cumprec <- qprec[k] 
        defout[k] <- petlag
        if(cumprec >= qpet[k] + petlag) { 
          qin[k] <- max(b2*amax-bigmalarea(e0+5),0)*(cumprec-qpet[k]-petlag)
          petlag <- 0
        } else if (cumprec >= qpet[k]) {
          qin[k] <- 0
          petlag <- petlag - (cumprec-qpet[k]) 
        } else {
          qin[k] <- 0
          petlag = min(petlag + 1.0*(qpet[k] - cumprec),sscap)
        }
        
        storx <- storx-cleakx*exp(-2*(1-storx/(maxstorx+1)))*storx+qin[k]
        if(storx>maxstorx) {
          qinlag[k] <- (storx-maxstorx) + cgainx*maxstorx
          storx <- storx-qinlag[k]
        } else {
          qinlag[k] <- cgainx*exp(-2*(1-storx/(maxstorx+1)))*storx
          storx <- storx-qinlag[k]
        }  
        
        qin[k] <- qinlag[k] +0*qin[k]
        qout[k] <- 1.0*max(0,v0+qin[k]-evest[k]-qgw[k]-vspl)
        vest[k] <- max(v0 + qin[k] - evest[k] - qgw[k] - qout[k],0)
        qout[k] <- qout[k]+qgw[k]
        elest[k] <- bigmalelev(vest[k]) ; aest[k] <- bigmalarea(elest[k])
        v0 <- vest[k]; e0 <- bigmalelev(v0); a0 <- bigmalarea(e0)
        
      }
      
    }
  }
  
  #  For model diagnostic plots, save the quarterly water balance variables in a
  #  dataframe called watbalxx. Note the double assignment arrow, which means the
  #  dataframe is saved in the workspace when the function is done.
  
  dqtrxx <- 1896+seq(0,nyrs*4-1,1)/4
  yrxx <- floor(dqtrxx)
  qtrxx <- (dqtrxx-floor(dqtrxx))*4 +1
  watbalxx <<- data.frame(yrxx,qtrxx,dqtrxx,elest,aest,vest,qin,qout,evest,qgw,defout)
  
  # Function output consists of the fitted annual lake volume changes for each water year

  volannual <- c(v0save,vest[seq(4,nyrs*4,4)])
  volannual <- volannual[-1]
  volannual
  
}

#################################################################################
###  These are the functions for computing lake volumes and areas for
###  Big Mallard Marsh
#################################################################################

bigmalvol <- function(elev) {
  #  computes lake volume (million cubic ft) given elevation
  elev <- max(elev,1820.7)
  vol <- 6.52*(elev-1820.7)^1.8
  vol 
}

bigmalarea <- function(elev) {
  #  computes lake area (million sq. ft) given elevation
  elev <- max(elev,1820.7)
  area <- 11.74*(elev-1820.7)^0.8
  area
}

bigmalelev <- function(vol) {
  # computes lake elevation given volume 
  vol <- max(vol,0)
  elev=1820.7+(vol/6.52)^(1/1.8)
  elev
}