diff --git a/R/plot.OutputsModel.R b/R/plot.OutputsModel.R
index 619e5ea9618613258f6f2d0928361468b89a29bd..86e3cd67a8cef69feb8da218a600ecc00d4b259e 100644
--- a/R/plot.OutputsModel.R
+++ b/R/plot.OutputsModel.R
@@ -1,50 +1,55 @@
 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,
+                              AxisTS = function(x) axis.POSIXct(side = 1, x = x$DatesR, ...),
                               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
-  
-  
-  
+  if (!is.null(IndPeriod_Plot)) {
+    OutputsModel <- .IndexOutputsModel(OutputsModel, IndPeriod_Plot)
+    Qobs <- Qobs[IndPeriod_Plot]
+    IndPeriod_Plot <- seq_along(IndPeriod_Plot)
+  }
+
+
   ## ---------- 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
   if ("Precip" %in% names(OutputsModel)) {
     BOOL_Pobs <- TRUE
   }
-  
+
   BOOL_EPobs <- FALSE
   if ("PotEvap" %in% names(OutputsModel)) {
     BOOL_EPobs <- TRUE
   }
-  
+
   BOOL_EAobs <- FALSE
   if ("AE" %in% names(OutputsModel)) {
     BOOL_EAobs <- TRUE
   }
-  
+
   BOOL_Qsim <- FALSE
   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)) {
@@ -53,27 +58,27 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   } else if (inherits(OutputsModel, "GR") & !is.null(Qobs)) {
     warning("incorrect length of 'Qobs'. Time series of observed flow not drawn")
   }
-  
+
   BOOL_Error <- FALSE
   if (BOOL_Qsim & BOOL_Qobs) {
     BOOL_Error <- TRUE
   }
-  
+
   BOOL_Snow <- FALSE
   if ("CemaNeigeLayers" %in% names(OutputsModel)) {
     if ("SnowPack" %in% names(OutputsModel$CemaNeigeLayers[[1]])) {
       BOOL_Snow <- TRUE
     }
   }
-  
+
   BOOL_Psol <- FALSE
   if ("CemaNeigeLayers" %in% names(OutputsModel)) {
     if ("Psol"     %in% names(OutputsModel$CemaNeigeLayers[[1]])) {
       BOOL_Psol <- TRUE
     }
   }
-  
-  
+
+
   ## check 'which'
   whichNeedQobs  <- c("Error", "CorQQ")
   whichDashboard <- c("all", "synth", "ts", "perf")
@@ -139,7 +144,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")
@@ -147,7 +152,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   if (inherits(OutputsModel, "GR") & !BOOL_Qsim) {
     stop("'OutputsModel' must contain at least 'Qsim' to allow plotting")
   }
-  
+
   if (BOOL_Dates) {
     # MyRollMean1 <- function(x, n) {
     #   return(filter(x, rep(1 / n, n), sides = 2))
@@ -157,7 +162,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     # }
     MyRollMean3 <- function(x, n) {
       return(filter(x, filter = rep(1 / n, n), sides = 2, circular = TRUE))
-    }    
+    }
     BOOL_TS  <- FALSE
     if (inherits(OutputsModel, "hourly")) {
       BOOL_TS <- TRUE
@@ -232,11 +237,11 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     warning("zeroes detected in 'Qsim': some plots in the log space will not be created using all time-steps")
   }
   BOOL_FilterZero <- TRUE
-  
-  
-  
+
+
+
   ## ---------- plot
-  
+
   ## plot choices
   BOOLPLOT_Precip   <- "Precip"   %in% which &  BOOL_Pobs
   BOOLPLOT_PotEvap  <- "PotEvap"  %in% which &  BOOL_EPobs
