Line data Source code
1 : !> \file mo_read_timeseries.f90
2 : !> \brief \copybrief mo_read_timeseries
3 : !> \details \copydetails mo_read_timeseries
4 :
5 : !> \brief Routines to read files containing timeseries data.
6 : !> \details This routine is reading time series input data for a particular time period. The files need to have a
7 : !! specific header specified in the different routines.
8 : !> \authors Matthias Zink, Juliane Mai
9 : !> \date Jan 2013
10 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
11 : !! mHM is released under the LGPLv3+ license \license_note
12 : !> \ingroup f_common
13 : MODULE mo_read_timeseries
14 :
15 : USE mo_kind, ONLY : i4, dp
16 : USE mo_os, ONLY : check_path_isfile
17 :
18 : IMPLICIT NONE
19 :
20 : PUBLIC :: read_timeseries
21 :
22 : CONTAINS
23 :
24 : ! -----------------------------------------------------------------
25 : ! NAME
26 : ! read_timeseries
27 :
28 : ! PURPOSE
29 : !> \brief Reads time series in ASCII format.
30 :
31 : !> \details Reads time series in ASCII format.
32 : !> Needs specific header lines:
33 : !> \verbatim
34 : !> <description>
35 : !> nodata <nodata value>
36 : !> n <number of measurements per day> measurements per day [1, 1440]
37 : !> start <YYYY_i4> <MM_i4> <DD_i4> <HH_i4> <MM_i4> (YYYY MM DD HH MM)
38 : !> end <YYYY_i4> <MM_i4> <DD_i4> <HH_i4> <MM_i4> (YYYY MM DD HH MM)
39 : !> \endverbatim
40 : !> Line 6 is the first line with data in the following format:
41 : !> \verbatim
42 : !> <YYYY_i4> <MM_i4> <DD_i4> <HH_i4> <MM_i4> <data_dp>
43 : !> \endverbatim
44 : !> The routine checks for missing data points and if data points are equal distanced.
45 : !> The first data point at each day has to be at HH:MM = 00:00.
46 :
47 : ! INTENT(IN)
48 : !> \param[in] "character(len = *) :: filename" File name
49 : !> \param[in] "integer(i4) :: fileunit" Unit to open file
50 : !> \param[in] "integer(i4), dimension(3) :: periodStart" Start day of reading (YYYY,MM,DD)
51 : !> \param[in] "integer(i4), dimension(3) :: periodEnd" End day of reading (YYYY,MM,DD)
52 : !> \param[in] "logical :: optimize" optimization flag
53 : !> \param[in] "integer(i4) :: opti_function"
54 :
55 : ! INTENT(OUT)
56 : !> \param[out] "real(dp), dimension(:) :: data" Data vector
57 :
58 : ! INTENT(OUT), OPTIONAL
59 : !> \param[out] "logical, dimension(:), optional :: mask" Mask for nodata values in data
60 : !> \param[out] "integer(i4), optional :: nMeasPerDay" Number of data points per day
61 :
62 : ! HISTORY
63 : !> \authors Matthias Zink, Juliane Mai
64 :
65 : !> \date Jan 2013
66 :
67 : ! Modifications:
68 : ! Stephan Thober Mar 2014 - : read data even if shorter than model period
69 : ! MatthiasZink Mar 2014 - : enable read in of nodata periods, e.g. forecast mode
70 : ! Matthias Zink, Juliane Mai Jan 2015 - : corrected read in of nodata periods, e.g. forecast mode
71 :
72 24 : subroutine read_timeseries(filename, fileunit, periodStart, periodEnd, optimize, opti_function, data, mask, &
73 : nMeasPerDay)
74 :
75 : use mo_julian, only : julday
76 : use mo_message, only : error_message
77 : use mo_string_utils, only : num2str
78 :
79 : implicit none
80 :
81 : ! File name
82 : character(len = *), intent(in) :: filename
83 :
84 : ! Unit to open file
85 : integer(i4), intent(in) :: fileunit
86 :
87 : ! Start day of reading (YYYY,MM,DD)
88 : integer(i4), dimension(3), intent(in) :: periodStart
89 :
90 : ! End day of reading (YYYY,MM,DD)
91 : integer(i4), dimension(3), intent(in) :: periodEnd
92 :
93 : ! optimization flag
94 : logical, intent(in) :: optimize
95 :
96 : integer(i4), intent(in) :: opti_function
97 :
98 : ! Data vector
99 : real(dp), dimension(:), allocatable, intent(out) :: data
100 :
101 : ! Mask for nodata values in data
102 : logical, dimension(:), allocatable, optional, intent(out) :: mask
103 :
104 : ! Number of data points per day
105 : integer(i4), optional, intent(out) :: nMeasPerDay
106 :
107 : ! no data value of data
108 24 : real(dp) :: nodata_file
109 :
110 : ! [h] time resolution of input data
111 : integer(i4) :: timestep_file
112 :
113 : ! starting date of timeseries data within file
114 : integer(i4), dimension(3) :: periodStart_file
115 :
116 : ! ending date of timeseries data within file
117 : integer(i4), dimension(3) :: periodEnd_file
118 :
119 : ! year, month, day, hour, minute
120 : integer(i4), dimension(5) :: time_file
121 :
122 : integer(i4) :: i, j, i_max, ios
123 :
124 : ! index to put data from file to data
125 : integer(i4) :: idx_st_period
126 :
127 : ! index to put data from file to data
128 : integer(i4) :: idx_en_period
129 :
130 : ! index to put data from file to data
131 : integer(i4) :: idx_st_file
132 :
133 : ! index to put data from file to data
134 : integer(i4) :: idx_en_file
135 :
136 : ! start julian day of available data
137 : integer(i4) :: startJul_file
138 :
139 : ! end julian day of available data
140 : integer(i4) :: endJul_file
141 :
142 : ! start julian day of needed data
143 : integer(i4) :: startJul_period
144 :
145 : ! end julian day of needed data
146 : integer(i4) :: endJul_period
147 :
148 : ! number of days in file
149 : integer(i4) :: length_file
150 :
151 : ! number of days in period
152 : integer(i4) :: length_period
153 :
154 : ! time series output (fileStart:fileEnd)
155 : ! points --> 0.25 [d-1]
156 24 : real(dp), dimension(:), allocatable :: data_file
157 :
158 : ! dummy for char read in
159 : character(256) :: dummy
160 :
161 :
162 : !checking whether the file exists
163 24 : call check_path_isfile(path = filename, raise=.true.)
164 24 : open(unit = fileunit, file = filename, action = 'read', status = 'old')
165 : ! read header
166 24 : read(fileunit, '(a256)') dummy
167 24 : read(fileunit, *) dummy, nodata_file
168 24 : read(fileunit, *) dummy, timestep_file
169 96 : read(fileunit, *) dummy, (periodStart_file(i), i = 1, 3)
170 96 : read(fileunit, *) dummy, (periodEnd_file(i), i = 1, 3)
171 24 : dummy = dummy // '' ! only to avoid warning
172 24 : if ((timestep_file .lt. 1_i4) .or. (timestep_file .gt. 1440_i4)) then
173 : call error_message('***ERROR: Number of measurements per day has to be between 1 (daily) and 1440 (every minute)! ', &
174 0 : trim(filename))
175 : end if
176 :
177 : ! checking if period is covered by data in file
178 24 : startJul_period = julday(periodStart(3), periodStart(2), periodStart(1))
179 24 : endJul_period = julday(periodEnd(3), periodEnd(2), periodEnd(1))
180 24 : startJul_file = julday(periodStart_file(3), periodStart_file(2), periodStart_file(1))
181 24 : endJul_file = julday(periodEnd_file(3), periodEnd_file(2), periodEnd_file(1))
182 :
183 : if (((startJul_period < startJul_file) .OR. (endJul_period > endJul_file)) &
184 24 : .AND. optimize .and. ((opti_function .le. 9_i4) .or. &
185 : (opti_function .eq. 14_i4) .or. &
186 : (opti_function .eq. 31_i4) .or. &
187 : (opti_function .eq. 33_i4))) then
188 : ! adjust this whenever a new opti function on discharge is added to mhm!
189 0 : call error_message('***ERROR: Simulation period is not covered by observations! ', trim(filename))
190 : end if
191 :
192 : ! allocation of arrays
193 72 : allocate(data ((endJul_period - startJul_period + 1_i4) * timestep_file))
194 10612 : data = nodata_file
195 72 : allocate(data_file((endJul_file - startJul_file + 1_i4) * timestep_file))
196 47082 : data_file = nodata_file
197 :
198 24 : if (present(mask)) then
199 72 : allocate(mask((endJul_period - startJul_period + 1_i4) * timestep_file))
200 10612 : mask = .true.
201 : end if
202 :
203 : ! read data from file to temporal array
204 24 : i_max = (endJul_file - startJul_file + 1_i4) * timestep_file
205 47082 : do i = 1, i_max
206 : ! read date and data
207 282348 : read(fileunit, *, iostat=ios) (time_file(j), j = 1, 5), data_file(i)
208 47058 : if ( ios /= 0 ) call error_message( &
209 : "ERROR: wanted to read ", num2str(i_max), " time-steps from ", &
210 : filename, ", but reached end-of-file at ", num2str(i) &
211 24 : )
212 : end do
213 24 : time_file(1) = time_file(2) + 1 ! only to avoid warning ! PKS - why do we need this?
214 :
215 24 : length_file = (endJul_file - startJul_file + 1)
216 24 : length_period = (endJul_period - startJul_period + 1)
217 :
218 : ! |---------------------------------| FILE
219 : ! |--------------| PERIOD
220 24 : if ((startJul_period .ge. startJul_file) .and. (endJul_period .le. endJul_file)) then
221 23 : idx_st_period = 1
222 23 : idx_en_period = length_period * timestep_file
223 23 : idx_st_file = (startJul_period - startJul_file) * timestep_file + 1
224 23 : idx_en_file = (idx_st_file + length_period - 1) * timestep_file
225 : end if
226 :
227 : ! |--------------| FILE
228 : ! |---------------------------------| PERIOD
229 24 : if ((startJul_period .lt. startJul_file) .and. (endJul_period .gt. endJul_file)) then
230 0 : idx_st_period = (startJul_file - startJul_period) * timestep_file + 1
231 0 : idx_en_period = (idx_st_period + length_file - 1) * timestep_file
232 0 : idx_st_file = 1
233 0 : idx_en_file = length_file * timestep_file
234 : end if
235 :
236 : ! |--------------| FILE
237 : ! |---------------------------------| PERIOD
238 24 : if ((startJul_period .ge. startJul_file) .and. (endJul_period .gt. endJul_file)) then
239 0 : idx_st_period = 1
240 0 : idx_en_period = (endJul_file - startJul_period) * timestep_file + 1
241 0 : idx_st_file = (startJul_period - startJul_file) * timestep_file + 1
242 0 : idx_en_file = length_file * timestep_file
243 : end if
244 :
245 : ! |--------------| FILE
246 : ! |---------------------------------| PERIOD
247 24 : if ((startJul_period .lt. startJul_file) .and. (endJul_period .le. endJul_file)) then
248 1 : idx_st_period = (startJul_file - startJul_period) * timestep_file + 1
249 1 : idx_en_period = (length_period) * timestep_file
250 1 : idx_st_file = 1
251 1 : idx_en_file = (endJul_period - startJul_file) * timestep_file + 1
252 : end if
253 :
254 10247 : data(idx_st_period : idx_en_period) = data_file(idx_st_file : idx_en_file)
255 :
256 24 : if (present(mask)) then
257 10612 : where (abs(data - nodata_file) .lt. tiny(1.0_dp))
258 : mask = .false.
259 : end where
260 : end if
261 :
262 24 : if (present(nMeasPerDay)) then
263 24 : nMeasPerDay = timestep_file
264 : end if
265 :
266 24 : close(fileunit)
267 :
268 24 : end subroutine read_timeseries
269 :
270 : END MODULE mo_read_timeseries
|
