Imax.R

Imax <- function(InputsModel,
                 IndPeriod_Run,
                 TestedValues = seq(from = 0.1, to = 3, by = 0.1)) {


  ## ---------- check arguments

  ##  InputsModel
  if (!inherits(InputsModel, "InputsModel")) {
    stop("'InputsModel' must be of class 'InputsModel'")
  }
  if (!inherits(InputsModel, "hourly")) {
    stop("'InputsModel' must be of class 'hourly'")
  }

  ## IndPeriod_Run
  if (!is.vector(IndPeriod_Run)) {
    stop("'IndPeriod_Run' must be a vector of numeric values")
  }
  if (!inherits(IndPeriod_Run, "integer")) {
    stop("'IndPeriod_Run' must be of type integer")
  }
  if (!identical(as.integer(IndPeriod_Run), IndPeriod_Run[1]:IndPeriod_Run[length(IndPeriod_Run)])) {
    stop("'IndPeriod_Run' must be a continuous sequence of integers")
  }

  ## TestedValues
  if (!(is.numeric(TestedValues))) {
    stop("'TestedValues' must be 'numeric'")
  }


  ## ---------- hourly inputs aggregation

  ## aggregate data at the daily time step
  daily_data <- SeriesAggreg(InputsModel[IndPeriod_Run], Format = "%Y%m%d")


  ## ---------- calculate interception

  ## calculate total interception of daily GR models on the period
  cum_daily <- sum(pmin(daily_data$Precip, daily_data$PotEvap))
  if (anyNA(cum_daily)) {
    stop("'IndPeriod_Run' must be set to select 24 hours by day")
  }

  ## calculate total interception of the GR5H interception store on the period
  ## and compute difference with daily values
  differences <- array(NA, c(length(TestedValues)))
  for (Imax in TestedValues) {
    C0 <- 0
    cum_hourly <- 0
    for (i in IndPeriod_Run) {
      Ec  <- min(InputsModel$PotEvap[i], InputsModel$Precip[i] + C0)
      Pth <- max(0, InputsModel$Precip[i] - (Imax-C0) - Ec)
      C0  <- C0 + InputsModel$Precip[i] - Ec - Pth
      cum_hourly <- cum_hourly + Ec
    }
    differences[which(Imax == TestedValues)] <- abs(cum_hourly - cum_daily)
  }

  ## return the Imax value that minimises the difference
  return(TestedValues[which.min(differences)])

}