Utils.R

## =================================================================================
## commands to avoid warnings during package checking when global variables are used
## =================================================================================

if (getRversion() >= "2.15.1") {
  utils::globalVariables(c(".ShinyGR.args"))
  utils::suppressForeignCheck(c(".ShinyGR.args"))
  utils::globalVariables(c(".ShinyGR.hist"))
  utils::suppressForeignCheck(c(".ShinyGR.hist"))
}




## =================================================================================
## function to manage the model units
## =================================================================================

.TypeModelGR <- function(x) {
  
  if (!is.list(x)) {
    x <- list(TypeModel = x)
  }
  if (any(class(x) %in% c("PrepGR", "CalGR", "SimGR")) || names(x) %in% "TypeModel") {
    x <- x$TypeModel
  }
  
  StrName    <- "(.*)(GR)(\\d{1})(\\D{1})"
  NameModel  <- gsub(StrName, "\\2\\3\\4", x)
  TimeUnitFR <- gsub(StrName, "\\4", x)
  # TimeUnit   <- switch(TimeUnitFR, H = "hour", J = "day", M = "month", A = "year")
  # TimeLag    <- switch(TimeUnit, "hour" = 3600, "day" = 3600*24, "month" = 3600*24*31, "year" = 366)
  TimeUnit   <- switch(TimeUnitFR, H = "h", J = "d", M = "month", A = "y")
  TimeLag    <- switch(TimeUnit, "h" = 3600, "d" = 3600*24, "month" = 3600*24*31, "y" = 366)
  NbParam    <- gsub(StrName, "\\3", x)
  isCN       <- grepl("CemaNeige"  , x)
  
  res <- list(TypeModel = x, NameModel = NameModel, CemaNeige = isCN, 
              NbParam = as.numeric(NbParam),
              TimeUnit = TimeUnit, TimeLag = TimeLag)
  return(res)
}




## =================================================================================
## function to plot the gr models diagrams (only GR4J and GR5J)
## =================================================================================

