R/Evapotranspiration.R

Defines functions ET.McGuinnessBordne ET.JensenHaise ET.ChapmanAustralian ET.HargreavesSamani ET.Abtew ET.Romanenko ET.Linacre ET.Hamon ET.Turc ET.BlaneyCriddle ET.Makkink ET.SzilagyiJozsa ET.GrangerGray ET.BrutsaertStrickler ET.PriestleyTaylor ET.MattShuttleworth ET.PenmanMonteith ET.Penman ET.default ET

Documented in ET ET.Abtew ET.BlaneyCriddle ET.BrutsaertStrickler ET.ChapmanAustralian ET.default ET.GrangerGray ET.Hamon ET.HargreavesSamani ET.JensenHaise ET.Linacre ET.Makkink ET.MattShuttleworth ET.McGuinnessBordne ET.Penman ET.PenmanMonteith ET.PriestleyTaylor ET.Romanenko ET.SzilagyiJozsa ET.Turc

ET <- function(data, constants, ...) UseMethod("ET")

ET.default <- function(data, constants, crop=NULL, alpha=NULL, solar=NULL, wind=NULL, ...) {

  if (is.null(solar)) {
    if (!is.null(data$Rs)) {
      solar = "data"
    } else if (!is.null(data$n)) {
      solar = "sunshine hours"
    } else if (!is.null(data$Cd)) {
      solar = "cloud"
    } else if (!is.null(data$Precip)) {
      solar = "monthly precipitation"
    }
  }

  if (is.null(wind)) {
    if (!is.null(data$u2)|!is.null(data$uz)) {
      wind = "yes"
    } else {
      wind = "no"
    }
  }
  if ( all(any(is.null(data$RHmax),is.null(data$RHmin)),
           is.null(data$Rs),is.null(data$n),is.null(data$Cd),is.null(data$Precip),
           all(is.null(data$uz),is.null(data$u2)),
           any(is.null(data$Tmax),is.null(data$Tmin))) ) { # no data available
    stop("No ET model is suitable according to the data availability")

  } else if ( all(any(is.null(data$RHmax),is.null(data$RHmin)),
                  is.null(data$Rs),is.null(data$n),is.null(data$Cd),is.null(data$Precip)) &
              all(is.null(data$uz),is.null(data$u2)) &
              all(!is.null(data$Tmax),!is.null(data$Tmin)) ) { # Only Tmax/Tmin available
    message("No ET model specified, choose the Hargreaves-Samani model according to the data availability")
    class(data) = "HargreavesSamani"
    ET(data, constants, ...)
  } else if ( all(any(is.null(data$RHmax),is.null(data$RHmin)),all(is.null(data$uz),is.null(data$u2))) &
              any(!is.null(data$Rs),!is.null(data$n),!is.null(data$Cd),!is.null(data$Precip)) &
              all(!is.null(data$Tmax),!is.null(data$Tmin)) ) { # Tmax/Tmin & any Rs data available
    message("No ET model specified, choose the Makkink model according to the data availability")
    class(data) = "Makkink"

    ET(data, constants, solar=solar, ...)
  } else if ( all(is.null(data$uz),is.null(data$u2)) &
              any(!is.null(data$Rs),!is.null(data$n),!is.null(data$Cd),!is.null(data$Precip)) &
              all(!is.null(data$Tmax),!is.null(data$Tmin),!is.null(data$RHmax),!is.null(data$RHmin)) &
              !is.null(alpha) ) { # Tmax/Tmin & any Rs & RHmax/RHmin data available
    message("No ET model specified, choose the Priestley-Taylor model according to the data availability")
    class(data) = "PriestleyTaylor"

    ET(data, constants, solar=solar, ...)
  } else if ( any(!is.null(data$Rs),!is.null(data$n),!is.null(data$Cd),!is.null(data$Precip)) &
              all(!is.null(data$Tmax),!is.null(data$Tmin),!is.null(data$RHmax),!is.null(data$RHmin)) &
              any(!is.null(data$uz),!is.null(data$u2)) )  { # All data available & crop specified
    Flag = 1
  } else {
    "No ET model can be recommended according to the data availability"
  }

  if (exists('Flag')) {
    if (Flag == 1) { #Penman-Monteith or Penman

      if (is.null(crop) | all(crop != "short",crop != "tall")) {
        alpha=0.08
        z0=0.001
        message("No ET model specified and no valid evaporative surface specified, choose the Penman model for open-water evaporation according to the data availability")
        class(data) = "Penman"

        ET(data, constants, solar=solar, wind=wind, windfunction_ver=1948, alpha=alpha, z0=z0, ...)
      } else {

        message("No ET model specified, choose the Penman-Monteith model according to the data availability")
        class(data) = "PenmanMonteith"
        ET(data, constants, solar=solar, wind=wind, crop=crop, ...)
      }
    }
  }


}

#-------------------------------------------------------------------------------------

ET.Penman <- function(data, constants, ts="daily", solar="sunshine hours", wind="yes", windfunction_ver=1948, alpha = 0.08, z0 = 0.001,
                      message = "yes", AdditionalStats= "yes", save.csv = "no", ...) {
  #class(data) <- "funname"

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }

  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }

  if (wind == "yes") { # wind data is required
    if (is.null(data$uz) & is.null(data$u2)) {
      stop("Required data missing for 'uz' or 'u2'")
    }
  }

  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  if (wind != "yes" & wind != "no") {
    stop("Please choose if actual data will be used for wind speed from wind = 'yes' and wind = 'no'")
  }

  # check user-input albedo
  if (wind == "yes") {
    if (is.na(as.numeric(alpha))) {
      stop("Please use a numeric value for the alpha (albedo of evaporative surface)")
    }
    if (!is.na(as.numeric(alpha))) {
      if (as.numeric(alpha) < 0 | as.numeric(alpha) > 1) {
        stop("Please use a value between 0 and 1 for the alpha (albedo of evaporative surface)")
      }
    }
    if (is.na(as.numeric(z0))) {
      stop("Please use a numeric value for the z0 (roughness height)")
    }
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  # Calculations from data and constants for Penman

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar declination (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }

  if (wind == "yes") {
    # Wind speed at 2 meters
    if (is.null(data$u2)) {
      u2 <- data$uz * log(2/z0) / log(constants$z/z0) # Equation S4.4
    } else {
      u2 <- data$u2
    }

    if (!is.null(data$va) & !is.null(data$vs)) {
      vabar <- data$va
      vas <- data$vs
    } else {
      # Saturated vapour pressure
      vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
      vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
      vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

      # Vapour pressure
      vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7

    }

    R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
    R_ns <- (1 - alpha) * R_s # net incoming shortwave radiation - water or other evaporative surface with specified Albedo (S3.2)

    R_n = R_ns - R_nl # net radiation (S3.1)
    if (windfunction_ver == "1948") {
      f_u = 2.626 + 1.381 * u2 # wind function Penman 1948 (S4.11)
    } else if (windfunction_ver == "1956") {
      f_u = 1.313 + 1.381 * u2 # wind function Penman 1956 (S4.3)
    } else if (windfunction_ver != "1948" & windfunction_ver != "1956") {
      stop("Please select the version of wind function (1948 or 1956)")
    }

    Ea = f_u * (vas - vabar) # (S4.2)

    Epenman.Daily <-  delta / (delta +  gamma) * (R_n / constants$lambda) + gamma  / (delta + gamma) * Ea # Penman open-water evaporation (S4.1)
  } else {
    # mean relative humidity
    RHmean <- (data$RHmax + data$RHmin) / 2

    Epenman.Daily <-  0.047 * R_s * sqrt(Ta + 9.5) - 2.4 * (R_s/R_a)^2 + 0.09 * (Ta + 20) * (1 - RHmean/100) # Penman open-water evaporation without wind data by Valiantzas (2006) (S4.12)
  }

  ET.Daily <- Epenman.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }
  # Generate summary message for results
  ET_formulation <- "Penman"
  if (wind == "no") {
    ET_type <- "Open-water Evaporation"
    Surface <- paste("water, albedo =", alpha, "; roughness height =", z0, "m")
  } else {
    if (alpha != 0.08) {
      ET_type <- "Potential ET"
      Surface <- paste("user-defined, albedo =", alpha, "; roughness height =", z0, "m")
    } else if (alpha == 0.08) {
      ET_type <- "Open-water Evaporation"
      Surface <- paste("water, albedo =", alpha, "; roughness height =", z0, "m")
    }
  }
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }

  if (wind == "yes") {
    if (windfunction_ver == "1948") {
      message2 <- "Wind data have been used for calculating the Penman evaporation. Penman 1948 wind function has been used."
    } else if (windfunction_ver == "1956") {
      message2 <- "Wind data have been used for calculating the Penman evaporation. Penman 1956 wind function has been used."
    }
  } else {
    message2 <- "Alternative calculation for Penman evaporation without wind data have been performed"
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1, message2=message2)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {


    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    message(message1)
    message(message2)


    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      #message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_Penman.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_Penman.csv", col.names=F, append= T, sep=',' )
    }
    #class(results) <- funname
    invisible(results)
  } else {
    return(results)
  }
}

#-------------------------------------------------------------------------------------

