diff --git a/DESCRIPTION b/DESCRIPTION
index 783aa3553e07212941f30f8b669ed7b52d0dd044..aa8950b0408e46697b9e1b6b541a9266c652f9aa 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.2.57
+Version: 1.3.2.58
 Date: 2019-11-19
 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 640cfb3b16f9eae862a8676af53f283ffcd0b8e4..d5fe4de40df2e023d6736321e06bd7a83b7f4552 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 
 
 
-### 1.3.2.57 Release Notes (2019-11-19)
+### 1.3.2.58 Release Notes (2019-11-19)
 
 
 #### New features
@@ -14,6 +14,7 @@
 
 - Fixed bug in <code>TransfoParam_GR1A()</code>. The number of model parameters was wrong (2 instead of 1) which caused an error when the GR1 model during the calibration.
 - Fixed bug in <code>plot.OutputsModel()</code>. The function does not return any error message when <code>log_scale = TRUE</code>, <code>Qobs = NULL</code> and user want to draw flows time series.
+- Fixed bug in <code>RunModel_&#42;GR&#42;()</code>. The functions do not return any error message anymore due to slightly negative values returned by GR4H, GR4J, GR5J or GR6J Fortran codes (the message was returned by <code>CreateIniStates()</code> when the final states were created). The <code>RunModel_&#42;GR&#42;()</code> functions now returns zero instead of these slightly negative values, except for the ExpStore where negatives values are allowed.
 
 
 #### CRAN-compatibility updates
diff --git a/R/RunModel_CemaNeigeGR4H.R b/R/RunModel_CemaNeigeGR4H.R
index 7d62f57fb924afac196a03ed67c87e4f6a083b4e..7aecbaf11c5f7e0e7b701069f7f618ca5119a4a1 100644
--- a/R/RunModel_CemaNeigeGR4H.R
+++ b/R/RunModel_CemaNeigeGR4H.R
@@ -131,6 +131,7 @@ RunModel_CemaNeigeGR4H <- function(InputsModel,RunOptions,Param){
       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_CemaNeigeGR4H, InputsModel = InputsModel, IsHyst = IsHyst,
                                             ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
diff --git a/R/RunModel_CemaNeigeGR4J.R b/R/RunModel_CemaNeigeGR4J.R
index 9365769e56b1bdf637312a6d54918f4bbb526184..f3a7f8226511264f189df5346a32e29b8a8528d7 100644
--- a/R/RunModel_CemaNeigeGR4J.R
+++ b/R/RunModel_CemaNeigeGR4J.R
@@ -131,6 +131,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;
       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_CemaNeigeGR4J, InputsModel = InputsModel, IsHyst = IsHyst,
                                             ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
diff --git a/R/RunModel_CemaNeigeGR5J.R b/R/RunModel_CemaNeigeGR5J.R
index 529359df264ecb18540b8b11c36e2bddab804c02..2d41be506a20dd5cebe5ef3e0ca5006aa5ada343 100644
--- a/R/RunModel_CemaNeigeGR5J.R
+++ b/R/RunModel_CemaNeigeGR5J.R
@@ -128,7 +128,8 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
-      if (ExportStateEnd) { 
+      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_CemaNeigeGR5J, InputsModel = InputsModel, IsHyst = IsHyst,
                                             ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL,
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index 3baf57c336b4027b8e3ca60914945014479b8b35..54b4497f6f6a989952eb6978d37f9bc1e46af026 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -133,7 +133,8 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
-      if (ExportStateEnd) { 
+      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_CemaNeigeGR6J, InputsModel = InputsModel, IsHyst = IsHyst,
                                             ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = RESULTS$StateEnd[3L],
diff --git a/R/RunModel_GR4H.R b/R/RunModel_GR4H.R
index fc29cbd2d9da720765009158a22a98580e24bc29..71ecde1d427433921c2ef3c5ce42d0fca660303c 100644
--- a/R/RunModel_GR4H.R
+++ b/R/RunModel_GR4H.R
@@ -64,7 +64,8 @@ RunModel_GR4H <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
-      if (ExportStateEnd) { 
+      if (ExportStateEnd) {
+        RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
         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)],
diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index cb1fad6c523831b0c4b1d1198111300ad6e9eb74..8d383a18a17e16f0ee3aac77e952ebc96727b47f 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -63,7 +63,8 @@ RunModel_GR4J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
-      if (ExportStateEnd) { 
+      if (ExportStateEnd) {
+        RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
         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)],
diff --git a/R/RunModel_GR5J.R b/R/RunModel_GR5J.R
index 96fc0b320e1125a9fa4c378f51c958d7fecf7fe4..56e750ae73aef1affd9762460837d445f3dc1b5a 100644
--- a/R/RunModel_GR5J.R
+++ b/R/RunModel_GR5J.R
@@ -64,7 +64,8 @@ RunModel_GR5J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
-      if (ExportStateEnd) { 
+      if (ExportStateEnd) {
+        RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
         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)],
diff --git a/R/RunModel_GR6J.R b/R/RunModel_GR6J.R
index 7fff09370b0ce6dcd161ff71c11cfeaea9ca0497..15bd28423f7efa46e2962b9aa47dea9fddce0813 100644
--- a/R/RunModel_GR6J.R
+++ b/R/RunModel_GR6J.R
@@ -69,7 +69,8 @@ RunModel_GR6J <- function(InputsModel,RunOptions,Param){
                      )
       RESULTS$Outputs[ round(RESULTS$Outputs ,3)==(-999.999)] <- NA;
       RESULTS$StateEnd[round(RESULTS$StateEnd,3)==(-999.999)] <- NA;
-      if (ExportStateEnd) { 
+      if (ExportStateEnd) {
+        RESULTS$StateEnd[-3L] <- ifelse(RESULTS$StateEnd[-3L] < 0, 0, RESULTS$StateEnd[-3L]) ### remove negative values except for the ExpStore location
         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)],