diff --git a/.Rbuildignore b/.Rbuildignore
index 03d57e0e34158453fb11710577b435fedbba5b93..02bca46bffb309865243a7c0ee8e0bd876f74ed7 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,3 +1,4 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
 ^\.Rprofile$
+^packrat/
diff --git a/.gitignore b/.gitignore
index de2435597a53785e2ac420ba5c1973b355a8874f..3c605852332ec79777be9f733f93d0c252b81397 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 .Rhistory
 .RData
 airGR.Rproj
+packrat/lib*/
diff --git a/DESCRIPTION b/DESCRIPTION
index 26d9ac2579a10fa0e0cab41b0a6f449fcdf11a73..09f64021ab4cdf58f18fea7e70d65fbb6c1a7793 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.0.8.1
-Date: 2017-06-14
+Version: 1.0.8.2
+Date: 2017-06-21
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl")),
   person("Charles", "Perrin", role = c("aut", "ths")),
diff --git a/NAMESPACE b/NAMESPACE
index da72fe77336a757de68067cc7da2a45498596c20..b5a907419071c70bc39f67f647cd6749a172d5de 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -18,6 +18,7 @@ S3method("plot", "OutputsModel")
 export(Calibration)
 export(Calibration_Michel)
 export(CreateCalibOptions)
+export(CreateIniStates)
 export(CreateInputsCrit)
 export(CreateInputsModel)
 export(CreateRunOptions)
