plot.CalGR <- function(x, which = c("perf", "hist", "time"), ...) {
  
  if (! any(class(x) %in% "CalGR")) {
    stop("Non convenient data for x argument. Must be of class \"CalGR\"")
  }
  if (! any(which %in% c("perf", "hist", "time"))) {
    stop("Non convenient data for which argument. Must be of class \"perf\", \"hist\" or \"time\"")
  }  
  
  nbParamX <- .TypeModelGR(x)$NbParam #as.numeric(gsub("\\D", "", x$TypeModel))
  nbParamC <- ifelse(.TypeModelGR(x)$CemaNeige, 2, 0)
  nbParam <- nbParamX + nbParamC
  nmParam <- c(sprintf("X%i", 1:nbParamX), sprintf("C%i", seq_len(nbParamC)))
  
  if (any(which[1L] %in%  c("perf"))) {
    plot(x$OutputsModel, Qobs = x$Qobs, ...)
  }
  
  if (any(which[1L] %in%  c("hist"))) {
    layout.list <- list(matrix(c(1:2), ncol = 2),
                        matrix(c(1:3, 3), ncol = 2),
                        matrix(NA),
                        matrix(c(1:5, 5), ncol = 3),
                        matrix(c(1:5, 0, 6, 6), ncol = 4),
                        matrix(c(1:7, 7), ncol = 4),
                        matrix(c(1:5, 0, 6:8, 8), ncol = 5),
                        matrix(c(1:9, 9), ncol = 5))
    layout.mat <- layout.list[[nbParam]]
    layout.w <- c(rep(1, ncol(layout.mat)-1), 2)
    layout.h <- rep(1, nrow(layout.mat))
    layout(mat = layout.mat, widths = layout.w, heights = layout.h)
  }
  # if (which[1L] %in%  c("both")) {
  #   layout.mat <- rbind(layout.mat, max(layout.mat)+1, max(layout.mat)+2)
  #   layout(mat = layout.mat, widths = layout.w, heights = layout.h)
  # }
  if (which[1L] %in%  c("time")) {
    layout(mat = matrix(1:2), widths = c(1, 2), heights = c(1, 2))
  }
  
  if (any(which[1L] %in%  c("hist"))) {
    
    par(mar = c(2.5, 2.5, 4.5, 1))
    for (i in seq_len(ncol(x$OutputsCalib$HistParamR))) {
      
      plot(x$OutputsCalib$HistParamR[, i], col = ifelse(grepl("X", nmParam[i]), "darkorchid3", "darkolivegreen3"),
           type = "b", pch = 19, cex = 0.8,
           xlab = "", ylab = "", main = nmParam[i])
    }
    plot(x$OutputsCalib$HistCrit, col = "tomato",
         type = "b", pch = 19, cex = 1.0,
         xlab = "", ylab = "", main = "Efficiency criteria\nduring the steepest-descent step")
  }
  
  if (any(which[1L] %in%  c("time"))) {
    data <- data.frame(DatesR = x$OutputsModel$DatesR,
                       Precip = x$OutputsModel$Precip,
                       Qobs   = x$Qobs,
                       Qsim   = x$OutputsModel$Qsim)
    
    par(mar = c(0.1, 4, 4, 1), xaxt = "n")
    plot(Precip ~ DatesR, data = data, type = "h", col = "royalblue",
         xlab = "", ylab = "Precip.", main = "", ylim = rev(range(data$Precip)))
    
    par(mar = c(5, 4, 0.1, 1), xaxt = "s")
    matplot(data$DatesR, data[, c("Qobs", "Qsim")], type = "l", col = c("black", "orangered"), lty = 1,
         xlab = "Time", ylab = "Flow", main = "", xaxt = "n")
    axis.POSIXct(side = 1, x = data$DatesR)
    
    # if (plot.na) {
    #   axis(side = 1, at = as.POSIXct(data$DatesR[is.na(data$Qobs)]),
    #        labels = FALSE, lwd.ticks = 3, col.ticks = col.na, tck = 0.025, lend = "butt")	
    #   
    #   legend("topright", legend = "NA", pch = 15, col = col.na, bty = "n", cex = 0.8)
    # }
  }
}