#  File: /projects/dhsvm/uvm/scripts/plot_dhsvm1.r
#
#  source("/projects/dhsvm/uvm/scripts/plot_dhsvm1.r")  # PNNL
#  source("/users/b/w/bwemple/dhsvm1/scripts/plot_dhsvm1.r")  # UVM
#
#  R script to plot selected DHSVM output variables, and prepare tables.
#
#  The following R library packages must be installed for this script to work:
#  date, chron.
#
#  NOTES ON STATE OF CODE AND FURTHER DEVELOPMENT:
#  
#
#  Scott Waichler, 09/30/03
#  scott.waichler@pnl.gov
#  Battelle Pacific Northwest National Laboratory
#  "Serving the needs of the nation and humanity.
#   Putting technology to work."

# Set one of these to true (T), and the other to false
UVM <- T
PNNL <- F

doinit <- T  # initialization
doplot1 <- T  # plot HOURLY streamflow vs. time, both simulated and observed, 1 water year per graph,
              # and 1:1 plot
doplot2 <- T  # plot DAILY streamflow vs. time, both simulated and observed, 1 water year per graph,
              # and 1:1 plot
doplot3 <- T  # plot DAILY variables vs. time, both simulated and observed, 1 water year per graph
doplot4 <- T  # plot MONTHLY water balance

if(PNNL) {
  source("/home/waichler/R/functions1.R")  # PNNL, misc. functions
  source("/home/waichler/R/get_tickpoints.r")  # PNNL, misc. functions
  output.directory <- "/projects/dhsvm/uvm/output/"  # PNNL
  graphics.directory <- "/projects/dhsvm/uvm/graphics/ranch"  # PNNL
  obsfilen <- "/projects/dhsvm/uvm/data/flow-ranch.txt"
} else if (UVM) {
  source("/users/b/w/bwemple/dhsvm1/scripts/functions1.r")  # UVM, misc. functions
  source("/users/b/w/bwemple/dhsvm1/scripts/get_tickpoints.r")  # UVM, misc. functions
  output.directory <- "/users/b/w/bwemple/dhsvm1/outputmans/"  # UVM
  graphics.directory <- "/users/b/w/bwemple/dhsvm1/graphics/"  # UVM
  obsfilen <- "/users/b/w/bwemple/dhsvm1/flow-ranch.txt"
}

library(date)
library(chron)
 
