diff --git a/DESCRIPTION b/DESCRIPTION
index c0473074da506b3814b745ac3ee15fbc72cc7e0d..66ef4c4a7092944b9fa47954736804c609a653dc 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.3.0.18
+Version: 1.3.0.19
 Date: 2019-05-21
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
diff --git a/NEWS.rmd b/NEWS.rmd
index e5de3ce0523d7ecb539e35bab7ac0ff6d96632b2..a0930d178c60259b80fca7b9d53afdea18449b0a 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -14,7 +14,7 @@ output:
 
 
 
-### 1.3.0.18 Release Notes (2019-05-21)
+### 1.3.0.19 Release Notes (2019-05-21)
 
 
 #### New features
diff --git a/R/RunModel_CemaNeigeGR4H.R b/R/RunModel_CemaNeigeGR4H.R
index a9526bf67e27700ec7688f2326424f5e408756a7..ea9edfe5334508e9cc001d81457776ebd7b9453f 100644
--- a/R/RunModel_CemaNeigeGR4H.R
+++ b/R/RunModel_CemaNeigeGR4H.R
@@ -4,6 +4,7 @@ RunModel_CemaNeigeGR4H <- function(InputsModel,RunOptions,Param){
   ## Initialization of variables
   IsHyst <- inherits(RunOptions, "hysteresis")
   NParam <- ifelse(test = IsHyst, yes = 8L, no = 6L)
+  NParamCN <- NParam - 4L
   NStates <- 4L
   FortranOutputs <- .FortranOutputs(GR = "GR4H", isCN = TRUE)
 
@@ -70,16 +71,16 @@ RunModel_CemaNeigeGR4H <- function(InputsModel,RunOptions,Param){
                             InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
                             InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1],                       ### input series of air mean temperature [degC]
                             MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer],                       ### value of annual mean solid precip [mm/y]
-                            NParam=as.integer(NParam),                                                    ### number of model parameter = 2
+                            NParam=as.integer(NParamCN),                                                  ### number of model parameters = 2 or 4
                             Param=as.double(ParamCemaNeige),                                              ### parameter set
-                            NStates=as.integer(NStates),                                                  ### number of state variables used for model initialising = 2
+                            NStates=as.integer(NStates),                                                  ### number of state variables used for model initialisation = 4
                             StateStart=StateStartCemaNeige,                                               ### state variables used when the model run starts
                             IsHyst = as.integer(IsHyst),                                                  ### use of hysteresis
                             NOutputs=as.integer(length(IndOutputsCemaNeige)),                             ### number of output series
                             IndOutputs=IndOutputsCemaNeige,                                               ### indices of output series
                         ##outputs                                                               
                             Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
-                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                   ### state variables at the end of the model run (reservoir levels [mm] and HU)
+                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                   ### 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;
diff --git a/R/RunModel_CemaNeigeGR4J.R b/R/RunModel_CemaNeigeGR4J.R
index cdcbfcf4869db64ede60cb7f0e223da840cd8bd5..9254b597e048e3ff63da85a68503f7241398647b 100644
--- a/R/RunModel_CemaNeigeGR4J.R
+++ b/R/RunModel_CemaNeigeGR4J.R
@@ -4,6 +4,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
   ## Initialization of variables
   IsHyst <- inherits(RunOptions, "hysteresis")
   NParam <- ifelse(test = IsHyst, yes = 8L, no = 6L)
+  NParamCN <- NParam - 4L
   NStates <- 4L
   FortranOutputs <- .FortranOutputs(GR = "GR4J", isCN = TRUE)
 
@@ -70,16 +71,16 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param){
                             InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
                             InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1],                       ### input series of air mean temperature [degC]
                             MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer],                       ### value of annual mean solid precip [mm/y]
-                            NParam=as.integer(NParam),                                                    ### number of model parameter = 2
+                            NParam=as.integer(NParamCN),                                                  ### number of model parameters = 2 or 4
                             Param=as.double(ParamCemaNeige),                                              ### parameter set
-                            NStates=as.integer(NStates),                                                  ### number of state variables used for model initialising = 2
+                            NStates=as.integer(NStates),                                                  ### number of state variables used for model initialising = 4
                             StateStart=StateStartCemaNeige,                                               ### state variables used when the model run starts
                             IsHyst = as.integer(IsHyst),                                                  ### use of hysteresis
                             NOutputs=as.integer(length(IndOutputsCemaNeige)),                             ### number of output series
                             IndOutputs=IndOutputsCemaNeige,                                               ### indices of output series
                         ##outputs                                                               
                             Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