ET.PenmanMonteith <- function(data, constants, ts="daily", solar="sunshine hours", wind="yes", crop="short",
                              message = "yes", AdditionalStats= "yes", save.csv = "no", ...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }
  if (wind == "yes") { # wind data is required
    if (is.null(data$u2) & is.null(data$uz)) {
      stop("Required data missing for 'uz' or 'u2'")
    }
  }

  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  if (wind != "yes" & wind != "no") {
    stop("Please choose if actual data will be used for wind speed from wind = 'yes' and wind = 'no'")
  }
  # check user-input crop type and specify albedo
  if (wind == "yes") {
    if (crop != "short" & crop != "tall") {
      stop("Please enter 'short' or 'tall' for the desired reference crop type")
    } else {
      alpha <- 0.23 # albedo for both short and tall crop
      if (crop == "short") {
        z0 <- 0.02 # roughness height for short grass
      } else {
        z0 <- 0.1 # roughness height for tall grass
      }
    }
  } else {
    z0 <- 0.02 # roughness height for short grass
    alpha <- 0.25 # semi-desert short grass - will not be used for calculation - just informative
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  if (!is.null(data$va) & !is.null(data$vs)) {
    vabar <- data$va
    vas <- data$vs
  } else {
    # Saturated vapour pressure
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
    vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

    # Vapour pressure
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7

  }
  # Calculations from data and constants for Penman-Monteith Reference Crop

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  } else if (solar!="monthly precipitation") {
    # calculate R_s from sunshine hours - data or estimation using cloudness
    R_s <- (constants$as + constants$bs * (data$n/N))*R_a # estimated incoming solar radiation (S3.9)
  } else {
    # calculate R_s from cloudness estimated from monthly precipitation (#S3.14)
    R_s <- (0.85 - 0.047*data$Cd)*R_a 
  }

  R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
  R_nsg <- (1 - alpha) * R_s # net incoming shortwave radiation (S3.2)
  R_ng <- R_nsg - R_nl # net radiation (S3.1)

  if (wind == "yes") {
    # Wind speed
    if (is.null(data$u2)) {
      u2 <- data$uz * 4.87 / log(67.8*constants$z - 5.42) # Equation S5.20 for PET formulations other than Penman
    } else {
      u2 <- data$u2
    }

    if (crop == "short") {
      r_s <- 70 # will not be used for calculation - just informative
      CH <- 0.12 # will not be used for calculation - just informative
      ET_RC.Daily <- (0.408 * delta * (R_ng - constants$G) + gamma * 900 * u2 * (vas - vabar)/(Ta + 273)) / (delta + gamma * (1 + 0.34*u2)) # FAO-56 reference crop evapotranspiration from short grass (S5.18)
    } else {
      r_s <- 45 # will not be used for calculation - just informative
      CH <- 0.50 # will not be used for calculation - just informative
      ET_RC.Daily <- (0.408 * delta * (R_ng - constants$G) + gamma * 1600 * u2 * (vas - vabar)/(Ta + 273)) / (delta + gamma * (1 + 0.38*u2)) # ASCE-EWRI standardised Penman-Monteith for long grass (S5.19)
    }
    ET.Daily <- ET_RC.Daily
    ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
    ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)
  } else {
    # mean relative humidity
    RHmean <- (data$RHmax + data$RHmin) / 2

    R_s.Monthly <- aggregate(R_s, as.yearmon(data$Date.daily, "%m/%y"),mean)
    R_a.Monthly <- aggregate(R_a, as.yearmon(data$Date.daily, "%m/%y"),mean)
    Ta.Monthly <- aggregate(Ta, as.yearmon(data$Date.daily, "%m/%y"),mean)
    RHmean.Monthly <- aggregate(RHmean, as.yearmon(data$Date.daily, "%m/%y"),mean)
    #ET_RC.Daily <- matrix(NA,length(data$date.Daily),1)
    ET_RC.Monthly <- 0.038 * R_s.Monthly * sqrt(Ta.Monthly + 9.5) - 2.4 * (R_s.Monthly/R_a.Monthly)^2 + 0.075 * (Ta.Monthly + 20) * (1 - RHmean.Monthly/100) # Reference crop evapotranspiration without wind data by Valiantzas (2006) (S5.21)
    ET_RC.Daily <- data$Tmax
    for (cont in 1:length(data$Date.monthly)) {
      ET_RC.Daily[which(as.yearmon(time(ET_RC.Daily))==data$Date.monthly[cont])] <- ET_RC.Monthly[cont]
    }

    ET.Daily <- ET_RC.Daily
    ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)	
    ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily))), FUN = sum)	
  }

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  if (wind == "no") {
    ET_formulation <- "Penman-Monteith (without wind data)"
    ET_type <- "Reference Crop ET"
    Surface <- paste("short grass, albedo =", alpha, "; roughness height =", z0, "m")
  } else {
    if (crop == "short") {
      ET_formulation <- "Penman-Monteith FAO56"
      ET_type <- "Reference Crop ET"
      Surface <- paste("FAO-56 hypothetical short grass, albedo =", alpha, "; surface resistance =", r_s, "sm^-1; crop height =", CH, " m; roughness height =", z0, "m")
    } else {
      ET_formulation <- "Penman-Monteith ASCE-EWRI Standardised"
      ET_type <- "Reference Crop ET"
      Surface <- paste("ASCE-EWRI hypothetical tall grass, albedo =", alpha, "; surface resistance =", r_s, "sm^-1; crop height =", CH, " m; roughness height =", z0, "m")
    }
  }
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }

  if (wind == "yes") {
    message2 <- "Wind data have been used for calculating the reference crop evapotranspiration"
  } else {
    message2 <- "Alternative calculation for reference crop evapotranspiration without wind data have been performed"
  }

  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1, message2=message2)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }
  if (message == "yes") {

    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    message(message1)
    message(message2)


    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
  }

  #class(results) <- funname
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_PenmanMonteith.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_PenmanMonteith.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }


}

#-------------------------------------------------------------------------------------

ET.MattShuttleworth <- function(data, constants, ts="daily", solar="sunshine hours", alpha=0.23, r_s=70, CH=0.12,
                                message = "yes", AdditionalStats= "yes", save.csv = "no", ...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs.daily', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }
  if (is.null(data$u2) & is.null(data$uz)) {
    stop("Required data missing for 'uz' or 'u2'")
  }

  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  # check user-input albedo, surface resistance and crop height
  if (is.na(as.numeric(alpha))) {
    stop("Please use a numeric value for the alpha (albedo of evaporative surface)")
  }
  if (is.na(as.numeric(r_s))) {
    stop("Please use a numeric value for the r_s (surface resistance) in sm^-1")
  }
  if (is.na(as.numeric(CH))) {
    stop("Please use a numeric value for the CH (crop height) in m")
  }
  if (!is.na(as.numeric(alpha))) {
    if (as.numeric(alpha) < 0 | as.numeric(alpha) > 1) {
      stop("Please use a value between 0 and 1 for the alpha (albedo of evaporative surface)")
    }
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  if (!is.null(data$va) & !is.null(data$vs)) {
    vabar <- data$va
    vas <- data$vs
  } else {
    # Saturated vapour pressure
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
    vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

    # Vapour pressure
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7
  }

  # Calculations from data and constants for Matt-Shuttleworth Reference Crop

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }

  R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
  # For short grass
  R_nsg <- (1 - alpha) * R_s # net incoming shortwave radiation (S3.2)
  R_ng <- R_nsg - R_nl # net radiation (S3.1)

  # Wind speed
  if (is.null(data$u2)) {
    u2 <- data$uz * 4.87 / log(67.8*constants$z - 5.42) # Equation S5.20 for PET formulations other than Penman
  } else {
    u2 <- data$u2
  }

  r_clim <- 86400 * constants$Roua * constants$Ca * (vas - vabar) / (delta * R_ng) # clinmatological resistance (s*m^-1) (S5.34)
  r_clim[r_clim == 0] <- 0.1 # correction for r_clim = 0
  u2[u2 == 0] <- 0.1 # correction for u2 = 0
  VPD50toVPD2 <- (302 * (delta + gamma) + 70 * gamma * u2) / (208 * (delta + gamma) + 70 * gamma * u2) + 1/r_clim * ((302 * (delta + gamma) + 70 * gamma * u2) / (208 * (delta + gamma) + 70 * gamma * u2) * (208 / u2) - (302 / u2)) # ratio of vapour pressure deficits at 50m to vapour pressure deficits at 2m heights (S5.35)
  r_c50 <- 1 / ((0.41)^2) * log((50 - 0.67 * CH) / (0.123 * CH)) * log((50 - 0.67 * CH) / (0.0123 * CH)) * log((2 - 0.08) / 0.0148) / log((50 - 0.08) / 0.0148) # aerodynamic coefficient (s*m^-1) (S5.36)

  E_Tc.Daily <- 1/constants$lambda * (delta * R_ng + (constants$Roua * constants$Ca * u2 * (vas - vabar)) / r_c50 * VPD50toVPD2) / (delta + gamma * (1 + r_s * u2 / r_c50)) # well-watered crop evapotranspiration in a semi-arid and windy location (S5.37)

  ET.Daily <- E_Tc.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }
  # Generate summary message for results
  ET_formulation <- "Matt-Shuttleworth"
  ET_type <- "Reference Crop ET"
  Surface <- paste("user-defined, albedo =", alpha, "; surface resistance =", r_s, "sm^-1; crop height =", CH, "m")
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {


    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    message(message1)


    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
  }

  if (save.csv == "yes") {
    #class(results) <- funname
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_MattShuttleworth.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_MattShuttleworth.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
    # write to csv file
  } else {
    return(results)
  }


}

#-------------------------------------------------------------------------------------

ET.PriestleyTaylor <- function(data, constants, ts="daily", solar="sunshine hours", alpha=0.23,
                               message = "yes", AdditionalStats= "yes", save.csv = "no", ...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }

  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  # check user-input albedo
  if (is.na(as.numeric(alpha))) {
    stop("Please use a numeric value for the alpha (albedo of evaporative surface)")
  }
  if (!is.na(as.numeric(alpha))) {
    if (as.numeric(alpha) < 0 | as.numeric(alpha) > 1) {
      stop("Please use a value between 0 and 1 for the alpha (albedo of evaporative surface)")
    }
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  if (!is.null(data$va) & !is.null(data$vs)) {
    vabar <- data$va
    vas <- data$vs
  } else {
    # Saturated vapour pressure
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
    vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

    # Vapour pressure
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7
  }

  # Calculations from data and constants for Matt-Shuttleworth Reference Crop

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }

  R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
  # For short grass
  R_nsg <- (1 - alpha) * R_s # net incoming shortwave radiation (S3.2)
  R_ng <- R_nsg - R_nl # net radiation (S3.1)

  E_PT.Daily <- constants$alphaPT * (delta/(delta + gamma) * R_ng / constants$lambda - constants$G / constants$lambda) # well-watered crop evapotranspiration in a semi-arid and windy location (S5.37)

  ET.Daily <- E_PT.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "Priestley-Taylor"
  ET_type <- "Potential ET"
  if (alpha != 0.08) {
    Surface <- paste("user-defined, albedo =", alpha)
  } else if (alpha == 0.08) {
    Surface <- paste("water, albedo =", alpha)
  }
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {


    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    message(message1)


    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
  }

  if (save.csv == "yes") {
    #class(results) <- funname
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_PriestleyTaylor.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_PriestleyTaylor.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.PenPan <- function (data, constants, ts="daily", solar="sunshine hours", alpha=0.23, est="potential ET", pan_coeff=0.71, overest=F, message="yes",AdditionalStats= "yes",  save.csv="no", ...) {
  #class(data) <- funname
  if (is.null(data$Tmax) | is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }
  if (is.null(data$u2) & is.null(data$uz)) {
    stop("Required data missing for 'uz' or 'u2'")
  }
  if (solar == "data" & is.null(data$Rs)) {
    stop("Required data missing for 'Rs'")
  }
  else if (solar == "sunshine hours" & is.null(data$n)) {
    stop("Required data missing for 'n'")
  }
  else if (solar == "cloud" & is.null(data$Cd)) {
    stop("Required data missing for 'Cd'")
  }
  else if (solar == "monthly precipitation" & is.null(data$Precip)) {
    stop("Required data missing for 'Precip'")
  }
  if (is.na(as.numeric(alpha))) {
    stop("Please use a numeric value for the alpha (albedo of evaporative surface)")
  }
  if (!is.na(as.numeric(alpha))) {
    if (as.numeric(alpha) < 0 | as.numeric(alpha) > 1) {
      stop("Please use a value between 0 and 1 for the alpha (albedo of evaporative surface)")
    }
  }
  Ta <- (data$Tmax + data$Tmin)/2

  if (!is.null(data$va) & !is.null(data$vs)) {
    vabar <- data$va
    vas <- data$vs
  } else {
    # Saturated vapour pressure
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
    vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

    # Vapour pressure
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7
  }

  P <- 101.3 * ((293 - 0.0065 * constants$Elev)/293)^5.26
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta + 237.3)))/((Ta +
                                                                237.3)^2)
  gamma <- 0.00163 * P/constants$lambda
  d_r2 <- 1 + 0.033 * cos(2 * pi/365 * data$J)
  delta2 <- 0.409 * sin(2 * pi/365 * data$J - 1.39)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))
  N <- 24/pi * w_s
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) *
                                               sin(delta2) + cos(constants$lat_rad) * cos(delta2) *
                                               sin(w_s))
  R_so <- (0.75 + (2 * 10^-5) * constants$Elev) * R_a

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }
  if (is.null(data$u2)) {
    u2 <- data$uz * 4.87/log(67.8 * constants$z - 5.42)
  }
  else {
    u2 <- data$u2
  }
  R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) *
    ((data$Tmax + 273.2)^4 + (data$Tmin + 273.2)^4)/2 *
    (1.35 * R_s/R_so - 0.35)
  P_rad <- 1.32 + 4 * 10^(-4) * abs(constants$lat) + 8 * 10^(-5) *
    (constants$lat)^2
  f_dir <- -0.11 + 1.31 * R_s/R_a
  R_span <- (f_dir * P_rad + 1.42 * (1 - f_dir) + 0.42 * alpha) *
    R_s
  R_npan <- (1 - constants$alphaA) * R_span - R_nl
  f_pan_u <- 1.201 + 1.621 * u2
  Epenpan.Daily <- delta/(delta + constants$ap * gamma) *
    R_npan/constants$lambda + constants$ap * gamma/(delta +
                                                      constants$ap * gamma) * f_pan_u * (vas - vabar)
  if (overest == TRUE) {
    if (est == "potential ET") {
      Epenpan.Daily <- Epenpan.Daily/1.078*pan_coeff
    } else {
      Epenpan.Daily <- Epenpan.Daily/1.078
    }
  }

  ET.Daily <- Epenpan.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily,
                                               "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily,
                                                               "%m/%y"))), FUN = sum)
  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }
  ET_formulation <- "PenPan"
  if (est == "potential ET") {
    ET_type <- "potential ET"
  } else if (est == "pan") {
    ET_type <- "Class-A Pan Evaporation"
  }
  Surface <- paste("user-defined, albedo =", alpha)
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  }
  else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  }
  else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  }
  else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }
  results <- list(ET.Daily = ET.Daily, ET.Monthly = ET.Monthly,
                  ET.Annual = ET.Annual, ET.MonthlyAve = ET.MonthlyAve,
                  ET.AnnualAve = ET.AnnualAve, ET_formulation = ET_formulation,
                  ET_type = ET_type, message1 = message1)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {

    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    if (ET_type == "potential ET") {
      message("Pan coeffcient: ", pan_coeff)
    }
    message(message1)


    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
  }

  if (save.csv == "yes") {
    #class(results) <- funname
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_PenPan.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_PenPan.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}


