Commit 5107b586 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

Merge branch '62-embed-oudin-s-fortran-code' into 'dev'

Resolve "Embed Oudin's Fortran code"

Closes #62

See merge request !19
parents c5410c88 1d087445
Pipeline #17775 passed with stages
in 35 minutes and 3 seconds
Package: airGR
Type: Package
Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
Version: 1.6.3.68
Version: 1.6.3.72
Date: 2020-11-23
Authors@R: c(
person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
......
......@@ -4,7 +4,7 @@
### 1.6.3.68 Release Notes (2020-11-23)
### 1.6.3.72 Release Notes (2020-11-23)
#### New features
......
......@@ -28,10 +28,10 @@ PE_Oudin <- function(JD, Temp,
if (LatUnit[1L] == "deg" & (all(Lat >= 90) | all(Lat <= -90))) {
stop("'Lat' must be comprised between -90 and +90 degrees")
}
if (!RunFortran & LatUnit[1L] == "rad") {
if (LatUnit[1L] == "rad") {
FI <- Lat
}
if (!RunFortran & LatUnit[1L] == "deg") {
if (LatUnit[1L] == "deg") {
FI <- Lat / (180 / pi)
}
if (any(JD < 0) | any(JD > 366)) {
......@@ -61,9 +61,6 @@ PE_Oudin <- function(JD, Temp,
## ---------- Oudin's formula
if (RunFortran) {
if (LatUnit[1L] == "rad") {
Lat = Lat * 180 / pi
}
LInputs = as.integer(length(Temp))
......@@ -71,7 +68,7 @@ PE_Oudin <- function(JD, Temp,
Lat = rep(Lat, LInputs)
}
RESULTS <- .Fortran("frun_etp_oudin", PACKAGE = "airGR",
RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR",
##inputs
LInputs = LInputs,
InputsLAT = as.double(Lat),
......@@ -86,7 +83,6 @@ PE_Oudin <- function(JD, Temp,
PE_Oudin_D <- rep(NA, length(Temp))
COSFI <- cos(FI)
AFI <- abs(FI / 42)
for (k in seq_along(Temp)) {
......
......@@ -15,7 +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 *);
extern void F77_NAME(frun_pe_oudin)(int *, double *, double *, double *, double *);
static const R_FortranMethodDef FortranEntries[] = {
{"frun_cemaneige", (DL_FUNC) &F77_NAME(frun_cemaneige), 14},
......@@ -26,7 +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},
{"frun_pe_oudin", (DL_FUNC) &F77_NAME(frun_pe_oudin), 5},
{NULL, NULL, 0}
};
......
!------------------------------------------------------------------------------
! Subroutines relative to the Oudin potential evapotranspiration (PE) formula
!------------------------------------------------------------------------------
! TITLE : airGR
! PROJECT : airGR
! FILE : frun_PE.f90
!------------------------------------------------------------------------------
! AUTHORS
! Original code: L. Oudin
! Cleaning and formatting for airGR: Fr. Bourgin
! Further cleaning: O. Delaigue, G. Thirel
!------------------------------------------------------------------------------
! Creation date: 2004
! Last modified: 20/10/2020
!------------------------------------------------------------------------------
! REFERENCES
! 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.
!------------------------------------------------------------------------------
! Quick description of public procedures:
! 1. frun_pe_oudin
! 2. PE_OUDIN
!------------------------------------------------------------------------------
!*******************************************************************************
SUBROUTINE frun_etp_oudin(LInputs,InputsLAT,InputsTT,InputsJJ,OutputsETP)
SUBROUTINE frun_pe_oudin(LInputs,InputsLAT,InputsTemp,InputsJJ,OutputsPE)
!*******************************************************************************
! Subroutine that performs the call to the PE_OUDIN subroutine at each time step,
! and stores the final values
! Inputs
! LInputs ! Integer, length of input and output series
! InputsLAT ! Vector of real, input series of latitude [rad]
! InputsTemp ! Vector of real, input series of air mean temperature [degC]
! InputsJJ ! Vector of real, input series of Julian day [-]
! Outputs
! OutputsPE ! Vector of real, output series of potential evapotranspiration (PE) [mm/time step]
!DEC$ ATTRIBUTES DLLEXPORT :: frun_etp_oudin
!DEC$ ATTRIBUTES DLLEXPORT :: frun_pe_oudin
Implicit None
......@@ -10,27 +49,27 @@
! in
integer, intent(in) :: LInputs
doubleprecision, dimension(LInputs), intent(in) :: InputsLAT
doubleprecision, dimension(LInputs), intent(in) :: InputsTT
doubleprecision, dimension(LInputs), intent(in) :: InputsTemp
doubleprecision, dimension(LInputs), intent(in) :: InputsJJ
! out
doubleprecision, dimension(LInputs), intent(out) :: OutputsETP
doubleprecision, dimension(LInputs), intent(out) :: OutputsPE
!! locals
integer :: k
real :: FI, tt, jj, ETPoud
integer :: k
doubleprecision :: FI, tt, jj, PEoud
!--------------------------------------------------------------
! Time loop
!--------------------------------------------------------------
DO k=1,LInputs
tt = InputsTT(k)
DO k = 1, LInputs
tt = InputsTemp(k)
jj = InputsJJ(k)
FI = InputsLAT(k) / 57.296
FI = InputsLAT(k)!
!model run on one time step
CALL ETP_OUDIN(FI,tt,jj,ETPoud)
CALL PE_OUDIN(FI, tt, jj, PEoud)
!storage of outputs
OutputsETP(k) = ETPoud
OutputsPE(k) = PEoud
ENDDO
RETURN
......@@ -38,14 +77,22 @@
ENDSUBROUTINE
!################################################################################################################################
!*******************************************************************************
SUBROUTINE ETP_OUDIN(FI,DT,JD,DPE)
SUBROUTINE PE_OUDIN(FI,DT,JD,DPE)
!*******************************************************************************
! This subroutine calculates daily potential evapotranspiration (DPE)
! using daily temperature and daily extra-atmospheric global radiation
! Calculation of potential evapotranspiration (DPE) on a single time step
! using air temperature and daily extra-atmospheric global radiation
! (that depends only on Julian day)
!
! The PE formula is is that described in:
! The PE formula is 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
......@@ -58,67 +105,68 @@
! of hydrology. Journal of Hydrology 66 (1/4), 1-76.
!
!***************************************************************
! Inputs:
! xLAT: Latitude in decimal degrees
! DT: Temperature in degree C
! JD: Julian day
! Inputs
! FI ! Latitude [rad]
! DT ! Air Temperature [degC]
! JD ! Julian day [-]
!
! Output:
! DPE: Daily potential evapotranspiration in mm
! Outputs
! DPE ! Potential evapotranspiration [mm/time step]
!***************************************************************
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/
!! dummies
! in
doubleprecision, intent(in) :: FI, DT, JD
! out
doubleprecision, intent(out) :: DPE
!! locals
doubleprecision :: COSFI, TETA, COSTETA, COSGZ, GZ, COSGZ2
doubleprecision :: SINGZ, COSOM, COSOM2, SINOM, COSPZ, OM, GE
doubleprecision :: ETA
! 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))
COSGZ=MAX(0.001d0,COS(FI-TETA))
! GZ: Noon angular zenith distance of the sun
GZ=ACOS(COSGZ)
COSGZ2=COSGZ*COSGZ
IF(COSGZ2.GE.1.)THEN
SINGZ=0.
IF(COSGZ2.GE.1.) THEN
SINGZ=0.
ELSE
SINGZ=SQRT(1.-COSGZ2)
SINGZ=SQRT(1.-COSGZ2)
ENDIF
COSOM=1.-COSGZ/COSFI/COSTETA
IF(COSOM.LT.-1.)COSOM=-1.
IF(COSOM.GT.1.)COSOM=1.
IF(COSOM.LT.-1.) COSOM=-1.
IF(COSOM.GT.1.) COSOM=1.
COSOM2=COSOM*COSOM
IF(COSOM2.GE.1.)THEN
SINOM=0.
IF(COSOM2.GE.1.) THEN
SINOM=0.
ELSE
SINOM=SQRT(1.-COSOM2)
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
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)
DPE=MAX(0.d0,GE*(DT+5.)/100./28.5)
RETURN
END SUBROUTINE ETP_OUDIN
END SUBROUTINE PE_OUDIN
!*******************************************************************************
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment