diff --git a/NAMESPACE b/NAMESPACE
index c63da5247d87d1bb4c397d034ff35faf24281f0e..8348ba43e310606ba7183ff78740ef1d789dbab3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,79 +1,81 @@
-#####################################
-##              Load DLL           ##
-#####################################
-useDynLib(airGR, .registration = TRUE)
-
-
-
-#####################################
-##            S3 methods           ##
-#####################################
-S3method('[', InputsModel)
-#S3method('[', OutputsModel) ### to add in version 2.0
-S3method(plot, OutputsModel)
-S3method(SeriesAggreg, data.frame)
-S3method(SeriesAggreg, list)
-S3method(SeriesAggreg, InputsModel)
-S3method(SeriesAggreg, OutputsModel)
-
-
-
-#####################################
-##               Export            ##
-#####################################
-export(Calibration)
-export(Calibration_Michel)
-export(CreateCalibOptions)
-export(CreateErrorCrit_GAPX)
-export(CreateIniStates)
-export(CreateInputsCrit)
-export(CreateInputsCrit_Lavenne)
-export(CreateInputsModel)
-export(CreateRunOptions)
-export(DataAltiExtrapolation_Valery)
-export(ErrorCrit)
-export(ErrorCrit_KGE)
-export(ErrorCrit_KGE2)
-export(ErrorCrit_NSE)
-export(ErrorCrit_RMSE)
-export(Imax)
-export(PE_Oudin)
-export(plot.OutputsModel) ### to remove from version 2.0
-export(RunModel)
-export(RunModel_CemaNeige)
-export(RunModel_CemaNeigeGR4H)
-export(RunModel_CemaNeigeGR5H)
-export(RunModel_CemaNeigeGR4J)
-export(RunModel_CemaNeigeGR5J)
-export(RunModel_CemaNeigeGR6J)
-export(RunModel_GR1A)
-export(RunModel_GR2M)
-export(RunModel_GR4H)
-export(RunModel_GR5H)
-export(RunModel_GR4J)
-export(RunModel_GR5J)
-export(RunModel_GR6J)
-export(RunModel_Lag)
-export(SeriesAggreg)
-export(TransfoParam)
-export(TransfoParam_CemaNeige)
-export(TransfoParam_CemaNeigeHyst)
-export(TransfoParam_GR1A)
-export(TransfoParam_GR2M)
-export(TransfoParam_GR4H)
-export(TransfoParam_GR5H)
-export(TransfoParam_GR4J)
-export(TransfoParam_GR5J)
-export(TransfoParam_GR6J)
-export(TransfoParam_Lag)
-export(.ErrorCrit)
-export(.FeatModels)
-
-
-#####################################
-##               Import            ##
-#####################################
-import(stats)
-import(graphics)
-import(grDevices)
-import(utils)
+#####################################
+##              Load DLL           ##
+#####################################
+useDynLib(airGR, .registration = TRUE)
+
+
+
+#####################################
+##            S3 methods           ##
+#####################################
+S3method('[', InputsModel)
+#S3method('[', OutputsModel) ### to add in version 2.0
+S3method(plot, OutputsModel)
+S3method(SeriesAggreg, data.frame)
+S3method(SeriesAggreg, list)
+S3method(SeriesAggreg, InputsModel)
+S3method(SeriesAggreg, OutputsModel)
+
+
+
+#####################################
+##               Export            ##
+#####################################
+export(Calibration)
+export(Calibration_Michel)
+export(CreateCalibOptions)
+export(CreateErrorCrit_GAPX)
+export(CreateIniStates)
+export(CreateInputsCrit)
+export(CreateInputsCrit_Lavenne)
+export(CreateInputsModel)
+export(CreateRunOptions)
+export(DataAltiExtrapolation_Valery)
+export(ErrorCrit)
+export(ErrorCrit_KGE)
+export(ErrorCrit_KGE2)
+export(ErrorCrit_NSE)
+export(ErrorCrit_RMSE)
+export(Imax)
+export(PE_Oudin)
+export(plot.OutputsModel) ### to remove from version 2.0
+export(RunModel)
+export(RunModel_CemaNeige)
+export(RunModel_CemaNeigeGR4H)
+export(RunModel_CemaNeigeGR5H)
+export(RunModel_CemaNeigeGR4J)
+export(RunModel_CemaNeigeGR5J)
+export(RunModel_CemaNeigeGR6J)
+export(RunModel_GR1A)
+export(RunModel_GR2M)
+export(RunModel_GR4H)
+export(RunModel_GR5H)
+export(RunModel_GR4J)
+export(RunModel_GR5J)
+export(RunModel_GR6J)
+export(RunModel_Lag)
+export(RunModel_LLR)
+export(SeriesAggreg)
+export(TransfoParam)
+export(TransfoParam_CemaNeige)
+export(TransfoParam_CemaNeigeHyst)
+export(TransfoParam_GR1A)
+export(TransfoParam_GR2M)
+export(TransfoParam_GR4H)
+export(TransfoParam_GR5H)
+export(TransfoParam_GR4J)
+export(TransfoParam_GR5J)
+export(TransfoParam_GR6J)
+export(TransfoParam_Lag)
+export(TransfoParam_LLR)
+export(.ErrorCrit)
+export(.FeatModels)
+
+
+#####################################
+##               Import            ##
+#####################################
+import(stats)
+import(graphics)
+import(grDevices)
+import(utils)
diff --git a/R/Calibration_Michel.R b/R/Calibration_Michel.R
index c7119cba1a00fbcef54e38d7a499b46418a28f75..b5ae1e00f6c4247bcf04b415f24664077e0937f5 100644
--- a/R/Calibration_Michel.R
+++ b/R/Calibration_Michel.R
@@ -296,7 +296,7 @@ Calibration_Michel <- function(InputsModel,
     for (iNew in 1:nrow(CandidatesParamR)) {
       ##Model_run
       Param <- CandidatesParamR[iNew, ]
-      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, FUN_SD = RunOptions$FUN_SD, ...)
       ##Calibration_criterion_computation
       OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
       if (!is.na(OutputsCrit$CritValue)) {
@@ -357,7 +357,7 @@ Calibration_Michel <- function(InputsModel,
       CandidatesParamR <- FUN_TRANSFO(CandidatesParamT, "TR")
       ##Model_run
       Param <- CandidatesParamR[iNew, ]
-      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, ...)
+      OutputsModel <- RunModel(InputsModel, RunOptions, Param, FUN_MOD = FUN_MOD, FUN_SD = RunOptions$FUN_SD, ...)
       ##Calibration_criterion_computation
       OutputsCrit <- ErrorCrit(InputsCrit, OutputsModel, verbose = FALSE)
       if (OutputsCrit$CritValue * OutputsCrit$Multiplier < CritOptim) {
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index a7b0810feaecf4b4737a0d661b0d8fe9cddbb6cb..a33384b267a2fcc4fb3a00ed94ad9fe5ecd71fa5 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -205,7 +205,7 @@ CreateCalibOptions <- function(FUN_MOD,
         ParamT <- cbind(ParamTSD, ParamT)
       }
       if (identical(RunModel_LLR, FUN_SD)){
-        ParamTSD <- matrix(c(+1.25, +1.25
+        ParamTSD <- matrix(c(+1.25, +1.25,
                              +2.50, +2.50,
                              +5.00, +5.00), ncol = 2, byrow = TRUE)
         ParamT <- cbind(ParamTSD, ParamT)
diff --git a/R/RunModel.R b/R/RunModel.R
index 685cfbe9f2c96f731e7372a2f5101f691da4a1b9..be495e067343c210f489ecd7c356f49b4834b296 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -14,7 +14,7 @@ RunModel <- function(InputsModel, RunOptions, Param, FUN_MOD, FUN_SD = RunModel_
                           Param = Param[iFirstParamRunOffModel:length(Param)], ...)
 
   if (inherits(InputsModel, "SD") && !identical(FUN_MOD, FUN_SD)) {
-    OutputsModel <- RunModel_Lag(InputsModel, RunOptions, Param[1:RunOptions$FeatFUN_MOD$NParamSD],
+    OutputsModel <- FUN_SD(InputsModel, RunOptions, Param[1:RunOptions$FeatFUN_MOD$NParamSD],
                                  OutputsModel)
   }
   return(OutputsModel)
diff --git a/R/RunModel_LLR.R b/R/RunModel_LLR.R
index f5dff3161c35775e4d7aff98ca6be4b561fcea42..7c0f6eaf9a65be2ce3fe270a0dab82f6cd22bb34 100644
--- a/R/RunModel_LLR.R
+++ b/R/RunModel_LLR.R
@@ -63,71 +63,38 @@ RunModel_LLR <- function(InputsModel, RunOptions, Param, QcontribDown) {
                   QcontribDown)
   }
 
+
   ## Parameters set up
 
   TParam <- Param[1L]
   KParam <- Param[2L]
 
-  C0 <- exp(-1/KParam)
-  C1 <- ((KParam/1) * (1-C0))-C0
-  C2 <- 1 - (KParam/1) * (1-C0)
-
-
   ## propagation time from upstream meshes to outlet
-  PT <- InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep
-  print(PT)
-  HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
-  ## set up initial states
-  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)
-      }
-      as.vector(IniSD[iStart:(iStart + PT[x])])
-    })
-  } else {
-    IniStates <- lapply(
-      seq_len(NbUpBasins),
-      function(iUpBasins) {
-        iWarmUp <- seq.int(
-          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)) {
-          # If warm-up period is not enough long complete beginning with first value
-          ini <- c(rep(ini[1], floor(PT[iUpBasins] + 1) - length(ini)), ini)
-        }
-        return(as.vector(ini))
-      }
-    )
-  }
-  # message("IniStates: ", paste(IniStates, collapse = ", "))
+  PT <- round(InputsModel$LengthHydro * 1e3 / TParam / RunOptions$FeatFUN_MOD$TimeStep)
+  PK <- sqrt(InputsModel$LengthHydro/max(InputsModel$LengthHydro)) * KParam
+
+  C0 <- exp(-1/PK)
+  C1 <- ((PK/1) * (1-C0))-C0
+  C2 <- 1 - (PK/1) * (1-C0)
 
   ## Lag model computation