if(doinit) {
  
  # Read in the data (observed) file
  obs <- scan(file=obsfilen, list(date="", time="", flow=0), skip=1)
  obs$datetime <- chron(dates=dates(obs$date), times=times(paste(sep="", obs$time, ":00")))
  num.cells <- 10639  # number of DHSVM cells in watershed
  cell.area <- 900    # area per DHSVM grid cell (m2)
  # watershed conversion factor to go from streamflow in cfs to mm/hr
  millimeter.factor <- (0.0283196 * 1000 * 3600) / (num.cells * cell.area)
  obs$flow <- obs$flow * millimeter.factor
  cat(file="", "Finished reading obs file", obsfilen, "\n")

  # Read in the DHSVM Stream.Flow output file
  stream.file <- paste(sep="", output.directory, "Stream.Flow")
  input <- scan(stream.file, fill=T,
                list(dt="", dum1=0, dum2=0, dum3=0, flow=0, dum4=0, gage="", dum5=""),
                skip=0)
  input$flow <- (input$flow * 1000) / (num.cells * cell.area)  # convert from m3/hr to mm/hr
  gages <- c("Gage1") # list the names of the gage(s) that are output in Stream.Flow
  num.years <- 1  # number of years in target period
  stream <- list()  # initialize a list to hold the streamflow data for gage(s)

  for(i in 1:length(gages)) {
    stream[[i]] <- list()
    stream[[i]]$ind <- which(input$gage == gages[i])
    stream[[i]]$hourly.flow <- input$flow[stream[[i]]$ind]
    stream[[i]]$daily.flow <- as.vector(aggregate.ts(input$flow[stream[[i]]$ind], ndeltat = 24, FUN=sum))
    stream[[i]]$mean.annual.flow <- sum(stream[[i]]$flow) / num.years
  }
  cat(file="", "Finished reading the Stream.Flow file \n")

  # Read in the DHSVM Aggregated.Values output file
  # First read in the complete list of the variables in that file,
  # then construct and execute the scan() command.
  output.list.file <- paste(sep="", output.directory, "Output.list")
  aggvars <- scan(file=output.list.file, quiet=T, list(name="", ind=0), skip=1)
  aggvars$ind <- aggvars$ind + 1 # R indexing starts with 1, not 0
  number.vars <- max(aggvars$ind)
  # Change variable names by subsituting . for underscore _ so R isn't confused.
  for (i in 2:number.vars) {
    aggvars$name[i] <- gsub("_", ".", aggvars$name[i], ignore.case=FALSE, extended=TRUE)
  }
  # Construct the command to read in the output file
  output.agg.file <- paste(sep="", output.directory, "Aggregated.Values")
  read.agg.cmd <- paste(sep="", 'scan(file="', output.agg.file, '", list(datedashtime=""')
  for (i in 2:number.vars) read.agg.cmd <- paste(sep="", read.agg.cmd,  ", ", aggvars$name[i], "=0")
  read.agg.cmd <- paste(read.agg.cmd,  "), skip=0)")
  sim <- eval(parse(text=read.agg.cmd))  # execute the read command
  sim$date <- substr(sim$datedashtime, 1, 10)
  sim$time <- substr(sim$datedashtime, 12, 20) 
  sim$datetime <- chron(sim$date, sim$time, format=c("m/d/y", "h:m:s"))
#  sim$datetime <- sim$datetime + 0.5/24.0/3600.0 # add a half-second to round to whole hours
  sim$flow <- stream[[1]]$hourly.flow  # put the streamflow in the same list as other variables
  
}  # end doinit


