diff --git a/DESCRIPTION b/DESCRIPTION
index d2d491493f047629dbb001538651c898f25caace..28d6655405fec6cf60b250661a87473352a6c032 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.4.3.91
-Date: 2020-08-26
+Version: 1.4.4.2
+Date: 2020-10-14
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"),
@@ -10,7 +10,7 @@ Authors@R: c(
person("Charles", "Perrin", role = c("aut", "ths"), comment = c(ORCID = "0000-0001-8552-1881")),
person("Claude", "Michel", role = c("aut", "ths")),
person("Vazken", "Andréassian", role = c("ctb", "ths"), comment = c(ORCID = "0000-0001-7124-9303")),
- person("François", "Bourgin", role = c("ctb"), comment = c(ORCID = "0000-0002-2820-7260", vignette = "'Parameter estimation' vignettes")),
+ person("François", "Bourgin", role = c("ctb"), comment = c(ORCID = "0000-0002-2820-7260")),
person("Pierre", "Brigode", role = c("ctb"), comment = c(ORCID = "0000-0001-8257-0741")),
person("Nicolas", "Le Moine", role = c("ctb")),
person("Thibaut", "Mathevet", role = c("ctb"), comment = c(ORCID = "0000-0002-4142-4454")),
diff --git a/NEWS.md b/NEWS.md
index 6862db949f2e3498c450757325949d5cf850e99e..51cdf0480e7dc889511adb2f5fcfd21a9561beac 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,7 +4,12 @@
-### 1.4.3.91 Release Notes (2020-08-26)
+### 1.4.4.2 Release Notes (2020-10-14)
+
+#### New features
+
+- <code>PE_Oudin()</code> now presents a <code>RunFortran</code> argument to run the code in Fortran or in R. The Fortran mode is the fastest.
+
#### Version control and issue tracking
@@ -25,6 +30,8 @@
- Added output to <code>RunModel_GR2M()</code> function (Ps).
+- <code>PE_Oudin()</code> can now run for several locations (i.e. several latitudes) in the Fortran mode (<code>RunFortran = TRUE</code>). In this case <code>Lat</code> must be of the same length as <code>Temp</code>.
+
#### Minor user-visible changes
diff --git a/R/PE_Oudin.R b/R/PE_Oudin.R
index 4113fc63b099aa1ad063a59a9d3c23f88347b74f..07b64ea1f3ae7d5b02a39cdba49398f5bbefd55c 100644
--- a/R/PE_Oudin.R
+++ b/R/PE_Oudin.R
@@ -1,6 +1,7 @@
PE_Oudin <- function(JD, Temp,
Lat, LatUnit = c("rad", "deg"),
- TimeStepIn = "daily", TimeStepOut = "daily") {
+ TimeStepIn = "daily", TimeStepOut = "daily",
+ RunFortran = FALSE) {
## ---------- check arguments
@@ -18,22 +19,26 @@ PE_Oudin <- function(JD, Temp,
stop("'Temp' must be of class 'numeric'")
}
if (length(JD) != length(Temp)) {
- stop("'Temp' and 'LatUnit' must have the same length")
+ stop("'JD' and 'Temp' must have the same length")
}
- LatUnit <- match.arg(LatUnit, choices = c("rad", "deg"))
- if (!inherits(Lat, "numeric") | length(Lat) != 1) {
+
+ if (!RunFortran & (!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))) {
+ if (RunFortran & (!inherits(Lat, "numeric") | (!length(Lat) %in% c(1, length(Temp))))) {
+ stop("'Lat' must be a 'numeric' of length one or of the same length as 'Temp'")
+ }
+ LatUnit <- match.arg(LatUnit, choices = c("rad", "deg"))
+ if (LatUnit[1L] == "rad" & (all(Lat >= pi/2) | all(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] == "deg" & (all(Lat >= 90) | all(Lat <= -90))) {
+ stop("'Lat' must be comprised between -90 and +90 degrees")
}
- if (LatUnit[1L] == "rad") {
+ if (!RunFortran & LatUnit[1L] == "rad") {
FI <- Lat
}
- if (LatUnit[1L] == "deg") {
+ if (!RunFortran & LatUnit[1L] == "deg") {
FI <- Lat / (180 / pi)
}
if (any(JD < 0) | any(JD > 366)) {
@@ -62,6 +67,31 @@ PE_Oudin <- function(JD, Temp,
## ---------- Oudin's formula
+ if (RunFortran)
+ {
+ if (LatUnit[1L] == "rad") {
+ Lat = Lat * 180 / pi
+ }
+
+ LInputs = as.integer(length(Temp))
+
+ if (length(Lat) == 1) {
+ Lat = rep(Lat, LInputs)
+ }
+
+ RESULTS <- .Fortran("frun_etp_oudin", PACKAGE="airGR",
+ ##inputs
+ LInputs = LInputs,
+ InputsLAT = as.double(Lat),
+ InputsTT = as.double(Temp),
+ InputsJJ = as.double(JD),
+ ##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)
AFI <- abs(FI / 42)
@@ -111,6 +141,7 @@ PE_Oudin <- function(JD, Temp,
}
+ }
## ---------- disaggregate PE from daily to hourly
diff --git a/man/PE_Oudin.Rd b/man/PE_Oudin.Rd
index 031d8d1090315ba3ba79485665e784f61275eaaa..6760cbe1ec9b44a1a3c18228aa876cb18712f172 100644
--- a/man/PE_Oudin.Rd
+++ b/man/PE_Oudin.Rd
@@ -10,8 +10,9 @@
\usage{
-PE_Oudin(JD, Temp, Lat, LatUnit,
- TimeStepIn = "daily", TimeStepOut = "daily")
+PE_Oudin(JD, Temp, Lat, LatUnit = c("rad", "deg"),
+ TimeStepIn = "daily", TimeStepOut = "daily",
+ RunFortran = FALSE)
## deprectated function
PEdaily_Oudin(JD, Temp, LatRad, Lat, LatUnit)
@@ -25,13 +26,15 @@ PEdaily_Oudin(JD, Temp, LatRad, Lat, LatUnit)
\item{LatRad}{(deprecated)[numeric] latitude of measurement for the temperature series [radians]. Please use \code{Lat} instead}
-\item{Lat}{[numeric] latitude of measurement for the temperature series [radians or degrees]}
+\item{Lat}{[numeric] latitude of measurement for the temperature series [radians or degrees]. Atomic vector, except if \code{RunFortran = TRUE}, it can be a vector of the same length as \code{Temp}}
\item{LatUnit}{[character] latitude unit (default = \code{"rad"} or \code{"deg"})}
\item{TimeStepIn}{[character] time step of inputs (e.g. \code{"daily"} or \code{"hourly"}, default = \code{"daily"})}
\item{TimeStepOut}{[character] time step of outputs (e.g. \code{"daily"} or \code{"hourly"}, default = \code{"daily"})}
+
+\item{RunFortran}{[boolean] to run the code in the Fortran mode or in the R mode (default)}
}
@@ -70,7 +73,7 @@ PotEvap <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
\author{
-Laurent Coron, Ludovic Oudin, Olivier Delaigue, Guillaume Thirel
+Laurent Coron, Ludovic Oudin, Olivier Delaigue, Guillaume Thirel, François Bourgin
}
diff --git a/src/airGR.c b/src/airGR.c
index 639675a4d94eed29148afc175092f2c1ef95edd7..b60f212f8e3f9ff51abd9f333b05bd4db03816a3 100644
--- a/src/airGR.c
+++ b/src/airGR.c
@@ -15,6 +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 *);
static const R_FortranMethodDef FortranEntries[] = {
{"frun_cemaneige", (DL_FUNC) &F77_NAME(frun_cemaneige), 14},
@@ -25,6 +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},
{NULL, NULL, 0}
};
diff --git a/src/frun_ETP.f90 b/src/frun_ETP.f90
new file mode 100644
index 0000000000000000000000000000000000000000..c0c6f20ddf9dd5936066c43ab5b9ef34ae9e09d0
--- /dev/null
+++ b/src/frun_ETP.f90
@@ -0,0 +1,124 @@
+!*******************************************************************************
+ 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/tests/testthat/test-evap.R b/tests/testthat/test-evap.R
new file mode 100644
index 0000000000000000000000000000000000000000..aac8122d134e540315c3f4c986defad81de0d452
--- /dev/null
+++ b/tests/testthat/test-evap.R
@@ -0,0 +1,62 @@
+context("Test evaporation")
+
+comp_evap <- function(BasinObs,
+ Lat, LatUnit,
+ TimeStepIn = "daily",
+ TimeStepOut = "daily") {
+ PotEvap <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+ Temp = BasinObs$T,
+ Lat = Lat, LatUnit = LatUnit,
+ TimeStepIn = TimeStepIn, TimeStepOut = TimeStepOut)
+ PotEvapFor <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+ Temp = BasinObs$T,
+ Lat = Lat, LatUnit = LatUnit,
+ TimeStepIn = TimeStepIn, TimeStepOut = TimeStepOut,
+ RunFortran = TRUE)
+ all(range(PotEvap - PotEvapFor) < 0.000001)
+}
+
+test_that("PE_Oudin works", {
+ skip_on_cran()
+ rm(list = ls())
+
+ data(L0123001); BasinObs_L0123001 <- BasinObs
+ data(L0123002); BasinObs_L0123002 <- BasinObs
+
+ expect_true(comp_evap(BasinObs = BasinObs_L0123001,
+ Lat = 0.8, LatUnit = "rad",
+ TimeStepIn = "daily", TimeStepOut = "daily"))
+ expect_true(comp_evap(BasinObs = BasinObs_L0123001,
+ Lat = 0.8, LatUnit = "rad",
+ TimeStepIn = "daily", TimeStepOut = "hourly"))
+ expect_true(comp_evap(BasinObs = BasinObs_L0123002,
+ Lat = 0.9, LatUnit = "rad",
+ TimeStepIn = "daily", TimeStepOut = "daily"))
+ expect_true(comp_evap(BasinObs = BasinObs_L0123002,
+ Lat = 0.9, LatUnit = "rad",
+ TimeStepIn = "daily", TimeStepOut = "hourly"))
+
+ ## check with several catchments using different values for Lat
+
+ ## one by one
+ PotEvapFor1 <- PE_Oudin(JD = as.POSIXlt(BasinObs_L0123001$DatesR)$yday + 1,
+ Temp = BasinObs_L0123001$T,
+ Lat = 0.8, LatUnit = "rad",
+ RunFortran = TRUE)
+ PotEvapFor2 <- PE_Oudin(JD = as.POSIXlt(BasinObs_L0123002$DatesR)$yday + 1,
+ Temp = BasinObs_L0123002$T,
+ Lat = 0.9, LatUnit = "rad",
+ RunFortran = TRUE)
+
+ ## all in one
+ BasinObs_L0123001$Lat <- 0.8
+ BasinObs_L0123002$Lat <- 0.9
+ BasinObs <- rbind(BasinObs_L0123001, BasinObs_L0123002)
+ PotEvapFor <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+ Temp = BasinObs$T,
+ Lat = BasinObs$Lat, LatUnit = "rad",
+ RunFortran = TRUE)
+
+ expect_equal(PotEvapFor, c(PotEvapFor1, PotEvapFor2))
+
+})