diff --git a/R/RunModel_GR4J.R b/R/RunModel_GR4J.R
index 0af262be0d397dcdd5a81533ac8e38bdf4746c22..c3b3dbd8a753644b089c30160ff2ebf207dec7b5 100644
--- a/R/RunModel_GR4J.R
+++ b/R/RunModel_GR4J.R
@@ -3,8 +3,6 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
 
   ## Initialization of variables
   NParam <- 4
-  FortranOutputs <- .FortranOutputs(GR = "GR4J")$GR
-
 
   ## Arguments check
   if (!inherits(InputsModel, "InputsModel")) {
@@ -52,16 +50,11 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
   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))
+    IndOutputs <- as.integer(1:length(RunOptions$FortranOutputs$GR))
   } else {
-    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+    IndOutputs <- which(RunOptions$FortranOutputs$GR %in% RunOptions$Outputs_Sim)
   }
 
-  ## Output data preparation
-  IndPeriod2     <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
-  ExportDatesR   <- "DatesR"   %in% RunOptions$Outputs_Sim
-  ExportStateEnd <- "StateEnd" %in% RunOptions$Outputs_Sim
-
   ## Use of IniResLevels
   if (!is.null(RunOptions$IniResLevels)) {
     RunOptions$IniStates[1] <- RunOptions$IniResLevels[1] * Param[1] ### production store level (mm)
@@ -86,7 +79,7 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
   )
   RESULTS$Outputs[RESULTS$Outputs   <= -99e8] <- NA
   RESULTS$StateEnd[RESULTS$StateEnd <= -99e8] <- NA
-  if (ExportStateEnd) {
+  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
     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,
@@ -97,30 +90,10 @@ RunModel_GR4J <- function(InputsModel, RunOptions, Param) {
   }
 
   ## Output data preparation
-  ## OutputsModel only
-  if (!ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i])
-    names(OutputsModel) <- FortranOutputs[IndOutputs]
-  }
-  ## DatesR and OutputsModel only
-  if (ExportDatesR & !ExportStateEnd) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs])
-  }
-  ## OutputsModel and StateEnd only
-  if (!ExportDatesR & ExportStateEnd) {
-    OutputsModel <- c(lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c(FortranOutputs[IndOutputs], "StateEnd")
-  }
-  ## DatesR and OutputsModel and StateEnd
-  if ((ExportDatesR & ExportStateEnd) | "all" %in% RunOptions$Outputs_Sim) {
-    OutputsModel <- c(list(InputsModel$DatesR[RunOptions$IndPeriod_Run]),
-                      lapply(seq_len(RESULTS$NOutputs), function(i) RESULTS$Outputs[IndPeriod2, i]),
-                      list(RESULTS$StateEnd))
-    names(OutputsModel) <- c("DatesR", FortranOutputs[IndOutputs], "StateEnd")
-  }
+  OutputsModel <- .GetOutputsModel(InputsModel,
+                                   RunOptions,
+                                   RESULTS,
+                                   LInputSeries)
 
   ## End
   rm(RESULTS)
diff --git a/R/UtilsRunModel.R b/R/UtilsRunModel.R
new file mode 100644
index 0000000000000000000000000000000000000000..fb0d45248637dd3a0ffe19f3b4dc3e8a2ce6b21a
--- /dev/null
+++ b/R/UtilsRunModel.R
@@ -0,0 +1,48 @@
+#' Create `OutputsModel` for GR non-Cemaneige models
+#'
+#' @param InputsModel output of [CreateInputsModel]
+#' @param RunOptions output of [CreateRunOptions]
+#' @param RESULTS outputs of [.Fortran]
+#' @param LInputSeries number of time steps of warm-up + run periods
+#' @param CemaNeigeLayers outputs of Cemaneige pre-process
+#'
+#' @return
+#' @noRd
+#'
+.GetOutputsModel <- function(InputsModel,
+                                RunOptions,
+                                RESULTS,
+                                LInputSeries,
+                                CemaNeigeLayers = NULL) {
+
+  IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
+  FortranOutputs <- RunOptions$FortranOutputs$GR
+
+  if ("all" %in% RunOptions$Outputs_Sim) {
+    IndOutputs <- as.integer(1:length(FortranOutputs))
+  } else {
+    IndOutputs <- which(FortranOutputs %in% RunOptions$Outputs_Sim)
+  }
+
+  OutputsModel <- list()
+
+  if("DatesR" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$DatesR <- list(InputsModel$DatesR[RunOptions$IndPeriod_Run])
+  }
+
+  seqOutputs <- seq_len(RESULTS$NOutputs)
+  names(seqOutputs) <- FortranOutputs[IndOutputs]
+
+  OutputsModel <- c(OutputsModel,
+                    lapply(seqOutputs, function(i) RESULTS$Outputs[IndPeriod2, i]))
+
+  if(!is.null(CemaNeigeLayers)) {
+    OutputsModel$CemaNeigeLayers <- CemaNeigeLayers
+  }
+
+  if("StateEnd" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$StateEnd <- RESULTS$StateEnd
+  }
+
+  return(OutputsModel)
+}