diff --git a/R/CreateIniStates.R b/R/CreateIniStates.R
new file mode 100644
index 0000000000000000000000000000000000000000..4540c537e2b65d18a873725eb192fe4a453b7d1b
--- /dev/null
+++ b/R/CreateIniStates.R
@@ -0,0 +1,227 @@
+CreateIniStates <- function(FUN_MOD, InputsModel,
+                            ProdStore = 0.3, RoutStore = 0.5, ExpStore = NULL,
+                            UH1 = NULL, UH2 = NULL,
+                            GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+                            verbose = TRUE) {
+  
+  
+  ObjectClass <- NULL
+  
+  UH1n <- 20L
+  UH2n <- UH1n * 2L
+  
+  
+  ## check FUN_MOD
+  BOOL <- FALSE
+  if (identical(FUN_MOD, RunModel_GR4H)) {
+    ObjectClass <- c(ObjectClass, "GR", "hourly")
+    BOOL <- TRUE
+  }
+  if (identical(FUN_MOD, RunModel_GR4J) |
+      identical(FUN_MOD, RunModel_GR5J) |
+      identical(FUN_MOD, RunModel_GR6J)) {
+    ObjectClass <- c(ObjectClass, "GR", "daily")
+    BOOL <- TRUE
+  }
+  if (identical(FUN_MOD, RunModel_GR2M)) {
+    ObjectClass <- c(ObjectClass, "GR", "monthly")
+    BOOL <- TRUE
+  }
+  if (identical(FUN_MOD, RunModel_GR1A)) {
+    stop("'RunModel_GR1A' does not require 'IniStates' object")
+  }
+  if (identical(FUN_MOD, RunModel_CemaNeige)) {
+    ObjectClass <- c(ObjectClass, "CemaNeige", "daily")
+    BOOL <- TRUE
+  }
+  if (identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
+      identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
+      identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
+    ObjectClass <- c(ObjectClass, "GR", "CemaNeige", "daily")
+    BOOL <- TRUE
+  }
+  if (!BOOL) {
+    stop("Incorrect 'FUN_MOD' for use in 'CreateIniStates'")
+    return(NULL)
+  }
+  
+  ## check InputsModel
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+    return(NULL)
+  }
+  if ("GR" %in% ObjectClass & !inherits(InputsModel, "GR")) {
+    stop("'InputsModel' must be of class 'GR'")
+    return(NULL)
+  }
+  if ("CemaNeige" %in% ObjectClass &
+      !inherits(InputsModel, "CemaNeige")) {
+    stop("'RunModel_CemaNeigeGR*' must be of class 'CemaNeige'")
+    return(NULL)
+  }
+  
+  
+  ## check states
+  if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
+    if (is.null(ExpStore)) {
+      stop("'RunModel_*GR6J' need an 'ExpStore' value")
+      return(NULL)
+    }
+  } else if (!is.null(ExpStore) & verbose) {
+    warning("This 'FUN_MOD' does not require 'ExpStore'. Value set to NA")
+    ExpStore <- Inf
+  }
+  
+  if (identical(FUN_MOD, RunModel_GR2M) & verbose) {
+    if (!is.null(UH1)) {
+      warning("This 'FUN_MOD' does not require 'UH1'. Values set to NA")
+      UH1 <- rep(Inf, UH1n)
+    }
+    if (!is.null(UH2)) {
+      warning("This 'FUN_MOD' does not require 'UH2'. Values set to NA")
+      UH2 <- rep(Inf, UH2n)
+    }
+    
+  }
+  
+  if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !is.null(UH1) & verbose) {
+    warning("This 'FUN_MOD' does not require 'UH1'. Values set to NA")
+    UH1 <- rep(Inf, UH1n)
+  }
+  
+  if("CemaNeige" %in% ObjectClass &
+     (is.null(GCemaNeigeLayers) | is.null(eTGCemaNeigeLayers))) {
+    stop("'RunModel_CemaNeigeGR*' need values for 'GCemaNeigeLayers' and 'GCemaNeigeLayers'")
+    return(NULL)
+  }
+  if(!"CemaNeige" %in% ObjectClass &
+     (!is.null(GCemaNeigeLayers) | !is.null(eTGCemaNeigeLayers)) & verbose) {
+    warning("This 'FUN_MOD' does not require 'GCemaNeigeLayers' and 'GCemaNeigeLayers'. Values set to NA")
+    GCemaNeigeLayers   <- Inf
+    eTGCemaNeigeLayers <- Inf
+  }
+  
+  
+  ## set states
+  if("CemaNeige" %in% ObjectClass) {
+    NLayers <- length(InputsModel$LayerPrecip)
+  } else {
+    NLayers <- 1
+  }
+  
+  
+  ## manage NULL values
+  if (is.null(ExpStore)) {
+    ExpStore <- Inf 
+  }
+  if (is.null(UH1)) {
+    if ("hourly"  %in% ObjectClass) {
+      k <- 24
+    } else {
+      k <- 1
+    }
+    UH1 <- rep(Inf, 20 * k)
+  }
+  if (is.null(UH2)) {
+    if ("hourly"  %in% ObjectClass) {
+      k <- 24
+    } else {
+      k <- 1
+    }
+    UH2 <- rep(Inf, 40 * k)
+  }
+  if (is.null(GCemaNeigeLayers)) {
+    GCemaNeigeLayers <- rep(Inf, NLayers)
+  }
+  if (is.null(eTGCemaNeigeLayers)) {
+    eTGCemaNeigeLayers <- rep(Inf, NLayers)
+  } 
+  
+  
+  ## check length
+  if (!is.numeric(ProdStore) || length(ProdStore) != 1L) {
+    stop("'ProdStore' must be numeric of length one")
+  }
+  if (!is.numeric(RoutStore) || length(RoutStore) != 1L) {
+    stop("'RoutStore' must be numeric of length one")
+  }
+  if (!is.numeric(ExpStore) || length(ExpStore) != 1L) {
+    stop("'ExpStore' must be numeric of length one")
+  }
+  if ( "hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != 480L)) {
+    stop(sprintf("'UH1' must be numeric of length 480 (%i * 24)", UH1n))
+  }
+  if (!"hourly" %in% ObjectClass & (!is.numeric(UH1) || length(UH1) != 20L)) {
+    stop(sprintf("'UH1' must be numeric of length %i", UH1n))
+  }
+  if ( "hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != 960L)) {
+    stop(sprintf("'UH2' must be numeric of length 960 (%i * 24)", UH2n))
+  }
+  if (!"hourly" %in% ObjectClass & (!is.numeric(UH2) || length(UH2) != 40L)) {
+    stop(sprintf("'UH2' must be numeric of length %i (2 * %i)", UH2n, UH1n))
+  }
+  if (!is.numeric(GCemaNeigeLayers) || length(GCemaNeigeLayers) != NLayers) {
+    stop(sprintf("'GCemaNeigeLayers' must be numeric of length %i", NLayers))
+  }
+  if (!is.numeric(eTGCemaNeigeLayers) || length(eTGCemaNeigeLayers) != NLayers) {
+    stop(sprintf("'eTGCemaNeigeLayers' must be numeric of length %i", NLayers))
+  }
+  
+  
+  # if ("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass) {
+  #   if ("hourly"  %in% ObjectClass) {
+  #     NState <- 3 * 24 * 20 + 7
+  #   }
+  #   if ("daily"   %in% ObjectClass) {
+  #     if (identical(FUN_MOD, RunModel_GR5J)) {
+  #       NState <-
+  #         2 * 20 + 2 * NLayers + 7
+  #     } else {
+  #       NState <- 3 * 20 + 2 * NLayers + 7
+  #     }
+  #   }
+  #   if ("monthly" %in% ObjectClass) {
+  #     NState <- 2
+  #   }
+  #   if ("yearly"  %in% ObjectClass) {
+  #     NState <- 1
+  #   }
+  # }
+  # if (!is.null(IniStates)) {
+  #   if (!is.vector(IniStates) | !is.numeric(IniStates)) {
+  #     stop("IniStates must be a vector of numeric values")
+  #     return(NULL)
+  #   }
+  #   if (length(IniStates) != NState) {
+  #     stop(paste0(
+  #       "The length of IniStates must be ",
+  #       NState,
+  #       " for the chosen FUN_MOD"
+  #     ))
+  #     return(NULL)
+  #   }
+  # } else {
+  #   IniStates <- as.double(rep(0.0, NState))
+  #   IniStates[1:3] <- NA
+  # }
+  
+  # if ("yearly"  %in% ObjectClass) {
+  #   IniStates <- c(ProdStore)
+  # }
+  # else if ("monthly"  %in% ObjectClass) {
+  #   IniStates <- c(ProdStore, RoutStore)
+  # }
+  
+  
+  IniStates <- list(Store = list(Prod = ProdStore, Rout = RoutStore, Exp = ExpStore),
+                    UH = list(UH1 = UH1, UH2 = UH2),
+                    CemaNeigeLayers = list(G = GCemaNeigeLayers, eTG = eTGCemaNeigeLayers))
+  IniStatesNA <- unlist(IniStates)
+  IniStatesNA[is.infinite(IniStatesNA)] <- NA
+  IniStatesNA <- relist(IniStatesNA, skeleton = IniStates)
+  
+  class(IniStatesNA) <- c("IniStates", ObjectClass)
+  return(IniStatesNA)
+  
+  
+}
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index ced3cc714e71f7853fd928799b0ea25a010c7a04..3fb6c9204deaae33fca7bc66bc9e3b107717b3ac 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -1,8 +1,8 @@
-CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod_Run,IniStates=NULL,IniResLevels=NULL,
-                             Outputs_Cal=NULL,Outputs_Sim="all",RunSnowModule=TRUE,MeanAnSolidPrecip=NULL, verbose = TRUE){
+CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, IniStates = NULL, IniResLevels = NULL, 
+                             Outputs_Cal = NULL, Outputs_Sim = "all", RunSnowModule = TRUE, MeanAnSolidPrecip = NULL,  verbose = TRUE) {
 
 
-  ObjectClass <- NULL;
+  ObjectClass <- NULL
 
   ##check_FUN_MOD
     BOOL <- FALSE;
@@ -93,37 +93,158 @@ CreateRunOptions <- function(FUN_MOD,InputsModel,IndPeriod_WarmUp=NULL,IndPeriod
     }
     if(!is.null(WTxt) & verbose){ warning(WTxt); }
 
-
-  ##check_IniStates_and_IniResLevels
-    if(is.null(IniStates) & is.null(IniResLevels) & verbose){ 
-      warning("\t Model states initialisation not defined -> default configuration used \n"); }
+    
+    ## check IniResLevels    
+    if ("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)) {
+      if (!is.null(IniResLevels)) {
+        if (!is.vector(IniResLevels) | !is.numeric(IniResLevels) | any(is.na(IniResLevels))) {
+          stop("IniResLevels must be a vector of numeric values \n")
+          return(NULL)
+        }
+        if ((identical(FUN_MOD, RunModel_GR4H) |
+             identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
+             identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
+             identical(FUN_MOD, RunModel_GR2M)) &
+            length(IniResLevels) != 2) {
+          stop("The length of IniStates must be 2 for the chosen FUN_MOD \n")
+          return(NULL)
+        }
+        if ((identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) &
+            length(IniResLevels) != 3) {
+          stop("The length of IniStates must be 3 for the chosen FUN_MOD \n")
+          return(NULL)
+        }
+      } else if (is.null(IniStates)) {
+        if (identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) {
+          IniResLevels <- as.double(c(0.3, 0.5, 0))
+        } else {
+          IniResLevels <- as.double(c(0.3, 0.5, NA))
+        }
+      }
+    } else {
+      if (!is.null(IniResLevels)) {
+        stop("IniResLevels can only be used with monthly or daily or hourly GR models \n")
+      }
+    }
+  ## check IniStates
+    if (is.null(IniStates) & is.null(IniResLevels) & verbose) {
+      warning("\t Model states initialisation not defined -> default configuration used \n")
+    }
+    if (!is.null(IniStates) & !is.null(IniResLevels) & verbose) {
+      warning("\t IniStates and IniResLevels are both defined -> Store levels are taken from IniResLevels \n")
+    }
     if("CemaNeige" %in% ObjectClass){ NLayers <- length(InputsModel$LayerPrecip); } else { NLayers <- 0; }
     NState <- NULL;
-    if("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass){ 
-      if("hourly"  %in% ObjectClass){ NState <- 3*24*20 + 7; }
-      if("daily"   %in% ObjectClass){ if (identical(FUN_MOD,RunModel_GR5J)){NState <- 2*20 + 2*NLayers + 7; } else {NState <- 3*20 + 2*NLayers + 7; }}
+    if("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass){
+      if("hourly"  %in% ObjectClass){ NState <- 7 + 3*24*20 }
+      if("daily"   %in% ObjectClass){ NState <- 7 + 3*20 + 2*NLayers }
       if("monthly" %in% ObjectClass){ NState <- 2; }
       if("yearly"  %in% ObjectClass){ NState <- 1; }
     }
-    if(!is.null(IniStates)){
-      if(!is.vector( IniStates)   ){ stop("IniStates must be a vector of numeric values \n"); return(NULL);  }
-      if(!is.numeric(IniStates)   ){ stop("IniStates must be a vector of numeric values \n"); return(NULL);  }
-      if(length(IniStates)!=NState){ stop(paste("The length of IniStates must be ",NState," for the chosen FUN_MOD \n",sep="")); return(NULL);  }
-    } else {
-      IniStates <- as.double(rep(0.0,NState));
-    }
-    if("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)){
-      if(!is.null(IniResLevels)){
-        if(!is.vector(IniResLevels) ){ stop("IniResLevels must be a vector of numeric values \n"); return(NULL);  }
-        if(!is.numeric(IniResLevels)){ stop("IniResLevels must be a vector of numeric values \n"); return(NULL);  }
-        if(length(IniResLevels)!=2 ) { stop("The length of IniStates must be 2 for the chosen FUN_MOD \n"); return(NULL);  }
-      } else {
-        IniResLevels <- as.double(c(0.3,0.5));
+    if (!is.null(IniStates)) {
+      if (!inherits(IniStates, "IniStates")) {
+        stop("IniStates must be an object of class IniStates\n")
+        return(NULL)
+      }
+      # print(str(IniStates))
+      # print(ObjectClass)
+      # print(class(IniStates))
+      # print(sum(ObjectClass %in% class(IniStates)))
+      if (sum(ObjectClass %in% class(IniStates)) < 2) {
+        stop(paste0("Non convenient IniStates for this FUN_MOD\n"))
+        return(NULL)
+      }
+      if (identical(FUN_MOD, RunModel_GR1A) & !is.null(IniStates)) { ## GR1A
+        stop(paste0("IniStates is not available for this FUN_MOD\n"))
+        return(NULL)
+      }
+      if ((identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J)) & !all(is.na(IniStates$UH$UH1))) { ## GR5J
+        stop(paste0("Non convenient IniStates for this FUN_MOD. In IniStates, UH1 have to be a vector of NA for GR5J \n"))
+        return(NULL)
       }
+      if ((identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & is.na(IniStates$Store$Exp)) { ## GR6J
+        stop(paste0("Non convenient IniStates for this FUN_MOD. GR6J need an exponential store value in IniStates\n"))
+        return(NULL)
+      }
+      if (!(identical(FUN_MOD, RunModel_GR6J) | identical(FUN_MOD, RunModel_CemaNeigeGR6J)) & !is.na(IniStates$Store$Exp)) { ## except GR6J
+        stop(paste0("Non convenient IniStates for this FUN_MOD\n"))
+        return(NULL)
+      }
+      # if (length(na.omit(unlist(IniStates))) != NState) {
+      #   stop(paste0("The length of IniStates must be ", NState, " for the chosen FUN_MOD \n"))
+      #   return(NULL)
+      # }
+      if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$G  ))) {
+        IniStates$CemaNeigeLayers$G   <- NULL
+      }
+      if (!"CemaNeige" %in% ObjectClass & any(is.na(IniStates$CemaNeigeLayers$eTG))) {
+        IniStates$CemaNeigeLayers$eTG <- NULL
+      }
+      IniStates$Store$Rest <- rep(NA, 4)
+      IniStates <- unlist(IniStates)
+      IniStates[is.na(IniStates)] <- 0
+      # print(IniStates)
+      if ("monthly" %in% ObjectClass) {
+        IniStates <- IniStates[seq_len(NState)]
+        # print(NState)
+        # print(IniStates)
+        }
     } else {
-      if(!is.null(IniResLevels)){ stop("IniResLevels can only be used with monthly or daily or hourly GR models \n") }
+      IniStates <- as.double(rep(0.0, NState))
     }
 
+    # if (is.null(IniStates) & is.null(IniResLevels) & verbose) { 
+    #   warning("\t Model states initialisation not defined -> default configuration used \n")
+    # }
+    # if("CemaNeige" %in% ObjectClass){ NLayers <- length(InputsModel$LayerPrecip); } else { NLayers <- 0; }
+    # NState <- NULL;
+    # if("GR" %in% ObjectClass | "CemaNeige" %in% ObjectClass){ 
+    #   if("hourly"  %in% ObjectClass){ NState <- 3*24*20 + 7; }
+    #   if("daily"   %in% ObjectClass){ if (identical(FUN_MOD,RunModel_GR5J)){NState <- 2*20 + 2*NLayers + 7; } else {NState <- 3*20 + 2*NLayers + 7; }}
+    #   if("monthly" %in% ObjectClass){ NState <- 2; }
+    #   if("yearly"  %in% ObjectClass){ NState <- 1; }
+    # }
+    # if (!is.null(IniStates)) {
+    #   if (!is.vector( IniStates) | !is.numeric(IniStates)) {
+    #     stop("IniStates must be a vector of numeric values \n")
+    #     return(NULL)
+    #   }
+    #   if (length(IniStates) != NState) {
+    #     stop(paste0("The length of IniStates must be ", NState, " for the chosen FUN_MOD \n"))
+    #     return(NULL)
+    #   }
+    # } else {
+    #   IniStates <- as.double(rep(0.0, NState))
+    #   IniStates[1:3] <- NA
+    # }
+    # if ("GR" %in% ObjectClass & ("monthly" %in% ObjectClass | "daily" %in% ObjectClass | "hourly" %in% ObjectClass)) {
+    #   if (!is.null(IniResLevels)) {
+    #     if (!is.vector(IniResLevels) | !is.numeric(IniResLevels)) {
+    #       stop("IniResLevels must be a vector of numeric values \n")
+    #       return(NULL)
+    #     }
+    #     if ((identical(FUN_MOD, RunModel_GR4H) |
+    #          identical(FUN_MOD, RunModel_GR4J) | identical(FUN_MOD, RunModel_CemaNeigeGR4J) |
+    #          identical(FUN_MOD, RunModel_GR5J) | identical(FUN_MOD, RunModel_CemaNeigeGR5J) |
+    #          identical(FUN_MOD, RunModel_GR2M)) &
+    #         length(IniResLevels) != 2) {
+    #       stop("The length of IniStates must be 2 for the chosen FUN_MOD \n")
+    #       return(NULL)
+    #     }
+    #     if ((identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) &
+    #         length(IniResLevels) != 3) {
+    #       stop("The length of IniStates must be 3 for the chosen FUN_MOD \n")
+    #       return(NULL)
+    #     }
+    #   } else {
+    #     IniResLevels <- as.double(c(0.3, 0.5, 0))
+    #   }
+    # } else {
+    #   if (!is.null(IniResLevels)) {
+    #     stop("IniResLevels can only be used with monthly or daily or hourly GR models \n")
+    #   }
+    # }
+    
 
   ##check_Outputs_Cal_and_Sim
 
diff --git a/R/RunModel_CemaNeigeGR4J.R b/R/RunModel_CemaNeigeGR4J.R
index 1413f96acdd72b707118bdec654fffdc9e5586d4..01f1ecfb4934fe10caa84680f6ed9744ce6884f7 100644
--- a/R/RunModel_CemaNeigeGR4J.R
+++ b/R/RunModel_CemaNeigeGR4J.R
@@ -54,7 +54,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
 
     ##Call_DLL_CemaNeige_________________________
       for(iLayer in 1:NLayers){
-        StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
+        StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
         RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
                         ##inputs
                             LInputs=LInputSeries,                                                          ### length of input and output series
@@ -74,6 +74,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
                          )
         RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
         RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+        
 
         ##Data_storage
         CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
@@ -97,7 +98,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
       } else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim);  }
 
     ##Use_of_IniResLevels
