diff --git a/DESCRIPTION b/DESCRIPTION
index 158bc726f5cc006b94855035f20daf65b49d3468..0974026234fe67a478be4ff6727a93d987eb7fd7 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.6.3.57
+Version: 1.6.3.58
 Date: 2020-11-16
 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 8bf504ebb398e93132914b261bfca39ec29b2230..1307e01ada64410a17953a16f119b6a54b20f615 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,7 +4,7 @@
 
 
 
-### 1.6.3.57 Release Notes (2020-11-16)
+### 1.6.3.58 Release Notes (2020-11-16)
 
 #### New features
 
diff --git a/R/RunModel_GR1A.R b/R/RunModel_GR1A.R
index 8e973fcc59309636f2025dd59c932dc288c45c62..2470a7fd9b96a1f1e8a462364e087264eda5f568 100644
--- a/R/RunModel_GR1A.R
+++ b/R/RunModel_GR1A.R
@@ -66,8 +66,8 @@ RunModel_GR1A <- function(InputsModel, RunOptions, Param) {
                       Outputs = matrix(as.double(-999.999), nrow = LInputSeries, ncol=length(IndOutputs)), ### output series [mm]
                       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
+  RESULTS$Outputs[ round(RESULTS$Outputs , 3) == -999.999] <- NA
+  RESULTS$StateEnd[round(RESULTS$StateEnd, 3) == -999.999] <- NA
   
   
   ## Output data preparation
diff --git a/R/RunModel_GR2M.R b/R/RunModel_GR2M.R
index 8586fd890cb53c2365e6779f9a5e2aa31fb8790e..ab6ab67be91245c540e2412945a00e3c3b0543cf 100644
--- a/R/RunModel_GR2M.R
+++ b/R/RunModel_GR2M.R
@@ -79,8 +79,8 @@ RunModel_GR2M <- function(InputsModel, RunOptions, Param) {
                       Outputs = matrix(as.double(-999.999), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm]
                       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
+  RESULTS$Outputs [round(RESULTS$Outputs , 3) == -999.999] <- NA
+  RESULTS$StateEnd[round(RESULTS$StateEnd, 3) == -999.999] <- NA
   if (ExportStateEnd) { 
     RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel, 
                                         ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL, 
diff --git a/R/RunModel_GR4H.R b/R/RunModel_GR4H.R
index 57870e06d0de1c904f6e16dc187b1d17e9caa325..38d0cfd559fa7c044fb565b3f3a1c47b85b7f9a9 100644
--- a/R/RunModel_GR4H.R
+++ b/R/RunModel_GR4H.R
@@ -84,8 +84,8 @@ RunModel_GR4H <- function(InputsModel, RunOptions, Param) {
                       Outputs = matrix(as.double(-999.999), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm]
                       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
+  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
     RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4H, InputsModel = InputsModel, 
diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index 464ab9a5acc561920d49f75fa9bfe565a5331390..96e8cc715cf3b2c2112c7169dbf6c4a01daf40c9 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -83,8 +83,8 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
                       Outputs = matrix(as.double(-999.999), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm]
                       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
+  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
     RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel, 
diff --git a/R/RunModel_GR5H.R b/R/RunModel_GR5H.R
index f683c3f61451a5a279f8cddbefc895ece06931b4..cf2cce4ce720ccceb73a239cd636981d6c964a7e 100644
--- a/R/RunModel_GR5H.R
+++ b/R/RunModel_GR5H.R
@@ -93,8 +93,8 @@ RunModel_GR5H <- function(InputsModel, RunOptions, Param) {
                       Outputs = matrix(as.double(-999.999), nrow = LInputSeries, ncol = length(IndOutputs)), ### 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
+  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
     RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5H, InputsModel = InputsModel, 
diff --git a/R/RunModel_GR5J.R b/R/RunModel_GR5J.R
index 06a62030369cccb419cdc4aa96b5ec6b20286a74..edd304c2305243d6836cdb2a6f88c410df58b2ef 100644
--- a/R/RunModel_GR5J.R
+++ b/R/RunModel_GR5J.R
@@ -84,8 +84,8 @@ RunModel_GR5J <- function(InputsModel, RunOptions, Param) {
                       Outputs = matrix(as.double(-999.999), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm]
                       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
+  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
     RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR5J, InputsModel = InputsModel, 
diff --git a/R/RunModel_GR6J.R b/R/RunModel_GR6J.R
index fb79a3979b25e07a1e12f297fd0e883044e47cab..376604a62644659cc99a9e81d94d0802a388fd85 100644
--- a/R/RunModel_GR6J.R
+++ b/R/RunModel_GR6J.R
@@ -89,8 +89,8 @@ RunModel_GR6J <- function(InputsModel, RunOptions, Param) {
                       Outputs = matrix(as.double(-999.999), nrow = LInputSeries, ncol = length(IndOutputs)), ### output series [mm]
                       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
+  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
     RESULTS$StateEnd <- CreateIniStates(FUN_MOD = RunModel_GR6J, InputsModel = InputsModel,