# Plot streamflow timeseries and 1:1 graph
if(doplot1) {
 
#  Uncomment this section to make a Postscript file
# ps.options(onefile=F, print.it=F, paper="special", width=8.5, height=5.0,
#             horizontal=F, pointsize=12, family="Helvetica")
#  pfilen <- paste(sep="", graphics.directory, "_streamf_ts.ps")  # streamflow timeseries file
#  postscript(file=pfilen)
  
#  Uncomment this section to make a PNG file
  pfilen <- paste(sep="", graphics.directory, "ranch_streamf_ts.png")  # streamflow timeseries file#
  bitmap(pfilen, type = "png256", height = 5, width = 5, res = 200)
  
  par(omi=c(0,0,0,0))
  par(mar=c(2, 2, 1, 1))

  # get start and end dates for each year 
  begindate <- "10/01/2001"
#  begindatetime <- chron(dates=dates(begindate), times=times("00:00:00"))
  begindatetime <- chron(dates=dates(begindate), times=times("01:00:00"))
#  enddate <- "09/30/2002"
  enddate <- "09/10/2002"
  enddatetime <- chron(dates=dates(enddate), times=times("23:59:59"))
  x.ticks <- as.chron(seq.dates(begindatetime, enddatetime, by="months")) # monthly spacing
  x.tick.labels <- substr(as.character(dates(x.ticks)), 1, 5)
  x.tick.labels <- c("Oct1", "Nov1", "Dec1", "Jan1", "Feb1", "Mar1",  "Apr1", "May1",
                     "Jun1", "Jul1", "Aug1", "Sep1")

  # Get the indices for the right dates
  obs.ind <- which(obs$datetime >= begindatetime & obs$datetime <= enddatetime)
  sim.ind <- which(sim$datetime >= begindatetime & sim$datetime <= enddatetime)

  hydrograph.datetime <- obs$datetime[obs.ind]
  
  lw <- 1.0
  size <- 1.0
  obs1.color <- "blue"
  sim1.color <- "magenta"
  obs2.color <- "blue"
  sim2.color <- "red"        
  ymin1 <- min(obs$flow[obs.ind], sim$flow[sim.ind])
  ymax1 <- max(obs$flow[obs.ind], sim$flow[sim.ind])
  xlabel <- "WY 2002" 
  wy <- 2002
  
  plot(hydrograph.datetime, obs$flow[obs.ind],
       type="l", lty=1, col=obs2.color, lwd=lw, ylim=c(0, ymax1), xaxt="n", 
       xlim=c(begindatetime, enddatetime), ylab = "Flow (mm/hr)", xlab="")
  axis(1, x.ticks, x.tick.labels) 
  lines(hydrograph.datetime, sim$flow[sim.ind], type="l", lty=2, lwd=lw, col=sim1.color)
  lines(hydrograph.datetime, 1000*sim[["ChannelInt"]][sim.ind], type="l", lty=2, lwd=lw, col=sim2.color)

  note.x.origin <- paste("04/01/", as.character(wy), sep="")  
  note.x.origin <- chron(dates=(note.x.origin), times=times("00:00:00"))
  text(note.x.origin, ymax1, xlabel, cex=1.0)

  legend.x.origin <- paste("06/01/", as.character(wy), sep="")  
  legend.x.origin <- chron(dates=(legend.x.origin), times=times("00:00:00"))
  legend(legend.x.origin, 0.9*ymax1,
     c("Observed", "Simulated", "Channel Int"),
     col=c(obs1.color, sim1.color, sim2.color), bty="n",lty=c(1, 2, 2), lwd=lw, cex=size)

  #note <- paste(sep="", "E=", stats[1], ", E1prime=", stats[2], ", d1prime=", stats[3])        
  #note.x.origin <- paste("06/01/", as.character(wy), sep="")  
  #note.x.origin <- chron(dates=(note.x.origin), times=times("00:00:00"))
  #text(note.x.origin, 0.5*ymax1, note, cex=size, adj=0)

  dev.off()


  # 1:1 plot of same data as above
  ps.options(onefile=F, print.it=F, paper="special", width=5.0, height=5.0,
             horizontal=F, pointsize=12, family="Helvetica")
  pfilen <- paste(sep="", graphics.directory, "ranch_streamf_1to1.ps")  # streamflow timeseries file
  postscript(file=pfilen)
  par(pty="s")    
  par(omi=c(0,0,0,0))
  par(mar=c(4, 4, 0.2, 0.2))

  lw <- 1.0
  size <- 1.0
  obs1.color <- "black"
  sim1.color <- "black"
  obs2.color <- "black"
  sim2.color <- "black"        
  ymin1 <- min(obs$flow[obs.ind], sim$flow[sim.ind])
  ymax1 <- max(obs$flow[obs.ind], sim$flow[sim.ind])
  xlabel <- "WY 2002" 

  plot(obs$flow[obs.ind], sim$flow[sim.ind],
       type="p", col=obs2.color, pch=".", ylim=c(0, ymax1), 
       xlim=c(0, ymax1), ylab = "Simulated Flow (mm/hr)",
       xlab="Observed Flow (mm/hr)")
  abline(0,1)

  dev.off()       

}  # end doplot1


