calculate_Imax <- function(InputsModel, 
                           IndPeriod_Run, 
                           tested_values = seq(0.1, 3, 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'")
  }  
  if (!(class(tested_values) == "numeric")){
    stop("'tested_values' must be 'numeric'")
  }
  
  ##check_IndPeriod_Run
  if (!is.vector(IndPeriod_Run)) {
    stop("'IndPeriod_Run' must be a vector of numeric values")
  }
  if (!is.numeric(IndPeriod_Run)) {
    stop("'IndPeriod_Run' must be a vector of numeric values")
  }
  if (!identical(as.integer(IndPeriod_Run), as.integer(seq(from = IndPeriod_Run[1], to = tail(IndPeriod_Run, 1), by = 1)))) {
    stop("'IndPeriod_Run' must be a continuous sequence of integers")
  }
  if (storage.mode(IndPeriod_Run) != "integer") {
    stop("'IndPeriod_Run' should be of type integer")
  }
  
  ##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(tested_values)))
  for (Imax in tested_values) {
    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 == tested_values)] <- abs(cum_hourly - cum_daily)
  }
  
  ##return the Imax value that minimises the difference
  return(tested_values[which.min(differences)])
  
}