diff --git a/.Rbuildignore b/.Rbuildignore
index 0e731ea9e35d14b012ea4b826fb7e4c70a139325..88a083344290395ebb1b3f767616fcb2bdc82cf4 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -3,6 +3,7 @@
 ^\.Rprofile$
 ^packrat/
 ^tests/tmp/
+^\.gitlab-ci.yml$
 ^\.regressionignore$
 ^\.gitlab-ci\.yml$
 ^\.vscode$
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 4b1828ba0745a86928138167e60cd1cbb30f4681..599ec847b10e5b32ad7701752284926fa025ca7d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -33,17 +33,11 @@ default:
 
 .regression:
   stage: regression
-  only:
-    refs:
-      - schedules
   script:
-    - Rscript -e 'source("tests/testthat/store_examples.R"); StoreRefExampleResults("airGR");'
+    - Rscript tests/testthat/regression_tests.R stable
     - R CMD INSTALL .
-    - Rscript -e 'source("tests/testthat/store_examples.R"); StoreTestExampleResults("airGR");'
-  artifacts:
-    paths:
-      - tests/tmp/
-    expire_in: 1 week
+    - Rscript tests/testthat/regression_tests.R dev
+    - Rscript tests/testthat/regression_tests.R compare
 
 .check_not_cran:
   stage: tests
@@ -73,11 +67,17 @@ regression_patched:
   extends: .regression
 
 regression_devel:
+  only:
+    refs:
+      - schedules
   variables:
     R_VERSION: "devel"
   extends: .regression
 
 regression_oldrel:
+  only:
+    refs:
+      - schedules
   variables:
     R_VERSION: "oldrel"
   extends: .regression
diff --git a/.regressionignore b/.regressionignore
index f045b6921506c22997c1b574c55b606d36a1c05f..e88367fe53557b59629d2765ff2b5a90e3b2df1d 100644
--- a/.regressionignore
+++ b/.regressionignore
@@ -5,3 +5,4 @@
 # Example for ignoring OutputsModel variable produced by example("RunModel_GR2M"): RunModel_GR2M OutputsModel
 RunModel_GR2M OutputsModel
 RunModel_GR2M RunOptions
+RunModel_GR1A OutputsModel
diff --git a/NAMESPACE b/NAMESPACE
index 0c106ded63b9cfe7c35db6723b515d7933231b52..7119024b0009bfe35e663a2b7907c5cf1d112a5a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -63,8 +63,8 @@ export(TransfoParam_GR5J)
 export(TransfoParam_GR6J)
 export(TransfoParam_Lag)
 export(plot)
-exportPattern(".FortranOutputs")
-exportPattern(".ErrorCrit")
+export(plot.OutputsModel)
+export(.ErrorCrit)
 
 
 
diff --git a/NEWS.md b/NEWS.md
index 6fdd93f079db1a45f4e07800481dc3f8f578f8e4..83b55d04b2bd5444030795f48012588db2eced46 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -16,7 +16,7 @@
 ____________________________________________________________________________________
 
 
-### 1.6.3.65 Release Notes (2020-11-17)
+### 1.6.3.73 Release Notes (2020-11-24)
 
 #### New features
 
@@ -46,7 +46,11 @@ ________________________________________________________________________________
 
 #### Minor user-visible changes
 