# Plot DAILY streamflow timeseries and 1:1 graph
if(doplot2) {
 
#  Uncomment this section to make a Postscript file
#  ps.options(onefile=F, print.it=F, paper="special", width=8.5, height=5.0,
#             horizontal=F, pointsize=12, family="Helvetica")
#  pfilen <- paste(sep="", graphics.directory, "_streamf_daily.ps")  # streamflow timeseries file
#  postscript(file=pfilen)
  
#  Uncomment this section to make a PNG file
  pfilen <- paste(sep="", graphics.directory, "ranch_streamf_daily.png")  # streamflow timeseries file
  bitmap(pfilen, type = "png256", height = 5, width = 5, res = 200)
    
  par(omi=c(0,0,0,0))
  par(mar=c(2, 2, 1, 1))

  # get start and end dates for each year 
  begindate <- "10/02/2001"
#  begindatetime <- chron(dates=dates(begindate), times=times("00:00:00"))
  begindatetime <- chron(dates=dates(begindate), times=times("00:00:00"))
#  enddate <- "09/30/2002"
  enddate <- "09/10/2002"
  enddatetime <- chron(dates=dates(enddate), times=times("23:59:59"))
  x.ticks <- as.chron(seq.dates(begindatetime, enddatetime, by="months")) # monthly spacing
  x.tick.labels <- substr(as.character(dates(x.ticks)), 1, 5)
  x.tick.labels <- c("Oct1", "Nov1", "Dec1", "Jan1", "Feb1", "Mar1",  "Apr1", "May1",
                     "Jun1", "Jul1", "Aug1", "Sep1")

  # Get the indices for the right dates
  obs.ind <- which(obs$datetime >= begindatetime & obs$datetime <= enddatetime)
  sim.ind <- which(sim$datetime >= begindatetime & sim$datetime <= enddatetime)

  # aggregate to daily
  sim.flow <- as.vector(aggregate.ts(sim$flow[sim.ind], ndeltat = 24, FUN=sum))
  sim.channelint <- as.vector(aggregate.ts(1000*sim[["ChannelInt"]][sim.ind], ndeltat = 24, FUN=sum))
  obs.flow <- as.vector(aggregate.ts(obs$flow[obs.ind], ndeltat = 24, FUN=sum))
  hydrograph.datetime <- as.vector(aggregate.ts(obs$datetime[obs.ind], ndeltat = 24, FUN=mean))
  
  lw <- 1.0
  size <- 1.0
  obs1.color <- "blue"
  sim1.color <- "red"
  obs2.color <- "blue"
  sim2.color <- "black"        
  ymin1 <- min(obs.flow, sim.flow)
  ymax1 <- max(obs.flow, sim.flow)
  xlabel <- "Daily Flow WY 2002" 
  wy <- 2002
  
  plot(hydrograph.datetime, obs.flow,
       type="l", lty=1, col=obs2.color, lwd=lw, ylim=c(0, ymax1), xaxt="n", 
       xlim=c(begindatetime, enddatetime), ylab = "Flow (mm/day)", xlab="")
  axis(1, x.ticks, x.tick.labels) 
  lines(hydrograph.datetime, sim.flow, type="l", lty=1, lwd=lw, col=sim1.color)
  lines(hydrograph.datetime, sim.channelint, type="l", lty=3, lwd=lw, col=sim2.color)

  note.x.origin <- paste("04/01/", as.character(wy), sep="")  
  note.x.origin <- chron(dates=(note.x.origin), times=times("00:00:00"))
  text(note.x.origin, ymax1, xlabel, cex=1.0)

  legend.x.origin <- paste("06/01/", as.character(wy), sep="")  
  legend.x.origin <- chron(dates=(legend.x.origin), times=times("00:00:00"))
  legend(legend.x.origin, 0.9*ymax1,
     c("Observed", "Simulated", "Channel Int"),
     col=c(obs1.color, sim1.color, sim2.color), bty="n",lty=c(1, 1, 3), lwd=lw, cex=size)

  #note <- paste(sep="", "E=", stats[1], ", E1prime=", stats[2], ", d1prime=", stats[3])        
  #note.x.origin <- paste("06/01/", as.character(wy), sep="")  
  #note.x.origin <- chron(dates=(note.x.origin), times=times("00:00:00"))
  #text(note.x.origin, 0.5*ymax1, note, cex=size, adj=0)

  dev.off()


  # 1:1 plot of same data as above
  ps.options(onefile=F, print.it=F, paper="special", width=5.0, height=5.0,
             horizontal=F, pointsize=12, family="Helvetica")
  pfilen <- paste(sep="", graphics.directory, "ranch_streamf_1to1.ps")  # streamflow timeseries file
  postscript(file=pfilen)
  par(pty="s")    
  par(omi=c(0,0,0,0))
  par(mar=c(4, 4, 0.2, 0.2))

  lw <- 1.0
  size <- 1.0
  obs1.color <- "black"
  sim1.color <- "black"
  obs2.color <- "black"
  sim2.color <- "black"        
  ymin1 <- min(obs.flow, sim.flow)
  ymax1 <- max(obs.flow, sim.flow)
  xlabel <- "WY 2002" 

  plot(obs.flow, sim.flow,
       type="p", col=obs2.color, pch=".", ylim=c(0, ymax1), 
       xlim=c(0, ymax1), ylab = "Simulated Flow (mm/hr)",
       xlab="Observed Flow (mm/hr)")
  abline(0,1)

  dev.off()       

}  # end doplot2


