diff --git a/.Rbuildignore b/.Rbuildignore
index 047e4e576fe142755c31c6616bb60f7d15e1c37b..a85fa700638207047334f50b0c444225db5e9cd9 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,4 +1,4 @@
-^griwrm\.Rproj$
+^airGRiwrm\.Rproj$
 ^\.Rproj\.user$
 ^LICENSE\.md$
 ^\.gitlab-ci\.yml$
diff --git a/DESCRIPTION b/DESCRIPTION
index 250528a8613b02a1edb250faf6f35cbdee125d66..c388a6073f54d7657bdbf661caa5ea67e0282e15 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,21 +1,23 @@
 Package: airGRiwrm
-Title: airGR Integrated Water Resource Management
-Version: 0.4.0
+Title: 'airGR' Integrated Water Resource Management
+Version: 0.5.0
 Authors@R: c(
     person("David", "Dorchies", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-6595-7984"), email = "david.dorchies@inrae.fr"),
     person("Olivier", "Delaigue", role = c("ctb"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
     person("Guillaume", "Thirel", role = c("ctb"), comment = c(ORCID = "0000-0002-1444-1830"))
     )
-Description: This R package aims to model water basin using 'airGR' based semi-distributive model with the integration of human infrastructures and their management.
+Description: Semi-distributive Precipitation-Runoff Modelling based on 'airGR' package models integrating human infrastructures and their managements.
 License: AGPL-3
 Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.1.0
+RoxygenNote: 7.1.1
 Imports:
     dplyr,
     utils,
-    airGR (>= 1.6.1.11)
+    airGR (>= 1.6.1.11),
+    grDevices,
+    graphics
 Suggests:
     knitr,
     rmarkdown,
diff --git a/R/Calibration.GRiwrmInputsModel.R b/R/Calibration.GRiwrmInputsModel.R
index 3dfbe106e0bbf59801dc637b59c2c00e775d34f7..3b62b7761caa3166ec8000ece5591de648991702 100644
--- a/R/Calibration.GRiwrmInputsModel.R
+++ b/R/Calibration.GRiwrmInputsModel.R
@@ -30,7 +30,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
       IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]]$IndPeriod_Run, OutputsModel)
     }
 
