From 83bf349f09e3889bc65a38baf1f454d56aa23abd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Bourgin?= <francois.bourgin@inrae.fr>
Date: Fri, 17 Apr 2020 18:10:13 +0200
Subject: [PATCH] NEW: ETP Oudin daily in Fortran
---
R/PEdaily_Oudin.R | 122 +++++++++++++++++++++++++-------------------
src/airGR.c | 2 +
src/frun_ETP.f90 | 126 ++++++++++++++++++++++++++++++++++++++++++++++
3 files changed, 199 insertions(+), 51 deletions(-)
create mode 100644 src/frun_ETP.f90
diff --git a/R/PEdaily_Oudin.R b/R/PEdaily_Oudin.R
index b71ac999..b8fcef19 100644
--- a/R/PEdaily_Oudin.R
+++ b/R/PEdaily_Oudin.R
@@ -49,67 +49,87 @@ 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
+ run_fortran = T
+ if (run_fortran)
+ {
+ if (LatUnit[1L] == "rad") {
+ Lat = Lat * 180 / pi
}
+ LInputs = as.integer(length(Temp))
+
+ RESULTS <- .Fortran("frun_etp_oudin", PACKAGE="airGR",
+ ##inputs
+ LInputs = LInputs,
+ xLAT = as.double(Lat),
+ InputsTT = as.double(Temp),
+ InputsJJ = as.double(JD),
+ ##outputs
+ PE_Oudin_D = rep(as.double(-999.999), LInputs)
+ )
+ } else {
- COSOM2 <- COSOM * COSOM
+ PE_Oudin_D <- rep(NA, length(Temp))
+ COSFI <- cos(FI)
+ AFI <- abs(FI / 42)
- if (COSOM2 >= 1) {
- SINOM <- 0
- } else {
- SINOM <- sqrt(1 - COSOM2)
+ 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
+ }
+ }
+
}
- 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 (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)
}
diff --git a/src/airGR.c b/src/airGR.c
index 639675a4..b60f212f 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 00000000..5bde3e9c
--- /dev/null
+++ b/src/frun_ETP.f90
@@ -0,0 +1,126 @@
+!*******************************************************************************
+ SUBROUTINE frun_etp_oudin(LInputs,LAT,InputsTT,InputsJJ,OutputsETP)
+!*******************************************************************************
+
+ !DEC$ ATTRIBUTES DLLEXPORT :: frun_etp_oudin
+
+ Implicit None
+
+ !! dummies
+ ! in
+ integer, intent(in) :: LInputs
+ doubleprecision, intent(in) :: LAT
+ 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
+
+ FI = LAT / 57.296
+
+ !--------------------------------------------------------------
+ ! Time loop
+ !--------------------------------------------------------------
+ DO k=1,LInputs
+ tt = InputsTT(k)
+ jj = InputsJJ(k)
+ ETPoud = -99.9
+ !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
+!*******************************************************************************
--
GitLab