#-------------------------------------------------------------------------------------

ET.BrutsaertStrickler <- function(data, constants, ts="daily", solar="sunshine hours", alpha=0.23,
                                  message = "yes", AdditionalStats= "yes", save.csv = "no", ...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }
  if (is.null(data$u2) & is.null(data$uz)) {
    stop("Required data missing for 'uz' or 'u2'")
  }

  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  # check user-input albedo
  if (is.na(as.numeric(alpha))) {
    stop("Please use a numeric value for the alpha (albedo of evaporative surface)")
  }
  if (!is.na(as.numeric(alpha))) {
    if (as.numeric(alpha) < 0 | as.numeric(alpha) > 1) {
      stop("Please use a value between 0 and 1 for the alpha (albedo of evaporative surface)")
    }
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  if (!is.null(data$va) & !is.null(data$vs)) {
    vabar <- data$va
    vas <- data$vs
  } else {
    # Saturated vapour pressure
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
    vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

    # Vapour pressure
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7
  }

  # update alphaPT according to Brutsaert and Strickler (1979)
  constants$alphaPT <- 1.28

  # Calculations from data and constants for Brutsaert and Strickler actual evapotranspiration

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }

  # Wind speed
  if (is.null(data$u2)) {
    u2 <- data$uz * 4.87 / log(67.8*constants$z - 5.42) # Equation S5.20 for PET formulations other than Penman
  } else {
    u2 <- data$u2
  }

  R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
  # For short grass
  R_nsg <- (1 - alpha) * R_s # net incoming shortwave radiation (S3.2)
  R_ng <- R_nsg - R_nl # net radiation (S3.1)
  f_u2 <- 2.626 + 1.381 * u2 # Penman's wind function adopted with Priestley and Tayloy constant alphaPT by Brutsaert and Strickler (1979) (S8.3)

  ET_BS_Act.Daily <- (2 * constants$alphaPT - 1) * (delta / (delta + gamma)) * R_ng / constants$lambda - gamma / (delta + gamma) * f_u2 * (vas - vabar) # Brutsaert and Strickler actual areal evapotranspiration (mm.day^-1) Brutsaert and Strickler (1979) (S8.2)

  ET.Daily <- ET_BS_Act.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "Brutsaert-Strickler"
  ET_type <- "Actual Areal ET"
  Surface <- paste("user-defined, albedo =", alpha)
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {


    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    message(message1)


    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }

  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_BrutsaertStrickler.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_BrutsaertStrickler.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.GrangerGray <- function(data, constants, ts="daily", solar="sunshine hours", windfunction_ver=1948, alpha=0.23,
                           message = "yes", AdditionalStats= "yes", save.csv = "no",...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }
  if (is.null(data$u2) & is.null(data$uz)) {
    stop("Required data missing for 'uz' or 'u2'")
  }

  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  # check user-input albedo
  if (is.na(as.numeric(alpha))) {
    stop("Please use a numeric value for the alpha (albedo of evaporative surface)")
  }
  if (!is.na(as.numeric(alpha))) {
    if (as.numeric(alpha) < 0 | as.numeric(alpha) > 1) {
      stop("Please use a value between 0 and 1 for the alpha (albedo of evaporative surface)")
    }
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  if (!is.null(data$va) & !is.null(data$vs)) {
    vabar <- data$va
    vas <- data$vs
  } else {
    # Saturated vapour pressure
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
    vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

    # Vapour pressure
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7
  }
  # Calculations from data and constants for Granger and Gray actual evapotranspiration

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }
  R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
  # For short grass
  R_nsg <- (1 - alpha) * R_s # net incoming shortwave radiation (S3.2)
  R_ng <- R_nsg - R_nl # net radiation (S3.1)

  # Wind speed
  if (is.null(data$u2)) {
    u2 <- data$uz * 4.87 / log(67.8*constants$z - 5.42) # Equation S5.20 for PET formulations other than Penman
  } else {
    u2 <- data$u2
  }

  if (windfunction_ver == "1948") {
    f_u = 2.626 + 1.381 * u2 # wind function Penman 1948 (S4.11)
  } else if (windfunction_ver == "1956") {
    f_u = 1.313 + 1.381 * u2 # wind function Penman 1956 (S4.3)
  } else if (windfunction_ver != "1948" & windfunction_ver != "1956") {
    stop("Please select the version of wind function (1948 or 1956)")
  }
  Ea = f_u * (vas - vabar) # (S4.2)
  D_p <- Ea / (Ea + (R_ng - constants$G) / constants$lambda) # dimensionless relative drying power (S8.6)
  G_g <- 1 / (0.793 + 0.20 * exp(4.902 * D_p)) + 0.006 * D_p # dimensionless evaporation parameter (S8.5)
  ET_GG_Act.Daily <- delta * G_g / (delta * G_g + gamma) * (R_ng - constants$G) / constants$lambda + gamma * G_g / (delta * G_g + gamma) * Ea # Granger and Gray actual areal evapotranspiration (mm.day^-1) Granger and Gray (1989) (S8.4)

  ET.Daily <- ET_GG_Act.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "Granger-Gray"
  ET_type <- "Actual Areal ET"
  Surface <- paste("user-defined, albedo =", alpha)
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }

  if (windfunction_ver == "1948") {
    message2 <- "Wind data have been used for the calculation of the drying power of air, using Penman 1948 wind function."
  } else if (windfunction_ver == "1956") {
    message2 <- "Wind data have been used for the calculation of the drying power of air, using Penman 1956 wind function."
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1, message2=message2)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }
  if (message == "yes") {


    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    message(message1)
    message(message2)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }

  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_GrangerGray.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_GrangerGray.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }


}

#-------------------------------------------------------------------------------------

