diff --git a/DESCRIPTION b/DESCRIPTION
index 77af64aff30b4796be79ee59310b9f39569c8db5..763dff675d4f612685d559554c7dfdbf5c16ae2e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.2.15.3
+Version: 1.2.15.4
 Date: 2019-05-02
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
diff --git a/NEWS.rmd b/NEWS.rmd
index 8c6ace522086fb2e1ec70e368e8e91c8d5c09ef3..f59e604bbf49419a66b74f1113c0bab09f3047ef 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -14,7 +14,7 @@ output:
 
 
 
-### 1.2.15.3 Release Notes (2019-05-02)
+### 1.2.15.34 Release Notes (2019-05-02)
 
 
 #### New features
diff --git a/R/plot.OutputsModel.R b/R/plot.OutputsModel.R
index 7a7ca69d26e513e235c772c0dc0079c057e273cd..3e0af2c1359c08fe4b849993e65a61ff59352d8f 100644
--- a/R/plot.OutputsModel.R
+++ b/R/plot.OutputsModel.R
@@ -315,7 +315,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
     
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     ylim1 <- range(OutputsModel$Precip[IndPeriod_Plot], na.rm = TRUE)
     ylim2 <- ylim1 * c(1.0, 1.1)
     ylim2 <- rev(ylim2)
@@ -329,9 +329,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
          col = "royalblue", lwd = lwdP * lwdk, lend = 1,
          xlab = "", ylab = "", ...)
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    par(las = 0)
     mtext(side = 2, paste("precip.", plotunit), cex = cex.lab, adj = 1, line = line)
-    par(las = 0)
     
     if (BOOL_Psol) {
       legend("bottomright", legend = c("solid","liquid"),
@@ -357,7 +355,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
     
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     ylim1 <- range(OutputsModel$PotEvap[IndPeriod_Plot], na.rm = TRUE)
     ylim2 <- ylim1 #* c(1.0, 1.1)
     
@@ -366,9 +364,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
          col = "green3", lwd = lwd * lwdk,
          xlab = "", ylab = "", ...)
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    par(las = 0)
+    
     mtext(side = 2, paste("pot. evap.", plotunit), cex = cex.lab, line = line)
-    par(las = 0)
+    
     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, ...)
@@ -382,7 +380,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   if (BOOLPLOT_Temp) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     ylim1 <- c(+99999, -99999)
     for(iLayer in 1:NLayers) {
       ylim1[1] <- min(ylim1[1], OutputsModel$CemaNeigeLayers[[iLayer]]$Temp)
@@ -400,9 +398,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     abline(h = 0, col = "grey", lty = 2)
     lines(SnowPackLayerMean[IndPeriod_Plot], type = "l", lwd = lwd * lwdk *1.0, col = "darkorchid4")
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    par(las = 0)
+    
     mtext(side = 2, expression(paste("temp. [", degree, "C]"), sep = ""),  padj = 0.2, line = line, cex = cex.lab)
-    par(las = 0)
+    
     legend("topright", legend = c("mean", "layers"), col = c("darkorchid4", "orchid"),
            lty = c(1, 3), lwd = c(lwd*1.0, lwd*0.8),
            bty = "o", bg = bg, box.col = bg, cex = cex.leg)
@@ -420,7 +418,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   if (BOOLPLOT_SnowPack) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     ylim1 <- c(+99999, -99999)
     for(iLayer in 1:NLayers) {
       ylim1[1] <- min(ylim1[1], OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack)
@@ -436,9 +434,8 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       lines(OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack[IndPeriod_Plot], lty = 3, col = "royalblue", lwd = lwd * lwdk *0.8)
     }
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    par(las = 0)
+    
     mtext(side = 2, paste("snow pack", "[mm]"), line = line, cex = cex.lab)
-    par(las = 0)
     legend("topright", legend = c("mean", "layers"), col = c("royalblue", "royalblue"),
            lty = c(1, 3), lwd = c(lwd*1.2, lwd*0.8),
            bty = "o", bg = bg, box.col = bg, cex = cex.leg)
@@ -456,7 +453,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   if (BOOLPLOT_Flows & log_scale) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     
     DATA2 <- Qobs
     DATA2[!SelectQobsNotZero] <- mean(Qobs, na.rm = TRUE) / 10000
@@ -485,15 +482,11 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       colleg <- c(colleg, "orangered")
     }
     axis(side = 2, at = seqDATA1, labels = seqDATA2, cex.axis = cex.axis, ...)
