diff --git a/DESCRIPTION b/DESCRIPTION
index 09f5f69e02807a774e4ee3ee1a2b90b4e71834fe..481de5209390f56875800948d061aefb75420269 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.6.8.44
+Version: 1.6.9.9
 Date: 2021-01-08
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
diff --git a/NAMESPACE b/NAMESPACE
index 5d211da10e84c0fd65df7e835c17e466ee96d4a7..bf71216c1b7e82680ccdaf969983b105c814200d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,6 +8,8 @@ useDynLib(airGR, .registration = TRUE)
 #####################################
 ##            S3 methods           ##
 #####################################
+S3method('[', InputsModel)
+S3method('[', OutputsModel)
 S3method(plot, OutputsModel)
 S3method(SeriesAggreg, data.frame)
 S3method(SeriesAggreg, list)
diff --git a/NEWS.md b/NEWS.md
index b9d37506e33fb29edbcce3663c249500cf827147..145e53bd95478c7dcf339e76889f0ca70b2a48e6 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 
 
 
-### 1.6.8.44 Release Notes (2021-01-08)
+### 1.6.9.9 Release Notes (2021-01-08)
 
 #### New features
 
@@ -12,6 +12,7 @@
 - `PE_Oudin()` now presents a `RunFortran` argument to run the code in Fortran or in R. The Fortran mode is the fastest. ([#62](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/62))
 - Added `RunModel_Lag()` which allows to perform a single run for the Lag model over the test period in order to run semi-distributed GR models. ([#34](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/34))
 - Added the 'sd_model' vignette to explain how to manage the use of semi-distributed GR models. ([#34](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/34))
+- Added `[` S3 method for `InputsModel` and, `OutputsModel` class objects in order to extract subsets of them. ([#67](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/67))
 
 
 #### Deprecated and defunct
diff --git a/R/SeriesAggreg.list.R b/R/SeriesAggreg.list.R
index 8834347dbc5ac40fa35a5c19b074699173c34f15..ab7df82d2cfe03ed205d8a47c31231077bb11be8 100644
--- a/R/SeriesAggreg.list.R
+++ b/R/SeriesAggreg.list.R
@@ -6,6 +6,10 @@ SeriesAggreg.list <- function(x,
                               except = NULL,
                               recursive = TRUE,
                               ...) {
+
+  classIni <- class(x)
+  class(x) <- "list" # in order to avoid the use of '['.InputsModel' or '['.OutputsModel' when x[i] is used
+
   if (missing(Format)) {
     Format <- .GetSeriesAggregFormat(NewTimeFormat)
   } else if (!is.null(NewTimeFormat)) {
@@ -143,7 +147,7 @@ SeriesAggreg.list <- function(x,
 
     class(res) <- gsub("hourly|daily|monthly|yearly",
                        .GetSeriesAggregClass(Format),
-                       class(x))
+                       classIni)
 
     return(res)
 
diff --git a/R/Utils.R b/R/Utils.R
index d1cde5eb0dc10d9f3f686ddf03dab9c24671d850..fe496be1d0a2e23f92ccca6b1813ed18c5e63a07 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -11,6 +11,8 @@
 #   }
 # }
 
+
+
 ## =================================================================================
 ## function to manage Fortran outputs
 ## =================================================================================
@@ -74,3 +76,93 @@
   res <- list(GR = outGR, CN = outCN)
 
 }
+
+
+
+## =================================================================================
+## functions to extract parts of InputsModel or OutputsModel objects
+## =================================================================================
+
+## InputsModel
+
+.ExtractInputsModel <- function(x, i) {
+  res <- lapply(x, function(x) {
+    if (is.matrix(x)) {
+      res0 <- x[i, , drop = FALSE]
+    }
+    if (is.vector(x) | inherits(x, "POSIXt")) {
+      res0 <- x[i]
+    }
+    if (is.list(x) & !inherits(x, "POSIXt")) {
+      if (inherits(x, "OutputsModel")) {
+        res0 <- .ExtractOutputsModel(x = x, i = i)
+      } else {
+        res0 <- .ExtractInputsModel(x = x, i = i)
+      }
+    }
+    return(res0)
+  })
+  if (!is.null(x$ZLayers)) {
+    res$ZLayers <- x$ZLayers
+  }
+  if (inherits(x, "SD")) {
+    res$LengthHydro <- x$LengthHydro
+    res$BasinAreas  <- x$BasinAreas
+  }
+  class(res) <- class(x)
+  res
+}
+
+'[.InputsModel' <- function(x, i) {
+  if (!inherits(x, "InputsModel")) {
+    stop("'x' must be of class 'InputsModel'")
+  }
+  if (is.factor(i)) {
+    i <- as.character(i)
+  }
+  if (is.numeric(i)) {
+    .ExtractInputsModel(x, i)
+  } else {
+    NextMethod()
+  }
+}
+
+
+## InputsModel
+
+.ExtractOutputsModel <- function(x, i) {
+  res <- lapply(x, function(x) {
+    if (is.matrix(x)  && length(dim(x)) == 2L) {
+      res0 <- x[i, ]
+    }
+    if (is.array(x) && length(dim(x)) == 3L) {
+      res0 <- x[i, , ]
+    }
+    if (is.vector(x) | inherits(x, "POSIXt")) {
+      res0 <- x[i]
+    }
+    if (is.list(x) & !inherits(x, "POSIXt")) {
+      res0 <- .ExtractOutputsModel(x = x, i = i)
+    }
+    return(res0)
+  })
+  if (!is.null(x$StateEnd)) {
+    res$StateEnd <- x$StateEnd
+  }
+  class(res) <- class(x)
+  res
+}
+
+'[.OutputsModel' <- function(x, i) {
+  if (!inherits(x, "OutputsModel")) {
+    stop("'x' must be of class 'OutputsModel'")
+  }
+  if (is.factor(i)) {
+    i <- as.character(i)
+  }
+  if (is.numeric(i)) {
+    .ExtractOutputsModel(x, i)
+  } else {
+    NextMethod()
+  }
+}
diff --git a/man/CreateInputsModel.Rd b/man/CreateInputsModel.Rd
index 3716e10108f338ae2be103361067424244ed5cef..b94b675d2d42c6f3a52ab27c6288339f5e923fed 100644
--- a/man/CreateInputsModel.Rd
+++ b/man/CreateInputsModel.Rd
@@ -3,6 +3,7 @@
 
 \name{CreateInputsModel}
 \alias{CreateInputsModel}
+\alias{[.InputsModel}
 
 
 \title{Creation of the InputsModel object required to the RunModel functions}
@@ -19,6 +20,8 @@ CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale = TRUE, PotEvap = NULL,
                   ZInputs = NULL, HypsoData = NULL, NLayers = 5,
                   Qupstream = NULL, LengthHydro = NULL, BasinAreas = NULL,
                   verbose = TRUE)
+
+\method{[}{InputsModel}(x, i)
 }
 
 
@@ -53,6 +56,9 @@ CreateInputsModel(FUN_MOD, DatesR, Precip, PrecipScale = TRUE, PotEvap = NULL,
 
 \item{BasinAreas}{(optional) [numeric] real giving the area of each upstream sub-catchment [km2] and the area of the downstream sub-catchment in the last item, required to create the SD model inputs (see details)}
 
+\item{x}{[InputsModel] object of class InputsModel}
+
+\item{i}{[integer] of the indices to subset a time series or [character] names of the elements to extract}
 }
 
 
diff --git a/man/RunModel.Rd b/man/RunModel.Rd
index ceae5f1dc480d3ef4516b31ad145db9d769749a0..efead0f00cb88c0444224b2586647b6dbb4da1ad 100644
--- a/man/RunModel.Rd
+++ b/man/RunModel.Rd
@@ -3,6 +3,7 @@
 
 \name{RunModel}
 \alias{RunModel}
+\alias{[.OutputsModel}
 
 
 \title{Run with the provided hydrological model function}
@@ -15,6 +16,8 @@ Function which performs a single model run with the provided function over the s
 
 \usage{
 RunModel(InputsModel, RunOptions, Param, FUN_MOD)
+
+\method{[}{OutputsModel}(x, i)
 }
 
 
@@ -26,6 +29,10 @@ RunModel(InputsModel, RunOptions, Param, FUN_MOD)
 \item{Param}{[numeric] vector of model parameters (See details for SD lag model)}
 
 \item{FUN_MOD}{[function] hydrological model function (e.g. \code{\link{RunModel_GR4J}}, \code{\link{RunModel_CemaNeigeGR4J}})}
+
+\item{x}{[InputsModel] object of class InputsModel}
+
+\item{i}{[integer] of the indices to subset a time series or [character] names of the elements to extract}
 }