diff --git a/DESCRIPTION b/DESCRIPTION
index ea9b15ef120b2bd281d9fac9058535ce4b4dcaaa..ebfc7ab9f7bb303a1c096188e7475b4431282b1d 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.2.6.0
-Date: 2019-02-25
+Version: 1.2.7.0
+Date: 2019-02-28
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
   person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
diff --git a/NAMESPACE b/NAMESPACE
index 031f782463a365c03c27be6f04c2bdd97d228b39..b33fec1818ee3b60d2a4fc7d1ca561c562c5bfcd 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -52,6 +52,7 @@ export(TransfoParam_GR5J)
 export(TransfoParam_GR6J)
 export(plot.OutputsModel)
 export(plot_OutputsModel)
+exportPattern(".FortranOutputs")
 
 
 
diff --git a/NEWS.rmd b/NEWS.rmd
index e8c94884a0e04f5d8c6cbff166c7fd556c9dd749..6a0db407e475da988edb2ad76ec6f07955094861 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -13,7 +13,7 @@ output:
 
 
 
-### 1.2.6.0 Release Notes (2019-02-25) 
+### 1.2.7.0 Release Notes (2019-02-28) 
 
 
 
@@ -51,6 +51,8 @@ output:
 #### Minor user-visible changes
 
 - <code>ErrorCrit_&#42;()</code> functions now return objects of class <code>ErrorCrit</code> and <code>NSE</code>, <code>KGE</code>, <code>KGE2</code> or <code>RMSE</code>.
+- <code>.FortranOutputs()</code> private function added to manage Fortran outputs
+- Outputs of frun_GR2M Fortran subroutine were reordered
 
 ____________________________________________________________________________________
 
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index b01041181219206e56feccffb86dd79d8e7a026f..e11071681c0fca51cfff9df650e5b1170fcd2fe3 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -269,28 +269,25 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
   ##Outputs_all
   Outputs_all <- NULL
   if (identical(FUN_MOD,RunModel_GR4H)) {
-    Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Prod", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch", "AExch", "QR", "QD", "Qsim")
+    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4H")$GR)
   }
   if (identical(FUN_MOD,RunModel_GR4J) | identical(FUN_MOD,RunModel_CemaNeigeGR4J)) {
-    Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
-                     "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim")
+    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR4J")$GR)
   }
   if (identical(FUN_MOD,RunModel_GR5J) | identical(FUN_MOD,RunModel_CemaNeigeGR5J)) {
-    Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
-                     "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim")
+    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR5J")$GR)
   }
   if (identical(FUN_MOD,RunModel_GR6J) | identical(FUN_MOD,RunModel_CemaNeigeGR6J)) {
-    Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
-                     "AExch1", "AExch2", "AExch", "QR", "QRExp", "Exp", "QD", "Qsim")
+    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR6J")$GR)
   }
   if (identical(FUN_MOD,RunModel_GR2M)) {
-    Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "AE", "Pn", "Perc", "PR", "Exch", "Prod", "Rout", "Qsim")
+    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR2M")$GR)
   }
   if (identical(FUN_MOD,RunModel_GR1A)) {
-    Outputs_all <- c(Outputs_all,"PotEvap", "Precip", "Qsim")
+    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = "GR1A")$GR)
   }
   if ("CemaNeige" %in% ObjectClass) {
-    Outputs_all <- c(Outputs_all,"Pliq", "Psol", "SnowPack", "ThermalState", "Gratio", "PotMelt", "Melt", "PliqAndMelt", "Temp", "Gthreshold", "Glocalmax")
+    Outputs_all <- c(Outputs_all, .FortranOutputs(GR = NULL, isCN = TRUE)$CN)
   }
   
   ##check_Outputs_Sim
@@ -358,7 +355,9 @@ CreateRunOptions <- function(FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndP
     
   }
   Outputs_Cal <- Outputs_Cal[!duplicated(Outputs_Cal)]
-  
+  Outputs_Calxxx <- unique(Outputs_Cal[!duplicated(Outputs_Cal)])
+  print(Outputs_Cal)
+  print(Outputs_Calxxx)
   
   
   ##check_MeanAnSolidPrecip
