diff --git a/DESCRIPTION b/DESCRIPTION
index 33733c865c4072f85e88dcafe82b7c51468ab725..bafe11af9e50da32ee1cb51e12c67a053bbc6fca 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: airGRteaching
 Type: Package
 Title: Tools to Simplify the Use of the airGR Hydrological Package for Education (Including a Shiny Application)
-Version: 0.1.2.4
+Version: 0.1.2.5
 Date: 2017-04-04
 Authors@R: c(person("Olivier", "Delaigue", role = c("aut", "cre"), email = "airGR@irstea.fr"), person("Laurent", "Coron", role = c("aut")), person("Pierre", "Brigode", role = c("aut")), person("Guillaume", "Thirel", role = c("ctb")))
 Depends: airGR (>= 1.0.5.22)
diff --git a/R/Utils.R b/R/Utils.R
index 9a875b2199449bbc15c43ab10d1a187bbc2fe370..e68d77e061f8cba03475579086ebf9cae459e217 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -27,10 +27,13 @@ if (getRversion() >= "2.15.1") {
   StrName <- "(.*)(GR)(\\d{1})(\\D{1})"
   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)
   NbParam  <- gsub(StrName, "\\3", x)
   isCN     <- grepl("CemaNeige"  , x)
   
-  res <- list(TypeModel = x, NbParam = as.numeric(NbParam), TimeUnit = TimeUnit, CemaNeige = isCN)
+  res <- list(TypeModel = x, NbParam = as.numeric(NbParam),
+              TimeUnit = TimeUnit, TimeLag = TimeLag,
+              CemaNeige = isCN)
   return(res)
 }
 
@@ -60,14 +63,14 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   nom_modele      <- "GR4J"
   xy_E            <- c(250, 970)
   xy_P            <- c(600, 970)
-  xy_Q            <- c(750, 30)
+  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)
+  xy_Q9           <- c(400, 310)
+  xy_Q1           <- c(800, 310)
   y_routage       <- 100
   fact_res        <- 200/800/3
   base_res        <- 300
@@ -75,20 +78,24 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   xy_min_PROD     <- c(200, 610)
   xy_min_ROUT     <- c(250, 150)
   y_entreeHU      <- 500
-  xy_HU1          <- c(500, 420)
-  xy_HU2          <- c(900, 420)
-  y_Ech_Q1        <- 200
-  y_Ech_Q9        <- 180
+  xy_UH1          <- c(500, 420)
+  xy_UH2          <- c(900, 420)
+  y_Ech_Q1        <- 170 #200
+  y_Ech_Q9        <- 150 #180
   y_out_aubes     <- 400
   D               <- 5/2
   xpad            <- 1.5
   ypad            <- 1.5
   max_triangle    <- max(unlist(OutputsModel[c("Perc", "PR", "Q9", "Q1", "QR", "QD")]))
   fact_var        <- 40
-  fact_triangle   <- 25
+  fact_triangle   <- 100#25
+  cex_tri         <- function(cex, fact = 25, max) suppressWarnings(log(cex * fact + 1) / max)
   radius1         <- 0
   radius2         <- 60
-  cex_max_poly    <- 0.005
+  cex_max_poly    <- 2 # 0.005
+  tri_R           <- -0x25BA
+  tri_B           <- -0x25BC
+  tri_L           <- -0x25C4
   
   par(col.axis = par("fg"))
   