-      if("IniResLevels" %in% names(RunOptions)){
+      if(!is.null(RunOptions$IniResLevels)){
         RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*ParamMod[1];  ### production store level (mm)
         RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*ParamMod[3];  ### routing store level (mm)
       }
@@ -120,6 +121,13 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+      RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel,
+                                          ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
+                                          UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
+                                          GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %%2 == 1]],
+                                          eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %%2 == 0]],
+                                          verbose = FALSE)
+      
       if(RunOptions$RunSnowModule & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
 
     ##Output_data_preparation
@@ -138,14 +146,14 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
       if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
-                           list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
+                           list(RESULTS$StateEnd) );
         names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
       ##DatesR_and_OutputsModel_and_SateEnd
       if(ExportDatesR==TRUE & ExportStateEnd==TRUE){
         OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                            lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
-                           list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
+                           list(RESULTS$StateEnd) );
         names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
 
     ##End
diff --git a/R/RunModel_CemaNeigeGR5J.R b/R/RunModel_CemaNeigeGR5J.R
index a703657b3f698970ecc6c4ef2731cda6e9fec782..9e9816d15dd4c4ca30f87f3904a9b36d6f86d62a 100644
--- a/R/RunModel_CemaNeigeGR5J.R
+++ b/R/RunModel_CemaNeigeGR5J.R
@@ -54,7 +54,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
 
     ##Call_DLL_CemaNeige_________________________
       for(iLayer in 1:NLayers){
-        StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
+        StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
         RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
                         ##inputs
                             LInputs=LInputSeries,                                                          ### length of input and output series
@@ -97,7 +97,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
       } else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim);  }
 
     ##Use_of_IniResLevels