-    par(las = 0)
     mtext(side = 2, paste("flow", plotunit), line = line, cex = cex.lab)
-    par(las = 0)
     if (!is.null(BasinArea)) {
       Factor <- Factor_UNIT_M3S
       axis(side = 4, at = seqDATA1ba, labels = seqDATA2ba, cex.axis = cex.axis, ...)
-      par(las = 0)
       mtext(side = 4, paste("flow", "[m3/s]"), line = line, cex = cex.lab)
-      par(las = 0)
     }
     if (BOOL_Dates) {
       axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cex.axis, ...)
@@ -508,7 +501,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   if (BOOLPLOT_Flows & !log_scale) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     ylim1 <- range(OutputsModel$Qsim[IndPeriod_Plot], na.rm = TRUE)
     if (BOOL_Qobs) {
       ylim1 <- range(c(ylim1, Qobs[IndPeriod_Plot]), na.rm = TRUE)
@@ -528,15 +521,11 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       colleg <- c(colleg, "orangered")
     }
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    par(las = 0)
     mtext(side = 2, paste("flow", plotunit), line = line, cex = cex.lab)
-    par(las = 0)
     if (!is.null(BasinArea)) {
       Factor <- Factor_UNIT_M3S
       axis(side = 4, at = pretty(ylim1*Factor)/Factor, labels = pretty(ylim1*Factor), cex.axis = cex.axis, ...)
-      par(las = 0)
       mtext(side = 4, paste("flow", "[m3/s]"), line = line, cex = cex.lab)
-      par(las = 0)
     }
     if (BOOL_Dates) {
       axis(side = 1, at = Seq1, labels = FALSE, cex.axis = cex.axis, ...)
@@ -555,7 +544,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     mar <- c(3, 5, 1, 5)
     
     errorQ <- OutputsModel$Qsim / Qobs
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     ylim1 <- range(errorQ[IndPeriod_Plot], na.rm = TRUE)
     
     plot(Xaxis, errorQ[IndPeriod_Plot],
@@ -564,9 +553,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
          xlab = "", ylab = "", log = ifelse(log_scale, "y", ""),
          panel.first = abline(h = 1, col = "grey", lty = 2), ...)
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    par(las = 0)
     mtext(side = 2, paste("flow err.", plotunit), cex = cex.lab, line = line)
-    par(las = 0)
     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, ...)
@@ -581,25 +568,20 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     kPlot <- kPlot + 1
     mar <- c(6, 5, 1, 5)
     plotunitregime <- "[mm/month]"
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     ## Empty plot
     if ((NameTS == "hour"  & length(IndPeriod_Plot) < 697) |
         (NameTS == "day"   & length(IndPeriod_Plot) <  30) |
         (NameTS == "month" & length(IndPeriod_Plot) <   2) |
         (NameTS == "year"  & length(IndPeriod_Plot) <   2)) {
       plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
-      par(las = 0)
       mtext(side = 1, text = "", line = line, cex = cex.lab)
-      par(las = 0)
       text(0, 0, labels = "NO ENOUGH VALUES", col = "grey40")
-      par(las = 0)
       txtlab <- "flow regime"
       if (BOOL_Pobs) {
         txtlab <- "precip. & flow regime"
       }
-      par(las = 0)
       mtext(side = 2, paste(txtlab, plotunitregime), line = line, cex = cex.lab)
-      par(las = 0)
     } else {
       ## Data_formating_as_table
       DataModel <- as.data.frame(matrix(as.numeric(NA), nrow = length(IndPeriod_Plot), ncol = 5))
@@ -734,24 +716,18 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       axis(side = 1, at = SeqX1, tick = TRUE , labels = xLabels1, cex.axis = cex.axis, ...)
       axis(side = 1, at = SeqX2, tick = FALSE, labels = xLabels2, cex.axis = cex.axis, ...)
       axis(side = 2, at = pretty(ylimQ), labels = pretty(ylimQ), cex.axis = cex.axis, ...)
-      par(las = 0)
       mtext(side = 1, labX, line = line, cex = cex.lab)
-      par(las = 0)
       posleg <- "topright"
       txtlab <- "flow regime"
       if (BOOL_Pobs) {
         posleg <- "right"
         txtlab <- "precip. & flow regime"
       }
-      par(las = 0)
       mtext(side = 2, paste(txtlab, plotunitregime), line = line, cex = cex.lab)
-      par(las = 0)
       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, ...)