diff --git a/R/RunModel_CemaNeige.R b/R/RunModel_CemaNeige.R
index 13e6321b4c79b3e9dd143ae518a63a5a6e94c5b2..fb1517a5566e0a91ba28a45ff99ef9de565ff28e 100644
--- a/R/RunModel_CemaNeige.R
+++ b/R/RunModel_CemaNeige.R
@@ -11,12 +11,8 @@ RunModel_CemaNeige <- function(InputsModel, RunOptions, Param, IsHyst = FALSE) {
   ## Initialization of variables
   NParam <- ifelse(IsHyst, 4L, 2L)
   NStates <- 4L
-  FortranOutputsCemaNeige <- c("Pliq", "Psol",
-                               "SnowPack", "ThermalState", "Gratio",
-                               "PotMelt", "Melt", "PliqAndMelt", "Temp",
-                               "Gthreshold", "Glocalmax")
-  
-  
+  FortranOutputsCemaNeige <- .FortranOutputs(GR = NULL, isCN = TRUE)$CN
+  print(FortranOutputsCemaNeige)
   
   ## Arguments_check
   if (!inherits(InputsModel, "InputsModel")) {
diff --git a/R/RunModel_CemaNeigeGR4J.R b/R/RunModel_CemaNeigeGR4J.R
index c426f90575c74440b9283d76f9a7cce58cf35e90..6a8533c728aaac53faeab1a0930faff38561378e 100644
--- a/R/RunModel_CemaNeigeGR4J.R
+++ b/R/RunModel_CemaNeigeGR4J.R
@@ -9,9 +9,8 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
   
   NParam <- ifelse(IsHyst, 8L, 6L)
   NStates <- 4L
-    FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp", "Gthreshold", "Glocalmax");
-    FortranOutputsMod <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
-                           "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
+  FortranOutputs <- .FortranOutputs(GR = "GR4J", isCN = TRUE)
+
     
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
@@ -56,8 +55,8 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
       
     ##SNOW_MODULE________________________________________________________________________________##
     if(inherits(RunOptions,"CemaNeige")==TRUE){
-      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputsCemaNeige)); 
-      } else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim);  }
+      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN)); 
+      } else { IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim);  }
       CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
 
 
@@ -87,7 +86,7 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
 
         ##Data_storage
         CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
-        names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
+        names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige];
         IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
         if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
         if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
@@ -103,8 +102,8 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
 
 
     ##MODEL______________________________________________________________________________________##
-      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputsMod)); 
-      } else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim);  }
+      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)){
@@ -142,33 +141,33 @@ RunModel_CemaNeigeGR4J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
                                             verbose = FALSE)
       }
       
-      if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
+      if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputs$GR[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
 
     ##Output_data_preparation
       ##OutputsModel_only
       if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers) );
-        names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); }
+        names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
       ##DatesR_and_OutputsModel_only
       if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
         OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                            lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers) );
-        names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers);      }
+        names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers);      }
       ##OutputsModel_and_SateEnd_only
       if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
                            list(RESULTS$StateEnd) );
-        names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
+        names(OutputsModel) <- c(FortranOutputs$GR[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(RESULTS$StateEnd) );
-        names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
+        names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
 
     ##End
       rm(RESULTS); 
diff --git a/R/RunModel_CemaNeigeGR5J.R b/R/RunModel_CemaNeigeGR5J.R
index fc58b9c9a32668f33d92ebd0350f2647a5b830a2..de0cc94a46ec36e07f33c6e1d5cb04934a019a3c 100644
--- a/R/RunModel_CemaNeigeGR5J.R
+++ b/R/RunModel_CemaNeigeGR5J.R
@@ -9,9 +9,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
   
   NParam <- ifelse(IsHyst, 9L, 7L)
   NStates <- 4L
-    FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp", "Gthreshold", "Glocalmax");
-    FortranOutputsMod <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
-                           "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
+  FortranOutputs <- .FortranOutputs(GR = "GR5J", isCN = TRUE)
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
@@ -56,8 +54,8 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
 
     ##SNOW_MODULE________________________________________________________________________________##
     if(inherits(RunOptions,"CemaNeige")==TRUE){
-      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputsCemaNeige)); 
-      } else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim);  }
+      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN)); 
+      } else { IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim);  }
       CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
 
       
@@ -87,7 +85,7 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
 
         ##Data_storage
         CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
-        names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
+        names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige];
         IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
         if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
         if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
@@ -103,8 +101,8 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
 
 
     ##MODEL______________________________________________________________________________________##
-      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputsMod)); 
-      } else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim);  }
+      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)){
@@ -142,33 +140,33 @@ RunModel_CemaNeigeGR5J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
                                             verbose = FALSE)
       }
       
-      if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
+      if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputs$GR[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
 
     ##Output_data_preparation
       ##OutputsModel_only
       if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers) );