-    OutputsCalib[[IM$id]] <- Calibration.InputsModel(
+    OutputsCalib[[IM$id]] <- Calibration(
       InputsModel = IM,
       RunOptions = RunOptions[[IM$id]],
       InputsCrit = InputsCrit[[IM$id]],
@@ -41,7 +41,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
     if(useUpstreamQsim) {
       # Run the model for the sub-basin
       OutputsModel[[IM$id]] <- RunModel(
-        InputsModel = IM,
+        x = IM,
         RunOptions = RunOptions[[IM$id]],
         Param = OutputsCalib[[IM$id]]$ParamFinalR
       )
diff --git a/R/CreateInputsCrit.GRiwrmInputsModel.R b/R/CreateInputsCrit.GRiwrmInputsModel.R
index 56af942644ae37dcfba799109b45c1215fe1395c..889600818f88276667a706f849023ef5b7f5fd4e 100644
--- a/R/CreateInputsCrit.GRiwrmInputsModel.R
+++ b/R/CreateInputsCrit.GRiwrmInputsModel.R
@@ -1,11 +1,14 @@
 #' Create \emph{GRiwrmInputsCrit} object for **airGRiwrm**.
+#'
+#' This function does the same operations as [airGR::CreateInputsCrit] for all sub-basins of the GRiwrm model.
+#'
 #' @param InputsModel  object of class \emph{GRiwrmInputsModel}, see \code{\link{CreateInputsModel.GRiwrm}} for details.
 #' @param FUN_CRIT \[function (atomic or list)\] error criterion function (e.g. \code{\link[airGR]{ErrorCrit_RMSE}}, \code{\link[airGR]{ErrorCrit_NSE}})
 #' @param RunOptions object of class \emph{GRiwrmRunOptions}, see \code{[CreateRunOptions.GRiwrm]} for details.
 #' @param Qobs matrix or data frame containing observed flows. Column names correspond to nodes ID
-#' @param ... further arguments passed to \code{\link[airGR]{CreateInputsCrit}}.
+#' @param ... further arguments passed to [airGR::CreateInputsCrit].
 #'
-#' @return Object of class \emph{GRiwrmInputsCrit}
+#' @return Object of class \emph{GRiwrmInputsCrit} which is a list of `airGR::InputsCrit` objects (See [airGR::CreateInputsCrit])
 #' @export
 CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel,
                                                FUN_CRIT = airGR::ErrorCrit_NSE,
diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index 2b541e7975aa5c2b67da4adb473e65cb5d4618e5..78287a8b21d13919f2d4fa036a9ffde62199400d 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -11,8 +11,8 @@
 #' @export
 CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
 
-  InputsModel <- CreateEmptyGRiwrmInputsModel()
-  Qobs[is.na(Qobs)] <- -99 # airGRCreateInputsModel doesn't accept NA values
+  InputsModel <- CreateEmptyGRiwrmInputsModel(x)
+  Qobs[is.na(Qobs)] <- -99 # airGR::CreateInputsModel doesn't accept NA values
 
   for(id in getNodeRanking(x)) {
     message("CreateInputsModel.GRiwrm: Treating sub-basin ", id, "...")
@@ -20,16 +20,20 @@ CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
       id, x, DatesR,Precip[,id], PotEvap[,id], Qobs, ...
     )
   }
+  attr(InputsModel, "TimeStep") <- getModelTimeStep(InputsModel)
   return(InputsModel)
 }
 
 
 #' Create an empty InputsModel object for **airGRiwrm** nodes
 #'
+#' @param griwrm a `GRiwrm` object (See [GRiwrm])
+#'
 #' @return \emph{GRiwrmInputsModel} empty object
-CreateEmptyGRiwrmInputsModel <- function() {
+CreateEmptyGRiwrmInputsModel <- function(griwrm) {
   InputsModel <- list()
-  class(InputsModel) <- append(class(InputsModel), "GRiwrmInputsModel")
+  class(InputsModel) <- c("GRiwrmInputsModel", class(InputsModel))
+  attr(InputsModel, "GRiwrm") <- griwrm
   return(InputsModel)
 }
 
@@ -79,6 +83,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs
 
   # Add Identifiers of connected nodes in order to be able to update them with simulated flows
   InputsModel$id <- id
+  InputsModel$down <- node$down
   if(length(UpstreamNodes) > 0) {
     InputsModel$UpstreamNodes <- UpstreamNodes
     InputsModel$UpstreamIsRunoff <- !is.na(griwrm$model[match(UpstreamNodes, griwrm$id)])
@@ -88,3 +93,28 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs
 
   return(InputsModel)
 }
+
+
+#' Check time steps of the model of all the nodes and return the time step in seconds
+#'
+#' This function is called inside [CreateInputsModel.GRiwrm] for defining the time step of the big model.
+#'
+#' @param InputsModel a `GRiwrmInputsModel`
+#'
+#' @return A [numeric] representing the time step in seconds
+#'
+getModelTimeStep <- function(InputsModel) {
+  TS <- sapply(InputsModel, function(x) {
+    if (inherits(x, "hourly")) {
+      TimeStep <- 60 * 60
+    } else if (inherits(x, "daily")) {
+      TimeStep <- 60 * 60 * 24
+    } else {
+      stop("All models should be at hourly or daily time step")
+    }
+  })
+  if(length(unique(TS)) != 1) {
+    stop("Time steps of the model of all nodes should be identical")
+  }
+  return(unique(TS))
+}
diff --git a/R/CreateInputsModel.R b/R/CreateInputsModel.R
index 618463f62379c4fd3a641bb64a95124700bd8790..6d17ddabaec16154568cc04297fcb6508ade9ad0 100644
--- a/R/CreateInputsModel.R
+++ b/R/CreateInputsModel.R
@@ -1,4 +1,6 @@
-#' Create InputsModel object for either **airGR** or **airGRiwrm**
+#' Generic function for creating `InputsModel` object for either **airGR** or **airGRiwrm**
+#'
+#' See the methods [CreateInputsModel.GRiwrm] for **airGRiwrm** and [CreateInputsModel.default] for **airGR**.
 #'
 #' @param x First parameter determining which InputsModel object is created
 #' @param ... further arguments passed to or from other methods.
diff --git a/R/CreateSupervisor.R b/R/CreateSupervisor.R
new file mode 100644
index 0000000000000000000000000000000000000000..d0be4935cbd7e1d915ac8d2bb5cca3e84f74a5f2
--- /dev/null
+++ b/R/CreateSupervisor.R
@@ -0,0 +1,71 @@
+#' Create a Supervisor for handling regulation in a model
+#'
+#' @param InputsModel `GRiwrmInputsModel` The inputs of the basin model
+#' @param TimeStep [integer] The number of time steps between each supervision
+#'
+#' @return `Supervisor` object
+#' @export
+#'
+#' @examples
+#' data(Severn)
+#' nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+#' nodes$distance_downstream <- nodes$distance_downstream * 1000 # Conversion km -> m
+#' nodes$model <- "RunModel_GR4J"
+#' griwrm <- GRiwrm(nodes,
+#'                  list(id = "gauge_id",
+#'                       down = "downstream_id",
+#'                       length = "distance_downstream"))
+#' BasinsObs <- Severn$BasinsObs
+#' DatesR <- BasinsObs[[1]]$DatesR
+#' PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+#' PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+#' Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
+#' Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+#' PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+#' InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+#' sv <- CreateSupervisor(InputsModel)
+CreateSupervisor <- function(InputsModel, TimeStep = 1L) {
+  if(!inherits(InputsModel, "GRiwrmInputsModel")) {
+    stop("`InputsModel` parameter must of class 'GRiwrmInputsModel' (See ?CreateInputsModel.GRiwrm)")
+  }
+  if(!is.integer(TimeStep)) stop("`TimeStep` parameter must be an integer")
+
+  # Create Supervisor environment from the parent of GlobalEnv
+  e <- new.env(parent = parent.env(globalenv()))
+  class(e) <- c("Supervisor", class(e))
+
+  # Hidden variable to detect which environment it is
+  e$.isSupervisor <- "3FJKmDcJ4snDbVBg"
+
+  # Add pointer to itself in order to assign variable from function environment
+  e$supervisor <- e
+
+  # Copy of InputsModel, griwrm and prepare OutputsModel
+  e$InputsModel <- InputsModel
+  e$griwrm <- attr(InputsModel, "GRiwrm")
+  e$OutputsModel <- list()
+  e$.TimeStep <- TimeStep
+
+  # Controller list
+  e$controllers <- list()
+  class(e$controllers) <- c("Controllers", class(e$controllers))
+
+  # Copy functions to be used enclosed in the Supervisor environment
+  e$CreateController <- CreateController
+  environment(e$CreateController) <- e
+
+  # Time steps handling: these data are provided by RunModel
+  # Index of the current time steps in the modelled time series between 1 and length(RunOptions$Ind_Period)
+  e$ts.index <- NA
+  # Index of the previous time steps in the modelled time series
+  e$ts.previous <- NA
+  # Index of the time step preceding RunOptions$Ind_Period
+  e$ts.index0 <- NA
+  # Date/Time of the current time step (For controller calculations based on date)
+  e$ts.date <- NULL
+
+  # Current Controller ID (Updated in doSupervision)
+  e$controller.id <- NULL
+
+  return(e)
+}
diff --git a/R/DiagramGRiwrm.R b/R/DiagramGRiwrm.R
index 1930c6c20d130bc6deb129b643536afe0678358c..dad08f9bef26ee9e67636e8a951d8901d7d8320f 100644
--- a/R/DiagramGRiwrm.R
+++ b/R/DiagramGRiwrm.R
@@ -23,7 +23,7 @@ DiagramGRiwrm <- function(griwrm, display = TRUE, orientation = "LR") {
   if(Sys.getenv("RSTUDIO") != "1") {
     return()
   }
-  if(!"DiagrammeR" %in% rownames(installed.packages())) {
+  if(!"DiagrammeR" %in% rownames(utils::installed.packages())) {
     stop("The 'DiagrammeR' package should be installed. Type: install.packages('DiagrammeR')")
   }
   g2 <- griwrm[!is.na(griwrm$down),]
diff --git a/R/GRiwrm.R b/R/GRiwrm.R
index 6ee029720809f731753bcffc46c66953b6f8ec04..34f6bcd12f5226a00ff71a0ae234be60878d328f 100644
--- a/R/GRiwrm.R
+++ b/R/GRiwrm.R
@@ -42,7 +42,6 @@ GRiwrm <- function(db,
                         length = "double",
                         model = "character",
                         area = "double"))
-  rownames(db) <- db$id
   class(db) <- c("GRiwrm", class(db))
   db
 }
diff --git a/R/RunModel.GR.R b/R/RunModel.GR.R
new file mode 100644
index 0000000000000000000000000000000000000000..da2bc072a62b7dc6f720c4df406734898db192fa
--- /dev/null
+++ b/R/RunModel.GR.R
@@ -0,0 +1,22 @@
+#' Run rainfall-runoff part of a sub-basin model
+#'
+#' @inherit airGR::RunModel
+#' @param x `InputsModel` used as `InputsModel` parameter for [airGR::RunModel]
+#' @param ... further arguments passed to or from other methods
+#'
+#' @export
+#'
+RunModel.GR <- function(x, RunOptions, Param, ...) {
+
+  if (inherits(x, "SD")) {
+    # Lag model take one parameter at the beginning of the vector
+    iFirstParamRunOffModel <- 2
+  } else {
+    # All parameters
+    iFirstParamRunOffModel <- 1
+  }
+
+  FUN_MOD <- match.fun(x$FUN_MOD)
+  FUN_MOD(x, RunOptions = RunOptions,
+          Param = Param[iFirstParamRunOffModel:length(Param)])
+}
diff --git a/R/RunModel.GRiwrmInputsModel.R b/R/RunModel.GRiwrmInputsModel.R
index 87f00c2144e440e9715abaafe4b925a3e17233f5..c40140600fcddb3144ef7e807a7602ad7d0ecaec 100644
--- a/R/RunModel.GRiwrmInputsModel.R
+++ b/R/RunModel.GRiwrmInputsModel.R
@@ -1,32 +1,34 @@
-#' Title
+#' RunModel function for GRiwrmInputsModel object
 #'
-#' @param InputsModel object of class \emph{GRiwrmInputsModel}, see \code{[CreateInputsModel.GRiwrm]} for details.
+#' @param x object of class \emph{GRiwrmInputsModel}, see \code{[CreateInputsModel.GRiwrm]} for details.
 #' @param RunOptions object of class \emph{GRiwrmRunOptions}, see \code{[CreateRunOptions.GRiwrm]} for details.
 #' @param Param list of parameter. The list item names are the IDs of the sub-basins. Each item is a vector of numerical parameters.
 #' @param ... Mandatory for S3 method signature function compatibility with generic.
 #'
 #' @return \emph{GRiwrmOutputsModel} object which is a list of \emph{OutputsModel} objects (See \code{\link[airGR]{RunModel}}) for each node of the semi-distributed model.
 #' @export
-RunModel.GRiwrmInputsModel <- function(InputsModel, RunOptions, Param, ...) {
+RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
+
+  checkRunModelParameters(x, RunOptions, Param)
 
   OutputsModel <- list()
-  class(OutputsModel) <- append(class(OutputsModel), "GRiwrmOutputsModel")
+  class(OutputsModel) <- c("GRiwrmOutputsModel", class(OutputsModel))
 
-  for(IM in InputsModel) {
-    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", IM$id, "...")
+  for(id in names(x)) {
+    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", x[[id]]$id, "...")
 
-    # Update InputsModel$Qupstream with simulated upstream flows
-    if(any(IM$UpstreamIsRunoff)) {
-      IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]]$IndPeriod_Run, OutputsModel)
+    # Update x[[id]]$Qupstream with simulated upstream flows
+    if(any(x[[id]]$UpstreamIsRunoff)) {
+      x[[id]] <- UpdateQsimUpstream(x[[id]], RunOptions[[id]]$IndPeriod_Run, OutputsModel)
     }
 
     # Run the model for the sub-basin
-    OutputsModel[[IM$id]] <- RunModel(
-      InputsModel = IM,
-      RunOptions = RunOptions[[IM$id]],
-      Param = Param[[IM$id]]
+    OutputsModel[[id]] <- RunModel.InputsModel(
+      x[[id]],
+      RunOptions = RunOptions[[id]],
+      Param = Param[[id]]
     )
-
   }
+  attr(OutputsModel, "Qm3s") <- OutputsModelQsim(x, OutputsModel, RunOptions[[1]]$IndPeriod_Run)
   return(OutputsModel)
 }
diff --git a/R/RunModel.InputsModel.R b/R/RunModel.InputsModel.R
index 83affc57243f387b1818f4a9b0c367634c5eb81d..c2c290aabbb66188eea542c927b8c4682f803d86 100644
--- a/R/RunModel.InputsModel.R
+++ b/R/RunModel.InputsModel.R
@@ -1,11 +1,12 @@
 #' Wrapper for \code{\link[airGR]{RunModel}} for one sub-basin.
 #'
 #' @inherit airGR::RunModel
+#' @param x `InputsModel` used as `InputsModel` parameter for [airGR::RunModel]
 #' @param ... Further arguments for compatibility with S3 method
 #' @export
-RunModel.InputsModel <- function(InputsModel, RunOptions, Param, FUN_MOD = NULL, ...) {
+RunModel.InputsModel <- function(x, RunOptions, Param, FUN_MOD = NULL, ...) {
   if(is.null(FUN_MOD)) {
-    FUN_MOD <- InputsModel$FUN_MOD
+    FUN_MOD <- x$FUN_MOD
   }
-  airGR::RunModel(InputsModel, RunOptions, Param, FUN_MOD)
+  airGR::RunModel(x, RunOptions, Param, FUN_MOD)
 }
diff --git a/R/RunModel.R b/R/RunModel.R
index 48eb869b2e34f923398fe72301ef4b76df2bf5c7..0ac7f23a9c50d35ca70a0877a1f8f1d3693f3727 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -1,10 +1,10 @@
 #' RunModel function for both **airGR** InputsModel and GRiwrmInputsModel object
 #'
-#' @param InputsModel object of class \emph{InputsModel} or \emph{GRiwrmInputsModel}. See \code{\link{CreateInputsModel}} for details
+#' @param x object of class \emph{InputsModel} or \emph{GRiwrmInputsModel}. See \code{\link{CreateInputsModel}} for details
 #' @param ... further arguments passed to or from other methods
 #'
 #' @return Either a [list] of OutputsModel object (for GRiwrmInputsModel) or an OutputsModel object (for InputsModel)
 #' @export
-RunModel <- function(InputsModel, ...) {
-  UseMethod("RunModel", InputsModel)
+RunModel <- function(x, ...) {
+  UseMethod("RunModel", x)
 }
diff --git a/R/RunModel.SD.R b/R/RunModel.SD.R
new file mode 100644
index 0000000000000000000000000000000000000000..284a85a2533ee48b74ce291303bb37f9d23c9661
--- /dev/null
+++ b/R/RunModel.SD.R
@@ -0,0 +1,14 @@
+#' Run SD Model from run-off model outputs
+#'
+#' @inheritParams airGR::RunModel_Lag
+#' @param x `InputsModel` used as `InputsModel` parameter for [airGR::RunModel_Lag]
+#' @param QsimDown a [numeric] corresponding to the runoff of the sub-basin (Typically the `Qsim` outputs of the GR model)
+#' @param ... further arguments passed to or from other methods
+#'
+#' @return `OutputsModel` object. See [airGR::RunModel_Lag]
+#' @export
+#'
+RunModel.SD <- function(x, RunOptions, Param, QsimDown, ...) {
+  x$OutputsModel <- list(Qsim = QsimDown)
+  RunModel_Lag(x, RunOptions = RunOptions, Param = Param[1])
+}
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
new file mode 100644
index 0000000000000000000000000000000000000000..c1b7802c8f66b1cf61f04d087024df29e0522c7b
--- /dev/null
+++ b/R/RunModel.Supervisor.R
@@ -0,0 +1,86 @@
+#' RunModel function for GRiwrmInputsModel object
+#'
+#' @param x object of class `Supervisor`, see [CreateSupervisor] for details.
+#' @param RunOptions object of class \emph{GRiwrmRunOptions}, see \code{[CreateRunOptions.GRiwrm]} for details.
+#' @param Param list of parameter. The list item names are the IDs of the sub-basins. Each item is a vector of numerical parameters.
+#' @param ... Mandatory for S3 method signature function compatibility with generic.
+#'
+#' @return \emph{GRiwrmOutputsModel} object which is a list of \emph{OutputsModel} objects (See \code{\link[airGR]{RunModel}}) for each node of the semi-distributed model.
+#' @export
+RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
+
+  # Time steps handling
+  x$ts.index0 <- RunOptions[[1]]$IndPeriod_Run[1] - 1
+  ts.start <- RunOptions[[1]]$IndPeriod_Run[1]
+  ts.end <- RunOptions[[1]]$IndPeriod_Run[length(RunOptions[[1]]$IndPeriod_Run)]
+  superTSstarts <- seq(ts.start, ts.end, x$.TimeStep)
+  lSuperTS <- lapply(
+    superTSstarts, function(x, TS, xMax) {
+      seq(x, min(x + TS - 1, xMax))
+    },
+    TS = x$.TimeStep,
+    xMax = ts.end
+  )
+
+  # Run runoff model for each sub-basin
+  x$OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
+    RunModel.GR(IM,
+                RunOptions = RunOptions[[IM$id]],
+                Param = Param[[IM$id]])
+    })
+  class(x$OutputsModel) <- append(class(x$OutputsModel), "GRiwrmOutputsModel")
+  # Copy simulated pure runoff flows (no SD nodes) to Qupstream in downstream SD nodes
+  for(id in getNoSD_Ids(x$InputsModel)) {
+    downId <- x$InputsModel[[id]]$down
+    x$InputsModel[[downId]]$Qupstream[RunOptions[[downId]]$IndPeriod_Run, id] <-
+      x$OutputsModel[[id]]$Qsim
+  }
+
+  # Save Qsim for step by step simulation
+  Qsim <- lapply(x$OutputsModel, function(OM) {
+    OM$Qsim
+  })
+
+  # Adapt RunOptions to step by step simulation
+  for(id in getSD_Ids(x$InputsModel)) {
+    RunOptions[[id]]$IndPeriod_WarmUp <- 0L
+    RunOptions[[id]]$Outputs_Sim <- "StateEnd"
+  }
+
+  # Loop over time steps with a step equal to the supervision time step
+  for(iTS in lSuperTS) {
+    # Run regulation on the whole basin for the current time step
+    x$ts.index <- iTS - x$ts.index0
+    x$ts.date <- x$InputsModel[[1]]$DatesR[iTS]
+    # Regulation occurs from second time step
+    if(iTS[1] > ts.start) {
+      doSupervision(x)
+    }
+
+    # Loop over sub-basin using SD model
+    for(id in getSD_Ids(x$InputsModel)) {
+      # Run the SD model for the sub-basin and one time step
+      RunOptions[[id]]$IndPeriod_Run <- iTS
+      RunOptions[[id]]$IniStates <- unlist(x$OutputsModel[[id]]$StateEnd)
+      x$OutputsModel[[id]] <- RunModel.SD(
+        x$InputsModel[[id]],
+        RunOptions = RunOptions[[id]],
+        Param = Param[[id]],
+        QsimDown = Qsim[[id]][x$ts.index]
+      )
+      # Storing Qsim in the data.frame Qsim
+      Qsim[[id]][x$ts.index] <- x$OutputsModel[[id]]$Qsim
+      # Routing Qsim to the downstream node
+      if(!is.na(x$InputsModel[[id]]$down)) {
+        x$InputsModel[[x$InputsModel[[id]]$down]]$Qupstream[iTS, id] <-
+          x$OutputsModel[[id]]$Qsim
+      }
+    }
+    x$ts.previous <- x$ts.index
+  }
+  for(id in getSD_Ids(x$InputsModel)) {
+    x$OutputsModel[[id]]$Qsim <- Qsim[[id]]
+  }
+  attr(x$OutputsModel, "Qm3s") <- OutputsModelQsim(x$InputsModel, x$OutputsModel, RunOptions[[1]]$IndPeriod_Run)
+  return(x$OutputsModel)
+}
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
new file mode 100644
index 0000000000000000000000000000000000000000..b246b1f4599f9b126fff7d82f1a808e2be0ce659
--- /dev/null
+++ b/R/RunModel_Lag.R
@@ -0,0 +1,111 @@
+RunModel_Lag <- function(InputsModel, RunOptions, Param) {
+  NParam <- 1
+
+  ##Arguments_check
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }
+  if (!inherits(InputsModel, "SD")) {
+    stop("'InputsModel' must be of class 'SD'")
+  }
+  if (!inherits(RunOptions, "RunOptions")) {
+    stop("'RunOptions' must be of class 'RunOptions'")
+  }
+  if (!is.vector(Param) | !is.numeric(Param)) {
+    stop("'Param' must be a numeric vector")
+  }
+  if (sum(!is.na(Param)) != NParam) {
+    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
+  }
+  if (is.null(InputsModel$OutputsModel)) {
+    stop(
+      "'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment"
+    )
+  }
+  if (is.null(InputsModel$OutputsModel$Qsim)) {
+    stop(
+      "'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment"
+    )
+  }
+  if (sum(!is.na(InputsModel$OutputsModel$Qsim)) != length(RunOptions$IndPeriod_Run)) {
+    stop(
+      "'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA"
+    )
+  }
+
+  OutputsModel <- InputsModel$OutputsModel
+  OutputsModel$QsimDown <- OutputsModel$Qsim
+
+  if (inherits(InputsModel, "hourly")) {
+    TimeStep <- 60 * 60
+  } else if (inherits(InputsModel, "daily")) {
+    TimeStep <- 60 * 60 * 24
+  } else {
+    stop("'InputsModel' should be of class \"daily\" or \"hourly\"")
+  }
+
+  # propagation time from upstream meshes to outlet
+  PT <- InputsModel$LengthHydro / Param[1L] / TimeStep
+  HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
+
+  NbUpBasins <- length(InputsModel$LengthHydro)
+  LengthTs <- length(OutputsModel$QsimDown)
+  OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+
+  IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
+  if (length(IniSD) > 0) {
+    if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
+      stop(
+        sprintf(
+          "SD initial states has a length of %i and a length of %i is required",
+          length(IniSD),
+          sum(floor(PT)) + NbUpBasins
+        )
+      )
+    }
+    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+      iStart <- 1
+      if (x > 1) {
+        iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
+      }
+      IniSD[iStart:(iStart + PT[x])]
+    })
+  } else {
+    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+      rep(0, floor(PT[x] + 1))
+    })
+  }
+  #message("Initstates: ",paste(IniStates, collapse = ", "))
+
+  for (upstream_basin in seq_len(NbUpBasins)) {
+    Qupstream <- c(IniStates[[upstream_basin]],
+                   InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin])
+    if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
+      # Upstream flow with area needs to be converted to m3 by time step
+      Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
+    }
+    #message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
+    OutputsModel$Qsim <-
+      OutputsModel$Qsim +
+      Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
+      Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
+  }
+  # Warning for negative flows
+  if (any(OutputsModel$Qsim < 0)) {
+    warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
+    OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
+  }
+  # Convert back Qsim to mm
+  OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+  #message("Qsim: ",paste(OutputsModel$Qsim, collapse = ", "))
+
+  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
+      lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
+      InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
+    })
+    #message("StateEnd: ",paste(OutputsModel$StateEnd$SD, collapse = ", "))
+  }
+
+  return(OutputsModel)
+}
diff --git a/R/createController.R b/R/createController.R
new file mode 100644
index 0000000000000000000000000000000000000000..bb58dd9699d6f9b40ff751bfa8991f8c9645109c
--- /dev/null
+++ b/R/createController.R
@@ -0,0 +1,71 @@
+#' Create and add a controller in a supervisor
+#'
+#' @details
+#' `ctrl.id` parameter is a unique id for finding the controller in the supervisor.
+#' If a controller with the same id already exists, it is overwritten by this new one.
+#'
+#' `FUN` parameter should be a function with one [numeric] parameter.
+#' This parameter will receive the measured values of at `Y` locations as input
+#' for the previous time step and returns calculated `U`. These `U` will then be applied
+#' at their location for the current time step of calculation of the model.
+#'
+#' @param supervisor `Supervisor` object, see [CreateSupervisor]
+#' @param ctrl.id [character] id of the controller (see Details)
+#' @param Y [character] location of the controlled and/or measured variables in the model. See [createControl]
+#' @param U [character] location of the command variables in the model. See [createControl]
+#' @param FUN [function] controller logic which calculates `U` from `Y` (see Details)
+#'
+#' @return `Controller`
+#' @export
+#'
+#' @examples
+#' # First create a Supervisor from a model
+#' example("CreateSupervisor")
+#' # A controller which usually releases 0.1 m3/s and provides
+#' # extra release if the downstream flow is below 0.5 m3/s
+#' logicDamRelease <- function(Y) max(0.5 - Y[1], 0.1)
+#' CreateController(sv, "DamRelease", Y = c("54001"), U = c("54095"), FUN = logicDamRelease)
+CreateController <- function(supervisor, ctrl.id, Y, U, FUN){
+
+  if(!is.character(ctrl.id)) stop("Parameter `ctrl.id` should be character")
+
+  FUN <- match.fun(FUN)
+
+  ctrlr <- list(
+    id = ctrl.id,
+    U = createControl(U),
+    Unames = U,
+    Y = createControl(Y),
+    Ynames = Y,
+    FUN = FUN
+  )
+  class(ctrlr) <- c("Controller", class(ctrlr))
+
+  # Function called from Supervisor environment
+  #environment(ctrlr$FUN) <- supervisor
+  if(!is.null(supervisor$controllers[[ctrl.id]])) {
+    warning("Controller '", ctrl.id, "' already exists in the supervisor: overwriting")
+  }
+  supervisor$controllers[[ctrl.id]] <- ctrlr
+  message("The controller has been added to the supervisor")
+  invisible(ctrlr)
+}
+
+#' Create a list of controls for command (U) and controlled variables (Y)
+#'
+#' @param locations vector of [character] containing the location of the variable in the model (see details)
+#'
+#' @return a [matrix] with each column is the location of the variables and the rows
+#' the values of the variable for the current time steps (empty by default)
+#' @export
+#'
+#' @examples
+#' # For pointing the discharge at the oulet of basins "54095" and "54002"
+#' createControl(c("54095", "54002"))
+createControl <- function(locations) {
+  if(!is.character(locations)) {
+    stop("Parameter `locations` should be character")
+  }
+  m <- matrix(NA, ncol = length(locations), nrow = 0)
+  return(m)
+}
diff --git a/R/plot.Qm3s.R b/R/plot.Qm3s.R
new file mode 100644
index 0000000000000000000000000000000000000000..80793de81ea3541f3ae7c2d69e581418d2dff098
--- /dev/null
+++ b/R/plot.Qm3s.R
@@ -0,0 +1,45 @@
+#' Plot a `Qm3s` object (time series of simulated flows)
+#'
+#' @param x a [data.frame] with a first column with [POSIXt] dates and followings columns with flows at each node of the network
+#' @param type 1-character string (See [plot.default]), default "l"
+#' @param xlab a label for the x axis, defaults to "Date"
+#' @param ylab a label for the y axis, defaults to "Flow (m3/s)"
+#' @param main a main title for the plot, defaults to "Simulated flows"
+#' @param col plotting color (See [par]), defaults to rainbow colors
+#' @param legend See parameter `legend` of [legend]. Set to [NULL] to not display the legend
+#' @param legend.cex `cex` parameter for the text of the legend (See [par])
+#' @param lty The line type (See [par])
+#' @param ... Further arguments to pass to the [matplot] functions
+#'
+#' @importFrom grDevices rainbow
+#' @importFrom graphics matplot
+#' @export
+#'
+plot.Qm3s <- function(x,
+                      type = 'l',
+                      xlab = "Date",
+                      ylab = "Flow (m3/s)",
+                      main = "Simulated flows",
+                      col = rainbow(ncol(x) - 1),
+                      legend = colnames(x)[-1],
+                      legend.cex = 0.7,
+                      lty = 1,
+                      ...) {
+  matplot(
+    x$DatesR,
+    x[, -1],
+    type = type,
+    lty = lty,
+    xlab = xlab,
+    ylab = ylab,
+    main = main,
+    col = col, ...
+  )
+  if(!is.null(legend)) {
+    legend('topright',
+           legend = legend,
+           cex = legend.cex,
+           lty = lty,
+           col = col)
+  }
+}
diff --git a/R/utils.R b/R/utils.R
new file mode 100644
index 0000000000000000000000000000000000000000..d98dd6a4fce8c2a83fbf81809ac94cf7ef5c90db
--- /dev/null
+++ b/R/utils.R
@@ -0,0 +1,155 @@
+#' Id of sub-basins using SD model
+#'
+#' @param InputsModel `GRiwrmInputsModel` object
+#'
+#' @return [character] IDs of the sub-basins using SD model
+#'
+getSD_Ids <- function(InputsModel) {
+  if (!inherits(InputsModel, "GRiwrmInputsModel")) {
+    stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
+  }
+  bSDs <- sapply(InputsModel, function (IM) {
+    inherits(IM, "SD")
+  })
+  names(InputsModel)[bSDs]
+}
+
+#' Id of sub-basins not using SD model
+#'
+#' @param InputsModel `GRiwrmInputsModel` object
+#'
+#' @return [character] IDs of the sub-basins not using SD model
+#'
+getNoSD_Ids <- function(InputsModel) {
+  if (!inherits(InputsModel, "GRiwrmInputsModel")) {
+    stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
+  }
+  bSDs <- sapply(InputsModel, function (IM) {
+    !inherits(IM, "SD")
+  })
+  names(InputsModel)[bSDs]
+}
+
+
+#' Retrieve data in the model for the current time steps
+#'
+#' This function should be call inside a Supervisor
+#'
+#' @param loc location of the data
+#' @param sv a `Supervisor` (See [CreateSupervisor])
+#'
+#' @return [numeric] retrieved data at the location
+getDataFromLocation <- function(loc, sv) {
+  if (length(grep("\\[[0-9]+\\]$", loc)) > 0) {
+    stop("Reaching output of other controller is not implemented yet")
+  } else {
+    node <- sv$griwrm$down[sv$griwrm$id == loc]
+    sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.previous, loc]
+  }
+}
+
+
+#' Write data to model input for the current time step
+#'
+#' @param ctrlr a `Controller` object (See [CreateController])
+#' @param sv `Supervisor` (See [CreateSupervisor])
+#'
+#' @return [NULL]
+setDataToLocation <- function(ctrlr, sv) {
+  l <- lapply(seq(length(ctrlr$Unames)), function(i) {
+    node <- sv$griwrm$down[sv$griwrm$id == ctrlr$Unames[i]]
+    # limit U size to the number of simulation time steps of the current supervision time step
+    U <- ctrlr$U[seq.int(length(sv$ts.index)),i]
+    # ! Qupstream contains warm up period and run period => the index is shifted
+    sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index, ctrlr$Unames[i]] <- U
+  })
+}
+
+
+#' Do the supervision for the current time step
+#'
+#' @param supervisor `Supervisor` (See [CreateSupervisor])
+#'
+doSupervision <- function(supervisor) {
+  for (id in names(supervisor$controllers)) {
+    supervisor$controller.id <- id
+    # Read Y from locations in the model
+    supervisor$controllers[[id]]$Y <- do.call(
+      cbind,
+      lapply(supervisor$controllers[[id]]$Ynames, getDataFromLocation, sv = supervisor)
+    )
+    # Run logic
+    supervisor$controllers[[id]]$U <-
+      supervisor$controllers[[id]]$FUN(supervisor$controllers[[id]]$Y)
+    # Write U to locations in the model
+    setDataToLocation(supervisor$controllers[[id]], sv = supervisor)
+  }
+}
+
+
+#' Check the parameters of RunModel methods
+#'
+#' Stop the execution if an error is detected.
+#'
+#' @param InputsModel a `GRiwrmInputsModel` object (See [CreateInputsModel.GRiwrm])
+#' @param RunOptions a `GRiwrmRunOptions` object (See [CreateRunOptions.GRiwrmInputsModel])
+#' @param Param a [list] of [numeric] containing model parameters of each node of the network
+#'
+checkRunModelParameters <- function(InputsModel, RunOptions, Param) {
+  if(!inherits(InputsModel, "GRiwrmInputsModel")) stop("`InputsModel` parameter must of class 'GRiwrmRunoptions' (See ?CreateRunOptions.GRiwrmInputsModel)")
+  if(!inherits(RunOptions, "GRiwrmRunOptions")) stop("Argument `RunOptions` parameter must of class 'GRiwrmRunOptions' (See ?CreateRunOptions.GRiwrmInputsModel)")
+  if(!is.list(Param) || !all(names(InputsModel) %in% names(Param))) stop("Argument `Param` must be a list with names equal to nodes IDs")
+}
+
+
+#' Create a data.frame with simulated flows at each nodes of the [GRiwrm] object
+#'
+#' @details
+#' This function can only be called inside [RunModel.GRiwrmInputsModel] or [RunModel.Supervisor]
+#' because it needs a `GRiwrmInputsModel` object internally modified by these functions
+#' (`Qupstream` updated with simulated flows).
+#'
+#' @param InputsModel a `GRiwrmInputsModel` object created by [CreateInputsModel.GRiwrm]
+#' @param OutputsModel a `GRiwrmOutputsModel` object created by [RunModel.GRiwrmInputsModel] or [RunModel.Supervisor]
+#' @param IndPeriod_Run an [integer] vector (See [airGR::CreateRunOptions])
+#'
+#' @return a [data.frame] containing the simulated flows (in m3/time step) structured with the following columns:
+#' - 'DatesR' containing the timestamps of the time series
+#' - one column by node with the simulated flows
+#'
+OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
+  griwrm <- attr(InputsModel, "GRiwrm")
+  # Get simulated flow for each node
+  # Flow for each node is available in InputsModel$Qupstream except for the downstream node
+  upperNodes <- griwrm$id[!is.na(griwrm$down)]
+  lQsim <- lapply(
+    upperNodes,
+    function(x, griwrm, IndPeriod_Run) {
+      node <- griwrm$down[griwrm$id == x]
+      InputsModel[[node]]$Qupstream[IndPeriod_Run, x]
+    },
+    griwrm = griwrm, IndPeriod_Run = IndPeriod_Run
+  )
+  names(lQsim) <- upperNodes
+  # Flow of the downstream node is only available in OutputsModel[[node]]$Qsim
+  downNode <- names(InputsModel)[length(InputsModel)]
+  lQsim[[downNode]] <- OutputsModel[[downNode]]$Qsim
+
+  # Conversion to m3/s
+  lQsim <- lapply(
+    names(lQsim),
+    function(x) {
+      i <- which(griwrm$id == x)
+      if(is.na(griwrm$area[i])) { # m3/time step => m3/s
+        return(lQsim[[x]] / attr(InputsModel, "TimeStep"))
+      } else { # mm/time step => m3/s
+        return(lQsim[[x]] * griwrm$area[i] * 1E3 / attr(InputsModel, "TimeStep"))
+      }
+    }
+  )
+  names(lQsim) <- c(upperNodes, downNode)
+  dfQsim <- cbind(data.frame(DatesR = as.POSIXct(InputsModel[[1]]$DatesR[IndPeriod_Run])),
+                  do.call(cbind,lQsim))
+  class(dfQsim) <- c("Qm3s", class(dfQsim)) # For S3 methods
+  return(dfQsim)
+}
diff --git a/tests/testthat/test-RunModel.R b/tests/testthat/test-RunModel.R
new file mode 100644
index 0000000000000000000000000000000000000000..a5f5a720038957e3d5e86324a8df33e720be4cf6
--- /dev/null
+++ b/tests/testthat/test-RunModel.R
@@ -0,0 +1,99 @@
+context("RunModel.Supervisor")
+
+# Load data
+data(Severn)
+
+# Network configuration
+nodes <- Severn$BasinsInfo[c(1,2,5), c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes$distance_downstream <- nodes$distance_downstream * 1000 # Conversion km -> m
+nodes$model <- NA
+nodes$model[1] <- "RunModel_GR4J"
+griwrm <- GRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
+
+# InputsModel
+DatesR <- Severn$BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$peti}))
+Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
+InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+
+# RunOptions
+nTS <- 365
+IndPeriod_Run <- seq(
+  length(InputsModel[[1]]$DatesR) - nTS + 1,
+  length(InputsModel[[1]]$DatesR)
+)
+IndPeriod_WarmUp = seq(IndPeriod_Run[1]-366,IndPeriod_Run[1]-1)
+RunOptions <- CreateRunOptions(
+  InputsModel = InputsModel,
+  IndPeriod_WarmUp = IndPeriod_WarmUp,
+  IndPeriod_Run = IndPeriod_Run
+)
+
+# RunModel.GRiwrmInputsModel
+Param <- list("54057" = c(0.727,  175.493,   -0.082,    0.029,    4.654))
+OM_GriwrmInputs <- RunModel(
+  InputsModel,
+  RunOptions = RunOptions,
+  Param = Param
+)
+
+test_that("RunModelSupervisor with no regulation should returns same results as RunModel.GRiwrmInputsModel", {
+  sv <- CreateSupervisor(InputsModel)
+  OM_Supervisor <- RunModel(
+    sv,
+    RunOptions = RunOptions,
+    Param = Param
+  )
+  expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})
+
+# Add 2 nodes to the network
+griwrm2 <- rbind(griwrm,
+                 data.frame(
+                   id = c("R1", "R2"),
+                   down = "54057",
+                   length = 100000,
+                   area = NA,
+                   model = NA
+                 ))
+# Add Qobs for the 2 new nodes and create InputsModel
+Qobs2 <- cbind(Qobs, matrix(data = rep(0, 2*nrow(Qobs)), ncol = 2))
+colnames(Qobs2) <- c(colnames(Qobs2)[1:6], "R1", "R2")
+InputsModel <- CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qobs2)
+
+test_that("RunModelSupervisor with two regulations that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
+  # Create Supervisor
+  sv <- CreateSupervisor(InputsModel)
+  # Function to withdraw half of the measured flow
+  fWithdrawal <- function(y) { -y/2 }
+  # Function to release half of the the measured flow
+  fRelease <- function(y) { y/2 }
+  # Controller that withdraw half of the flow measured at node "54002" at location "R1"
+  CreateController(sv, "Withdrawal", Y = c("54002"), U = c("R1"), FUN = fWithdrawal)
+  # Controller that release half of the flow measured at node "54002" at location "R2"
+  CreateController(sv, "Release", Y = c("54002"), U = c("R2"), FUN = fRelease)
+
+  OM_Supervisor <- RunModel(
+    sv,
+    RunOptions = RunOptions,
+    Param = Param
+  )
+  expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})
+
+test_that("RunModelSupervisor with multi time steps controller, two regulations in 1 centralised controller that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
+  sv <- CreateSupervisor(InputsModel, TimeStep = 10L)
+  fEverything <- function(y) {
+    matrix(c(y[,1]/2, -y[,1]/2), ncol = 2)
+  }
+  CreateController(sv, "Everything", Y = c("54002", "54032"), U = c("R1", "R2"), FUN = fEverything)
+  OM_Supervisor <- RunModel(
+    sv,
+    RunOptions = RunOptions,
+    Param = Param
+  )
+  expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})
diff --git a/vignettes/V01_Structure_SD_model.Rmd b/vignettes/V01_Structure_SD_model.Rmd
index 4761eecd89c78e9ec7a97016766635a178266dc5..4b733900b3e431cd4ef7f3254f35565c9f714a9b 100644
--- a/vignettes/V01_Structure_SD_model.Rmd
+++ b/vignettes/V01_Structure_SD_model.Rmd
@@ -16,7 +16,7 @@ knitr::opts_chunk$set(echo = TRUE)
 library(airGRiwrm)
 ```
 
-The package **airGRiwrm** is a modeling tool for integrated water resource management based on the package **airGR** package [See @coron_airgr_2020].
+The package **airGRiwrm** is a modeling tool for integrated water resource management based on the package **airGR** package [See @coronSuiteLumpedGR2017].
 
 In a semi-distributive model, the catchment is divided into several sub-catchments. Each sub-catchment is an hydrological entity where a runfall-runoff model produces a flow time series at the outlet of the sub-catchment. Then a hydraulic link is set between sub-catchment outlets to model the flow at the outlet of the whole catchment. The aim of **airGRiwrm** is to organise the structure and schedule the execution of the hydrological and hydraulic sub-models contained in the semi-distributive model.
 
@@ -24,7 +24,7 @@ In this vignette, we show how to prepare observation data for the model.
 
 ## Description of the example used in this tutorial
 
-The example of this tutorial takes place on the Severn River in United Kingdom. The data set comes from the CAMEL GB database [See @coxon_catchment_2020].
+The example of this tutorial takes place on the Severn River in United Kingdom. The data set comes from the CAMEL GB database [See @coxonCatchmentAttributesHydrometeorological2020].
 
 ```{r}
 data(Severn)