-      if("IniResLevels" %in% names(RunOptions)){
+      if(!is.null(RunOptions$IniResLevels)){
         RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*ParamMod[1];  ### production store level (mm)
         RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*ParamMod[3];  ### routing store level (mm)
       }
@@ -120,6 +120,13 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+      RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR5J, InputsModel = InputsModel,
+                                          ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
+                                          UH1 = NULL, UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
+                                          GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %%2 == 1]],
+                                          eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %%2 == 0]],
+                                          verbose = FALSE)
+      
       if(RunOptions$RunSnowModule & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
 
     ##Output_data_preparation
@@ -138,14 +145,14 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
       if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
-                           list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
+                           list(RESULTS$StateEnd) );
         names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
       ##DatesR_and_OutputsModel_and_SateEnd
       if(ExportDatesR==TRUE & ExportStateEnd==TRUE){
         OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                            lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
-                           list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
+                           list(RESULTS$StateEnd) );
         names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
 
     ##End
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index 4cf4f11724b35d662b2e5ffddb30efc3f388c4b1..50d87cf7d342b07f2242a4bcffbbd9264031d9ba 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -59,7 +59,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
 
     ##Call_DLL_CemaNeige_________________________
       for(iLayer in 1:NLayers){
-        StateStartCemaNeige <- RunOptions$IniStates[ (NStatesMod+2*(iLayer-1)+1):(NStatesMod+2*(iLayer-1)+2) ];
+        StateStartCemaNeige <- RunOptions$IniStates[(7+20+40) + c(iLayer, iLayer+NLayers)]
         RESULTS <- .Fortran("frun_CemaNeige",PACKAGE="airGR",
                         ##inputs
                             LInputs=LInputSeries,                                                          ### length of input and output series
@@ -102,9 +102,10 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
       } else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim);  }
 
     ##Use_of_IniResLevels