ET.SzilagyiJozsa <- function(data, constants, ts="daily", solar="sunshine hours", wind="yes", windfunction_ver=1948, alpha=0.23, z0=0.2,
                             message = "yes", AdditionalStats= "yes", save.csv = "no", ...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }
  if (wind == "yes") { # wind data is required
    if (is.null(data$u2) & is.null(data$uz)) {
      stop("Required data missing for 'uz' or 'u2'")
    }
  }

  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  if (wind != "yes" & wind != "no") {
    stop("Please choose if actual data will be used for wind speed from wind = 'yes' and wind = 'no'")
  }

  # check user-input albedo
  if (wind == "yes") {
    if (is.na(as.numeric(alpha))) {
      stop("Please use a numeric value for the alpha (albedo of evaporative surface)")
    }
    if (!is.na(as.numeric(alpha))) {
      if (as.numeric(alpha) < 0 | as.numeric(alpha) > 1) {
        stop("Please use a value between 0 and 1 for the alpha (albedo of evaporative surface)")
      }
    }
    if (is.na(as.numeric(z0))) {
      stop("Please use a numeric value for the z0 (roughness height)")
    }
  }

  # update alphaPT according to Szilagyi and Jozsa (2008)
  constants$alphaPT <- 1.31

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  if (!is.null(data$va) & !is.null(data$vs)) {
    vabar <- data$va
    vas <- data$vs
  } else {
    # Saturated vapour pressure
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
    vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

    # Vapour pressure
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7
  }

  # Calculations from data and constants for Penman

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }
  if (wind == "yes") {
    # Wind speed at 2 meters
    if (is.null(data$u2)) {
      u2 <- data$uz * log(2/z0) / log(constants$z/z0) # Equation S4.4
    } else {
      u2 <- data$u2
    }


    R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
    # For vegetated surface
    R_nsg <- (1 - alpha) * R_s # net incoming shortwave radiation (S3.2)
    R_ng = R_nsg - R_nl # net radiation (S3.1)
    if (windfunction_ver == "1948") {
      f_u = 2.626 + 1.381 * u2 # wind function Penman 1948 (S4.11)
    } else if (windfunction_ver == "1956") {
      f_u = 1.313 + 1.381 * u2 # wind function Penman 1956 (S4.3)
    } else if (windfunction_ver != "1948" & windfunction_ver != "1956") {
      stop("Please select the version of wind function (1948 or 1956)")
    }
    Ea = f_u * (vas - vabar) # (S4.2)

    Epenman.Daily <-  delta / (delta +  gamma) * (R_ng / constants$lambda) + gamma  / (delta + gamma) * Ea # Penman open-water evaporation (S4.1)
  } else {
    R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
    # For vegetated surface
    R_nsg <- (1 - alpha) * R_s # net incoming shortwave radiation (S3.2)
    R_ng = R_nsg - R_nl # net radiation (S3.1)
    # mean relative humidity
    RHmean <- (data$RHmax + data$RHmin) / 2

    Epenman.Daily <-  0.047 * R_s * sqrt(Ta + 9.5) - 2.4 * (R_s/R_a)^2 + 0.09 * (Ta + 20) * (1 - RHmean/100) # Penman open-water evaporation without wind data by Valiantzas (2006) (S4.12)
  }

  # Iteration for equilibrium temperature T_e
  T_e <- Ta
  for (i in 1:99999) {
    v_e <- 0.6108 * exp(17.27 * T_e/(T_e + 237.3)) # saturated vapour pressure at T_e (S2.5)
    T_enew <- Ta - 1 / gamma * (1 - R_ng / (constants$lambda * Epenman.Daily)) * (v_e - vabar) # rearranged from S8.8
    deltaT_e <- na.omit(T_enew - T_e)
    maxdeltaT_e <- abs(max(deltaT_e))
    T_e <- T_enew
    if (maxdeltaT_e < 0.01) break
  }
  deltaT_e <- 4098 * (0.6108 * exp((17.27 * T_e)/(T_e+237.3))) / ((T_e + 237.3)^2)  # slope of vapour pressure curve (S2.4)
  E_PT_T_e <- constants$alphaPT * (deltaT_e / (deltaT_e + gamma) * R_ng / constants$lambda) # Priestley-Taylor evapotranspiration at T_e
  E_SJ_Act.Daily <- 2 * E_PT_T_e - Epenman.Daily # actual evapotranspiration by Szilagyi and Jozsa (2008) (S8.7)

  ET.Daily <- E_SJ_Act.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "Szilagyi-Jozsa"
  ET_type <- "Actual ET"
  Surface <- paste("user-defined, albedo =", alpha, "; roughness height", z0, "m")

  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }

  if (wind == "yes") {
    if (windfunction_ver == "1948") {
      message2 <- "Wind data have been used for calculating the Penman evaporation. Penman 1948 wind function has been used."
    } else if (windfunction_ver == "1956") {
      message2 <- "Wind data have been used for calculating the Penman evaporation. Penman 1956 wind function has been used."
    }
  } else {
    message2 <- "Alternative calculation for Penman evaporation without wind data have been performed"
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1, message2=message2)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }


  if (message == "yes") {


    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    message(message1)
    message(message2)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }

  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_SzilagyiJozsa.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_SzilagyiJozsa.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.Makkink <- function(data, constants, ts="daily", solar="sunshine hours",
                       message = "yes", AdditionalStats= "yes", save.csv = "no", ...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  # Calculations from data and constants for Makkink

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }

  Emakkink.Daily <- 0.61 * (delta / (delta + gamma) * R_s/2.45) - 0.12  # potential evapotranspiration by Bruin (1981) (S9.6)

  ET.Daily <- Emakkink.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "Makkink"
  ET_type <- "Reference crop ET"
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }

  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {

    message(ET_formulation, " ", ET_type)
    message(message1)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }

  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_Makkink.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_Makkink.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.BlaneyCriddle <- function(data, constants, ts="daily", solar="sunshine hours", height=F,
                             message = "yes", AdditionalStats= "yes", save.csv = "no",...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (solar == "sunshine hours" & is.null(data$n)) { # sunshine hour data is required
    stop("Required data missing for 'n'")
  } else if (solar == "cloud") {
    if (is.null(data$n)) { # for alternative calculation of sunshine hours using cloud cover
      stop("Required data missing for 'Cd'")
    }
    if (is.null(data$u2) & is.null(data$uz)) {
      stop("Required data missing for 'uz' or 'u2'")
    }
    if (is.null(data$RHmin)) {
      stop("Required data missing for 'RHmin'")
    }
  }
  if (solar == "data" | solar == "monthly precipitation") {
    stop("Only 'sunshine hours' and 'cloud' are accepted because estimations of sunshine hours is required")
  }
  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  # Calculations from data and constants for Blaney and Criddle
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values

  # Wind speed
  if (is.null(data$u2)) {
    u2 <- data$uz * 4.87 / log(67.8*constants$z - 5.42) # Equation S5.20 for PET formulations other than Penman
  } else {
    u2 <- data$u2
  }

  bvar <- constants$e0 + constants$e1 * data$RHmin + constants$e2 * data$n/N + constants$e3 * u2 + constants$e4 * data$RHmin * data$n/N + constants$e5 * data$RHmin * u2 # undefined working variable (Allena and Pruitt, 1986; Shuttleworth, 1992) (S9.8)
  N.annual <- ave(N, format(time(N), "%y"), FUN = sum) # Annual sum of maximum sunshine hours
  # chech if data from first/last years is incomplete, and adjust N.annual values for incomplete years
  if(data$J[1]!=1 & is.integer(floor(as.numeric(as.yearmon(data$Date.daily)))[1]/4)==FALSE) { # first year a normal year
    N.annual[floor(as.numeric(as.yearmon(data$Date.daily)))==floor(as.numeric(as.yearmon(data$Date.daily)))[1]] <-
      sum(24/pi * acos(-tan(constants$lat_rad) * tan(0.409 * sin(2*pi/365 * c(1:365) - 1.39))))
  }
  if(data$J[1]!=1 & is.integer(floor(as.numeric(as.yearmon(data$Date.daily)))[1]/4)==TRUE) { # first year a leap year
    N.annual[floor(as.numeric(as.yearmon(data$Date.daily)))==floor(as.numeric(as.yearmon(data$Date.daily)))[1]] <-
      sum(24/pi * acos(-tan(constants$lat_rad) * tan(0.409 * sin(2*pi/365 * c(1:366) - 1.39))))
  }
  if (data$J[length(data$J)]!=365 & is.integer(floor(as.numeric(as.yearmon(data$Date.daily)))[length(data$J)]/4)==FALSE) { # last year a normal year
    N.annual[floor(as.numeric(as.yearmon(data$Date.daily)))==floor(as.numeric(as.yearmon(data$Date.daily)))[length(data$J)]] <-
      sum(24/pi * acos(-tan(constants$lat_rad) * tan(0.409 * sin(2*pi/365 * c(1:365) - 1.39))))
  }
  if (data$J[length(data$J)]!=366 & is.integer(floor(as.numeric(as.yearmon(data$Date.daily)))[length(data$J)]/4)==TRUE) { # first year a leap year
    N.annual[floor(as.numeric(as.yearmon(data$Date.daily)))==floor(as.numeric(as.yearmon(data$Date.daily)))[length(data$J)]] <-
      sum(24/pi * acos(-tan(constants$lat_rad) * tan(0.409 * sin(2*pi/366 * c(1:366) - 1.39))))
  }
  p_y <- 100 * data$n/N.annual # percentage of actual daytime hours for the day comparing to the annual sum of maximum sunshine hours

  ET_BC.Daily <- (0.0043 * data$RHmin - data$n/N - 1.41) + bvar * p_y * (0.46 * Ta +8.13) # Blaney-Criddle Reference Crop evapotranspiration (mm.day^-1) (S9.7)

  ET.Daily <- ET_BC.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  if (height == T) {
    ET_BC.Daily = ET_BC.Daily * (1 + 0.1 * constants$Elev/1000) # with adjustment for site elevation by Allen and Pruitt (1986) (S9.9)
  }

  # Generate summary message for results
  ET_formulation <- "Blaney-Criddle"
  ET_type <- "Reference Crop ET"
  if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }

  if (height == T) {
    message3 <- "Height adjustment has been applied to calculated Blaney-Criddle reference crop evapotranspiration"
  } else {
    message3 <- "No height adjustment has been applied to calculated Blaney-Criddle reference crop evapotranspiration"
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1, message3=message3)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }
  if (message == "yes") {


    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: reference crop")
    message(message1)
    message(message3)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_BlaneyCriddle.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_BlaneyCriddle.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.Turc <- function(data, constants, ts="daily", solar="sunshine hours", humid=F,
                    message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }

  if (humid == TRUE & (is.null(data$RHmax)|is.null(data$RHmin))) { # for adjustment for non-humid conditions
    stop("Required data missing for 'RHmax' and 'RHmin', or 'RH'")
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  # Calculations from data and constants for Turc

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }
  ET_Turc.Daily <- 0.013 * (23.88 * R_s + 50) * Ta / (Ta + 15) # reference crop evapotranspiration by Turc (1961) (S9.10)

  if (humid == TRUE) {
    # mean relative humidity
    RHmean <- (data$RHmax + data$RHmin) / 2

    ET_Turc.Daily[RHmean < 50] <- 0.013 * (23.88 * R_s + 50) * Ta[RHmean < 50] / (Ta[RHmean < 50] + 15) * (1 + (50 - RHmean[RHmean < 50]) / 70) # Turc reference crop evapotranspiration adjusted for non-humid conditions (RH < 50) by Alexandris et al., (S9.11)
  }

  ET.Daily <- ET_Turc.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "Turc"
  ET_type <- "Reference Crop ET"

  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }

  if (humid == TRUE) {
    message4 <- "Adjustment for non-humid conditions has been applied to calculated Turc reference crop evapotranspiration"
  } else {
    message4 <- "No adjustment for non-humid conditions has been applied to calculated Turc reference crop evapotranspiration"
  }


  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1, message4=message4)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {
    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: reference crop")
    message(message1)
    message(message4)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv =="yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_Turc.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_Turc.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.Hamon <- function(data, constants = NULL, ts="daily",
                     message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (is.null(data$n)) {
    stop("Required data missing for 'n'")
  }
  Ta <- (data$Tmax + data$Tmin) / 2

  # Saturated vapour pressure
  vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
  vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
  vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

  ET_Hamon.Daily <- 0.55 * 25.4 * (data$n/12)^2 * (216.7 * vas * 10 / (Ta + 273.3))/100 # Rosenberry et al., 2004

  ET.Daily <- ET_Hamon.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "Hamon"
  ET_type <- "Potential ET"

  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type)

  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {
    message(ET_formulation, " ", ET_type)
    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_Hamon.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_Hamon.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}
# Oudin et al., 2005

#-------------------------------------------------------------------------------------

ET.Linacre <- function(data, constants, ts="daily",
                       message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  # Check of specific data requirement
  if (is.null(data$Tdew)) {
    stop("Required data missing for 'Tdew' or 'Tdew'")
  }

  Ta <- (data$Tmax + data$Tmin) / 2
  T_m <- Ta + 0.006 * constants$Elev
  ET_Linacre.Daily <- (500 * T_m /(100 - constants$lat)+15*(Ta - data$Tdew))/(80 - Ta)

  ET.Daily <- ET_Linacre.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "Linacre"
  ET_type <- "Actual ET"

  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type)
  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }
  if (message == "yes") {
    message(ET_formulation, " ", ET_type)
    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_Linacre.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_Linacre.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}
# Linacre ET. 1977. A simple formula for estimating evaporation rates in various climates, using temperature data alone. AgriculturalMeteorology 18: 409-424.

#-------------------------------------------------------------------------------------

ET.Romanenko <- function(data, constants = NULL, ts="daily",
                         message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  # Check of specific data requirement
  if (is.null(data$va)|is.null(data$vs)) {
    if (is.null(data$RHmax)|is.null(data$RHmin)) {
      stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
    }
  }
  Ta <- (data$Tmax + data$Tmin) / 2

  if (!is.null(data$va) & !is.null(data$vs)) {
    vabar <- data$va
    vas <- data$vs
  } else {
    # Saturated vapour pressure
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
    vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

    # Vapour pressure
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7
  }

  ET_Romanenko.Daily <- 4.5 * (1 + Ta / 25)^2 * (1 - vabar/vas)
  ET.Daily <- ET_Romanenko.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }
  # Generate summary message for results
  ET_formulation <- "Romanenko"
  ET_type <- "Actual ET"

  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type)
  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {
    message(ET_formulation, " ", ET_type)
    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_Romanenko.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_Romanenko.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}
# Oudin et al., 2005

#-------------------------------------------------------------------------------------

