Merge branch 'OudinFortran' into 'dev'

Resolve "Embed Oudin's Fortran code"

Closes #62

See merge request !14

DESCRIPTION

 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")),

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
 

R/PE_Oudin.R

 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
   

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
 }
 
 

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}
 };
 

src/frun_ETP.f90

+!*******************************************************************************
+      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
+!*******************************************************************************

tests/testthat/test-evap.R

+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))
+  
+})