Commit 7d0ed682 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

v1.6.9.20 fix(PET): replace misuse of Lat by FI in PE_Oudin

Reviewed-by: @francois.bourgin
Refs #62, d6f27bc0
Showing with 39 additions and 39 deletions
+39 -39
Package: airGR Package: airGR
Type: Package Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.6.9.19 Version: 1.6.9.20
Date: 2021-01-12 Date: 2021-01-13
Authors@R: c( Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")), person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"), person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@inrae.fr"),
......
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
### 1.6.9.19 Release Notes (2021-01-12) ### 1.6.9.20 Release Notes (2021-01-13)
#### New features #### New features
......
...@@ -2,10 +2,10 @@ PE_Oudin <- function(JD, Temp, ...@@ -2,10 +2,10 @@ PE_Oudin <- function(JD, Temp,
Lat, LatUnit = c("rad", "deg"), Lat, LatUnit = c("rad", "deg"),
TimeStepIn = "daily", TimeStepOut = "daily", TimeStepIn = "daily", TimeStepOut = "daily",
RunFortran = FALSE) { RunFortran = FALSE) {
## ---------- check arguments ## ---------- check arguments
if (!(inherits(JD, "numeric") | inherits(JD, "integer"))) { if (!(inherits(JD, "numeric") | inherits(JD, "integer"))) {
stop("'JD' must be of class 'numeric'") stop("'JD' must be of class 'numeric'")
} }
...@@ -47,76 +47,76 @@ PE_Oudin <- function(JD, Temp, ...@@ -47,76 +47,76 @@ PE_Oudin <- function(JD, Temp,
if (TimeStepIn == "hourly" & any(rleJD$lengths != 24)) { if (TimeStepIn == "hourly" & any(rleJD$lengths != 24)) {
stop("each day must have 24 identical values of julian days (one for each hour)") stop("each day must have 24 identical values of julian days (one for each hour)")
} }
## ---------- hourly inputs aggregation ## ---------- hourly inputs aggregation
if (TimeStepIn == "hourly") { if (TimeStepIn == "hourly") {
JD <- rleJD$values JD <- rleJD$values
idJD <- rep(seq_along(JD), each = rleJD$lengths[1L]) idJD <- rep(seq_along(JD), each = rleJD$lengths[1L])
Temp <- as.vector(tapply(X = Temp, INDEX = idJD, FUN = mean)) Temp <- as.vector(tapply(X = Temp, INDEX = idJD, FUN = mean))
} }
## ---------- Oudin's formula ## ---------- Oudin's formula
if (RunFortran) { if (RunFortran) {
LInputs = as.integer(length(Temp)) LInputs = as.integer(length(Temp))
if (length(Lat) == 1) { if (length(FI) == 1) {
Lat = rep(Lat, LInputs) FI <- rep(FI, LInputs)
} }
RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR", RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR",
##inputs ##inputs
LInputs = LInputs, LInputs = LInputs,
InputsLAT = as.double(Lat), InputsLAT = as.double(FI),
InputsTT = as.double(Temp), InputsTT = as.double(Temp),
InputsJJ = as.double(JD), InputsJJ = as.double(JD),
##outputs ##outputs
PE_Oudin_D = rep(as.double(-999.999), LInputs) PE_Oudin_D = rep(as.double(-999.999), LInputs)
) )
PE_Oudin_D = RESULTS$PE_Oudin_D PE_Oudin_D = RESULTS$PE_Oudin_D
} else { } else {
PE_Oudin_D <- rep(NA, length(Temp)) PE_Oudin_D <- rep(NA, length(Temp))
COSFI <- cos(FI) COSFI <- cos(FI)
for (k in seq_along(Temp)) { for (k in seq_along(Temp)) {
TETA <- 0.4093 * sin(JD[k] / 58.1 - 1.405) TETA <- 0.4093 * sin(JD[k] / 58.1 - 1.405)
COSTETA <- cos(TETA) COSTETA <- cos(TETA)
COSGZ <- max(0.001, cos(FI - TETA)) COSGZ <- max(0.001, cos(FI - TETA))
GZ <- acos(COSGZ) GZ <- acos(COSGZ)
COSOM <- 1 - COSGZ / COSFI / COSTETA COSOM <- 1 - COSGZ / COSFI / COSTETA
if (COSOM < -1) { if (COSOM < -1) {
COSOM <- -1 COSOM <- -1
} }
if (COSOM > 1) { if (COSOM > 1) {
COSOM <- 1 COSOM <- 1
} }
COSOM2 <- COSOM * COSOM COSOM2 <- COSOM * COSOM
if (COSOM2 >= 1) { if (COSOM2 >= 1) {
SINOM <- 0 SINOM <- 0
} else { } else {
SINOM <- sqrt(1 - COSOM2) SINOM <- sqrt(1 - COSOM2)
} }
OM <- acos(COSOM) OM <- acos(COSOM)
COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1) COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1)
if (COSPZ < 0.001) { if (COSPZ < 0.001) {
COSPZ <- 0.001 COSPZ <- 0.001
} }
ETA <- 1 + cos(JD[k] / 58.1) / 30 ETA <- 1 + cos(JD[k] / 58.1) / 30
GE <- 446 * OM * COSPZ * ETA GE <- 446 * OM * COSPZ * ETA
if (is.na(Temp[k])) { if (is.na(Temp[k])) {
PE_Oudin_D[k] <- NA PE_Oudin_D[k] <- NA
} else { } else {
...@@ -126,13 +126,13 @@ PE_Oudin <- function(JD, Temp, ...@@ -126,13 +126,13 @@ PE_Oudin <- function(JD, Temp,
PE_Oudin_D[k] <- 0 PE_Oudin_D[k] <- 0
} }
} }
} }
} }
## ---------- disaggregate PE from daily to hourly ## ---------- disaggregate PE from daily to hourly
if (TimeStepOut == "hourly") { if (TimeStepOut == "hourly") {
parab_D2H <- c(0.000, 0.000, 0.000, 0.000, 0.000, 0.000, parab_D2H <- c(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.035, 0.062, 0.079, 0.097, 0.110, 0.117,
...@@ -140,10 +140,10 @@ PE_Oudin <- function(JD, Temp, ...@@ -140,10 +140,10 @@ PE_Oudin <- function(JD, Temp,
0.000, 0.000, 0.000, 0.000, 0.000, 0.000) 0.000, 0.000, 0.000, 0.000, 0.000, 0.000)
PE_Oudin_H <- rep(PE_Oudin_D, each = 24) * rep(parab_D2H, times = length(PE_Oudin_D)) PE_Oudin_H <- rep(PE_Oudin_D, each = 24) * rep(parab_D2H, times = length(PE_Oudin_D))
} }
## ---------- outputs warnings ## ---------- outputs warnings
if (any(is.na(Temp))) { if (any(is.na(Temp))) {
if (any(is.na(PE_Oudin_D))) { if (any(is.na(PE_Oudin_D))) {
warning("'Temp' time series, and therefore the returned 'PE' time series, contain missing value(s)") warning("'Temp' time series, and therefore the returned 'PE' time series, contain missing value(s)")
...@@ -154,12 +154,12 @@ PE_Oudin <- function(JD, Temp, ...@@ -154,12 +154,12 @@ PE_Oudin <- function(JD, Temp,
if (!any(is.na(Temp)) & any(is.na(PE_Oudin_D))) { if (!any(is.na(Temp)) & any(is.na(PE_Oudin_D))) {
warning("returned 'PE' time series contains missing value(s)") warning("returned 'PE' time series contains missing value(s)")
} }
if (TimeStepOut == "daily") { if (TimeStepOut == "daily") {
PE_Oudin <- PE_Oudin_D PE_Oudin <- PE_Oudin_D
} else { } else {
PE_Oudin <- PE_Oudin_H PE_Oudin <- PE_Oudin_H
} }
return(PE_Oudin) return(PE_Oudin)
} }
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment