diff --git a/.gitignore b/.gitignore
index 72342a65dbae583572a67520e9566b7e5bd68583..188eac083f4542383a37a30f2515c869bc20cc09 100644
--- a/.gitignore
+++ b/.gitignore
@@ -34,7 +34,7 @@ packrat/lib*/
 /*.Rcheck/
 
 # RStudio files
-.Rproj.user/
+.Rproj.user
 
 # produced vignettes
 vignettes/*.html
diff --git a/R/CreateInputsModel.R b/R/CreateInputsModel.R
index 34c85db755340d68940bffc4fad0787e4a36aed9..c55c928265642d630c60d0e033f99a0f9c988d50 100644
--- a/R/CreateInputsModel.R
+++ b/R/CreateInputsModel.R
@@ -207,6 +207,9 @@ CreateInputsModel <- function(FUN_MOD,
       if (ncol(QobsUpstr)+1 != ncol(BasinAreas)) {
         stop("'BasinAreas' must have one column more than 'QobsUpstr' and 'LengthHydro'")
       }
+      if (nrow(QobsUpstr) != LLL) {
+        stop("'QobsUpstr' must have same number of rows as 'DatesR' length")
+      }
       if (nrow(LengthHydro) != 1 | nrow(BasinAreas) != 1) {
         stop("'LengthHydro' and 'BasinAreas' must have only one row")
       }
diff --git a/R/RunModel.R b/R/RunModel.R
index cd3013d9388c1c1758c27cba1ec2d0738cdd449b..496a3d533d870e5b609facbbf3dd3502e6fb31fd 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -26,12 +26,13 @@ RunModel <- function (InputsModel, RunOptions, Param, FUN_MOD) {
     OutputsModelDown$Qsim <- OutputsModelDown$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] / AreaTot
     
     for (upstream_basin in seq_len(NbUpBasins)) {
+      QobsUpstr <- InputsModel$QobsUpstr[RunOptions$IndPeriod_Run, upstream_basin]
       OutputsModelDown$Qsim <- OutputsModelDown$Qsim + 
         c(rep(0, floor(PT[upstream_basin])), 
-          InputsModel$QobsUpstr[(1 + floor(PT[upstream_basin])):LengthTs, upstream_basin]) *
+          QobsUpstr[(1 + floor(PT[upstream_basin])):LengthTs]) *
         HUTRANS[1, upstream_basin] * InputsModel$BasinAreas[upstream_basin] / AreaTot + 
         c(rep(0, floor(PT[upstream_basin] + 1)), 
-          InputsModel$QobsUpstr[(2 + floor(PT[upstream_basin])):LengthTs, upstream_basin]) *
+          QobsUpstr[(2 + floor(PT[upstream_basin])):LengthTs]) *
         HUTRANS[2, upstream_basin] * InputsModel$BasinAreas[upstream_basin] / AreaTot
     }
     
diff --git a/tests/testthat/test-RunModel_LAG.R b/tests/testthat/test-RunModel_LAG.R
new file mode 100644
index 0000000000000000000000000000000000000000..9254cb4bee4a885d166175a2a54ba7bbc71391b7
--- /dev/null
+++ b/tests/testthat/test-RunModel_LAG.R
@@ -0,0 +1,85 @@
+context("RunModel_LAG")
+
+data(L0123001)
+
+
+test_that("'BasinAreas' must have one column more than 'QobsUpstr' and 'LengthHydro'", {
+    expect_error(
+        InputsModel <- CreateInputsModel(
+            FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+            Precip = BasinObs$P, PotEvap = BasinObs$E,
+            QobsUpstr = matrix(BasinObs$Qmm * 2,ncol=1),
+            LengthHydro = matrix(c(1),nrow=1),
+            BasinAreas = matrix(c(1),nrow=1)
+        ),
+        regexp = "'BasinAreas' must have one column more than 'QobsUpstr' and 'LengthHydro'"
+    )
+})
+
+test_that("'BasinAreas' must have one column more than 'QobsUpstr' and 'LengthHydro'", {
+    expect_error(
+        InputsModel <- CreateInputsModel(
+            FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+            Precip = BasinObs$P, PotEvap = BasinObs$E,
+            QobsUpstr = matrix(BasinObs$Qmm, ncol = 1),
+            LengthHydro = matrix(c(1), nrow = 1),
+            BasinAreas = matrix(c(1, 2),nrow = 1)
+        ),
+        regexp = "'QobsUpstr' cannot contain any NA value"
+    )
+})
+
+QobsUpstr = BasinObs$Qmm
+QobsUpstr[is.na(QobsUpstr)] = mean(QobsUpstr, na.rm = TRUE)
+
+InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P, PotEvap = BasinObs$E,
+    QobsUpstr = matrix(QobsUpstr,ncol=1),
+    LengthHydro = matrix(c(1),nrow=1),
+    BasinAreas = matrix(c(BasinInfo$BasinArea,BasinInfo$BasinArea),nrow=1)
+)
+
+Ind_Run <- seq(
+    which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1990-01-01"),
+    which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="1999-12-31")
+)
+
+RunOptions <- CreateRunOptions(
+    FUN_MOD = RunModel_GR4J,
+    InputsModel = InputsModel, IndPeriod_Run = Ind_Run
+)
+
+Param = c(257.237556, 1.012237, 88.234673, 2.207958) #  From vignettes/V01_get_started
+
+OutputsGR4JOnly <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
+
+
+test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
+  UpstBasinArea = InputsModel$BasinAreas[1,1]
+  InputsModel$BasinAreas[1,1] <<- 0
+  OutputsSD <- RunModel(InputsModel, RunOptions, Param = c(Param, 1), FUN_MOD = RunModel_GR4J)
+  expect_equal(OutputsGR4JOnly$Qsim, OutputsSD$Qsim)
+  InputsModel$BasinAreas[1,1] <<- UpstBasinArea
+})
+
+test_that("Downstream basin with nil area and nul upstream length should return same Qdown as Qupstream alone", {
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P, PotEvap = BasinObs$E,
+    QobsUpstr = matrix(QobsUpstr,ncol=1),
+    LengthHydro = matrix(c(0),nrow=1),
+    BasinAreas = matrix(c(BasinInfo$BasinArea,0),nrow=1)
+  )
+  OutputsSD <- RunModel(InputsModel, RunOptions, Param = c(Param, 1), FUN_MOD = RunModel_GR4J)
+  expect_equal(OutputsSD$Qsim, QobsUpstr[Ind_Run])
+})
+
+test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
+    QlsGR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2] * 1E6 / 86400
+    ParamSD = c(Param, InputsModel$LengthHydro * 1000 / (24 * 60 * 60))
+    OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
+    QlsSdSim <- OutputsSD$Qsim * sum(InputsModel$BasinAreas) * 1E6 / 86400
+    QlsUpstLagObs <- c(0, QobsUpstr[Ind_Run[1:(length(Ind_Run)-1)] + 1]) * InputsModel$BasinAreas[1] * 1E6 / 86400
+    expect_equal(QlsSdSim - QlsGR4Only, QlsUpstLagObs)
+})