-      if("IniResLevels" %in% names(RunOptions)){
-        RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*ParamMod[1];  ### production store level (mm)
-        RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*ParamMod[3];  ### routing store level (mm)
+      if(!is.null(RunOptions$IniResLevels)){
+        RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * ParamMod[1] ### production store level (mm)
+        RunOptions$IniStates[2] <- RunOptions$IniResLevels[2] * ParamMod[3] ### routing store level (mm)
+        RunOptions$IniStates[3] <- RunOptions$IniResLevels[3]               ### exponential store level (mm)
       }
 
     ##Call_fortan
@@ -125,6 +126,13 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+      RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_CemaNeigeGR6J, InputsModel = InputsModel,
+                                          ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
+                                          UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
+                                          GCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %%2 == 1]],
+                                          eTGCemaNeigeLayers = CemaNeigeStateEnd[seq_len(2*NLayers)[seq_len(2*NLayers) %%2 == 0]],
+                                          verbose = FALSE)
+      
       if(RunOptions$RunSnowModule & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
 
     ##Output_data_preparation
@@ -143,14 +151,14 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
       if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
-                           list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
+                           list(RESULTS$StateEnd) );
         names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
       ##DatesR_and_OutputsModel_and_SateEnd
       if(ExportDatesR==TRUE & ExportStateEnd==TRUE){
         OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                            lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
-                           list(c(RESULTS$StateEnd,CemaNeigeStateEnd)) );
+                           list(RESULTS$StateEnd) );
         names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
 
     ##End
diff --git a/R/RunModel_GR2M.R b/R/RunModel_GR2M.R
index d0b1784f981e4a1486b68beb5bad0b7982a63eb5..d6e029a5e4d19002219919997507167ed770e86c 100644
--- a/R/RunModel_GR2M.R
+++ b/R/RunModel_GR2M.R
@@ -27,8 +27,9 @@ RunModel_GR2M <- function(InputsModel,RunOptions,Param){
       } else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim);  }
 
     ##Use_of_IniResLevels
