diff --git a/.gitignore b/.gitignore
index 72342a65dbae583572a67520e9566b7e5bd68583..f10d27f9f1b02ad12198493cb5b9ecba808ae567 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,15 +1,3 @@
-# Specific files for airGR
-packrat/lib*/
-
-# Compiled files
-/src/*.o
-/src/*.dll
-/src-*
-
-# Test temporary files
-/tests/tmp/
-/tests/testthat/*.pdf
-
 ######################################################################################################
 ### Generic .gitignore for R (source: https://github.com/github/gitignore/blob/master/R.gitignore) ###
 ######################################################################################################
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index bc8fed2b4448287abe5d2028d4d374113314d29b..fc25b6bb26cceabfc087f54e18d3e59910e115f1 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -13,6 +13,7 @@ cache:
     - R_LIBS/
 
 before_script:
+  - apt-get update && apt-get -y install libxml2-dev
   - echo "R_LIBS='$(pwd)/R_LIBS'" > .Renviron
 
 update_packages:
@@ -20,7 +21,6 @@ update_packages:
   script:
     - mkdir -p R_LIBS
     - Rscript -e 'if(!dir.exists("R_LIBS/remotes")) install.packages("remotes", lib = "R_LIBS")'
-    - apt-get update && apt-get -y install libxml2-dev
     - Rscript -e 'remotes::update_packages(c("knitr", "rmarkdown"), lib = "R_LIBS")'
     - Rscript -e 'remotes::install_gitlab("HYCAR-Hydro/airgr@dev", host = "gitlab.irstea.fr", lib = "R_LIBS")'
 
diff --git a/DESCRIPTION b/DESCRIPTION
index 46c3dfef00e9ec0633efb027ec84e49a703ff90e..c05b36498c68cda6d8b742fc0df641112382617a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,17 +1,24 @@
 Package: airGRdatassim
 Type: Package
-Title: Suite of Tools to Perform Ensemble-Based Data Assimilation with GR Hydrological Models
-Version: 0.0.2.44
-Date: 2020-12-04
+Title: Ensemble-Based Data Assimilation with GR Hydrological Models
+Version: 0.1.3
+Date: 2021-02-08
 Authors@R: c(
   person("Gaia", "Piazzi", role = c("aut"), comment = c(ORCID = "0000-0002-4167-0794")),
   person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
   person("Guillaume", "Thirel", role = c("ctb"), comment = c(ORCID = "0000-0002-1444-1830")),
   person("Maxime", "Logez", role = c("ctb"), comment = c(ORCID = "0000-0001-9843-0495"))
   )
-Depends: airGR (>= 1.4.3.65)
-Suggests: knitr, rmarkdown
-Description: Add-on to the 'airGR' package which provides the tools to assimilate observed discharges in daily GR hydrological model. The package consists in two functions allowing to perform the assimilation of observed discharges via Ensemble Kalman filter or Particle filter.
+Depends: airGR (>= 1.6.9.27)
+Imports: 
+  graphics,
+  grDevices,
+  stats,
+  utils
+Suggests:
+  knitr,
+  rmarkdown
+Description: Add-on to the 'airGR' package which provides the tools to assimilate observed discharges in daily GR hydrological models. The package consists in two functions allowing to perform the assimilation of observed discharges via the Ensemble Kalman filter or the Particle filter as described in Piazzi et al. (2021) <doi:10.1029/2020WR028390>.
 License: GPL-2
 BugReports: https://gitlab.irstea.fr/HYCAR-Hydro/airgrdatassim/-/issues
 NeedsCompilation: no
diff --git a/NAMESPACE b/NAMESPACE
index 9e4ceb0b8e6516746a1a11664c7146886a0fb396..c1642d0e088037c997eae4b94ae17f1dc221b019 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,8 @@
 #####################################
 ##            S3 methods           ##
 #####################################
-S3method('[', InputsModel)
+S3method('[', InputsPert)
+S3method('[', OutputsModelDA)
 S3method(plot, InputsPert)
 S3method(plot, OutputsModelDA)
 
diff --git a/R/CreateInputsPert.R b/R/CreateInputsPert.R
index a518c15ba5f36aff066653b0683f3102d1d18cb4..ae7932527da15eca26d89368498cd169ebd55da1 100644
--- a/R/CreateInputsPert.R
+++ b/R/CreateInputsPert.R
@@ -33,7 +33,7 @@ CreateInputsPert <- function(FUN_MOD, DatesR, Precip = NULL, PotEvap = NULL, Tem
   MeteoNames <- c("Precip", "PotEvap")
 
   # length of the input vectors
-  Nt <- length(DatesR)
+  NbTime <- length(DatesR)
 
   # model time step [day]
   Dt <- 1
@@ -57,11 +57,11 @@ CreateInputsPert <- function(FUN_MOD, DatesR, Precip = NULL, PotEvap = NULL, Tem
   skipVarMeteo <- NULL
   if (is.null(Precip)) {
     skipVarMeteo <- c(skipVarMeteo, "Precip")
-    Precip <- rep(0, times = Nt)
+    Precip <- rep(0, times = NbTime)
   }
   if (is.null(PotEvap)) {
     skipVarMeteo <- c(skipVarMeteo, "PotEvap")
-    PotEvap <- rep(0, times = Nt)
+    PotEvap <- rep(0, times = NbTime)
   }
 
   # select MeteoNames and Eps in function of the provided variables
@@ -94,21 +94,21 @@ CreateInputsPert <- function(FUN_MOD, DatesR, Precip = NULL, PotEvap = NULL, Tem
   MbrNames <- sprintf("Mbr_%s", seq_len(NbMbr))
 
   # time names
-  TimeNames <- sprintf("Time_%s", seq_len(Nt))
+  TimeNames <- sprintf("Time_%s", seq_len(NbTime))
 
 
   # error time evolution
-  S <- array(data = rep(NaN, times = NbMeteo*NbMbr*Nt),
-             dim = c(NbMeteo, NbMbr, Nt),
+  S <- array(data = rep(NaN, times = NbMeteo*NbMbr*NbTime),
+             dim = c(NbMeteo, NbMbr, NbTime),
              dimnames = list(MeteoNames, MbrNames, TimeNames))
   # uniform random numbers
-  U <- array(data = rep(NaN, times = NbMeteo*NbMbr*Nt),
-             dim = c(NbMeteo, NbMbr, Nt),
+  U <- array(data = rep(NaN, times = NbMeteo*NbMbr*NbTime),
+             dim = c(NbMeteo, NbMbr, NbTime),
              dimnames = list(MeteoNames, MbrNames, TimeNames))
 
   # multiplicative perturbations of model inputs
-  Fi <- array(data = rep(NaN, times = NbMeteo*NbMbr*Nt),
-              dim = c(NbMeteo, NbMbr, Nt),
+  Fi <- array(data = rep(NaN, times = NbMeteo*NbMbr*NbTime),
+              dim = c(NbMeteo, NbMbr, NbTime),
               dimnames = list(MeteoNames, MbrNames, TimeNames))
 
 
@@ -117,7 +117,7 @@ CreateInputsPert <- function(FUN_MOD, DatesR, Precip = NULL, PotEvap = NULL, Tem
 
   MeteoEns <- sapply(MeteoNames, function(iMeteo) {
     matrix(data = rep(InputsPert[[iMeteo]], each = NbMbr),
-           nrow = Nt, byrow = TRUE,
+           nrow = NbTime, byrow = TRUE,
            dimnames = list(TimeNames, MbrNames))
   }, simplify = "array")
 
@@ -128,7 +128,7 @@ CreateInputsPert <- function(FUN_MOD, DatesR, Precip = NULL, PotEvap = NULL, Tem
 
   # ------ Perturbation of the time series of model forcings (i.e., precipitation and potential evapotranspiration)
 
-  for (iTime in seq_len(Nt)) { # olivier, can we do apply here?
+  for (iTime in seq_len(NbTime)) { # olivier, can we do apply here?
 
     # white noise
     if (!is.null(Seed)) {
@@ -182,6 +182,6 @@ CreateInputsPert <- function(FUN_MOD, DatesR, Precip = NULL, PotEvap = NULL, Tem
 
   # ------ Class
 
-  class(InputsPert) <- c(class(InputsPert), "InputsPert")
+  class(InputsPert) <- c("InputsPert", class(InputsPert))
   return(InputsPert)
 }
diff --git a/R/RunModel_DA.R b/R/RunModel_DA.R
index 1c0d1d58122cf78d1929c4e8792068a7c09d6fd9..05c183296c70f2e304f17d4a4f2d0d8187b7ac32 100644
--- a/R/RunModel_DA.R
+++ b/R/RunModel_DA.R
@@ -20,7 +20,7 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
                        "RunModel_CemaNeigeGR6J")
 
   FUN_MOD <- match.fun(FUN_MOD)
-  if (!any(sapply(c(FUN_MODList,FUN_MODSnowList), function(x) identical(FUN_MOD, match.fun(x))))) {
+  if (!any(sapply(c(FUN_MODList, FUN_MODSnowList), function(x) identical(FUN_MOD, match.fun(x))))) {
     stop(sprintf("incorrect 'FUN_MOD' for use in 'CreateInputsPerturb'. Only %s can be used",
                  paste(c(FUN_MODList, FUN_MODSnowList), collapse = ", ")))
   }
@@ -142,7 +142,7 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
   # data assimilation method used (not open-loop simulation)
   IsDa <- DaMethod != "none"
 
-  Nt <- length(IndRun)
+  NbTime <- length(IndRun)
 
   NbState <- length(StateNames)
 
@@ -153,7 +153,7 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
   MbrNames <- sprintf("Mbr_%s", seq_len(NbMbr))
 
   # time names
-  TimeNames <- sprintf("Time_%s", seq_len(Nt))
+  TimeNames <- sprintf("Time_%s", seq_len(NbTime))
 
   # InputsModel
   InputsModel <- InputsModel[IndRun]
@@ -169,27 +169,22 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
   # ------ Ensemble initializations
 
   ObsPert <- matrix(data = NA,
-                    nrow = NbMbr, ncol = Nt,
+                    nrow = NbMbr, ncol = NbTime,
                     dimnames = list(MbrNames,
                                     TimeNames))
 
   QsimEns <- ObsPert
 
   IniStatesEns   <- list()
-  IniStatesEnsNt <- list()
+  IniStatesEnsNbTime <- list()
 
-  OutputStates <- matrix(data = NaN,
-                         nrow = NbState, ncol = NbMbr, byrow = TRUE,
-                         dimnames = list(StateNames,
-                                         MbrNames))
-
-  EnsStateBkg <- array(data = rep(NaN, times = NbState*NbMbr*Nt),
-                       dim = c(NbState, NbMbr, Nt),
+  EnsStateBkg <- array(data = rep(NaN, times = NbState*NbMbr*NbTime),
+                       dim = c(NbState, NbMbr, NbTime),
                        dimnames = list(StateNames,
                                        MbrNames,
                                        TimeNames))
 
-  EnsStateA   <- EnsStateBkg
+  EnsStateA <- EnsStateBkg
 
   ItAssim <- 0
 
@@ -230,7 +225,7 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
     }
     for (iMbr in seq_len(NbMbr)) {
 
-      if (iTime == 1) { #default (one year by default) warmup
+      if (iTime == 1) { # default (one year by default) warmup
 
         RunOptionsIni$IndPeriod_Run <- iTime
         OutputsModel <- FUN_MOD(InputsModel = InputsModel,
@@ -245,7 +240,7 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
         RunOptionsIter$IniStates <- IniStates
         RunOptionsIter$IniResLevels <- NULL
 
-        # Definition of run options
+        # definition of run options
         if (IsMeteo) {
           InputsPertMbr <- InputsPert
           InputsPertMbr$Precip <- InputsPert$Precip[, iMbr]
@@ -264,9 +259,9 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
 
       EnsStateBkg["Prod", iMbr, iTime] <- OutputsModel$Prod
       EnsStateBkg["Rout", iMbr, iTime] <- OutputsModel$Rout
-      EnsStateBkg["UH2"  , iMbr, iTime] <- OutputsModel$StateEnd$UH$UH2[1]
+      EnsStateBkg["UH2" , iMbr, iTime] <- OutputsModel$StateEnd$UH$UH2[1]
       if ("UH1" %in% StateNames) {
-        EnsStateBkg["UH1"  , iMbr, iTime] <- OutputsModel$StateEnd$UH$UH1[1]
+        EnsStateBkg["UH1", iMbr, iTime] <- OutputsModel$StateEnd$UH$UH1[1]
       }
 
       QsimEns[iMbr, iTime]  <- OutputsModel$Qsim
@@ -297,9 +292,9 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
           }
         }
 
-        if (iTime < Nt) {
-          IniStatesEnsNt[[iTime+1]] <- IniStatesEns
-          names(IniStatesEnsNt)[iTime+1] <- sprintf("Time_%s",iTime+1)
+        if (iTime < NbTime) {
+          IniStatesEnsNbTime[[iTime+1]] <- IniStatesEns
+          names(IniStatesEnsNbTime)[iTime+1] <- sprintf("Time_%s",iTime+1)
         }
 
         if (!is.null(StatePert)) {
@@ -315,7 +310,7 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
 
         EnsStateA[, , iTime] <- ans$EnsStateEnkf
 
-        if (iTime < Nt) {
+        if (iTime < NbTime) {
           if (!is.null(StatePert)) {
             EnsStateBkg[, , iTime+1] <- ans$EnsStatePert
           } else {
@@ -340,20 +335,20 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
         EnsStateA["Rout", , iTime] <- sapply(seq_along(ans$EnsStatePf), function(x) ans$EnsStatePf[[x]]$Store$Rout)
         EnsStateA["UH2" , , iTime] <- sapply(seq_along(ans$EnsStatePf), function(x) ans$EnsStatePf[[x]]$UH$UH2[1])
         if ("UH1" %in% StateNames) {
-          EnsStateA["UH1" , , iTime] <- sapply(seq_along(ans$EnsStatePf), function(x) ans$EnsStatePf[[x]]$UH$UH1[1])
+          EnsStateA["UH1", , iTime] <- sapply(seq_along(ans$EnsStatePf), function(x) ans$EnsStatePf[[x]]$UH$UH1[1])
         }
 
-        if (iTime < Nt) { # olivier?
-          IniStatesEnsNt[[iTime+1]] <- ans$EnsStatePf
-          names(IniStatesEnsNt)[iTime+1] <- sprintf("Time_%s", iTime+1)
+        if (iTime < NbTime) { # olivier?
+          IniStatesEnsNbTime[[iTime+1]] <- ans$EnsStatePf
+          names(IniStatesEnsNbTime)[iTime+1] <- sprintf("Time_%s", iTime+1)
         }
       }
 
     } else { # IF no assimilation
 
-      if (iTime < Nt) {
-        IniStatesEnsNt[[iTime+1]] <- IniStatesEns
-        names(IniStatesEnsNt)[iTime+1] <- sprintf("Time_%s", iTime+1)
+      if (iTime < NbTime) {
+        IniStatesEnsNbTime[[iTime+1]] <- IniStatesEns
+        names(IniStatesEnsNbTime)[iTime+1] <- sprintf("Time_%s", iTime+1)
 
         EnsStateBkg[, , iTime+1] <- EnsStateBkg[, , iTime]
       }
@@ -372,13 +367,13 @@ RunModel_DA <- function(InputsModel, InputsPert = NULL, Qobs = NULL,
   # ------ Outputs and class
 
   res <- list(DatesR = InputsModel$DatesR,
-              QsimEns = QsimEns,
-              EnsStateBkg = EnsStateBkg,
-              EnsStateA = EnsStateA,
+              QsimEns = t(QsimEns),
+              EnsStateBkg = aperm(EnsStateBkg),
+              EnsStateA = aperm(EnsStateA),
+              NbTime = NbTime,
               NbMbr = NbMbr,
-              Nt = Nt,
               NbState = NbState)
-  class(res) <- c("OutputsModelDA", DaMethod, TimeUnit)
+  class(res) <- c("OutputsModelDA", "OutputsModel", DaMethod, TimeUnit)
   return(res)
 
 }
diff --git a/R/Utils.R b/R/Utils.R
index 3a3c556ee25368e3495ff719164f52a5d22a4c4a..5c63f932e7b4ebecef4a5f84a9d8f7cf3c127968 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -1,27 +1,24 @@
-.ExtractInputsModel <- function(Inputs, IndRun) {
-  res <- lapply(Inputs, function(x) {
-    if (is.matrix(x)) {
-      res0 <- x[IndRun, ]
-    }
-    if (is.vector(x) | inherits(x, "POSIXt")) {
-      res0 <- x[IndRun]
-    }
-    if (is.list(x) & !inherits(x, "POSIXt")) {
-      res0 <- .ExtractInputsModel(Inputs = x, IndRun = IndRun)
-    }
-    return(res0)
-  })
-  if (!is.null(Inputs$ZLayers)) {
-    res$ZLayers <- Inputs$ZLayers
+
+## =================================================================================
+## function to extract parts of InputsPert or OutputsModelDA objects
+## =================================================================================
+
+'[.InputsPert' <-  function(x, i) {
+  res <- NextMethod()
+  if (is.numeric(i)) {
+    res$NbMbr <- x$NbMbr
   }
-  class(res) <- class(Inputs)
-  res
+  return(res)
 }
 
 
-'[.InputsModel' <- function(Inputs, IndRun) {
-  if (!inherits(Inputs, "InputsModel")) {
-    stop("'Inputs' must be of class 'InputsModel'")
+'[.OutputsModelDA' <- function(x, i) {
+  res <- NextMethod()
+  if (is.numeric(i)) {
+    res <- airGR:::.ExtractOutputsModel(x, i)
+    res$NbMbr   <- x$NbMbr
+    res$NbTime  <- x$NbTime
+    res$NbState <- x$NbState
   }
-  .ExtractInputsModel(Inputs, IndRun)
+  return(res)
 }
diff --git a/R/plot.InputsPert.R b/R/plot.InputsPert.R
index 4ef8d29506a59050b8f236d04cd5642e74d99b6f..fd071be5d73b8a5ff5490543713660d595498bd2 100644
--- a/R/plot.InputsPert.R
+++ b/R/plot.InputsPert.R
@@ -1,5 +1,7 @@
 plot.InputsPert <- function(x, which = "all", main = NULL,
-                            ColPrecip = "royalblue", ColPotEvap = "green3", ...) {
+                            ColPrecip = "royalblue", ColPotEvap = "green3",
+                            ask = prod(par("mfcol")) < length(which) && dev.interactive(),
+                            ...) {
 
 
   ## ---------- check arguments
@@ -21,6 +23,12 @@ plot.InputsPert <- function(x, which = "all", main = NULL,
     stop(sprintf("'%s' element not available in x", which))
   }
 
+  ## ask
+  if (length(NamesInputsPert) > 1L && ask) {
+    oask <- devAskNewPage(TRUE)
+    on.exit(devAskNewPage(oask))
+  }
+
 
   ## ---------- graphical variables
 
@@ -45,6 +53,7 @@ plot.InputsPert <- function(x, which = "all", main = NULL,
 
     ## ---------- plot
 
+    dev.hold()
     RangeInputsPert <- apply(x[[iName]], MARGIN = 1, FUN = range)
     plot(x = x$DatesR, y = rowMeans(x[[iName]]),
          ylim = range(RangeInputsPert),
@@ -54,6 +63,8 @@ plot.InputsPert <- function(x, which = "all", main = NULL,
          panel.first = polygon(x = c(as.numeric(x$DatesR), rev(as.numeric(x$DatesR))),
                                y = c(RangeInputsPert[1L, ], rev(RangeInputsPert[2L, ])),
                                col = adjustcolor(Col, alpha.f = 0.25), border = NA))
+    dev.flush()
+
   }
 
 
diff --git a/R/plot.OutputsModelDA.R b/R/plot.OutputsModelDA.R
index d3d56db16d6ddca325613f5b5a540829b0f961c7..99fd80137fcb1ffd58b2dd93c14077076297b1be 100644
--- a/R/plot.OutputsModelDA.R
+++ b/R/plot.OutputsModelDA.R
@@ -17,7 +17,7 @@ plot.OutputsModelDA <- function(x, Qobs = NULL, main = NULL,
 
   ## ---------- graphical variables
 
-  RangeQsimEns <- apply(x$QsimEns, MARGIN = 2, FUN = range)
+  RangeQsimEns <- apply(x$QsimEns, MARGIN = 1, FUN = range)
 
   if (!is.null(Qobs)) {
     LegObs <- "obs"
@@ -47,7 +47,7 @@ plot.OutputsModelDA <- function(x, Qobs = NULL, main = NULL,
 
   ## ---------- plot
 
-  plot(x = x$DatesR, y = colMeans(x$QsimEns),
+  plot(x = x$DatesR, y = rowMeans(x$QsimEns),
        ylim = range(RangeQsimEns, Qobs, na.rm = TRUE),
        type = "l", col = ColSim, lwd = 2,
        main = Main,
diff --git a/README.md b/README.md
index bc881e7653c7e7f9244e5ca7f5de608342c805fa..5c65606299946b9f53c2575d3da0bda064419d3c 100644
--- a/README.md
+++ b/README.md
@@ -1,18 +1,19 @@
 
-# airGRdatassim: Suite of tools to perform ensemble-based discharge assimilation in GR Hydrological Models
+# airGRdatassim: Suite of Tools to Perform Ensemble-Based Data Assimilation with GR Hydrological Models
 
 ## Overview
 
-airGRdatassim is a package based on the airGR hydrological modeling package. It provides the tools to assimilate observed discharges in the GR daily hydrological model (GR4J, GR5J and GR6J, with and without the CemaNeige snow model). The package is developed at INRAE-Antony ([Catchment Hydrology research group](https://webgr.inrae.fr/en/home/) of the HYCAR Research Unit, France). 
+airGRdatassim is a package based on the airGR hydrological modeling package. It provides the tools to assimilate observed discharges in the GR daily hydrological model (GR4J, GR5J and GR6J, with and without the CemaNeige snow model). The package is developed at INRAE-Antony ([Catchment Hydrology research group](https://webgr.inrae.fr/en/home/) of the HYCAR Research Unit, France). More information about the efficiency of these data assimilation schemes with GR5J can be found in Piazzi et al. (accepted)
 
 
 ## Installation
 
-To download the version of the airGRdatassim package that is on GitLab, you have first install the [Git software](https://git-scm.com/downloads). Then you can install the package in the R environment, using the following command lines:
+To download the version of the airGRdatassim package that is on GitLab, you have first install the [Git software](https://git-scm.com/downloads). Since you need the latest version of airGR (not yet on CRAN), you need to install [Rtools](https://cran.r-project.org/bin/windows/Rtools/history.html) in order to do so from its sources.  Then you can install the package in the R environment, using the following command lines:
 
 ``` r
 Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS="true")
 install.packages("remotes")
+remotes::install_git(url = "https://gitlab.irstea.fr/HYCAR-Hydro/airgr")
 remotes::install_git(url = "https://gitlab.irstea.fr/HYCAR-Hydro/airgrdatassim")
 ```
 
@@ -34,7 +35,7 @@ Consistently with the airGR package, both structure and class of function argume
 
 ## Hydrological model
 
-DA schemes are designed to be coupled with GR daily hydrological model, which is implemented in the airGR package. This model can be called within the airGRdatassim package using the following airGR functions (see the [airGR manual](https://cran.r-project.org/web/packages/airGR/airGR.pdf#page.4) to get the references of the GR models): 
+DA schemes are designed to be coupled with GR daily hydrological model, which is implemented in the airGR package. This model can be called within the airGRdatassim package using the following airGR functions (use the command `?airGR` to get the references of the GR models): 
 
   - `RunModel_GR4J()`: four-parameter daily lumped hydrological model
   - `RunModel_GR5J()`: five-parameter daily lumped hydrological model
@@ -59,4 +60,4 @@ For more information and to get started with the package, you can refer to the v
 
 ## References
 
-- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (in review). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model. Water Resources Research.
+- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (accepted). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model. Water Resources Research, doi: [10.1029/2020WR028390](https://doi.org/10.1029/2020WR028390).
diff --git a/inst/CITATION b/inst/CITATION
index f0261e270cd0bb106fbcf3ec519847ddd6cb4d20..14a531412329b43f868cd9d1fa5b62a6611810ef 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -5,31 +5,35 @@ note    <- sprintf("R package version %s", version)
 
 citHeader("To cite airGR in publications use these two references:")
 
-bibentry(bibtype="Article",
-         title = "Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertinties and to updated variables of a conceptual hydrological model",
-         author = personList(as.person("G. Piazzi"), as.person("G. Thirel"), as.person("C. Perrin"), as.person("O. Delaigue")),
-         journal = "Water Resources Research",
-         year = "in review",
-         textVersion = 
+bibentry(bibtype     = "Article",
+         title       = "Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertinties and to updated variables of a conceptual hydrological model",
+         author      = personList(as.person("G. Piazzi"), as.person("G. Thirel"), as.person("C. Perrin"), as.person("O. Delaigue")),
+         journal     = "Water Resources Research",
+         year        = "accepted",
+         doi         = "10.1029/2020WR028390",
+         textVersion =
            paste("Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O.",
-                 "(in review).", 
+                 "(accepted).",
                  "Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertinties and to updated variables of a conceptual hydrological model.",
-                 "Water Resources Research.",
+                 "Water Resources Research,",
+                 "doi: 10.1029/2020WR028390.",
                  sep = " ")
 )
 
-bibentry(bibtype="Manual",
-         title = "{airGRdatassim}: Suite of Tools to Perform Ensemble-Based Data Assimilation in {GR} Hydrological Models",
-         author = personList(as.person("G. Piazzi"), as.person("O. Delaigue")),
-         journal = "R News",
-         year = year,
-         note = note,       
-         url = "https://gitlab.irstea.fr/HYCAR-Hydro/airgrdatassim",
-         textVersion = 
-           paste0("Piazzi, G., Delaigue, O. (",
+bibentry(bibtype     = "Manual",
+         title       = "{airGRdatassim}: Suite of Tools to Perform Ensemble-Based Data Assimilation in {GR} Hydrological Models",
+         author      = personList(as.person("G. Piazzi"), as.person("O. Delaigue")),
+         journal     = "R News",
+         year        = year,
+         note        = note,
+         doi         = "10.15454/WEYYVZ",
+         url         = "https://gitlab.irstea.fr/HYCAR-Hydro/airgrdatassim",
+         textVersion =
+           paste0("Piazzi, G. and Delaigue, O. (",
                   year,
                   "). airGRdatassim: Suite of Tools to Perform Ensemble-Based Data Assimilation in GR Hydrological Models. ",
-                  note,
-                  ". URL: https://gitlab.irstea.fr/HYCAR-Hydro/airgrdatassim."
+                  note, ", ",
+                  "doi: 10.15454/WEYYVZ", ", ",
+                  "URL  : https://gitlab.irstea.fr/HYCAR-Hydro/airgrdatassim."
                   )
 )
diff --git a/man/CreateInputsPert.Rd b/man/CreateInputsPert.Rd
index 869c627f3d242abb202a96a86eb326c314d2314b..2f5e1c2fa2a2f77eee53bab2d64aaf0324ca10cd 100644
--- a/man/CreateInputsPert.Rd
+++ b/man/CreateInputsPert.Rd
@@ -3,19 +3,28 @@
 
 \name{CreateInputsPert}
 \alias{CreateInputsPert}
+\alias{[.InputsPert}
+\alias{plot.InputsPert}
 
 
 \title{Generation of ensemble model inputs for data assimilation}
 
+
 \usage{
 CreateInputsPert(FUN_MOD, DatesR, Precip = NULL,
                  PotEvap = NULL, TempMean = NULL,
                  ZInputs = NULL, HypsoData = NULL, NLayers = 5,
                  NbMbr = 50, Seed = NULL)
+
+\method{[}{InputsPert}(x, i)
+
+\method{plot}{InputsPert}(x, which = "all", main = NULL,
+     ColPrecip = "royalblue", ColPotEvap = "green3",
+     ask = prod(par("mfcol")) < length(which) && dev.interactive(), ...)
 }
 
-\arguments{
 
+\arguments{
 \item{FUN_MOD}{[function] hydrological model function (e.g. \code{\link[airGR]{RunModel_GR5J}}, \code{\link[airGR]{RunModel_CemaNeigeGR5J}})}
 
 \item{DatesR}{[POSIXct] vector of dates}
@@ -35,27 +44,42 @@ CreateInputsPert(FUN_MOD, DatesR, Precip = NULL,
 \item{NbMbr}{(optional) [numeric] number of ensemble members (minimum of 20 recommanded for the EnKF scheme and of 30 for the PF scheme)}
 
 \item{Seed}{(optional) [numeric] seed of random number generator}
+
+\item{x}{[InputsPert] containing the vector of dates (\emph{POSIXt}) and the time series of numeric values list perturbed}
+
+\item{i}{[integer] of the indices to subset a time series or [character] names of the elements to extract}
+
+\item{which}{(optional) [character] choice of plots (e.g. \code{"Precip"} or \code{"PotEvap"}, default = \code{"all"})}
+
+\item{main}{(optional) [character] an overall title for the plot (see \code{\link[graphics]{title}})}
+
+\item{ColPrecip, ColPotEvap}{(optional) [character] color to be used for perturbed precipitation and perturbed potential evapotransipration (in any format that \code{\link[grDevices]{col2rgb}} accepts)}
+
+\item{ask}{(optional) [logical] if \code{TRUE}, the user is asked before each plot, see \code{\link[graphics]{par}(ask = .)}}
+
+\item{...}{other parameters to be passed through to plotting functions}
 }
 
-\value{
 
+\value{
 \emph{InputsPert}{[list] object of class \emph{InputsModel} containing the ensembles of perturbed model inputs required to perform data assimilation:}
-
   \tabular{ll}{
      \emph{$DatesR } \tab [POSIXlt] vector of dates
                           \cr\cr
-     \emph{$Precip } \tab [numeric] matrix (dim(Nt, NbMbr)) of ensemblist time series of perturbed total precipitation [mm/d] (with Nt the length of the period; generated only if the precipitation time series is provided as a function argument)
+     \emph{$Precip } \tab [numeric] matrix (dim(NbTime, NbMbr)) of ensemblist time series of perturbed total precipitation [mm/d] (with NbTime the length of the period; generated only if the precipitation time series is provided as a function argument)
                           \cr\cr
-     \emph{$PotEvap} \tab [numeric] matrix (dim(Nt, NbMbr)) of ensemblist time series of perturbed potential evapotranspiration [mm/d] (generated only if the time series of potential evapotranspiration is provided as a function argument)
+     \emph{$PotEvap} \tab [numeric] matrix (dim(NbTime, NbMbr)) of ensemblist time series of perturbed potential evapotranspiration [mm/d] (generated only if the time series of potential evapotranspiration is provided as a function argument)
                           \cr\cr
      \emph{$NbMbr  } \tab [integer] atomic vector of number of ensemble members
   }
 }
 
+
 \description{
 Function which perturbs the model inputs to generate probabilistic meteorological forcings required to perform ensemble-based data assimilation.
 }
 
+
 \details{
 The function generates an ensemble of precipitation or/and potential evapotranspiration time series, from data provided as function argument/s. The mean air temperature required to create CemaNeige inputs is not perturbed.
 \cr\cr
@@ -68,8 +92,13 @@ In order to ensure reproducible results, \emph{Seed} can be set to fix the rando
 For further details, see the references section.
 \cr\cr
 Nota: The function can be applied when using GR4J, GR5J and GR6J models (i.e. daily model time step), with or withouth the CemaNeige module.
+
+On the graphical outputs: \cr
+  - solid line: medians of the input values \cr
+  - polygon: minima and maxima of the input values \cr
 }
 
+
 \examples{
 library(airGRdatassim)
 
@@ -92,21 +121,30 @@ InputsPert <- CreateInputsPert(FUN_MOD = RunModel_GR5J,
                                NbMbr = 100L)
 str(InputsPert)
 
-## results preview on a subset
-par(mfrow = c(2, 1))
-plot(InputsPert[IndRun])
+## results preview
+oldpar <- par(mfrow = c(2, 1))
+plot(InputsPert)
+par(oldpar)
+
+## results preview on a subset on one perturbed variable
+oldpar <- par(mfrow = c(1, 1))
+plot(InputsPert[IndRun], which = "PotEvap")
+par(oldpar)
 }
 
+
 \author{
 Gaia Piazzi, Olivier Delaigue
 }
 
+
 \references{
-- Clark, M. P., Rupp, D. E., Woods, R. A., Zheng, X., Ibbitt, R. P., Slater, A. G. et al. (2008). Hydrological data assimilation with the ensemble Kalman filter: Use of streamflow observations to update states in a distributed hydrological model. Advances in Water Resources, 31(10), 1309-1324, doi: \href{https://www.doi.org/10.1016/j.advwatres.2008.06.005}{10.1016/j.advwatres.2008.06.005}
+- Clark, M. P., Rupp, D. E., Woods, R. A., Zheng, X., Ibbitt, R. P., Slater, A. G. et al. (2008). Hydrological data assimilation with the ensemble Kalman filter: Use of streamflow observations to update states in a distributed hydrological model. Advances in Water Resources, 31(10), 1309-1324, \doi{10.1016/j.advwatres.2008.06.005}
 \cr\cr
-- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (in review). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model. Water Resources Research.
+- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (accepted). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model. Water Resources Research, \doi{10.1029/2020WR028390}.
 }
 
+
 \seealso{
 \code{\link{RunModel_DA}}
 }
diff --git a/man/RunModel_DA.Rd b/man/RunModel_DA.Rd
index cb393e634e918d8926834f45cf50ad93da0f3287..f828d1c42ab35b578ede25d9f489e159598f0300 100644
--- a/man/RunModel_DA.Rd
+++ b/man/RunModel_DA.Rd
@@ -3,9 +3,13 @@
 
 \name{RunModel_DA}
 \alias{RunModel_DA}
+\alias{[.OutputsModelDA}
+\alias{plot.OutputsModelDA}
+
 
 \title{Run discharge simulations with ensemble-based data assimilation}
 
+
 \usage{
 RunModel_DA(InputsModel, InputsPert = NULL, Qobs = NULL,
             IndRun,
@@ -13,20 +17,23 @@ RunModel_DA(InputsModel, InputsPert = NULL, Qobs = NULL,
             DaMethod = c("EnKF", "PF", "none"), NbMbr = NULL,
             StateEnKF = NULL, StatePert = NULL,
             Seed = NULL)
-}
 
-\arguments{
+\method{[}{OutputsModelDA}(x, i)
+
+\method{plot}{OutputsModelDA}(x, Qobs = NULL, main = NULL,
+     ColSim = "orangered", ColObs = par("fg"), ...)
+}
 
 
+\arguments{
 \item{InputsModel}{[list] object of class \emph{InputsModel} containing the data required to run the model  (see \code{\link[airGR]{CreateInputsModel}} for details)}
 
 \item{InputsPert}{(optional) [list] object of class \emph{InputsModel} containing the ensembles of perturbed data required to evaluate the model probabilistic outputs, if the uncertainty in meteorological forcings is taken into account. The variables data proposed in \emph{InputsPert} will erase the variables data given in \emph{InputsModel}}
 
-\item{Qobs}{(optional) [numeric] time series of observed discharges [mm/d]}
+\item{Qobs}{(optional) [numeric] time series of observed flow [mm/d]}
 
 \item{IndRun}{[numeric] index of period to be used for the model run [-]}
 
-
 \item{FUN_MOD}{[function] daily hydrological model functions (GR4J, GR5J or GR6J; e.g. \code{\link[airGR]{RunModel_GR5J}}, \code{\link[airGR]{RunModel_CemaNeigeGR5J}})}
 
 \item{Param}{[numeric] vector of model parameters (number of parameters depends on the used hydrological model, from 4 parameters in GR4J up to 10 parameters in GR6J with CemaNeige)}
@@ -58,38 +65,53 @@ RunModel_DA(InputsModel, InputsPert = NULL, Qobs = NULL,
 
 \item{Seed}{(optional) [numeric] seed of random number generator}
 
+\item{x}{[OutputsModelDA] containing the vector of dates (\emph{POSIXt}) and the time series of numeric values DA model outputs}
+
+\item{i}{[integer] of the indices to subset a time series or [character] names of the elements to extract}
+
+\item{main}{(optional) [character] an overall title for the plot (see \code{\link[graphics]{title}})}
+
+\item{ColSim, ColObs}{(optional) [character] color to be used for simulated flow and observed flow (in any format that \code{\link[grDevices]{col2rgb}} accepts)}
+
+\item{...}{other parameters to be passed through to plotting functions}
 }
 
+
 \value{
 [list] runModel_DA function provides a list containing the outputs organised as follows:
-
   \tabular{ll}{
     \emph{$DatesR     } \tab [POSIXlt] series of dates (length of IndRun)
             \cr
             \cr
-    \emph{$QsimEns    } \tab [numeric] matrix (dim(NbMbr, Nt)) of ensemble discharge simulations
+    \emph{$QsimEns    } \tab [numeric] matrix (dim(NbTime, NbMbr)) of ensemble discharge simulations
             \cr
             \cr
-    \emph{$EnsStateBkg} \tab [numeric] array (dim(Nstate, NbMbr, Nt)) of ensemble values of background model states (before the filter update)
+    \emph{$EnsStateBkg} \tab [numeric] array (dim(NbTime, NbMbr, Nstate)) of ensemble values of background model states (before the filter update)
             \cr
             \cr
-    \emph{$EnsStateA  } \tab [numeric] array (dim(Nstate, NbMbr, Nt)) of ensemble values of analysis model states (after the filter update)
+    \emph{$EnsStateA  } \tab [numeric] array (dim(NbTime, NbMbr, Nstate)) of ensemble values of analysis model states (after the filter update)
             \cr
             \cr
-    \emph{$NbMbr      } \tab [integer] atomic vector of number of ensemble members
+    \emph{$NbTime     } \tab [integer] atomic vector of length of IndRun
             \cr
             \cr
-    \emph{$Nt         } \tab [integer] atomic vector of length of IndRun
+    \emph{$NbMbr      } \tab [integer] atomic vector of number of ensemble members
             \cr
             \cr
     \emph{$NbState    } \tab [integer] atomic vector of number of of model states
          }
+
+On the graphical outputs: \cr
+  - solid line: medians of the input values \cr
+  - polygon: minima and maxima of the input values \cr
 }
 
+
 \description{
-Function which performs discharge ensemble simulations with the assimilation of observed discharges through the Ensemble Kalman filter (EnKF) or the Particle filter (PF) schemes. More information about the efficiency of these data assimilation schemes with GR4J can be found in Piazzi et al. (in review).
+Function which performs discharge ensemble simulations with the assimilation of observed discharges through the Ensemble Kalman filter (EnKF) or the Particle filter (PF) schemes. More information about the efficiency of these data assimilation schemes with GR4J can be found in Piazzi et al. (accepted).
 }
 
+
 \details{
 Discharge observations are sequentially assimilated at each time step (if \emph{DaMethod} != "none") via the EnKF or PF schemes. Because of the sequential approach, the analysis states resulting from the assimilation procedure at time step \emph{t} are used as initial states at the following prediction time step \emph{t + 1}.
 \cr\cr
@@ -106,6 +128,7 @@ For further details and guidelines on the choice of the DA technique, see the re
 Nota: The function can be applied when using GR4J, GR5J and GR6J models (i.e. daily model time step), with or withouth the CemaNeige module (see \code{\link[airGR]{airGR}} package).
 }
 
+
 \examples{
 library(airGRdatassim)
 
@@ -144,18 +167,23 @@ OutputsModelDA <- RunModel_DA(InputsModel = InputsModel,
 
 ## results preview
 plot(OutputsModelDA, Qobs = BasinObs$Qmm[IndRun])
+
+## results preview on a subset
+plot(OutputsModelDA[1:10], Qobs = BasinObs$Qmm[IndRun][1:10])
 }
 
+
 \author{
 Gaia Piazzi, Olivier Delaigue
 }
 
+
 \references{
-- Clark, M. P., Rupp, D. E., Woods, R. A., Zheng, X., Ibbitt, R. P., Slater, A. G. et al. (2008). Hydrological data assimilation with the ensemble Kalman filter: Use of streamflow observations to update states in a distributed hydrological model. Advances in Water Resources, 31(10), 1309-1324, doi: \href{https://www.doi.org/10.1016/j.advwatres.2008.06.005}{10.1016/j.advwatres.2008.06.005}
+- Clark, M. P., Rupp, D. E., Woods, R. A., Zheng, X., Ibbitt, R. P., Slater, A. G. et al. (2008). Hydrological data assimilation with the ensemble Kalman filter: Use of streamflow observations to update states in a distributed hydrological model. Advances in Water Resources, 31(10), 1309-1324, \doi{10.1016/j.advwatres.2008.06.005}
 \cr\cr
-- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (in review). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model. Water Resources Research.
+- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (accepted). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model. Water Resources Research, \doi{10.1029/2020WR028390}.
 \cr\cr
-- Salamon, P. and Feyen, L. (2009). Assessing parameter, precipitation, and predictive uncertainty in a distributed hydrological model using sequential data assimilation with the particle filter. Journal of Hydrology, 376(3-4), 428-442, doi: \href{https://www.doi.org/10.1016/j.jhydrol.2009.07.051}{10.1016/j.jhydrol.2009.07.051}
+- Salamon, P. and Feyen, L. (2009). Assessing parameter, precipitation, and predictive uncertainty in a distributed hydrological model using sequential data assimilation with the particle filter. Journal of Hydrology, 376(3-4), 428-442, \doi{10.1016/j.jhydrol.2009.07.051}
 }
 
 
diff --git a/man/airGRdatassim.Rd b/man/airGRdatassim.Rd
index b69c7e5dafafc1b37416d1bdf3d92003353b1983..e65eae9851f6bdd55851ad5a40fc860011a53fdd 100644
--- a/man/airGRdatassim.Rd
+++ b/man/airGRdatassim.Rd
@@ -2,9 +2,9 @@
 \alias{airGRdatassim}
 \docType{package}
 \encoding{UTF-8}
-\title{Suite of tools to perform ensemble-based discharge assimilation with GR Hydrological Models}
+\title{\packageTitle{airGRdatassim}}
 \description{
-airGRdatassim is a package based on the airGR hydrological modelling package. It provides the tools to assimilate observed discharges in the GR daily hydrological models (GR4J, GR5J and GR6J, with and without the CemaNeige snow model). The package is developed at INRAE-Antony (\href{https://webgr.inrae.fr/en/}{Catchment Hydrology research group} of the HYCAR Research Unit, France). More information about the efficiency of these data assimilation schemes with GR4J can be found in Piazzi et al. (in review). \cr\cr
+\packageDescription{airGRdatassim} airGRdatassim currently runs on daily hydrological models (GR4J, GR5J and GR6J, with and without the CemaNeige snow model). The package is developed at INRAE-Antony (\href{https://webgr.inrae.fr/en/home/}{Catchment Hydrology research group} of the HYCAR Research Unit, France). More information about the efficiency of these data assimilation schemes with GR5J can be found in Piazzi et al. (accepted). \cr\cr
 
 
 ## --- Functions and objects
@@ -22,7 +22,7 @@ Consistently with the airGR package, both structure and class of function argume
 
 ## --- Models
 
-DA schemes are designed to be coupled with GR daily hydrological models, which are implemented in the airGR package. These models can be called within the airGRdatassim package using the following airGR functions (see the \code{\link[airGR]{airGR}} manual to get the references of the GR models): \cr
+DA schemes are designed to be coupled with GR daily hydrological models, which are implemented in the airGR package. These models can be called within the airGRdatassim package using the following airGR functions (use the command `?airGR` to get the references of the GR models): \cr
   - \code{\link[airGR]{RunModel_GR4J}}: four-parameter daily lumped hydrological model \cr
   - \code{\link[airGR]{RunModel_GR5J}}: five-parameter daily lumped hydrological model \cr
   - \code{\link[airGR]{RunModel_GR6J}}: six-parameter daily lumped hydrological model  \cr
@@ -46,7 +46,7 @@ For more information and to get started with the package, you can refer to the v
 
 ## --- References
 
-- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (in review). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model. Water Resources Research. \cr
+- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (accepted). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model. Water Resources Research, \doi{10.1029/2020WR028390}.
 
 }
 
diff --git a/man/plot.Rd b/man/plot.Rd
deleted file mode 100644
index 42f542fbb7637124f53f9831073fe9f14a8ef565..0000000000000000000000000000000000000000
--- a/man/plot.Rd
+++ /dev/null
@@ -1,54 +0,0 @@
-\encoding{UTF-8}
-
-
-\name{plot}
-\alias{plot.InputsPert}
-\alias{plot.OutputsModelDA}
-\alias{plot}
-
-
-\title{Default preview of perturbed inputs and DA model outputs}
-
-
-\description{
-Function which creates time series plots for perturbed inputs and DA model outputs.
-}
-
-
-\usage{
-\method{plot}{InputsPert}(x, which = "all", main = NULL,
-     ColPrecip = "royalblue", ColPotEvap = "green3", ...)
-\method{plot}{OutputsModelDA}(x, Qobs = NULL, main = NULL,
-     ColSim = "orangered", ColObs = par("fg"), ...)
-}
-
-
-\arguments{
-\item{x}{[InputsPert] or [OutputsModelDA] containing the vector of dates (\emph{POSIXt}) and the time series of numeric values list perturbed inputs and DA model outputs (see \code{\link{CreateInputsPert}} and \code{\link{RunModel_DA}})}
-
-\item{Qobs}{(optional) [numeric] time series of observed flow (for the same time steps than simulated) [mm/time step]}
-
-\item{which}{(optional) [character] choice of plots (e.g. \code{"Precip"} or \code{"PotEvap"}, default = \code{"all"})}
-
-\item{main}{(optional) [character] an overall title for the plot (see \code{\link[graphics]{title}})}
-
-\item{ColPrecip}{(optional) [character] color to be used for perturbed precipitation (in any format that \code{\link[grDevices]{col2rgb}} accepts)}
-
-\item{ColPotEvap}{(optional) [character] color to be used for perturbed potential evapotransipration (in any format that \code{\link[grDevices]{col2rgb}} accepts)}
-
-\item{ColSim}{(optional) [character] color to be used for simulated flow (in any format that \code{\link[grDevices]{col2rgb}} accepts)}
-
-\item{ColObs}{(optional) [character] color to be used for observed flow}
-
-\item{...}{other parameters to be passed through to plotting functions}
-}
-
-
-\author{
-Olivier Delaigue, Gaia Piazzi
-}
-
-
-\seealso{
-\code{\link{CreateInputsPert}}, \code{\link{RunModel_DA}}
-}
diff --git a/vignettes/get_started.Rmd b/vignettes/get_started.Rmd
index dcdc29e22e8efd7f1c2fe28e849b37f263ac6488..8ad58e5a8ec34632b5fb3c195770285233a06ccf 100644
--- a/vignettes/get_started.Rmd
+++ b/vignettes/get_started.Rmd
@@ -18,7 +18,7 @@ airGR_bib <- toBibtex(citation("airGR"))
 airGR_bib <- gsub("@Article\\{", "@Article{airGR2017", airGR_bib)
 airGR_bib <- gsub("@Manual\\{", "@Manual{R-airGR", airGR_bib)
 airGRdatassim_bib <- toBibtex(citation("airGRdatassim")) #toBibtex(citation("airGRdatassim")[1])
-airGRdatassim_bib <- gsub("@Article\\{", "@Article{piazzi_sequential_2020", airGRdatassim_bib)
+airGRdatassim_bib <- gsub("@Article\\{", "@Article{piazzi_sequential_2021", airGRdatassim_bib)
 airGRdatassim_bib <- gsub("@Manual\\{", "@Manual{R-airGRdatassim", airGRdatassim_bib)
 options(encoding = "UTF-8")
 writeLines(text = airGR_bib, con = "airgr_ref.bib")
@@ -35,7 +35,7 @@ Discharge simulations are affected by several sources of uncertainty (e.g. error
 
 Data assimilation (DA) allows to combine measurements and model simulations under the assumption that both supply useful information to obtain the most likely estimate of the system state. Among the sequential methods, the ensemble-based techniques allow to explicitly handle different sources of uncertainty and to quantify the unknown errors of both model states and observations.
 
-After the study of @piazzi_sequential_2020, airGRdatassim package [@R-airGRdatassim] provides functions allowing to perform DA-based ensemble discharge simulations. airGRdatassim is a package based on the [airGR](https://CRAN.R-project.org/package=airGR) hydrological modeling package [@airGR2017; @R-airGR]. This package is designed to flexibly use the Ensemble Kalman filter (EnKF) or Particle filter (PF) technique to assimilate daily discharge observations into GR4J, GR5J, and GR6J models  (see the [airGR manual](https://cran.r-project.org/web/packages/airGR/airGR.pdf#page.4) to get the references of the GR models). 
+After the study of @piazzi_sequential_2021, airGRdatassim package [@R-airGRdatassim] provides functions allowing to perform DA-based ensemble discharge simulations. airGRdatassim is a package based on the [airGR](https://CRAN.R-project.org/package=airGR) hydrological modeling package [@airGR2017; @R-airGR]. This package is designed to flexibly use the Ensemble Kalman filter (EnKF) or Particle filter (PF) technique to assimilate daily discharge observations into GR4J, GR5J, and GR6J models (use the command `?airGR` to get the references of the GR models). 
 
 The uncertainty in both the meteorological forcings (i.e. precipitation and potential evapotranspiration) and model states (i.e. production store level, routing store level and unit hydrographs) can be easily taken into account by enabling specific perturbation procedures.
 
@@ -81,7 +81,7 @@ If `DaMethod = "EnKF"`, the EnKF performs the update of the state variables spec
 
 If `DaMethod = "none"`, discharge assimilation is not performed (i.e., open-loop simulation).
 
-Among the model states, the update of the routing store level ensures the most benefit from the assimilation of observed discharges. Conversely, the update of the unit hydrograph only has a slight impact on the accuracy of the DA-based discharge simulations [@piazzi_sequential_2020].
+Among the model states, the update of the routing store level ensures the most benefit from the assimilation of observed discharges. Conversely, the update of the unit hydrograph only has a slight impact on the accuracy of the DA-based discharge simulations [@piazzi_sequential_2021].
 
 
 ## Model uncertainties
@@ -95,7 +95,7 @@ The perturbation of model states is performed through a normally-distributed nul
 When considering the uncertainty in the forcings, the `CreateInputsPert()` function generates an ensemble for precipitation or/and potential evapotranspiration, provided as function argument/s.
 It is noteworthy that a significant reduction in the ensemble spread may occur especially in no-rain periods, when accounting only for meteorological uncertainty.
 
-Even though a more comprehensive representation of the state uncertainty can prevent the ensemble shrinkage, it has a different impact on the performance of the PF and the EnKF schemes [@piazzi_sequential_2020]. Indeed, while a higher model error covariance can lead to an overweighting of observations in the EnKF-based analysis procedure, a larger spread of the ensemble simulations allows for a more straightforward weighting and thus a more efficient resampling of particles in the PF scheme. 
+Even though a more comprehensive representation of the state uncertainty can prevent the ensemble shrinkage, it has a different impact on the performance of the PF and the EnKF schemes [@piazzi_sequential_2021]. Indeed, while a higher model error covariance can lead to an overweighting of observations in the EnKF-based analysis procedure, a larger spread of the ensemble simulations allows for a more straightforward weighting and thus a more efficient resampling of particles in the PF scheme. 
 
 
 ## Definition of DA settings
@@ -251,14 +251,14 @@ CritOL <- ErrorCrit(InputsCrit = InputsCritMulti, OutputsModel = OutputsModel_OL
 ```{r, warning=FALSE, eval=TRUE, message=FALSE}
 # EnKF run
 OutputsModel_EnKF <- OutputsModel_OL
-OutputsModel_EnKF$Qsim <- colMeans(ResEnKF$QsimEns)
+OutputsModel_EnKF$Qsim <- rowMeans(ResEnKF$QsimEns)
 CritEnKF <- ErrorCrit(InputsCrit = InputsCritMulti, OutputsModel = OutputsModel_EnKF)
 ```
 
 ```{r, warning=FALSE, eval=TRUE, message=FALSE}
 # PF run
 OutputsModel_PF <- OutputsModel_OL
-OutputsModel_PF$Qsim <- colMeans(ResPF$QsimEns)
+OutputsModel_PF$Qsim <- rowMeans(ResPF$QsimEns)
 CritPF <- ErrorCrit(InputsCrit = InputsCritMulti, OutputsModel = OutputsModel_PF)
 ```
 
@@ -279,24 +279,25 @@ CritVal
 
 
 ```{r, fig.width=7, fig.height=5, warning=FALSE, echo=FALSE}
-ylimMin <- min(colMeans(ResEnKF$QsimEns), colMeans(ResPF$QsimEns), OutputsModel_OL$Qsim, Qobs)
-ylimMax <- max(colMeans(ResEnKF$QsimEns), colMeans(ResPF$QsimEns), OutputsModel_OL$Qsim, Qobs)
+ylimMin <- min(rowMeans(ResEnKF$QsimEns), rowMeans(ResPF$QsimEns), OutputsModel_OL$Qsim, Qobs)
+ylimMax <- max(rowMeans(ResEnKF$QsimEns), rowMeans(ResPF$QsimEns), OutputsModel_OL$Qsim, Qobs)
 
-par(mar = c(5, 4, 2, 2))
+oldpar <- par(mar = c(5, 4, 2, 2))
 matplot(x = BasinObs$DatesR[IndRun],
         y = cbind(OutputsModel_OL$Qsim,
-                  colMeans(ResEnKF$QsimEns),
-                  colMeans(ResPF$QsimEns),
+                  rowMeans(ResEnKF$QsimEns),
+                  rowMeans(ResPF$QsimEns),
                   Qobs),
         type = "l", col = c("red", "green", "blue", "black"),
         lty = c(1, 1, 1, 3), lwd = c(1.5, 1.5, 1.5, 1.5),
-        xlab = "Time", ylab = "Discharge [mm/day]",
+        xlab = "Time", ylab = "Discharge [mm/d]",
         xaxt = "n"
 )
 axis.POSIXct(side = 1, x = BasinObs$DatesR[IndRun])
 legend("topleft", legend = c("OL", "EnKF", "PF", "Obs"),
        col = c("red", "green", "blue", "black"),
        lty = c(1, 1, 1, 3), lwd = rep(1.5, 4))
+par(oldpar)
 ```