ET.Abtew <- function(data, constants, ts="daily", solar="sunshine hours",
                     message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs.daily'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n.daily'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd.daily'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip.daily'")
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  # Calculations from data and constants for Penman

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
  R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }

  ET_Abtew.Daily <- 0.52 * R_s/constants$lambda

  ET.Daily <- ET_Abtew.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }
  }
  # Generate summary message for results
  ET_formulation <- "Abtew"
  ET_type <- "Actual ET"
  if (solar == "data") {
    message1 <- "Solar radiation data have been used directly for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1)
  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {


    message(ET_formulation, " ", ET_type)
    message(message1)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_Abtew.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_Abtew.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}
# Abtew, W. (1996), EVAPOTRANSPIRATION MEASUREMENTS AND MODELING FOR THREE WETLAND SYSTEMS IN SOUTH FLORIDA. JAWRA Journal of the American Water Resources Association, 32: 465-473. doi: 10.1111/j.1752-1688.1996.tb04044.x
#-------------------------------------------------------------------------------------

ET.HargreavesSamani <- function(data, constants, ts="daily",
                                message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  # Calculations from data and constants for Hargreaves-Samani

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)

  C_HS <- 0.00185 * (data$Tmax - data$Tmin)^2 - 0.0433 * (data$Tmax - data$Tmin) + 0.4023 # empirical coefficient by Hargreaves and Samani (1985) (S9.13)
  ET_HS.Daily <- 0.0135 * C_HS * R_a / constants$lambda * (data$Tmax - data$Tmin)^0.5 * (Ta + 17.8) # reference crop evapotranspiration by Hargreaves and Samani (1985) (S9.12)

  ET.Daily <- ET_HS.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }
  # Generate summary message for results
  ET_formulation <- "Hargreaves-Samani"
  ET_type <- "Reference Crop ET"

  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type)
  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }
  if (message == "yes") {
    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: reference crop")

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_HargreavesSamani.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_HargreavesSamani.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.ChapmanAustralian <- function(data, constants, ts="daily",PenPan=T, solar="sunshine hours", alpha=0.23,
                                 message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (PenPan == TRUE) { # Calculate Class-A pan evaporation using PenPan formula
    if (is.null(data$Tmax)|is.null(data$Tmin)) {
      stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
    }
    if (is.null(data$va)|is.null(data$vs)) {
      if (is.null(data$RHmax)|is.null(data$RHmin)) {
        stop("Required data missing: need either 'va' and 'vs', or 'RHmax' and 'RHmin' (or 'RH')")
      }
    }
    if (is.null(data$u2) & is.null(data$uz)) {
      stop("Required data missing for 'uz' or 'u2'")
    }
    if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
      stop("Required data missing for 'Rs'")
    } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
      stop("Required data missing for 'n'")
    } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
      stop("Required data missing for 'Cd'")
    } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
      stop("Required data missing for 'Precip'")
    }
  }
  if (PenPan == FALSE & is.null(data$Epan)) { # for using Class-A pan evaporation data
    stop("Required data missing for 'Epan'")
  }
  # check user-input albedo
  if (PenPan == TRUE) {
    if (is.na(as.numeric(alpha))) {
      stop("Please use a numeric value for the alpha (albedo of evaporative surface)")
    }
    if (!is.na(as.numeric(alpha))) {
      if (as.numeric(alpha) < 0 | as.numeric(alpha) > 1) {
        stop("Please use a value between 0 and 1 for the alpha (albedo of evaporative surface)")
      }
    }
  }

  # Calculations from data and constants for daily equivalent Penman-Monteith potential evaporation
  A_p <- 0.17 + 0.011 * abs(constants$lat) # constant (S13.2)
  B_p <- 10 ^ (0.66 - 0.211 * abs(constants$lat)) # constants (S13.3)

  if (PenPan == TRUE) {
    # Calculating mean temperature
    Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

    if (!is.null(data$va) & !is.null(data$vs)) {
      vabar <- data$va
      vas <- data$vs
    } else {
      # Saturated vapour pressure
      vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax / (data$Tmax + 237.3)) # Equation S2.5
      vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin / (data$Tmin + 237.3)) # Equation S2.5
      vas <- (vs_Tmax + vs_Tmin)/2 # Equation S2.6

      # Vapour pressure
      vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2 # Equation S2.7
    }
    # estimating class-A pan evaporation using PenPan model
    P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
    delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
    gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
    d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
    delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
    w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
    N <- 24/pi * w_s # calculating daily values
    R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)
    R_so <- (0.75 + (2*10^-5)*constants$Elev) * R_a # clear sky radiation (S3.4)

    if (solar == "data") {
      R_s <- data$Rs
    }else if (solar != "monthly precipitation"&
              solar != "cloud") {
      R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
    }else {
      R_s <- (0.85 - 0.047 * data$Cd) * R_a
    }
    # Wind speed
    if (is.null(data$u2)) {
      u2 <- data$uz * 4.87 / log(67.8*constants$z - 5.42) # Equation S5.20 for PET formulations other than Penman
    } else {
      u2 <- data$u2
    }

    R_nl <- constants$sigma * (0.34 - 0.14 * sqrt(vabar)) * ((data$Tmax+273.2)^4 + (data$Tmin+273.2)^4)/2  * (1.35 * R_s / R_so - 0.35) # estimated net outgoing longwave radiation (S3.3)
    # For short grass
    P_rad <- 1.32 + 4 * 10^(-4) * abs(constants$lat) + 8 * 10^(-5) * (constants$lat)^2 # pan radiation factor (S6.6)
    f_dir <- -0.11 + 1.31 * R_s / R_a # fraction of R_S that os direct (S6.5)
    R_span <- (f_dir * P_rad + 1.42 * (1 - f_dir) + 0.42 * alpha) * R_s # total shortwave radiation received (S6.4)
    R_npan <-(1 - constants$alphaA) * R_span - R_nl # net radiation at the pan (S6.3)
    f_pan_u <-1.201 + 1.621 * u2 # (S6.2)

    Epan <- delta / (delta + constants$ap * gamma) * R_npan / constants$lambda + constants$ap * gamma / (delta + constants$ap * gamma) * f_pan_u * (vas - vabar) # PenPan estimation of Class-A pan evaporation (S6.1)
    ET_eqPM.Daily <- A_p * Epan + B_p # daily equivalent Penman-Monteith potential evaporation (mm.day^-1)
  } else if (PenPan == FALSE & is.null(data$Epan)) {
    stop("No data available for Class-A pan evaporation ")
  } else {
    ET_eqPM.Daily <- A_p * data$Epan + B_p # daily equivalent Penman-Monteith potential evaporation (mm.day^-1)
  }

  ET.Daily <- ET_eqPM.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL

  if (AdditionalStats=="yes") {

    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }


  }
  # Generate summary message for results
  ET_formulation <- "Chapman"
  ET_type <- "Potential ET"
  if (PenPan == TRUE) {
    Surface <- paste("user-defined, albedo =", alpha)
  } else {
    Surface <- paste("not specified, actual Class-A pan evaporation data is used")
  }

  if (solar == "data") {
    message1 <- "Solar radiation data have been used for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }

  if (PenPan == TRUE) {
    message5 <- "PenPan formulation has been used to estimate Class-A pan evaporation for the calculation of potential evapotranspiration"
  } else {
    message5 <- "Class-A pan evaporation has been used for the calculation of potential evapotranspiration"
  }

  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type, message1=message1, message5=message5)
  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {
    # Generate summary message for results

    message(ET_formulation, " ", ET_type)
    message("Evaporative surface: ", Surface)
    message(message1)
    message(message5)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_ChapmanAustralian.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_ChapmanAustralian.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.JensenHaise <- function(data, constants, ts="daily",solar="sunshine hours",
                           message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }
  if (solar == "data" & is.null(data$Rs)) { # solar radiation data is required
    stop("Required data missing for 'Rs'")
  } else if (solar == "sunshine hours" & is.null(data$n)) { # for alternative calculation of solar radiation with sunshine hour
    stop("Required data missing for 'n'")
  } else if (solar == "cloud" & is.null(data$Cd)) { # for alternative calculation of sunshine hours using cloud cover
    stop("Required data missing for 'Cd'")
  } else if (solar == "monthly precipitation" & is.null(data$Precip)) { # for alternative calculation of cloudiness using monthly precipitation
    stop("Required data missing for 'Precip'")
  }


  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)

  if (solar == "data") {
    R_s <- data$Rs
  }else if (solar != "monthly precipitation"&
            solar != "cloud") {
    R_s <- (constants$as + constants$bs * (data$n/N)) * R_a
  }else {
    R_s <- (0.85 - 0.047 * data$Cd) * R_a
  }

  # estimating evapotranspiration using Jensen-Haise

  ET_JH.Daily <- 0.025* (Ta + 3) * R_s / constants$lambda # Jensen-Haise daily evapotranspiration by Prudhomme & Williams 2013
  ET.Daily <- ET_JH.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)


  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {

    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  ET_formulation <- "Jensen-Haise"
  ET_type <- "Potential ET"


  if (solar == "data") {
    message1 <- "Solar radiation data have been used for calculating evapotranspiration"
  } else if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  } else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  } else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }
  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type)
  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }

  if (message =="yes") {
    # Generate summary message for results

    message(ET_formulation, " ", ET_type)
    message(message1)
    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_JensenHaise.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_JensenHaise.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-------------------------------------------------------------------------------------

ET.McGuinnessBordne <- function(data, constants, ts="daily",
                                message = "yes",AdditionalStats= "yes", save.csv ="no",...) {
  #class(data) <- funname

  # Check of specific data requirement
  if (is.null(data$Tmax)|is.null(data$Tmin)) {
    stop("Required data missing for 'Tmax' and 'Tmin', or 'Temp'")
  }

  # Calculating mean temperature
  Ta <- (data$Tmax + data$Tmin) / 2   # Equation S2.1 in Tom McMahon's HESS 2013 paper, which in turn was based on Equation 9 in Allen et al, 1998.

  P <- 101.3 * ((293 - 0.0065 * constants$Elev) / 293)^5.26 # atmospheric pressure (S2.10)
  delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta+237.3))) / ((Ta + 237.3)^2) # slope of vapour pressure curve (S2.4)
  gamma <- 0.00163 * P / constants$lambda # psychrometric constant (S2.9)
  d_r2 <- 1 + 0.033*cos(2*pi/365 * data$J) # dr is the inverse relative distance Earth-Sun (S3.6)
  delta2 <- 0.409 * sin(2*pi/365 * data$J - 1.39) # solar dedication (S3.7)
  w_s <- acos(-tan(constants$lat_rad) * tan(delta2))  # sunset hour angle (S3.8)
  N <- 24/pi * w_s # calculating daily values
  R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) * sin(delta2) + cos(constants$lat_rad) * cos(delta2) * sin(w_s)) # extraterristrial radiation (S3.5)

  # estimating evapotranspiration using McGuinness-Bordne

  ET_MB.Daily <- R_a * (Ta + 5)/ (constants$lambda*68)  # McGuinness-Bordne daily evapotranspiration by McGuinness-Bordne  (1972) (mm.day^-1) (Oudin et al., 2005

  ET.Daily <- ET_MB.Daily
  ET.Monthly <- aggregate(ET.Daily, as.yearmon(data$Date.daily, "%m/%y"), FUN = sum)
  ET.Annual <- aggregate(ET.Daily, floor(as.numeric(as.yearmon(data$Date.daily, "%m/%y"))), FUN = sum)

  ET.MonthlyAve <- ET.AnnualAve <- NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.daily)$mon):max(as.POSIXlt(data$Date.daily)$mon)){
      i = mon - min(as.POSIXlt(data$Date.daily)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$mon== mon])
    }
    for (year in min(as.POSIXlt(data$Date.daily)$year):max(as.POSIXlt(data$Date.daily)$year)){
      i = year - min(as.POSIXlt(data$Date.daily)$year) + 1
      ET.AnnualAve[i] <- mean(ET.Daily[as.POSIXlt(data$Date.daily)$year== year])
    }

  }

  # Generate summary message for results
  ET_formulation <- "McGuinness-Bordne"
  ET_type <- "Potential ET"

  results <- list(ET.Daily=ET.Daily, ET.Monthly=ET.Monthly, ET.Annual=ET.Annual, ET.MonthlyAve=ET.MonthlyAve, ET.AnnualAve=ET.AnnualAve, ET_formulation=ET_formulation, ET_type=ET_type)
  if (ts=="daily") {
    res_ts <-  ET.Daily
  } else if (ts=="monthly") {
    res_ts <- ET.Monthly
  } else if (ts=="annual") {
    res_ts <- ET.Annual
  }
  if (message =="yes") {
    message(ET_formulation, " ", ET_type)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ", length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ",round(mean(res_ts,na.rm=T),digits=2))
      message("Max: ",round(max(res_ts,na.rm=T),digits=2))
      message("Min: ",round(min(res_ts,na.rm=T),digits=2))
    } else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ",round(mean(res_ts),digits=2))
      message("Max: ",round(max(res_ts),digits=2))
      message("Min: ",round(min(res_ts),digits=2))
    }
    #class(results) <- funname
  }
  if (save.csv == "yes") {
    # write to csv file
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file="ET_McGuinnessBordne.csv", dec=".",
                  quote=FALSE, col.names=FALSE, row.names=F, append=TRUE, sep=",")
      write.table(data.frame(get(namer,results)), file="ET_McGuinnessBordne.csv", col.names=F, append= T, sep=',' )
    }
    invisible(results)
  } else {
    return(results)
  }

}
#####################################################

