diff --git a/DESCRIPTION b/DESCRIPTION
index 8b8cf3daeaad4245b8ac80daea3fa546452dbac2..2b7b80351894a393d99e9c10101dc2ad63d1600e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.3.2.16
-Date: 2019-06-04
+Version: 1.3.2.17
+Date: 2019-06-12
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
   person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@irstea.fr"),
diff --git a/NEWS.rmd b/NEWS.rmd
index ec3dfb91d2b38baddc24a775a6aadece5bb7daa7..a4bd9a5fd48fcb0bf05ef84b97f8808e94bf784e 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -14,7 +14,7 @@ output:
 
 
 
-### 1.3.2.16 Release Notes (2019-06-04)
+### 1.3.2.17 Release Notes (2019-06-12)
 
 
 #### New features
@@ -30,6 +30,12 @@ output:
 
 - The <code>PEdaily_Oudin()</code> function is deprecated and his use has been replaced by the use of <code>PE_Oudin()</code>.
 
+- <code>plot.OutputsModel()</code> now presents a <code>LayoutMat</code> argument (and additionnal related argument: <code>LayoutWidths</code>, <code>LayoutHeights</code>) to specify complex plot arrangements.
+
+#### Bug fixes
+
+- Fixed bug in <code>plot.OutputsModel()</code>. The function now runs correctly when the <code>which</code> argument contains the <code>"CorQQ"</code> value without <code>"CumFreq"</code>.
+
 
 #### Major user-visible changes
 
@@ -40,8 +46,12 @@ output:
 
 #### Minor user-visible changes
 
+- <code>.ErrorCrit()</code> private function added to check inputs into <code>ErrorCrit_&#42;()</code> functions. The <code>ErrorCrit_&#42;()</code> functions were simplified accordingly.
+
 - <code>CreateInputsCrit()</code> now returns <code>FUN_CRIT</code> as a character string.
 
+- An example is addeed to illustred the use of the <code>plot.OutputsModel()</code> function.
+
 ____________________________________________________________________________________
 
 
diff --git a/R/plot.OutputsModel.R b/R/plot.OutputsModel.R
index d97feb100ab4ecc4abd726d2bbb82eb0c4358a86..deca5dc582c0d3b5491c10bb10021581941a9f3a 100644
--- a/R/plot.OutputsModel.R
+++ b/R/plot.OutputsModel.R
@@ -1,20 +1,31 @@
 plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea = NULL, which = "synth", log_scale = FALSE,