@@ -109,20 +116,20 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
     cex <- log(cex * fact_triangle + 1)
     
     if (dir == "D") {
-      xy[,1] <- c(x_center-cex/2, x_center-cex/2, x_center+cex/2)
-      xy[,2] <- c(y_center-cex/2, y_center+cex/2, y_center)
+      xy[, 1] <- c(x_center-cex/2, x_center-cex/2, x_center+cex/2)
+      xy[, 2] <- c(y_center-cex/2, y_center+cex/2, y_center)
     }
     if (dir == "G") {
-      xy[,1] <- c(x_center+cex/2, x_center+cex/2, x_center-cex/2)
-      xy[,2] <- c(y_center-cex/2, y_center+cex/2, y_center)
+      xy[, 1] <- c(x_center+cex/2, x_center+cex/2, x_center-cex/2)
+      xy[, 2] <- c(y_center-cex/2, y_center+cex/2, y_center)
     }
     if (dir == "H") {
-      xy[,1] <- c(x_center-cex/2, x_center+cex/2, x_center)
-      xy[,2] <- c(y_center-cex/2, y_center-cex/2, y_center+cex/2)
+      xy[, 1] <- c(x_center-cex/2, x_center+cex/2, x_center)
+      xy[, 2] <- c(y_center-cex/2, y_center-cex/2, y_center+cex/2)
     }
     if (dir == "B") {
-      xy[,1] <- c(x_center-cex/2, x_center+cex/2, x_center)
-      xy[,2] <- c(y_center+cex/2, y_center+cex/2, y_center-cex/2)
+      xy[, 1] <- c(x_center-cex/2, x_center+cex/2, x_center)
+      xy[, 2] <- c(y_center+cex/2, y_center+cex/2, y_center-cex/2)
     }
     
     return(xy) 
@@ -140,7 +147,7 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   # HUs
   ##################################################################################
   
-  # Calcul des ordonnees SH1 de l' "hydrogramme unitaire cumule" HU1
+  # 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
@@ -148,7 +155,7 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
     if (i >= Param[4])           SH1[i] <- 1
   }
   
-  # Calcul des ordonnees UH1 de l' "hydrogramme unitaire discret" HU1
+  # Calcul des ordonnees UH1 de l' "hydrogramme unitaire discret" UH1
   UH1     <- array(NA, NH)
   for (j in 1:NH) {
     if (j == 1) {
@@ -158,7 +165,7 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
     }  
   }
   
-  # Calcul des ordonnees SH2 de l' "hydrogramme unitaire cumule" HU2
+  # 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
@@ -167,7 +174,7 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
     if (i >= 2*Param[4])                  SH2[i] <- 1
   }
   