-  Qsim_m3 <- QsimDown *
-    InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+  Qsim_m3 <- matrix(QsimDown *
+                      InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+                    , nrow = length(IndPeriod1), ncol = NbUpBasins)
+
 
   for (upstream_basin in seq_len(NbUpBasins)) {
-    Qupstream <- c(IniStates[[upstream_basin]],
-                   InputsModel$Qupstream[IndPeriod1, upstream_basin])
-    # message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
-    Qsim_m3 <- Qsim_m3 +
-      Qupstream[2:(1 + LInputSeries)] * HUTRANS[1, upstream_basin] +
-      Qupstream[1:LInputSeries] * HUTRANS[2, upstream_basin]
+    Qupstream_m3 <- InputsModel$Qupstream[IndPeriod1, upstream_basin]
+    Qroute <- Qllr <- Qupstream_m3
+    for (q_upstream in seq_along(Qupstream_m3)[2:max(seq_along(Qupstream_m3))]){
+      Qroute[q_upstream] <- C0 * Qroute[q_upstream - 1] + C1 * Qupstream_m3[q_upstream - 1] + C2 * Qupstream_m3[q_upstream]
+      Qllr[q_upstream + PT[upstream_basin]] <- Qroute[q_upstream]
+    }
+    Qsim_m3[, upstream_basin] <- Qllr[IndPeriod1]
   }
 