.DiagramGR <- function(OutputsModel, Param, SimPer, EventDate, HydroModel, CemaNeige) {
  
  
  # --------------------------------------------------------------------------------
  # PARAMETRES
  # --------------------------------------------------------------------------------
  
  # Parametres
  mgp             <- c(0, 0.75, 0)
  col_P           <- rgb(066, 139, 202, maxColorValue = 255) #"royalblue"
  col_E           <- rgb(164, 196, 000, maxColorValue = 255) #"forestgreen"
  col_Q           <- "orangered"
  col_SP          <- adjustcolor("cyan4"   , alpha.f = 0.60)
  col_SR          <- adjustcolor("darkblue", alpha.f = 0.60)
  col_R           <- rgb(066, 139, 202, maxColorValue = 255) #rgb(037, 155, 210, maxColorValue = 255)
  col_mod_bg      <- rgb(245, 245, 245, maxColorValue = 255)
  col_mod_bd      <- rgb(231, 231, 231, maxColorValue = 255)
  xy_E            <- c(250, 980)
  xy_PE           <- c(250, 940)
  xy_AE           <- c(250, 860)
  xy_P            <- c(600, 980)
  xy_Precip       <- c(600, 950)
  xy_Q            <- c(700,  30)
  x_Ps            <- 440
  x_PnPs          <- 700
  y_interception  <- 900
  y_rendement     <- 815
  y_percolation   <- 575
  xy_Q9           <- c(400, 310)
  xy_Q1           <- c(800, 310)
  y_routage       <- 100
  fact_res        <- 200/800/3
  fact_resExp     <- 1
  base_res        <- 300
  NH              <- 10
  xy_min_PROD     <- c(200, 610)
  xy_min_ROUT     <- c(250, 150)
  xy_min_EXPO     <- c(200, 250)
  y_entreeUH      <- 500
  xy_UH1          <- c(500, 420)
  xy_UH2          <- c(900, 420)
  y_Ech_Q1        <- 170 #200
  y_Ech_Q9        <- 150 #180
  D               <- 5/2
  xpad            <- 1.5
  ypad            <- 1.5
  max_triangle    <- max(unlist(OutputsModel[c("Perc", "PR", "Q9", "Q1", "QR", "QD", "Pn", "Ps", 
                                               "AE", "Precip", "PotEvap", "AExch1", "AExch2")]))
  fact_var        <- 40
  fact_triangle   <- 100#25
  cex_max_poly    <- 2 # 0.005
  cex_tri         <- function(cex, fact = 25, max) suppressWarnings(log(abs(cex) * fact + 1) / max)
  radius1         <- 0
  radius2         <- 60
  tri_R           <- -0x25BA
  tri_B           <- -0x25BC
  tri_L           <- -0x25C4
  tri_T           <- -0x25B2
  
  par(col.axis = par("fg"), cex.axis = 1.3, cex.lab = 1.3, cex = 0.7, mgp = mgp)
  
  if (.GlobalEnv$.ShinyGR.args$theme == "Cyborg") {
    col_mod_bg    <- rgb(255-245, 255-245, 255-245, maxColorValue = 255)
    col_mod_bd    <- rgb(255-231, 255-231, 255-231, maxColorValue = 255)
  }
  if (.GlobalEnv$.ShinyGR.args$theme == "Flatly") {
    col_mod_bg    <- "#ECF0F1"
    col_mod_bd    <- "#ECF0F1"
  }
  
  # Pas de temps
  dates_deb       <- EventDate
  n_pdt           <- length(which(OutputsModel$DatesR >= EventDate & OutputsModel$DatesR <= SimPer[2L]))
  i_pdt           <- which(format(OutputsModel$DatesR, "%Y%m%d") == format(EventDate, "%Y%m%d"))
  
  
  
  # --------------------------------------------------------------------------------
  # UH 1 & 2
  # --------------------------------------------------------------------------------
  
  if (HydroModel %in% c("GR4J", "GR6J")) {
    # Calcul des ordonnees SH1 de l' "hydrogramme unitaire cumule" UH1
    SH1     <- array(NA, NH)
    for (i in 1:NH) {
      if (i <= 0)                  SH1[i] <- 0
      if (i > 0 & i < Param[4])    SH1[i] <- (i/Param[4])^(D)
      if (i >= Param[4])           SH1[i] <- 1
    }
    
    # Calcul des ordonnees UH1 de l' "hydrogramme unitaire discret" UH1
    UH1     <- array(NA, NH)
    for (j in 1:NH) {
      if (j == 1) {
        UH1[j] <- SH1[j]
      } else {
        UH1[j] <- SH1[j] - SH1[j-1]
      }  
    }
    
    # Parametres
    max_UH1 <- log(sqrt(max(max(UH1)*OutputsModel$PR*0.9))+1)
    
  }
  
  
  # Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" UH2
  SH2     <- array(NA, 2*NH)
  for (i in 1:(2*NH)) {
    if (i <= 0)                           SH2[i] <- 0
    if (i > 0 & i < Param[4])             SH2[i] <- 0.5*(i/Param[4])^(D)
    if (i >= Param[4] & i < 2*Param[4])   SH2[i] <- 1 - (0.5*(2-i/Param[4])^(D))
    if (i >= 2*Param[4])                  SH2[i] <- 1
  }
  
  # Calcul des ordonnees UH2 de l' "hydrogramme unitaire discret" UH2
  UH2     <- array(NA, 2*NH)
  for (j in 1:(2*NH)) {
    if (j == 1) {
      UH2[j] <- SH2[j]
    } else {
      UH2[j] <- SH2[j] - SH2[j-1]
    }  
  }
  
  # Parametres
  max_UH2 <- log(sqrt(max(max(UH2)*OutputsModel$PR*0.1))+1) 
  
  
  # --------------------------------------------------------------------------------
  # PARTITIONNEMENT FENETRE GRAPHIQUE
  # --------------------------------------------------------------------------------
  
  # layout(matrix(c(1:4, 4, 4), nrow = 3, ncol = 2, byrow = FALSE), widths = c(1.0, 0.6))
  
  # --------------------------------------------------------------------------------
  # PLUIE ET ETP
  # --------------------------------------------------------------------------------
  
  # P
  # par(mar = c(2, 4, 1, 1), mgp = mgp)
  # plot(OutputsModel$Dates, OutputsModel$Precip, type = "h", col = col_P, ylim = rev(range(OutputsModel$Precip)), xaxt = "n", ylab = "precip. [mm/d]")
  # rect(xleft = EventDate, ybottom = par("usr")[3], xright =  par("usr")[2], ytop =  par("usr")[4], col = adjustcolor(par("bg"), alpha.f = 0.75), border = NA)
  # abline(v = EventDate, col = "grey", lwd = 2, lty = 2)
  # box()
  
  # ETP
  # par(mar = c(2, 4, 1, 1), mgp = mgp)
  # plot(OutputsModel$Dates, OutputsModel$PotEvap, pch = 19, col = col_E, xaxt = "n", ylab = "evapo. [mm/d]")
  # rect(xleft = EventDate, ybottom = par("usr")[3], xright =  par("usr")[2], ytop =  par("usr")[4], col = adjustcolor(par("bg"), alpha.f = 0.75), border = NA)
  # abline(v = EventDate, col = "grey", lwd = 2, lty = 2)
  # box()  
  
  
  # --------------------------------------------------------------------------------
  # DEBIT
  # --------------------------------------------------------------------------------
  
  # Q
  # par(mar = c(2, 4, 1, 1), mgp = mgp)
  # plot(OutputsModel$Dates, OutputsModel$Qobs, type = "l", ylab = "flow [mm/d]")
  # lines(OutputsModel$Dates[1:i_pdt], OutputsModel$Qsim[1:i_pdt], type = "l", col = "orangered")
  # rect(xleft = EventDate, ybottom = par("usr")[3], xright =  par("usr")[2], ytop =  par("usr")[4], col = adjustcolor(par("bg"), alpha.f = 0.75), border = NA)
  # abline(v = EventDate, col = "grey", lwd = 2, lty = 2)
  # box()  
  
  
  # --------------------------------------------------------------------------------
  # SCHEMAS MODELES
  # --------------------------------------------------------------------------------
  
  # Cadre
  par(mar = rep(0.2, 4))
  par(fg = par("fg"))
  plot(x = 0, type = "n", xlab = "", ylab = "", axes = FALSE, ylim = c(0, 1000), xlim = c(0, 1000))
  
  # Le modele
  rect(xleft = 0, xright = 1000, ybottom = 50, ytop = 970, col = col_mod_bg, border = col_mod_bd)
  
  
  # --------------------------------------------------------------------------------
  # ENTREES / SORTIES
  # --------------------------------------------------------------------------------
  
  # Entrees P et ETP
  if (CemaNeige) {
    text(x = xy_P[1]*1.65, y = xy_P[2]*0.98, labels = "+ CemaNeige", adj = c(1, 1), font = 2, col = "grey40", cex = 1.6)
  }
  text(x = xy_P[1], y = xy_P[2], labels = "P", pos = 3, font = 2, col = col_P, cex = 1.8)
  text(x = xy_E[1], y = xy_E[2], labels = "E", pos = 3, font = 2, col = col_E, cex = 1.8)
  
  # Sorties Q
  text(x = xy_Q[1], y = xy_Q[2], labels = "Q", pos = 1, font = 2, col = col_Q, cex = 1.8)
  
  # Parametres
  tmp_decal   <- 20
  
  # --------------------------------------------------------------------------------
  # NEUTRALISATION DE P
  # --------------------------------------------------------------------------------
  
  # Interception
  segments(x0 = xy_E[1]-50, x1 = xy_P[1]+50,
           y0 = y_interception+tmp_decal, y1 = y_interception+tmp_decal)
  text(x = xy_P[1]+50, y = y_interception+20, labels = "Interception", pos = 4, font = 1, cex = 1.4)
  
  # P vers Pn
  segments(x0 = xy_P[1], x1 = xy_P[1], y0 = xy_P[2], y1 = y_interception+tmp_decal)
  
  # Pn vers Ps
  segments(x0 = xy_P[1], x1 = xy_P[1],
           y0 = y_interception, y1 = y_rendement+2*tmp_decal)
  segments(x0 = x_Ps, x1 = xy_P[1],
           y0 = y_rendement+2*tmp_decal, y1 = y_rendement+2*tmp_decal)
  segments(x0 = x_Ps, x1 = x_Ps,
           y0 = y_rendement+2*tmp_decal, y1 = y_rendement)
  
  # Pn vers Pn - Ps
  segments(x0 = xy_P[1], x1 = x_PnPs,
           y0 = y_rendement+2*tmp_decal, y1 = y_rendement+2*tmp_decal)
  segments(x0 = x_PnPs , x1 = x_PnPs,
           y0 = y_rendement+2*tmp_decal, y1 = y_rendement)
  
  # Pn - Ps vers Pr
  segments(x0 = x_PnPs, x1 = x_PnPs,
           y0 = y_rendement, y1 = y_percolation)
  
  # E vers En puis Es
  segments(x0 = xy_E[1], x1 = xy_E[1],
           y0 = xy_E[2], y1 = y_interception+tmp_decal)
  segments(x0 = xy_E[1], x1 = xy_E[1],
           y0 = y_interception, y1 = y_rendement)
  
  # Ecriture
  plotrix::boxed.labels(x = xy_P[1], y = y_interception, labels = "Pn",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  plotrix::boxed.labels(x = xy_E[1], y = y_interception, labels = "En",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  
  # ETP
  if (OutputsModel$PotEvap[i_pdt] != 0) {
    points(x = xy_PE[1], y =  xy_PE[2],
           type = "p", pch = tri_T, col = col_E,
           cex = cex_tri(OutputsModel$PotEvap[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
  # Precipitation
  if (OutputsModel$Precip[i_pdt] != 0) {
    points(x = xy_Precip[1], y =  xy_Precip[2],
           type = "p", pch = tri_B, col = col_P,
           cex = cex_tri(OutputsModel$Precip[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
  # Pn et Ps
  
  points(x = x_Ps, y = y_rendement+1.2*tmp_decal,
         type = "p", pch = tri_B, col = col_P,
         cex = cex_tri(OutputsModel$Ps[i_pdt], fact = fact_triangle, max = cex_max_poly))
  
  points(x = x_PnPs, y = y_rendement+1.2*tmp_decal,
         type = "p", pch = tri_B, col = col_P,
         cex = cex_tri(OutputsModel$Pn[i_pdt] - OutputsModel$Ps[i_pdt], fact = fact_triangle, max = cex_max_poly))
  
  
  # --------------------------------------------------------------------------------
  # FONCTION DE RENDEMENT
  # --------------------------------------------------------------------------------
  
  # Es
  plotrix::boxed.labels(x = xy_E[1], y = y_rendement, labels = "Es",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  
  # Evaporation reelle
  if (OutputsModel$AE[i_pdt] != 0) {
    points(x = xy_AE[1], y =  xy_AE[2],
           type = "p", pch = tri_T, col = col_P,
           cex = cex_tri(OutputsModel$AE[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
  # Ps et Pn - Ps
  plotrix::boxed.labels(x = x_Ps  , y = y_rendement, labels = "Ps"   ,
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  plotrix::boxed.labels(x = x_PnPs, y = y_rendement, labels = "Pn - Ps",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  
  # Reservoir de production
  rect(xleft   = xy_min_PROD[1], xright = xy_min_PROD[1]+base_res,
       ybottom = xy_min_PROD[2], ytop   = xy_min_PROD[2]+OutputsModel$Prod[i_pdt]*fact_res,
       col = col_SP, border = NA)
  segments(x0 = xy_min_PROD[1], x1 = xy_min_PROD[1]+base_res,
           y0 = xy_min_PROD[2], y1 = xy_min_PROD[2])
  segments(x0 = xy_min_PROD[1], x1 = xy_min_PROD[1],
           y0 = xy_min_PROD[2], y1 = xy_min_PROD[2]+Param[1]*fact_res)
  segments(x0 = xy_min_PROD[1]+base_res, x1 = xy_min_PROD[1]+base_res,
           y0 = xy_min_PROD[2], y1 = xy_min_PROD[2]+Param[1]*fact_res)
  text(x = 30, y = xy_min_PROD[2]+15, labels = "Prod.\nstore", cex = 1.4, pos = 4)
  
  
  # --------------------------------------------------------------------------------
  # PERCOLATION
  # --------------------------------------------------------------------------------
  
  # Reservoir de production vers Pr
  segments(x0 = xy_min_PROD[1]+base_res/2, x1 = xy_min_PROD[1]+base_res/2,
           y0 = xy_min_PROD[2], y1 = y_percolation)
  segments(x0 = xy_min_PROD[1]+base_res/2, x1 = x_PnPs,
           y0 = y_percolation, y1 = y_percolation)
  
  # Perc
  plotrix::boxed.labels(x = xy_min_PROD[1]+base_res/2, y = y_percolation, labels = "Perc.",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  
  # Valeur de Perc
  if (OutputsModel$Perc[i_pdt] != 0) {
    points(x = xy_min_PROD[1]+base_res+75, y = y_percolation,
           type = "p", pch = tri_R, col = col_P,
           cex = cex_tri(OutputsModel$Perc[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
  # parametres
  tmp_decal   <- (y_percolation - y_entreeUH) / 2 
  
  # Pr vers UH
  segments(x0 = x_PnPs, x1 = x_PnPs,
           y0 = y_percolation, y1 = y_entreeUH+tmp_decal/2)
  
  if (HydroModel %in% c("GR4J", "GR6J")) {
    
    # --------------------------------------------------------------------------------
    # SEPARATION DE PR
    # --------------------------------------------------------------------------------
    
    # Pr vers UH1
    segments(x0 = xy_Q9[1], x1 = x_PnPs,
             y0 = y_entreeUH+tmp_decal/2, y1 = y_entreeUH+tmp_decal/2)
    segments(x0 = xy_Q9[1], x1 = xy_Q9[1],
             y0 = y_entreeUH+tmp_decal/2, y1 = xy_Q9[2])
    
    # Pr vers UH2
    segments(x0 = x_PnPs, x1 = xy_Q1[1],
             y0 = y_entreeUH+tmp_decal/2, y1 = y_entreeUH+tmp_decal/2)
    segments(x0 = xy_Q1[1], x1 = xy_Q1[1],
             y0 = y_entreeUH+tmp_decal/2, y1 = y_routage)
    
    # Pr
    plotrix::boxed.labels(x = x_PnPs, y = y_percolation, labels = "Pr",
                          bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
    
    # Pr
    if (OutputsModel$PR[i_pdt] != 0) {
      points(x = x_PnPs[1], y = y_entreeUH+tmp_decal,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(OutputsModel$PR[i_pdt], fact = fact_triangle, max = cex_max_poly))
    }
    
    # --------------------------------------------------------------------------------
    # HYDROGRAMME UNITAIRE 1
    # --------------------------------------------------------------------------------
    
    # Entree de UH1
    if (OutputsModel$PR[i_pdt] != 0) {
      points(x = xy_Q9[1], y =y_entreeUH,
             type = "p", pch = tri_B, col = col_P,
             cex =  cex_tri(OutputsModel$PR[i_pdt]*0.9, fact = fact_triangle, max = cex_max_poly))
    }
    
    # Remplissage de UH1
    PR_mat_UH1_lg <- ceiling(Param[4])
    PR_mat_UH1_id <- max(i_pdt-PR_mat_UH1_lg+1, 1):i_pdt
    PR_mat_UH1 <- matrix(rep(c(rep(0, times = PR_mat_UH1_lg-length(PR_mat_UH1_id)+1),
                               OutputsModel$PR[PR_mat_UH1_id]), times = PR_mat_UH1_lg),
                         ncol = PR_mat_UH1_lg+1)[, -1L]
    PR_mat_UH1[lower.tri(PR_mat_UH1)] <- 0
    
    # --------------------------------------------------------------------------------
    # HYDROGRAMME UNITAIRE 2
    # --------------------------------------------------------------------------------
    
    # Entree de UH2
    if (OutputsModel$PR[i_pdt] != 0) {
      points(x = xy_Q1[1], y = y_entreeUH,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(OutputsModel$PR[i_pdt]*0.1, fact = fact_triangle, max = cex_max_poly))
    }
    
    
    # Remplissage de UH2
    PR_mat_UH2_lg <- ceiling(Param[4]*2)
    PR_mat_UH2_id <- max(i_pdt-PR_mat_UH2_lg+1, 1):i_pdt
    PR_mat_UH2 <- matrix(rep(c(rep(0, times = PR_mat_UH2_lg-length(PR_mat_UH2_id)+1),
                               OutputsModel$PR[PR_mat_UH2_id]), times = PR_mat_UH2_lg),
                         ncol = PR_mat_UH2_lg+1)[, -1L]
    PR_mat_UH2[lower.tri(PR_mat_UH2)] <- 0
    
  } 
  
  if (HydroModel == "GR5J") {
    
    # --------------------------------------------------------------------------------
    # SEPARATION DE PR
    # --------------------------------------------------------------------------------
    
    # sortie UH
    segments(x0 = x_PnPs, x1 = x_PnPs,
             y0 = y_entreeUH-2*tmp_decal, y1 = y_entreeUH-3*tmp_decal)
    
    # sortie UH vers branche 1
    segments(x0 = xy_Q9[1], x1 = x_PnPs,
             y0 = y_entreeUH-3*tmp_decal, y1 = y_entreeUH-3*tmp_decal)
    segments(x0 = xy_Q9[1], x1 = xy_Q9[1],
             y0 = y_entreeUH-3*tmp_decal, y1 = xy_Q9[2])
    
    # sortie UH vers branche 2
    segments(x0 = x_PnPs, x1 = xy_Q1[1],
             y0 = y_entreeUH-3*tmp_decal, y1 = y_entreeUH-3*tmp_decal)
    segments(x0 = xy_Q1[1], x1 = xy_Q1[1],
             y0 = y_entreeUH-3*tmp_decal, y1 = y_routage)
    
    # Pr
    plotrix::boxed.labels(x = x_PnPs, y = y_percolation, labels = "Pr",
                          bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
    
    
    # --------------------------------------------------------------------------------
    # HYDROGRAMME UNITAIRE
    # --------------------------------------------------------------------------------
    
    # Entree de UH (PR)
    if (OutputsModel$PR[i_pdt] != 0) {
      points(x = x_PnPs[1], y = y_entreeUH+tmp_decal,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(OutputsModel$PR[i_pdt], fact = fact_triangle, max = cex_max_poly))
    }
    
    
    # Remplissage de UH2
    PR_mat_UH2_lg <- ceiling(Param[4]*2)
    PR_mat_UH2_id <- max(i_pdt-PR_mat_UH2_lg+1, 1):i_pdt
    PR_mat_UH2 <- matrix(rep(c(rep(0, times = PR_mat_UH2_lg-length(PR_mat_UH2_id)+1),
                               OutputsModel$PR[PR_mat_UH2_id]), times = PR_mat_UH2_lg),
                         ncol = PR_mat_UH2_lg+1)[, -1L]
    PR_mat_UH2[lower.tri(PR_mat_UH2)] <- 0
    
    
    # Sorties de UH
    if (PR_mat_UH2[1] != 0) {
      points(x = x_PnPs[1], y =  y_entreeUH-5*tmp_decal/2,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(PR_mat_UH2[1], fact = fact_triangle, max = cex_max_poly))
    }
    
  }   
  
  # sortie de UH1 vers reservoirs exponentiel et de routage
  if (HydroModel == "GR6J") {
    segments(x0 = xy_Q9[1], x1 = xy_Q9[1],
             y0 = y_entreeUH-3*tmp_decal, y1 = xy_Q9[2])
    segments(x0 = xy_Q9[1]*0.80, x1 = xy_Q9[1]*1.30,
             y0 = xy_Q9[2], y1 = xy_Q9[2])
    segments(x0 = xy_Q9[1]*1.30, x1 = xy_Q9[1]*1.30,
             y0 = xy_Q9[2], y1 = xy_Q9[2]*0.65)    
    segments(x0 = xy_Q9[1]*0.80, x1 = xy_Q9[1]*0.80,
             y0 = xy_Q9[2], y1 = xy_Q9[2]*0.90)
    segments(x0 = xy_Q9[1]*0.55, x1 = xy_Q9[1]*0.55,
             y0 = xy_Q9[2]*0.70, y1 = y_routage)
    segments(x0 = xy_Q9[1]*0.55, x1 = xy_min_ROUT[1]+base_res/2,
             y0 = y_routage, y1 = y_routage)
  }
  
  # Q9
  if (OutputsModel$Q9[i_pdt] != 0) {
    points(x = xy_Q9[1], y = xy_Q9[2]+tmp_decal,
           type = "p", pch = tri_B, col = col_P,
           cex = cex_tri(OutputsModel$Q9[i_pdt], fact = fact_triangle, max = cex_max_poly))
    if (HydroModel == "GR6J") {
      # Q9 exp
      points(x = xy_Q9[1]*0.80, y = xy_Q9[1]*0.73,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(OutputsModel$Q9[i_pdt]*0.4, fact = fact_triangle, max = cex_max_poly))
      # Q9 rout      
      points(x = xy_Q9[1]*1.30, y = xy_Q9[1]*0.73,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(OutputsModel$Q9[i_pdt]*0.6, fact = fact_triangle, max = cex_max_poly))
      # QrExp
      plotrix::boxed.labels(x = xy_Q9[1]*0.55, y = y_routage, labels = "QrExp", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
    }
  }
  
  plotrix::boxed.labels(x = xy_Q9[1], y = xy_Q9[2], labels = "Q9",
                        bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  
  # Q1
  if (OutputsModel$Q1[i_pdt] != 0) {
    points(x = xy_Q1[1], y =  xy_Q1[2]+tmp_decal,
           type = "p", pch = tri_B, col = col_P,
           cex = cex_tri(OutputsModel$Q1[i_pdt], fact = fact_triangle, max = cex_max_poly))
    segments(x0 = xy_Q[1], x1 = xy_Q1[1], y0 = y_routage, y1 = y_routage)
  }
  
  plotrix::boxed.labels(x = xy_Q1[1], y = xy_Q1[2], labels = "Q1", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  
  # Valeur de Qd
  if (OutputsModel$QD[i_pdt] != 0) {
    points(x = xy_Q[1]+30, y =  y_routage,
           type = "p", pch = tri_L, col = col_P,
           cex = cex_tri(OutputsModel$QD[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
  # Qd
  plotrix::boxed.labels(x = xy_Q1[1], y = y_routage, labels = "Qd", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  
  
  # --------------------------------------------------------------------------------
  # RESERVOIR DE ROUTAGE
  # --------------------------------------------------------------------------------
  
  # Triche pour la taille du reservoire de routage
  tmp_triche   <- 0#80
  
  # Reservoir de routage
  rect(xleft = xy_min_ROUT[1], xright = xy_min_ROUT[1]+base_res,
       ybottom = xy_min_ROUT[2], ytop = xy_min_ROUT[2]+OutputsModel$Rout[i_pdt]*fact_res+tmp_triche,
       col = col_SR, border = NA)
  segments(x0 = xy_min_ROUT[1], x1 = xy_min_ROUT[1]+base_res, 
           y0 = xy_min_ROUT[2], y1 = xy_min_ROUT[2])
  segments(x0 = xy_min_ROUT[1], x1 = xy_min_ROUT[1],
           y0 = xy_min_ROUT[2], y1 = xy_min_ROUT[2]+Param[3]*fact_res+tmp_triche)
  segments(x0 = xy_min_ROUT[1]+base_res, x1 = xy_min_ROUT[1]+base_res,
           y0 = xy_min_ROUT[2], y1 = xy_min_ROUT[2]+Param[3]*fact_res+tmp_triche)
  text(x = 30, y = xy_min_ROUT[2]+15, labels = "Routing\nstore", cex = 1.4, pos = 4)
  
  # Sorties du reservoir
  segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = xy_min_ROUT[1]+base_res/2,
           y0 = xy_min_ROUT[2], y1 = y_routage)
  segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = xy_Q[1],
           y0 = y_routage, y1 = y_routage)
  
  # Qr
  if (HydroModel != "GR6J") {
    plotrix::boxed.labels(x = xy_min_ROUT[1]+base_res/2, y = y_routage, labels = "Qr",
                          bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  }
  if (HydroModel == "GR6J") {
    plotrix::boxed.labels(x = xy_min_ROUT[1]+base_res/1.5, y = (xy_min_ROUT[2]+y_routage)/2, labels = "Qr",
                          bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
  }
  
  # Valeur de Qr
  if (OutputsModel$QR[i_pdt] != 0) {
    if (HydroModel != "GR6J") {
    points(x = xy_Q[1]-100, y = y_routage,
           type = "p", pch = tri_R, col = col_P,
           cex = cex_tri(OutputsModel$QR[i_pdt], fact = fact_triangle, max = cex_max_poly))
    } else {
      points(x = xy_min_ROUT[1]+base_res/2, y = (xy_min_ROUT[2]+y_routage)/2,
             type = "p", pch = tri_B, col = col_P,
             cex = cex_tri(OutputsModel$QR[i_pdt], fact = fact_triangle, max = cex_max_poly))
    }
  }
  

  
  # --------------------------------------------------------------------------------
  # RESERVOIR EXPONENTIEL
  # --------------------------------------------------------------------------------

  if (HydroModel == "GR6J") {
    
    # Triche pour la taille du reservoire exponentiel
    tmp_triche   <- 0#80
    
    # Exp en log
    signExp <- ifelse(OutputsModel$Exp[i_pdt] > 0, +1, -1)
    logExpIpdt <- log(abs(OutputsModel$Exp[i_pdt])+1e-6) * signExp
    logExpMax  <- log(max(abs(OutputsModel$Exp  ))+1e-6)

    # Reservoir exponentiel
    rect(xleft = xy_min_EXPO[1], xright = xy_min_EXPO[1]+base_res,
         ybottom = xy_min_EXPO[2], ytop = xy_min_EXPO[2]+logExpIpdt*fact_resExp+tmp_triche,
         col = ifelse(OutputsModel$Exp[i_pdt] > 0, "#10B510", "#FF0303"), border = NA)
    segments(x0 = xy_min_EXPO[1], x1 = xy_min_EXPO[1]+base_res, 
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2])
    segments(x0 = xy_min_EXPO[1], x1 = xy_min_EXPO[1],
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]+logExpMax*fact_resExp+tmp_triche)
    segments(x0 = xy_min_EXPO[1]+base_res, x1 = xy_min_EXPO[1]+base_res,
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]+logExpMax*fact_resExp+tmp_triche)
    segments(x0 = xy_min_EXPO[1], x1 = xy_min_EXPO[1],
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]-logExpMax*fact_resExp-tmp_triche)
    segments(x0 = xy_min_EXPO[1]+base_res, x1 = xy_min_EXPO[1]+base_res,
             y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2]-logExpMax*fact_resExp-tmp_triche)
    text(x = 30, y = xy_min_EXPO[2]+00, labels = "Exp.\nstore", cex = 1.4, pos = 4)
    points(x = 180, y = xy_min_EXPO[2]+20, pch = 43, # +
         cex = 2.0, col = "#10B510")
    points(x = 178, y = xy_min_EXPO[2]-20, pch = 95, # -
         cex = 1.6, col = "#FF0303")

    # Valeur de QrExp
    if (OutputsModel$QR[i_pdt] != 0) {
      points(x = xy_Q[1]-350, y = y_routage,
             type = "p", pch = tri_R, col = col_P,
             cex = cex_tri(OutputsModel$QRExp[i_pdt], fact = fact_triangle, max = cex_max_poly))
    }

    # Valeur de Qr + QrExp
    if (OutputsModel$QR[i_pdt] != 0) {    
      points(x = xy_Q[1]-100, y = y_routage,
             type = "p", pch = tri_R, col = col_P,
             cex = cex_tri(OutputsModel$QR[i_pdt]+OutputsModel$QRExp[i_pdt], fact = fact_triangle, max = cex_max_poly))
    }
  }
  
  
  # --------------------------------------------------------------------------------
  # Q FINAL
  # -------------------------------------------------------------------------------
  
  # Q final
  segments(x0 = xy_Q[1], x1 = xy_Q[1], y0 = y_routage, y1 = xy_Q[2]+10)
  if (OutputsModel$Qsim[i_pdt] != 0) {
    points(x = xy_Q[1], y =  mean(c(y_routage, xy_Q[2]+10)),
           type = "p", pch = tri_B, col = col_Q,
           cex = cex_tri(OutputsModel$Qsim[i_pdt], fact = fact_triangle, max = cex_max_poly))
  }
  
  
  # --------------------------------------------------------------------------------
  # EXCHANGE
  # --------------------------------------------------------------------------------
  
  # Actual exchange Q9
  arrows(x0 = xy_min_ROUT[1]+base_res, x1 = 1025,
         y0 = y_Ech_Q9               , y1 = y_Ech_Q9,
         length = 0.075, angle = 20, code = 3)
  pch <- ifelse(OutputsModel$AExch1[i_pdt] < 0, tri_R, tri_L)
  points(x = xy_min_ROUT[1]+base_res+130, y = y_Ech_Q9,
         type = "p", pch = pch, col = col_P,
         cex = cex_tri(OutputsModel$AExch1[i_pdt], fact = fact_triangle, max = cex_max_poly))
  
  # Actual exchange Q1
  arrows(x0 = xy_Q1[1], x1 = 1025,
         y0 = y_Ech_Q1, y1 = y_Ech_Q1,
         length = 0.075, angle = 20, code = 3)
  pch <- ifelse(OutputsModel$AExch2[i_pdt] < 0, tri_R, tri_L)
  points(x = xy_Q1[1]+100, y = y_Ech_Q1,
         type = "p", pch = pch, col = col_P,
         cex = cex_tri(OutputsModel$AExch2[i_pdt], fact = fact_triangle, max = cex_max_poly)) 
  
  # Actual exchange Q9 exp.
  if (HydroModel%in% c("GR6J")) {
    arrows(x0 = xy_min_EXPO[1]+base_res[1], x1 = 1025,
           y0 = xy_min_EXPO[2], y1 = xy_min_EXPO[2],
           length = 0.075, angle = 20, code = 3)
    pch <- ifelse(OutputsModel$Exch[i_pdt] < 0, tri_R, tri_L)
    points(x = xy_Q1[1]+100, y = xy_min_EXPO[2],
           type = "p", pch = pch, col = col_P,
           cex = cex_tri(OutputsModel$Exch[i_pdt], fact = fact_triangle, max = cex_max_poly)) 
  }
  
  if (HydroModel%in% c("GR4J", "GR6J")) {
    
    # --------------------------------------------------------------------------------
    # UH 1 & 2 PLOT
    # --------------------------------------------------------------------------------
    
    tmpUH1 <- PR_mat_UH1 %*% UH1[seq_len(PR_mat_UH1_lg)] * 0.75 # 0.75 au lieu de 0.90 pour reduire dif. visuelle
    tmpUH2 <- PR_mat_UH2 %*% UH2[seq_len(PR_mat_UH2_lg)] * 0.25 # 0.25 au lieu de 0.10 pour reduire dif. visuelle
    palUH0 <- colorRampPalette(c("#428BCA", "#B5D2EA"))
    par(mar = rep(1, 4))
    par(new = TRUE, plt = c(0.35, 0.46,  0.40, 0.45), yaxt = "n", bg = "red")
    rect(xleft = 325, xright = 475,  ybottom = 390, ytop = 445, col = col_mod_bg, border = NA)
    par(new = TRUE)
    barplot(tmpUH1, beside = TRUE, space = 0, col = palUH0(length(tmpUH1)),
            ylim = range(c(0, 5/length(PR_mat_UH2_lg)), na.rm = TRUE))
    par(new = TRUE, plt = c(0.67, 0.89,  0.40, 0.45), yaxt = "n")
    rect(xleft = 0, xright = 900,  ybottom = 0, ytop = 470, col = col_mod_bg, border = NA)
    par(new = TRUE)
    barplot(tmpUH2, beside = TRUE, space = 0, col = palUH0(length(tmpUH2)),
            ylim = range(c(0, 5/length(PR_mat_UH2_lg)), na.rm = TRUE))
    
  }
  
  
  if (HydroModel == "GR5J") {
    
    # --------------------------------------------------------------------------------
    # UH PLOT
    # --------------------------------------------------------------------------------
    
    tmpUH2 <- PR_mat_UH2 %*% UH2[seq_len(PR_mat_UH2_lg)] * 1 / 2 # * 1/2 car seul HU et pour que ca rentre
    palUH0 <- colorRampPalette(c("#428BCA", "#B5D2EA"))
    par(mar = rep(1, 4))
    par(new = TRUE, plt = c(0.57, 0.79,  0.45, 0.5), yaxt = "n")
    par(new = TRUE)
    barplot(tmpUH2, beside = TRUE, space = 0, col = palUH0(length(tmpUH2)),
            ylim = range(c(0, 5/length(PR_mat_UH2_lg)), na.rm = TRUE))
    
  }
  
}