-  # Calcul des ordonnees UH2 de l' "hydrogramme unitaire discret" HU2
+  # Calcul des ordonnees UH2 de l' "hydrogramme unitaire discret" UH2
   UH2     <- array(NA, 2*NH)
   for (j in 1:(2*NH)) {
     if (j == 1) {
@@ -303,31 +310,42 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   tmp_decal   <- 20
   
   # 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, col = par("fg"), cex = 1.4)
+  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 Ps
-  segments(x0 = x_PnPs, x1 = x_PnPs, y0 = y_rendement, y1 = y_percolation)
+  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 Ps
+  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)
+  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
-  boxed.labels(x = xy_P[1], y = y_interception, labels = "Pn", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
-  boxed.labels(x = xy_E[1], y = y_interception, labels = "En", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  boxed.labels(x = xy_P[1], y = y_interception, labels = "Pn",
+               bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  boxed.labels(x = xy_E[1], y = y_interception, labels = "En",
+               bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
   
   
   ##################################################################################
@@ -335,20 +353,27 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   ##################################################################################
   
   # Es
-  boxed.labels(x = xy_E[1], y = y_rendement, labels = "Es", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  boxed.labels(x = xy_E[1], y = y_rendement, labels = "Es",
+               bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
   
-  # Ps et Pn-Ps
-  boxed.labels(x = x_Ps  , y = y_rendement, labels = "Ps"   , col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
-  boxed.labels(x = x_PnPs, y = y_rendement, labels = "Pn-Ps", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  # Ps et Pn - Ps
+  boxed.labels(x = x_Ps  , y = y_rendement, labels = "Ps"   ,
+               bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  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_R, 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], col = par("fg"))
-  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, col = par("fg"))
-  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, col = par("fg"))
-  text(x = 30, y = xy_min_PROD[2]+Param[1]*fact_res/3, "Prod. store", cex = 1.4, col = par("fg"), pos = 4)
+  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]+Param[1]*fact_res/3, labels = "Prod. store",
+       cex = 1.4, pos = 4)
   
   
   ##################################################################################
@@ -356,17 +381,23 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   ##################################################################################
   
   # 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)
+  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
-  boxed.labels(x = xy_min_PROD[1]+base_res/2, y = y_percolation, labels = "Perc.", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  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) {
     tmp_triangle  <- create_polygon(x_center = xy_min_PROD[1]+base_res+75, y_center = y_percolation, 
                                     cex = OutputsModel$Perc[i_pdt]*max_triangle / cex_max_poly, dir = "D")
     polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
+    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))
   }
   
   
@@ -377,34 +408,44 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   # parametres
   tmp_decal   <- (y_percolation - y_entreeHU) / 2 
   
-  # Pr vers HU1
-  segments(x0 = x_PnPs, x1 = x_PnPs, y0 = y_percolation, y1 = y_entreeHU+tmp_decal)
-  segments(x0 = xy_Q9[1], x1 = x_PnPs, y0 = y_entreeHU+tmp_decal, y1 = y_entreeHU+tmp_decal)
-  segments(x0 = xy_Q9[1], x1 = xy_Q9[1], y0 = y_entreeHU+tmp_decal, y1 = xy_Q9[2])
+  # Pr vers UH1
+  segments(x0 = x_PnPs, x1 = x_PnPs,
+           y0 = y_percolation, y1 = y_entreeHU+tmp_decal)
+  segments(x0 = xy_Q9[1], x1 = x_PnPs,
+           y0 = y_entreeHU+tmp_decal, y1 = y_entreeHU+tmp_decal)
+  segments(x0 = xy_Q9[1], x1 = xy_Q9[1],
+           y0 = y_entreeHU+tmp_decal, y1 = xy_Q9[2])
   
-  # Pr vers HU2
-  segments(x0 = x_PnPs, x1 = xy_Q1[1], y0 = y_entreeHU+tmp_decal, y1 = y_entreeHU+tmp_decal)
-  segments(x0 = xy_Q1[1], x1 = xy_Q1[1], y0 = y_entreeHU+tmp_decal, y1 = y_routage)
+  # Pr vers UH2
+  segments(x0 = x_PnPs, x1 = xy_Q1[1],
+           y0 = y_entreeHU+tmp_decal, y1 = y_entreeHU+tmp_decal)
+  segments(x0 = xy_Q1[1], x1 = xy_Q1[1],
+           y0 = y_entreeHU+tmp_decal, y1 = y_routage)
   
   # Pr
-  boxed.labels(x = x_PnPs, y = y_percolation, labels = "Pr", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  boxed.labels(x = x_PnPs, y = y_percolation, labels = "Pr",
+               bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
   
   
   ##################################################################################
   # HYDROGRAMME UNITAIRE 1 (AUBE)
   ##################################################################################
   
-  # Entree de HU1
+  # Entree de UH1
   if (OutputsModel$PR[i_pdt] != 0) {
     tmp_triangle  <- create_polygon(x_center = xy_Q9[1], y_center = y_entreeHU+tmp_decal/2, 
                                     cex = OutputsModel$PR[i_pdt]*max_triangle / cex_max_poly, dir = "B")
     polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
+    points(x = xy_Q9[1], y =y_entreeHU+tmp_decal/2,
+           type = "p", pch = tri_B, col = col_P,
+           cex =  cex_tri(OutputsModel$PR[i_pdt]*0.9, fact = fact_triangle, max = cex_max_poly))
   }
+
   
-  # Fleche vers HU1
+  # Fleche vers UH1
   # arrows(x0 = xy_Q9[1], y0 = y_entreeHU, x1 = xy_Q9[1]+30, y1 = y_entreeHU-10, length = 0.075, angle = 20, col = "red")
   
-  # Remplissage de HU1
+  # Remplissage de UH1
   # tmp_UH1      <- OutputsModel$PR[i_pdt]*0.9*fact_var*UH1
   # tmp_UH1[1]   <- 0
   # tmp_UH1      <- suppressWarnings(log(sqrt(OutputsModel$PR[i_pdt]*0.9*UH1) +1)) * 60 / max_UH1
@@ -427,8 +468,8 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   #   }
   #   angles     <- seq(angle1, angle2, by = 0.01)
   #   angles[length(angles)] <- angle2
-  #   xpos <- c(cos(angles) * radius1, cos(rev(angles)) * radius2) + xy_HU1[1]
-  #   ypos <- c(sin(angles) * radius1, sin(rev(angles)) * radius2) + xy_HU1[2]
+  #   xpos <- c(cos(angles) * radius1, cos(rev(angles)) * radius2) + xy_UH1[1]
+  #   ypos <- c(sin(angles) * radius1, sin(rev(angles)) * radius2) + xy_UH1[2]
   #   polygon(xpos, ypos, col = col_R, border = NA)
   # }
   
@@ -449,39 +490,48 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   #     radius2 <- temp
   #   }
   #   
-  #   segments(x0 = xy_HU1[1], y0 = xy_HU1[2], x1 = (cos(angle1)*radius2) + xy_HU1[1], y1 = (sin(angle1)*radius2) + xy_HU1[2])
-  #   segments(x0 = xy_HU1[1], y0 = xy_HU1[2], x1 = (cos(angle2)*radius2) + xy_HU1[1], y1 = (sin(angle2)*radius2) + xy_HU1[2])
+  #   segments(x0 = xy_UH1[1], y0 = xy_UH1[2], x1 = (cos(angle1)*radius2) + xy_UH1[1], y1 = (sin(angle1)*radius2) + xy_UH1[2])
+  #   segments(x0 = xy_UH1[1], y0 = xy_UH1[2], x1 = (cos(angle2)*radius2) + xy_UH1[1], y1 = (sin(angle2)*radius2) + xy_UH1[2])
   # }
   
-  # Fleche sortant de HU1
+  # Fleche sortant de UH1
   # arrows(x0 = xy_Q9[1]+30, y0 = y_out_aubes+10, x1 = xy_Q9[1], y1 = y_out_aubes, length = 0.075, angle = 20, col = "red")
   
-  # Sorties de HU1
-  # if (OutputsModel$Q9[i_pdt] != 0) {
-  #   tmp_triangle  <- create_polygon(x_center = xy_Q9[1], y_center = xy_Q9[2]+tmp_decal,
-  #                                   cex = OutputsModel$Q9[i_pdt]*max_triangle / cex_max_poly, dir = "B")
-  #   polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
-  # }
+  # Sorties de UH1
+  if (OutputsModel$Q9[i_pdt] != 0) {
+    tmp_triangle  <- create_polygon(x_center = xy_Q9[1], y_center = xy_Q9[2]+tmp_decal,
+                                    cex = OutputsModel$Q9[i_pdt]*max_triangle / cex_max_poly, dir = "B")
+    polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
+    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))
+  }
+
   
   # Q9
-  boxed.labels(x = xy_Q9[1], y = xy_Q9[2], labels = "Q9", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  boxed.labels(x = xy_Q9[1], y = xy_Q9[2], labels = "Q9",
+               bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
   
   
   ##################################################################################
   # HYDROGRAMME UNITAIRE 2
   ##################################################################################
   
-  # Entree de HU2
+  # Entree de UH2
   if (OutputsModel$PR[i_pdt] != 0) {
     tmp_triangle  <- create_polygon(x_center = xy_Q1[1], y_center = y_entreeHU+tmp_decal/2,
                                     cex = OutputsModel$PR[i_pdt]*0.1*max_triangle / cex_max_poly, dir = "B")
     polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
+    points(x = xy_Q1[1], y = y_entreeHU+tmp_decal/2,
+           type = "p", pch = tri_B, col = col_P,
+           cex = cex_tri(OutputsModel$PR[i_pdt]*0.1, fact = fact_triangle, max = cex_max_poly))
   }
+
   
-  # Fleche vers HU2
+  # Fleche vers UH2
   # arrows(x0 = xy_Q1[1], y0 = y_entreeHU, x1 = xy_Q1[1]+30, y1 = y_entreeHU-10, length = 0.075, angle = 20, col = "red")
   
-  # Remplissage de HU2
+  # Remplissage de UH2
   # tmp_UH2      <- OutputsModel$PR[i_pdt]*0.1*fact_var*UH2
   # tmp_UH2[1]   <- 0
   # tmp_UH1      <- suppressWarnings(log((sqrt(OutputsModel$PR[i_pdt]*0.1*UH2)) +1)) * 60 / max_UH2
@@ -504,8 +554,8 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   #   # }
   #   angles     <- seq(angle1, angle2, by = 0.01)
   #   angles[length(angles)] <- angle2
-  #   xpos <- c(cos(angles) * radius1, cos(rev(angles)) * radius2) + xy_HU2[1]
-  #   ypos <- c(sin(angles) * radius1, sin(rev(angles)) * radius2) + xy_HU2[2]
+  #   xpos <- c(cos(angles) * radius1, cos(rev(angles)) * radius2) + xy_UH2[1]
+  #   ypos <- c(sin(angles) * radius1, sin(rev(angles)) * radius2) + xy_UH2[2]
   #   polygon(xpos, ypos, col = col_R, border = NA)
   # }
   
@@ -526,37 +576,55 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   #     radius2 <- temp
   #   }
   #   
-  #   segments(x0 = xy_HU2[1], y0 = xy_HU2[2], x1 = (cos(angle1)*radius2) + xy_HU2[1], y1 = (sin(angle1)*radius2) + xy_HU2[2])
-  #   segments(x0 = xy_HU2[1], y0 = xy_HU2[2], x1 = (cos(angle2)*radius2) + xy_HU2[1], y1 = (sin(angle2)*radius2) + xy_HU2[2])
+  #   segments(x0 = xy_UH2[1], y0 = xy_UH2[2], x1 = (cos(angle1)*radius2) + xy_UH2[1], y1 = (sin(angle1)*radius2) + xy_UH2[2])
+  #   segments(x0 = xy_UH2[1], y0 = xy_UH2[2], x1 = (cos(angle2)*radius2) + xy_UH2[1], y1 = (sin(angle2)*radius2) + xy_UH2[2])
   # }
   
-  # Fleche sortant de HU2
+  # Fleche sortant de UH2
   # arrows(x0 = xy_Q1[1]+30, y0 = y_out_aubes+10, x1 = xy_Q1[1], y1 = y_out_aubes, length = 0.075, angle = 20, col = "red")
   
-  # Sorties de HU2
+  # Sorties de UH2
   if (OutputsModel$Q1[i_pdt] != 0) {
     tmp_triangle  <- create_polygon(x_center = xy_Q1[1], y_center = xy_Q1[2]+tmp_decal,
                                     cex = OutputsModel$Q1[i_pdt]*max_triangle / cex_max_poly, dir = "B")
     polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
+    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)
   }
-  segments(x0 = xy_Q[1], x1 = xy_Q1[1], y0 = y_routage, y1 = y_routage)
+
   
   # Q1
-  boxed.labels(x = xy_Q1[1], y = xy_Q1[2], labels = "Q1", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  boxed.labels(x = xy_Q1[1], y = xy_Q1[2], labels = "Q1", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
   
   # Echange
-  arrows(x0 = xy_Q1[1], y0 = y_Ech_Q1, x1 = 1025    , y1 = y_Ech_Q1, length = 0.075, angle = 20)
-  arrows(x0 = 1025    , y0 = y_Ech_Q1, x1 = xy_Q1[1], y1 = y_Ech_Q1, length = 0.075, angle = 20)
+  tmp_Exch    <- Param[2]*(OutputsModel$Rout[i_pdt]/Param[3])^(7/2)
+  if (tmp_Exch < 0) {
+    tmp_sens  <- "G"
+  } else {
+    tmp_sens  <- "D"
+  }
+  arrows(x0 = xy_Q1[1], x1 = 1025,
+         y0 = y_Ech_Q1, y1 = y_Ech_Q1,
+         length = 0.075, angle = 20, code = 3)
+  points(x = xy_Q1[1]+100, y = y_Ech_Q1,
+         type = "p", pch = ifelse(tmp_Exch < 0, tri_R, tri_L), col = col_P,
+         cex = cex_tri(tmp_Exch, fact = fact_triangle, max = cex_max_poly))  
   
   # Valeur de Qd
   if (OutputsModel$QD[i_pdt] != 0) {
     tmp_triangle  <- create_polygon(x_center = xy_Q[1]+10, y_center = y_routage,
                                     cex = OutputsModel$QD[i_pdt]*max_triangle / cex_max_poly, dir = "G")
     polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
+    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))
+    # print(cex_tri(OutputsModel$QD[i_pdt], fact = fact_triangle, max = cex_max_poly))
   }
-  
+
   # Qd
-  boxed.labels(x = xy_Q1[1], y = y_routage, labels = "Qd", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  boxed.labels(x = xy_Q1[1], y = y_routage, labels = "Qd", bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
   
   
   ##################################################################################
@@ -564,69 +632,98 @@ DiagramGR4J <- function(OutputsModel, Param, SimPer, EventDate) {
   ##################################################################################
   
   # Triche pour la taille du reservoire de routage
-  tmp_triche   <- 80
+  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_R, 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], col = par("fg"))
-  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, col = par("fg"))
-  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, col = par("fg"))
-  text(x = 50, y = xy_min_ROUT[2]+Param[1]*fact_res/3, "Routing store", cex = 1.4, col = par("fg"), pos = 4)
+  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_R, 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 = 50, y = xy_min_ROUT[2]+Param[3]*fact_res/3, labels = "Routing store",
+       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, col = par("fg"))
-  segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = xy_Q[1], y0 = y_routage, y1 = y_routage)
+  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)
   
   # Echange
