PE_Oudin.R 4.78 KB
Newer Older
1
PE_Oudin <- function(JD, Temp,
2
                     Lat, LatUnit = c("rad", "deg"),
3
                     TimeStepIn = "daily", TimeStepOut = "daily",
4
                     RunFortran = FALSE) {
5
6


7
  ## ---------- check arguments
8

9
10
11
12
13
14
15
  if (!(inherits(JD, "numeric") | inherits(JD, "integer"))) {
    stop("'JD' must be of class 'numeric'")
  }
  if (!(inherits(Temp, "numeric") | inherits(Temp, "integer"))) {
    stop("'Temp' must be of class 'numeric'")
  }
  if (length(JD) != length(Temp)) {
16
    stop("'JD' and 'Temp' must have the same length")
17
  }
18
  if (!RunFortran & (!inherits(Lat, "numeric") | length(Lat) != 1)) {
19
20
    stop("'Lat' must be a 'numeric' of length one")
  }
21
22
23
24
25
  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))) {
26
27
    stop("'Lat' must be comprised between -pi/2 and +pi/2 degrees")
  }
28
  if (LatUnit[1L] == "deg" & (all(Lat >= 90) | all(Lat <= -90))) {
29
    stop("'Lat' must be comprised between -90 and +90 degrees")
30
  }
31
  if (LatUnit[1L] == "rad") {
32
33
    FI <- Lat
  }
34
  if (LatUnit[1L] == "deg") {
35
36
37
38
39
    FI <- Lat / (180 / pi)
  }
  if (any(JD < 0) | any(JD > 366)) {
    stop("'JD' must only contain integers from 1 to 366")
  }
40
41
42
  TimeStep <- c("daily", "hourly")
  TimeStepIn  <- match.arg(TimeStepIn , choices = TimeStep)
  TimeStepOut <- match.arg(TimeStepOut, choices = TimeStep)
43
  rleJD <- rle(JD)
44
45
  msgDaliy <- "each day should have only one identical value of julian days. The time series is not sorted, or contains duplicate or missing dates"
  msgHourly <- "each day must have 24 identical values of julian days (one for each hour). The time series is not sorted, or contains duplicate or missing dates"
46
  if (TimeStepIn == "daily" & any(rleJD$lengths != 1)) {
47
    warning(msgDaliy)
48
  }
49
50


51
  ## ---------- hourly inputs aggregation
52

53
54
55
  if (TimeStepIn == "hourly") {
    JD <- rleJD$values
    idJD <- rep(seq_along(JD), each = rleJD$lengths[1L])
56
57
58
59
60
61
62
63
    if (length(Temp) != length(idJD)) {
      stop(msgHourly)
    } else {
      Temp <- as.vector(tapply(X = Temp, INDEX = idJD, FUN = mean))
    }
  }
  if (TimeStepIn == "hourly" & any(rleJD$lengths != 24)) {
    warning(msgHourly)
64
  }
65
66


67
  ## ---------- Oudin's formula
68

69
  if (RunFortran) {
70

71
    LInputs = as.integer(length(Temp))
72
73

    if (length(FI) == 1) {
74
      FI <- rep(FI, LInputs)
75
    }
76

77
    RESULTS <- .Fortran("frun_pe_oudin", PACKAGE = "airGR",
78
79
                        ##inputs
                        LInputs = LInputs,
80
                        InputsLAT = as.double(FI),
81
82
                        InputsTT = as.double(Temp),
                        InputsJJ = as.double(JD),
83
                        ##outputs
84
                        PE_Oudin_D = rep(as.double(-99e9), LInputs)
85
86
    )
    PE_Oudin_D = RESULTS$PE_Oudin_D
87

88
  } else {
89

90
91
    PE_Oudin_D <- rep(NA, length(Temp))
    COSFI <- cos(FI)
92

93
    for (k in seq_along(Temp)) {
94

95
96
97
98
99
      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
100

101
102
103
104
      if (COSOM < -1) {
        COSOM <- -1
      }
      if (COSOM > 1) {
105
        COSOM <- 1
106
      }
107

108
      COSOM2 <- COSOM * COSOM
109

110
111
      if (COSOM2 >= 1) {
        SINOM <- 0
112
      } else {
113
        SINOM <- sqrt(1 - COSOM2)
114
      }
115

116
117
      OM <- acos(COSOM)
      COSPZ <- COSGZ + COSFI * COSTETA * (SINOM/OM - 1)
118

119
120
121
      if (COSPZ < 0.001) {
        COSPZ <- 0.001
      }
122

123
124
      ETA <- 1 + cos(JD[k] / 58.1) / 30
      GE <- 446 * OM * COSPZ * ETA
125

126
127
128
129
130
131
132
133
134
      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
        }
      }
135

136
    }
137

138
  }
139

140
  ## ---------- disaggregate PE from daily to hourly
141

142
  if (TimeStepOut == "hourly") {
143
    parab_D2H <- c(0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
144
145
                   0.035, 0.062, 0.079, 0.097, 0.110, 0.117,
                   0.117, 0.110, 0.097, 0.079, 0.062, 0.035,
146
                   0.000, 0.000, 0.000, 0.000, 0.000, 0.000)
147
    PE_Oudin_H <- rep(PE_Oudin_D, each = 24) * rep(parab_D2H, times = length(PE_Oudin_D))
148
  }
149
150


151
  ## ---------- outputs warnings
152

153
154
155
156
157
158
159
160
161
162
  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)")
  }
163

164
165
166
167
168
169
  if (TimeStepOut == "daily") {
    PE_Oudin <- PE_Oudin_D
  } else {
    PE_Oudin <- PE_Oudin_H
  }
  return(PE_Oudin)
170

171
}