@@ -248,30 +253,30 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
   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
   BLOC <- TRUE
   if (BLOC) {
     lwdk <- 1.8
     line <- 2.6
     bg   <- NA
-    
+
     ## Set plot arrangement
-    if (is.null(LayoutMat)) {     
+    if (is.null(LayoutMat)) {
       matlayout <- NULL
       iPlot <- 0
       iHght <- NULL
-      
+
       listBOOLPLOT1 <- c(Precip = BOOLPLOT_Precip,
-                         PotEvap  = BOOLPLOT_PotEvap | BOOLPLOT_ActuEvap,
+                         PotEvap = BOOLPLOT_PotEvap | BOOLPLOT_ActuEvap,
                          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
@@ -305,24 +310,25 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       matlayout <- LayoutMat
     }
     layout(matlayout, widths = LayoutWidths, heights = LayoutHeights)
-    
 
-    
-    
-    Xaxis <- 1:length(IndPeriod_Plot)
-    if (BOOL_Dates) {
-      if (NameTS %in% c("hour", "day", "month")) {
-        Seq1 <- which(OutputsModel$DatesR[IndPeriod_Plot]$mday == 1 & OutputsModel$DatesR[IndPeriod_Plot]$mon %in% c(0, 3, 6, 9))
-        Seq2 <- which(OutputsModel$DatesR[IndPeriod_Plot]$mday == 1 & OutputsModel$DatesR[IndPeriod_Plot]$mon == 0)
-        Labels2 <- format(OutputsModel$DatesR[IndPeriod_Plot], format = formatAxis)[Seq2]
-      }
-      if (NameTS %in% c("year")) {
-        Seq1 <- 1:length(OutputsModel$DatesR[IndPeriod_Plot])
-        Seq2 <- 1:length(OutputsModel$DatesR[IndPeriod_Plot])
-        Labels2 <- format(OutputsModel$DatesR[IndPeriod_Plot], format = formatAxis)[Seq2]
-      }
-    }
-    
+
+
+
+    Xaxis <- as.POSIXct(OutputsModel$DatesR[IndPeriod_Plot])
+    # Xaxis <- 1:length(IndPeriod_Plot)
+    # if (BOOL_Dates) {
+    #   if (NameTS %in% c("hour", "day", "month")) {
+    #     Seq1 <- which(OutputsModel$DatesR[IndPeriod_Plot]$mday == 1 & OutputsModel$DatesR[IndPeriod_Plot]$mon %in% c(0, 3, 6, 9))
+    #     Seq2 <- which(OutputsModel$DatesR[IndPeriod_Plot]$mday == 1 & OutputsModel$DatesR[IndPeriod_Plot]$mon == 0)
+    #     Labels2 <- format(OutputsModel$DatesR[IndPeriod_Plot], format = formatAxis)[Seq2]
+    #   }
+    #   if (NameTS %in% c("year")) {
+    #     Seq1 <- 1:length(OutputsModel$DatesR[IndPeriod_Plot])
+    #     Seq2 <- 1:length(OutputsModel$DatesR[IndPeriod_Plot])
+    #     Labels2 <- format(OutputsModel$DatesR[IndPeriod_Plot], format = formatAxis)[Seq2]
+    #   }
+    # }
+
     if (!is.null(BasinArea)) {
       Factor_UNIT_M3S <- switch(NameTS,
                                 hour  = 60 * 60,
@@ -332,9 +338,9 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       Factor_UNIT_M3S <- BasinArea / (Factor_UNIT_M3S / 1000)
     }
   }
-  
+
   kPlot <- 0
-  
+
   ## vector of Q values for the y-axis when it is expressed in
   Factor <- ifelse(!is.null(BasinArea), Factor_UNIT_M3S, 1)
   seqDATA0 <- c(0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000)
@@ -344,17 +350,17 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     seqDATA1ba <- log(seqDATA0 * Factor_UNIT_M3S)
     seqDATA2ba <- round(exp(seqDATA1ba), digits = 2)
   }
-  
+
   ## Precip
   if (BOOLPLOT_Precip) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
-    
+
     par(new = FALSE, mar = mar)
     ylim1 <- range(OutputsModel$Precip[IndPeriod_Plot], na.rm = TRUE)
     ylim2 <- ylim1 * c(1.0, 1.1)
     ylim2 <- rev(ylim2)
-    
+
     lwdP <- lwd * 0.7
     if (NameTS %in% c("month", "year")) {
       lwdP <- lwd * 2
@@ -365,31 +371,33 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
          xlab = "", ylab = "", ...)
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
     mtext(side = 2, paste("precip.", plotunit), cex = cex.lab, adj = 1, line = line)
-    
+
     if (BOOL_Psol) {
       legend("bottomright", legend = c("solid","liquid"),
              col = c("lightblue", "royalblue"), lty = c(1, 1), lwd = c(lwd, lwd),
              bty = "o", bg = bg, box.col = bg, cex = cex.leg)
-      par(new = TRUE)
-      plot(Xaxis, PsolLayerMean[IndPeriod_Plot],
+      # par(new = TRUE)
+      # plot(Xaxis, PsolLayerMean[IndPeriod_Plot],
+      points(Xaxis, PsolLayerMean[IndPeriod_Plot],
            type = "h", xaxt = "n", yaxt = "n", yaxs = "i", ylim = ylim2,
            col = "lightblue", lwd = lwdP * lwdk, lend = 1,
            xlab = "", ylab = "", ...)
     }
-    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 (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, ...)
+    # }
+    AxisTS(OutputsModel)
   }
-  
-  
+
+
   ## PotEvap
   if (BOOLPLOT_PotEvap) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
-    
+
     par(new = FALSE, mar = mar)
     if (!BOOLPLOT_ActuEvap) {
       ylim1 <- range(OutputsModel$PotEvap[IndPeriod_Plot], na.rm = TRUE)
@@ -401,8 +409,8 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       xlabE <- "evap."
     }
     ylim2 <- ylim1 #* c(1.0, 1.1)
-    
-    
+
+
     plot(Xaxis, OutputsModel$PotEvap[IndPeriod_Plot],
          type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
          col = "green3", lwd = lwd * lwdk,
@@ -410,49 +418,51 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     if (BOOLPLOT_ActuEvap) {
       lines(Xaxis, OutputsModel$AE[IndPeriod_Plot],
            type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
-           col = "green4", lwd = lwd * lwdk, lty = 3) 
+           col = "green4", lwd = lwd * lwdk, lty = 3)
       legend("topright", legend = c("pot.", "actu."), col = c("green3", "green4"),
              lty = c(1, 3), lwd = c(lwd*1.0, lwd*0.8),
              bty = "o", bg = bg, box.col = bg, cex = cex.leg)
     }
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    
+
     mtext(side = 2, paste(xlabE, plotunit), cex = cex.lab, line = line)
-    
-    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 (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, ...)
+    # }
+    AxisTS(OutputsModel)
   }
-  
+
   ## ActuEvap
   if (BOOLPLOT_ActuEvap & !BOOLPLOT_PotEvap) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
-    
+
     par(new = FALSE, mar = mar)
     ylim1 <- range(OutputsModel$AE[IndPeriod_Plot], na.rm = TRUE)
     ylim2 <- ylim1 #* c(1.0, 1.1)
-    
+
     plot(Xaxis, OutputsModel$AE[IndPeriod_Plot],
          type = "l", xaxt = "n", yaxt = "n", ylim = ylim2,
          col = "green4", lwd = lwd * lwdk,
          xlab = "", ylab = "", ...)
     axis(side = 2, at = pretty(ylim1), labels = pretty(ylim1), cex.axis = cex.axis, ...)
-    
+
     mtext(side = 2, paste("actu. evap.", plotunit), cex = cex.lab, line = line)
-    
-    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 (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, ...)
+    # }
+    AxisTS(OutputsModel)
   }
-  
-  
+
+
   ## Temp
   if (BOOLPLOT_Temp) {
     kPlot <- kPlot + 1
@@ -468,29 +478,30 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
         SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$Temp/NLayers
       }
     }
-    plot(SnowPackLayerMean[IndPeriod_Plot], type = "n", ylim = ylim1, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
+    plot(Xaxis, SnowPackLayerMean[IndPeriod_Plot], type = "n", ylim = ylim1, xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
     for(iLayer in 1:NLayers) {
-      lines(OutputsModel$CemaNeigeLayers[[iLayer]]$Temp[IndPeriod_Plot], lty = 3, col = "orchid", lwd = lwd * lwdk * 0.8)
+      lines(Xaxis, OutputsModel$CemaNeigeLayers[[iLayer]]$Temp[IndPeriod_Plot], lty = 3, col = "orchid", lwd = lwd * lwdk * 0.8)
     }
     abline(h = 0, col = "grey", lty = 2)
-    lines(SnowPackLayerMean[IndPeriod_Plot], type = "l", lwd = lwd * lwdk *1.0, col = "darkorchid4")
+    lines(Xaxis, 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, ...)
-    
+
     mtext(side = 2, expression(paste("temp. [", degree, "C]"), sep = ""),  padj = 0.2, line = line, cex = cex.lab)
-    
+
     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)
     box()
-    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 (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, ...)
+    # }
+    AxisTS(OutputsModel)
   }
-  
-  
+
+
   ## SnowPack
   if (BOOLPLOT_SnowPack) {
     kPlot <- kPlot + 1
@@ -506,33 +517,34 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
         SnowPackLayerMean <- SnowPackLayerMean + OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack/NLayers
       }
     }
-    plot(SnowPackLayerMean[IndPeriod_Plot], type = "l", ylim = ylim1, lwd = lwd * lwdk *1.2, col = "royalblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
+    plot(Xaxis, SnowPackLayerMean[IndPeriod_Plot], type = "l", ylim = ylim1, lwd = lwd * lwdk *1.2, col = "royalblue", xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
     for(iLayer in 1:NLayers) {
-      lines(OutputsModel$CemaNeigeLayers[[iLayer]]$SnowPack[IndPeriod_Plot], lty = 3, col = "royalblue", lwd = lwd * lwdk *0.8)
+      lines(Xaxis, 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, ...)
-    
+
     mtext(side = 2, paste("snow pack", "[mm]"), line = line, cex = cex.lab)
     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)
     box()
-    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 (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, ...)
+    # }
+    AxisTS(OutputsModel)
   }
-  
-  
+
+
   ## Flows
   if (BOOLPLOT_Flows & log_scale) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
     par(new = FALSE, mar = mar)
-    
-    if (BOOL_Qobs) {    
+
+    if (BOOL_Qobs) {
       DATA2 <- Qobs
       DATA2[!SelectQobsNotZero] <- mean(Qobs, na.rm = TRUE) / 10000
       DATA2 <- log(DATA2)
@@ -540,7 +552,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     DATA3 <- OutputsModel$Qsim
     DATA3[!SelectQsimNotZero] <- mean(OutputsModel$Qsim, na.rm = TRUE) / 10000
     DATA3 <- log(DATA3)
-    
+
     ylim1 <- range(DATA3[IndPeriod_Plot], na.rm = TRUE)
     if (BOOL_Qobs) {
       ylim1 <- range(c(ylim1, DATA2[IndPeriod_Plot]), na.rm = TRUE)
@@ -566,12 +578,13 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       axis(side = 4, at = seqDATA1ba, labels = seqDATA2ba, cex.axis = cex.axis, ...)
       mtext(side = 4, paste("flow", "[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 (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, ...)
+    # }
+    AxisTS(OutputsModel)
     legend("topright", txtleg, col = colleg, lty = 1, lwd = lwd * lwdk , 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)
     box()
@@ -605,22 +618,23 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       axis(side = 4, at = pretty(ylim1*Factor)/Factor, labels = pretty(ylim1*Factor), cex.axis = cex.axis, ...)
       mtext(side = 4, paste("flow", "[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 (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, ...)
+    # }
+    AxisTS(OutputsModel)
     legend("topright", txtleg, col = colleg, lty = 1, lwd = lwd * lwdk , bty = "o", bg = bg, box.col = bg, cex = cex.leg)
     box()
   }
-  
-  
+
+
   ## Error
   if (BOOLPLOT_Error) {
     kPlot <- kPlot + 1
     mar <- c(3, 5, 1, 5)
-    
+
     if (log_scale) {
       errorQ <- log(OutputsModel$Qsim[IndPeriod_Plot]) - log(Qobs[IndPeriod_Plot])
     } else {
@@ -640,18 +654,19 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       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 (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, ...)
+    # }
+    AxisTS(OutputsModel)
     if (log_scale) {
       legend("bottomright", "log scale", lty = 1, col = NA, bty = "o", bg = bg, box.col = bg, cex = cex.leg)
     }
   }
-  
-  
+
+
   ## Regime
   if (BOOLPLOT_Regime) {
     kPlot <- kPlot + 1
@@ -722,20 +737,20 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
                                                      format = "%m%d")))
         DataDailyInterAn <- as.data.frame(aggregate(DataDaily[, 2:5], by = list(as.numeric(substr(TxtDatesDataDaily , 5, 8))), FUN = mean, na.rm = TRUE))
         colnames(DataDailyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
-        DataDailyInterAn <- merge(SeqY, DataDailyInterAn, by = "Dates", all.x = TRUE, all.y = FALSE) 
+        DataDailyInterAn <- merge(SeqY, DataDailyInterAn, by = "Dates", all.x = TRUE, all.y = FALSE)
       }
       if (!is.null(DataMonthly)) {
         SeqM <- data.frame(Dates = 1:12)
         DataMonthlyInterAn <- as.data.frame(aggregate(DataMonthly[, 2:5], by = list(as.numeric(substr(TxtDatesDataMonthly, 5, 6))), FUN = mean, na.rm = TRUE))
         colnames(DataMonthlyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
-        DataMonthlyInterAn <- merge(SeqM, DataMonthlyInterAn, by = "Dates", all.x = TRUE, all.y = FALSE) 
+        DataMonthlyInterAn <- merge(SeqM, DataMonthlyInterAn, by = "Dates", all.x = TRUE, all.y = FALSE)
       }
       ## Smoothing_of_daily_series_and_scale_conversion_to_make_them_become_a_monthly_regime
       if (!is.null(DataDaily)) {
         ## Smoothing
         NDaysWindow <- 30
-        DataDailyInterAn <- as.data.frame(cbind(DataDailyInterAn$Dates, 
-                                                MyRollMean3(DataDailyInterAn$Precip, NDaysWindow), MyRollMean3(DataDailyInterAn$Psol, NDaysWindow), 
+        DataDailyInterAn <- as.data.frame(cbind(DataDailyInterAn$Dates,
+                                                MyRollMean3(DataDailyInterAn$Precip, NDaysWindow), MyRollMean3(DataDailyInterAn$Psol, NDaysWindow),
                                                 MyRollMean3(DataDailyInterAn$Qobs  , NDaysWindow), MyRollMean3(DataDailyInterAn$Qsim, NDaysWindow)))
         colnames(DataDailyInterAn) <- c("Dates", "Precip", "Psol", "Qobs", "Qsim")
         ## Scale_conversion_to_make_them_become_a_monthly_regime
@@ -823,8 +838,8 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       box()
     }
   }
-  
-  
+
+
   ## Cumulative_frequency
   if (BOOLPLOT_CumFreq) {
     kPlot <- kPlot + 1
@@ -851,7 +866,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
            xlab = "", ylab = "", ...)
       axis(side = 1, at = pretty(xlim), labels = pretty(xlim), cex.axis = cex.axis, ...)
       mtext(side = 1, text = "non-exceedance prob. [-]", line = line, cex = cex.lab)
-      axis(side = 2, at = seqDATA1, labels = seqDATA2, cex.axis = cex.axis, ...) 
+      axis(side = 2, at = seqDATA1, labels = seqDATA2, cex.axis = cex.axis, ...)
       mtext(side = 2, text = paste("flow", plotunit), line = line, cex = cex.lab)
       txtleg <- NULL
       colleg <- NULL
@@ -896,8 +911,8 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     }
     box()
   }
-  
-  
+
+
   ## Correlation_QQ
   if (BOOLPLOT_CorQQ) {
     kPlot <- kPlot + 1
@@ -928,7 +943,7 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
     }
     box()
   }
-  
+
   ## Empty_plots
   if (exists("iPlotMax")) {
     while (kPlot < iPlotMax) {
@@ -937,8 +952,8 @@ plot.OutputsModel <- function(x, Qobs = NULL, IndPeriod_Plot = NULL, BasinArea =
       plot(0, 0, type = "n", xlab = "", ylab = "", axes = FALSE, ...)
     }
   }
-  
+
   ## Restoring_layout_options
   # layout(1)
-  
+
 }
diff --git a/man/plot.OutputsModel.Rd b/man/plot.OutputsModel.Rd
index e9a9042a7279fce9baa7bc6a3492a9489d5ec046..10c033c9dffe48498a9a36c06c842031454aefc5 100644
--- a/man/plot.OutputsModel.Rd
+++ b/man/plot.OutputsModel.Rd
@@ -21,6 +21,7 @@ Function which creates a screen plot giving an overview of the model outputs.
 \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,
+     AxisTS = function(x) axis.POSIXct(side = 1, x = x$DatesR, ...),
      LayoutMat = NULL,
      LayoutWidths = rep.int(1, ncol(LayoutMat)),
      LayoutHeights = rep.int(1, nrow(LayoutMat)),
@@ -49,6 +50,8 @@ Function which creates a screen plot giving an overview of the model outputs.
 
 \item{lwd}{(optional) [numeric] the line width (a positive number)}
 
+\item{AxisTS}{(optional) [function] to manage x-axis representing calendar dates and times on time series plots (see \code{\link{axis}} and \code{\link{axis.POSIXct}})}
+
 \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}})}