diff --git a/DESCRIPTION b/DESCRIPTION
index 345698f56c7735c350ac620c99175c56e96888cd..1544339a4a97627ae16ba313965d55be5a5fe310 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.4.3.5
+Version: 1.4.3.6
 Date: 2019-12-12
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
diff --git a/NEWS.md b/NEWS.md
index 4dd5a64f9e6eb18da5ce3d0c14d71db2815ab946..cbb4d11274b0fee3b0d86e6078c3e49902c5595a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 
 
 
-### 1.4.3.5 Release Notes (2019-12-12)
+### 1.4.3.6 Release Notes (2019-12-12)
 
 
 #### New features
diff --git a/R/RunModel_CemaNeigeGR5H.R b/R/RunModel_CemaNeigeGR5H.R
index 53a3c7a2bb4fbd0e4039637358e03f1f6d4c4f40..0b59a5883bbe043daee24e3776c0176295d20ef2 100644
--- a/R/RunModel_CemaNeigeGR5H.R
+++ b/R/RunModel_CemaNeigeGR5H.R
@@ -43,8 +43,8 @@ RunModel_CemaNeigeGR5H <- function(InputsModel,RunOptions,Param){
   if(identical(RunOptions$IndPeriod_WarmUp,as.integer(0))){ RunOptions$IndPeriod_WarmUp <- NULL; }
   IndPeriod1   <- c(RunOptions$IndPeriod_WarmUp,RunOptions$IndPeriod_Run);
   LInputSeries <- as.integer(length(IndPeriod1))
-  if("all" %in% RunOptions$Outputs_Sim){ IndOutputs <- as.integer(1:length(FortranOutputs)); 
-  } else { IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim);  }
+  if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputs)); 
+  } else { IndOutputsMod <- which(FortranOutputs %in% RunOptions$Outputs_Sim);  }
   
   ParamCemaNeige <- Param[(length(Param)-1-2*as.integer(IsHyst)):length(Param)];
   NParamMod      <- as.integer(length(Param)-(2+2*as.integer(IsHyst)));
@@ -108,7 +108,10 @@ RunModel_CemaNeigeGR5H <- function(InputsModel,RunOptions,Param){
     CatchMeltAndPliq  <- InputsModel$Precip[IndPeriod1]; }
   
   
+  
   ##MODEL______________________________________________________________________________________##
+  if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputs$GR)); 
+  } else { IndOutputsMod <- which(FortranOutputs$GR %in% RunOptions$Outputs_Sim);  }
   
   ##Use_of_IniResLevels
   if(!is.null(RunOptions$IniResLevels)){
@@ -130,16 +133,17 @@ RunModel_CemaNeigeGR5H <- function(InputsModel,RunOptions,Param){
                       NStates=as.integer(length(RunOptions$IniStates)), ### number of state variables used for model initialising
                       StateStart=RunOptions$IniStates,                  ### state variables used when the model run starts
                       Imax=Imax,                                        ### maximal capacity of interception store
-                      NOutputs=as.integer(length(IndOutputs)),          ### number of output series
-                      IndOutputs=IndOutputs,                            ### indices of output series
+                      NOutputs=as.integer(length(IndOutputsMod)),          ### number of output series
+                      IndOutputs=IndOutputsMod,                            ### indices of output series
                       ##outputs
-                      Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputs)),  ### output series [mm or mm/h]
+                      Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsMod)),  ### output series [mm or mm/h]
                       StateEnd=rep(as.double(-999.999),length(RunOptions$IniStates))                  ### state variables at the end of the model run
   )
   RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
   RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
   if (ExportStateEnd) {
     RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
+    idNStates <- seq_len(NStates*NLayers) %% NStates
     RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel, IsHyst = IsHyst,
                                         ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
                                         IntStore = RESULTS$StateEnd[4L],
@@ -150,7 +154,7 @@ RunModel_CemaNeigeGR5H <- function(InputsModel,RunOptions,Param){
                                         GlocmaxCemaNeigeLayers = CemaNeigeStateEnd[seq_len(NStates*NLayers)[idNStates == 0]],
                                         verbose = FALSE)
   }
-  
+
   if(inherits(RunOptions,"CemaNeige") & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputs$GR[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
   
   ##Output_data_preparation