diff --git a/NAMESPACE b/NAMESPACE
index dc19affe404c304a5ab2f64f7d103b86dfeae1dc..a6e3b9618b3e40ac965d2462cf710b9d5b24c35a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -38,6 +38,7 @@ export(RunModel)
 export(getNoSD_Ids)
 export(getSD_Ids)
 export(isNodeDownstream)
+export(plot.Qm3s)
 import(airGR)
 importFrom(grDevices,rainbow)
 importFrom(graphics,matplot)
diff --git a/R/plot.Qm3s.R b/R/plot.Qm3s.R
index 6b74e610a5c75924d79d0fa9d7adb1fe2007d919..0b378cacdbfde6546b236de4378ddff7c77cc1d1 100644
--- a/R/plot.Qm3s.R
+++ b/R/plot.Qm3s.R
@@ -1,12 +1,18 @@
 #' Plot of a `Qm3s` object (time series of simulated flows)
 #'
-#' @param x [data.frame] with a first column with [POSIXt] dates and followings columns with flows at each node of the network
+#' This function plot time series of flow rate in m3/s. It's a method for object
+#' of class "Qm3s" which can be directly called by `plot`. It can also be called
+#' as a function `plot.Qm3s` if the first parameter has the good format.
+#'
+#' @param x [data.frame] with a first column with [POSIXt] dates and followings
+#'        columns with flows at each node of the network
 #' @param type [character] plot type (See [plot.default]), default "l"
 #' @param xlab [character] label for the x axis, default to "Date"
 #' @param ylab [character] label for the y axis, default to "Flow (m3/s)"
 #' @param main [character] main title for the plot, default to "Simulated flows"
 #' @param col [character] plotting color (See [par]), default to rainbow colors
-#' @param legend [character] see parameter `legend` of [legend]. Set to [NULL] if display of the legend is not wanted
+#' @param legend [character] see parameter `legend` of [legend]. Set it to [NULL]
+#'        to hide the legend
 #' @param legend.cex [character] `cex` parameter for the text of the legend (See [par])
 #' @param legend.x,legend.y Legend position, see `x` and `y` parameters in [graphics::legend]
 #' @param lty [character] or [numeric] The line type (See [par])
@@ -18,12 +24,13 @@
 #'
 #' @importFrom grDevices rainbow
 #' @importFrom graphics matplot
+#' @export plot.Qm3s
 #' @export
 #'
 plot.Qm3s <- function(x,
                       type = "l",
                       xlab = "Date",
-                      ylab = expression("Flow (m"^"3"*"/s)"),
+                      ylab = expression("Flow rate (m"^"3"*"/s)"),
                       main = "Simulated flows",
                       col = rainbow(ncol(x) - 1),
                       legend = colnames(x)[-1],
@@ -32,8 +39,12 @@ plot.Qm3s <- function(x,
                       legend.y = NULL,
                       lty = 1,
                       ...) {
+
+  stopifnot(is.data.frame(x),
+            inherits(x[, 1], "POSIXct"))
+
   matplot(
-    x$DatesR,
+    x[, 1],
     x[, -1],
     type = type,
     lty = lty,
diff --git a/man-examples/RunModel.GRiwrmInputsModel.R b/man-examples/RunModel.GRiwrmInputsModel.R
index aec25ca1ee33c6bcd6c2b11e94b7d05eae677d25..9aafca1133e02849179b920c0f46fe4c3548c5c0 100644
--- a/man-examples/RunModel.GRiwrmInputsModel.R
+++ b/man-examples/RunModel.GRiwrmInputsModel.R
@@ -197,15 +197,14 @@ Ind_Plot <- which(
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
-class(dfQdiv) <- c("Qm3s", class(dfQdiv))
+
 oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
-plot(dfQdiv)
+plot.Qm3s(dfQdiv)
 
 # Plot natural and influenced flow at station "54001"
 df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
                  attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
 names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
-class(df54001) <- c("Qm3s", class(df54001))
-plot(df54001, ylim = c(0,70))
+plot.Qm3s(df54001, ylim = c(0,70))
 abline(h = 12, col = "green")
 par(oldpar)
diff --git a/man-examples/RunModel.Supervisor.R b/man-examples/RunModel.Supervisor.R
index 9db04cd28d97cb408a594ceed68800cac4edc4ee..d1aa0f1871052e07f55634bad8d67f126388d306 100644
--- a/man-examples/RunModel.Supervisor.R
+++ b/man-examples/RunModel.Supervisor.R
@@ -120,6 +120,5 @@ plot(attr(OM_dam, "Qm3s")[, c("DatesR", "54095", "Dam", "54057")],
      ylim = c(0, 200))
 Vres <- data.frame(DatesR = attr(OM_dam, "Qm3s")$DatesR,
                    V_reservoir = sv$Vres)
-class(Vres) <- c("Qm3s", class(Vres))
-plot(Vres)
+plot.Qm3s(Vres)
 par(oldpar)
diff --git a/man/CreateController.Rd b/man/CreateController.Rd
index fcd57d65ee81ad1a56e6e26d59fe4ab6a1272a4a..d474e737805e487df315f47ec164880a254e8cbd 100644
--- a/man/CreateController.Rd
+++ b/man/CreateController.Rd
@@ -165,7 +165,6 @@ plot(attr(OM_dam, "Qm3s")[, c("DatesR", "54095", "Dam", "54057")],
      ylim = c(0, 200))
 Vres <- data.frame(DatesR = attr(OM_dam, "Qm3s")$DatesR,
                    V_reservoir = sv$Vres)
-class(Vres) <- c("Qm3s", class(Vres))
-plot(Vres)
+plot.Qm3s(Vres)
 par(oldpar)
 }
diff --git a/man/CreateInputsModel.GRiwrm.Rd b/man/CreateInputsModel.GRiwrm.Rd
index 6e66cf02c4e2b6259de991daa7dc1e093bfb82ed..608f953b2828bddf235714cd12848cb27a501696 100644
--- a/man/CreateInputsModel.GRiwrm.Rd
+++ b/man/CreateInputsModel.GRiwrm.Rd
@@ -300,16 +300,15 @@ Ind_Plot <- which(
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
-class(dfQdiv) <- c("Qm3s", class(dfQdiv))
+
 oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
-plot(dfQdiv)
+plot.Qm3s(dfQdiv)
 
 # Plot natural and influenced flow at station "54001"
 df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
                  attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
 names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
-class(df54001) <- c("Qm3s", class(df54001))
-plot(df54001, ylim = c(0,70))
+plot.Qm3s(df54001, ylim = c(0,70))
 abline(h = 12, col = "green")
 par(oldpar)
 }
diff --git a/man/CreateRunOptions.Rd b/man/CreateRunOptions.Rd
index a25a127aff7e3b8b7f0abc0ea143d040b8129ae4..56ed9eb1c341e5f00ed43cdb418f9b8b12b94510 100644
--- a/man/CreateRunOptions.Rd
+++ b/man/CreateRunOptions.Rd
@@ -243,16 +243,15 @@ Ind_Plot <- which(
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
-class(dfQdiv) <- c("Qm3s", class(dfQdiv))
+
 oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
-plot(dfQdiv)
+plot.Qm3s(dfQdiv)
 
 # Plot natural and influenced flow at station "54001"
 df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
                  attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
 names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
-class(df54001) <- c("Qm3s", class(df54001))
-plot(df54001, ylim = c(0,70))
+plot.Qm3s(df54001, ylim = c(0,70))
 abline(h = 12, col = "green")
 par(oldpar)
 }
diff --git a/man/CreateSupervisor.Rd b/man/CreateSupervisor.Rd
index 6d28f25d15e738c236913f139af28c75a6ef1350..d68a3c417056a7af209568e5a14cff412a8de83e 100644
--- a/man/CreateSupervisor.Rd
+++ b/man/CreateSupervisor.Rd
@@ -147,7 +147,6 @@ plot(attr(OM_dam, "Qm3s")[, c("DatesR", "54095", "Dam", "54057")],
      ylim = c(0, 200))
 Vres <- data.frame(DatesR = attr(OM_dam, "Qm3s")$DatesR,
                    V_reservoir = sv$Vres)
-class(Vres) <- c("Qm3s", class(Vres))
-plot(Vres)
+plot.Qm3s(Vres)
 par(oldpar)
 }
diff --git a/man/RunModel.GRiwrmInputsModel.Rd b/man/RunModel.GRiwrmInputsModel.Rd
index d3e75da557470151142a9b17e5f68bc44c7a7d24..097b74dc95bb08e264716e443825d597c1c22039 100644
--- a/man/RunModel.GRiwrmInputsModel.Rd
+++ b/man/RunModel.GRiwrmInputsModel.Rd
@@ -221,16 +221,15 @@ Ind_Plot <- which(
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
-class(dfQdiv) <- c("Qm3s", class(dfQdiv))
+
 oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
-plot(dfQdiv)
+plot.Qm3s(dfQdiv)
 
 # Plot natural and influenced flow at station "54001"
 df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
                  attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
 names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
-class(df54001) <- c("Qm3s", class(df54001))
-plot(df54001, ylim = c(0,70))
+plot.Qm3s(df54001, ylim = c(0,70))
 abline(h = 12, col = "green")
 par(oldpar)
 }
diff --git a/man/RunModel.Supervisor.Rd b/man/RunModel.Supervisor.Rd
index 9bb72a134ffdfebd8d9d83b4b924eb8b94970ea5..667be3f8e295c1037529dd946745d93b120c2261 100644
--- a/man/RunModel.Supervisor.Rd
+++ b/man/RunModel.Supervisor.Rd
@@ -144,7 +144,6 @@ plot(attr(OM_dam, "Qm3s")[, c("DatesR", "54095", "Dam", "54057")],
      ylim = c(0, 200))
 Vres <- data.frame(DatesR = attr(OM_dam, "Qm3s")$DatesR,
                    V_reservoir = sv$Vres)
-class(Vres) <- c("Qm3s", class(Vres))
-plot(Vres)
+plot.Qm3s(Vres)
 par(oldpar)
 }
diff --git a/man/plot.GRiwrmOutputsModel.Rd b/man/plot.GRiwrmOutputsModel.Rd
index 1bf143c9222bb8c3af2e7906de22aa0398cb4bcd..60162ec0d3652076fb6744a7b3df8e67857e934e 100644
--- a/man/plot.GRiwrmOutputsModel.Rd
+++ b/man/plot.GRiwrmOutputsModel.Rd
@@ -221,16 +221,15 @@ Ind_Plot <- which(
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
-class(dfQdiv) <- c("Qm3s", class(dfQdiv))
+
 oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
-plot(dfQdiv)
+plot.Qm3s(dfQdiv)
 
 # Plot natural and influenced flow at station "54001"
 df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
                  attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
 names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
-class(df54001) <- c("Qm3s", class(df54001))
-plot(df54001, ylim = c(0,70))
+plot.Qm3s(df54001, ylim = c(0,70))
 abline(h = 12, col = "green")
 par(oldpar)
 }
diff --git a/man/plot.Qm3s.Rd b/man/plot.Qm3s.Rd
index 34682b1bddd62b4d35f107f52bd972081bc46759..acd5b4fa4e04bf30543c0783f48b7bacc4ffbad5 100644
--- a/man/plot.Qm3s.Rd
+++ b/man/plot.Qm3s.Rd
@@ -8,7 +8,7 @@
   x,
   type = "l",
   xlab = "Date",
-  ylab = expression("Flow (m"^"3" * "/s)"),
+  ylab = expression("Flow rate (m"^"3" * "/s)"),
   main = "Simulated flows",
   col = rainbow(ncol(x) - 1),
   legend = colnames(x)[-1],
@@ -20,7 +20,8 @@
 )
 }
 \arguments{
-\item{x}{\link{data.frame} with a first column with \link{POSIXt} dates and followings columns with flows at each node of the network}
+\item{x}{\link{data.frame} with a first column with \link{POSIXt} dates and followings
+columns with flows at each node of the network}
 
 \item{type}{\link{character} plot type (See \link{plot.default}), default "l"}
 
@@ -32,7 +33,8 @@
 
 \item{col}{\link{character} plotting color (See \link{par}), default to rainbow colors}
 
-\item{legend}{\link{character} see parameter \code{legend} of \link{legend}. Set to \link{NULL} if display of the legend is not wanted}
+\item{legend}{\link{character} see parameter \code{legend} of \link{legend}. Set it to \link{NULL}
+to hide the legend}
 
 \item{legend.cex}{\link{character} \code{cex} parameter for the text of the legend (See \link{par})}
 
@@ -46,7 +48,9 @@
 Screen plot window.
 }
 \description{
-Plot of a \code{Qm3s} object (time series of simulated flows)
+This function plot time series of flow rate in m3/s. It's a method for object
+of class "Qm3s" which can be directly called by \code{plot}. It can also be called
+as a function \code{plot.Qm3s} if the first parameter has the good format.
 }
 \examples{
 ###################################################################
@@ -248,16 +252,15 @@ Ind_Plot <- which(
 )
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR[Ind_Plot],
                      Diverted_flow = Qdiv_m3s[Ind_Plot])
-class(dfQdiv) <- c("Qm3s", class(dfQdiv))
+
 oldpar <- par(mfrow=c(2,1), mar = c(2.5,4,1,1))
-plot(dfQdiv)
+plot.Qm3s(dfQdiv)
 
 # Plot natural and influenced flow at station "54001"
 df54001 <- cbind(attr(OM_div, "Qm3s")[Ind_Plot, c("DatesR", "54001")],
                  attr(OM_severn, "Qm3s")[Ind_Plot, "54001"])
 names(df54001) <- c("DatesR", "54001 with irrigation", "54001 natural flow")
-class(df54001) <- c("Qm3s", class(df54001))
-plot(df54001, ylim = c(0,70))
+plot.Qm3s(df54001, ylim = c(0,70))
 abline(h = 12, col = "green")
 par(oldpar)
 }