-      if("IniResLevels" %in% names(RunOptions)){
+      if(!is.null(RunOptions$IniResLevels)){
         RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1];  ### production store level (mm)
+		RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[2];  ### routing store level (mm)
       }
 
     ##Call_fortan
@@ -49,6 +50,12 @@ RunModel_GR2M <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs [round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+      RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel,
+                                          ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
+                                          UH1 = NULL, UH2 = NULL,
+                                          GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+                                          verbose = FALSE)
+      
 
     ##Output_data_preparation
       IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
diff --git a/R/RunModel_GR4H.R b/R/RunModel_GR4H.R
index c914468ff9cb1b0445523c6ec0ee0be68ec4ea8b..0333c972b5ea5c16f03e6768cd10dd789b6ea496 100644
--- a/R/RunModel_GR4H.R
+++ b/R/RunModel_GR4H.R
@@ -36,7 +36,7 @@ RunModel_GR4H <- function(InputsModel,RunOptions,Param){
       } else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim);  }
 
     ##Use_of_IniResLevels
-      if("IniResLevels" %in% names(RunOptions)){
+      if(!is.null(RunOptions$IniResLevels)){
         RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1];  ### production store level (mm)
         RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3];  ### routing store level (mm)
       }
@@ -59,6 +59,11 @@ RunModel_GR4H <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+      RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4H, InputsModel = InputsModel,
+                                          ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
+                                          UH1 = RESULTS$StateEnd[(1:(20*24))+7], UH2 = RESULTS$StateEnd[(1:(40*24))+(7+20*24)],
+                                          GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+                                          verbose = FALSE)
 
     ##Output_data_preparation
       IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index 4577fa851b15cf87399307370c224423d8d45c0e..2e522ffe69ffda0646aa2ed8dfaf1b02143821ac 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -37,7 +37,7 @@ RunModel_GR4J <- function(InputsModel,RunOptions,Param){
       } else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim);  }
 
     ##Use_of_IniResLevels
-      if("IniResLevels" %in% names(RunOptions)){
+      if(!is.null(RunOptions$IniResLevels)){
         RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1];  ### production store level (mm)
         RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3];  ### routing store level (mm)
       }
@@ -60,6 +60,11 @@ RunModel_GR4J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+      RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
+                                          ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
+                                          UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
+                                          GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+                                          verbose = FALSE)
 
     ##Output_data_preparation
       IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
@@ -74,12 +79,12 @@ RunModel_GR4J <- function(InputsModel,RunOptions,Param){
         OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                            lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]) );
         names(OutputsModel) <- c("DatesR",FortranOutputs[IndOutputs]);      }
-      ##OutputsModel_and_SateEnd_only
+      ##OutputsModel_and_StateEnd_only
       if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(RESULTS$StateEnd) );
         names(OutputsModel) <- c(FortranOutputs[IndOutputs],"StateEnd");      }
-      ##DatesR_and_OutputsModel_and_SateEnd
+      ##DatesR_and_OutputsModel_and_StateEnd
       if((ExportDatesR==TRUE & ExportStateEnd==TRUE) | "all" %in% RunOptions$Outputs_Sim){
         OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                            lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
diff --git a/R/RunModel_GR5J.R b/R/RunModel_GR5J.R
index 6cbf01ff5c36593b70620cc0f3f4d936058f7bff..489f3410570ee749be21f06bcc38bb334256a879 100644
--- a/R/RunModel_GR5J.R
+++ b/R/RunModel_GR5J.R
@@ -37,7 +37,7 @@ RunModel_GR5J <- function(InputsModel,RunOptions,Param){
       } else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim);  }
 
     ##Use_of_IniResLevels
-      if("IniResLevels" %in% names(RunOptions)){
+      if(!is.null(RunOptions$IniResLevels)){
         RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1];  ### production store level (mm)
         RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3];  ### routing store level (mm)
       }
@@ -60,6 +60,11 @@ RunModel_GR5J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+      RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5J, InputsModel = InputsModel,
+                                          ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
+                                          UH1 = NULL, UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
+                                          GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+                                          verbose = FALSE)
 
     ##Output_data_preparation
       IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
diff --git a/R/RunModel_GR6J.R b/R/RunModel_GR6J.R
index b7582eb8a231dd586a08fedd9400c9f58ed95c69..c2ca3b384b1a8ff635f56867369cb551f102ad45 100644
--- a/R/RunModel_GR6J.R
+++ b/R/RunModel_GR6J.R
@@ -41,9 +41,10 @@ RunModel_GR6J <- function(InputsModel,RunOptions,Param){
       } else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim);  }
 
     ##Use_of_IniResLevels
-      if("IniResLevels" %in% names(RunOptions)){
+      if(!is.null(RunOptions$IniResLevels)){
         RunOptions$IniStates[1] <- RunOptions$IniResLevels[1]*Param[1];  ### production store level (mm)
         RunOptions$IniStates[2] <- RunOptions$IniResLevels[2]*Param[3];  ### routing store level (mm)
+        RunOptions$IniStates[3] <- RunOptions$IniResLevels[3]            ### exponential store level (mm)
       }
 
     ##Call_fortan
@@ -64,6 +65,11 @@ RunModel_GR6J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
+      RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR6J, InputsModel = InputsModel,
+                                          ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
+                                          UH1 = RESULTS$StateEnd[(1:20)+7], UH2 = RESULTS$StateEnd[(1:40)+(7+20)],
+                                          GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+                                          verbose = FALSE)
 
     ##Output_data_preparation
       IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries;
