diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index be71d890fc117ec75e3a6c0461bd937eb7911ef8..dea36093c0ba94564c16d8e721963847bc2751b8 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -1,7 +1,7 @@
 RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
   NParam <- 1
 
-  ##Arguments_check
+  ## argument check
   if (!inherits(InputsModel, "InputsModel")) {
     stop("'InputsModel' must be of class 'InputsModel'")
   }
@@ -21,35 +21,49 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
     if (is.null(QcontribDown$Qsim)) {
       stop("'QcontribDown' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
     }
-    OutputsModel <- QcontribDown
-    OutputsModel$QsimDown <- OutputsModel$Qsim
+    if (length(QcontribDown$Qsim) != length(RunOptions$IndPeriod_Run)) {
+      stop("Time series Qsim in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
+    }
+    if (!is.null(QcontribDown$WarmUpQsim) &&
+        length(QcontribDown$WarmUpQsim) != length(RunOptions$IndPeriod_WarmUp) &&
+        RunOptions$IndPeriod_WarmUp != 0L) {
+      stop("Time series WarmUpQsim in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_WarmUp'")
+    }
   } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
-    OutputsModel <- list()
-    class(OutputsModel) <- c("OutputsModel", class(OutputsModel))
-    OutputsModel$QsimDown <- QcontribDown
+    if (length(QcontribDown) != length(RunOptions$IndPeriod_Run)) {
+      stop("'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
+    }
   } else {
     stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
   }
-  if (length(OutputsModel$QsimDown) != length(RunOptions$IndPeriod_Run)) {
-    stop("Time series in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
+
+  # data set up
+  NbUpBasins <- length(InputsModel$LengthHydro)
+  if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
+    RunOptions$IndPeriod_WarmUp <- NULL
   }
+  IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
+  LInputSeries <- as.integer(length(IndPeriod1))
+  IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
 
-  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\"")
+  if (inherits(QcontribDown, "OutputsModel")) {
+    OutputsModel <- QcontribDown
+    if (is.null(OutputsModel$WarmUpQsim)) {
+      OutputsModel$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
+    }
+    QsimDown <- c(OutputsModel$WarmUpQsim, OutputsModel$Qsim)
+  } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
+    OutputsModel <- list()
+    class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
+    QsimDown <- c(rep(NA, length(RunOptions$IndPeriod_WarmUp)),
+                  QcontribDown)
   }
 
-  # propagation time from upstream meshes to outlet
-  PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / TimeStep
+  ## propagation time from upstream meshes to outlet
+  PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / RunOptions$FeatFUN_MOD$TimeStep
   HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
 
-  NbUpBasins <- length(InputsModel$LengthHydro)
-  LengthTs <- length(OutputsModel$QsimDown)
-  OutputsModel$Qsim_m3 <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
-
+  ## set up initial states
   IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
   if (length(IniSD) > 0) {
     if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
@@ -66,15 +80,15 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
       if (x > 1) {
         iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
       }
-      IniSD[iStart:(iStart + PT[x])]
+      as.vector(IniSD[iStart:(iStart + PT[x])])
     })
   } else {
     IniStates <- lapply(
       seq_len(NbUpBasins),
       function(iUpBasins) {
         iWarmUp <- seq.int(
-          from = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)] - floor(PT[iUpBasins])),
-          to   = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)])
+          from = max(1, IndPeriod1[1] - floor(PT[iUpBasins]) - 1),
+          to   = max(1, IndPeriod1[1] - 1)
         )
         ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
         if (length(ini) != floor(PT[iUpBasins] + 1)) {
@@ -85,19 +99,30 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
       }
     )
   }
-  # message("Initstates: ", paste(IniStates, collapse = ", "))
+  # message("IniStates: ", paste(IniStates, collapse = ", "))
+
+  ## Lag model computation
+  Qsim_m3 <- QsimDown *
+    InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
 
   for (upstream_basin in seq_len(NbUpBasins)) {
     Qupstream <- c(IniStates[[upstream_basin]],
-                   InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin])
+                   InputsModel$Qupstream[IndPeriod1, upstream_basin])
     # message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
-    OutputsModel$Qsim_m3 <- OutputsModel$Qsim_m3 +
-      Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
-      Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
+    Qsim_m3 <- Qsim_m3 +
+      Qupstream[2:(1 + LInputSeries)] * HUTRANS[1, upstream_basin] +
+      Qupstream[1:LInputSeries] * HUTRANS[2, upstream_basin]
+  }
+
+  ## OutputsModel
+
+  OutputsModel$Qsim_m3 <- Qsim_m3[IndPeriod2]
+
+  if ("Qsim" %in% RunOptions$Outputs_Sim) {
+    # Convert back Qsim to mm
+    OutputsModel$Qsim <- OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+    # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
   }
-  # Convert back Qsim to mm
-  OutputsModel$Qsim <- OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
-  # message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
 
   # Warning for negative flows or NAs only in extended outputs
   if (length(RunOptions$Outputs_Sim) > 2) {
@@ -112,12 +137,20 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
   }
 
   if ("StateEnd" %in% RunOptions$Outputs_Sim) {
-    OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
+    SD <- lapply(seq(NbUpBasins), function(x) {
       lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
       InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
     })
+    if(is.null(OutputsModel$StateEnd)) {
+      OutputsModel$StateEnd <- CreateIniStates(RunModel_Lag, InputsModel, SD = SD)
+    } else {
+      OutputsModel$StateEnd$SD <- SD
+    }
     # message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
   }
+  if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))]
+  }
 
   if ("Param" %in% RunOptions$Outputs_Sim) {
     OutputsModel$Param <- c(Param, OutputsModel$Param)