PE_Oudin.R 3.89 KB
Newer Older
1
PE_Oudin <- function(JD, Temp,
2
                     Lat, LatUnit = c("rad", "deg"),
3
4
5
6
7
                     TimeStepIn = "daily", TimeStepOut = "daily") {
  
  
  ## ---------- check arguments
  
8
9
10
11
12
13
  # if (!missing(LatRad)) {
  #   warning("Deprecated 'LatRad' argument. Please, use 'Lat' instead.")
  #   if (missing(Lat)) {
  #     Lat <- LatRad
  #   }
  # }
14
15
16
17
18
19
20
21
22
  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)) {
    stop("'Temp' and 'LatUnit' must have the same length")
  }
23
  LatUnit <- match.arg(LatUnit, choices = c("rad", "deg"))
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
  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")
  }
  if (LatUnit[1L] == "deg" & ((Lat >= 90) | (Lat <= -90))) {
    stop("'Lat' must be  comprised between -90 and +90 degrees")
  }
  if (LatUnit[1L] == "rad") {
    FI <- Lat
  }
  if (LatUnit[1L] == "deg") {
    FI <- Lat / (180 / pi)
  }
  if (any(JD < 0) | any(JD > 366)) {
    stop("'JD' must only contain integers from 1 to 366")
  }
42
43
44
  TimeStep <- c("daily", "hourly")
  TimeStepIn  <- match.arg(TimeStepIn , choices = TimeStep)
  TimeStepOut <- match.arg(TimeStepOut, choices = TimeStep)
45
46
47
48
49
50
  rleJD <- rle(JD)
  if (TimeStepIn == "daily" & any(rleJD$lengths != 1)) {
    stop("each day must have only one identical value of julian days")
  }
  if (TimeStepIn == "hourly" & any(rleJD$lengths != 24)) {
    stop("each day must have 24 identical values of julian days (one for each hour)")
51
  }
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
  
  
  ## ---------- hourly inputs aggregation
  
  if (TimeStepIn == "hourly") {
    JD <- rleJD$values
    idJD <- rep(seq_along(JD), each = rleJD$lengths[1L])
    Temp <- as.vector(tapply(X = Temp, INDEX = idJD, FUN = mean))
  }
  
  
  ## ---------- 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
    }
    
    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
      }
    }
    
  }
  
  
  ## ---------- disaggregate PE from daily to hourly
  
  if (TimeStepOut == "hourly") {
118
    parab_D2H <- c(0.000, 0.000, 0.000, 0.000, 0.000, 0.000,
119
120
                   0.035, 0.062, 0.079, 0.097, 0.110, 0.117,
                   0.117, 0.110, 0.097, 0.079, 0.062, 0.035,
121
                   0.000, 0.000, 0.000, 0.000, 0.000, 0.000)
122
    PE_Oudin_H <- rep(PE_Oudin_D, each = 24) * rep(parab_D2H, times = length(PE_Oudin_D))
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
  }
  
  
  ## ---------- outputs warnings
  
  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 (TimeStepOut == "daily") {
    PE_Oudin <- PE_Oudin_D
  } else {
    PE_Oudin <- PE_Oudin_H
  }
  return(PE_Oudin)
  
}