-        names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); }
+        names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
       ##DatesR_and_OutputsModel_only
       if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
         OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                            lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers) );
-        names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers);      }
+        names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers);      }
       ##OutputsModel_and_SateEnd_only
       if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
                            list(RESULTS$StateEnd) );
-        names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
+        names(OutputsModel) <- c(FortranOutputs$GR[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(RESULTS$StateEnd) );
-        names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
+        names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
 
     ##End
       rm(RESULTS); 
diff --git a/R/RunModel_CemaNeigeGR6J.R b/R/RunModel_CemaNeigeGR6J.R
index 862dad0eccb3a5e016a7714ffe3ba00470b3e6fc..006ef0faa7768e1498ede5d0fdf0c5d6326e5d0f 100644
--- a/R/RunModel_CemaNeigeGR6J.R
+++ b/R/RunModel_CemaNeigeGR6J.R
@@ -9,9 +9,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
   
   NParam <- ifelse(IsHyst, 10L, 8L)
   NStates <- 4L
-    FortranOutputsCemaNeige <- c("Pliq","Psol","SnowPack","ThermalState","Gratio","PotMelt","Melt","PliqAndMelt", "Temp", "Gthreshold", "Glocalmax");
-    FortranOutputsMod       <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1",
-			"Rout", "Exch", "AExch1", "AExch2", "AExch", "QR", "QRExp", "Exp", "QD", "Qsim");
+  FortranOutputs <- .FortranOutputs(GR = "GR6J", isCN = TRUE)
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
@@ -60,8 +58,8 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
 
     ##SNOW_MODULE________________________________________________________________________________##
     if(inherits(RunOptions,"CemaNeige")==TRUE){
-      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputsCemaNeige)); 
-      } else { IndOutputsCemaNeige <- which(FortranOutputsCemaNeige %in% RunOptions$Outputs_Sim);  }
+      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsCemaNeige <- as.integer(1:length(FortranOutputs$CN)); 
+      } else { IndOutputsCemaNeige <- which(FortranOutputs$CN %in% RunOptions$Outputs_Sim);  }
       CemaNeigeLayers <- list(); CemaNeigeStateEnd <- NULL; NameCemaNeigeLayers <- "CemaNeigeLayers";
 
       
@@ -91,7 +89,7 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
 
         ##Data_storage
         CemaNeigeLayers[[iLayer]] <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]);
-        names(CemaNeigeLayers[[iLayer]]) <- FortranOutputsCemaNeige[IndOutputsCemaNeige];
+        names(CemaNeigeLayers[[iLayer]]) <- FortranOutputs$CN[IndOutputsCemaNeige];
         IndPliqAndMelt <- which(names(CemaNeigeLayers[[iLayer]]) == "PliqAndMelt");
         if(iLayer==1){ CatchMeltAndPliq <- RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
         if(iLayer >1){ CatchMeltAndPliq <- CatchMeltAndPliq + RESULTS$Outputs[,IndPliqAndMelt]/NLayers; }
@@ -107,8 +105,8 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
 
 
     ##MODEL______________________________________________________________________________________##
-      if("all" %in% RunOptions$Outputs_Sim){ IndOutputsMod <- as.integer(1:length(FortranOutputsMod)); 
-      } else { IndOutputsMod <- which(FortranOutputsMod %in% RunOptions$Outputs_Sim);  }
+      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)){
@@ -147,33 +145,33 @@ RunModel_CemaNeigeGR6J <- function(InputsModel,RunOptions,Param, IsHyst = FALSE)
                                             verbose = FALSE)
       }
       
-      if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputsMod[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
+      if(inherits(RunOptions,"CemaNeige")==TRUE & "Precip" %in% RunOptions$Outputs_Sim){ RESULTS$Outputs[,which(FortranOutputs$GR[IndOutputsMod]=="Precip")] <- InputsModel$Precip[IndPeriod1]; }
 
     ##Output_data_preparation
       ##OutputsModel_only
       if(ExportDatesR==FALSE & ExportStateEnd==FALSE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers) );
-        names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers); }
+        names(OutputsModel) <- c(FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers); }
       ##DatesR_and_OutputsModel_only
       if(ExportDatesR==TRUE & ExportStateEnd==FALSE){
         OutputsModel <- c( list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
                            lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers) );
