diff --git a/DESCRIPTION b/DESCRIPTION
index c7daa70b755713420ba59666d37703d46594305e..8dfab7e1c23c5de48c76b42226866c4d20ae8a29 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.1.4
+Version: 1.3.2.0
Date: 2019-05-22
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
diff --git a/NEWS.rmd b/NEWS.rmd
index da3f2b091044209dfd563daaef054417fa7e9dd8..805f06f416e43c80345f0e14dad81623fef4f2de 100644
--- a/NEWS.rmd
+++ b/NEWS.rmd
@@ -14,7 +14,7 @@ output:
-### 1.3.1.4 Release Notes (2019-05-22)
+### 1.3.2.0 Release Notes (2019-05-22)
#### New features
@@ -23,6 +23,8 @@ output:
- Added <code>RunModel_CemaNeigeGR4H()</code> function to run the hourly model GR4H with the CemaNeige module.
+- Added <code>PE_Oudin()</code> function to compute Oudin's potential evapotranspiration for hourly or daily tim steps.
+
#### Major user-visible changes
diff --git a/R/PE_Oudin.R b/R/PE_Oudin.R
new file mode 100644
index 0000000000000000000000000000000000000000..e02bfa4751ac7d99acefb85d9f3eed1732a977b5
--- /dev/null
+++ b/R/PE_Oudin.R
@@ -0,0 +1,154 @@
+PE_Oudin <- function(JD, Temp,
+ LatRad, Lat, LatUnit = c("rad", "deg"),
+ TimeStepIn = "daily", TimeStepOut = "daily") {
+
+
+ ## ---------- check arguments
+
+ if (!missing(LatRad)) {
+ warning("Deprecated 'LatRad' argument. Please, use 'Lat' instead.")
+ if (missing(Lat)) {
+ Lat <- LatRad
+ }
+ }
+ if (!(inherits(JD, "numeric") | inherits(JD, "integer"))) {
+ stop("'JD' must be of class 'numeric'")
+ }
+ if (!(inherits(Temp, "numeric") | inherits(Temp, "integer"))) {
+ stop("'Temp' must be of class 'numeric'")
+ }
+ if (length(JD) != length(Temp)) {
+ stop("'Temp' and 'LatUnit' must have the same length")
+ }
+ if (!any(LatUnit %in% c("rad", "deg"))) {
+ stop("'LatUnit' must be one of \"rad\" or \"deg\"")
+ }
+ if (!inherits(Lat, "numeric") | length(Lat) != 1) {
+ stop("'Lat' must be a 'numeric' of length one")
+ }
+ if (LatUnit[1L] == "rad" & ((Lat >= pi/2) | (Lat <= -pi/2))) {
+ stop("'Lat' must be comprised between -pi/2 and +pi/2 degrees")
+ }
+ if (LatUnit[1L] == "deg" & ((Lat >= 90) | (Lat <= -90))) {
+ stop("'Lat' must be comprised between -90 and +90 degrees")
+ }
+ if (LatUnit[1L] == "rad") {
+ FI <- Lat
+ }
+ if (LatUnit[1L] == "deg") {
+ FI <- Lat / (180 / pi)
+ }
+ if (any(JD < 0) | any(JD > 366)) {
+ stop("'JD' must only contain integers from 1 to 366")
+ }
+ if (!inherits(TimeStepIn, "character") | length(TimeStepIn) != 1) {
+ stop("'TimeStepIn' must be a 'character' of length one")
+ }
+ if (!inherits(TimeStepOut, "character") | length(TimeStepOut) != 1) {
+ stop("'TimeStepIn' must be a 'character' of length one")
+ }
+ if (!(TimeStepIn %in% c("daily", "hourly"))) {
+ stop("'TimeStepIn' must be one of \"daily\" or \"hourly\"")
+ }
+ if (!(TimeStepOut %in% c("daily", "hourly"))) {
+ stop("'TimeStepOut' must be one of \"daily\" or \"hourly\"")
+ }
+
+
+ ## ---------- hourly inputs aggregation
+
+ if (TimeStepIn == "hourly") {
+ rleJD <- rle(JD)
+ JD <- rleJD$values
+ if (any(rleJD$lengths != 24)) {
+ stop("each days must have 24 identical Julian day values (one for each hour)")
+ }
+ idJD <- rep(seq_along(JD), each = rleJD$lengths[1L])
+ Temp <- as.vector(tapply(X = Temp, INDEX = idJD, FUN = mean))
+ }
+
+
+ ## ---------- Oudin's formula
+
+ PE_Oudin_D <- rep(NA, length(Temp))
+ COSFI <- cos(FI)
+ AFI <- abs(FI / 42)
+
+ for (k in seq_along(Temp)) {
+
+ TETA <- 0.4093 * sin(JD[k] / 58.1 - 1.405)
+ COSTETA <- cos(TETA)
+ COSGZ <- max(0.001, cos(FI - TETA))
+ GZ <- acos(COSGZ)
+ COSOM <- 1 - COSGZ / COSFI / COSTETA
+
+ if (COSOM < -1) {
+ COSOM <- -1
+ }
+ if (COSOM > 1) {
+ COSOM <- 1
+ }
+
+ COSOM2 <- COSOM * COSOM
+
+ if (COSOM2 >= 1) {
+ SINOM <- 0
+ } else {
+ SINOM <- sqrt(1 - COSOM2)
+ }
+
+ OM <- acos(COSOM)
+ COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1)
+
+ if (COSPZ < 0.001) {
+ COSPZ <- 0.001
+ }
+
+ ETA <- 1 + cos(JD[k] / 58.1) / 30
+ GE <- 446 * OM * COSPZ * ETA
+
+ if (is.na(Temp[k])) {
+ PE_Oudin_D[k] <- NA
+ } else {
+ if (Temp[k] >= -5.0) {
+ PE_Oudin_D[k] <- GE * (Temp[k] + 5) / 100 / 28.5
+ } else {
+ PE_Oudin_D[k] <- 0
+ }
+ }
+
+ }
+
+
+ ## ---------- disaggregate PE from daily to hourly
+
+ if (TimeStepOut == "hourly") {
+ sinus_D2H <- c(0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
+ 0.035, 0.062, 0.079, 0.097, 0.110, 0.117,
+ 0.117, 0.110, 0.097, 0.079, 0.062, 0.035,
+ 0.000, 0.000, 0.000, 0.000, 0.000)
+ PE_Oudin_H <- rep(PE_Oudin_D, each = 24) * rep(sinus_D2H, times = length(PE_Oudin_D))
+ }
+
+
+ ## ---------- outputs warnings
+
+ if (any(is.na(Temp))) {
+ if (any(is.na(PE_Oudin_D))) {
+ warning("'Temp' time series, and therefore the returned 'PE' time series, contain missing value(s)")
+ } else {
+ warning("'Temp' time series contains missing value(s)")
+ }
+ }
+ if (!any(is.na(Temp)) & any(is.na(PE_Oudin_D))) {
+ warning("returned 'PE' time series contains missing value(s)")
+ }
+
+ if (TimeStepOut == "daily") {
+ PE_Oudin <- PE_Oudin_D
+ } else {
+ PE_Oudin <- PE_Oudin_H
+ }
+ return(PE_Oudin)
+
+}