#  File: /projects/dhsvm/uvm/scripts/make_dhsvm_met.r.
#
#  Last changed: Tue Dec  2 09:56:20 2003 by Scott Waichler <waichler@snow.pnl.gov>
#
#  source("/projects/dhsvm/uvm/scripts/make_dhsvm_met.r")
#
#  This script reads in Beverly's supplied file and outputs a DHSVM met input file.
#  Incoming longwave radiation is generated using the Stefan-Boltzmann equation with
#  a constant emissivity.   Soil temperature is set equal to air temperature.
#  Output file has no header lines, with line format similar to following:
#  05/31/1999-21 22.0  2  92   0.0 327.8 22.0 22.0 22.0 0.0002 -999.9
#
#  The following R library packages must be installed for this script to work: (none).
#
#
#  Scott Waichler, 09/29/03.   
#  scott.waichler@pnl.gov
#  Battelle Pacific Northwest National Laboratory
#  "Serving the needs of the nation and humanity.
#   Putting technology to work."

doinit <- T
doprocess <- T

library(date)
library(chron)

#####################################################################################
# 1.  PRELIMINARIES
#####################################################################################
if(doinit) {
  
  #infile1 <- "/projects/dhsvm/uvm/met/mansfield-met.txt"  # WY 2001

  #outfile1 <- "/projects/dhsvm/uvm/met/met_dhsvm_control_wy01.txt" # WY 2001
  #errfile <- "/projects/dhsvm/uvm/met/met_nonzero_sw_wy01.txt"


  #infile1 <- "mansfield-met-adj.txt"  # WY 2002

  #outfile1 <- "mansfield-met-wy02.txt" # WY 2002

  #errfile <- "mansfield-met_nonzero_sw_wy02.txt"

  infile1 <- "met-wy03-adj.txt"  # WY 2003

  outfile1 <- "mansfield-met-wy03-adj.txt" # WY 2003

  errfile <- "mansfield-met_nonzero_sw_wy03.txt"




  nodata <- -999.9
  StefanBoltz <- 5.67e-8

  # 2001
  input01 <- scan(file=infile1, list(d="", t="", tair=0, wind=0, rh=0, sw=0, precip=0), skip=1)
  num.timesteps01 <- length(input01$d)

  # Get the date-time string in the proper format for DHSVM
  date.components <- unlist(strsplit(input01$d, "/")) # split the date string, then flatten list to vector
  vec.len <- length(date.components)
  month.ind <- seq(1, vec.len, by=3)
  day.ind <- seq(2, vec.len, by=3)
  year.ind <- seq(3, vec.len, by=3)
  mymonths <- as.integer(date.components[month.ind])
  mydays <- as.integer(date.components[day.ind])
  years.char <-  date.components[year.ind]
  months.char <- ifelse(mymonths < 10,
                        paste(sep="", "0", as.character(mymonths)),
                        as.character(mymonths)
                       )
  days.char <- ifelse(mydays < 10,
                      paste(sep="", "0", as.character(mydays)),
                      as.character(mydays)
                     )
  time.components <- unlist(strsplit(input01$t, ":")) # split the date string, then flatten list to vector
  vec.len <- length(time.components)
  hour.ind <- seq(1, vec.len, by=2)
  myhours <-  as.integer(time.components[hour.ind])
  hours.char <- ifelse(myhours < 10,
                        paste(sep="", "0", as.character(myhours)),
                        as.character(myhours)
                       )
  input01$date <- paste(sep="", months.char, "/", days.char, "/", years.char)
  input01$time <- paste(sep="", hours.char, ":00:00")
  input01$datetime <- chron(dates=input01$date,  times=input01$time)
  
  # string for dhsvm file
  input01$dt <- paste(sep="", months.char, "/", days.char, "/", years.char, "-",  hours.char)
  

}  # end doinit


#####################################################################################
# 4.  DO BUILDMET PROCESSING
#####################################################################################
if(doprocess) {

  # Downward shortwave radiation.  Input file apparently has offset of 28.43 for some solar
  # radiation, so I find out the extent of the problem:
  ind.problem <- which((myhours > 21 | myhours < 5) &
                       input01$sw != 0.0
                      )
  cat("The problem datetimes for non-zero solar radiation are: \n")
  cat(file=errfile, append=F, "")
  for(i in ind.problem) {
    cat(input01$dt[i], "\n")
    cat(file=errfile, append=T, input01$dt[i], "\n")
  }

  sw.err.begin <-  chron(dates=c("10/01/2001", "04/08/2002", "06/13/2002"),
                         times=c("00:00:00",   "02:00:00",   "11:00:00"))
  sw.err.end <-  chron(dates=c("10/11/2001", "06/04/2002", "07/10/2002"),
                         times=c("09:00:00",   "16:00:00",   "12:00:00", ))

  ind <- NULL
  input01$sw.corr <- input01$sw # init
  for(i in 1:length(sw.err.begin)) {
    ind.corr <- which(input01$datetime >= sw.err.begin[i] & input01$datetime <= sw.err.end[i]) 
    input01$sw.corr[ind.corr] <- input01$sw[ind.corr] - 28.43
  }

  # Downward longwave radiation.  Assume constant emissivity value of 0.876, which is from HUBEX
  input01$lw <- round(0.876 * StefanBoltz * ((input01$tair+273.3)^4), 1)


  met01 <-  data.frame(input01$dt, input01$tair, input01$wind,
                      input01$rh, input01$sw.corr, input01$lw,
                      input01$tair, input01$tair, input01$tair,
                      input01$precip)

  write.table(met01,  file=outfile1, append=F, quote=F, row.names=F, col.names=F)

}  # end doprocess


# end of R script