diff --git a/vignettes/V02_Calibration_SD_model.Rmd b/vignettes/V02_Calibration_SD_model.Rmd
index 319e3dc32282aae5281c0eb41f8a01bbfa177654..d15f9846e15c34b6dcfcfa7a5a075e35a03a32d6 100644
--- a/vignettes/V02_Calibration_SD_model.Rmd
+++ b/vignettes/V02_Calibration_SD_model.Rmd
@@ -128,7 +128,7 @@ ParamMichel <- sapply(griwrm$id, function(x) {OutputsCalib[[x]]$Param})
 
 ```{r RunModel}
 OutputsModels <- RunModel(
-  InputsModel = InputsModel,
+  InputsModel,
   RunOptions = RunOptions,
   Param = ParamMichel
 )
@@ -145,8 +145,9 @@ htmltools::tagList(lapply(
 ))
 ```
 
+The resulting flow of each node in m3/s is directly available and can be plotted with these commands:
 
-
-
-
-
+```{r}
+Qm3s <- attr(OutputsModels, "Qm3s")
+plot(Qm3s[1:100,])
+```
diff --git a/vignettes/V03_Open-loop_influenced_flow.Rmd b/vignettes/V03_Open-loop_influenced_flow.Rmd
index d98c37bbd68f21933ba35c3eb3f80b8b2fd0a663..463bbe62858d06425e8be802d1a245ef785f63a7 100644
--- a/vignettes/V03_Open-loop_influenced_flow.Rmd
+++ b/vignettes/V03_Open-loop_influenced_flow.Rmd
@@ -23,15 +23,15 @@ library(airGRiwrm)
 
 
 