# Calculate radiation variables
Radiation <- function (data, constants, ts = "monthly", solar = "sunshine hours",
                       Tdew = T, alpha = NULL, model = "CRAE")
{
  if (model != "CRWE" & model != "CRAE" & model != "CRLE") {
    stop("Error: please choose either 'CRWE' or 'CRAE' for the model to use")
  }
  if (ts == "daily") {
    stop("Error: Morton models are not available for daily time step")
  }
  if (is.null(data$Tmax)) {
    stop("Required data missing for 'Tmax' or 'Temp'")
  }
  if (is.null(data$Tmin)) {
    stop("Required data missing for 'Tmin' or 'Temp'")
  }
  if (Tdew == TRUE & is.null(data$Tdew)) {
    stop("Required data missing for 'Tdew'")
  }
  if (Tdew == FALSE) {
    if (is.null(data$va)) {
      if (is.null(data$RHmax) | is.null(data$RHmin)) {
        stop("Required data missing for 'va' , or 'RHmax' and 'RHmin' (or 'RH')")
      }
    }
  }

  if (solar == "sunshine hours" & is.null(data$n)) {
    stop("Required data missing for 'n'")
  }
  if (solar == "monthly precipitation") {
    stop("Only 'data', 'sunshine hours' and 'cloud' are accepted because estimations of sunshine hours is required")
  }
  if (solar == "sunshine hours" & is.null(data$Precip)) {
    if ("PA" %in% names(constants) == FALSE) {
      stop("Required data missing for 'Precip.daily' or required constant missing for 'PA'")
    }
  }
  T_Mo.temp <- (data$Tmax + data$Tmin)/2
  T_Mo <- aggregate(T_Mo.temp, as.yearmon(data$Date.daily,
                                          "%m/%y"), FUN = mean)
  # calculate Tdew
  if (Tdew == TRUE) {
    Tdew_Mo <- aggregate(data$Tdew, as.yearmon(data$Date.daily,
                                               "%d/%m/%y"), FUN = mean)
  } else if (!is.null(data$va)) {
    vabar_Mo <- aggregate(data$va, as.yearmon(data$Date.daily,
                                              "%m/%y"), FUN = mean)
    Tdew_Mo <- (116.9 + 237.3 * log(vabar_Mo))/(16.78 -
                                                  log(vabar_Mo))
  } else {
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax/(data$Tmax + 
                                                 237.3))
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin/(data$Tmin + 
                                                 237.3))
    vabar <- (vs_Tmin * data$RHmax/100 + vs_Tmax * data$RHmin/100)/2
    vabar_Mo <- aggregate(vabar, as.yearmon(data$Date.daily,
                                            "%m/%y"), FUN = mean)
    Tdew_Mo <- (116.9 + 237.3 * log(vabar_Mo))/(16.78 -
                                                  log(vabar_Mo))
  }

  # calculate vas
  if (is.null(data$vs)) {
    vas <- data$vs
  } else {
    vs_Tmax <- 0.6108 * exp(17.27 * data$Tmax/(data$Tmax +
                                                 237.3))
    vs_Tmin <- 0.6108 * exp(17.27 * data$Tmin/(data$Tmin +
                                                 237.3))
    vas <- (vs_Tmax + vs_Tmin)/2
  }


  delta <- 4098 * (0.6108 * exp(17.27 * T_Mo/(T_Mo + 237.3)))/(T_Mo +
                                                                 237.3)^2
  deltas <- 0.409 * sin(2 * pi/365 * data$J - 1.39)
  omegas <- acos(-tan(constants$lat_rad) * tan(deltas))
  if (!is.null(data$Precip)) {
    PA <- mean(aggregate(data$Precip, floor(as.numeric(as.yearmon(data$Date.daily,
                                                                  "%m/%y"))), FUN = sum))
  }
  else {
    PA <- constants$PA
  }
  if (model == "CRLE" | model == "CRWE") {
    constants$epsilonMo <- 0.97
    constants$fz <- 25
    constants$b0 <- 1.12
    constants$b1 <- 13
    constants$b2 <- 1.12
  }
  ptops <- ((288 - 0.0065 * constants$Elev)/288)^5.256
  alpha_zd <- 0.26 - 0.00012 * PA * sqrt(ptops) * (1 +
                                                     abs(constants$lat/42) + (constants$lat/42)^2)
  if (alpha_zd < 0.11) {
    alpha_zd <- 0.11
  }
  else if (alpha_zd > 0.17) {
    alpha_zd <- 0.17
  }

  if (solar == "sunshine hours" | solar == "cloud") {
    N <- 24/pi * omegas
    # S_daily <- data$n/N
    # for (i in 1:length(S_daily)) {
    #   if (S_daily[i] > 1) {
    #     S_daily[i] <- 1
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    S_daily <- pmin(data$n/N, 1)
    S <- mean(S_daily)


    vD_Mo <- 6.11 * exp(constants$alphaMo * Tdew_Mo/(Tdew_Mo +
                                                       constants$betaMo))
    v_Mo <- 6.11 * exp(constants$alphaMo * T_Mo/(T_Mo +
                                                   constants$betaMo))
    deltaMo <- constants$alphaMo * constants$betaMo * v_Mo/((T_Mo +
                                                               constants$betaMo)^2)
    thetaMo <- (23.2 * sin((29.5 * data$i - 94) * pi/180)) *
      pi/180
    Z_Mo <- acos(cos(constants$lat_rad - thetaMo))
    # for (i in 1:length(Z_Mo)) {
    #   if (cos(Z_Mo[i]) < 0.001) {
    #     Z_Mo[i] <- acos(0.001)
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    filt <- which(cos(Z_Mo) < 0.001)
    if (length(filt)>0) {
      Z_Mo[filt] <- acos(0.001)
    }


    omegaMo <- acos(1 - cos(Z_Mo)/(cos(constants$lat_rad) *
                                     cos(thetaMo)))
    cosz <- cos(Z_Mo) + (sin(omegaMo)/omegaMo - 1) * cos(constants$lat_rad) *
      cos(thetaMo)
    etaMo <- 1 + 1/60 * sin((29.5 * data$i - 106) * pi/180)
    G_E <- 1354/(etaMo^2) * omegaMo/pi * cosz
    alpha_zz <- matrix(NA, length(v_Mo), 1)

    if (model == "CRLE" | model == "CRWE") {
      alpha_zz[1:length(v_Mo)] <- 0.05
    } else {
      alpha_zz[1:length(v_Mo)] <- alpha_zd
      # for (i in 1:length(v_Mo)) {
      #   if (alpha_zz[i] < 0.11) {
      #     alpha_zz[i] <- 0.11
      #   }
      #   else {
      #     if (alpha_zz[i] > 0.5 * (0.91 - vD_Mo[i]/v_Mo[i])) {
      #       alpha_zz[i] <- 0.91 - vD_Mo[i]/v_Mo[i]
      #     }
      #     else {
      #     }
      #   }
      # }
      # Vectorised version : Tim Peterson 23 May 2019
      # IS 0.91 CORRECT?? 'DATA' VERSION BELOW USES 0.5
      # DG - should have the 0.5
      alpha_zz <- pmax(alpha_zz, 0.11)
      filt <- which(alpha_zz > (0.5 * (0.91 - vD_Mo/v_Mo)))
      if (length(filt)>0) {
        alpha_zz[filt] <- 0.5 * (0.91 - vD_Mo[filt]/v_Mo[filt])
      }
    }
    c_0 <- as.vector(v_Mo - vD_Mo)
    # for (i in 1:length(c_0)) {
    #   if (c_0[i] < 0) {
    #     c_0[i] <- 0
    #   }
    #   else {
    #     if (c_0[i] > 1) {
    #       c_0[i] <- 1
    #     }
    #     else {
    #       c_0[i] <- c_0[i]
    #     }
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    c_0 <- pmax(c_0, 0)
    c_0 <- pmin(c_0, 1)

    alpha_z <- alpha_zz + (1 - c_0^2) * (0.34 - alpha_zz)
    alpha_0 <- alpha_z * (exp(1.08) - ((2.16 * cos(Z_Mo))/pi +
                                         sin(Z_Mo)) * exp(0.012 * Z_Mo * 180/pi))/(1.473 *
                                                                                     (1 - sin(Z_Mo)))
    W_Mo <- vD_Mo/(0.49 + T_Mo/129)
    c_1 <- as.vector(21 - T_Mo)
    # for (i in 1:length(c_1)) {
    #   if (c_1[i] < 0) {
    #     c_1[i] <- 0
    #   }
    #   else {
    #     if (c_1[i] > 5) {
    #       c_1[i] <- 5
    #     }
    #     else {
    #       c_1[i] <- c_1[i]
    #     }
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    c_1 <- pmax(c_1, 0)
    c_1 <- pmin(c_1, 5)

    j_Mo <- (0.5 + 2.5 * (cosz)^2) * exp(c_1 * (ptops -
                                                  1))
    tauMo <- exp(-0.089 * (ptops * 1/cosz)^0.75 - 0.083 *
                   (j_Mo/cosz)^0.9 - 0.029 * (W_Mo/cosz)^0.6)
    tauaMo <- as.vector(exp(-0.0415 * (j_Mo/cosz)^0.9 -
                              (0.0029)^0.5 * (W_Mo/cosz)^0.3))
    for (i in 1:length(tauaMo)) {
      if (tauaMo[i] < exp(-0.0415 * (as.matrix(j_Mo/cosz)[i])^0.9 -
                          0.029 * (as.matrix(W_Mo/cosz)[i])^0.6)) {
        tauaMo[i] <- exp(-0.0415 * (as.matrix(j_Mo/cosz)[i])^0.9 -
                           0.029 * (as.matrix(W_Mo/cosz)[i])^0.6)
      }
      else {
        tauaMo[i] <- tauaMo[i]
      }
    }
    G_0 <- G_E * tauMo * (1 + (1 - tauMo/tauaMo) * (1 +
                                                      alpha_0 * tauMo))
    G_Mo <- S * G_0 + (0.08 + 0.3 * S) * (1 - S) * G_E
    alpha_Mo <- alpha_0 * (S + (1 - S) * (1 - Z_Mo/330 *
                                            180/pi))
    c_2 <- as.vector(10 * (vD_Mo/v_Mo - S - 0.42))
    # for (i in 1:length(c_2)) {
    #   if (c_2[i] < 0) {
    #     c_2[i] <- 0
    #   }
    #   else {
    #     if (c_2[i] > 1) {
    #       c_2[i] <- 1
    #     }
    #     else {
    #       c_2[i] <- c_2[i]
    #     }
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    c_2 <- pmax(c_2, 0)
    c_2 <- pmin(c_2, 1)

    rouMo <- 0.18 * ((1 - c_2) * (1 - S)^2 + c_2 * (1 -
                                                      S)^0.5) * 1/ptops
    B_Mo <- as.vector(constants$epsilonMo * constants$sigmaMo *
                        (T_Mo + 273)^4 * (1 - (0.71 + 0.007 * vD_Mo * ptops) *
                                            (1 + rouMo)))
    # for (i in 1:length(B_Mo)) {
    #   if (B_Mo[i] < 0.05 * constants$epsilonMo * constants$sigmaMo *
    #       (T_Mo[i] + 274)^4) {
    #     B_Mo[i] <- 0.05 * constants$epsilonMo * constants$sigmaMo *
    #       (T_Mo[i] + 274)^4
    #   }
    #   else {
    #     B_Mo[i] <- B_Mo[i]
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    B_Mo <- pmax(B_Mo, 0.05 * constants$epsilonMo * constants$sigmaMo * (T_Mo + 274)^4)


    R_T <- (1 - alpha_Mo) * G_Mo - B_Mo
  }
  else if (solar == "data") {
    vD_Mo <- 6.11 * exp(constants$alphaMo * Tdew_Mo/(Tdew_Mo +
                                                       constants$betaMo))
    v_Mo <- 6.11 * exp(constants$alphaMo * T_Mo/(T_Mo +
                                                   constants$betaMo))
    ptops <- ((288 - 0.0065 * constants$Elev)/288)^5.256
    deltaMo <- constants$alphaMo * constants$betaMo * v_Mo/((T_Mo +
                                                               constants$betaMo)^2)
    thetaMo <- (23.2 * sin((29.5 * data$i - 94) * pi/180)) *
      pi/180
    Z_Mo <- acos(cos(constants$lat_rad - thetaMo))
    # for (i in 1:length(Z_Mo)) {
    #   if (cos(Z_Mo[i]) < 0.001) {
    #     Z_Mo[i] <- acos(0.001)
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    filt <- which(cos(Z_Mo) < 0.001)
    if (length(filt)>0) {
      Z_Mo[filt] <- acos(0.001)
    }

    omegaMo <- acos(1 - cos(Z_Mo)/(cos(constants$lat_rad) *
                                     cos(thetaMo)))
    cosz <- cos(Z_Mo) + (sin(omegaMo)/omegaMo - 1) * cos(constants$lat_rad) *
      cos(thetaMo)
    etaMo <- 1 + 1/60 * sin((29.5 * data$i - 106) * pi/180)
    G_E <- 1354/(etaMo^2) * omegaMo/pi * cosz
    G_Mo <- aggregate(data$Rs * 10^6/86400, as.yearmon(data$Date.daily,
                                                       "%m/%y"), FUN = mean)
    #S = NULL
    Ta <- (data$Tmax + data$Tmin)/2
    P <- 101.3 * ((293 - 0.0065 * constants$Elev)/293)^5.26
    delta <- 4098 * (0.6108 * exp((17.27 * Ta)/(Ta + 237.3)))/((Ta +
                                                                  237.3)^2)
    gamma <- 0.00163 * P/constants$lambda
    d_r2 <- 1 + 0.033 * cos(2 * pi/365 * data$J)
    delta2 <- 0.409 * sin(2 * pi/365 * data$J - 1.39)
    w_s <- acos(-tan(constants$lat_rad) * tan(delta2))
    N <- 24/pi * w_s
    R_a <- (1440/pi) * d_r2 * constants$Gsc * (w_s * sin(constants$lat_rad) *
                                                 sin(delta2) + cos(constants$lat_rad) * cos(delta2) *
                                                 sin(w_s))
    R_so <- (0.75 + (2 * 10^-5) * constants$Elev) * R_a

    alpha_zz <- matrix(NA, length(v_Mo), 1)
    if (model == "CRLE" | model == "CRWE") {
      alpha_zz[1:length(v_Mo)] <- 0.05
    } else {
      alpha_zz[1:length(v_Mo)] <- alpha_zd
      #for (i in 1:length(v_Mo)) {
      #   if (alpha_zz[i] < 0.11) {
      #     alpha_zz[i] <- 0.11
      #   }
      #   else {
      #     if (alpha_zz[i] > 0.5 * (0.91 - vD_Mo[i]/v_Mo[i])) {
      #       alpha_zz[i] <- 0.5 * (0.91 - vD_Mo[i]/v_Mo[i])
      #     }
      #     else {
      #     }
      #   }
      # }
      # Vectorised version : Tim Peterson 23 May 2019
      alpha_zz <- pmax(alpha_zz, 0.11)
      filt <- which(alpha_zz > (0.5 * (0.91 - vD_Mo/v_Mo)))
      if (length(filt)>0) {
        alpha_zz[filt] <- 0.5 * (0.91 - vD_Mo[filt]/v_Mo[filt])
      }


    }
    c_0 <- as.vector(v_Mo - vD_Mo)
    # for (i in 1:length(c_0)) {
    #   if (c_0[i] < 0) {
    #     c_0[i] <- 0
    #   }
    #   else {
    #     if (c_0[i] > 1) {
    #       c_0[i] <- 1
    #     }
    #     else {
    #       c_0[i] <- c_0[i]
    #     }
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    c_0 <- pmax(c_0, 0)
    c_0 <- pmin(c_0, 1)

    alpha_z <- alpha_zz + (1 - c_0^2) * (0.34 - alpha_zz)
    alpha_0 <- alpha_z * (exp(1.08) - ((2.16 * cos(Z_Mo))/pi +
                                         sin(Z_Mo)) * exp(0.012 * Z_Mo * 180/pi))/(1.473 *
                                                                                     (1 - sin(Z_Mo)))
    W_Mo <- vD_Mo/(0.49 + T_Mo/129)
    c_1 <- as.vector(21 - T_Mo)
    # for (i in 1:length(c_1)) {
    #   if (c_1[i] < 0) {
    #     c_1[i] <- 0
    #   }
    #   else {
    #     if (c_1[i] > 5) {
    #       c_1[i] <- 5
    #     }
    #     else {
    #       c_1[i] <- c_1[i]
    #     }
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    c_1 <- pmax(c_1, 0)
    c_1 <- pmin(c_1, 5)

    j_Mo <- (0.5 + 2.5 * (cosz)^2) * exp(c_1 * (ptops -
                                                  1))
    tauMo <- exp(-0.089 * (ptops * 1/cosz)^0.75 - 0.083 *
                   (j_Mo/cosz)^0.9 - 0.029 * (W_Mo/cosz)^0.6)
    tauaMo <- as.vector(exp(-0.0415 * (j_Mo/cosz)^0.9 -
                              (0.0029)^0.5 * (W_Mo/cosz)^0.3))
    # for (i in 1:length(tauaMo)) {
    #   if (tauaMo[i] < exp(-0.0415 * (as.matrix(j_Mo/cosz)[i])^0.9 -
    #                       0.029 * (as.matrix(W_Mo/cosz)[i])^0.6)) {
    #     tauaMo[i] <- exp(-0.0415 * (as.matrix(j_Mo/cosz)[i])^0.9 -
    #                        0.029 * (as.matrix(W_Mo/cosz)[i])^0.6)
    #   }
    #   else {
    #     tauaMo[i] <- tauaMo[i]
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    tauaMo <- pmax(tauaMo, exp(-0.0415 * (as.matrix(j_Mo/cosz))^0.9 -
                                 0.029 * (as.matrix(W_Mo/cosz))^0.6))

    G_0 <- G_E * tauMo * (1 + (1 - tauMo/tauaMo) * (1 +
                                                      alpha_0 * tauMo))


    S <- 0.53*G_Mo / (G_0 - 0.47*G_Mo)
    S[which(S<0)] <- 0
    S[which(S>1)] <- 1

    alpha_Mo <- alpha_0 * (S + (1 - S) * (1 - Z_Mo/330 *
                                            180/pi))
    c_2 <- as.vector(10 * (vD_Mo/v_Mo - S - 0.42))
    # for (i in 1:length(c_2)) {
    #   if (c_2[i] < 0) {
    #     c_2[i] <- 0
    #   }
    #   else {
    #     if (c_2[i] > 1) {
    #       c_2[i] <- 1
    #     }
    #     else {
    #       c_2[i] <- c_2[i]
    #     }
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    c_2 <- pmax(c_2, 0)
    c_2 <- pmin(c_2, 1)
    rouMo <- 0.18 * ((1 - c_2) * (1 - S)^2 + c_2 * (1 -
                                                      S)^0.5) * 1/ptops
    B_Mo <- as.vector(constants$epsilonMo * constants$sigmaMo *
                        (T_Mo + 273)^4 * (1 - (0.71 + 0.007 * vD_Mo * ptops) *
                                            (1 + rouMo)))
    # for (i in 1:length(B_Mo)) {
    #   if (B_Mo[i] < 0.05 * constants$epsilonMo * constants$sigmaMo *
    #       (T_Mo[i] + 274)^4) {
    #     B_Mo[i] <- 0.05 * constants$epsilonMo * constants$sigmaMo *
    #       (T_Mo[i] + 274)^4
    #   }
    #   else {
    #     B_Mo[i] <- B_Mo[i]
    #   }
    # }
    # Vectorised version : Tim Peterson 23 May 2019
    B_Mo <- pmax(B_Mo, 0.05 * constants$epsilonMo * constants$sigmaMo * (T_Mo + 274)^4)

    R_T <- (1 - alpha_Mo) * G_Mo - B_Mo
  }
  if (solar == "sunshine hours") {
    message1 <- "Sunshine hour data have been used for calculating incoming solar radiation"
  }
  else if (solar == "cloud") {
    message1 <- "Cloudiness data have been used for calculating sunshine hour and thus incoming solar radiation"
  }
  else {
    message1 <- "Monthly precipitation data have been used for calculating incoming solar radiation"
  }
  if (Tdew == TRUE) {
    message6 <- "Data of dew point temperature has been used"
  }
  else {
    message6 <- "Data of average vapour pressure has been used to estimate dew point pressure"
  }
  variables <- list(T_Mo = T_Mo, Tdew_Mo = Tdew_Mo, S = S,
                    R_T = R_T, ptops = ptops, vD_Mo = vD_Mo, v_Mo = v_Mo,
                    deltaMo = deltaMo, G_E = G_E, G_Mo = G_Mo, alpha_Mo = alpha_Mo,
                    B_Mo = B_Mo, message1 = message1, message6 = message6)
  return(variables)
}


