Line data Source code
1 : !> \file mo_meteo_helper.f90
2 : !> \brief \copybrief mo_meteo_helper
3 : !> \details \copydetails mo_meteo_helper
4 :
5 : !> \brief Prepare meteorological forcings data for mHM.
6 : !> \details Prepare meteorological forcings data for mHM.
7 : !> \authors Rohini Kumar
8 : !> \date Jan 2012
9 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
10 : !! mHM is released under the LGPLv3+ license \license_note
11 : !> \ingroup f_meteo
12 : MODULE mo_meteo_helper
13 :
14 : ! This module provides routines to read meteorological data.
15 :
16 : ! Written Rohini Kumar, Jan 2013
17 : ! Modified
18 :
19 : USE mo_kind, ONLY : i4, dp
20 : use mo_message, only: message, error_message
21 :
22 : IMPLICIT NONE
23 :
24 : PRIVATE
25 :
26 : public :: meteo_forcings_wrapper
27 : public :: meteo_weights_wrapper
28 : public :: chunk_config
29 :
30 : CONTAINS
31 :
32 : !> \brief Prepare meteorological forcings data for mHM at Level-1
33 : !> \details Prepare meteorological forcings data for mHM, which include
34 : !! 1) Reading meteo. datasets at their native resolution for every Domain
35 : !! 2) Perform aggregation or disaggregation of meteo. datasets from their
36 : !! native resolution (level-2) to the required hydrologic resolution (level-1)
37 : !! 3) Pad the above datasets of every Domain to their respective global ones
38 : !> \changelog
39 : !! - Stephan Thober Jun 2014
40 : !! - changed to readPer
41 : !! - Stephan Thober Feb 2016
42 : !! - refactored deallocate statements
43 : !! - Robert Schweppe Jun 2018
44 : !! - refactoring and reformatting
45 : !! - Sebastian Müller Mar 2023
46 : !! - used by meteo-handler
47 : !! - now independent of global variables
48 : !> \authors Rohini Kumar
49 : !> \date Jan 2013
50 391 : subroutine meteo_forcings_wrapper( &
51 391 : iDomain, dataPath, inputFormat, dataOut1, readPer, nTstepForcingDay, level1, level2, lower, upper, ncvarName, bound_error)
52 :
53 : use mo_append, only : append
54 : use mo_common_constants, only : nodata_dp
55 : use mo_common_types, only : period, grid
56 : use mo_read_nc, only : read_nc
57 : use mo_meteo_spatial_tools, only : spatial_aggregation, spatial_disaggregation
58 :
59 : implicit none
60 :
61 : !> Domain Id
62 : integer(i4), intent(in) :: iDomain
63 : !> Data path where a given meteo. variable is stored
64 : character(len = *), intent(in) :: dataPath
65 : !> only 'nc' possible at the moment
66 : character(len = *), intent(in) :: inputFormat
67 : !> Packed meterological variable for the whole simulation period
68 : real(dp), dimension(:, :), allocatable, intent(inout) :: dataOut1
69 : !> start and end dates of read period
70 : type(period), intent(in) :: readPer
71 : !> Number of forcing intervals per day
72 : integer(i4), intent(inout) :: nTstepForcingDay
73 : !> grid information at hydrologic level
74 : type(Grid), dimension(:), intent(in) :: level1
75 : !> Reference of the metereological variables
76 : type(Grid), dimension(:), intent(in) :: level2
77 : !> Lower bound for check of validity of data values
78 : real(dp), optional, intent(in) :: lower
79 : !> Upper bound for check of validity of data values
80 : real(dp), optional, intent(in) :: upper
81 : !> name of the variable (for .nc files)
82 : character(len = *), optional, intent(in) :: ncvarName
83 : !> .FALSE. to only warn about bound (lower, upper) violations, default = .TRUE. - raise an error
84 : logical, optional, intent(in) :: bound_error
85 :
86 391 : logical, dimension(:, :), allocatable :: mask1
87 : integer(i4) :: ncells1
88 : integer(i4) :: nrows2, ncols2
89 391 : logical, dimension(:, :), allocatable :: mask2
90 : ! meteo data at level-2
91 391 : real(dp), dimension(:, :, :), allocatable :: L2_data
92 : ! meteo data at level-1
93 391 : real(dp), dimension(:, :, :), allocatable :: L1_data
94 : ! packed meteo data at level-1 from 3D to 2D
95 391 : real(dp), dimension(:, :), allocatable :: L1_data_packed
96 : integer(i4) :: nTimeSteps
97 : ! level-1_resolution/level-2_resolution
98 391 : real(dp) :: cellFactorHbyM
99 : integer(i4) :: t
100 :
101 : ! get basic Domain information at level-1
102 391 : nCells1 = level1(iDomain)%nCells
103 27147 : mask1 = level1(iDomain)%mask
104 :
105 : ! make basic Domain information at level-2
106 391 : nrows2 = level2(iDomain)%nrows
107 391 : ncols2 = level2(iDomain)%ncols
108 25999 : mask2 = level2(iDomain)%mask
109 :
110 :
111 782 : select case (trim(inputFormat))
112 : case('nc')
113 : ! Fortran 2008 allowes to pass optional dummy arguments to another routine (when they are optional there too)
114 : CALL read_nc(dataPath, nRows2, nCols2, ncvarName, mask2, L2_data, target_period=readPer, &
115 391 : lower=lower, upper=upper, is_meteo=.True., bound_error=bound_error, nTstepForcingDay=nTstepForcingDay)
116 : case DEFAULT
117 782 : call error_message('***ERROR: meteo_forcings_wrapper: Not recognized input format')
118 : end select
119 : ! cellfactor to decide on the upscaling or downscaling of meteo. fields
120 391 : cellFactorHbyM = level1(iDomain)%cellsize / level2(iDomain)%cellsize
121 :
122 : ! upscaling & packing
123 391 : if(cellFactorHbyM .gt. 1.0_dp) then
124 0 : call spatial_aggregation(L2_data, level2(iDomain)%cellsize, level1(iDomain)%cellsize, mask1, mask2, L1_data)
125 : ! downscaling
126 400 : elseif(cellFactorHbyM .lt. 1.0_dp) then
127 9 : call spatial_disaggregation(L2_data, level2(iDomain)%cellsize, level1(iDomain)%cellsize, mask1, mask2, L1_data)
128 : ! nothing
129 : else
130 1910 : allocate(L1_data(size(L2_data, 1), size(L2_data, 2), size(L2_data, 3)))
131 3585940 : L1_data(:, :, :) = L2_data(:, :, :)
132 : end if
133 : ! free memory immediately
134 391 : deallocate(L2_data)
135 :
136 : ! pack variables
137 391 : nTimeSteps = size(L1_data, 3)
138 1564 : allocate(L1_data_packed(nCells1, nTimeSteps))
139 63442 : do t = 1, nTimeSteps
140 :
141 63442 : L1_data_packed(:, t) = pack(L1_data(:, :, t), MASK = mask1(:, :))
142 :
143 : end do
144 :
145 : ! free memory immediately
146 391 : deallocate(L1_data)
147 :
148 : ! append
149 391 : call append(dataOut1, L1_data_packed(:, :), nodata_dp)
150 :
151 : !free space
152 391 : deallocate(L1_data_packed)
153 :
154 391 : end subroutine meteo_forcings_wrapper
155 :
156 :
157 : !> \brief Prepare weights for meteorological forcings data for mHM at Level-1
158 : !> \details Prepare meteorological weights data for mHM, which include
159 : !! 1) Reading meteo. weights datasets at their native resolution for every Domain
160 : !! 2) Perform aggregation or disaggregation of meteo. weights datasets from their
161 : !! native resolution (level-2) to the required hydrologic resolution (level-1)
162 : !! 3) Pad the above datasets of every Domain to their respective global ones
163 : !> \changelog
164 : !! - Stephan Thober May 2017
165 : !! - updated documentation
166 : !! - Robert Schweppe Jun 2018
167 : !! - refactoring and reformatting
168 : !! - Oldrich Rakovec Aug 2018
169 : !! - adding message about reading meteo_weights
170 : !! - Sebastian Müller Mar 2023
171 : !! - used by meteo-handler
172 : !! - now independent of global variables
173 : !> \authors Stephan Thober & Rohini Kumar
174 : !> \date Jan 2017
175 156 : subroutine meteo_weights_wrapper(iDomain, read_meteo_weights, dataPath, dataOut1, level1, level2, lower, upper, ncvarName)
176 :
177 391 : use mo_append, only : append
178 : use mo_common_constants, only : nodata_dp
179 : use mo_common_types, only : grid
180 : use mo_read_nc, only : read_weights_nc
181 : use mo_meteo_spatial_tools, only : spatial_aggregation, spatial_disaggregation
182 :
183 : implicit none
184 :
185 : !> Domain Id
186 : integer(i4), intent(in) :: iDomain
187 : !> Flag for reading meteo weights
188 : logical, intent(in) :: read_meteo_weights
189 : !> Data path where a given meteo. variable is stored
190 : character(len = *), intent(in) :: dataPath
191 : !> Packed meterological variable for the whole simulation period
192 : real(dp), dimension(:, :, :), allocatable, intent(inout) :: dataOut1
193 : !> grid information at hydrologic level
194 : type(Grid), dimension(:), intent(in) :: level1
195 : !> Reference of the metereological variables
196 : type(Grid), dimension(:), intent(in) :: level2
197 : !> Lower bound for check of validity of data values
198 : real(dp), optional, intent(in) :: lower
199 : !> Upper bound for check of validity of data values
200 : real(dp), optional, intent(in) :: upper
201 : !> name of the variable (for .nc files)
202 : character(len = *), optional, intent(in) :: ncvarName
203 :
204 78 : logical, dimension(:, :), allocatable :: mask1
205 : integer(i4) :: ncells1
206 : integer(i4) :: nrows2, ncols2
207 78 : logical, dimension(:, :), allocatable :: mask2
208 : ! meteo weights data at level-2
209 78 : real(dp), dimension(:, :, :, :), allocatable :: L2_data
210 : ! meteo weights data at level-1
211 78 : real(dp), dimension(:, :, :, :), allocatable :: L1_data
212 : ! packed meteo weights data at level-1 from 4D to 3D
213 78 : real(dp), dimension(:, :, :), allocatable :: L1_data_packed
214 : integer(i4) :: nMonths, nHours
215 : ! level-1_resolution/level-2_resolution
216 78 : real(dp) :: cellFactorHbyM
217 : integer(i4) :: t, j
218 :
219 : ! get basic Domain information at level-1
220 78 : nCells1 = level1(iDomain)%nCells
221 6489 : mask1 = level1(iDomain)%mask
222 :
223 : ! make basic Domain information at level-2
224 78 : nrows2 = level2(iDomain)%nrows
225 78 : ncols2 = level2(iDomain)%ncols
226 5028 : mask2 = level2(iDomain)%mask
227 :
228 78 : if (read_meteo_weights) then
229 0 : call message(' read_meteo_weights = .TRUE. ... Reading meteo weights ... ')
230 0 : CALL read_weights_nc(dataPath, nRows2, nCols2, ncvarName, L2_data, mask2, lower=lower, upper=upper)
231 :
232 : ! cellfactor to decide on the upscaling or downscaling of meteo. fields
233 0 : cellFactorHbyM = level1(iDomain)%cellsize / level2(iDomain)%cellsize
234 :
235 : ! upscaling & packing
236 0 : if(cellFactorHbyM .gt. 1.0_dp) then
237 0 : call spatial_aggregation(L2_data, level2(iDomain)%cellsize, level1(iDomain)%cellsize, mask1, mask2, L1_data)
238 : ! downscaling
239 0 : elseif(cellFactorHbyM .lt. 1.0_dp) then
240 0 : call spatial_disaggregation(L2_data, level2(iDomain)%cellsize, level1(iDomain)%cellsize, mask1, mask2, L1_data)
241 : ! nothing
242 : else
243 0 : L1_data = L2_data
244 : end if
245 : ! free memory immediately
246 0 : deallocate(L2_data)
247 :
248 : ! pack variables
249 0 : nMonths = size(L1_data, dim = 3)
250 0 : nHours = size(L1_data, dim = 4)
251 0 : allocate(L1_data_packed(nCells1, nMonths, nHours))
252 0 : do t = 1, nMonths
253 0 : do j = 1, nHours
254 0 : L1_data_packed(:, t, j) = pack(L1_data(:, :, t, j), MASK = mask1(:, :))
255 : end do
256 : end do
257 : ! free memory immediately
258 0 : deallocate(L1_data)
259 : else
260 : ! dummy allocation
261 312 : allocate(L1_data_packed(nCells1, 12, 24))
262 956748 : L1_data_packed = nodata_dp
263 : end if
264 :
265 : ! append
266 78 : call append(dataOut1, L1_data_packed)
267 :
268 : !free space
269 78 : deallocate(L1_data_packed)
270 :
271 78 : end subroutine meteo_weights_wrapper
272 :
273 :
274 : !> \brief determines the start date, end date, and read_flag given Domain id and current timestep
275 : !> \changelog
276 : !! - Robert Schweppe Jun 2018
277 : !! - refactoring and reformatting
278 : !! - Sebastian Müller Mar 2023
279 : !! - used by meteo-handler
280 : !! - now independent of global variables
281 : !> \authors Stephan Thober
282 : !> \date Jun 2014
283 131082 : subroutine chunk_config(iDomain, tt, nTstepDay, simPer, timestep, timeStep_model_inputs, read_flag, readPer)
284 :
285 78 : use mo_common_constants, only : nodata_dp
286 : use mo_common_types, only: period
287 :
288 : implicit none
289 :
290 : !> current Domain
291 : integer(i4), intent(in) :: iDomain
292 : !> current timestep
293 : integer(i4), intent(in) :: tt
294 : !> Number of time intervals per day
295 : integer(i4), intent(in) :: nTstepDay
296 : !> warmPer + evalPer
297 : type(period), dimension(:), intent(in) :: simPer
298 : !> [h] simulation time step (= TS) in [h]
299 : integer(i4), intent(in) :: timeStep
300 : !> frequency for reading meteo input
301 : integer(i4), dimension(:), intent(in) :: timeStep_model_inputs
302 : !> indicate whether reading data should be read
303 : logical, intent(out) :: read_flag
304 : ! start and end dates of reading Period
305 : type(period), intent(out) :: readPer
306 :
307 : ! initialize
308 65541 : read_flag = .false.
309 65541 : if (tt .eq. 1_i4) then
310 26 : readPer%julStart = int(nodata_dp, i4)
311 26 : readPer%julEnd = int(nodata_dp, i4)
312 26 : readPer%dstart = int(nodata_dp, i4)
313 26 : readPer%mstart = int(nodata_dp, i4)
314 26 : readPer%ystart = int(nodata_dp, i4)
315 26 : readPer%dend = int(nodata_dp, i4)
316 26 : readper%mend = int(nodata_dp, i4)
317 26 : readper%yend = int(nodata_dp, i4)
318 26 : readPer%Nobs = int(nodata_dp, i4)
319 : end if
320 :
321 : ! evaluate date and timeStep_model_inputs to get read_flag -------
322 65541 : read_flag = is_read(iDomain, tt, nTstepDay, simPer, timestep, timeStep_model_inputs)
323 : !
324 : ! determine start and end date of chunk to read
325 65541 : if (read_flag) call chunk_size(iDomain, tt, nTstepDay, simPer, timeStep_model_inputs, readPer)
326 : !
327 65541 : end subroutine chunk_config
328 :
329 :
330 : !> \brief evaluate whether new chunk should be read at this timestep
331 : !> \changelog
332 : !! - Robert Schweppe Jun 2018
333 : !! - refactoring and reformatting
334 : !! - Sebastian Müller Mar 2023
335 : !! - used by meteo-handler
336 : !! - now independent of global variables
337 : !> \authors Stephan Thober
338 : !> \date Jun 2014
339 131082 : function is_read(iDomain, tt, nTstepDay, simPer, timestep, timeStep_model_inputs)
340 :
341 65541 : use mo_common_types, only : period
342 : use mo_julian, only : caldat
343 :
344 : implicit none
345 :
346 : !> current Domain
347 : integer(i4), intent(in) :: iDomain
348 : !> current time step
349 : integer(i4), intent(in) :: tt
350 : !> Number of time intervals per day
351 : integer(i4), intent(in) :: nTstepDay
352 : !> warmPer + evalPer
353 : type(period), dimension(:), intent(in) :: simPer
354 : !> [h] simulation time step (= TS) in [h]
355 : integer(i4), intent(in) :: timeStep
356 : !> frequency for reading meteo input
357 : integer(i4), dimension(:), intent(in) :: timeStep_model_inputs
358 :
359 : logical :: is_read
360 :
361 : ! number of simulated days
362 : integer(i4) :: Ndays
363 : ! day
364 : integer(i4) :: day
365 : ! months
366 : integer(i4) :: month
367 : ! years
368 : integer(i4) :: year
369 : ! number of simulated days one timestep before
370 : integer(i4) :: Ndays_before
371 : ! day one simulated timestep before
372 : integer(i4) :: day_before
373 : ! month one simulated timestep before
374 : integer(i4) :: month_before
375 : ! year one simulated timestep before
376 : integer(i4) :: year_before
377 :
378 : ! initialize
379 65541 : is_read = .false.
380 :
381 : ! special case for first timestep
382 65541 : if (tt .eq. 1_i4) then
383 : is_read = .true.
384 : else
385 : ! check if a new day started by comparing the day of the current time step (Ndays)
386 : ! with the one before (Ndays_before)
387 65515 : Ndays = ceiling((real(tt, dp)) / real(nTstepDay, dp))
388 65515 : Ndays_before = ceiling((real(tt, dp) - 1.0_dp) / real(nTstepDay, dp))
389 :
390 : ! evaluate cases of given timeStep_model_inputs
391 65515 : select case(timeStep_model_inputs(iDomain))
392 : case(0) ! only at the beginning of the period
393 : if (tt .eq. 1_i4) is_read = .true.
394 : case(1 :) ! every timestep with frequency timeStep_model_inputs
395 0 : if (mod((tt - 1) * timestep, timeStep_model_inputs(iDomain) * 24) .eq. 0_i4) is_read = .true.
396 : case(-1) ! every day
397 0 : if (Ndays .ne. Ndays_before) is_read = .true.
398 : case(-2) ! every month
399 65515 : if (Ndays .ne. Ndays_before) then
400 : ! calculate months
401 2725 : call caldat(simPer(iDomain)%julStart + Ndays - 1, dd = day, mm = month, yy = year)
402 2725 : call caldat(simPer(iDomain)%julStart + Ndays_before - 1, dd = day_before, mm = month_before, yy = year_before)
403 2725 : if (month .ne. month_before) is_read = .true.
404 : end if
405 : case(-3) ! every year
406 0 : if (Ndays .ne. Ndays_before) then
407 : ! calculate months
408 0 : call caldat(simPer(iDomain)%julStart + Ndays - 1, dd = day, mm = month, yy = year)
409 0 : call caldat(simPer(iDomain)%julStart + Ndays_before - 1, dd = day_before, mm = month_before, yy = year_before)
410 0 : if (year .ne. year_before) is_read = .true.
411 : end if
412 : case default ! not specified correctly
413 65515 : call error_message('ERROR*** mo_meteo_helper: function is_read: timStep_model_inputs not specified correctly!')
414 : end select
415 : end if
416 :
417 65541 : end function is_read
418 :
419 :
420 : !> \brief calculate beginning and end of read Period, i.e. that is length of current chunk to read
421 : !> \changelog
422 : !! - Stephan Thober Jan 2015
423 : !! - added iDomain
424 : !! - Robert Schweppe Jun 2018
425 : !! - refactoring and reformatting
426 : !! - Sebastian Müller Mar 2023
427 : !! - used by meteo-handler
428 : !! - now independent of global variables
429 : !> \authors Stephan Thober
430 : !> \date Jun 2014
431 111 : subroutine chunk_size(iDomain, tt, nTstepDay, simPer, timeStep_model_inputs, readPer)
432 :
433 65541 : use mo_common_types, only: period
434 : use mo_julian, only : caldat, julday
435 :
436 : implicit none
437 :
438 : !> current Domain to process
439 : integer(i4), intent(in) :: iDomain
440 : !> current time step
441 : integer(i4), intent(in) :: tt
442 : !> Number of time intervals per day
443 : integer(i4), intent(in) :: nTstepDay
444 : !> warmPer + evalPer
445 : type(period), dimension(:), intent(in) :: simPer
446 : !> frequency for reading meteo input
447 : integer(i4), dimension(:), intent(in) :: timeStep_model_inputs
448 : !> start and end dates of read Period
449 : type(period), intent(out) :: readPer
450 :
451 : ! number of simulated days
452 : integer(i4) :: Ndays
453 : ! day
454 : integer(i4) :: day
455 : ! months
456 : integer(i4) :: month
457 : ! years
458 : integer(i4) :: year
459 :
460 : ! calculate date of start date
461 111 : Ndays = ceiling(real(tt, dp) / real(nTstepDay, dp))
462 :
463 : ! get start date
464 111 : readPer%julStart = simPer(iDomain)%julStart + Ndays - 1
465 :
466 : ! calculate end date according to specified frequency
467 111 : select case (timeStep_model_inputs(iDomain))
468 : case(0) ! length of chunk has to cover whole period
469 21 : readPer%julEnd = simPer(iDomain)%julEnd
470 : case(1 :) ! every timestep with frequency timeStep_model_inputs
471 0 : readPer%julEnd = readPer%julStart + timeStep_model_inputs(iDomain) - 1
472 : case(-1) ! every day
473 0 : readPer%julEnd = readPer%julStart
474 : case(-2) ! every month
475 : ! calculate date
476 90 : call caldat(simPer(iDomain)%julStart + Ndays, dd = day, mm = month, yy = year)
477 : ! increment month
478 90 : if (month .eq. 12) then
479 5 : month = 1
480 5 : year = year + 1
481 : else
482 85 : month = month + 1
483 : end if
484 90 : readPer%julEnd = julday(dd = 1, mm = month, yy = year) - 1
485 : case(-3) ! every year
486 : ! calculate date
487 0 : call caldat(simPer(iDomain)%julStart + Ndays, dd = day, mm = month, yy = year)
488 0 : readPer%julEnd = julday(dd = 31, mm = 12, yy = year)
489 : case default ! not specified correctly
490 111 : call error_message('ERROR*** mo_meteo_helper: chunk_size: timStep_model_inputs not specified correctly!')
491 : end select
492 :
493 : ! end date should not be greater than end of simulation period
494 111 : readPer%julEnd = min(readPer%julEnd, simPer(iDomain)%julEnd)
495 :
496 : ! calculate the dates of the start and end dates
497 111 : call caldat(readPer%julStart, dd = readPer%dstart, mm = readPer%mstart, yy = readPer%ystart)
498 111 : call caldat(readPer%julEnd, dd = readPer%dEnd, mm = readPer%mend, yy = readPer%yend)
499 111 : readPer%Nobs = readPer%julEnd - readPer%julstart + 1
500 :
501 111 : end subroutine chunk_size
502 :
503 : END MODULE mo_meteo_helper
|