-We can see in vignette 'V02_Calibration_SD_model' that calibration result for the discharge simulated on the Avon at Evesham (Gauging station '54002') and on the Severn at Buildwas (Gauging station '54095') is not satisfactory. These upper basins are heavily influenced by impoundments and inter-basin transfers [@higgs_hydrological_1988].
+We can see in vignette 'V02_Calibration_SD_model' that calibration result for the flow simulated on the Avon at Evesham (Gauging station '54002') and on the Severn at Buildwas (Gauging station '54095') is not satisfactory. These upper basins are heavily influenced by impoundments and inter-basin transfers [@higgsHydrologicalChangesRiver1988].
 
-To cope with this influenced flow, we won't try to simulate it with an hydrological model and we will directly inject this discharge measurement in the model.
+To cope with this influenced flow, we won't try to simulate it with a hydrological model and we will directly inject in the model the observed flow at these nodes.
 
 ## Convert a gauging station into a release spot
 
 ### Modifying the GRiwrm object
 
-The production method of the `GRiwrm` object is detailed in the vignette "01_First_network" of the package.
+The production method of the `GRiwrm` object is detailed in the vignette "01_First_network" of the package (and also in the documentation, type `?GRiwrm` to access to it).
 
 Run `vignette("01_Structure_SD_model", package = "airGRiwrm")` before this one in order to create the Rdata file containing the `GRiwrmInputsModel` object:
 
