diff --git a/DESCRIPTION b/DESCRIPTION
index 42fc032d6b981b1a06a293314f7fa6dc45dd2b8f..5ba39d1c43ea91f7f7422273545e9575d588a203 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.2.15.7
-Date: 2019-05-03
+Version: 1.2.16.0
+Date: 2019-05-20
 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 7a0eb0ac57a29570e889222c53cb8c7c4d1c2f24..f96e4a973c6a17a9f4983bb1f300289bfea3b60b 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -14,13 +14,15 @@ output:
 
 
 
-### 1.2.15.7 Release Notes (2019-05-03)
+### 1.2.16.0 Release Notes (2019-05-20)
 
 
 #### New features
 
 - <code>CreateInputsCrit()</code> now allows power (as numeric or character values) and the Box-Cox transformations in the <code>transfo</code> argument.
 
+- <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.
+
 
 #### Major user-visible changes
 
diff --git a/R/plot.OutputsModel.R b/R/plot.OutputsModel.R
index d97feb100ab4ecc4abd726d2bbb82eb0c4358a86..9ac287b145b0c53e9b77e8c137859ebfc45c303b 100644
--- a/R/plot.OutputsModel.R
+++ b/R/plot.OutputsModel.R
@@ -1,20 +1,30 @@
 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, verbose = TRUE,
+                              LayoutMat = NULL, LayoutWidths = rep.int(1, ncol(LayoutMat)), LayoutHeights = rep.int(1, nrow(LayoutMat)), ...) {
+  
+  
+  ## 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 +38,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,6 +68,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   }
   
   
+  ## check 'which'
   if (is.null(which)) {
     stop("'which' must be a vector of character")
   }
@@ -76,7 +88,6 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     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")
   }
@@ -90,7 +101,6 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     which <- c("Error", "Regime", "CumFreq", "CorQQ")
   }
   
-  
   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'")
@@ -101,7 +111,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     }
   }
   
-  
+  ## check dates
   if (!BOOL_Dates) {
     stop("'OutputsModel' must contain at least 'DatesR' to allow plotting")
   }
@@ -190,7 +200,11 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   }
   BOOL_FilterZero <- TRUE
   
-  ## Plots_choices
+  
+  
+  ## ---------- plot
+  
+  ## plot choices
   BOOLPLOT_Precip   <- "Precip"   %in% which & BOOL_Pobs
   BOOLPLOT_PotEvap  <- "PotEvap"  %in% which & BOOL_Eobs
   BOOLPLOT_Temp     <- "Temp"     %in% which & BOOL_Snow
@@ -202,74 +216,82 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   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
-    }
-    if (Sum1 == 0 & Sum2 == 1) {
-      matlayout <- rbind(matlayout, iPlot + 1)
-      iPlot <- iPlot + 1
-    }
-    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)
-    #}
+    ## Set plot arrangement
+    if (is.null(LayoutMat)) {     
+      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
+      }
+      if (Sum1 == 0 & Sum2 == 1) {
+        matlayout <- rbind(matlayout, iPlot + 1)
+        iPlot <- iPlot + 1
+      }
+      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)
+      #}
+      
+      LayoutWidths  <- rep.int(1, ncol(matlayout))
+      LayoutHeights <- rep.int(1, nrow(matlayout))
+    }
+    if (!is.null(LayoutMat)) {
+      matlayout <- LayoutMat
+    }
+    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 +431,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) {
@@ -839,13 +861,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)
   
 }