-        par(las = 0)
         mtext(side = 4, paste("flow regime", "[m3/s]"), line = line, cex = cex.lab)
-        par(las = 0)
       }
       ### 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)
@@ -765,7 +741,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   if (BOOLPLOT_CumFreq) {
     kPlot <- kPlot + 1
     mar <- c(6, 5, 1, 5)
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     xlim <- c(0, 1)
     if ( BOOL_Qobs & !BOOL_Qsim) {
       SelectNotZero <- SelectQobsNotZero
@@ -786,13 +762,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
            xaxt = "n", yaxt = "n",
            xlab = "", ylab = "", ...)
       axis(side = 1, at = pretty(xlim), labels = pretty(xlim), cex.axis = cex.axis, ...)
-      par(las = 0)
       mtext(side = 1, text = "non-exceedance prob. [-]", line = line, cex = cex.lab)
-      par(las = 0)
       axis(side = 2, at = seqDATA1, labels = seqDATA2, cex.axis = cex.axis, ...) 
-      par(las = 0)
       mtext(side = 2, text = paste("flow", plotunit), line = line, cex = cex.lab)
-      par(las = 0)
       txtleg <- NULL
       colleg <- NULL
       if (BOOL_Qobs) {
@@ -824,20 +796,14 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       if (!is.null(BasinArea)) {
         Factor <- Factor_UNIT_M3S
         axis(side = 4, at = seqDATA1, labels = round(seqDATA2*Factor, digits = 2), cex.axis = cex.axis, ...)
-        par(las = 0)
         mtext(side = 4, paste("flow", "[m3/s]"), line = line, cex = cex.lab)
-        par(las = 0)
       }
       legend("topleft", txtleg, col = colleg, lty = 1, lwd = lwd, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
       legend("bottomright", "log scale", lty = 1, col = NA, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
     } else {
       plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
-      par(las = 0)
       mtext(side = 1, text = "non-exceedance prob. [-]", line = line, cex = cex.lab)
-      par(las = 0)
-      par(las = 0)
       mtext(side = 2, text = paste("flow", plotunit), line = line, cex = cex.lab)
-      par(las = 0)
       text(0, 0, labels = "NO COMMON DATA", col = "grey40")
     }
     box()
@@ -848,7 +814,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   if (BOOLPLOT_CorQQ) {
     kPlot <- kPlot + 1
     mar <- c(6, 5, 1, 5)
-    par(new = FALSE, mar = mar, las = 0)
+    par(new = FALSE, mar = mar)
     if (any(SelectNotZero)) {
       ylim <- log(range(c(Qobs[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero], OutputsModel$Qsim[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero]), na.rm = TRUE))
       plot(log(Qobs[IndPeriod_Plot][SelectQobsNotZero & SelectQsimNotZero]),
@@ -858,28 +824,18 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       abline(a = 0, b = 1, col = "royalblue", lwd = lwd)
       axis(side = 1, at = seqDATA1, labels = seqDATA2, cex = cex.leg, cex.axis = cex.axis, ...)
       axis(side = 2, at = seqDATA1, labels = seqDATA2, cex = cex.leg, cex.axis = cex.axis, ...)
-      par(las = 0)
       mtext(side = 1, paste("observed flow", plotunit), line = line, cex = cex.lab)
-      par(las = 0)
-      par(las = 0)
       mtext(side = 2, paste("simulated flow", plotunit), line = line, cex = cex.lab)
-      par(las = 0)
       if (!is.null(BasinArea)) {
         Factor <- Factor_UNIT_M3S
         axis(side = 4, at = seqDATA1, labels = round(seqDATA2*Factor, digits = 2), cex.axis = cex.axis, ...)
-        par(las = 0)
         mtext(side = 4, paste("simulated flow", "[m3/s]"), line = line, cex = cex.lab)
-        par(las = 0)
       }
       legend("bottomright", "log scale", lty = 1, col = NA, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
     } else {
       plot(0, 0, type = "n", xaxt = "n", yaxt = "n", xlab = "", ylab = "", ...)
-      par(las = 0)
       mtext(side = 1, paste("observed flow", plotunit), line = line, cex = cex.lab)
-      par(las = 0)
-      par(las = 0)
       mtext(side = 2, paste("simulated flow", plotunit), line = line, cex = cex.lab)
-      par(las = 0)
       text(0, 0, labels = "NO COMMON DATA", col = "grey40")
     }
     box()