#-------------------------------------------------------------------------------------
ET.MortonCRAE <- function (data, constants, ts = "monthly", est = "potential ET",
                           solar = "sunshine hours", Tdew = T, alpha = NULL,
                           message = "yes",AdditionalStats= "yes", save.csv ="no",...)
{
  variables <- Radiation(data, constants, ts, solar, Tdew,
                         alpha)
  R_T <- (1 - variables$alpha_Mo) * variables$G_Mo - variables$B_Mo
  R_TC <- as.vector(R_T)
  for (i in 1:length(R_TC)) {
    if (R_TC[i] < 0) {
      R_TC[i] <- 0
    }
    else {
      R_TC[i] <- R_TC[i]
    }
  }
  xiMo <- 1/(0.28 * (1 + variables$vD_Mo/variables$v_Mo) +
               R_TC * variables$deltaMo/(variables$ptops * constants$gammaps *
                                           (1/variables$ptops)^0.5 * constants$b0 * constants$fz *
                                           (variables$v_Mo - variables$vD_Mo)))
  # for (i in 1:length(xiMo)) {
  #   if (xiMo[i] < 1) {
  #     xiMo[i] <- 1
  #   }
  #   else {
  #     xiMo[i] <- xiMo[i]
  #   }
  # }
  # Vectorised version : Tim Peterson 23 May 2019
  xiMo <- pmax(xiMo, 1)

  f_T <- (1/variables$ptops)^0.5 * constants$fz/xiMo
  lambdaMo1 <- constants$gammaps * variables$ptops + 4 * constants$epsilonMo *
    constants$sigmaMo * (variables$T_Mo + 274)^3/f_T
  T_p <- variables$T_Mo
  for (i in 1:99999) {
    v_p <- 6.11 * exp((constants$alphaMo * T_p)/(T_p + constants$betaMo))
    delta_p <- constants$alphaMo * constants$betaMo * v_p/((T_p +
                                                              constants$betaMo)^2)
    delta_T_p <- (R_T/f_T + variables$vD_Mo - v_p + lambdaMo1 *
                    (variables$T_Mo - T_p))/(delta_p + lambdaMo1)
    T_p <- T_p + delta_T_p
    if (abs(max(na.omit(delta_T_p))) < 0.01)
      break
  }
  v_p <- 6.11 * exp((constants$alphaMo * T_p)/(T_p + constants$betaMo))
  delta_p <- constants$alphaMo * constants$betaMo * v_p/((T_p +
                                                            constants$betaMo)^2)
  E_TP.temp <- R_T - lambdaMo1 * f_T * (T_p - variables$T_Mo)
  R_TP <- E_TP.temp + variables$ptops * constants$gammaps *
    f_T * (T_p - variables$T_Mo)
  E_TW.temp <- constants$b1 + constants$b2 * R_TP/(1 + variables$ptops *
                                                     constants$gammaps/delta_p)
  E_T_Mo.temp <- 2 * E_TW.temp - E_TP.temp
  E_TP.temp <- 1/(constants$lambdaMo) * E_TP.temp
  E_TW.temp <- 1/(constants$lambdaMo) * E_TW.temp
  E_T_Mo.temp <- 1/(constants$lambdaMo) * E_T_Mo.temp
  E_TP <- E_TP.temp * data$Ndays
  E_TW <- E_TW.temp * data$Ndays
  E_T_Mo <- E_T_Mo.temp * data$Ndays
  if (est == "potential ET") {
    ET_Mo.Monthly <- E_TP
    ET_Mo.Average <- E_TP.temp
    ET_type <- "Potential ET"
  }
  else if (est == "wet areal ET") {
    ET_Mo.Monthly <- E_TW
    ET_Mo.Average <- E_TW.temp
    ET_type <- "Wet-environment Areal ET"
  }
  else if (est == "actual areal ET") {
    ET_Mo.Monthly <- E_T_Mo
    ET_Mo.Average <- E_T_Mo.temp
    ET_type <- "Actual Areal ET"
  }
  ET.Daily <- NULL
  ET.Monthly <- ET_Mo.Monthly
  ET.Annual <- aggregate(ET.Monthly, floor(as.numeric(as.yearmon(data$Date.monthly,
                                                                 "%m/%y"))), FUN = sum)
  ET.MonthlyAve = ET.AnnualAve = NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.monthly)$mon):max(as.POSIXlt(data$Date.monthly)$mon)) {
      i = mon - min(as.POSIXlt(data$Date.monthly)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET_Mo.Average[as.POSIXlt(data$Date.monthly)$mon ==
                                               mon])
    }
    for (year in min(as.POSIXlt(data$Date.monthly)$year):max(as.POSIXlt(data$Date.monthly)$year)) {
      i = year - min(as.POSIXlt(data$Date.monthly)$year) +
        1
      ET.AnnualAve[i] <- mean(ET_Mo.Average[as.POSIXlt(data$Date.monthly)$year ==
                                              year])
    }

  }

  ET_formulation <- "Morton CRAE"

  results <- list(ET.Daily = ET.Daily, ET.Monthly = ET.Monthly,
                  ET.Annual = ET.Annual, ET.MonthlyAve = ET.MonthlyAve,
                  ET.AnnualAve = ET.AnnualAve, ET_formulation = ET_formulation,
                  ET_type = ET_type, message1 = variables$message1, message6 = variables$message6)
  if (ts == "monthly") {
    res_ts <- ET.Monthly
  }
  else if (ts == "annual") {
    res_ts <- ET.Annual
  }

  if (message == "yes") {
    message(ET_formulation, " ", ET_type)
    message(variables$message1)
    message(variables$message6)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ",
              length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ", round(mean(res_ts, na.rm = T), digits = 2))
      message("Max: ", round(max(res_ts, na.rm = T), digits = 2))
      message("Min: ", round(min(res_ts, na.rm = T), digits = 2))
    }
    else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ", round(mean(res_ts), digits = 2))
      message("Max: ", round(max(res_ts), digits = 2))
      message("Min: ", round(min(res_ts), digits = 2))
    }
  }
  if (save.csv== "yes") {
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file = "ET_MortonCRAE.csv",
                  dec = ".", quote = FALSE, col.names = FALSE, row.names = F,
                  append = TRUE, sep = ",")
      write.table(data.frame(get(namer, results)), file = "ET_MortonCRAE.csv",
                  col.names = F, append = T, sep = ",")
    }
    invisible(results)
  } else {
    return(results)
  }

}

