diff --git a/DESCRIPTION b/DESCRIPTION index 29584a97d7995587cfd8833c6a9b7eb46cf28e32..6c7c5d7294c27c3f5359075c0a459d8c3dc76d63 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.6.9.19 -Date: 2021-01-12 +Version: 1.6.9.20 +Date: 2021-01-13 Authors@R: c( 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"), diff --git a/NEWS.md b/NEWS.md index e1c3dbea7b84bd793bbf70a16c813d9b8c9ec736..c9374a6d13b470846e00961d4f6eb6deecb59eef 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,7 +2,7 @@ -### 1.6.9.19 Release Notes (2021-01-12) +### 1.6.9.20 Release Notes (2021-01-13) #### New features diff --git a/R/PE_Oudin.R b/R/PE_Oudin.R index a156f3b20a64a9ecb91ea57ee3e77289c4ed0554..2ce49b334736ff8ab69e7ac6823be62d7b4f3684 100644 --- a/R/PE_Oudin.R +++ b/R/PE_Oudin.R @@ -2,10 +2,10 @@ PE_Oudin <- function(JD, Temp, Lat, LatUnit = c("rad", "deg"), TimeStepIn = "daily", TimeStepOut = "daily", RunFortran = FALSE) { - - + + ## ---------- check arguments - + if (!(inherits(JD, "numeric") | inherits(JD, "integer"))) { stop("'JD' must be of class 'numeric'") } @@ -47,76 +47,76 @@ PE_Oudin <- function(JD, Temp, if (TimeStepIn == "hourly" & any(rleJD$lengths != 24)) { stop("each day must have 24 identical values of julian days (one for each hour)") } - - + + ## ---------- hourly inputs aggregation - + if (TimeStepIn == "hourly") { JD <- rleJD$values idJD <- rep(seq_along(JD), each = rleJD$lengths[1L]) Temp <- as.vector(tapply(X = Temp, INDEX = idJD, FUN = mean)) } - - + + ## ---------- Oudin's formula - + if (RunFortran) { - + LInputs = as.integer(length(Temp)) - - if (length(Lat) == 1) { - Lat = rep(Lat, LInputs) + + if (length(FI) == 1) { + FI <- rep(FI, LInputs) } - + RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR", ##inputs LInputs = LInputs, - InputsLAT = as.double(Lat), + InputsLAT = as.double(FI), InputsTT = as.double(Temp), InputsJJ = as.double(JD), - ##outputs + ##outputs PE_Oudin_D = rep(as.double(-999.999), LInputs) ) PE_Oudin_D = RESULTS$PE_Oudin_D - + } else { - + PE_Oudin_D <- rep(NA, length(Temp)) COSFI <- cos(FI) - + 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 { @@ -126,13 +126,13 @@ PE_Oudin <- function(JD, Temp, PE_Oudin_D[k] <- 0 } } - + } - + } - + ## ---------- disaggregate PE from daily to hourly - + if (TimeStepOut == "hourly") { 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, @@ -140,10 +140,10 @@ PE_Oudin <- function(JD, Temp, 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)) } - - + + ## ---------- 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)") @@ -154,12 +154,12 @@ PE_Oudin <- function(JD, Temp, 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) - + }