-  arrows(x0 = xy_min_ROUT[1]+base_res, y0 = y_Ech_Q9, x1 = 1025,y1 = y_Ech_Q9, length = 0.075, angle = 20)
-  arrows(x0 = 1025, y0 = y_Ech_Q9, x1 = xy_min_ROUT[1]+base_res,y1 = y_Ech_Q9, length = 0.075, angle = 20)
+  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)
   
   # Qr
-  boxed.labels(x = xy_min_ROUT[1]+base_res/2, y = y_routage, labels = "Qr", col = par("fg"), bg = col_mod_bg, border = NA, xpad = xpad, ypad = ypad)
+  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)
   
   # Valeur de Qr
   if (OutputsModel$QR[i_pdt] != 0) {
     tmp_triangle  <- create_polygon(x_center = xy_Q[1]-50, y_center = y_routage,
                                     cex = OutputsModel$QR[i_pdt]*max_triangle / cex_max_poly, dir = "D")
     polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
+    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))
   }
   
   # Q final
   segments(x0 = xy_Q[1], x1 = xy_Q[1], y0 = y_routage, y1 = xy_Q[2]+10)
   
-  
-  
-  
-  
-  
   # Echange
-  tmp_Exch    <- Param[2]*(OutputsModel$Rout[i_pdt]/Param[3])^(7/2)
-  if (tmp_Exch < 0) {
-    tmp_sens  <- "G"
-  } else {
-    tmp_sens  <- "D"
-  }
-  # segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = -10, y0 = y_routage+40, y1 = y_routage+40, col = par("fg"))
+  # tmp_Exch    <- Param[2]*(OutputsModel$Rout[i_pdt]/Param[3])^(7/2)
+  # if (tmp_Exch < 0) {
+  #   tmp_sens  <- "G"
+  # } else {
+  #   tmp_sens  <- "D"
+  # }
+  # segments(x0 = xy_min_ROUT[1]+base_res/2, x1 = -10, y0 = y_routage+40, y1 = y_routage+40)
   tmp_triangle  <- create_polygon(x_center = xy_min_ROUT[1]-base_res/2, y_center = y_routage+40, cex = tmp_Exch*50, dir = tmp_sens)
   polygon(x = tmp_triangle[,1], y = tmp_triangle[,2], border = NA, col = col_R)