-                              cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1, verbose = TRUE, ...) {
+                              cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1,
+                              LayoutMat = NULL, LayoutWidths = rep.int(1, ncol(LayoutMat)), LayoutHeights = rep.int(1, nrow(LayoutMat)),
+                              verbose = TRUE, ...) {
+  
+  
+  ## save default graphical parameters and resetting on exit
+  opar <- par(no.readonly = TRUE)
+  on.exit(par(opar))
+  
   
   OutputsModel <- x
   
+  
+  
+  ## ---------- check arguments
+  
   if (!inherits(OutputsModel, "GR") & !inherits(OutputsModel, "CemaNeige")) {
     stop("'OutputsModel' not in the correct format for default plotting")
   }
   
-  
+  ## check 'OutputsModel'
   BOOL_Dates <- FALSE
-  
   if ("DatesR" %in% names(OutputsModel)) {
     BOOL_Dates <- TRUE
   }
-  BOOL_Pobs <- FALSE
   
+  BOOL_Pobs <- FALSE
   if ("Precip" %in% names(OutputsModel)) {
     BOOL_Pobs <- TRUE
   }
@@ -28,6 +39,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   if ("Qsim"   %in% names(OutputsModel)) {
     BOOL_Qsim <- TRUE
   }
+  
   BOOL_Qobs <- FALSE
   if (BOOL_Qsim & length(Qobs) == length(OutputsModel$Qsim)) {
     if (sum(is.na(Qobs)) != length(Qobs)) {
@@ -57,51 +69,72 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   }
   
   
+  ## check 'which'
+  whichNeedQobs  <- c("Error", "CorQQ")
+  whichDashboard <- c("all", "synth", "ts", "perf")
+  whichAll   <- c("Precip", "PotEvap", "Temp", "SnowPack", "Flows", "Error", "Regime", "CumFreq", "CorQQ")
+  whichSynth <- c("Precip"           , "Temp", "SnowPack", "Flows"         , "Regime", "CumFreq", "CorQQ")
+  whichTS    <- c("Precip", "PotEvap", "Temp", "SnowPack", "Flows"                                       )
+  whichPerf  <- c(                                                  "Error", "Regime", "CumFreq", "CorQQ")
+  whichCN    <- c(                     "Temp", "SnowPack"                                                )
+  warnMsgWhich   <- "'which' must be a vector of character"
+  warnMsgNoQobs  <- "the %s plot(s) cannot be drawn if there is no 'Qobs'"
+  warnMsgWhichCN <- sprintf("incorrect element found in argument 'which':\n\twithout CemaNeige, %s are not available \n\tit can only contain %s",
+                            paste0(shQuote(whichCN), collapse = " and "),
+                            paste0(shQuote(c(whichDashboard, whichAll[!whichAll %in% whichCN])), collapse = ", "))
   if (is.null(which)) {
-    stop("'which' must be a vector of character")
+    stop(warnMsgWhich)
   }
   if (!is.vector(which)) {
-    stop("'which' must be a vector of character")
+    stop(warnMsgWhich)
   }
   if (!is.character(which)) {
-    stop("'which' must be a vector of character")
+    stop(warnMsgWhich)
   }
-  if (any(!which %in% c("all", "synth", "ts", "perf", "PotEvap", "Precip", 'Temp', "SnowPack", "Flows", "Error", "Regime", "CumFreq", "CorQQ"))) {
-    stop("incorrect element found in argument 'which':\nit can only contain 'all', 'synth', 'ts', 'perf', 'Precip', 'PotEvap', 'Temp', 'SnowPack', 'Error', 'Flows', 'Regime', 'CumFreq' or 'CorQQ'")
+  if (any(!which %in% c(whichDashboard, whichAll))) {
+    stop(sprintf("incorrect element found in argument 'which': %s\nit can only contain %s",
+                 paste0(shQuote(which[!which %in% c(whichDashboard, whichAll)])),
+                 paste0(shQuote(c(whichDashboard, whichAll)), collapse = ", ")))
   }
-  if (all(which %in% c("Temp", "SnowPack")) & !inherits(OutputsModel, "CemaNeige")) {
-    stop("Incorrect element found in argument 'which':\nwithout CemaNeige it can only contain 'all', 'synth', 'ts', 'perf', 'Precip', 'PotEvap', 'Flows', 'Error', 'Regime', 'CumFreq' or 'CorQQ'")
+  if (all(which %in% whichCN) & !inherits(OutputsModel, "CemaNeige")) {
+    stop(warnMsgWhichCN)
   }
-  if (length(unique(which %in% c("Temp", "SnowPack"))) == 2 & !inherits(OutputsModel, "CemaNeige")) {
-    warning("Incorrect element found in argument 'which':\nit can only contain 'all', 'synth', 'ts', 'perf', 'Precip', 'PotEvap', 'Flows', 'Error', 'Regime', 'CumFreq' or 'CorQQ'\nwithout CemaNeige 'Temp' and 'SnowPack' are not available")
-  }  
-  
-  
-  if ("all" %in% which) {
-    which <- c("Precip", "PotEvap", "Temp", "SnowPack", "Flows", "Error", "Regime", "CumFreq", "CorQQ")
+  if (length(unique(which %in% whichCN)) == 2 & !inherits(OutputsModel, "CemaNeige")) {
+    warning(warnMsgWhichCN)
   }
-  if ("synth" %in% which) {
-    which <- c("Precip", "Temp", "SnowPack", "Flows", "Regime", "CumFreq", "CorQQ")
+  if (all(!which %in% c("all", "synth", "ts", whichCN)) & !inherits(OutputsModel, "GR")) {
+    stop(sprintf("incorrect element found in argument 'which': \nwith CemaNeige alone, only %s are available",
+                 paste0(shQuote(c("all", "synth", "ts", "Temp", "SnowPack")), collapse = ", ")))
   }
-  if ("ts" %in% which) {
-    which <- c("Precip", "PotEvap", "Temp", "SnowPack", "Flows")
+  if (any(!which %in% c("all", "synth", "ts", whichCN)) & !inherits(OutputsModel, "GR")) {
+    warning(sprintf("incorrect element found in argument 'which': \nwith CemaNeige alone, only %s are available",
+                    paste0(shQuote(c("all", "synth", "ts", "Temp", "SnowPack")), collapse = ", ")))
   }
   if ("perf" %in% which) {
-    which <- c("Error", "Regime", "CumFreq", "CorQQ")
+    which <- c(which, whichPerf)
   }
-  
-  
-  if (is.null(Qobs)) {
-    if (length(which) == 1 & any(which %in% "Error")) {
-      stop("the 'Error' time srie cannot be draw if there is no 'Qobs'")
+  if ("ts" %in% which) {
+    which <- c(which, whichTS)
+  }
+  if ("synth" %in% which) {
+    which <- c(which, whichSynth)
+  }
+  if ("all" %in% which) {
+    which <- c(which, whichAll)
+  }
+  if (is.null(Qobs) & inherits(OutputsModel, "GR")) {
+    if (length(which) == 1 & (any(which %in% whichNeedQobs))) {
+      stop(sprintf(warnMsgNoQobs, shQuote(which)))
     }
-    if (length(which) != 1 & any(which %in% c("Error", "all"))) {
+    if (length(which) != 1 & any(which %in% whichNeedQobs)) {
+      BOOL_CorQQ <- FALSE
       BOOL_Error <- FALSE
-      warning("the 'Error' time serie cannot be draw if there is no 'Qobs'")
+      warning(sprintf(warnMsgNoQobs, paste0(shQuote(whichNeedQobs), collapse = " and ")))
     }
   }
+
   
-  
+  ## check dates
   if (!BOOL_Dates) {
     stop("'OutputsModel' must contain at least 'DatesR' to allow plotting")
   }
@@ -182,6 +215,15 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     SelectQsimNotZero <- round(OutputsModel$Qsim[IndPeriod_Plot], 4) != 0
     BOOL_QsimZero <- sum(!SelectQsimNotZero, na.rm = TRUE) > 0
   }
+  if ( BOOL_Qobs & !BOOL_Qsim) {
+    SelectNotZero <- SelectQobsNotZero
+  }
+  if (!BOOL_Qobs &  BOOL_Qsim) {
+    SelectNotZero <- SelectQsimNotZero
+  }
+  if ( BOOL_Qobs &  BOOL_Qsim) {
+    SelectNotZero <- SelectQobsNotZero & SelectQsimNotZero
+  }
   if (BOOL_QobsZero & verbose) {
     warning("zeroes detected in 'Qobs': some plots in the log space will not be created using all time-steps")
   }
@@ -190,86 +232,79 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   }
   BOOL_FilterZero <- TRUE
   
-  ## Plots_choices
-  BOOLPLOT_Precip   <- "Precip"   %in% which & BOOL_Pobs
-  BOOLPLOT_PotEvap  <- "PotEvap"  %in% which & BOOL_Eobs
-  BOOLPLOT_Temp     <- "Temp"     %in% which & BOOL_Snow
-  BOOLPLOT_SnowPack <- "SnowPack" %in% which & BOOL_Snow
+  
+  
+  ## ---------- plot
+  
+  ## plot choices
+  BOOLPLOT_Precip   <- "Precip"   %in% which &  BOOL_Pobs
+  BOOLPLOT_PotEvap  <- "PotEvap"  %in% which &  BOOL_Eobs
+  BOOLPLOT_Temp     <- "Temp"     %in% which &  BOOL_Snow
+  BOOLPLOT_SnowPack <- "SnowPack" %in% which &  BOOL_Snow
   BOOLPLOT_Flows    <- "Flows"    %in% which & (BOOL_Qsim | BOOL_Qobs)
-  BOOLPLOT_Error    <- "Error"    %in% which & BOOL_Error
-  BOOLPLOT_Regime   <- "Regime"   %in% which & BOOL_TS & BOOL_Qsim & (NameTS %in% c("hour", "day", "month"))
+  BOOLPLOT_Error    <- "Error"    %in% which &  BOOL_Error
+  BOOLPLOT_Regime   <- "Regime"   %in% which &  BOOL_Qsim & BOOL_TS & (NameTS %in% c("hour", "day", "month"))
   BOOLPLOT_CumFreq  <- "CumFreq"  %in% which & (BOOL_Qsim | BOOL_Qobs) & BOOL_FilterZero
   BOOLPLOT_CorQQ    <- "CorQQ"    %in% which & (BOOL_Qsim & BOOL_Qobs) & BOOL_FilterZero
   
   
-  ## Options
+  ## options
   BLOC <- TRUE
   if (BLOC) {
     lwdk <- 1.8
     line <- 2.6
     bg   <- NA
     
-    matlayout <- NULL
-    iPlot <- 0
-    
-    Sum1 <- sum(c(BOOLPLOT_Precip, BOOLPLOT_SnowPack, BOOLPLOT_Flows))
-    Sum2 <- sum(c(BOOLPLOT_Regime, BOOLPLOT_CumFreq, BOOLPLOT_CorQQ))
-    if (BOOLPLOT_Precip) {
-      matlayout <- rbind(matlayout, c(iPlot + 1, iPlot + 1, iPlot + 1))
-      iPlot <- iPlot + 1
-    }
-    if (BOOLPLOT_PotEvap) {
-      matlayout <- rbind(matlayout, c(iPlot + 1, iPlot + 1, iPlot + 1), c(iPlot + 1, iPlot + 1, iPlot + 1))
-      iPlot <- iPlot + 1
-    } 
-    if (BOOLPLOT_Temp) {
-      matlayout <- rbind(matlayout, c(iPlot + 1, iPlot + 1, iPlot + 1), c(iPlot + 1, iPlot + 1, iPlot + 1))
-      iPlot <- iPlot + 1
-    }      
-    if (BOOLPLOT_SnowPack) {
-      matlayout <- rbind(matlayout, c(iPlot + 1, iPlot + 1, iPlot + 1), c(iPlot + 1, iPlot + 1, iPlot + 1))
-      iPlot <- iPlot + 1
-    }
-    if (BOOLPLOT_Flows) {
-      matlayout <- rbind(matlayout, c(iPlot + 1, iPlot + 1, iPlot + 1), c(iPlot + 1, iPlot + 1, iPlot + 1))
-      iPlot <- iPlot + 1
-    }
-    if (BOOLPLOT_Error) {
-      matlayout <- rbind(matlayout, c(iPlot + 1, iPlot + 1, iPlot + 1), c(iPlot + 1, iPlot + 1, iPlot + 1))
-      iPlot <- iPlot + 1
-    }    
-    if ((Sum1 >= 1 & Sum2 != 0) | (Sum1 == 0 & Sum2 == 3)) {
-      matlayout <- rbind(matlayout, c(iPlot + 1, iPlot + 2, iPlot + 3), c(iPlot + 1, iPlot + 2, iPlot + 3))
-      iPlot <- iPlot + 3
-    }
-    if (Sum1 == 0 & Sum2 == 2) {
-      matlayout <- rbind(matlayout, c(iPlot + 1, iPlot + 2))
-      iPlot <- iPlot + 2
+    ## Set plot arrangement
+    if (is.null(LayoutMat)) {     
+      matlayout <- NULL
+      iPlot <- 0
+      iHght <- NULL
+      
+      listBOOLPLOT1 <- c(Precip = BOOLPLOT_Precip, PotEvap  = BOOLPLOT_PotEvap,
+                         Temp   = BOOLPLOT_Temp  , SnowPack = BOOLPLOT_SnowPack,
+                         Flows  = BOOLPLOT_Flows , Error    = BOOLPLOT_Error)
+      listBOOLPLOT2 <- c(Regime = BOOLPLOT_Regime, CumFreq  = BOOLPLOT_CumFreq,
+                         CorQQ = BOOLPLOT_CorQQ)
+      Sum1 <- sum(listBOOLPLOT1)
+      Sum2 <- sum(listBOOLPLOT2)
+      
+      for (k in seq_len(Sum1)) {
+        matlayout <- rbind(matlayout, iPlot + c(1, 1, 1))
+        iPlot <- iPlot + 1
+        iHght <- c(iHght, 0.7)
+      }
+      ## Flows plot is higher than the other TS
+      listBOOLPLOT1 <- listBOOLPLOT1[listBOOLPLOT1]
+      listBOOLPLOTF <- (names(listBOOLPLOT1) == "Flows") * BOOLPLOT_Flows
+      iHght <- iHght + listBOOLPLOTF * listBOOLPLOT1 * 0.3
+      if (Sum2 >= 1) {
+        iHght <- c(iHght, 1.0)
+      }
+      if ((Sum1 >= 1 & Sum2 != 0) | (Sum1 == 0 & Sum2 == 3)) {
+        matlayout <- rbind(matlayout, iPlot + c(1, 2, 3))
+        iPlot <- iPlot + 3
+      }
+      if (Sum1 == 0 & Sum2 == 2) {
+        matlayout <- rbind(matlayout, iPlot + c(1, 2))
+        iPlot <- iPlot + 2
+      }
+      if (Sum1 == 0 & Sum2 == 1) {
+        matlayout <- rbind(matlayout, iPlot + 1)
+        iPlot <- iPlot + 1
+      }
+
+      iPlotMax <- iPlot
+      LayoutWidths  <- rep.int(1, ncol(matlayout))
+      LayoutHeights <- iHght #rep.int(1, nrow(matlayout))
     }
-    if (Sum1 == 0 & Sum2 == 1) {
-      matlayout <- rbind(matlayout, iPlot + 1)
-      iPlot <- iPlot + 1
+    if (!is.null(LayoutMat)) {
+      matlayout <- LayoutMat
     }
-    iPlotMax <- iPlot
-    
-    # isRStudio <- Sys.getenv("RSTUDIO") == "1";
-    # if (!isRStudio) {
-    #   if (Sum1 == 1 & Sum2 == 0) {width = 10; height = 05;}
-    #   if (Sum1 == 1 & Sum2 != 0) {width = 10; height = 07;}
-    #   if (Sum1 == 2 & Sum2 == 0) {width = 10; height = 05;}
-    #   if (Sum1 == 2 & Sum2 != 0) {width = 10; height = 07;}
-    #   if (Sum1 == 3 & Sum2 == 0) {width = 10; height = 07;}
-    #   if (Sum1 == 3 & Sum2 != 0) {width = 10; height = 10;}
-    #   if (Sum1 == 0 & Sum2 == 1) {width = 05; height = 05;}
-    #   if (Sum1 == 0 & Sum2 == 2) {width = 10; height = 04;}
-    #   if (Sum1 == 0 & Sum2 == 3) {width = 10; height = 03;}
-    #   dev.new(width = width, height = height)
-    #}
+    layout(matlayout, widths = LayoutWidths, heights = LayoutHeights)
     
-    opar <- par(no.readonly = TRUE)
-    on.exit(par(opar))
+
     
-    layout(matlayout)
     
     Xaxis <- 1:length(IndPeriod_Plot)
     if (BOOL_Dates) {
@@ -409,7 +444,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cex.axis, ...)
     }
   }
-
+  
   
   ## SnowPack
   if (BOOLPLOT_SnowPack) {
@@ -540,23 +575,34 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
     
-    errorQ <- OutputsModel$Qsim / Qobs
+    if (log_scale) {
+      errorQ <- log(OutputsModel$Qsim[IndPeriod_Plot]) - log(Qobs[IndPeriod_Plot])
+    } else {
+      errorQ <- OutputsModel$Qsim[IndPeriod_Plot] - Qobs[IndPeriod_Plot]
+    }
     par(new = FALSE, mar = mar)
-    ylim1 <- range(errorQ[IndPeriod_Plot], na.rm = TRUE)
-    
-    plot(Xaxis, errorQ[IndPeriod_Plot],
+    ylim1 <- range(errorQ[SelectNotZero], na.rm = TRUE)
+    plot(Xaxis, errorQ,
          type = "l", xaxt = "n", yaxt = "n", ylim = ylim1,
-         col = "grey50", lwd = lwd * lwdk,
-         xlab = "", ylab = "", log = ifelse(log_scale, "y", ""),
-         panel.first = abline(h = 1, col = "grey", lty = 2), ...)
+         col = par("fg"), lwd = lwd * lwdk,
+         xlab = "", ylab = "",
+         panel.first = abline(h = 0, col = "royalblue"), ...)
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    mtext(side = 2, paste("flow err.", plotunit), cex = cex.lab, line = line)
+    mtext(side = 2, paste("flow error", plotunit), cex = cex.lab, line = line)
+    if (!is.null(BasinArea)) {
+      Factor <- Factor_UNIT_M3S
+      axis(side = 4, at = pretty(ylim1*Factor)/Factor, labels = pretty(ylim1*Factor), cex.axis = cex.axis, ...)
+      mtext(side = 4, paste("flow error", "[m3/s]"), line = line, cex = cex.lab)
+    }
     if (BOOL_Dates) {
       axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cex.axis, ...)
       axis(side = 1, at = Seq2, labels = Labels2, lwd.ticks = 1.5, cex.axis = cex.axis, ...)
     } else {
       axis(side = 1, at = pretty(Xaxis), labels = pretty(Xaxis), cex.axis = cex.axis, ...)
     }
+    if (log_scale) {
+      legend("bottomright", "log scale", lty = 1, col = NA, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
+    }
   }
   
   
@@ -574,9 +620,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
       mtext(side = 1, text = "", line = line, cex = cex.lab)
       text(0, 0, labels = "NO ENOUGH VALUES", col = "grey40")
-      txtlab <- "flow regime"
+      txtlab <- "flow"
       if (BOOL_Pobs) {
-        txtlab <- "precip. & flow regime"
+        txtlab <- "precip. & flow"
       }
       mtext(side = 2, paste(txtlab, plotunitregime), line = line, cex = cex.lab)
     } else {
@@ -715,16 +761,16 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       axis(side = 2, at = pretty(ylimQ), labels = pretty(ylimQ), cex.axis = cex.axis, ...)
       mtext(side = 1, labX, line = line, cex = cex.lab)
       posleg <- "topright"
-      txtlab <- "flow regime"
+      txtlab <- "flow"
       if (BOOL_Pobs) {
         posleg <- "right"
-        txtlab <- "precip. & flow regime"
+        txtlab <- "precip. & flow"
       }
       mtext(side = 2, paste(txtlab, plotunitregime), line = line, cex = cex.lab)
       if (!is.null(BasinArea)) {
         Factor <- Factor_UNIT_M3S / (365.25 / 12)
         axis(side = 4, at = pretty(ylimQ*Factor)/Factor, labels = pretty(ylimQ*Factor), cex.axis = cex.axis, ...)
-        mtext(side = 4, paste("flow regime", "[m3/s]"), line = line, cex = cex.lab)
+        mtext(side = 4, paste("flow", "[m3/s]"), line = line, cex = cex.lab)
       }
       ### posleg <- "topright"; if (BOOL_Pobs) {posleg <- "right";}
       ### legend(posleg, txtleg, col = colleg, lty = 1, lwd = lwdleg, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
@@ -733,7 +779,6 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   }
   
   
-  
   ## Cumulative_frequency
   if (BOOLPLOT_CumFreq) {
     kPlot <- kPlot + 1
@@ -741,15 +786,15 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     par(new = FALSE, mar = mar)
     xlim <- c(0, 1)
     if ( BOOL_Qobs & !BOOL_Qsim) {
-      SelectNotZero <- SelectQobsNotZero
+      # SelectNotZero <- SelectQobsNotZero
       ylim <- range(log(Qobs[IndPeriod_Plot][SelectNotZero]), na.rm = TRUE)
     }
     if (!BOOL_Qobs &  BOOL_Qsim) {
-      SelectNotZero <- SelectQsimNotZero
+      # SelectNotZero <- SelectQsimNotZero
       ylim <- range(log(OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero]), na.rm = TRUE)
     }
     if ( BOOL_Qobs &  BOOL_Qsim) {
-      SelectNotZero <- SelectQobsNotZero & SelectQsimNotZero
+      # SelectNotZero <- SelectQobsNotZero & SelectQsimNotZero
       ylim <- range(log(c(Qobs[IndPeriod_Plot][SelectNotZero], OutputsModel$Qsim[IndPeriod_Plot][SelectNotZero])), na.rm = TRUE)
     }
     SelectNotZero <- ifelse(is.na(SelectNotZero), FALSE, SelectNotZero)
@@ -839,13 +884,15 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   }
   
   ## Empty_plots
-  while (kPlot < iPlotMax) {
-    kPlot <- kPlot + 1
-    par(new = FALSE)
-    plot(0, 0, type = "n", xlab = "", ylab = "", axes = FALSE, ...)
+  if (exists("iPlotMax")) {
+    while (kPlot < iPlotMax) {
+      kPlot <- kPlot + 1
+      par(new = FALSE)
+      plot(0, 0, type = "n", xlab = "", ylab = "", axes = FALSE, ...)
+    }
   }
   
   ## Restoring_layout_options
-  layout(1)
+  # layout(1)
   
 }
diff --git a/data/exampleSimPlot.rda b/data/exampleSimPlot.rda
new file mode 100644
index 0000000000000000000000000000000000000000..fd930f46a4e74046d5cecc7f5c556a8176aff2a0
Binary files /dev/null and b/data/exampleSimPlot.rda differ
diff --git a/man/plot.OutputsModel.Rd b/man/plot.OutputsModel.Rd
index c2fd5e27c255b8a29fcd7c3c02575e804dad52a6..bf88ca091f4aaed07f688fa37ce71621f9a578ea 100644
--- a/man/plot.OutputsModel.Rd
+++ b/man/plot.OutputsModel.Rd
@@ -4,6 +4,9 @@
 \name{plot.OutputsModel}
 \alias{plot.OutputsModel}
 \alias{plot}
+\alias{exampleSimPlot}
+\alias{simGR4J}
+\alias{simCNGR4J}
 
 
 \title{Default preview of model outputs}
@@ -12,12 +15,16 @@
 \usage{
 \method{plot}{OutputsModel}(x, Qobs = NULL, IndPeriod_Plot = NULL,
   BasinArea = NULL, which = "synth", log_scale = FALSE,
-  cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1, verbose = TRUE, ...)
+  cex.axis = 1, cex.lab = 0.9, cex.leg = 0.9, lwd = 1,
+  LayoutMat = NULL,
+  LayoutWidths = rep.int(1, ncol(LayoutMat)),
+  LayoutHeights = rep.int(1, nrow(LayoutMat)),
+  verbose = TRUE, ...)
 }
 
 
 \arguments{
-\item{x}{[object of class \emph{OutputsModel}] list of model outputs (which must at least include DatesR, Precip and Qsim) [POSIXlt, mm, mm]}
+\item{x}{[object of class \emph{OutputsModel}] list of model outputs (which must at least include DatesR, Precip and Qsim) [POSIXlt, mm/time step, mm/time step]}
 
 \item{Qobs}{(optional) [numeric] time series of observed flow (for the same time steps than simulated) [mm/time step]}
 
@@ -25,9 +32,9 @@
 
 \item{BasinArea}{(optional) [numeric] basin area [km2], used to plot flow axes in m3/s}
 
-\item{which}{(optional) [character] choice of plots \cr (e.g. c(\code{"Precip"}, \code{"Temp"}, \code{"SnowPack"}, \code{"Flows"}, \code{"Regime"}, \code{"CumFreq"}, \code{"CorQQ"})), default = \code{"all"}}
+\item{which}{(optional) [character] choice of plots \cr (e.g. c(\code{"Precip"}, \code{"Temp"}, \code{"SnowPack"}, \code{"Flows"}, \code{"Regime"}, \code{"CumFreq"}, \code{"CorQQ"})), default = \code{"synth"}, see details below}
 
-\item{log_scale}{(optional) [boolean] indicating if the flow axis is to be logarithmic, default = \code{FALSE}}
+\item{log_scale}{(optional) [boolean] indicating if the flow time series axis and the flow error time series axis are to be logarithmic, default = \code{FALSE}}
 
 \item{cex.axis}{(optional) [numeric] the magnification to be used for axis annotation relative to the current setting of \code{cex}}
 
@@ -37,6 +44,12 @@
 
 \item{lwd}{(optional) [numeric] the line width (a positive number)}
 
+\item{LayoutMat}{(optional) [numeric] a matrix object specifying the location of the next N figures on the output device. Each value in the matrix must be 0 or a positive integer. If N is the largest positive integer in the matrix, then the integers {1, \dots, N-1} must also appear at least once in the matrix (see \code{\link{layout}})}
+
+\item{LayoutWidths}{(optional) [numeric] a vector of values for the widths of columns on the device (see \code{\link{layout}})}
+
+\item{LayoutHeights}{(optional) [numeric] a vector of values for the heights of rows on the device (see \code{\link{layout}})}
+
 \item{verbose}{(optional) [boolean] indicating if the function is run in verbose mode or not, default = \code{TRUE}}
 
 \item{...}{other parameters to be passed through to plotting functions}
@@ -54,16 +67,28 @@ Function which creates a screen plot giving an overview of the model outputs.
 
 
 \details{
-Dashboard of results including various graphs (depending on the model):\cr
-  (1) time series of total precipitation\cr
-  (2) time series of temperature (plotted only if CemaNeige is used)\cr
-  (3) time series of snow pack (plotted only if CemaNeige is used)\cr
-  (4) time series of simulated flows (and observed flows if provided)\cr
-  (5) interannual median monthly simulated flow (and observed flows if provided)\cr
-  (6) correlation plot between simulated and observed flows (if observed flows provided)\cr
-  (7) cumulative frequency plot for simulated flows (and observed flows if provided)
+Different types of independent graphs are available (depending on the model, but always drawn in this order):
+\itemize{
+  \item \code{"Precip"}: time series of total precipitation
+  \item \code{"PotEvap"}: time series of potential evapotranspiration
+  \item \code{"Temp"}: time series of temperature (plotted only if CemaNeige is used)
+  \item \code{"SnowPack"}: time series of snow water equivalent (plotted only if CemaNeige is used)
+  \item \code{"Flows"}: time series of simulated flows (and observed flows if provided)
+  \item \code{"Regime"}: interannual median monthly simulated flow (and observed flows if provided)
+  \item \code{"CorQQ"}: correlation plot between simulated and observed flows (only if observed flows provided)
+  \item \code{"CumFreq"}: cumulative frequency plot for simulated flows (and observed flows if provided)
 }
 
+Different dashboards of results including various graphs are available:
+\itemize{
+  \item \code{"perf"}: corresponds to \code{"Error"}, \code{"Regime"}, \code{"CumFreq"} and \code{"CorQQ"}
+  \item \code{"ts"}: corresponds to \code{"Precip"}, \code{"PotEvap"}, \code{"Temp"}, \code{"SnowPack"} and \code{"Flows"}
+  \item \code{"synth"}: corresponds to \code{"Precip"}, \code{"Temp"}, \code{"SnowPack"}, \code{"Flows"}, \code{"Regime"}, \code{"CumFreq"} and \code{"CorQQ"}
+  \item \code{"all"}: corresponds to \code{"Precip"}, \code{"PotEvap"}, \code{"Temp"}, \code{"SnowPack"}, \code{"Flows"}, \code{"Error"}, \code{"Regime"}, \code{"CumFreq"} and \code{"CorQQ"}
+  }
+  
+If several dashboards are selected, or if an independent graph is called with a dashboard, the graphical device will include all the requested graphs without redundancy.
+}
 
 \author{
 Laurent Coron, Olivier Delaigue, Guillaume Thirel
@@ -71,5 +96,51 @@ Laurent Coron, Olivier Delaigue, Guillaume Thirel
 
 
 \examples{
-### See examples of RunModel_GR4J or RunModel_CemaNeigeGR4J functions
+### see examples of RunModel_GR4J or RunModel_CemaNeigeGR4J functions
+### to understand how the example datasets have been prepared
+
+
+## loading examples dataset for GR4J and GR4J + CemaNeige
+data(exampleSimPlot)
+
+
+### Qobs and outputs from GR4J and GR4J + CemaNeige models
+str(simGR4J, max.level = 1)
+str(simCNGR4J, max.level = 1)
+
+
+### default dashboard (which = "synth")
+
+## GR models whithout CemaNeige
+plot(simGR4J$OutputsModel, Qobs = simGR4J$Qobs)
+
+## GR models whith CemaNeige ("Temp" and "SnowPack" added)
+plot(simCNGR4J$OutputsModel, Qobs = simCNGR4J$Qobs)
+
+
+### "Error" and "CorQQ" plots cannot be display whithout Qobs
+plot(simGR4J$OutputsModel, which = "all", Qobs = simGR4J$Qobs)
+plot(simGR4J$OutputsModel, which = "all", Qobs = NULL)
+
+
+### complex plot arrangements
+plot(simGR4J$OutputsModel, Qobs = simGR4J$Qobs,
+     which = c("Flows", "Regime", "CumFreq", "CorQQ"),
+     LayoutMat = matrix(c(1, 2, 3, 1, 4, 4), ncol = 2),
+     LayoutWidths  = c(1.5, 1),
+     LayoutHeights = c(0.5, 1, 1))
+    
+    
+### add a main title
+
+## the whole list of settable par's
+opar <- par(no.readonly = TRUE)
+
+## define outer margins and a title inside it
+par(oma = c(0, 0, 3, 0))
+plot(simGR4J$OutputsModel, Qobs = simGR4J$Qobs)
+title(main = "GR4J outputs", outer = TRUE, line = 1.2, cex.main = 1.4)
+
+## reset original par
+par(opar)
 }