# Daily fluxes
if(doplot3) {

#  Uncomment this section to make a Postscript file
#  ps.options(onefile=F, print.it=F, paper="special", width=8.5, height=5.0,
#             horizontal=F, pointsize=12, family="Helvetica")
#  ps.options(onefile=F, print.it=F, paper="letter", width=11.0, height=8.5,
#             horizontal=T, pointsize=12, family="Helvetica")
#  pfilen <- paste(sep="", graphics.directory, "_daily_vars.ps") # streamflow timeseries file
#  postscript(file=pfilen)
  
#  Uncomment this section to make a PNG file
   pfilen <- paste(sep="", graphics.directory, "ranch_daily_vars.png")  # streamflow timeseries file
   bitmap(pfilen, type = "png256", height = 5, width = 8.5, res = 200)

  par(pty="m") # take up the maximum size
  par(mfrow=c(2,1))  # two plots per page
  par(omi=c(0.2,0.2,0.2,0.2))
  par(mar=c(2, 3, 1, 3))
  
  # get start and end dates for each year 
  begindate <- "10/02/2001"
  begindatetime <- chron(dates=dates(begindate), times=times("00:00:00"))
  enddate <- "09/30/2002"
  enddatetime <- chron(dates=dates(enddate), times=times("23:59:59"))
  x.ticks <- as.chron(seq.dates(begindatetime, enddatetime, by="months")) # monthly spacing
  x.tick.labels <- substr(as.character(dates(x.ticks)), 1, 5)
  x.tick.labels <- c("Oct1", "Nov1", "Dec1", "Jan1", "Feb1", "Mar1",  "Apr1", "May1",
                     "Jun1", "Jul1", "Aug1", "Sep1")
  # Get the indices for the right dates
  obs.ind <- which(obs$datetime >= begindatetime & obs$datetime <= enddatetime)
  sim.ind <- which(sim$datetime >= begindatetime & sim$datetime <= enddatetime)

  daily <- list()
  daily$datetime <- as.vector(aggregate.ts(obs$datetime[obs.ind], ndeltat = 24, FUN=mean)) 
  daily$precip <- as.vector(aggregate.ts(1000*sim[["Precip"]][sim.ind], ndeltat = 24, FUN=sum))
  daily$obs.stream <- as.vector(aggregate.ts(obs$flow[obs.ind], ndeltat = 24, FUN=sum))
  daily$stream <- as.vector(aggregate.ts(sim$flow[sim.ind], ndeltat = 24, FUN=sum))
  daily$gwrecharge <- as.vector(aggregate.ts(1000*sim[["GwRecharge"]][sim.ind], ndeltat = 24, FUN=sum))
  daily$tair <- as.vector(aggregate.ts(sim[["MeanTair"]][sim.ind], ndeltat = 24, FUN=mean))
  daily$swq <- as.vector(aggregate.ts(1000*sim[["Swq"]][sim.ind], ndeltat = 24, FUN=mean))
  daily$et <- as.vector(aggregate.ts(1000*sim[["ETot"]][sim.ind], ndeltat = 24, FUN=sum))

  
  # First plot:  precip, swe, tair
  # left vertical axis
  y <- c(daily$precip, daily$swq, daily$et)
  num.tick.marks <- 8
  axis.min <- 0
  axis.max<- "float"
  y1.tickvalues <- GET.FIRST.AXIS(y, num.tick.marks, axis.min, axis.max)
  y1labs <- as.character(y1.tickvalues)
  miny1 <- range(y1.tickvalues)[1]
  maxy1 <- range(y1.tickvalues)[2]
#  y1.tickvalues <- seq(10, 70, by=10)

  # right vertical axis
  y2 <- daily$tair
  axis.min <- "float"
  axis.max<- "float"
  y2.tickvalues <- GET.SECOND.AXIS(y1.tickvalues, y2, axis.min, axis.max)
  y2labs <- as.character(y2.tickvalues)
  miny2 <- range(y2.tickvalues)[1]
  maxy2 <- range(y2.tickvalues)[2]
  # R only knows how to plot against the left axis scale,
  # so we trick it through scaling the variables to be plotted on the right axis.
  scaledy2 <- (y2 - miny2) / (maxy2 - miny2) 
  scaledy2 <- (scaledy2 * (maxy1 - miny1)) + miny1
    
  plot(daily$datetime, daily$precip,   # precip
         xlab = "Date",
         ylab = "Flux (mm)", axes=F, ylim=range(y1.tickvalues),
         type="h", lty=1, cex=1.2, col="blue"
      )

  lines(daily$datetime, daily$swq, lty=1, lwd=1, col="green") # Snow water equivalent
  lines(daily$datetime, scaledy2, lty=1, lwd=1, col="purple") # Air temperature

  format = c("m/d/y","h:m:s")

  axis(1, x.ticks, x.tick.labels, cex=1.2) 
  axis(2, at=y1.tickvalues, labels=y1labs, cex=1.2) 
  axis(4, at=y1.tickvalues, labels=y2labs, cex=1.2) 

  box()
  legend(begindatetime+42, max(y1.tickvalues),
         c("Precip", "Tair", "SWE"),
         lty=c(1, 1, 1),
         lwd=c(1, 1, 1),
         col=c("blue", "purple", "green"),
         bty="n", cex=1)
  
  mtext("Air Temp (C)", side=4, line=2, cex=1, outer=F)

  
  # Second plot:  streamflow and et
  # left vertical axis
  y <- c(daily$obs.stream, daily$stream, daily$et)
  num.tick.marks <- 6
  axis.min <- 0
  axis.max<- "float"
  y1.tickvalues <- GET.FIRST.AXIS(y, num.tick.marks, axis.min, axis.max)
  #y1.tickvalues <- c(0,10,20,30,40,50,60,70)  
  y1labs <- as.character(y1.tickvalues)
  # y1.tickvalues <- seq(10, 70, by=10)

  # right vertical axis
#  y2 <- daily$tair
#  axis.min <- "float"
#  axis.max<- "float"
#  y2.tickvalues <- GET.SECOND.AXIS(y1.tickvalues, y2, axis.min, axis.max)
#  y2labs <- as.character(y2.tickvalues)
#  miny2 <- range(y2.tickvalues)[1]
#  maxy2 <- range(y2.tickvalues)[2]
  # R only knows how to plot against the left axis scale,
  # so we trick it through scaling the variables to be plotted on the right axis.
#  scaledy2 <- (y2 - miny2) / (maxy2 - miny2) 
#  scaledy2 <- (scaledy2 * (maxy1 - miny1)) + miny1
    
  plot(daily$datetime, daily$obs.stream,   # precip
         xlab = "Date",
         ylab = "Flux (mm)", axes=F, ylim=range(y1.tickvalues),
         type="l", lty=1, cex=1.2, col="blue"
      )
  lines(daily$datetime, daily$stream, lty=1, lwd=1, col="red") # simulated flow
  lines(daily$datetime, daily$et, lty=1, lwd=1, col="orange") # simulated flow

  format = c("m/d/y","h:m:s")

  axis(1, x.ticks, x.tick.labels, cex=1.2) 
  axis(2, at=y1.tickvalues, labels=y1labs, cex=1.2) 
#  axis(4, at=y1.tickvalues, labels=y2labs, cex=1.2) 

  box()
  legend(begindatetime+42, max(y1.tickvalues),
         c("Obs Flow", "Sim Flow", "ET"),
         lty=c(1, 1, 1),
         lwd=c(1, 1, 1),
         col=c("blue", "red", "orange"),
         bty="n", cex=1)
  
  dev.off()
  
} # end doplot3


