Merge branch 'etp' into OudinFortran

R/PEdaily_Oudin.R

-PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg")) {
+PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg"), run_fortran = F) {
   
   ## ---------- deprecated function
   
@@ -28,9 +28,9 @@ PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg")) {
   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 (!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")
   }
@@ -49,67 +49,93 @@ PEdaily_Oudin <- function(JD, Temp, LatRad, Lat, LatUnit = c("rad", "deg")) {
   
   
   ## ---------- 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
+
+  if (run_fortran)
+  {
+    if (LatUnit[1L] == "rad") {
+      Lat = Lat * 180 / pi
     }
     
-    COSOM2 <- COSOM * COSOM
+    LInputs = as.integer(length(Temp))
     
-    if (COSOM2 >= 1) {
-      SINOM <- 0
-    } else {
-      SINOM <- sqrt(1 - COSOM2)
+    if (length(Lat) == 1) {
+      Lat = rep(Lat, LInputs)
     }
     
-    OM <- acos(COSOM)
-    COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1)
+    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 {
     
-    if (COSPZ < 0.001) {
-      COSPZ <- 0.001
-    }
+    PE_Oudin_D <- rep(NA, length(Temp))
+    COSFI <- cos(FI)
+    AFI <- abs(FI / 42) 
     
-    ETA <- 1 + cos(JD[k] / 58.1) / 30
-    GE <- 446 * OM * COSPZ * ETA
+    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
+        }
+      }
+      
+    }
     
-    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
+    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 (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)")
-  }
-  
   return(PE_Oudin_D)
   
 }

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