-                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                   ### state variables at the end of the model run (reservoir levels [mm] and HU)
+                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                   ### 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;
diff --git a/R/RunModel_CemaNeigeGR5J.R b/R/RunModel_CemaNeigeGR5J.R
index 55fd7d2d5763e19bce9e76fa69986714ff56653d..699117ab9c77f488acde82e3210b0ded06c74e64 100644
--- a/R/RunModel_CemaNeigeGR5J.R
+++ b/R/RunModel_CemaNeigeGR5J.R
@@ -3,6 +3,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
   
   IsHyst <- inherits(RunOptions, "hysteresis")
   NParam <- ifelse(test = IsHyst, yes = 9L, no = 7L)
+  NParamCN <- NParam - 5L
   NStates <- 4L
   FortranOutputs <- .FortranOutputs(GR = "GR5J", isCN = TRUE)
 
@@ -68,16 +69,16 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
                             InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1], ### input series of fraction of solid precipitation [0-1]
                             InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1],                       ### input series of air mean temperature [degC]
                             MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer],                       ### value of annual mean solid precip [mm/y]
-                            NParam=as.integer(NParam),                                                    ### number of model parameter = 2
+                            NParam=as.integer(NParamCN),                                                  ### number of model parameters = 2 or 4
                             Param=as.double(ParamCemaNeige),                                              ### parameter set
-                            NStates=as.integer(NStates),                                                  ### number of state variables used for model initialising = 2
+                            NStates=as.integer(NStates),                                                  ### number of state variables used for model initialising = 4
                             StateStart=StateStartCemaNeige,                                               ### state variables used when the model run starts
                             IsHyst = as.integer(IsHyst),                                                  ### use of hysteresis
                             NOutputs=as.integer(length(IndOutputsCemaNeige)),                             ### number of output series
                             IndOutputs=IndOutputsCemaNeige,                                               ### indices of output series
                         ##outputs                                                               
                             Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)), ### output series [mm]
-                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                   ### state variables at the end of the model run (reservoir levels [mm] and HU)
+                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                   ### 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;
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index 49bb37cb0cfcbefc2eff426efcce511b6973ee18..a288a008bc4e05946d2bc791bec6c1069219f657 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -3,6 +3,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
   
   IsHyst <- inherits(RunOptions, "hysteresis")
   NParam <- ifelse(test = IsHyst, yes = 10L, no = 8L)
+  NParamCN <- NParam - 6L
   NStates <- 4L
   FortranOutputs <- .FortranOutputs(GR = "GR6J", isCN = TRUE)
 
@@ -72,16 +73,16 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
                             InputsFracSolidPrecip=InputsModel$LayerFracSolidPrecip[[iLayer]][IndPeriod1],  ### input series of fraction of solid precipitation [0-1]
                             InputsTemp=InputsModel$LayerTemp[[iLayer]][IndPeriod1],                        ### input series of air mean temperature [degC]
                             MeanAnSolidPrecip=RunOptions$MeanAnSolidPrecip[iLayer],                        ### value of annual mean solid precip [mm/y]
-                            NParam=as.integer(NParam),                                                     ### number of model parameter = 2
+                            NParam=as.integer(NParamCN),                                                   ### number of model parameters = 2 or 4
                             Param=as.double(ParamCemaNeige),                                               ### parameter set
-                            NStates=as.integer(NStates),                                                   ### number of state variables used for model initialising = 2
+                            NStates=as.integer(NStates),                                                   ### number of state variables used for model initialising = 4
                             StateStart=StateStartCemaNeige,                                                ### state variables used when the model run starts
                             IsHyst = as.integer(IsHyst),                       ### use of hysteresis                        
                             NOutputs=as.integer(length(IndOutputsCemaNeige)),                              ### number of output series
                             IndOutputs=IndOutputsCemaNeige,                                                ### indices of output series
                         ##outputs                                                               
                             Outputs=matrix(as.double(-999.999),nrow=LInputSeries,ncol=length(IndOutputsCemaNeige)),  ### output series [mm]
-                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                    ### state variables at the end of the model run (reservoir levels [mm] and HU)
+                            StateEnd=rep(as.double(-999.999),as.integer(NStates))                                    ### 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;