-  
-  
-  tmpHU1 <- OutputsModel$PR[i_pdt]*0.75*UH1 #*0.9*fact_var*UH1
-  tmpHU2 <- OutputsModel$PR[i_pdt]*0.25*UH2 #*0.1*fact_var*UH2
-  colHU1 <- rgb(066, 139, 202, maxColorValue = 255, alpha = seq(255, 100, length.out = length(tmpHU1)))
-  colHU2 <- rgb(066, 139, 202, maxColorValue = 255, alpha = seq(255, 100, length.out = length(tmpHU2)))
+  points(x = xy_min_ROUT[1]+base_res+100, y = y_Ech_Q9,
+         type = "p", pch = ifelse(tmp_Exch < 0, tri_R, tri_L), col = col_P,
+         cex = cex_tri(tmp_Exch, fact = fact_triangle, max = cex_max_poly))
 
+  max_UH1 <- max(max(UH1)*OutputsModel$PR*0.75)
+  max_UH2 <- max(max(UH2)*OutputsModel$PR*0.25)
+  
+  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
+  
+  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
+  
+  
+  # tmpUH1 <- OutputsModel$PR[i_pdt]*0.75*UH1 #*0.9*fact_var*UH1
+  # tmpUH2 <- OutputsModel$PR[i_pdt]*0.25*UH2 #*0.1*fact_var*UH2
+  tmpUH1 <- PR_mat_UH1 %*% UH1[seq_len(PR_mat_UH1_lg)] * 0.75 #*0.9*fact_var*UH1
+  tmpUH2 <- PR_mat_UH2 %*% UH2[seq_len(PR_mat_UH2_lg)] * 0.25 #*0.1*fact_var*UH2
+  # colUH1 <- rgb(066, 139, 202, maxColorValue = 255, alpha = seq(255, 100, length.out = length(tmpUH1)))
+  # colUH2 <- rgb(066, 139, 202, maxColorValue = 255, alpha = seq(255, 100, length.out = length(tmpUH2)))
+  palUH0 <- colorRampPalette(c("#428BCA", "#B5D2EA"))
+  colUH1 <- palUH0(length(tmpUH1))
+  colUH2 <- palUH0(length(tmpUH1))
+  
   par(mar = rep(1, 4))