-        names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers);      }
+        names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers);      }
       ##OutputsModel_and_SateEnd_only
       if(ExportDatesR==FALSE & ExportStateEnd==TRUE){
         OutputsModel <- c( lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2,i]),
                            list(CemaNeigeLayers),
                            list(RESULTS$StateEnd) );
-        names(OutputsModel) <- c(FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
+        names(OutputsModel) <- c(FortranOutputs$GR[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(RESULTS$StateEnd) );
-        names(OutputsModel) <- c("DatesR",FortranOutputsMod[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
+        names(OutputsModel) <- c("DatesR",FortranOutputs$GR[IndOutputsMod],NameCemaNeigeLayers,"StateEnd");      }
 
     ##End
       rm(RESULTS); 
diff --git a/R/RunModel_GR1A.R b/R/RunModel_GR1A.R
index 81f69207f7544fefac330656f1c130ff3d32011f..b05757b650ec2c999f960a573f35410202b14b7c 100644
--- a/R/RunModel_GR1A.R
+++ b/R/RunModel_GR1A.R
@@ -1,7 +1,7 @@
 RunModel_GR1A <- function(InputsModel,RunOptions,Param){
 
     NParam <- 1;
-    FortranOutputs <- c("PotEvap","Precip","Qsim");
+    FortranOutputs <- .FortranOutputs(GR = "GR1A")$GR
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_GR2M.R b/R/RunModel_GR2M.R
index 83fc84de2d346df457c0b10ba4e2628d78cb3416..0dcb2ddfbf7cebc28a222eadefa6d0e0af19f38d 100644
--- a/R/RunModel_GR2M.R
+++ b/R/RunModel_GR2M.R
@@ -1,7 +1,7 @@
 RunModel_GR2M <- function(InputsModel,RunOptions,Param){
 
     NParam <- 2;
-    FortranOutputs <- c("PotEvap", "Precip", "AE", "Pn", "Perc", "PR", "Exch", "Prod", "Rout", "Qsim")
+    FortranOutputs <- .FortranOutputs(GR = "GR2M")$GR
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_GR4H.R b/R/RunModel_GR4H.R
index 81324211d4d42c3941f68c8739f418a98bc0a1af..46cfdc6249f9394c85195e77c84eb92eb23f452a 100644
--- a/R/RunModel_GR4H.R
+++ b/R/RunModel_GR4H.R
@@ -1,7 +1,7 @@
 RunModel_GR4H <- function(InputsModel,RunOptions,Param){
 
     NParam <- 4;
-    FortranOutputs <- c("PotEvap","Precip","Prod","AE","Perc","PR","Q9","Q1","Rout","Exch","AExch","QR","QD","Qsim");
+    FortranOutputs <- .FortranOutputs(GR = "GR4H")$GR
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index c7701e6000529331e193c0564976fd13e71fbcda..8934fa9d415400b60a8f23a0625db575a523a2cf 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -1,8 +1,7 @@
 RunModel_GR4J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 4;
-    FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
-                        "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
+    FortranOutputs <- .FortranOutputs(GR = "GR4J")$GR
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_GR5J.R b/R/RunModel_GR5J.R
index d12be2557c5bbb33e6ea97f635c9eae120239fe3..275d9d7e32506b4541ecf40e656df1c6f6d47758 100644
--- a/R/RunModel_GR5J.R
+++ b/R/RunModel_GR5J.R
@@ -1,8 +1,7 @@
 RunModel_GR5J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 5;
-    FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1", "Rout", "Exch",
-                        "AExch1", "AExch2", "AExch", "QR", "QD", "Qsim");
+    FortranOutputs <- .FortranOutputs(GR = "GR5J")$GR
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/RunModel_GR6J.R b/R/RunModel_GR6J.R
index 5014ba4807db1f51a0a4be98a460a17679eca812..644ab273dd2cc5cec44168403820c2b8b6b5539b 100644
--- a/R/RunModel_GR6J.R
+++ b/R/RunModel_GR6J.R
@@ -1,8 +1,7 @@
 RunModel_GR6J <- function(InputsModel,RunOptions,Param){
 
     NParam <- 6;
-    FortranOutputs <- c("PotEvap", "Precip", "Prod", "Pn", "Ps", "AE", "Perc", "PR", "Q9", "Q1",
-			"Rout", "Exch", "AExch1", "AExch2", "AExch", "QR", "QRExp", "Exp", "QD", "Qsim");
+    FortranOutputs <- .FortranOutputs(GR = "GR6J")$GR
 
     ##Arguments_check
       if(inherits(InputsModel,"InputsModel")==FALSE){ stop("InputsModel must be of class 'InputsModel' \n"); return(NULL); }  
diff --git a/R/Utils.R b/R/Utils.R
new file mode 100644
index 0000000000000000000000000000000000000000..e8c2ef4c6048f93d4376c9bda34ad212ea22c641
--- /dev/null
+++ b/R/Utils.R
@@ -0,0 +1,74 @@
+
+
+## =================================================================================
+## function to manage Fortran outputs
+## =================================================================================
+
+.FortranOutputs <- function(GR = NULL, isCN = FALSE) {
+  
+  outGR <- NULL
+  outCN <- NULL
+  
+  if (is.null(GR)) {
+    GR <- ""
+  }
+  if (GR == "GR1A") {
+    outGR <- c("PotEvap", "Precip",
+               "Qsim")
+  } else if (GR == "GR2M") {
+    outGR <- c("PotEvap", "Precip", "Prod", "Pn",
+               "AE",
+               "Perc", "PR",
+               "Rout", "Exch",
+               "Qsim")
+  } else if (GR == "GR4H") {
+    outGR <- c("PotEvap", "Precip", "Prod",
+               "AE",
+               "Perc", "PR",
+               "Q9", "Q1",
+               "Rout", "Exch",
+               "AExch", "QR",
+               "QD",
+               "Qsim")
+  } else if (GR %in% c("GR4J", "GR5J")) {
+    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
+               "AE",
+               "Perc", "PR",
+               "Q9", "Q1",
+               "Rout", "Exch", 
+               "AExch1", "AExch2",
+               "AExch", "QR",
+               "QD",
+               "Qsim")
+  } else if (GR == "GR6J") {
+    outGR <- c("PotEvap", "Precip", "Prod", "Pn", "Ps",
+               "AE",
+               "Perc", "PR",
+               "Q9", "Q1",
+               "Rout", "Exch",
+               "AExch1", "AExch2",
+               "AExch", "QR",
+               "QRExp", "Exp",
+               "QD",
+               "Qsim")
+  }
+  if (isCN) {
+    outCN <- c("Pliq", "Psol", 
+               "SnowPack", "ThermalState", "Gratio", 
+               "PotMelt", "Melt", "PliqAndMelt", "Temp", 
+               "Gthreshold", "Glocalmax")
+  }
+  
+  res <- list(GR = outGR, CN = outCN)
+  
+}
+
+
+
+
+
+
+
+
+
+
diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f
index ed2bb0fe06af9175dd7b53e3064db9d937c8e5ef..7a2102c15e9872085b928cda173102bc3ec63cd5 100644
--- a/src/frun_GR2M.f
+++ b/src/frun_GR2M.f
@@ -168,13 +168,13 @@ C Updating store level
 C Variables storage
       MISC( 1)=E             ! PE     ! [numeric] observed potential evapotranspiration [mm/month]
       MISC( 2)=P             ! Precip ! [numeric] observed total precipitation  [mm/month]
-      MISC( 3)=AE            ! AE     ! [numeric] actual evapotranspiration [mm/month]
+      MISC( 3)=St(1)         ! Prod   ! [numeric] production store level (St(1)) [mm]
       MISC( 4)=P1            ! Pn     ! [numeric] net rainfall (P1) [mm/month]
-      MISC( 5)=P2            ! Perc   ! [numeric] percolation (P2) [mm/month]
-      MISC( 6)=P3            ! PR     ! [numeric] P3=P1+P2 [mm/month]
-      MISC( 7)=EXCH          ! EXCH   ! [numeric] groundwater exchange (EXCH) [mm/month]
-      MISC( 8)=St(1)         ! Prod   ! [numeric] production store level (St(1)) [mm]
-      MISC( 9)=St(2)         ! Rout   ! [numeric] routing store level (St(2)) [mm]
+      MISC( 5)=AE            ! AE     ! [numeric] actual evapotranspiration [mm/month]
+      MISC( 6)=P2            ! Perc   ! [numeric] percolation (P2) [mm/month]
+      MISC( 7)=P3            ! PR     ! [numeric] P3=P1+P2 [mm/month]
+      MISC( 8)=St(2)         ! Rout   ! [numeric] routing store level (St(2)) [mm]
+      MISC( 9)=EXCH          ! EXCH   ! [numeric] groundwater exchange (EXCH) [mm/month]
       MISC(10)=Q             ! Qsim   ! [numeric] simulated outflow at catchment outlet [mm/month]