diff --git a/man/CreateIniStates.Rd b/man/CreateIniStates.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..efb91556b576236397e01bcf90f1530a60771723
--- /dev/null
+++ b/man/CreateIniStates.Rd
@@ -0,0 +1,101 @@
+\encoding{UTF-8}
+
+
+\name{CreateIniStates}
+\alias{CreateIniStates}
+
+
+\title{Creation of the IniStates object possible required to the CreateRunOptions functions}
+
+
+\usage{
+CreateIniStates(FUN_MOD, InputsModel,
+  ProdStore = 0.3, RoutStore = 0.5, ExpStore = NULL,
+  UH1 = NULL, UH2 = NULL,
+  GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL,
+  verbose = TRUE)
+}
+
+
+\arguments{
+\item{FUN_MOD}{[function] hydrological model function (e.g. RunModel_GR4J, RunModel_CemaNeigeGR4J)}
+
+\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link{CreateInputsModel}} for details}
+
+\item{ProdStore}{[numeric]}
+
+\item{RoutStore}{[numeric]}
+
+\item{ExpStore}{[numeric]}
+
+\item{UH1}{[numeric]}
+
+\item{UH2}{[numeric]}
+
+\item{GCemaNeigeLayers}{[numeric]}
+
+\item{eTGCemaNeigeLayers}{[numeric]}
+
+\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, \code{default = TRUE}}
+
+}
+
+
+\value{
+[list] object of class \emph{IniStates} containing the initial model internal states; it includes the following:
+         \itemize{
+         \item{Store          }{}
+         \item{UH             }{}
+         \item{CemaNeigeLayers}{}
+         }
+}
+
+
+\description{
+Creation of the IniStates object required to XXXX
+}
+
+\details{
+XXXX
+}
+
+
+\examples{
+library(airGR)
+
+## loading catchment data
+data(L0123001)
+
+## preparation of the InputsModel object
+InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR, 
+                                 Precip = BasinObs$P, PotEvap = BasinObs$E)
+
+## run period selection
+Ind_Run <- seq(which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="01/01/1990 00:00"), 
+               which(format(BasinObs$DatesR, format = "\%d/\%m/\%Y \%H:\%M")=="31/12/1999 00:00"))
+
+
+## preparation of the IniStates object
+IniStates <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
+                             ProdStore = 0.1, RoutStore = 0.6, ExpStore = NULL,
+                             UH1 = c(0.52, 0.54, 0.15, rep(0, 17)),
+                             UH2 = c(0.057, 0.042, 0.015, 0.005, rep(0, 36)),
+                             GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL)
+str(IniStates)  
+
+## preparation of the RunOptions object
+RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, 
+                               IndPeriod_Run = Ind_Run, IniStates = IniStates)
+                           
+}
+
+
+\author{
+Olivier Delaigue
+}
+
+
+\seealso{
+\code{\link{CreateRunOptions}}
+}
+
diff --git a/man/CreateRunOptions.Rd b/man/CreateRunOptions.Rd
index 1607956b391e127d4b8800a46714d8520c492a88..0d5346d27d0b20a9e9a2f5181d30f528373607b9 100644
--- a/man/CreateRunOptions.Rd
+++ b/man/CreateRunOptions.Rd
@@ -19,18 +19,18 @@ CreateRunOptions(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run,
 
 \item{IniStates}{(optional) [numeric] vector of initial model internal states [mm]}
 
-\item{IniResLevels}{(optional) [numeric] vector of initial filling rates for production and routing stores (2 values between 0 and 1) [-]}
+\item{IniResLevels}{(optional) [numeric] vector of initial fillings for the GR stores (2 or 3 values according to the model) [- and/or mm]; see details}
 
 \item{Outputs_Cal}{(optional) [character] vector giving the outputs needed for the calibration \cr (e.g. c("Qsim")), the fewer outputs
  the faster the calibration}
 
-\item{Outputs_Sim}{(optional) [character] vector giving the requested outputs \cr (e.g. c("DatesR","Qsim","SnowPack")), default="all"}
+\item{Outputs_Sim}{(optional) [character] vector giving the requested outputs \cr (e.g. c("DatesR", "Qsim", "SnowPack")), default = "all"}
 
-\item{RunSnowModule}{(optional) [boolean] option indicating whether CemaNeige should be activated, default=TRUE}
+\item{RunSnowModule}{(optional) [boolean] option indicating whether CemaNeige should be activated, \code{default = TRUE}}
 
 \item{MeanAnSolidPrecip}{(optional) [numeric] vector giving the annual mean of average solid precipitation for each layer (computed from InputsModel if not defined) [mm/y]}
 
-\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, default=TRUE}
+\item{verbose}{(optional) [boolean] boolean indicating if the function is run in verbose mode or not, \code{default = TRUE}}
 }
 \value{
 [list] object of class \emph{RunOptions} containing the data required to evaluate the model outputs; it can include the following:
@@ -49,36 +49,42 @@ CreateRunOptions(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run,
 Creation of the RunOptions object required to the RunModel functions.
 }
 \details{
-Users wanting to use FUN_MOD functions that are not included in 
-the package must create their own RunOptions object accordingly.
+Users wanting to use \code{FUN_MOD} functions that are not included in 
+the package must create their own \code{RunOptions} object accordingly.
 
 ##### Initialisation options #####
 
 The model initialisation options can either be set to a default configuration or be defined by the user.
 
-This is done via three vectors: \cr \emph{IndPeriod_WarmUp}, \emph{IniStates}, \emph{IniResLevels}. \cr
+This is done via three vectors: \cr \code{IndPeriod_WarmUp}, \code{IniStates}, \code{IniResLevels}. \cr
 A default configuration is used for initialisation if these vectors are not defined.
 
 (1) Default initialisation options:
 
 \itemize{
-\item \emph{IndPeriod_WarmUp} default setting ensures a one-year warm-up using the time steps preceding the \emph{IndPeriod_Run}. 
+\item \code{IndPeriod_WarmUp} default setting ensures a one-year warm-up using the time steps preceding the \code{IndPeriod_Run}. 
 The actual length of this warm-up might be shorter depending on data availability (no missing value of climate inputs being allowed in model input series).
 
-\item \emph{IniStates} and \emph{IniResLevels} are automatically set to initialise all the model states at 0, except for the production and routing stores which are respectively initialised at 30 \% and 50 \% of their capacity. This initialisation is made at the very beginning of the model call (i.e. at the beginning of \emph{IndPeriod_WarmUp} or at the beginning of IndPeriod_Run if the warm-up period is disabled).
+\item \code{IniStates} and \code{IniResLevels} are automatically set to initialise all the model states at 0, except for the production and routing stores which are respectively initialised at 30 \% and 50 \% of their capacity. In case GR6J is used, the exponential store is initialised by default with 0 mm. This initialisation is made at the very beginning of the model call (i.e. at the beginning of \code{IndPeriod_WarmUp} or at the beginning of \code{IndPeriod_Run} if the warm-up period is disabled).
 }
 
 (2) Customisation of initialisation options:
 
 \itemize{
-\item \emph{IndPeriod_WarmUp} can be used to specify the indices of the warm-up period (within the time-series prepared in InputsModel). \cr
-- remark 1:	for most common cases, indices corresponding to one or several years preceding \emph{IndPeriod_Run} are used (e.g. \emph{IndPeriod_WarmUp = 1000:1365} and \emph{IndPeriod_Run = 1366:5000)}. \cr
-However, it is also possible to perform a long-term initialisation if other indices than the warm-up ones are set in \emph{IndPeriod_WarmUp} (e.g. \emph{IndPeriod_WarmUp <- c( 1:5000 , 1:5000 , 1:5000 ,1000:1365 )}). \cr
-- remark 2:	it is also possible to completely disable the warm-up period when using \emph{IndPeriod_WarmUp = 0L}.
-
-\item \emph{IniStates} and \emph{IniResLevels} can be used to specify the initial model states. \cr
-- remark 1:	if \emph{IniStates} is used, all model states must be provided (e.g. 60 floats [mm] are required for GR4J, GR5J and GR6J; 60+2*NLayers floats [mm] are required for CemaNeigeGR4J, CemaNeigeGR5J and CemaNeigeGR6J; see fortran source code for details). \cr
-- remark 2:	in addition to \emph{IniStates}, \emph{IniResLevels} allows to set the filling rate of the production and routing stores for the GR models. For instance for GR4J, GR5J and GR6J: \emph{IniResLevels <- c(0.3,0.5)} should be used to obtain initial fillings of 30\% and 50\% for the production and routing stores, respectively.  \emph{IniResLevels} is optional and can only be used if \emph{IniStates} is also defined (the state values corresponding to these two stores in \emph{IniStates} are not used in such case). \cr \cr
+\item \code{IndPeriod_WarmUp} can be used to specify the indices of the warm-up period (within the time-series prepared in InputsModel). \cr
+- remark 1:	for most common cases, indices corresponding to one or several years preceding \code{IndPeriod_Run} are used (e.g. \code{IndPeriod_WarmUp = 1000:1365} and \code{IndPeriod_Run = 1366:5000)}. \cr
+However, it is also possible to perform a long-term initialisation if other indices than the warm-up ones are set in \code{IndPeriod_WarmUp} (e.g. \code{IndPeriod_WarmUp <- c(1:5000 , 1:5000 , 1:5000 , 1000:1365)}). \cr
+- remark 2:	it is also possible to completely disable the warm-up period when using \code{IndPeriod_WarmUp = 0L}.
+
+\item \code{IniStates} and \code{IniResLevels} can be used to specify the initial model states. \cr
+- remark 1: \code{IniStates} and \code{IniResLevels} can not be used with GR1A.
+- remark 2:	\code{IniStates} can be set to the \code{$StateEnd} output of a previous \code{RunModel} call. \cr
+- remark 3:	if \code{IniStates} is used, all model states must be provided:
+  \itemize{
+  \item 7 values: production store level [mm],  routing store level [mm], exponential store level [mm] (only with GR6J, 0 otherwise), and the 4 other values have to be set to 0;
+  \item 60 floats [mm] are required for GR4J, GR5J and GR6J; 60+2*NLayers floats [mm] are required for \code{CemaNeigeGR4J}, \code{CemaNeigeGR5J} and \emph{CemaNeigeGR6J}; see fortran source code for details).
+  }
+- remark 4:	in addition to \code{IniStates}, \code{IniResLevels} allows to set the filling rate of the production and routing stores for the GR models. For instance for GR4J, GR5J and GR6J: \code{IniResLevels <- c(0.3, 0.5)} should be used to obtain initial fillings of 30\% and 50\% for the production and routing stores, respectively. \code{IniResLevels} is optional and can only be used if \code{IniStates} is also defined (the state values corresponding to these two stores in \code{IniStates} are not used in such case). \cr \cr
 }
 }
 \examples{
@@ -113,7 +119,7 @@ InputsCrit  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsMo
 OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
 }
 \author{
-Laurent Coron (June 2014)
+Laurent Coron, Olivier Delaigue, Guillaume Thirel
 }
 \seealso{
 \code{\link{RunModel}}, \code{\link{CreateInputsModel}}, \code{\link{CreateInputsCrit}}, \code{\link{CreateCalibOptions}}