Line data Source code
1 : !> \file mo_read_nc.f90
2 : !> \brief \copybrief mo_read_nc
3 : !> \details \copydetails mo_read_nc
4 :
5 : !> \brief Reads forcing input data.
6 : !> \details This module is to read forcing input data contained in netcdf files, e.g. temperature, precipitation,
7 : !! total_runoff, lai. Timesteps can be hourly, daily, monthly, and annual. The module provides a subroutine
8 : !! for NetCDF files only. First, the dimensions given are cross-checked with header.txt information. Second,
9 : !! the data of the specified period are read from the specified directory.
10 : !! If the optional lower and/or upper bound for the data values is given, the read data are checked for validity.
11 : !! The program is stopped if any value lies out of range.
12 : !> \changelog
13 : !! - Stephan Thober Sep 2015
14 : !! - separated routines for netcdf files from routines for binary files
15 : !! - Stephan Thober Jan 2017
16 : !! - added reading weights for disaggregation of daily meteorological values to hourly ones
17 : !! - Robert Schweppe Nov 2017
18 : !! - switched to mo_netcdf library and restuctured routines
19 : !! - Robert Schweppe Jun 2018
20 : !! - refactoring and reformatting
21 : !! - Sebastian Müller Mar 2023
22 : !! - documentation update
23 : !! - added bound_error to control if an error is raised if boundaries are violated
24 : !> \authors Juliane Mai
25 : !> \date Dec 2012
26 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
27 : !! mHM is released under the LGPLv3+ license \license_note
28 : !> \ingroup f_common
29 : module mo_read_nc
30 : use mo_message, only: error_message
31 :
32 : implicit none
33 : public :: read_nc
34 : public :: read_const_nc
35 : public :: read_weights_nc
36 : private
37 : !
38 : contains
39 :
40 : !> \brief Reads forcing input in NetCDF file format.
41 : !> \details Reads netCDF forcing files.
42 : !! First, the dimensions given are cross-checked with header.txt information. Second, the data of the
43 : !! specified period are read from the specified directory.
44 : !! If the optional lower and/or upper bound for the data values is given, the read data are checked for
45 : !! validity.
46 : !! The program is stopped if any value lies out of range.
47 : !! If the optinal argument nocheck is true, the data are not checked for coverage with the input mask.
48 : !! Additionally in this case an mask of vild data points can be received from the routine in maskout.
49 : !> \changelog
50 : !! - Stephan Thober Nov 2013
51 : !! - only read required chunk from nc file
52 : !! - Matthias Cuntz & Juliane Mai Nov 2014 - read daily, monthly or yearly files
53 : !! - Matthias Zink Mar 2014
54 : !! - added optional nocheck flag and optional maskout
55 : !! - Stephan Thober Sep 2015
56 : !! - added read for hourly data
57 : !! - Robert Schweppe Nov 2017
58 : !! - switched to mo_netcdf library and restuctured routines
59 : !! - Robert Schweppe Jun 2018
60 : !! - refactoring and reformatting
61 : !! - Sebastian Müller Mar 2023
62 : !! - added bound_error to control if an error is raised if boundaries are violated
63 : !! - add nTstepForcingDay as argument to be independent of global variables
64 : !> \authors Matthias Zink
65 : !> \date May 2013
66 802 : subroutine read_nc(folder, nRows, nCols, varName, mask, data, target_period, lower, upper, nctimestep, &
67 0 : fileName, nocheck, maskout, is_meteo, bound_error, nTstepForcingDay)
68 :
69 : use mo_constants, only : nodata_i4
70 : use mo_common_types, only: period
71 : use mo_kind, only : dp, i4
72 : use mo_netcdf, only : NcDataset, NcVariable
73 : use mo_string_utils, only : num2str
74 : use mo_utils, only : eq, ne
75 :
76 : implicit none
77 :
78 : !> Name of the folder where data are stored
79 : character(len = *), intent(in) :: folder
80 : !> Number of datapoints in longitudinal direction
81 : integer(i4), intent(in) :: nRows
82 : !> Number of datapoints in latitudinal direction
83 : integer(i4), intent(in) :: nCols
84 : !> Name of variable name to read
85 : character(len = *), intent(in) :: varName
86 : !> mask of valid data fields
87 : logical, dimension(:, :), intent(in) :: mask
88 : !> Data matrixdim_1 = longitude, dim_2 = latitude, dim_3 = time
89 : real(dp), dimension(:, :, :), allocatable, intent(out) :: data
90 : !> Period the data are needed for
91 : type(period), optional, intent(in) :: target_period
92 : !> Lower bound for check of validity of data values
93 : real(dp), optional, intent(in) :: lower
94 : !> Upper bound for check of validity of data values
95 : real(dp), optional, intent(in) :: upper
96 : !> timestep in netcdf file
97 : integer(i4), optional, intent(in) :: nctimestep
98 : !> name of file, defaults to varName
99 : character(256), optional, intent(in) :: fileName
100 : !> .TRUE. if check for nodata values deactivated, default = .FALSE. - check is done
101 : logical, optional, intent(in) :: nocheck
102 : !> mask of validdata points
103 : logical, dimension(:, :, :), allocatable, optional, intent(out) :: maskout
104 : !> logical whether meteorology is currently read
105 : logical, optional, intent(in) :: is_meteo
106 : !> .FALSE. to only warn about bound (lower, upper) violations, default = .TRUE. - raise an error
107 : logical, optional, intent(in) :: bound_error
108 : !> Number of datapoints in longitudinal direction
109 : integer(i4), optional, intent(inout) :: nTstepForcingDay
110 :
111 : ! netcdf file
112 : type(NcDataset) :: nc
113 : ! variables for data and time form netcdf
114 : type(NcVariable) :: var, time_var
115 : ! shape of NetCDF variable
116 401 : integer(i4), allocatable, dimension(:) :: var_shape
117 : ! index for selecting time vector
118 : integer(i4) :: time_start
119 : ! length of vector of selected time values
120 : integer(i4) :: time_cnt
121 : ! name of NetCDF file
122 : character(256) :: fName
123 : ! loop variable
124 : integer(i4) :: i
125 : ! data nodata value
126 401 : real(dp) :: nodata_value
127 : ! check if model domain is covered by data
128 : logical :: checking
129 : ! check if model domain is covered by data
130 : integer(i4) :: inctimestep
131 : ! raise error or only warn for bound violations
132 : logical :: bound_error_
133 :
134 : ! default value for performing checks on read input
135 401 : checking = .TRUE.
136 401 : if (present(nocheck)) checking = .NOT. nocheck
137 :
138 : ! default value for errors on bound violations
139 401 : bound_error_ = .TRUE.
140 401 : if (present(bound_error)) bound_error_ = bound_error
141 :
142 : ! default: fName = varname + '.nc'
143 : ! optional: fName = filename + '.nc'
144 401 : fName = varName
145 401 : if (present(fileName)) then
146 0 : fName = trim(folder) // trim(fileName) // '.nc'
147 : else
148 401 : fName = trim(folder) // trim(fName) // '.nc'
149 : end if
150 :
151 : ! read the Dataset
152 401 : nc = NcDataset(fname, "r")
153 : ! get the variable
154 401 : var = nc%getVariable(trim(varName))
155 :
156 : ! get dimensions and check if plane is correct
157 1604 : var_shape = var%getShape()
158 401 : if ((var_shape(1) .ne. nRows) .or. (var_shape(2) .ne. nCols)) then
159 0 : call error_message('***ERROR: read_nc: mHM generated x and y are not matching NetCDF dimensions')
160 : end if
161 :
162 : ! determine no data value, use _FillValue first, fall back to missing_value
163 401 : if (var%hasAttribute("_FillValue")) then
164 401 : call var%getAttribute('_FillValue', nodata_value)
165 0 : else if (var%hasAttribute("missing_value")) then
166 0 : call var%getAttribute('missing_value', nodata_value)
167 : else
168 0 : call error_message('***ERROR: read_nc: there must be either the attribute "missing_value" or "_FillValue"')
169 : end if
170 :
171 : ! get time variable
172 401 : time_var = nc%getVariable('time')
173 : ! read the time vector and get start index and count of selection
174 401 : call get_time_vector_and_select(time_var, fname, inctimestep, time_start, time_cnt, target_period)
175 :
176 401 : if (present(is_meteo)) then
177 391 : if (.not. present(nTstepForcingDay)) call error_message("read_nc: meteo file given but 'nTstepForcingDay' not passed")
178 391 : if (is_meteo) then
179 391 : select case(inctimestep)
180 : case(-1) ! daily
181 391 : if (nTstepForcingDay .eq. nodata_i4) then
182 14 : nTstepForcingDay = 1_i4
183 377 : else if (nTstepForcingDay .ne. 1_i4) then
184 : call error_message('***ERROR: read_forcing_nc: expected daily input forcing, but read something else. ' // &
185 0 : 'Ensure all input time steps have the same units across all input files')
186 : end if
187 : case(-4) ! hourly
188 0 : if (nTstepForcingDay .eq. nodata_i4) then
189 0 : nTstepForcingDay = 24_i4
190 0 : else if (nTstepForcingDay .ne. 24_i4) then
191 : call error_message('***ERROR: read_forcing_nc: expected hourly input forcing, but read something else. ' // &
192 0 : 'Ensure all input time steps have the same units across all input files')
193 : end if
194 : case default ! no output at all
195 391 : call error_message('***ERROR: read_nc: unknown nctimestep switch.')
196 : end select
197 : end if
198 : end if
199 :
200 : ! check optional nctimestep
201 401 : if (present(nctimestep)) then
202 10 : if (inctimestep .ne. nctimestep) then
203 : call error_message('***ERROR: provided timestep ' // num2str(nctimestep) //&
204 0 : ' does not match with the one in file ' // num2str(inctimestep), raise=.false.)
205 0 : call error_message('File: ' // trim(fname))
206 : end if
207 : end if
208 :
209 : ! extract data and select time slice
210 2807 : call var%getData(data, start = (/1, 1, time_start/), cnt = (/nRows, nCols, time_cnt/))
211 :
212 : ! save output mask if optional maskout is given
213 401 : if (present(maskout)) then
214 45 : allocate(maskout(var_shape(1), var_shape(2), var_shape(3)))
215 106529 : maskout = ne(data(:, :, :), nodata_value)
216 : end if
217 :
218 : ! start checking values
219 64159 : do i = 1, size(data, dim = 3)
220 : ! neglect checking for nodata values if optional nocheck is given
221 63758 : if (checking) then
222 6192456 : if (any(eq(data(:, :, i), nodata_value) .and. (mask))) then
223 0 : call error_message('***ERROR: read_nc: nodata value within domain ', raise=.false.)
224 0 : call error_message(' boundary in variable: ', trim(varName), raise=.false.)
225 0 : call error_message(' at timestep : ', trim(num2str(i)))
226 : end if
227 : end if
228 : ! optional check
229 63758 : if (present(lower)) then
230 6192456 : if (any((data(:, :, i) .lt. lower) .AND. mask(:, :))) then
231 : call error_message('***ERROR: read_nc: values in variable "', &
232 0 : trim(varName), '" are lower than ', trim(num2str(lower, '(F7.2)')), raise=.false.)
233 0 : call error_message(' at timestep : ', trim(num2str(i)), raise=.false.)
234 0 : call error_message('File: ', trim(fName), raise=.false.)
235 0 : call error_message('Minval at timestep: ', trim(num2str(minval(data(:, :, i)))), raise=.false.)
236 0 : call error_message('Total minval: ', trim(num2str(minval(data(:, :, :)))), raise=bound_error_)
237 : end if
238 : end if
239 :
240 64159 : if (present(upper)) then
241 6192456 : if (any((data(:, :, i) .gt. upper) .AND. mask(:, :))) then
242 : call error_message('***ERROR: read_nc: values in variable "', &
243 0 : trim(varName), '" are greater than ', trim(num2str(upper, '(F7.2)')), raise=.false.)
244 0 : call error_message(' at timestep : ', trim(num2str(i)), raise=.false.)
245 0 : call error_message('File: ', trim(fName), raise=.false.)
246 0 : call error_message('Maxval at timestep: ', trim(num2str(maxval(data(:, :, i)))), raise=.false.)
247 0 : call error_message('Total maxval: ', trim(num2str(maxval(data(:, :, :)))), raise=bound_error_)
248 : end if
249 : end if
250 :
251 : end do
252 :
253 802 : end subroutine read_nc
254 :
255 :
256 : !> \brief Reads time independent forcing input in NetCDF file format.
257 : !> \details Reads time independent netCDF forcing files. \n
258 : !! First, the dimensions given are cross-checked with header.txt information. Second, the data of the
259 : !! specified period are read from the specified directory.
260 : !! If the optional lower and/or upper bound for the data values is given, the read data are checked for validity.
261 : !! The program is stopped if any value lies out of range.\n
262 : !! If the optinal argument nocheck is true, the data are not checked for coverage with the input mask.
263 : !! Additionally in this case an mask of vild data points can be received from the routine in maskout.
264 : !!
265 : !! \note Files have to be called like defined in mo_files. Furthermore the variable names have to be called
266 : !! like they are defined in the declaration of this subroutine. The NetCDF file has to have 2 dimensions:
267 : !! 1. x, 2. y, It is expected that the variables (especially)within the NetCDF files contain an
268 : !! unit attribute. The timestep has to be equidistant.
269 : !!
270 : !> \author Lennart Schueler, heavily influenced by read_nc
271 : !> \date May 2018
272 1 : subroutine read_const_nc(folder, nRows, nCols, varName, data, fileName)
273 :
274 802 : use mo_kind, only: i4, dp
275 : use mo_netcdf, only: NcDataset, NcVariable
276 :
277 : implicit none
278 :
279 : character(len=*), intent(in) :: folder !< folder where data are stored
280 : integer(i4), intent(in) :: nRows !< number of rows of data fields:
281 : integer(i4), intent(in) :: nCols !< number of columns of data fields:
282 : character(len=*), intent(in) :: varName !< name of NetCDF variable
283 : real(dp), dimension(:,:), allocatable, intent(out) :: data !< data read in
284 : character(256), optional, intent(in) :: fileName !< name of file, defaults to varName
285 :
286 : ! local variables
287 : type(NcDataset) :: nc ! netcdf file
288 : type(NcVariable) :: var ! variables for data form netcdf
289 1 : integer(i4), allocatable, dimension(:) :: var_shape ! shape of NetCDF variable
290 :
291 : character(256) :: fName ! name of NetCDF file
292 1 : real(dp) :: nodata_value ! data nodata value
293 :
294 1 : fName = varName
295 1 : if (present(fileName)) then
296 1 : fName = trim(folder) // trim(fileName) // '.nc'
297 : else
298 0 : fName = trim(folder) // trim(fName) // '.nc'
299 : end if
300 :
301 : ! read the Dataset
302 1 : nc = NcDataset(fname, "r")
303 : ! get the variable
304 1 : var = nc%getVariable(trim(varName))
305 :
306 : ! get dimensions and check if plane is correct
307 3 : var_shape = var%getShape()
308 1 : if ( (var_shape(1) .ne. nRows) .or. (var_shape(2) .ne. nCols) ) then
309 0 : call error_message('***ERROR: read_const_nc: mHM generated x and y are not matching NetCDF dimensions')
310 : end if
311 :
312 : ! determine no data value, use _FillValue first, fall back to missing_value
313 1 : if (var%hasAttribute("_FillValue")) then
314 0 : call var%getAttribute('_FillValue', nodata_value)
315 1 : else if (var%hasAttribute("missing_value")) then
316 1 : call var%getAttribute('missing_value', nodata_value)
317 : else
318 0 : call error_message('***ERROR: read_const_nc: there must be either the attribute "missing_value" or "_FillValue"')
319 : end if
320 :
321 : ! extract data and select time slice
322 3 : call var%getData(data, start=(/1,1/), cnt=(/nRows,nCols/))
323 :
324 1 : end subroutine read_const_nc
325 :
326 :
327 : !> \brief Reads weights for meteo forcings input in NetCDF file format.
328 : !> \details Reads netCDF weight files.
329 : !! First, the dimensions given are cross-checked with header.txt information. If the optional lower
330 : !! and/or upper bound for the data values is given, the read data are checked for validity.
331 : !! The program is stopped if any value lies out of range.
332 : !! If the optinal argument nocheck is true, the data are not checked for coverage with the input mask.
333 : !! Additionally in this case an mask of vild data points can be received from the routine in maskout.
334 : !> \changelog
335 : !! - Robert Schweppe Nov 2017
336 : !! - switched to mo_netcdf library and restuctured routine
337 : !! - Robert Schweppe Jun 2018
338 : !! - refactoring and reformatting
339 : !! - Sebastian Müller Mar 2023
340 : !! - added bound_error to control if an error is raised if boundaries are violated
341 : !> \authors Stephan Thober & Matthias Zink
342 : !> \date Jan 2017
343 0 : subroutine read_weights_nc(folder, nRows, nCols, varName, data, mask, lower, upper, nocheck, maskout, fileName, bound_error)
344 :
345 1 : use mo_kind, only : dp, i4
346 : use mo_netcdf, only : NcDataset, NcVariable
347 : use mo_string_utils, only : num2str
348 : use mo_utils, only : eq, ne
349 :
350 : implicit none
351 :
352 : !> Name of the folder where data are stored
353 : character(len = *), intent(in) :: folder
354 : !> Number of datapoints in longitudinal direction
355 : integer(i4), intent(in) :: nRows
356 : !> Number of datapoints in latitudinal direction
357 : integer(i4), intent(in) :: nCols
358 : !> Name of variable name to read
359 : character(len = *), intent(in) :: varName
360 : !> Data matrixdim_1 = longitude, dim_2 = latitude, dim_3 = months, dim_4 = hours
361 : real(dp), dimension(:, :, :, :), allocatable, intent(out) :: data
362 : !> mask of valid data fields
363 : logical, dimension(:, :), intent(in) :: mask
364 : !> Lower bound for check of validity of data values
365 : real(dp), optional, intent(in) :: lower
366 : !> Upper bound for check of validity of data values
367 : real(dp), optional, intent(in) :: upper
368 : !> .TRUE. if check for nodata values deactivateddefault = .FALSE. - check is done
369 : logical, optional, intent(in) :: nocheck
370 : !> name of variable, defaults to fileName
371 : character(256), optional, intent(in) :: fileName
372 : !> ! mask of validdata points
373 : logical, dimension(:, :, :, :), allocatable, optional, intent(out) :: maskout
374 : !> .FALSE. to only warn about bound (lower, upper) violations, default = .TRUE. - raise an error
375 : logical, optional, intent(in) :: bound_error
376 :
377 : ! name of NetCDF file
378 : character(256) :: fName
379 : ! loop variable
380 : integer(i4) :: i
381 : ! loop variable
382 : integer(i4) :: j
383 : ! data nodata value
384 0 : real(dp) :: nodata_value
385 : ! check if model domain is covered by data
386 : logical :: checking
387 : ! container for Netcdf data
388 : type(NcDataset) :: nc
389 : ! container for Netcdf variable
390 : type(NcVariable) :: var
391 : ! shape of NetCDF variable
392 0 : integer(i4), allocatable, dimension(:) :: var_shape
393 : ! raise error or only warn for bound violations
394 : logical :: bound_error_
395 :
396 0 : checking = .TRUE.
397 0 : if (present(nocheck)) checking = .NOT. nocheck
398 :
399 : ! default value for errors on bound violations
400 0 : bound_error_ = .TRUE.
401 0 : if (present(bound_error)) bound_error_ = bound_error
402 :
403 0 : fName = varName
404 0 : if (present(fileName)) then
405 0 : fName = trim(folder) // fileName
406 : else
407 0 : fName = trim(folder) // trim(fName) // '.nc'
408 : end if
409 :
410 0 : nc = NcDataset(fname, "r")
411 0 : var = nc%getVariable(trim(varName))
412 :
413 : ! get dimensions
414 0 : var_shape = var%getShape()
415 0 : if ((var_shape(1) .ne. nRows) .or. (var_shape(2) .ne. nCols)) then
416 0 : call error_message('***ERROR: read_weights_nc: mHM generated x and y are not matching NetCDF dimensions')
417 : end if
418 :
419 : ! determine no data value
420 0 : call var%getAttribute('missing_value', nodata_value)
421 :
422 : ! extract data
423 0 : call var%getData(data)
424 :
425 : ! save output mask if optional maskout is given
426 0 : if (present(maskout)) then
427 0 : allocate(maskout(var_shape(1), var_shape(2), var_shape(3), var_shape(4)))
428 0 : maskout = ne(data(:, :, :, :), nodata_value)
429 : end if
430 :
431 : ! start checking values
432 0 : do i = 1, var_shape(3)
433 0 : do j = 1, var_shape(4)
434 : ! neglect checking for naodata values if optional nocheck is given
435 0 : if (checking) then
436 0 : if (any(eq(data(:, :, i, j), nodata_value) .and. (mask))) then
437 0 : call error_message('***ERROR: read_weights_nc: nodata value within domain ', raise=.false.)
438 0 : call error_message(' boundary in variable: ', trim(varName), raise=.false.)
439 0 : call error_message(' at hour : ', trim(num2str(i)))
440 : end if
441 : end if
442 : ! optional check
443 0 : if (present(lower)) then
444 0 : if (any((data(:, :, i, j) .lt. lower) .AND. mask(:, :))) then
445 : call error_message('***ERROR: read_weights_nc: values in variable "', &
446 0 : trim(varName), '" are lower than ', trim(num2str(lower, '(F7.2)')), raise=.false.)
447 0 : call error_message(' at hour : ', trim(num2str(i)), raise=.false.)
448 0 : call error_message('File: ', trim(fName), raise=.false.)
449 0 : call error_message('Minval at hour: ', trim(num2str(minval(data(:, :, i, j)), '(F7.2)')), raise=.false.)
450 0 : call error_message('Total minval: ', trim(num2str(minval(data(:, :, :, :)), '(F7.2)')), raise=bound_error_)
451 : end if
452 : end if
453 :
454 0 : if (present(upper)) then
455 0 : if (any((data(:, :, i, j) .gt. upper) .AND. mask(:, :))) then
456 : call error_message('***ERROR: read_weights_nc: values in variable "', &
457 0 : trim(varName), '" are greater than ', trim(num2str(upper, '(F7.2)')), raise=.false.)
458 0 : call error_message(' at hour : ', trim(num2str(i)), raise=.false.)
459 0 : call error_message('File: ', trim(fName), raise=.false.)
460 0 : call error_message('Maxval at hour: ', trim(num2str(maxval(data(:, :, i, j)), '(F7.2)')), raise=.false.)
461 0 : call error_message('Total maxval: ', trim(num2str(maxval(data(:, :, :, :)), '(F7.2)')), raise=bound_error_)
462 : end if
463 : end if
464 :
465 : end do
466 : end do
467 :
468 0 : end subroutine read_weights_nc
469 :
470 :
471 : !> \brief Extract time vector in unit julian hours and get supposed time step in hours
472 : !> \changelog
473 : !! - Matthias Cuntz & Juliane Mai Nov 2014
474 : !! - time int or double
475 : !! - Stephan Thober Sep 2015
476 : !! - added read for hourly data
477 : !! - Robert Schweppe Nov 2017
478 : !! - restructured routine, reads vector now
479 : !! - Maren Kaluza May 2018
480 : !! - fixed bug in time reading
481 : !! - Stephan Thober Aug 2020
482 : !! - fixed hourly reading
483 : !! - Stephan Thober Jan 2022
484 : !! - deactivated monthly and annual reading added nTstepForcingDay for hourly reading
485 : !> \authors Matthias Zink
486 : !> \date Oct 2012
487 401 : subroutine get_time_vector_and_select(var, fname, inctimestep, time_start, time_cnt, target_period)
488 :
489 0 : use mo_common_types, only: period
490 : use mo_constants, only : DayHours, DaySecs, YearDays
491 : use mo_julian, only : caldat, dec2date, julday
492 : use mo_kind, only : dp, i4, i8
493 : use mo_netcdf, only : NcVariable
494 : use mo_string_utils, only : DIVIDE_STRING
495 :
496 : implicit none
497 :
498 : !> variable of interest
499 : type(NcVariable), intent(in) :: var
500 : !> fname of ncfile for error message
501 : character(256), intent(in) :: fname
502 : !> flag for requested time step
503 : integer(i4), intent(out) :: inctimestep
504 : !> time_start index of time selection
505 : integer(i4), intent(out) :: time_start
506 : !> time_count of indexes of time selection
507 : integer(i4), intent(out) :: time_cnt
508 : !> reference period
509 : type(period), intent(in), optional :: target_period
510 :
511 : ! reference time of NetCDF
512 : integer(i4) :: yRef, dRef, mRef, hRef, jRef
513 : ! netcdf attribute values
514 : character(256) :: AttValues
515 : ! dummies for netcdf attribute handling
516 401 : character(256), dimension(:), allocatable :: strArr, date, time
517 : ! native time step converter in ncfile
518 : integer(i8) :: time_step_seconds
519 : ! time vector
520 401 : integer(i8), allocatable, dimension(:) :: time_data
521 401 : integer(i8), allocatable, dimension(:) :: time_diff
522 : ! period of ncfile, for clipping
523 : type(period) :: nc_period, clip_period
524 : integer(i4) :: ncJulSta1, dd, n_time
525 : integer(i4) :: mmcalstart, mmcalend, yycalstart, yycalend
526 : integer(i4) :: mmncstart, yyncstart
527 : ! helper variable for error output
528 : integer(i4) :: hstart_int, hend_int
529 :
530 401 : call var%getAttribute('units', AttValues)
531 : ! AttValues looks like "<unit> since YYYY-MM-DD[ HH:MM:SS]"
532 : ! split at space
533 401 : call DIVIDE_STRING(trim(AttValues), ' ', strArr)
534 :
535 : ! determine reference time at '-' and convert to integer
536 401 : call DIVIDE_STRING(trim(strArr(3)), '-', date)
537 401 : read(date(1), *) yRef
538 401 : read(date(2), *) mRef
539 401 : read(date(3), *) dRef
540 :
541 401 : jRef = julday(dd = dRef, mm = mRef, yy = yRef)
542 :
543 : ! if existing also read in the time (only hour so far)
544 401 : hRef = 0
545 401 : if(size(strArr) .gt. 3) then
546 347 : call DIVIDE_STRING(trim(strArr(4)), ':', time)
547 347 : read(time(1), *) hRef
548 : end if
549 :
550 : ! determine the step_size
551 401 : if (strArr(1) .EQ. 'days') then
552 : time_step_seconds = int(DaySecs)
553 8 : else if (strArr(1) .eq. 'hours') then
554 : time_step_seconds = int(DaySecs / DayHours)
555 0 : else if (strArr(1) .eq. 'minutes') then
556 : time_step_seconds = int(DaySecs / DayHours / 60._dp)
557 0 : else if (strArr(1) .eq. 'seconds') then
558 : time_step_seconds = 1_i8
559 : else
560 : call error_message('***ERROR: Please provide the input data in (days, hours, minutes, seconds) ', &
561 0 : 'since YYYY-MM-DD[ HH:MM:SS] in the netcdf file. Found: ', trim(AttValues))
562 : end if
563 :
564 : ! get the time vector
565 401 : call var%getData(time_data)
566 : ! convert array from units since to seconds
567 841186 : time_data = time_data * time_step_seconds
568 :
569 : ! check for length of time vector, needs to be at least of length 2, otherwise step width check fails
570 401 : if (size(time_data) .le. 1) then
571 0 : call error_message('***ERROR: length of time dimension needs to be at least 2 in file: ' // trim(fname))
572 : end if
573 :
574 : ! compare the read period from ncfile to the period required
575 : ! convert julian second information back to date via conversion to float
576 : ! the 0.5_dp is for the different reference of fractional julian days, hours are truncated
577 401 : n_time = size(time_data)
578 401 : call dec2date(time_data(1) / DaySecs - 0.5_dp + jRef + hRef / 24._dp, nc_period%dStart, nc_period%mStart, &
579 401 : nc_period%yStart, hstart_int)
580 401 : nc_period%julStart = int(time_data(1) / DaySecs + jRef + hRef / 24._dp)
581 401 : call dec2date(time_data(n_time) / DaySecs - 0.5_dp + jRef + hRef / 24._dp, nc_period%dEnd, nc_period%mEnd, &
582 802 : nc_period%yEnd, hend_int)
583 401 : nc_period%julEnd = int(time_data(n_time) / DaySecs + jRef + hRef / 24._dp)
584 :
585 : ! if no target period is present, use the whole time period
586 401 : if (present(target_period)) then
587 401 : clip_period = target_period
588 : else
589 0 : clip_period = nc_period
590 : end if
591 :
592 : ! calculate input resolution
593 1203 : allocate(time_diff(n_time - 1))
594 841186 : time_diff = (time_data(2 : n_time) - time_data(1 : n_time - 1)) / int(DaySecs, i8)
595 : ! difference must be 1 day
596 840350 : if (all(abs(time_diff - 1._dp) .lt. 1._dp)) then
597 392 : inctimestep = -1 ! daily
598 : ! difference must be between 28 and 31 days
599 888 : else if (all(abs(time_diff) .lt. 32._dp) .and. all(abs(time_diff) .gt. 27._dp)) then
600 9 : inctimestep = -2 ! monthly
601 : ! difference must be between 365 and 366 days
602 0 : else if ((all(abs(time_diff) .lt. YearDays + 2)) .and. (all(abs(time_diff) .gt. YearDays - 1._dp))) then
603 0 : inctimestep = -3 ! yearly
604 : ! difference must be 1 hour
605 0 : else if (all(abs((time_data(2 : n_time) - time_data(1 : n_time - 1)) / 3600._dp - 1._dp) .lt. 1.e-6)) then
606 0 : inctimestep = -4 ! hourly
607 : else
608 : call error_message('***ERROR: time_steps are not equal over all times in file and/or do not conform to' // &
609 0 : ' requested timestep in file (', trim(fname), ')')
610 : end if
611 :
612 : ! prepare the selection and check for required time_step
613 793 : select case(inctimestep)
614 : case(-1) ! daily
615 392 : ncJulSta1 = nc_period%julStart
616 392 : time_start = clip_period%julStart - ncJulSta1 + 1_i4
617 392 : time_cnt = clip_period%julEnd - clip_period%julStart + 1_i4
618 : case(-2) ! monthly
619 9 : call caldat(clip_period%julStart, dd, mmcalstart, yycalstart)
620 9 : call caldat(nc_period%julStart, dd, mmncstart, yyncstart)
621 : ! monthly timesteps are usually set by month end, so for beginning, we need 1st of month
622 9 : ncJulSta1 = julday(1, mmncstart, yyncstart)
623 9 : call caldat(clip_period%julEnd, dd, mmcalend, yycalend)
624 9 : time_start = (yycalstart * 12 + mmcalstart) - (yyncstart * 12 + mmncstart) + 1_i4
625 9 : time_cnt = (yycalend * 12 + mmcalend) - (yycalstart * 12 + mmcalstart) + 1_i4
626 : case(-3) ! yearly
627 0 : call caldat(clip_period%julStart, dd, mmcalstart, yycalstart)
628 0 : call caldat(nc_period%julStart, dd, mmncstart, yyncstart)
629 : ! yearly timesteps are usually set by year end, so for beginning, we need 1st of year
630 0 : ncJulSta1 = julday(1, 1, yyncstart)
631 0 : call caldat(clip_period%julEnd, dd, mmcalend, yycalend)
632 0 : time_start = yycalstart - yyncstart + 1_i4
633 0 : time_cnt = yycalend - yycalstart + 1_i4
634 : case(-4) ! hourly
635 0 : ncJulSta1 = nc_period%julStart
636 0 : time_start = (clip_period%julStart - ncJulSta1) * 24_i4 + 1_i4 ! convert to hours; always starts at one
637 0 : time_cnt = (clip_period%julEnd - clip_period%julStart + 1_i4) * 24_i4 ! convert to hours
638 : case default ! no output at all
639 401 : call error_message('***ERROR: read_nc: unknown nctimestep switch.')
640 : end select
641 :
642 : ! Check if time steps in file cover simulation period
643 401 : if (.not. ((ncJulSta1 .LE. clip_period%julStart) .AND. (nc_period%julEnd .GE. clip_period%julEnd))) then
644 : call error_message('***ERROR: read_nc: time period of input data: ', trim(fname), &
645 0 : ' is not matching modelling period.')
646 : end if
647 :
648 : ! free memory
649 401 : deallocate(time_diff)
650 :
651 401 : end subroutine get_time_vector_and_select
652 :
653 : end module mo_read_nc
|