-- <code>RunModel_GR1A()</code> now uses the Fortran version of the model code. This code is no longer duplicated: the R version which was used is removed.
+- The <code>.FortranOutputs()</code> function is no longer exported in the namespace.
+
+- <code>RunModel_GR1A()</code> now uses the Fortran version of the model code. This code is no longer duplicated: the R version which was used is removed. ([#65](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/65))
+
+
 - Character argument verification now use partial matching in <code>PE_Oudin()</code> and <code>SeriesAggreg()</code> functions. ([#37](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/37))
 
 
diff --git a/R/PE_Oudin.R b/R/PE_Oudin.R
index 54198f181dad65856df99b21d9ee01de9d0107b6..a156f3b20a64a9ecb91ea57ee3e77289c4ed0554 100644
--- a/R/PE_Oudin.R
+++ b/R/PE_Oudin.R
@@ -28,10 +28,10 @@ PE_Oudin <- function(JD, Temp,
   if (LatUnit[1L] == "deg" & (all(Lat >= 90) | all(Lat <= -90))) {
     stop("'Lat' must be comprised between -90 and +90 degrees")
   }
-  if (!RunFortran & LatUnit[1L] == "rad") {
+  if (LatUnit[1L] == "rad") {
     FI <- Lat
   }
-  if (!RunFortran & LatUnit[1L] == "deg") {
+  if (LatUnit[1L] == "deg") {
     FI <- Lat / (180 / pi)
   }
   if (any(JD < 0) | any(JD > 366)) {
@@ -61,9 +61,6 @@ PE_Oudin <- function(JD, Temp,
   ## ---------- Oudin's formula
   
   if (RunFortran) {
-    if (LatUnit[1L] == "rad") {
-      Lat = Lat * 180 / pi
-    }
     
     LInputs = as.integer(length(Temp))
     
@@ -71,7 +68,7 @@ PE_Oudin <- function(JD, Temp,
       Lat = rep(Lat, LInputs)
     }
     
-    RESULTS <- .Fortran("frun_etp_oudin", PACKAGE = "airGR",
+    RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR",
                         ##inputs
                         LInputs = LInputs,
                         InputsLAT = as.double(Lat),
@@ -86,7 +83,6 @@ PE_Oudin <- function(JD, Temp,
     
     PE_Oudin_D <- rep(NA, length(Temp))
     COSFI <- cos(FI)
-    AFI <- abs(FI / 42) 
     
     for (k in seq_along(Temp)) {
       
diff --git a/R/Utils.R b/R/Utils.R
index 13465c60afdcd53987a9895bc4dfc2ebe4d2e173..23eae05ef30aeb65e19164ce8ca7b3c9f62d475d 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -100,154 +100,3 @@
   res <- list(GR = outGR, CN = outCN)
 
 }
-
-
-
-## =================================================================================
-## function to manage inputs of specific ErrorCrit_*() functions
-## =================================================================================
-
-.ErrorCrit <- function(InputsCrit, Crit, OutputsModel, warnings) {
-
-  ## Arguments check
-  if (!inherits(InputsCrit, "InputsCrit")) {
-    stop("'InputsCrit' must be of class 'InputsCrit'", call. = FALSE)
-  }
-  if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
-    if (Crit == "RMSE") {
-      stop("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' with RMSE", call. = FALSE)
-    } else {
-      stop(paste0("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' or 'Compo' with ", Crit), call. = FALSE)
-    }
-  }
-
-
-  ## Initialisation
-  CritName <- NA
-  CritVar  <- InputsCrit$VarObs
-  if (InputsCrit$transfo == "") {
-    CritName <- paste0(Crit, "[CritVar]")
-  }
-  if (InputsCrit$transfo %in% c("sqrt", "log", "sort", "boxcox")) {
-    CritName <- paste0(Crit, "[", InputsCrit$transfo, "(CritVar)]")
-  }
-  if (InputsCrit$transfo == "inv") {
-    CritName <- paste0(Crit, "[1/CritVar]")
-  }
-  if (grepl("\\^", InputsCrit$transfo)) {
-    transfoPow <- suppressWarnings(as.numeric(gsub("\\^", "", InputsCrit$transfo)))
-    CritName <- paste0(Crit, "[CritVar^", transfoPow, "]")
-  }
-  CritName  <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName)
-  CritValue <- NA
-  if (Crit %in% c("RMSE")) {
-    CritBestValue <- +1
-    Multiplier    <- +1
-  }
-  if (Crit %in% c("NSE", "KGE", "KGE2")) {
-    CritBestValue <- +1
-    Multiplier    <- -1
-  }
-
-
-  ## Data preparation
-  VarObs <- InputsCrit$Obs
-  VarObs[!InputsCrit$BoolCrit] <- NA
-  if (InputsCrit$VarObs == "Q") {
-    VarSim <- OutputsModel$Qsim
-  }
-  if (InputsCrit$VarObs == "SCA") {
-    VarSim <- rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "Gratio"))
-  }
-  if (InputsCrit$VarObs == "SWE") {
-    VarSim <- rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "SnowPack"))
-  }
-  VarSim[!InputsCrit$BoolCrit] <- NA
-
-
-  ## Data transformation
-  if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) {
-    if (any(VarObs %in% 0)) {
-      warning("zeroes detected in 'Qobs': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
-    }
-    if (any(VarSim %in% 0)) {
-      warning("zeroes detected in 'Qsim': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
-    }
-  }
-  if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon) & !(InputsCrit$transfo == "boxcox")) {
-    VarObs <- VarObs + InputsCrit$epsilon
-    VarSim <- VarSim + InputsCrit$epsilon
-  }
-  if (InputsCrit$transfo == "sqrt") {
-    VarObs <- sqrt(VarObs)
-    VarSim <- sqrt(VarSim)
-  }
-  if (InputsCrit$transfo == "log") {
-    VarObs <- log(VarObs)
-    VarSim <- log(VarSim)
-    VarSim[VarSim      < -1e100] <- NA
-  }
-  if (InputsCrit$transfo == "inv") {
-    VarObs <- 1 / VarObs
-    VarSim <- 1 / VarSim
-    VarSim[abs(VarSim) > 1e+100] <- NA
-  }
-  if (InputsCrit$transfo == "sort") {
-    VarSim[is.na(VarObs)] <- NA
-    VarSim <- sort(VarSim, na.last = TRUE)
-    VarObs <- sort(VarObs, na.last = TRUE)
-    InputsCrit$BoolCrit <-  sort(InputsCrit$BoolCrit, decreasing = TRUE)
-  }
-  if (InputsCrit$transfo == "boxcox") {
-    muTransfoVarObs <- (0.01 * mean(VarObs, na.rm = TRUE))^0.25
-    VarSim <- (VarSim^0.25 - muTransfoVarObs) / 0.25
-    VarObs <- (VarObs^0.25 - muTransfoVarObs) / 0.25
-  }
-  if (grepl("\\^", InputsCrit$transfo)) {
-    VarObs <- VarObs^transfoPow
-    VarSim <- VarSim^transfoPow
-  }
-
-
-  ## TS_ignore
-  TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit
-  Ind_TS_ignore <- which(TS_ignore)
-  if (length(Ind_TS_ignore) == 0) {
-    Ind_TS_ignore <- NULL
-  }
-  if (sum(!TS_ignore) == 0 | (sum(!TS_ignore) == 1 & Crit %in% c("KGE", "KGE2"))) {
-    CritCompute <- FALSE
-  } else {
-    CritCompute <- TRUE
-  }
-  if (inherits(OutputsModel, "hourly")) {
-    WarningTS <- 365
-  }
-  if (inherits(OutputsModel, "daily")) {
-    WarningTS <- 365
-  }
-  if (inherits(OutputsModel, "monthly")) {
-    WarningTS <-  12
-  }
-  if (inherits(OutputsModel, "yearly")) {
-    WarningTS <-   3
-  }
-  if (sum(!TS_ignore) < WarningTS & warnings) {
-    warning("\t criterion computed on less than ", WarningTS, " time-steps", call. = FALSE)
-  }
-
-
-  ## Outputs
-  OutputsCritCheck <- list(WarningTS = WarningTS,
-                           VarObs = VarObs,
-                           VarSim = VarSim,
-                           CritBestValue = CritBestValue,
-                           Multiplier = Multiplier,
-                           CritName = CritName,
-                           CritVar = CritVar,
-                           CritCompute = CritCompute,
-                           TS_ignore = TS_ignore,
-                           Ind_TS_ignore = Ind_TS_ignore)
-}
-
-
diff --git a/R/UtilsErrorCrit.R b/R/UtilsErrorCrit.R
new file mode 100644
index 0000000000000000000000000000000000000000..d16a57e520c8e74a872c036c10fefcaa394d3f9b
--- /dev/null
+++ b/R/UtilsErrorCrit.R
@@ -0,0 +1,147 @@
+## =================================================================================
+## function to manage inputs of specific ErrorCrit_*() functions
+## =================================================================================
+
+.ErrorCrit <- function(InputsCrit, Crit, OutputsModel, warnings) {
+  
+  ## Arguments check
+  if (!inherits(InputsCrit, "InputsCrit")) {
+    stop("'InputsCrit' must be of class 'InputsCrit'", call. = FALSE)
+  }
+  if (inherits(InputsCrit, "Multi") | inherits(InputsCrit, "Compo")) {
+    if (Crit == "RMSE") {
+      stop("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' with RMSE", call. = FALSE)
+    } else {
+      stop(paste0("'InputsCrit' must be of class 'Single'. Use the 'ErrorCrit' function on objects of class 'Multi' or 'Compo' with ", Crit), call. = FALSE)
+    }
+  }
+  
+  
+  ## Initialisation
+  CritName <- NA
+  CritVar  <- InputsCrit$VarObs
+  if (InputsCrit$transfo == "") {
+    CritName <- paste0(Crit, "[CritVar]")
+  }
+  if (InputsCrit$transfo %in% c("sqrt", "log", "sort", "boxcox")) {
+    CritName <- paste0(Crit, "[", InputsCrit$transfo, "(CritVar)]")
+  }
+  if (InputsCrit$transfo == "inv") {
+    CritName <- paste0(Crit, "[1/CritVar]")
+  }
+  if (grepl("\\^", InputsCrit$transfo)) {
+    transfoPow <- suppressWarnings(as.numeric(gsub("\\^", "", InputsCrit$transfo)))
+    CritName <- paste0(Crit, "[CritVar^", transfoPow, "]")
+  }
+  CritName  <- gsub(pattern = "CritVar", replacement = CritVar, x = CritName)
+  CritValue <- NA
+  if (Crit %in% c("RMSE")) {
+    CritBestValue <- +1
+    Multiplier    <- +1
+  }
+  if (Crit %in% c("NSE", "KGE", "KGE2")) {
+    CritBestValue <- +1
+    Multiplier    <- -1
+  }
+  
+  
+  ## Data preparation
+  VarObs <- InputsCrit$Obs
+  VarObs[!InputsCrit$BoolCrit] <- NA
+  if (InputsCrit$VarObs == "Q") {
+    VarSim <- OutputsModel$Qsim
+  }
+  if (InputsCrit$VarObs == "SCA") {
+    VarSim <- rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "Gratio"))
+  }
+  if (InputsCrit$VarObs == "SWE") {
+    VarSim <- rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "SnowPack"))
+  }
+  VarSim[!InputsCrit$BoolCrit] <- NA
+  
+  
+  ## Data transformation
+  if (InputsCrit$transfo %in% c("log", "inv") & is.null(InputsCrit$epsilon) & warnings) {
+    if (any(VarObs %in% 0)) {
+      warning("zeroes detected in 'Qobs': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
+    }
+    if (any(VarSim %in% 0)) {
+      warning("zeroes detected in 'Qsim': the corresponding time-steps will be excluded from the criteria computation if the epsilon argument of 'CreateInputsCrit' = NULL", call. = FALSE)
+    }  
+  }
+  if ("epsilon" %in% names(InputsCrit) & !is.null(InputsCrit$epsilon) & !(InputsCrit$transfo == "boxcox")) {
+    VarObs <- VarObs + InputsCrit$epsilon
+    VarSim <- VarSim + InputsCrit$epsilon
+  }
+  if (InputsCrit$transfo == "sqrt") {
+    VarObs <- sqrt(VarObs)
+    VarSim <- sqrt(VarSim)
+  }
+  if (InputsCrit$transfo == "log") {
+    VarObs <- log(VarObs)
+    VarSim <- log(VarSim)
+    VarSim[VarSim      < -1e100] <- NA
+  }
+  if (InputsCrit$transfo == "inv") {
+    VarObs <- 1 / VarObs
+    VarSim <- 1 / VarSim
+    VarSim[abs(VarSim) > 1e+100] <- NA
+  }
+  if (InputsCrit$transfo == "sort") {
+    VarSim[is.na(VarObs)] <- NA
+    VarSim <- sort(VarSim, na.last = TRUE)
+    VarObs <- sort(VarObs, na.last = TRUE)
+    InputsCrit$BoolCrit <-  sort(InputsCrit$BoolCrit, decreasing = TRUE)
+  }
+  if (InputsCrit$transfo == "boxcox") {
+    muTransfoVarObs <- (0.01 * mean(VarObs, na.rm = TRUE))^0.25
+    VarSim <- (VarSim^0.25 - muTransfoVarObs) / 0.25
+    VarObs <- (VarObs^0.25 - muTransfoVarObs) / 0.25
+  }
+  if (grepl("\\^", InputsCrit$transfo)) {
+    VarObs <- VarObs^transfoPow
+    VarSim <- VarSim^transfoPow
+  }
+  
+  
+  ## TS_ignore
+  TS_ignore <- !is.finite(VarObs) | !is.finite(VarSim) | !InputsCrit$BoolCrit
+  Ind_TS_ignore <- which(TS_ignore)
+  if (length(Ind_TS_ignore) == 0) {
+    Ind_TS_ignore <- NULL
+  }
+  if (sum(!TS_ignore) == 0 | (sum(!TS_ignore) == 1 & Crit %in% c("KGE", "KGE2"))) { 
+    CritCompute <- FALSE
+  } else {
+    CritCompute <- TRUE
+  }
+  if (inherits(OutputsModel, "hourly")) {
+    WarningTS <- 365
+  }
+  if (inherits(OutputsModel, "daily")) {
+    WarningTS <- 365
+  }
+  if (inherits(OutputsModel, "monthly")) {
+    WarningTS <-  12
+  }
+  if (inherits(OutputsModel, "yearly")) {
+    WarningTS <-   3
+  }
+  if (sum(!TS_ignore) < WarningTS & warnings) {
+    warning("\t criterion computed on less than ", WarningTS, " time-steps", call. = FALSE)
+  }
+  
+  
+  ## Outputs
+  OutputsCritCheck <- list(WarningTS = WarningTS,
+                           VarObs = VarObs, 
+                           VarSim = VarSim, 
+                           CritBestValue = CritBestValue, 
+                           Multiplier = Multiplier, 
+                           CritName = CritName, 
+                           CritVar = CritVar, 
+                           CritCompute = CritCompute,
+                           TS_ignore = TS_ignore,
+                           Ind_TS_ignore = Ind_TS_ignore)
+}
+
diff --git a/src/airGR.c b/src/airGR.c
index b60f212f8e3f9ff51abd9f333b05bd4db03816a3..b4c256b242438470673b6a80c2f0b2ff26505b80 100644
--- a/src/airGR.c
+++ b/src/airGR.c
@@ -15,7 +15,7 @@ extern void F77_NAME(frun_gr5h)(int *, double *, double *, int *, double *, int
 extern void F77_NAME(frun_gr4j)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
 extern void F77_NAME(frun_gr5j)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
 extern void F77_NAME(frun_gr6j)(int *, double *, double *, int *, double *, int *, double *, int *, int *, double *, double *);
-extern void F77_NAME(frun_etp_oudin)(int *, double *, double *, double *, double *);
+extern void F77_NAME(frun_pe_oudin)(int *, double *, double *, double *, double *);
 
 static const R_FortranMethodDef FortranEntries[] = {
     {"frun_cemaneige", (DL_FUNC) &F77_NAME(frun_cemaneige), 14},
@@ -26,7 +26,7 @@ static const R_FortranMethodDef FortranEntries[] = {
     {"frun_gr4j",      (DL_FUNC) &F77_NAME(frun_gr4j),      11},
     {"frun_gr5j",      (DL_FUNC) &F77_NAME(frun_gr5j),      11},
     {"frun_gr6j",      (DL_FUNC) &F77_NAME(frun_gr6j),      11},
-    {"frun_etp_oudin", (DL_FUNC) &F77_NAME(frun_etp_oudin), 5},
+    {"frun_pe_oudin", (DL_FUNC) &F77_NAME(frun_pe_oudin), 5},
     {NULL, NULL, 0}
 };
 
diff --git a/src/frun_ETP.f90 b/src/frun_ETP.f90
deleted file mode 100644
index c0c6f20ddf9dd5936066c43ab5b9ef34ae9e09d0..0000000000000000000000000000000000000000
--- a/src/frun_ETP.f90
+++ /dev/null
@@ -1,124 +0,0 @@
-!*******************************************************************************
-      SUBROUTINE frun_etp_oudin(LInputs,InputsLAT,InputsTT,InputsJJ,OutputsETP)
-!*******************************************************************************
-
-      !DEC$ ATTRIBUTES DLLEXPORT :: frun_etp_oudin
-
-      Implicit None
-
-      !! dummies
-      ! in
-      integer, intent(in) :: LInputs
-      doubleprecision, dimension(LInputs), intent(in) :: InputsLAT
-      doubleprecision, dimension(LInputs), intent(in) :: InputsTT
-      doubleprecision, dimension(LInputs), intent(in) :: InputsJJ
-      
-      ! out
-      doubleprecision, dimension(LInputs), intent(out) :: OutputsETP
-    
-      !! locals
-      integer :: k
-      real :: FI, tt, jj, ETPoud
-      
-      !--------------------------------------------------------------
-      ! Time loop
-      !--------------------------------------------------------------
-      DO k=1,LInputs
-        tt = InputsTT(k)
-        jj = InputsJJ(k)
-        FI = InputsLAT(k) / 57.296
-        !model run on one time step
-        CALL ETP_OUDIN(FI,tt,jj,ETPoud)
-        !storage of outputs
-        OutputsETP(k) = ETPoud
-      ENDDO
-
-      RETURN
-
-      ENDSUBROUTINE
-
-
-!*******************************************************************************
-      SUBROUTINE ETP_OUDIN(FI,DT,JD,DPE)
-!*******************************************************************************
-!     This subroutine calculates daily potential evapotranspiration (DPE)
-!     using daily temperature and daily extra-atmospheric global radiation
-!     (that depends only on Julian day)
-!
-!     The PE formula is is that described in:
-!     Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V.,
-!     Anctil, F. and Loumagne, C., 2005. Which potential evapotranspiration
-!     input for a rainfall-runoff model? Part 2 - Towards a simple and
-!     efficient PE model for rainfall-runoff modelling. Journal of Hydrology
-!     303(1-4), 290-306.
-!
-!     For the calculation of extra-atmospheric global radiation, see Appendix C of
-!     the article by Morton, F.I., 1983. Operational estimates of areal
-!     evapotranspiration and their significance to the science and practice
-!     of hydrology. Journal of Hydrology 66 (1/4), 1-76.
-!
-!***************************************************************
-!     Inputs:
-!     xLAT: Latitude in decimal degrees
-!     DT: Temperature in degree C
-!     JD: Julian day
-!
-!     Output:
-!     DPE: Daily potential evapotranspiration in mm
-!***************************************************************
-      IMPLICIT NONE
-
-      REAL :: xLAT, FI, COSFI, TETA, COSTETA, COSGZ, GZ, COSGZ2
-      REAL :: SINGZ, COSOM, COSOM2, SINOM, COSPZ, OM, GE
-      REAL ::  ETA, DPE, DT, JD, RD
-
-!      DATA RD/57.296/
-
-! Calculation of extra-atmospheric global radiation (Appendix C in Morton
-! (1983), Eq. C-6 to C-11, p.60-61)
-! Converts latitude in radians
-!      FI=xLAT/RD
-      COSFI=COS(FI)
-!      AFI=ABS(xLAT/42.)
-
-! TETA: Declination of the sun in radians
-      TETA=0.4093*SIN(JD/58.1-1.405)
-      COSTETA=COS(TETA)
-      COSGZ=MAX(0.001,COS(FI-TETA))
-
-! GZ: Noon angular zenith distance of the sun
-      GZ=ACOS(COSGZ)
-      COSGZ2=COSGZ*COSGZ
-
-      IF(COSGZ2.GE.1.)THEN
-      SINGZ=0.
-      ELSE
-      SINGZ=SQRT(1.-COSGZ2)
-      ENDIF
-
-      COSOM=1.-COSGZ/COSFI/COSTETA
-      IF(COSOM.LT.-1.)COSOM=-1.
-      IF(COSOM.GT.1.)COSOM=1.
-      COSOM2=COSOM*COSOM
-      IF(COSOM2.GE.1.)THEN
-      SINOM=0.
-      ELSE
-      SINOM=SQRT(1.-COSOM2)
-      ENDIF
-
-      OM=ACOS(COSOM)
-! PZ: Average angular zenith distance of the sun
-      COSPZ=COSGZ+COSFI*COSTETA*(SINOM/OM-1.)
-      IF(COSPZ.LT.0.001)COSPZ=0.001
-! ETA: Radius vector of the sun
-      ETA=1.+COS(JD/58.1)/30.
-! GE: extra-atmospheric global radiation
-      GE=446.*OM*COSPZ*ETA
-
-! Daily PE by Oudin et al. (2006) formula:
-      DPE=MAX(0.,GE*(DT+5.)/100./28.5)
-
-      RETURN
-
-	END SUBROUTINE ETP_OUDIN
-!*******************************************************************************
diff --git a/src/frun_GR6J.f90 b/src/frun_GR6J.f90
index 9f52bffaa7835f5edb13597fb070fb81697e1e2a..0fa325770325b1905eb8fe1017c18761b51661cd 100644
--- a/src/frun_GR6J.f90
+++ b/src/frun_GR6J.f90
@@ -296,19 +296,14 @@
       IF(AR.GT.33.)  AR=33.
       IF(AR.LT.-33.) AR=-33.
 
-      IF(AR.GT.7.)THEN
+      IF(AR.GT.7.) THEN
         QRExp=St(3)+Param(6)/EXP(AR)
-      GOTO 3
-      ENDIF
-
-      IF(AR.LT.-7.)THEN
+      ELSEIF(AR.LT.-7.) THEN
         QRExp=Param(6)*EXP(AR)
-      GOTO 3
+      ELSE
+        QRExp=Param(6)*LOG(EXP(AR)+1.)
       ENDIF
 
-      QRExp=Param(6)*LOG(EXP(AR)+1.)
-    3 CONTINUE
-
       St(3)=St(3)-QRExp
 
 ! Runoff from direct branch QD
diff --git a/src/frun_PE.f90 b/src/frun_PE.f90
new file mode 100644
index 0000000000000000000000000000000000000000..be4e84a525dff970c19c1025fe01486d9b5bc298
--- /dev/null
+++ b/src/frun_PE.f90
@@ -0,0 +1,172 @@
+!------------------------------------------------------------------------------
+!    Subroutines relative to the Oudin potential evapotranspiration (PE) formula
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : frun_PE.f90
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: L. Oudin
+! Cleaning and formatting for airGR: Fr. Bourgin
+! Further cleaning: O. Delaigue, G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2004
+! Last modified: 20/10/2020
+!------------------------------------------------------------------------------
+! REFERENCES
+! Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V.,
+! Anctil, F. and Loumagne, C., 2005. Which potential evapotranspiration
+! input for a rainfall-runoff model? Part 2 - Towards a simple and
+! efficient PE model for rainfall-runoff modelling. Journal of Hydrology
+! 303(1-4), 290-306.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. frun_pe_oudin
+!         2. PE_OUDIN
+!------------------------------------------------------------------------------
+
+
+!*******************************************************************************
+      SUBROUTINE frun_pe_oudin(LInputs,InputsLAT,InputsTemp,InputsJJ,OutputsPE)
+!*******************************************************************************
+! Subroutine that performs the call to the PE_OUDIN subroutine at each time step, 
+! and stores the final values
+! Inputs
+!       LInputs      ! Integer, length of input and output series
+!       InputsLAT    ! Vector of real, input series of latitude [rad]
+!       InputsTemp   ! Vector of real, input series of air mean temperature [degC]
+!       InputsJJ     ! Vector of real, input series of Julian day [-]
+! Outputs
+!       OutputsPE    ! Vector of real, output series of potential evapotranspiration (PE) [mm/time step]
+
+
+
+      !DEC$ ATTRIBUTES DLLEXPORT :: frun_pe_oudin
+
+      Implicit None
+
+      !! dummies
+      ! in
+      integer, intent(in) :: LInputs
+      doubleprecision, dimension(LInputs), intent(in) :: InputsLAT
+      doubleprecision, dimension(LInputs), intent(in) :: InputsTemp
+      doubleprecision, dimension(LInputs), intent(in) :: InputsJJ
+      
+      ! out
+      doubleprecision, dimension(LInputs), intent(out) :: OutputsPE
+    
+      !! locals
+      integer         :: k
+      doubleprecision :: FI, tt, jj, PEoud
+      
+      !--------------------------------------------------------------
+      ! Time loop
+      !--------------------------------------------------------------
+      DO k = 1, LInputs
+        tt = InputsTemp(k)
+        jj = InputsJJ(k)
+        FI = InputsLAT(k)! 
+        !model run on one time step
+        CALL PE_OUDIN(FI, tt, jj, PEoud)
+        !storage of outputs
+        OutputsPE(k) = PEoud
+      ENDDO
+
+      RETURN
+
+      ENDSUBROUTINE
+
+
+
+
+
+!################################################################################################################################
+
+
+
+
+!*******************************************************************************
+      SUBROUTINE PE_OUDIN(FI,DT,JD,DPE)
+!*******************************************************************************
+!     Calculation of potential evapotranspiration (DPE) on a single time step
+!     using air temperature and daily extra-atmospheric global radiation
+!     (that depends only on Julian day)
+!
+!     The PE formula is described in:
+!     Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V.,
+!     Anctil, F. and Loumagne, C., 2005. Which potential evapotranspiration
+!     input for a rainfall-runoff model? Part 2 - Towards a simple and
+!     efficient PE model for rainfall-runoff modelling. Journal of Hydrology
+!     303(1-4), 290-306.
+!
+!     For the calculation of extra-atmospheric global radiation, see Appendix C of
+!     the article by Morton, F.I., 1983. Operational estimates of areal
+!     evapotranspiration and their significance to the science and practice
+!     of hydrology. Journal of Hydrology 66 (1/4), 1-76.
+!
+!***************************************************************
+! Inputs
+!       FI           ! Latitude [rad]
+!       DT           ! Air Temperature [degC]
+!       JD           ! Julian day [-]
+!
+! Outputs
+!       DPE          ! Potential evapotranspiration [mm/time step]
+!***************************************************************
+      IMPLICIT NONE
+
+      !! dummies
+      ! in
+      doubleprecision, intent(in) :: FI, DT, JD
+      ! out
+      doubleprecision, intent(out) :: DPE
+      !! locals
+      doubleprecision :: COSFI, TETA, COSTETA, COSGZ, GZ, COSGZ2
+      doubleprecision :: SINGZ, COSOM, COSOM2, SINOM, COSPZ, OM, GE
+      doubleprecision :: ETA
+
+! Calculation of extra-atmospheric global radiation (Appendix C in Morton
+! (1983), Eq. C-6 to C-11, p.60-61)
+      COSFI=COS(FI)
+
+! TETA: Declination of the sun in radians
+      TETA=0.4093*SIN(JD/58.1-1.405)
+      COSTETA=COS(TETA)
+      COSGZ=MAX(0.001d0,COS(FI-TETA))
+
+! GZ: Noon angular zenith distance of the sun
+      GZ=ACOS(COSGZ)
+      COSGZ2=COSGZ*COSGZ
+
+      IF(COSGZ2.GE.1.) THEN
+        SINGZ=0.
+      ELSE
+        SINGZ=SQRT(1.-COSGZ2)
+      ENDIF
+
+      COSOM=1.-COSGZ/COSFI/COSTETA
+      IF(COSOM.LT.-1.) COSOM=-1.
+      IF(COSOM.GT.1.) COSOM=1.
+      COSOM2=COSOM*COSOM
+      IF(COSOM2.GE.1.) THEN
+        SINOM=0.
+      ELSE
+        SINOM=SQRT(1.-COSOM2)
+      ENDIF
+
+      OM=ACOS(COSOM)
+! PZ: Average angular zenith distance of the sun
+      COSPZ=COSGZ+COSFI*COSTETA*(SINOM/OM-1.)
+      IF(COSPZ.LT.0.001) COSPZ=0.001
+! ETA: Radius vector of the sun
+      ETA=1.+COS(JD/58.1)/30.
+! GE: extra-atmospheric global radiation
+      GE=446.*OM*COSPZ*ETA
+
+! Daily PE by Oudin et al. (2006) formula:
+      DPE=MAX(0.d0,GE*(DT+5.)/100./28.5)
+
+      RETURN
+
+      END SUBROUTINE PE_OUDIN
+!*******************************************************************************
diff --git a/tests/testthat/store_examples.R b/tests/testthat/helper_regression.R
similarity index 76%
rename from tests/testthat/store_examples.R
rename to tests/testthat/helper_regression.R
index 24160040b04747114a25f5fc6006c727c5df1839..cbcf9dbfee09f0309b429ad54636123ba988a34a 100644
--- a/tests/testthat/store_examples.R
+++ b/tests/testthat/helper_regression.R
@@ -1,14 +1,16 @@
-StoreRefExampleResults <- function(package, ...) {
+StoreStableExampleResults <- function(
+  package = "airGR", 
+  path = file.path("tests/tmp", Sys.getenv("R_VERSION"), "stable"),
+  ...) {
   install.packages(package, repos = "http://cran.r-project.org")
-  StoreExampleResults(package = package, path = "tests/tmp/ref", ...)
+  StoreExampleResults(package = package, path = path, ...)
 }
 
-StoreTestExampleResults <- function(package, ...) {
-  StoreExampleResults(
-    package = package,
-    path = file.path("tests/tmp", Sys.getenv("R_VERSION", "test")),
-    ...
-  )
+StoreDevExampleResults <- function(
+  package = "airGR", 
+  path = file.path("tests/tmp", Sys.getenv("R_VERSION"), "dev"), 
+  ...) {
+  StoreExampleResults(package = package, path = path, ...)
 }
 
 #' Run examples of a package and store the output variables in RDS files for further testing.
@@ -72,3 +74,8 @@ StoreTopicResults <- function(topic, package, path, run.dontrun = TRUE, run.dont
 
 }
 
+CompareStableDev <- function() {
+  res = testthat::test_file("tests/testthat/regression.R")
+  dRes = as.data.frame(res)
+  if(any(dRes[,"failed"]>0) | any(dRes[,"error"])) quit(status = 1)
+}
diff --git a/tests/testthat/regression.R b/tests/testthat/regression.R
new file mode 100644
index 0000000000000000000000000000000000000000..2ccadba442cf1ab03b3834b7e28a6533e60f21b7
--- /dev/null
+++ b/tests/testthat/regression.R
@@ -0,0 +1,44 @@
+context("Compare example outputs with CRAN")
+
+CompareWithStable <- function(refVarFile, testDir, regIgnore) {
+  v <- data.frame(topic = basename(dirname(refVarFile)),
+                  var = gsub("\\.rds$", "", basename(refVarFile)))
+  if (is.null(regIgnore) || all(apply(regIgnore, 1, function(x) !all(x == v)))) {
+    test_that(paste("Compare", v$topic, v$var), {
+      testVarFile <- paste0(
+        file.path(testDir, v$topic, v$var),
+        ".rds"
+      )
+      expect_true(file.exists(testVarFile))
+      if (file.exists(testVarFile)) {
+        testVar <- readRDS(testVarFile)
+        refVar <- readRDS(refVarFile)
+        expect_equivalent(testVar, refVar)
+      }
+    })
+  }
+}
+
+tmp_path <- file.path("../tmp", Sys.getenv("R_VERSION"));
+
+if (dir.exists(file.path(tmp_path, "stable")) & dir.exists(file.path(tmp_path, "dev"))) {
+  refVarFiles <- list.files(file.path(tmp_path, "stable"), recursive = TRUE, full.names = TRUE)
+  regIgnoreFile <- "../../.regressionignore"
+  if (file.exists(regIgnoreFile)) {
+    message("Using .regressionignore file. The following variables are going to be skipped:")
+    regIgnore <- read.table(file = regIgnoreFile,
+                            sep = " ", header = FALSE, skip = 5,
+                            col.names = c("topic", "var"),
+                            stringsAsFactors = FALSE)
+    apply(regIgnore, 1, function(x) message(x[1], ": ", x[2]))
+  } else {
+    message("File ", file.path(getwd(), regIgnoreFile), " not found")
+    regIgnore <- NULL
+  }
+  lapply(X = refVarFiles, CompareWithStable, testDir = file.path(tmp_path, "dev"), regIgnore = regIgnore)
+} else {
+  stop("Regression tests compared to released version needs that you run the following instructions first:\n",
+       "Rscript tests/testthat/regression_tests.R stable\n",
+       "R CMD INSTALL .\n",
+       "Rscript tests/testthat/regression_tests.R dev")
+}
diff --git a/tests/testthat/regression_tests.R b/tests/testthat/regression_tests.R
new file mode 100644
index 0000000000000000000000000000000000000000..2b59829af88508a115abe5c20df9f0c8a7810107
--- /dev/null
+++ b/tests/testthat/regression_tests.R
@@ -0,0 +1,21 @@
+# Execute Regression test by comparing RD files stored in folders /tests/tmp/ref and /tests/tmp/test
+Args = commandArgs(trailingOnly=TRUE)
+
+source("tests/testthat/helper_regression.R")
+
+lActions = list(
+  stable = StoreStableExampleResults,
+  dev = StoreDevExampleResults,
+  compare = CompareStableDev
+)
+
+if(Args %in% names(lActions)) {
+  lActions[[Args]]()
+} else {
+  stop("This script should be run with one argument in the command line:\n",
+       "`Rscript tests/regression_tests.R [stable|dev|compare]`.\n",
+       "Available arguments are:\n",
+       "- stable: install stable version from CRAN, run and store examples\n", 
+       "- dev: install dev version from current directory, run and store examples\n",
+       "- compare: stored results of both versions")
+}
diff --git a/tests/testthat/test-regression.R b/tests/testthat/test-regression.R
deleted file mode 100644
index 432c561d644017d382aaacbd66b3753a54a33ec8..0000000000000000000000000000000000000000
--- a/tests/testthat/test-regression.R
+++ /dev/null
@@ -1,42 +0,0 @@
-context("Compare example outputs with CRAN")
-
-CompareWithRef <- function(refVarFile, testDir, regIgnore) {
-  v <- data.frame(topic = basename(dirname(refVarFile)),
-                  var = gsub("\\.rds$", "", basename(refVarFile)))
-  if (is.null(regIgnore) || all(apply(regIgnore, 1, function(x) !all(x == v)))) {
-    test_that(paste("Compare", v$topic, v$var), {
-      skip_on_cran()
-      testVarFile <- paste0(
-        file.path("../tmp", Sys.getenv("R_VERSION", "test"), v$topic, v$var),
-        ".rds"
-      )
-      expect_true(file.exists(testVarFile))
-      if (file.exists(testVarFile)) {
-        testVar <- readRDS(testVarFile)
-        refVar <- readRDS(refVarFile)
-        expect_equivalent(testVar, refVar)
-      }
-    })
-  }
-}
-
-if (dir.exists("../tmp/ref") & dir.exists("../tmp/test")) {
-  refVarFiles <- list.files("../tmp/ref", recursive = TRUE, full.names = TRUE)
-  regIgnoreFile <- "../../.regressionignore2"
-  if (file.exists(regIgnoreFile)) {
-    regIgnore <- read.table(file = regIgnoreFile,
-                            sep = " ", header = FALSE, skip = 5,
-                            col.names = c("topic", "var"),
-                            stringsAsFactors = FALSE)
-  } else {
-    regIgnore <- NULL
-  }
-  lapply(X = refVarFiles, CompareWithRef, testDir = "../tmp/test", regIgnore = regIgnore)
-} else {
-  warning("Regression tests compared to released version needs that you run the following instructions first:\n",
-          "Rscript -e 'source(\"tests/testthat/store_examples.R\"); StoreRefExampleResults(\"airGR\");'\n",
-          "R CMD INSTALL .\n",
-          "Rscript -e 'source(\"tests/testthat/store_examples.R\"); StoreTestExampleResults(\"airGR\");'\n")
-}
-
-