#-----------------------------------------------------------------------------------
ET.MortonCRWE <- function (data, constants, ts = "monthly", est = "potential ET",
                           solar = "sunshine hours", Tdew = T, alpha = NULL,
                           message = "yes",AdditionalStats= "yes", save.csv ="no",...)
{
  constants$epsilonMo <- 0.97
  constants$fz <- 25
  constants$b0 <- 1.12
  constants$b1 <- 13
  constants$b2 <- 1.12

  variables <- Radiation(data, constants, ts, solar, Tdew,
                         alpha, model = "CRWE")
  R_TC <- as.vector(variables$R_T)
  # for (i in 1:length(R_TC)) {
  #   if (R_TC[i] < 0) {
  #     R_TC[i] <- 0
  #   }
  #   else {
  #     R_TC[i] <- R_TC[i]
  #   }
  # }
  #
  # Vectorised version : Tim Peterson 23 May 2019
  R_TC = pmax(R_TC, 0)

  xiMo <- 1/(0.28 * (1 + variables$vD_Mo/variables$v_Mo) +
               R_TC * variables$deltaMo/(variables$ptops * constants$gammaps *
                                           (1/variables$ptops)^0.5 * constants$b0 * constants$fz *
                                           (variables$v_Mo - variables$vD_Mo)))
  # for (i in 1:length(xiMo)) {
  #   if (xiMo[i] < 1) {
  #     xiMo[i] <- 1
  #   }
  #   else {
  #     xiMo[i] <- xiMo[i]
  #   }
  # }
  # Vectorised version : Tim Peterson 23 May 2019
  xiMo <- pmax(xiMo, 1)

  f_T <- (1/variables$ptops)^0.5 * constants$fz/xiMo
  lambdaMo1 <- constants$gammaps * variables$ptops + 4 * constants$epsilonMo *
    constants$sigmaMo * (variables$T_Mo + 274)^3/f_T
  T_p <- variables$T_Mo
  for (i in 1:99999) {
    v_p <- 6.11 * exp((constants$alphaMo * T_p)/(T_p + constants$betaMo))
    delta_p <- constants$alphaMo * constants$betaMo * v_p/((T_p +
                                                              constants$betaMo)^2)
    delta_T_p <- (R_TC/f_T + variables$vD_Mo - v_p + lambdaMo1 *
                    (variables$T_Mo - T_p))/(delta_p + lambdaMo1)
    T_p <- T_p + delta_T_p
    if (abs(max(na.omit(delta_T_p))) < 0.01)
      break
  }
  v_p <- 6.11 * exp((constants$alphaMo * T_p)/(T_p + constants$betaMo))
  delta_p <- constants$alphaMo * constants$betaMo * v_p/((T_p +
                                                            constants$betaMo)^2)
  E_P.temp <- R_TC - lambdaMo1 * f_T * (T_p - variables$T_Mo)
  R_P <- E_P.temp + variables$ptops * constants$gammaps *
    f_T * (T_p - variables$T_Mo)
  E_W.temp <- constants$b1 + constants$b2 * R_P/(1 + variables$ptops *
                                                   constants$gammaps/delta_p)
  E_T_Mo.temp <- 2 * E_W.temp - E_P.temp
  E_P.temp <- 1/(constants$lambdaMo) * E_P.temp
  E_W.temp <- 1/(constants$lambdaMo) * E_W.temp
  E_T_Mo.temp <- 1/(constants$lambdaMo) * E_T_Mo.temp
  E_P <- E_P.temp * data$Ndays
  E_W <- E_W.temp * data$Ndays
  E_T_Mo <- E_T_Mo.temp * data$Ndays
  if (est == "potential ET") {
    ET_Mo.Monthly <- E_P
    ET_Mo.Average <- E_P.temp
    ET_type <- "Potential ET"
  }
  else if (est == "shallow lake ET") {
    ET_Mo.Monthly <- E_W
    ET_Mo.Average <- E_W.temp
    ET_type <- "Shallow Lake Evaporation"
  }
  ET.Daily <- NULL
  ET.Monthly <- ET_Mo.Monthly
  ET.Annual <- aggregate(ET.Monthly, floor(as.numeric(as.yearmon(data$Date.monthly,
                                                                 "%m/%y"))), FUN = sum)
  ET.MonthlyAve = ET.AnnualAve = NULL
  if (AdditionalStats=="yes") {
    for (mon in min(as.POSIXlt(data$Date.monthly)$mon):max(as.POSIXlt(data$Date.monthly)$mon)) {
      i = mon - min(as.POSIXlt(data$Date.monthly)$mon) + 1
      ET.MonthlyAve[i] <- mean(ET_Mo.Average[as.POSIXlt(data$Date.monthly)$mon ==
                                               mon])
    }
    for (year in min(as.POSIXlt(data$Date.monthly)$year):max(as.POSIXlt(data$Date.monthly)$year)) {
      i = year - min(as.POSIXlt(data$Date.monthly)$year) +
        1
      ET.AnnualAve[i] <- mean(ET_Mo.Average[as.POSIXlt(data$Date.monthly)$year ==
                                              year])
    }

  }

  ET_formulation <- "Morton CRWE"

  results <- list(ET.Daily = ET.Daily, ET.Monthly = ET.Monthly,
                  ET.Annual = ET.Annual, ET.MonthlyAve = ET.MonthlyAve,
                  ET.AnnualAve = ET.AnnualAve, ET_formulation = ET_formulation,
                  ET_type = ET_type, message1 = variables$message1, message6 = variables$message6)
  if (ts == "monthly") {
    res_ts <- ET.Monthly
  }
  else if (ts == "annual") {
    res_ts <- ET.Annual
  }
  if (message =="yes") {
    message(ET_formulation, " ", ET_type)
    message(variables$message1)
    message(variables$message6)

    message("Timestep: ", ts)
    message("Units: mm")
    message("Time duration: ", time(res_ts[1]), " to ", time(res_ts[length(res_ts)]))
    if (NA %in% res_ts) {
      message(length(res_ts), " ET estimates obtained; ",
              length(which(is.na(res_ts))), " NA output entries due to missing data")
      message("Basic stats (NA excluded)")
      message("Mean: ", round(mean(res_ts, na.rm = T), digits = 2))
      message("Max: ", round(max(res_ts, na.rm = T), digits = 2))
      message("Min: ", round(min(res_ts, na.rm = T), digits = 2))
    }
    else {
      message(length(res_ts), " ET estimates obtained")
      message("Basic stats")
      message("Mean: ", round(mean(res_ts), digits = 2))
      message("Max: ", round(max(res_ts), digits = 2))
      message("Min: ", round(min(res_ts), digits = 2))
    }
  }
  if (save.csv == "yes") {
    for (i in 1:length(results)) {
      namer <- names(results[i])
      write.table(as.character(namer), file = "ET_MortonCRWE.csv",
                  dec = ".", quote = FALSE, col.names = FALSE, row.names = F,
                  append = TRUE, sep = ",")
      write.table(data.frame(get(namer, results)), file = "ET_MortonCRWE.csv",
                  col.names = F, append = T, sep = ",")
    }
    invisible(results)
  } else {
    return(results)
  }

}

Try the Evapotranspiration package in your browser

Any scripts or data that you put into this service are public.

Evapotranspiration documentation built on Jan. 10, 2022, 9:06 a.m.