@@ -42,10 +42,10 @@ load("_cache/V01.RData")
 To notify the SD model that the provided flow on a node should be directly used instead of an hydrological model, you only need to declare its model as `NA`:
 
 ```{r}
-griwrm_OL <- griwrm
-griwrm_OL$model[griwrm$id == "54002"] <- NA
-griwrm_OL$model[griwrm$id == "54095"] <- NA
-griwrm_OL
+griwrmV03 <- griwrm
+griwrmV03$model[griwrm$id == "54002"] <- NA
+griwrmV03$model[griwrm$id == "54095"] <- NA
+griwrmV03
 ```
 
 Here, we keep the area of this basin which means that the discharge will be provided in mm per time step. If the discharge is provided in m<sup>3</sup>/s, then the area should be set to `NA` and downstream basin areas should be modified subsequently. 
@@ -57,7 +57,7 @@ The diagram of the network structure is represented below with:
   * in red the node with direct flow injection (no hydrological model)
 
 ```{r diagram}
-DiagramGRiwrm(griwrm_OL)
+DiagramGRiwrm(griwrmV03)
 ```
 
 ### Generate the GRiwrmInputsModel object
@@ -66,7 +66,7 @@ Since the network description has changed, the `GRiwrmInputsModel` should be gen
 
 ```{r}
 load("_cache/V01b.RData")
-IM_OL <- CreateInputsModel(griwrm_OL, DatesR, Precip, PotEvap, Qobs)
+IM_OL <- CreateInputsModel(griwrmV03, DatesR, Precip, PotEvap, Qobs)
 ```
 
 ## Calibration of the new model
@@ -82,16 +82,21 @@ The **airGR** calibration process is applied on each hydrological node of the `G
 ```{r Calibration}
 OC_OL <- suppressWarnings(
   Calibration(IM_OL, RunOptions, InputsCrit, CalibOptions))
-Param_OL <- sapply(griwrm$id, function(x) {OC_OL[[x]]$Param})
+ParamV03 <- sapply(griwrm$id, function(x) {OC_OL[[x]]$Param})
 ```
 
+```{r}
+save(griwrmV03, ParamV03, file = "_cache/V03.RData")
+```
+
+
 ## Run model with this newly calibrated parameters
 
 ```{r RunModel}
 OM_OL <- RunModel(
-  InputsModel = IM_OL,
+  IM_OL,
   RunOptions = RunOptions,
-  Param = Param_OL
+  Param = ParamV03
 )
 ```
 
@@ -108,4 +113,12 @@ htmltools::tagList(lapply(
 ))
 ```
 
+The resulting flow of each node in m3/s is directly available and can be plotted with these commands:
+
+```{r}
+Qm3s <- attr(OM_OL, "Qm3s")
+plot(Qm3s[1:100,])
+```
+
+
 # References
diff --git a/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd b/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..862652b18aa56e3b7ff811bd5c8ec0a7e60217ae
--- /dev/null
+++ b/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd
@@ -0,0 +1,336 @@
+---
+title: "Example of a regulated withdrawal (closed-loop control)"
+output: 
+  rmarkdown::html_vignette:
+    toc: true
+    toc_depth: 2
+vignette: >
+  %\VignetteIndexEntry{Example of a regulated withdrawal (closed-loop control)}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+bibliography: airGRiwrm.bib
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+library(airGRiwrm)
+```
+
+# Presentation of the case study
+
+Starting from the network and the calibration set in the vignette "V03_Open-loop_influenced_flow", we add 2 intake points for irrigation.
+
+Vignettes "V01_Structure_SD_model.Rmd" and  "V03_Open-loop_influenced_flow" should be run in order to load input data and parameters of the model:
+
+```{r load_cache}
+load("_cache/V01b.RData")
+load("_cache/V02.RData")
+load("_cache/V03.RData")
+```
+
+## Network configuration
+
+The intake points are located: 
+
+- on the Severn at 35 km upstream Bewdley (Gauging station '54001');
+- on the Severn at 10 km upstream Saxons Lode (Gauging station '54032').
+
+We have to add this 2 nodes in the `GRiwrm` object that describes the network:
+
+```{r griwrm}
+griwrmV04 <- rbind(
+  griwrmV03,
+  data.frame(
+    id = c("Irrigation1", "Irrigation2"),
+    down = c("54001", "54032"),
+    length = c(35000, 10000),
+    model = NA,
+    area = NA
+  )
+)
+DiagramGRiwrm(griwrmV04)
+```
+Nodes in red color are direct injection points (positive or negative flow) in the model, green nodes simulate either sub-basin rainfall-runoff and routing of upstream node flows and grey nodes simulate sub-basin rainfall-runoff flows only.
+
+It's important to notice that even if the points "Irrigation1" and "Irrigation2" are located on Severn such as gauging stations "54095", "54001", "54032", there respective nodes are not represented on the same path in the model. Consequently, with this network configuration, it is not possible to know the value of the flow in the Severn river at "Irrigation1" or "Irrigation2" nodes. These values are only available in nodes "54095", "54001", "54032".
+
+## Irrigation objectives and flow demand at intakes
+
+Irrigation perimeters #1 cover an area of 15 km² and the perimeter #2 an area of 30 km².
+
+The objective of these irrigation system is to cover the deficit in rain [@burtIrrigationPerformanceMeasures1997] with 80% of success. Here's below the calculation of the 8th decile of monthly water need given meteorological data of catchments "54001" and "54032" (unit mm/day) :
+
+```{r monthly_water_need}
+dailyWaterNeed <- PotEvap - Precip
+dailyWaterNeed <- cbind(as.data.frame(DatesR), dailyWaterNeed[,c("54001", "54032")])
+monthlyWaterNeed <- airGR::SeriesAggreg(dailyWaterNeed, "%m", rep("q80",2))
+monthlyWaterNeed$DatesR <- as.numeric(format(monthlyWaterNeed$DatesR,"%m"))
+names(monthlyWaterNeed)[1] <- "month"
+monthlyWaterNeed
+```
+
+We bound the irrigation season between April and August, then the need for the crops can the be expressed in m3/s as follow:
+
+```{r}
+irrigationObjective <- monthlyWaterNeed
+# Conversion in m3/day
+irrigationObjective$"54001" <- monthlyWaterNeed$"54001" * 15 * 1E3
+irrigationObjective$"54032" <- monthlyWaterNeed$"54032" * 30 * 1E3
+# Irrigation period between April and August
+irrigationObjective[-seq(4,8),-1] <- 0 
+# Conversion in m3/s
+irrigationObjective[,c(2,3)] <- round(irrigationObjective[,c(2,3)] / 86400, 1)
+irrigationObjective$total <- rowSums(irrigationObjective[,c(2,3)])
+irrigationObjective
+```
+
+With an hypothesis of an efficiency of the irrigation systems of 50% as proposed by @secklerConceptEfficiencyWaterresources2003, the flow demand at intake for each irrigation perimeter is as follow (unit: m3/s):
+
+```{r}
+# Apply 50% of efficiency on the water demand
+irrigationObjective[,seq(2,4)] <- irrigationObjective[,seq(2,4)] / 0.5
+# Display result in m3/s
+irrigationObjective
+```
+
+
+
+## Restriction of irrigation in case of water scarcity
+
+### Minimal environmental flow at the intakes
+
+In UK, abstraction restrictions are driven by Environmental Flow Indicator (EFI) supporting Good Ecological Status (GES) [@klaarDevelopingHydroecologicalModels2014]. Abstraction restriction consists in limiting the proportion of available flow for abstraction in function of the current flow regime (Reference taken for a river highly sensitive to abstraction classified "ASB3").
+
+```{r abstraction restriction rule}
+restriction_rule <- data.frame(quantile_natural_flow = c(.05, .3, 0.5, 0.7),
+           abstraction_rate = c(0.1, 0.15, 0.20, 0.24))
+```
+
+The control of the abstraction will be done at the gauging station downstream all the abstraction locations (node "54032"). So we need the flow corresponding to the quantiles of natural flow and flow available for abstraction in each case.
+
+```{r}
+quant_m3s32 <- quantile(
+  Qobs[,"54032"] * griwrmV04[griwrmV04$id == "54032", "area"] / 86.4,
+  restriction_rule$quantile_natural_flow,
+  na.rm = TRUE
+)
+restriction_rule_m3s <- data.frame(
+  threshold_natural_flow = quant_m3s32,
+  abstraction_rate = restriction_rule$abstraction_rate
+)
+
+matplot(restriction_rule$quantile_natural_flow, 
+        cbind(restriction_rule_m3s$threshold_natural_flow, 
+              restriction_rule$abstraction_rate * restriction_rule_m3s$threshold_natural_flow,
+              max(irrigationObjective$total)), 
+        log = "x", type = "l", 
+        main = "Quantiles of flow on the Severn at Saxons Lode (54032)",
+        xlab = "quantiles", ylab = "Flow (m3/s)",
+        lty = 1, col = rainbow(3, rev = TRUE)
+        )
+legend("topleft", legend = c("Natural flow", "Abstraction limit", "Irrigation max. objective"), 
+       col = rainbow(3, rev = TRUE), lty = 1)
+```
+The water availability or abstraction restriction depending on the natural flow is calculated with the function below:
+
+```{r}
+# A function to enclose the parameters in the function (See: http://adv-r.had.co.nz/Functional-programming.html#closures)
+getAvailableAbstractionEnclosed <- function(restriction_rule_m3s) {
+  function(Qnat) approx(restriction_rule_m3s$threshold_natural_flow,
+                        restriction_rule_m3s$abstraction_rate,
+                        Qnat, 
+                        rule = 2)
+}
+# The function with the parameters inside it :)
+getAvailableAbstraction <- getAvailableAbstractionEnclosed(restriction_rule_m3s)
+# You can check the storage of the parameters in the function with 
+as.list(environment(getAvailableAbstraction))
+```
+
+
+### Restriction rules
+
+The figure below shows that restrictions will be imposed to the irrigation perimeter if the natural flow at Saxons Lode is under around 20 m3/s.
+
+Applying restriction on the intake on a real field is always challenging since it's difficult to regulate day by day the flow at the intake. Policy makers often decide to close the irrigation abstraction points in turn several days a week based on the mean flow of the previous week.
+
+The number of authorized days per week for irrigation can be calculated as follow. All the calculation are based on the mean flow measured the week previous the current time step. First, the naturalised flow $N$ is equal to 
+
+$$ N = M + I_l $$
+with: 
+
+- $M$, the measured flow at the downstream gauging station
+- $I_l$, the total abstracted flow for irrigation for the last week
+
+Available flow for abstraction $A$ is: 
+
+$$A = f_{a}(N)$$
+
+with $f_a$ the availability function calculated from quantiles of natural flow and related restriction rates.
+
+The flow planned for irrigation $Ip$ is then:
+
+$$ I_p = \min (O, A)$$
+
+with $O$ the irrigation objective flow.
+
+The number of day for irrigation $n$ per week is then equal to:
+
+$$ n = \lfloor \frac{I_p}{O} \times 7 \rfloor$$
+with $\lfloor x \rfloor$ the function that returns the largest integers not greater than $x$
+
+The rotation of restriction days between the 2 irrigation perimeters is operated as follow:
+
+```{r}
+restriction_rotation <- matrix(c(5,7,6,4,2,1,3,3,1,2,4,6,7,5), ncol = 2)
+m <- do.call(
+  rbind, 
+  lapply(seq(0,7), function(x) {
+    b <- restriction_rotation <= x
+    rowSums(b)
+  })
+)
+# Display the planning of restriction
+image(1:ncol(m), 1:nrow(m), t(m), col = heat.colors(3, rev = TRUE), 
+      axes = FALSE, xlab = "week day", ylab = "number of restriction days", 
+      main = "Number of closed irrigation perimeters")
+axis(1, 1:ncol(m), unlist(strsplit("SMTWTFS", "")))
+axis(2, 1:nrow(m), seq(0,7))
+for (x in 1:ncol(m))
+  for (y in 1:nrow(m))
+    text(x, y, m[y,x])
+```
+
+# Implementation of the model
+
+As for the previous model, we need to set up an `GRiwrmInputsModel` object containing all the model inputs:
+
+```{r}
+# a flow is needed for all the nodes in the network
+# even if it is overwritten after by a controller
+QobsIrrig <- cbind(Qobs, Irrigation1 = 0, Irrigation2 = 0)
+IM_Irrig <- CreateInputsModel(griwrmV04, DatesR, Precip, PotEvap, QobsIrrig)
+```
+
+# Implementation of the regulation controller
+
+## The supervisor
+
+The simulation is piloted through a `Supervisor` which can contain one or more `Controller`. This supervisor will work with a cycle of 7 days: the measurement are taken on the last 7 days and decisions are taken for each time step for the next seven days.
+
+```{r}
+sv <- CreateSupervisor(IM_Irrig, TimeStep = 7L)
+```
+
+## The control logic function
+
+We need a controller that measures the flow at Saxon Lodes ("54032") and adapts weekly the abstracted flow at the two irrigation points. The supervisor will stop the simulation every 7 days and provides to the controller the last 7 simulated flow at Saxon Lodes ("54032") (measured variables) and controller should provide "command variables" for the next 7 days for the 2 irrigation points.
+
+A control logic function should be provide to the controller. This control logic function processes the logic of the regulation taking measured flows as input and returning the "command variables". Both measured variables and command variables are of type `matrix` with the variables in columns and the time steps in rows.
+
+In this example, the logic function must do the following tasks:
+
+1. Calculate the objective of irrigation according to the month of the current days of simulation
+2. calculate the naturalised flow from the measured flow and the abstracted flow of the previous week
+3. calculate the number of days of restriction for each irrigation point
+4. Return the abstracted flow for the next week taking into account restriction days
+
+```{r}
+fIrrigationFactory <- function(supervisor, 
+                               irrigationObjective, 
+                               restriction_rule_m3s, 
+                               restriction_rotation) {
+  function(Y) {
+    # Y is in mm/day and the basin's area is 6864.88 km2
+    # Calculate the objective of irrigation according to the month of the current days of simulation
+    month <- as.numeric(format(supervisor$ts.date, "%m"))
+    U <- irrigationObjective[month, c(2,3)] # m3/s
+    meanU <- mean(rowSums(U))
+    if(meanU > 0) {
+      # calculate the naturalised flow from the measured flow and the abstracted flow of the previous week
+      lastU <- supervisor$controllers[[supervisor$controller.id]]$U # m3/day
+      Qnat <- (Y * 6864.88 * 1E3 - rowSums(lastU)) / 86400 # m3/s
+      # Maximum abstracted flow available
+      Qrestricted <- mean(
+        approx(restriction_rule_m3s$threshold_natural_flow,
+               restriction_rule_m3s$abstraction_rate,
+               Qnat, 
+               rule = 2)$y * Qnat
+      )
+      # Total for irrigation
+      QIrrig <- min(meanU, Qrestricted)
+      # Number of days of irrigation
+      n <- floor(7 * (1 - QIrrig / meanU))
+      # Apply days off
+      U[restriction_rotation[seq(nrow(U)),] <= n] <- 0
+    }
+    return(-U * 86400) # withdrawal is a negative flow on an upstream node
+  }
+}
+```
+
+You can notice that the data required for processing the control logic are enclosed in the function `fIrrigationFactory` which takes the required data as arguments and return the control logic function.
+
+Creating `fIrrigation` by calling `fIrrigationFactory` with the arguments currently in memory save these variables in the environment of the function:
+
+```{r}
+fIrrigation <- fIrrigationFactory(supervisor = sv, 
+                                  irrigationObjective = irrigationObjective, 
+                                  restriction_rule_m3s = restriction_rule_m3s, 
+                                  restriction_rotation = restriction_rotation)
+```
+You can see what data is available in the environment of the function with:
+
+```{r}
+str(as.list(environment(fIrrigation)))
+```
+
+The `supervisor` variable is itself an environment which means that the variables contained inside it will be updated during the simulation. Some of them are useful for computing the control logic such as:
+
+- `supervisor$ts.index`: indexes of the current time steps of simulation (In `IndPeriod_Run`)
+- `supervisor$ts.date`: date/time of the current time steps of simulation
+- `supervisor$controller.id`: identifier of the current controller
+- `supervisor$controllers`: the `list` of `Controller`
+
+## The controller
+
+The controller contains:
+
+- the location of the measured flows
+- the location of the control commands
+- the logic control function
+
+```{r}
+CreateController(sv, 
+                 ctrl.id = "Irrigation", 
+                 Y = "54032", 
+                 U = c("Irrigation1", "Irrigation2"), 
+                 FUN = fIrrigation)
+```
+
+# Running the simulation
+
+For running a model with a supervision, you only need to substitute `InputsModel` by a `Supervisor` in the `RunModel` function call.
+
+```{r}
+OM_Irrig <- RunModel(sv, RunOptions = RunOptions, Param = ParamV03)
+```
+
+Simulated flows can be extracted and plot as follows:
+
+```{r}
+Qm3s <- attr(OM_Irrig, "Qm3s")
+Qm3s <- Qm3s[Qm3s$DatesR > "2003-03-25" & Qm3s$DatesR < "2003-09-05",]
+par(mfrow=c(2,1), mar = c(2.5,4,1,1))
+plot(Qm3s[, c("DatesR", "54095", "54001", "54032")], main = "", xlab = "", ylim = c(0,100))
+plot(Qm3s[, c("DatesR", "Irrigation1", "Irrigation2")], main = "", xlab = "")
+```
+We can observe  that the irrigations points are alternatively closed some days a week when the flow at node "54032" becomes low.
+
+# References
diff --git a/vignettes/airGRiwrm.bib b/vignettes/airGRiwrm.bib
index 439f05384cc1aee453d1b6a638695fd82e1759dc..8cdd6a2e17b59aa1979526af9d6d83ee3e44a803 100644
--- a/vignettes/airGRiwrm.bib
+++ b/vignettes/airGRiwrm.bib
@@ -1,105 +1,164 @@
 
-@misc{coxon_catchment_2020,
-	title = {Catchment attributes and hydro-meteorological timeseries for 671 catchments across Great Britain ({CAMELS}-{GB})},
-	rights = {This resource is available under the terms of the Open Government Licence},
-	url = {https://catalogue.ceh.ac.uk/id/8344e4f3-d2ea-44f5-8afa-86d2987543a9},
-	abstract = {This dataset provides hydro-meteorological timeseries and landscape attributes for 671 catchments across Great Britain. It collates river flows, catchment attributes and catchment boundaries from the {UK} National River Flow Archive together with a suite of new meteorological timeseries and catchment attributes. Daily timeseries for the time period 1st October 1970 to the 30th September 2015 are provided for a range of hydro-meteorological data (including rainfall, potential evapotranspiration, temperature, radiation, humidity and flow). A comprehensive set of catchment attributes are quantified describing a range of catchment characteristics including topography, climate, hydrology, land cover, soils, hydrogeology, human influences and discharge uncertainty. This dataset is intended for the community as a freely available, easily accessible dataset to use in a wide range of environmental data and modelling analyses. A research paper (Coxon et al, {CAMELS}-{GB}: Hydrometeorological time series and landscape attributes for 671 catchments in Great Britain) describing the dataset in detail will be made available in Earth System Science Data (https://www.earth-system-science-data.net/).},
-	publisher = {{NERC} Environmental Information Data Centre},
-	author = {Coxon, G. and Addor, N. and Bloomfield, J.P. and Freer, J. and Fry, M. and Hannaford, J. and Howden, N.J.K. and Lane, R. and Lewis, M. and Robinson, E.L. and Wagener, T. and Woods, R.},
-	editora = {Coxon, Gemma and Environmental Information Data Centre},
-	editoratype = {collaborator},
-	urldate = {2020-12-18},
-	date = {2020},
-	langid = {english},
-	doi = {10.5285/8344E4F3-D2EA-44F5-8AFA-86D2987543A9},
-	note = {Medium: text/csv Comma-separated values ({CSV}),Shapefile
-type: dataset},
-	keywords = {Hydrology}
+@article{burtIrrigationPerformanceMeasures1997,
+  title = {Irrigation {{Performance Measures}}: {{Efficiency}} and {{Uniformity}}},
+  shorttitle = {Irrigation {{Performance Measures}}},
+  author = {Burt, C. M. and Clemmens, A. J. and Strelkoff, T. S. and Solomon, K. H. and Bliesner, R. D. and Hardy, L. A. and Howell, T. A. and Eisenhauer, D. E.},
+  year = {1997},
+  month = nov,
+  volume = {123},
+  pages = {423--442},
+  publisher = {{American Society of Civil Engineers}},
+  issn = {0733-9437},
+  doi = {10.1061/(ASCE)0733-9437(1997)123:6(423)},
+  abstract = {It is essential to standardize the definitions and approaches to quantifying various irrigation performance measures. The ASCE Task Committee on Defining Irrigation Efficiency and Uniformity provides a comprehensive examination of various performance indices such as irrigation efficiency, application efficiency, irrigation sagacity, distribution uniformity, and others. Consistency is provided among different irrigation methods and different scales. Clarification of common points of confusion is provided, and methods are proposed whereby the accuracy of numerical values of the performance indicators can be assessed. This issue has two companion papers that provide more detailed information on statistical distribution uniformity and the accuracy of irrigation efficiency estimates.},
+  copyright = {Copyright \textcopyright{} 1997 American Society of Civil Engineers},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\D9Q9VV6J\\Burt et al. - 1997 - Irrigation Performance Measures Efficiency and Un.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\5NLYUVSN\\(ASCE)0733-9437(1997)1236(423).html},
+  journal = {Journal of Irrigation and Drainage Engineering},
+  language = {EN},
+  number = {6}
 }
 
-@article{coron_suite_2017,
-	title = {The suite of lumped {GR} hydrological models in an R package},
-	volume = {94},
-	issn = {13648152},
-	url = {https://linkinghub.elsevier.com/retrieve/pii/S1364815217300208},
-	doi = {10.1016/j.envsoft.2017.05.002},
-	pages = {166--171},
-	journaltitle = {Environmental Modelling \& Software},
-	shortjournal = {Environmental Modelling \& Software},
-	author = {Coron, L. and Thirel, G. and Delaigue, O. and Perrin, C. and Andréassian, V.},
-	urldate = {2020-12-18},
-	date = {2017-08},
-	langid = {english}
+@article{coronSuiteLumpedGR2017,
+  title = {The Suite of Lumped {{GR}} Hydrological Models in an {{R}} Package},
+  author = {Coron, L. and Thirel, G. and Delaigue, O. and Perrin, C. and Andr{\'e}assian, V.},
+  year = {2017},
+  month = aug,
+  volume = {94},
+  pages = {166--171},
+  issn = {1364-8152},
+  doi = {10.1016/j.envsoft.2017.05.002},
+  abstract = {Lumped hydrological models are catchment-scale representations of the transformation of precipitation into discharge. They are widely-used tools for real-time flow forecasting, flood design and climate change impact assessment, and they are often used for training and educational purposes. This article presents an R-package, airGR, to facilitate the implementation of the GR lumped hydrological models (including GR4J) and a snow-accumulation and melt model. The package allows users to calibrate and run hourly to annual models on catchment sets and to analyse their outputs. While the core of the models is implemented in Fortran, the user can manage the input/output data within R. A number of options and plotting functions are proposed to ease automate tests and analyses of the results. The codes are flexible enough to include external models, other calibration routines or efficiency criteria. To illustrate the features of airGR, we present one application example for a French mountainous catchment.},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\J94F6IEJ\\Coron et al. - 2017 - The suite of lumped GR hydrological models in an R.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\47JDT6WA\\S1364815217300208.html},
+  journal = {Environmental Modelling \& Software},
+  keywords = {Calibration,Flow simulation,Hydrological modelling,Lumped models,Parsimony,R package},
+  language = {en}
 }
 
-@misc{coron_airgr_2020,
-	title = {{airGR}: Suite of {GR} Hydrological Models for Precipitation-Runoff Modelling. R package version 1.4.3.65.},
-	url = {https://data.inra.fr/citation?persistentId=doi:10.15454/EX11NA},
-	shorttitle = {{airGR}},
-	abstract = {Hydrological modelling tools developed at {INRAE}-Antony ({HYCAR} Research Unit, France). The package includes several conceptual rainfall-runoff models ({GR}4H, {GR}5H, {GR}4J, {GR}5J, {GR}6J, {GR}2M, {GR}1A), a snow accumulation and melt model ({CemaNeige}) and the associated functions for their calibration and evaluation. Use help({airGR}) for package description and references.},
-	publisher = {Portail Data {INRAE}},
-	author = {Coron, Laurent and Delaigue, Olivier and Thirel, Guillaume and Perrin, Charles and Michel, Claude},
-	editora = {Delaigue, Olivier},
-	editoratype = {collaborator},
-	urldate = {2020-12-18},
-	date = {2020},
-	doi = {10.15454/EX11NA},
-	note = {type: dataset}
+@misc{coxonCatchmentAttributesHydrometeorological2020,
+  title = {Catchment Attributes and Hydro-Meteorological Timeseries for 671 Catchments across {{Great Britain}} ({{CAMELS}}-{{GB}})},
+  author = {Coxon, G. and Addor, N. and Bloomfield, J.P. and Freer, J. and Fry, M. and Hannaford, J. and Howden, N.J.K. and Lane, R. and Lewis, M. and Robinson, E.L. and Wagener, T. and Woods, R.},
+  year = {2020},
+  publisher = {{NERC Environmental Information Data Centre}},
+  doi = {10.5285/8344E4F3-D2EA-44F5-8AFA-86D2987543A9},
+  abstract = {This dataset provides hydro-meteorological timeseries and landscape attributes for 671 catchments across Great Britain. It collates river flows, catchment attributes and catchment boundaries from the UK National River Flow Archive together with a suite of new meteorological timeseries and catchment attributes. Daily timeseries for the time period 1st October 1970 to the 30th September 2015 are provided for a range of hydro-meteorological data (including rainfall, potential evapotranspiration, temperature, radiation, humidity and flow). A comprehensive set of catchment attributes are quantified describing a range of catchment characteristics including topography, climate, hydrology, land cover, soils, hydrogeology, human influences and discharge uncertainty. This dataset is intended for the community as a freely available, easily accessible dataset to use in a wide range of environmental data and modelling analyses. A research paper (Coxon et al, CAMELS-GB: Hydrometeorological time series and landscape attributes for 671 catchments in Great Britain) describing the dataset in detail will be made available in Earth System Science Data (https://www.earth-system-science-data.net/).},
+  collaborator = {Coxon, Gemma and Environmental Information Data Centre},
+  copyright = {This resource is available under the terms of the Open Government Licence},
+  keywords = {Hydrology},
+  language = {en}
 }
 
-@article{higgs_hydrological_1988,
-	title = {Hydrological changes and river regulation in the {UK}},
-	volume = {2},
-	issn = {1099-1646},
-	url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/rrr.3450020312},
-	doi = {https://doi.org/10.1002/rrr.3450020312},
-	abstract = {Water levels of streams and rivers in the United Kingdom have been regulated by weirs for more than one thousand years, but regulation of the flow regime by impoundments began in the latter half on the 19th Century. Organized river flow measurements were not undertaken until 1935, and today the average record length is about 20 years. Only three gauging stations have provided data suitable for pre- and post-impoundment comparisons. Other studies have relied on the comparison of regulated and naturalized discharges. In either case climate and land-use changes make evaluation of the hydrological effect of impoundments problematic. This paper reviews research on hydrological changes due to river regulation in the {UK}, and presents a case study of the River Severn to evaluate the influence of Clywedog Reservoir on flood magnitude and frequency. Consequent upon dam completion, on average, median flows have been reduced by about 50per cent; mean annual floods have been reduced by about 30per cent; and low flows have been maintained at about 22 per cent higher than the natural Q95 discharge. However, marked differences exist between rivers. The direct effect of reservoir compensation flows and the indirect effect of inter basin transfers for supply have significantly increased minimum flows in most rivers, although in the case of the latter this involves the discharge of treated effluents. In contrast, the effects of impoundments on flood magnitude and frequency is less clear and on the River Severn, at least, changes in flood hydrology during the past two decades are shown to be more related to climate change than to river regulation.},
-	pages = {349--368},
-	number = {3},
-	journaltitle = {Regulated Rivers: Research \& Management},
-	author = {Higgs, Gary and Petts, Geoff},
-	urldate = {2020-12-26},
-	date = {1988},
-	langid = {english},
-	note = {\_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/rrr.3450020312},
-	keywords = {Floods, Dry weather flow, Flow regime},
-	file = {Higgs et Petts - 1988 - Hydrological changes and river regulation in the U.pdf:C\:\\Users\\david.dorchies\\Zotero\\storage\\VIG8TKGU\\Higgs et Petts - 1988 - Hydrological changes and river regulation in the U.pdf:application/pdf;Snapshot:C\:\\Users\\david.dorchies\\Zotero\\storage\\74LGCRCE\\rrr.html:text/html}
+@article{graftonParadoxIrrigationEfficiency2018,
+  title = {The Paradox of Irrigation Efficiency},
+  author = {Grafton, R. Q. and Williams, J. and Perry, C. J. and Molle, F. and Ringler, C. and Steduto, P. and Udall, B. and Wheeler, S. A. and Wang, Y. and Garrick, D. and Allen, R. G.},
+  year = {2018},
+  month = aug,
+  volume = {361},
+  pages = {748--750},
+  publisher = {{American Association for the Advancement of Science}},
+  issn = {0036-8075, 1095-9203},
+  doi = {10.1126/science.aat9314},
+  abstract = {Higher efficiency rarely reduces water consumption Higher efficiency rarely reduces water consumption},
+  chapter = {Policy Forum},
+  copyright = {Copyright \textcopyright{} 2018, American Association for the Advancement of Science. http://www.sciencemag.org/about/science-licenses-journal-article-reuseThis is an article distributed under the terms of the Science Journals Default License.},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\VIHQIB84\\Grafton et al. - 2018 - The paradox of irrigation efficiency.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\XDKSCJIX\\748.html},
+  journal = {Science},
+  language = {en},
+  number = {6404},
+  pmid = {30139857}
 }
 
-@article{reynard_flood_2001,
-	title = {The Flood Characteristics of Large U.K. Rivers: Potential Effects of Changing Climate and Land Use},
-	volume = {48},
-	issn = {1573-1480},
-	url = {https://doi.org/10.1023/A:1010735726818},
-	doi = {10.1023/A:1010735726818},
-	shorttitle = {The Flood Characteristics of Large U.K. Rivers},
-	abstract = {A continuous flow simulation model({CLASSIC}) has been used to assess the potential impactof climate and land use changes on the flood regimesof large U.K. catchments. Climate change scenarios,based on the {HadCM}2 experiments from the {HadleyCentre}, are applied to the Severn and Thames rivers.The analysis shows that, for the 2050s, the climatechange scenarios result in an increase in both thefrequency and magnitude of flooding events in theserivers. The various ways of applying the rainfallscenario can have a significant effect on thesegeneral conclusions, although generally do not affecteither the direction or consistency of the changes.While ‘best guess’ land use changes show little impacton flood response, a 50\% increase in forest covercould counter-act the impact of climate change. Aswould be expected, a large change in the urban coverof the catchments does have a large effect on theflood regimes, increasing both the frequency andmagnitude of floods significantly beyond the changesdue to climate alone. Further research is requiredinto the potential impacts of seasonal changes in thedaily rainfall and potential evaporation regimes, landuse changes and the interaction between the two.},
-	pages = {343--359},
-	number = {2},
-	journaltitle = {Climatic Change},
-	shortjournal = {Climatic Change},
-	author = {Reynard, N. S. and Prudhomme, C. and Crooks, S. M.},
-	urldate = {2020-12-26},
-	date = {2001-02-01},
-	langid = {english},
-	file = {Reynard et al. - 2001 - The Flood Characteristics of Large U.K. Rivers Po.pdf:C\:\\Users\\david.dorchies\\Zotero\\storage\\6CMCZNJV\\Reynard et al. - 2001 - The Flood Characteristics of Large U.K. Rivers Po.pdf:application/pdf}
+@techreport{gustardStudyCompensationFlows1987,
+  title = {A Study of Compensation Flows in the {{UK}}},
+  author = {Gustard, Alan and Cole, Gwyneth and Marshall, David and Bayliss, Adrian},
+  year = {1987},
+  month = nov,
+  pages = {170},
+  institution = {{Institute of Hydrology}},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\Z3VSBM2Y\\Gustard et al. - 1987 - A study of compensation flows in the UK.pdf}
 }
 
-@article{jolley_large-scale_1996,
-	title = {A Large-Scale Grid-Based Hydrological Model of the Severn and Thames Catchments},
-	volume = {10},
-	issn = {1747-6593},
-	url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1747-6593.1996.tb00043.x},
-	doi = {https://doi.org/10.1111/j.1747-6593.1996.tb00043.x},
-	abstract = {This paper addresses the issues of scale and appropriate model complexity for large-scale hydrological models. A grid-based hydrological model, which employs the {UK} Meteorological Office Rainfall and Evaporation Calculation System, is applied to the Severn and Thames catchments using a grid scale of 40 km, and is shown to reproduce the observed mean annual runoff over a 10-year period to within 6\% with no prior calibration. The variation in the model performance is strongly correlated with the linearity of the annual rainfall/runoff relationship and a climate index. At the monthly scale, runoff routing becomes significant, and the introduction of a two-parameter routeing algorithm significantly improves the monthly runoff simulations giving efficiencies of 90\% and 88\% for the Severn and Thames respectively. The results provide guidance to climate modellers looking for efficient and robust land-surface parameterizations, and indicate the potential application of such a modelling scheme to water resource managers.},
-	pages = {253--262},
-	number = {4},
-	journaltitle = {Water and Environment Journal},
-	author = {Jolley, T. J. and Wheater, H. S.},
-	urldate = {2020-12-26},
-	date = {1996},
-	langid = {english},
-	note = {\_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1747-6593.1996.tb00043.x},
-	keywords = {Climate models, runoff, Severn, Thames, water resources, water-balance model},
-	file = {Jolley et Wheater - 1996 - A Large-Scale Grid-Based Hydrological Model of the.pdf:C\:\\Users\\david.dorchies\\Zotero\\storage\\AKJEZFAY\\Jolley et Wheater - 1996 - A Large-Scale Grid-Based Hydrological Model of the.pdf:application/pdf;Snapshot:C\:\\Users\\david.dorchies\\Zotero\\storage\\4G5L69X5\\j.1747-6593.1996.tb00043.html:text/html}
-}
\ No newline at end of file
+@article{higgsHydrologicalChangesRiver1988,
+  title = {Hydrological Changes and River Regulation in the {{UK}}},
+  author = {Higgs, Gary and Petts, Geoff},
+  year = {1988},
+  volume = {2},
+  pages = {349--368},
+  issn = {1099-1646},
+  doi = {10.1002/rrr.3450020312},
+  abstract = {Water levels of streams and rivers in the United Kingdom have been regulated by weirs for more than one thousand years, but regulation of the flow regime by impoundments began in the latter half on the 19th Century. Organized river flow measurements were not undertaken until 1935, and today the average record length is about 20 years. Only three gauging stations have provided data suitable for pre- and post-impoundment comparisons. Other studies have relied on the comparison of regulated and naturalized discharges. In either case climate and land-use changes make evaluation of the hydrological effect of impoundments problematic. This paper reviews research on hydrological changes due to river regulation in the UK, and presents a case study of the River Severn to evaluate the influence of Clywedog Reservoir on flood magnitude and frequency. Consequent upon dam completion, on average, median flows have been reduced by about 50per cent; mean annual floods have been reduced by about 30per cent; and low flows have been maintained at about 22 per cent higher than the natural Q95 discharge. However, marked differences exist between rivers. The direct effect of reservoir compensation flows and the indirect effect of inter basin transfers for supply have significantly increased minimum flows in most rivers, although in the case of the latter this involves the discharge of treated effluents. In contrast, the effects of impoundments on flood magnitude and frequency is less clear and on the River Severn, at least, changes in flood hydrology during the past two decades are shown to be more related to climate change than to river regulation.},
+  annotation = {\_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/rrr.3450020312},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\Q4EWWR7A\\Higgs et Petts - 1988 - Hydrological changes and river regulation in the U.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\CI9TWW3V\\rrr.html},
+  journal = {Regulated Rivers: Research \& Management},
+  keywords = {Dry weather flow,Floods,Flow regime},
+  language = {en},
+  number = {3}
+}
+
+@article{jolleyLargeScaleGridBasedHydrological1996,
+  title = {A {{Large}}-{{Scale Grid}}-{{Based Hydrological Model}} of the {{Severn}} and {{Thames Catchments}}},
+  author = {Jolley, T. J. and Wheater, H. S.},
+  year = {1996},
+  volume = {10},
+  pages = {253--262},
+  issn = {1747-6593},
+  doi = {10.1111/j.1747-6593.1996.tb00043.x},
+  abstract = {This paper addresses the issues of scale and appropriate model complexity for large-scale hydrological models. A grid-based hydrological model, which employs the UK Meteorological Office Rainfall and Evaporation Calculation System, is applied to the Severn and Thames catchments using a grid scale of 40 km, and is shown to reproduce the observed mean annual runoff over a 10-year period to within 6\% with no prior calibration. The variation in the model performance is strongly correlated with the linearity of the annual rainfall/runoff relationship and a climate index. At the monthly scale, runoff routing becomes significant, and the introduction of a two-parameter routeing algorithm significantly improves the monthly runoff simulations giving efficiencies of 90\% and 88\% for the Severn and Thames respectively. The results provide guidance to climate modellers looking for efficient and robust land-surface parameterizations, and indicate the potential application of such a modelling scheme to water resource managers.},
+  annotation = {\_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1747-6593.1996.tb00043.x},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\NIH9WIPD\\Jolley et Wheater - 1996 - A Large-Scale Grid-Based Hydrological Model of the.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\2T597AZG\\j.1747-6593.1996.tb00043.html},
+  journal = {Water and Environment Journal},
+  keywords = {Climate models,runoff,Severn,Thames,water resources,water-balance model},
+  language = {en},
+  number = {4}
+}
+
+@article{klaarDevelopingHydroecologicalModels2014,
+  title = {Developing Hydroecological Models to Inform Environmental Flow Standards: A Case Study from {{England}}},
+  shorttitle = {Developing Hydroecological Models to Inform Environmental Flow Standards},
+  author = {Klaar, Megan J. and Dunbar, Michael J. and Warren, Mark and Soley, Rob},
+  year = {2014},
+  volume = {1},
+  pages = {207--217},
+  issn = {2049-1948},
+  doi = {10.1002/wat2.1012},
+  abstract = {The concept of defining environmental flow regimes to balance the provision of water resources for both human and environmental needs has gained wide recognition. As the authority responsible for water resource management within England, the Environment Agency (EA) uses the Environmental Flow Indicator (EFI), which represents an allowable percentage deviation from the natural flow to determine where water may be available for new abstractions. In a simplified form, the EFI has been used as the hydrological supporting component of Water Framework Directive classification, to flag where hydrological alteration may be contributing to failure to achieve good ecological status, and to guide further ecological investigation. As the primary information source for the EFI was expert opinion, the EA aims to improve the evidence base linking flow alteration and ecological response, and to use this evidence to develop improved environmental flow criteria and implementation tools. Such tools will be required to make predictions at locations with no or limited ecological monitoring data. Hence empirical statistical models are required that provide a means to describe observed variation in ecological sensitivity to flow change. Models must also strike a balance between generic and local relationships. Multilevel (mixed effects) regression models provide a rich set of capabilities suitable for this purpose. Three brief examples of the application of these techniques in defining empirical relationships between flow alteration and ecological response are provided. Establishment of testable hydrological\textendash ecological relationships provides the framework for improving data collection, analysis, and ultimately water resources management models. This article is categorized under: Water and Life {$>$} Conservation, Management, and Awareness},
+  annotation = {\_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/wat2.1012},
+  copyright = {\textcopyright{} 2014 Wiley Periodicals, Inc.},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\RNDY5IW9\\wat2.html},
+  journal = {WIREs Water},
+  language = {en},
+  number = {2}
+}
+
+@article{reynardFloodCharacteristicsLarge2001,
+  title = {The {{Flood Characteristics}} of {{Large U}}.{{K}}. {{Rivers}}: {{Potential Effects}} of {{Changing Climate}} and {{Land Use}}},
+  shorttitle = {The {{Flood Characteristics}} of {{Large U}}.{{K}}. {{Rivers}}},
+  author = {Reynard, N. S. and Prudhomme, C. and Crooks, S. M.},
+  year = {2001},
+  month = feb,
+  volume = {48},
+  pages = {343--359},
+  issn = {1573-1480},
+  doi = {10.1023/A:1010735726818},
+  abstract = {A continuous flow simulation model(CLASSIC) has been used to assess the potential impactof climate and land use changes on the flood regimesof large U.K. catchments. Climate change scenarios,based on the HadCM2 experiments from the HadleyCentre, are applied to the Severn and Thames rivers.The analysis shows that, for the 2050s, the climatechange scenarios result in an increase in both thefrequency and magnitude of flooding events in theserivers. The various ways of applying the rainfallscenario can have a significant effect on thesegeneral conclusions, although generally do not affecteither the direction or consistency of the changes.While `best guess' land use changes show little impacton flood response, a 50\% increase in forest covercould counter-act the impact of climate change. Aswould be expected, a large change in the urban coverof the catchments does have a large effect on theflood regimes, increasing both the frequency andmagnitude of floods significantly beyond the changesdue to climate alone. Further research is requiredinto the potential impacts of seasonal changes in thedaily rainfall and potential evaporation regimes, landuse changes and the interaction between the two.},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\PZTI9ZQE\\Reynard et al. - 2001 - The Flood Characteristics of Large U.K. Rivers Po.pdf},
+  journal = {Climatic Change},
+  language = {en},
+  number = {2}
+}
+
+@incollection{secklerConceptEfficiencyWaterresources2003,
+  title = {The Concept of Efficiency in Water-Resources Management and Policy.},
+  booktitle = {Water Productivity in Agriculture: Limits and Opportunities for Improvement},
+  author = {Seckler, D. and Molden, D. and Sakthivadivel, R.},
+  editor = {Kijne, J. W. and Barker, R. and Molden, D.},
+  year = {2003},
+  pages = {37--51},
+  publisher = {{CABI}},
+  address = {{Wallingford}},
+  doi = {10.1079/9780851996691.0037},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\FM69CIK3\\Seckler et al. - 2003 - The concept of efficiency in water-resources manag.pdf},
+  isbn = {978-0-85199-669-1},
+  language = {en}
+}
+
+