-  par(new = TRUE, plt = c(0.35, 0.46,  0.40, 0.45), yaxt = "n") # c(0.43, 0.54,  0.42, 0.50)
+  par(new = TRUE, plt = c(0.35, 0.46,  0.40, 0.45), yaxt = "n", bg = "red") # c(0.43, 0.54,  0.42, 0.50)
   rect(xleft = 325, xright = 475,  ybottom = 390, ytop = 445, col = col_mod_bg, border = NA) #col_mod_bg
   par(new = TRUE)
-  barplot(tmpHU1, beside = TRUE, space = 0, col = colHU1, ylim = range(c(0, tmpHU1, tmpHU2), na.rm = TRUE))
+  #barplot(tmpUH1, beside = TRUE, space = 0, col = colUH1, ylim = range(c(0, tmpUH1, tmpUH2), na.rm = TRUE))
+  barplot(tmpUH1, beside = TRUE, space = 0, col = colUH1, ylim = range(c(0, 5/length(PR_mat_UH2_lg)), na.rm = TRUE))
   par(new = TRUE, plt = c(0.78-0.11, 0.89,  0.40, 0.45), yaxt = "n") # c(0.83, 0.94,  0.42, 0.50)
   rect(xleft = 0, xright = 900,  ybottom = 0, ytop = 470, col = col_mod_bg, border = NA)
   par(new = TRUE)
-  barplot(tmpHU2, beside = TRUE, space = 0, col = colHU2, ylim = range(c(0, tmpHU1, tmpHU2), na.rm = TRUE))
-
+  #barplot(tmpUH2, beside = TRUE, space = 0, col = colUH2, ylim = range(c(0, tmpUH1, tmpUH2), na.rm = TRUE))
+  barplot(tmpUH2, beside = TRUE, space = 0, col = colUH2, ylim = range(c(0, 5/length(PR_mat_UH2_lg)), na.rm = TRUE))
   
 }