Imax <- function(InputsModel, 
                 IndPeriod_Run, 
                 TestedValues = seq(from = 0.1, to = 3, seq = 0.1)) {
  
  ##_____Arguments_check_____________________________________________________________________
  if (!inherits(InputsModel, "InputsModel")) {
    stop("'InputsModel' must be of class 'InputsModel'")
  }  
  if (!inherits(InputsModel, "hourly")) {
    stop("'InputsModel' must be of class 'hourly'")
  }  
  
  ##check_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]:Ind_Run[length(IndPeriod_Run)])) {
    stop("'IndPeriod_Run' must be a continuous sequence of integers")
  }
  
  ##TestedValues  
  if (!(is.numeric(TestedValues))) {
    stop("'TestedValues' must be 'numeric'")
  }

  
  ##aggregate data at the daily time step  
  TabSeries <- data.frame(DatesR = InputsModel$DatesR[IndPeriod_Run], 
                          Precip = InputsModel$Precip[IndPeriod_Run], 
                          PotEvap = InputsModel$PotEvap[IndPeriod_Run])
  daily_data <- SeriesAggreg(TabSeries, "hourly", "daily", 
                             ConvertFun = c("sum", "sum"))
  
  ##calculate total interception of daily GR models on the period
  cum_daily <- sum(pmin(daily_data$Precip, daily_data$PotEvap))
  
  ##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)])
  
}