Line data Source code
1 : !> \file mo_temporal_aggregation.f90
2 : !> \brief \copybrief mo_temporal_aggregation
3 : !> \details \copydetails mo_temporal_aggregation
4 :
5 : !> \brief Temporal aggregation for time series (averaging)
6 : !> \details This module does temporal aggregation (averaging) of time series
7 : !> \changelog
8 : !! - Pallav Shrestha, Jun 2019
9 : !! - changed the output argument I/O from INOUT to OUT
10 : !> \authors Oldrich Rakovec, Rohini Kumar, Pallav Shrestha
11 : !> \date October 2015
12 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
13 : !! FORCES is released under the LGPLv3+ license \license_note
14 : MODULE mo_temporal_aggregation
15 :
16 : use mo_kind, ONLY : i4, dp
17 : use mo_julian, ONLY : julday, dec2date
18 : use mo_constants, ONLY : eps_dp
19 :
20 : IMPLICIT NONE
21 :
22 : PUBLIC :: day2mon_average ! converts daily time series to monthly
23 : PUBLIC :: hour2day_average ! converts hourly time series to daily
24 : PUBLIC :: day2mon_sum ! converts daily time series to monthly sums
25 :
26 : ! ------------------------------------------------------------------
27 :
28 : !> \brief Day-to-month average (day2mon_average)
29 :
30 : !> \details converts daily time series to monthly
31 :
32 : !> \param[in] "real(sp/dp) :: daily_data(:)" array of daily time series
33 : !> \param[in] "integer(i4) :: year" year of the starting time
34 : !> \param[in] "integer(i4) :: month" month of the starting time
35 : !> \param[in] "integer(i4) :: day" day of the starting time
36 : !> \param[out] "real(sp/dp) :: mon_average(:)" array of monthly averaged values
37 : !> \param[in] "real(sp/dp), optional :: misval" missing value definition
38 : !> \param[in] "logical, optional :: rm_misval" switch to exclude missing values
39 :
40 : !> \authors Oldrich Rakovec, Rohini Kumar
41 : !> \date Oct 2015
42 :
43 : INTERFACE day2mon_average
44 : MODULE PROCEDURE day2mon_average_dp
45 : END INTERFACE day2mon_average
46 :
47 : ! ------------------------------------------------------------------
48 :
49 : !> \brief Hour-to-day average (hour2day_average)
50 :
51 : !> \details converts hourly time series to daily
52 :
53 : !> \param[in] "real(sp/dp) :: hourly_data(:)" array of hourly time series
54 : !> \param[in] "integer(i4) :: year" year of the starting time
55 : !> \param[in] "integer(i4) :: month" month of the starting time
56 : !> \param[in] "integer(i4) :: day" day of the starting time
57 : !> \param[in] "integer(i4) :: hour" hour of the starting time
58 : !> \param[inout] "real(sp/dp) :: day_average(:)" array of daily averaged values
59 : !> \param[in] "real(sp/dp), optional :: misval" missing value definition
60 : !> \param[in] "logical, optional :: rm_misval" switch to exclude missing values
61 :
62 : !> \note Hours values should be from 0 to 23 (NOT from 1 to 24!)
63 :
64 : !> \author Oldrich Rakovec, Rohini Kumar
65 : !> \date Oct 2015
66 :
67 : INTERFACE hour2day_average
68 : MODULE PROCEDURE hour2day_average_dp
69 : END INTERFACE hour2day_average
70 :
71 : ! ------------------------------------------------------------------
72 :
73 : !> \brief Day-to-month sum (day2mon_sum)
74 :
75 : !> \details converts daily time series to monthly sums
76 :
77 : !> \param[in] "real(sp/dp) :: daily_data(:)" array of daily time series
78 : !> \param[in] "integer(i4) :: year" year of the starting time
79 : !> \param[in] "integer(i4) :: month" month of the starting time
80 : !> \param[in] "integer(i4) :: day" day of the starting time
81 : !> \param[out] "real(sp/dp) :: mon_sum(:)" array of monthly summed values
82 : !> \param[in] "real(sp/dp), optional :: misval" missing value definition
83 : !> \param[in] "logical, optional :: rm_misval" switch to exclude missing values
84 :
85 : !> \author Pallav Kumar Shrestha
86 : !> \date Apr 2019
87 :
88 : INTERFACE day2mon_sum
89 : MODULE PROCEDURE day2mon_sum_dp
90 : END INTERFACE day2mon_sum
91 :
92 : ! ------------------------------------------------------------------
93 :
94 : PRIVATE
95 :
96 : ! ------------------------------------------------------------------
97 :
98 : CONTAINS
99 :
100 2 : SUBROUTINE day2mon_average_dp(daily_data, yearS, monthS, dayS, mon_avg, misval, rm_misval)
101 :
102 : IMPLICIT NONE
103 :
104 : REAL(dp), dimension(:), INTENT(IN) :: daily_data ! array of daily data
105 : INTEGER(i4), INTENT(IN) :: yearS ! year of the initial time step
106 : INTEGER(i4), INTENT(IN) :: monthS ! month of the initial time step
107 : INTEGER(i4), INTENT(IN) :: dayS ! day of the initial time step
108 :
109 : REAL(dp), dimension(:), allocatable, INTENT(OUT) :: mon_avg ! array of the monthly averages
110 :
111 : REAL(dp), optional, INTENT(IN) :: misval ! missing value definition
112 : logical, optional, INTENT(IN) :: rm_misval ! switch to remove missing values
113 :
114 : ! local variables
115 : INTEGER(i4) :: ndays, tt, kk ! number of days, indices
116 : INTEGER(i4) :: start_day, end_day ! size of input array, size of days
117 : INTEGER(i4) :: y, m
118 : INTEGER(i4) :: year, month, day ! variables for date
119 : INTEGER(i4) :: yearE, monthE, dayE ! vatiables for End date
120 2 : REAL(dp) :: newTime
121 :
122 2 : REAL(dp), dimension(:, :), allocatable :: nCounter_m ! counter number of days in months (w/ data)
123 2 : REAL(dp), dimension(:, :), allocatable :: nCounter_m_full ! counter number of days in months (complete)
124 2 : REAL(dp), dimension(:, :), allocatable :: mon_sum ! monthly sum
125 :
126 : INTEGER(i4) :: nmonths ! number of days, number of months
127 : LOGICAL :: remove ! switch for considering missing data
128 2 : REAL(dp) :: missing ! switch for reading missing value or default -9999.
129 :
130 2 : if (present(misval)) then
131 0 : missing = misval
132 : else
133 : missing = -9999._dp
134 : end if
135 :
136 2 : if (present(rm_misval)) then
137 0 : remove = rm_misval
138 : else
139 : remove = .FALSE.
140 : end if
141 :
142 : ! get total number of days
143 2 : ndays = SIZE(daily_data)
144 :
145 : ! assign initial julian day
146 2 : start_day = julday(dayS, monthS, yearS)
147 :
148 : ! calculate last julian day
149 2 : end_day = start_day + ndays - 1_i4
150 :
151 : ! get year, month and day of the end date:
152 2 : call dec2date(real(end_day, dp), yy = yearE, mm = monthE, dd = dayE)
153 :
154 : ! get number of days with data for each month
155 8 : allocate(nCounter_m(yearS : yearE, 12))
156 4 : allocate(nCounter_m_full(yearS : yearE, 12))
157 4 : allocate(mon_sum(yearS : yearE, 12))
158 50 : nCounter_m(:, :) = 0
159 50 : nCounter_m_full(:, :) = 0
160 50 : mon_sum(:, :) = 0.0_dp
161 :
162 2 : newTime = real(start_day, dp)
163 : ! calculate monthly sums
164 733 : do tt = 1, (end_day - start_day + 1)
165 731 : call dec2date((newTime + tt - 1), yy = year, mm = month, dd = day)
166 731 : nCounter_m_full(year, month) = nCounter_m_full(year, month) + 1.0_dp
167 731 : if (abs(daily_data(tt) - missing) .lt. eps_dp) cycle
168 731 : mon_sum(year, month) = mon_sum(year, month) + daily_data(tt)
169 733 : nCounter_m(year, month) = nCounter_m(year, month) + 1.0_dp
170 : end do
171 :
172 : ! calculate number of months
173 2 : nmonths = 0
174 4 : do y = yearS, yearE
175 28 : do m = 1, 12
176 24 : if ((y .EQ. yearS) .AND. (m .LT. monthS)) cycle
177 24 : if ((y .EQ. yearE) .AND. (m .GT. monthE)) cycle
178 26 : nmonths = nmonths + 1
179 : end do
180 : end do
181 :
182 : ! calculate monthly averages
183 2 : if(allocated(mon_avg)) deallocate(mon_avg)
184 6 : allocate(mon_avg(nmonths))
185 26 : mon_avg(:) = missing
186 2 : kk = 0
187 4 : do y = yearS, yearE
188 28 : do m = 1, 12
189 24 : if ((y .EQ. yearS) .AND. (m .LT. monthS)) cycle
190 24 : if ((y .EQ. yearE) .AND. (m .GT. monthE)) cycle
191 24 : kk = kk + 1
192 24 : if ((nCounter_m(y, m) .GT. 0) .AND. &
193 2 : (abs(nCounter_m_full(y, m) - nCounter_m(y, m)) .LT. eps_dp)) then
194 24 : mon_avg(kk) = mon_sum(y, m) / real(nCounter_m(y, m), dp)
195 0 : else if ((nCounter_m(y, m) .GT. 0) .AND. remove) then
196 0 : mon_avg(kk) = mon_sum(y, m) / real(nCounter_m(y, m), dp)
197 : end if
198 : end do
199 : end do
200 :
201 2 : deallocate(nCounter_m_full)
202 2 : deallocate(nCounter_m)
203 2 : deallocate(mon_sum)
204 :
205 2 : END SUBROUTINE day2mon_average_dp
206 :
207 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
208 1 : SUBROUTINE hour2day_average_dp(hourly_data, yearS, monthS, dayS, hourS, day_avg, misval, rm_misval)
209 :
210 : IMPLICIT NONE
211 :
212 : REAL(dp), dimension(:), INTENT(IN) :: hourly_data ! array of hourly data
213 : INTEGER(i4), INTENT(IN) :: yearS ! year of the initial time step
214 : INTEGER(i4), INTENT(IN) :: monthS ! month of the initial time step
215 : INTEGER(i4), INTENT(IN) :: dayS ! day of the initial time step
216 : INTEGER(i4), INTENT(IN) :: hourS ! hour of the initial time step
217 :
218 : REAL(dp), dimension(:), allocatable, INTENT(INOUT) :: day_avg ! array of the daily averages
219 :
220 : REAL(dp), optional, INTENT(IN) :: misval ! missing value definition
221 : logical, optional, INTENT(IN) :: rm_misval ! switch to remove missing values
222 :
223 : ! local variables
224 : INTEGER(i4) :: nhours, ndays_dummy, tt, dd, kk
225 1 : REAL(dp) :: start_day, end_day ! assign julian values
226 : INTEGER(i4) :: yearE, monthE, dayE, hourE, hourEd ! End dates, incl. Dummy
227 :
228 1 : REAL(dp), dimension(:), allocatable :: nCounter_h ! counter number of hours in day (w/ data)
229 1 : REAL(dp), dimension(:), allocatable :: nCounter_h_full ! counter number of hours in day (complete)
230 1 : REAL(dp), dimension(:), allocatable :: day_sum ! daily sum
231 :
232 : LOGICAL :: remove ! switch for considering missing data
233 1 : REAL(dp) :: missing ! switch for reading missing value or default -9999.
234 :
235 1 : if (present(misval)) then
236 0 : missing = misval
237 : else
238 : missing = -9999._dp
239 : end if
240 :
241 1 : if (present(rm_misval)) then
242 0 : remove = rm_misval
243 : else
244 : remove = .FALSE.
245 : end if
246 :
247 : ! get total number of hours
248 1 : nhours = SIZE(hourly_data)
249 : ! assign initial julian day
250 1 : start_day = julday(dayS, monthS, yearS) - 0.5_dp + real(hourS, dp) / 24._dp
251 :
252 : ! calculate last julian day
253 1 : end_day = start_day + real(nhours - 1._dp, dp) / 24._dp
254 :
255 : ! get year, month and day of the end date
256 1 : call dec2date(end_day, yy = yearE, mm = monthE, dd = dayE, hh = hourE)
257 :
258 : ! get largerst possible number of calendar days
259 1 : ndays_dummy = ceiling(real(nhours, dp) / 24._dp + 2._dp)
260 :
261 3 : allocate(day_sum(ndays_dummy))
262 2 : allocate(nCounter_h(ndays_dummy))
263 2 : allocate(nCounter_h_full(ndays_dummy))
264 13 : day_sum(:) = 0.0_dp
265 13 : nCounter_h(:) = 0
266 13 : nCounter_h_full(:) = 0
267 :
268 : ! calculate daily sums
269 : dd = 1
270 241 : do tt = 1, nhours
271 240 : call dec2date(start_day + real(tt - 1, dp) / 24._dp, hh = hourEd)
272 240 : nCounter_h_full(dd) = nCounter_h_full(dd) + 1
273 240 : if (abs(hourly_data(tt) - missing) .lt. eps_dp) then
274 : day_sum(dd) = day_sum(dd)
275 : else
276 240 : day_sum(dd) = day_sum(dd) + hourly_data(tt)
277 240 : nCounter_h(dd) = nCounter_h(dd) + 1
278 : end if
279 241 : if ((hourEd .EQ. 23) .AND. (tt .LT. nhours)) dd = dd + 1
280 : end do
281 :
282 : ! dd is the total number of calendar days, between hourS and hourE
283 3 : allocate(day_avg(dd))
284 11 : day_avg(:) = missing
285 :
286 : ! calculate daily average
287 11 : do kk = 1, dd
288 10 : if ((nCounter_h(kk) .GT. 0) .AND. &
289 1 : (abs(nCounter_h_full(kk) - nCounter_h(kk)) .LT. eps_dp)) then
290 10 : day_avg(kk) = day_sum(kk) / real(nCounter_h(kk), dp)
291 0 : else if ((nCounter_h(kk) .GT. 0) .AND. remove) then
292 0 : day_avg(kk) = day_sum(kk) / real(nCounter_h(kk), dp)
293 : end if
294 : end do
295 :
296 1 : deallocate(nCounter_h_full)
297 1 : deallocate(nCounter_h)
298 1 : deallocate(day_sum)
299 :
300 2 : END SUBROUTINE hour2day_average_dp
301 :
302 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
303 :
304 0 : SUBROUTINE day2mon_sum_dp(daily_data, yearS, monthS, dayS, mon_sum, misval, rm_misval)
305 :
306 : IMPLICIT NONE
307 :
308 : REAL(dp), dimension(:), INTENT(IN) :: daily_data ! array of daily data
309 : INTEGER(i4), INTENT(IN) :: yearS ! year of the initial time step
310 : INTEGER(i4), INTENT(IN) :: monthS ! month of the initial time step
311 : INTEGER(i4), INTENT(IN) :: dayS ! day of the initial time step
312 :
313 : REAL(dp), dimension(:), allocatable, INTENT(OUT) :: mon_sum ! array of the monthly sums
314 :
315 : REAL(dp), optional, INTENT(IN) :: misval ! missing value definition
316 : logical, optional, INTENT(IN) :: rm_misval ! switch to remove missing values
317 :
318 : ! local variables
319 : INTEGER(i4) :: ndays, tt, kk ! number of days, indices
320 : INTEGER(i4) :: start_day, end_day ! size of input array, size of days
321 : INTEGER(i4) :: y, m
322 : INTEGER(i4) :: year, month, day ! variables for date
323 : INTEGER(i4) :: yearE, monthE, dayE ! vatiables for End date
324 0 : REAL(dp) :: newTime
325 :
326 0 : REAL(dp), dimension(:, :), allocatable :: nCounter_m ! counter number of days in months (w/ data)
327 0 : REAL(dp), dimension(:, :), allocatable :: nCounter_m_full ! counter number of days in months (complete)
328 0 : REAL(dp), dimension(:, :), allocatable :: mon_sum_temp_2d ! monthly sum temporary variable
329 0 : REAL(dp), dimension(:), allocatable :: mon_sum_temp_1d ! monthly sum temporary variable
330 :
331 : INTEGER(i4) :: nmonths ! number of days, number of months
332 : LOGICAL :: remove ! switch for considering missing data
333 0 : REAL(dp) :: missing ! switch for reading missing value or default -9999.
334 :
335 0 : if (present(misval)) then
336 0 : missing = misval
337 : else
338 : missing = -9999._dp
339 : end if
340 :
341 0 : if (present(rm_misval)) then
342 0 : remove = rm_misval
343 : else
344 : remove = .FALSE.
345 : end if
346 :
347 : ! get total number of days
348 0 : ndays = SIZE(daily_data)
349 :
350 : ! assign initial julian day
351 0 : start_day = julday(dayS, monthS, yearS)
352 :
353 : ! calculate last julian day
354 0 : end_day = start_day + ndays - 1_i4
355 :
356 : ! get year, month and day of the end date:
357 0 : call dec2date(real(end_day, dp), yy = yearE, mm = monthE, dd = dayE)
358 :
359 : ! get number of days with data for each month
360 0 : allocate(nCounter_m(yearS : yearE, 12))
361 0 : allocate(nCounter_m_full(yearS : yearE, 12))
362 0 : allocate(mon_sum_temp_2d(yearS : yearE, 12))
363 0 : nCounter_m(:, :) = 0
364 0 : nCounter_m_full(:, :) = 0
365 0 : mon_sum_temp_2d(:, :) = 0.0_dp
366 :
367 0 : newTime = real(start_day, dp)
368 : ! calculate monthly sums
369 0 : do tt = 1, (end_day - start_day + 1)
370 0 : call dec2date((newTime + tt - 1), yy = year, mm = month, dd = day)
371 0 : nCounter_m_full(year, month) = nCounter_m_full(year, month) + 1.0_dp
372 0 : if (abs(daily_data(tt) - missing) .lt. eps_dp) cycle
373 0 : mon_sum_temp_2d(year, month) = mon_sum_temp_2d(year, month) + daily_data(tt)
374 0 : nCounter_m(year, month) = nCounter_m(year, month) + 1.0_dp
375 : end do
376 :
377 : ! calculate number of months
378 0 : nmonths = 0
379 0 : do y = yearS, yearE
380 0 : do m = 1, 12
381 0 : if ((y .EQ. yearS) .AND. (m .LT. monthS)) cycle
382 0 : if ((y .EQ. yearE) .AND. (m .GT. monthE)) cycle
383 0 : nmonths = nmonths + 1
384 : end do
385 : end do
386 :
387 :
388 : ! store monthly sums
389 0 : allocate(mon_sum_temp_1d(nmonths))
390 0 : mon_sum_temp_1d(:) = missing
391 0 : kk = 0
392 0 : do y = yearS, yearE
393 0 : do m = 1, 12
394 0 : if ((y .EQ. yearS) .AND. (m .LT. monthS)) cycle
395 0 : if ((y .EQ. yearE) .AND. (m .GT. monthE)) cycle
396 0 : kk = kk + 1
397 0 : if ((nCounter_m(y, m) .GT. 0) .AND. &
398 0 : (abs(nCounter_m_full(y, m) - nCounter_m(y, m)) .LT. eps_dp)) then
399 0 : mon_sum_temp_1d(kk) = mon_sum_temp_2d(y, m)
400 0 : else if ((nCounter_m(y, m) .GT. 0) .AND. remove) then
401 0 : mon_sum_temp_1d(kk) = mon_sum_temp_2d(y, m)
402 : end if
403 : end do
404 : end do
405 :
406 0 : if(allocated(mon_sum)) deallocate(mon_sum)
407 0 : allocate(mon_sum(nmonths))
408 0 : mon_sum = mon_sum_temp_1d
409 :
410 0 : deallocate(nCounter_m_full)
411 0 : deallocate(nCounter_m)
412 0 : deallocate(mon_sum_temp_2d)
413 0 : deallocate(mon_sum_temp_1d)
414 :
415 1 : END SUBROUTINE day2mon_sum_dp
416 :
417 : END MODULE mo_temporal_aggregation
|