+
+  Qsim_m3 <- rowSums(Qsim_m3)
   ## OutputsModel
 
   if ("Qsim_m3" %in% RunOptions$Outputs_Sim) {
diff --git a/R/TransfoParam_LLR.R b/R/TransfoParam_LLR.R
index d791f8e4f11be93350ffa4a1f782f1a34bdc95da..675201e8dd746d59e6e667e1e5f6cccdc5ebc039 100644
--- a/R/TransfoParam_LLR.R
+++ b/R/TransfoParam_LLR.R
@@ -2,6 +2,7 @@ TransfoParam_LLR <- function(ParamIn, Direction) {
 
   ## number of model parameters
   NParam <- 2L
+
   ## check arguments
   isVecParamIn <- is.vector(ParamIn)
   if (isVecParamIn) {
@@ -20,10 +21,12 @@ TransfoParam_LLR <- function(ParamIn, Direction) {
 
   ## transformation
   if (Direction == "TR") {
+    ParamOut <- ParamIn
     ParamOut[, 1] <- 20 * (ParamIn[, 1] + 10) / 20.0
     ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
   }
   if (Direction == "RT") {
+    ParamOut <- ParamIn
     ParamOut[, 1] <- ParamIn[, 1] * 20.0 / 20 - 10
     ParamOut[, 2] <- 20 * (ParamIn[, 2] + 10) / 20.0
   }
diff --git a/R/TransfoParam_Lag.R b/R/TransfoParam_Lag.R
index 0ee13b0462e535df8da47ca60452115c610bc1f7..b286ef17ffc804129e913b80a3dd1234699add27 100644
--- a/R/TransfoParam_Lag.R
+++ b/R/TransfoParam_Lag.R
@@ -1,39 +1,38 @@
-TransfoParam_Lag <- function(ParamIn, Direction) {
-  
-  ## number of model parameters
-  NParam <- 1L
-  
-  
-  ## check arguments
-  isVecParamIn <- is.vector(ParamIn)
-  if (isVecParamIn) {
-    ParamIn <- matrix(ParamIn, nrow = 1)
-  }  
-  if (!inherits(ParamIn, "matrix")) {
-    stop("'ParamIn' must be of class 'matrix'")
-  }
-  if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
-    stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
-  }
-  if (ncol(ParamIn) != NParam) {
-    stop(sprintf("the Lag model requires %i parameters", NParam))
-  }
-  
-  
-  ## transformation
-  if (Direction == "TR") {
-    ParamOut <- 20 * (ParamIn + 10) / 20.0
-  }
-  if (Direction == "RT") {
-    ParamOut <- ParamIn * 20.0 / 20 - 10 
-  }
-  
-  
-  ## check output
-  if (isVecParamIn) {
-    ParamOut <- as.vector(ParamOut)
-  }
-  
-  return(ParamOut)
-  
-}
+TransfoParam_Lag <- function(ParamIn, Direction) {
+
+  ## number of model parameters
+  NParam <- 1L
+
+
+  ## check arguments
+  isVecParamIn <- is.vector(ParamIn)
+  if (isVecParamIn) {
+    ParamIn <- matrix(ParamIn, nrow = 1)
+  }
+  if (!inherits(ParamIn, "matrix")) {
+    stop("'ParamIn' must be of class 'matrix'")
+  }
+  if (!inherits(Direction, "character") | length(Direction) != 1 | any(!Direction %in% c("RT", "TR"))) {
+    stop("'Direction' must be a character vector of length 1 equal to 'RT' or 'TR'")
+  }
+  if (ncol(ParamIn) != NParam) {
+    stop(sprintf("the Lag model requires %i parameters", NParam))
+  }
+
+  ## transformation
+  if (Direction == "TR") {
+    ParamOut <- 20 * (ParamIn + 10) / 20.0
+  }
+  if (Direction == "RT") {
+    ParamOut <- ParamIn * 20.0 / 20 - 10
+  }
+
+
+  ## check output
+  if (isVecParamIn) {
+    ParamOut <- as.vector(ParamOut)
+  }
+
+  return(ParamOut)
+
+}
diff --git a/R/Utils.R b/R/Utils.R
index f3716df2bbdf5945ed4a8e22ab45dde021ac5c0b..87b91ce17a5458690d7fdd290f8ab39c8c265fa0 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -76,9 +76,10 @@
     if (!is.null(FUN_SD)){
       if (identical(RunModel_Lag, FUN_SD)){
         res$NParamSD <- FeatMod$NbParam[14]
-      }
-      if(identical(RunModel_LLR, FUN_SD)){
+        res$FUN_SD <- "Lag"
+      }else if(identical(RunModel_LLR, FUN_SD)){
         res$NParamSD <- FeatMod$NbParam[15]
+        res$FUN_SD <- "LLR"
       }
     }
     return(res)
diff --git a/R/UtilsCalibOptions.R b/R/UtilsCalibOptions.R
index a9ecfe084e28ad51b7738a62c5aa1e6e7dbe119e..1454498504572eb938fff0da3d199b81bca4af64 100644
--- a/R/UtilsCalibOptions.R
+++ b/R/UtilsCalibOptions.R
@@ -29,13 +29,11 @@
   if (IsSD) {
     if (identical(RunModel_Lag, FUN_SD)){
       FUN_ROUT <- TransfoParam_Lag
-    }
-    if (identical(RunModel_LLR, FUN_SD)){
+    } else if (identical(RunModel_LLR, FUN_SD)){
       FUN_ROUT <- TransfoParam_LLR
     }
   }
 
-
     ## set FUN_TRANSFO
   if (! "CemaNeige" %in% FeatFUN_MOD$Class) {
     if (!IsSD) {
@@ -46,10 +44,17 @@
         if (!Bool) {
           ParamIn <- rbind(ParamIn)
         }
+
         ParamOut <- NA * ParamIn
         NParam   <- ncol(ParamIn)
         ParamOut[, (FeatFUN_MOD$NParamSD + 1):NParam] <- FUN_GR(ParamIn[, (FeatFUN_MOD$NParamSD + 1):NParam], Direction)
-        ParamOut[, FeatFUN_MOD$NParamSD       ] <- FUN_ROUT(as.matrix(ParamIn[, FeatFUN_MOD$NParamSD]), Direction)
+
+        if (identical(RunModel_Lag, FUN_SD)){
+          ParamOut[, 1       ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
+        } else if (identical(RunModel_LLR, FUN_SD)){
+          ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
+        }
+
         if (!Bool) {
           ParamOut <- ParamOut[1, ]
         }
@@ -103,7 +108,7 @@
         NParam   <- ncol(ParamIn)
         ParamOut[, 2:(NParam - 4)     ] <- FUN_GR(ParamIn[, 2:(NParam - 4)], Direction)
         ParamOut[, (NParam - 3):NParam] <- FUN_SNOW(ParamIn[, (NParam - 3):NParam], Direction)
-        ParamOut[, 1                  ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
+        ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
         if (!Bool) {
           ParamOut <- ParamOut[1, ]
         }
@@ -124,7 +129,7 @@
           ParamOut[, 2:(NParam - 2)] <- FUN_GR(ParamIn[, 2:(NParam - 2)],  Direction)
         }
         ParamOut[, (NParam - 1):NParam] <- FUN_SNOW(ParamIn[, (NParam - 1):NParam], Direction)
-        ParamOut[, 1                  ] <- FUN_ROUT(as.matrix(ParamIn[, 1]), Direction)
+        ParamOut[, 1:FeatFUN_MOD$NParamSD] <- FUN_ROUT(ParamIn[, 1:FeatFUN_MOD$NParamSD], Direction)
         if (!Bool) {
           ParamOut <- ParamOut[1, ]
         }