# Monthly water balance
if(doplot4) {

#  Uncomment this section to make a Postscript file
#   ps.options(onefile=F, print.it=F, paper="special", width=8.5, height=5.0,
#             horizontal=F, pointsize=12, family="Helvetica")
#   pfilen <- paste(sep="", graphics.directory, "_monthly_bal.ps")  # streamflow timeseries file
#   postscript(file=pfilen)
  
# Uncomment this section to make a PNG file
  pfilen <- paste(sep="", graphics.directory, "ranch_monthly_bal.png")  # streamflow timeseries file
  bitmap(pfilen, type = "png256", height = 5, width = 5, res = 200)
  
  par(omi=c(0,0,0,0))
  par(mar=c(2, 2, 1, 1))

  mwb <- list()
  mwb$month <- date.mdy(as.date(obs$date[obs.ind], order="mdy"), weekday=FALSE)$month
  mwb$precip <- aggregate(1000*sim[["Precip"]][sim.ind], list(month = mwb$month), sum)$x
  mwb$stream <- aggregate(1000*sim[["ChannelInt"]][sim.ind], list(month = mwb$month), sum)$x
  mwb$gwrecharge <- aggregate(1000*sim[["GwRecharge"]][sim.ind], list(month = mwb$month), sum)$x
  mwb$et <- aggregate(1000*sim[["ETot"]][sim.ind], list(month = mwb$month), sum)$x
  mwb$swq <- aggregate(1000*sim[["Swq"]][sim.ind], list(month = mwb$month), mean)$x
  mwb$obs.stream <- aggregate(obs$flow[obs.ind], list(month = mwb$month), sum)$x

  mwb.remainder <- 100 * cumsum(mwb$precip - mwb$stream - mwb$et - mwb$gwrecharge) /
                           cumsum(mwb$precip)
  cy.month.names <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")  
  miny1 <- 0
  maxy1 <- 300
  miny2 <- min(mwb.remainder)
  maxy2 <- max(mwb.remainder)
  y1var <- mwb$precip
  y2var <- mwb.remainder
  scaledy2var <- (y2var - miny2) / (maxy2 - miny2) 
  scaledy2var <- (scaledy2var * (maxy1 - miny1)) + miny1
  num.yticks <- 5
  y1labs <- round(seq(miny1, maxy1, length=num.yticks), 0) 
  # convert these to the y2 axis scaling 
  y2labs <- (y1labs - miny1) / (maxy1 - miny1) 
  y2labs <- (y2labs * (maxy2 - miny2)) + miny2 
  y2labs <- round(y2labs, 2)
  # make up some labels and positions
  
  plot(mwb$precip,   #precip
           type="b", col="black", lty=1, lwd=1, pch=0, ylab = "Flux (mm)", xlab="Month",
           main="Monthly Water Balance", axes=F, ylim=c(0,400))
  axis(1, at=seq(1,12,1), labels=cy.month.names)
  axis(2, at=y1labs, labels=y1labs, cex=1.2) 
  axis(4, at=y1labs, labels=y2labs, cex=1.2) 
  box()
  mtext("Cumulative difference (%)", side=4, line=2.5, cex=1.2)
  
  lines(mwb$stream, type="b", lty=1, col="red", pch=5)  # runoff
  lines(mwb$et, type="b", lty=2, col="green", pch=1)  # et
  lines(mwb$gwrecharge, type="b", lty=2, col="brown", pch=2)  # gw recharge  
  lines(mwb$obs.stream, type="b", lty=1, col="blue", pch=5)  # obs runoff
  # lines(scaledy2var, type="b", lty=3, col="black", pch=3)  # cumulative error  
  legend(7, 400,
         c("Precip", "Obs Streamflow", "Streamflow", "ET", "GW recharge"),
         col=c("black", "blue", "red", "green", "brown"),
         bty="n",lty=c(1,1,1,2,2), pch=c(0,5,5,1,5), adj=0, cex=1)


  box()
  dev.off()

}  # end doplot4

