Line data Source code
1 : !> \file mo_read_wrapper.f90
2 : !> \brief \copybrief mo_read_wrapper
3 : !> \details \copydetails mo_read_wrapper
4 :
5 : !> \brief Wrapper for all reading routines.
6 : !> \details This module is to wrap up all reading routines.
7 : !! The general written reading routines are used to store now the read data into global variables.
8 : !> \authors Juliane Mai, Matthias Zink
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_mpr
13 : MODULE mo_read_wrapper
14 :
15 : USE mo_kind, ONLY : i4, dp
16 : use mo_common_constants, only : nodata_dp, nodata_i4
17 : use mo_message, only: message, error_message
18 :
19 : IMPLICIT NONE
20 :
21 : PUBLIC :: read_data ! reads all available data
22 :
23 : CONTAINS
24 :
25 : ! ------------------------------------------------------------------
26 :
27 : ! NAME
28 : ! read_data
29 :
30 : ! PURPOSE
31 : !> \brief Reads data.
32 :
33 : !> \details The namelists are already read by read_config call.
34 : !> All LUTs are read from their respective directory and information within those
35 : !> files are shared across all domains to be modeled.
36 :
37 : ! INTENT(IN), OPTIONAL
38 : !> \param[in] "type(period), dimension(:), optional :: LAIPer"
39 :
40 : ! HISTORY
41 : !> \authors Juliane Mai & Matthias Zink
42 :
43 : !> \date Feb 2013
44 :
45 : ! Modifications:
46 : ! Luis Samaniego Feb 2013 - rotate fdir variable to the new coordinate system
47 : ! Rohini Kumar Aug 2013 - name changed from "L0_LAI" to "L0_LCover_LAI"
48 : ! Rohini Kumar Aug 2013 - added dirSoil_LUT and dirGeology_LUT, and changed to
49 : ! read datapaths and variables made accordingly
50 : ! Rohini Kumar Aug 2013 - added iFlag_LAI_data_format to handle LAI options, and changed
51 : ! within the code made accordingly
52 : ! Rohini Kumar Sep 2013 - read input data for routing processes according to process_matrix flag
53 : ! Matthias Zink Mar 2014 - added inflow gauge
54 : ! Kumar & Schroen Apr 2014 - added check for consistency of L0 and L1 spatial resolution
55 : ! Stephan Thober Jun 2014 - added perform_mpr for omitting L0 read
56 : ! Matthias Cuntz & Juliane Mai Nov 2014 - LAI input from daily, monthly or yearly files
57 : ! Stephan Thober Aug 2015 - moved routing related variables and routines to mRM
58 : ! Rohini Kumar Mar 2016 - options to handle different soil databases
59 : ! Matthias Zink Mar 2014 - added subroutine for consistency check
60 : ! Stephan Thober Nov 2016 - moved processMatrix to common variables
61 : ! Rohini Kumar Dec 2016 - option to handle monthly mean gridded fields of LAI
62 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
63 :
64 26 : subroutine read_data(LAIPer)
65 :
66 : use mo_append, only : append, paste
67 : use mo_constants, only : YearMonths
68 : use mo_common_read_data, only : read_dem, read_lcover
69 : use mo_common_types, only: period, Grid
70 : use mo_common_variables, only : dirCommonFiles, dirMorpho, &
71 : global_parameters, level0, domainMeta, processMatrix
72 : use mo_mpr_file, only : file_aspect, file_geolut, file_hydrogeoclass, &
73 : file_laiclass, file_lailut, file_slope, file_soil_database, file_soil_database_1, &
74 : file_soilclass, uaspect, ugeolut, uhydrogeoclass, ulaiclass, ulailut, uslope, usoilclass
75 : use mo_mpr_global_variables, only : GeoUnitKar, &
76 : GeoUnitList, L0_asp, L0_geoUnit, L0_gridded_LAI, L0_slope, L0_soilId, LAILUT, &
77 : LAIUnitList, iFlag_soilDB, nGeoUnits, nLAI, nLAIclass, nSoilHorizons_mHM, soilDB, &
78 : timeStep_LAI_input, LAIBoundaries
79 : use mo_prepare_gridded_lai, only : prepare_gridded_daily_LAI_data, prepare_gridded_mean_monthly_LAI_data
80 : use mo_read_latlon, only : read_latlon
81 : use mo_read_lut, only : read_geoformation_lut, read_lai_lut
82 : use mo_read_spatial_data, only : read_spatial_data_ascii
83 : use mo_soil_database, only : read_soil_LUT
84 : use mo_string_utils, only : num2str
85 : use mo_timer, only : timer_get, timer_start, &
86 : timer_stop
87 :
88 : implicit none
89 :
90 : type(period), dimension(:), intent(in), optional :: LAIPer
91 :
92 : ! loop variables
93 : integer(i4) :: domainID, iDomain, iVar, iHorizon, iMon, itimer, ll
94 :
95 : ! dummy variable
96 : integer(i4) :: nH
97 :
98 : ! file unit of file to read
99 : integer(i4) :: nunit
100 :
101 : ! file name of file to read
102 : character(256) :: fName
103 :
104 13 : real(dp), dimension(:, :), allocatable :: data_dp_2d
105 :
106 13 : integer(i4), dimension(:, :), allocatable :: data_i4_2d
107 :
108 13 : integer(i4), dimension(:, :), allocatable :: dataMatrix_i4
109 :
110 13 : logical, dimension(:, :), allocatable :: mask_2d
111 :
112 13 : integer(i4), dimension(:), allocatable :: dummy_i4
113 :
114 : type(Grid), pointer :: level0_iDomain
115 :
116 : real(dp), parameter :: slope_minVal = 0.01_dp
117 :
118 : real(dp), parameter :: aspect_minVal = 1.00_dp
119 :
120 :
121 13 : call message(' Reading data ...')
122 13 : itimer = 1
123 13 : call timer_start(itimer)
124 :
125 13 : call message(' Reading dem and lcover ...')
126 13 : call read_dem()
127 13 : call read_lcover()
128 :
129 : ! ************************************************
130 : ! READ LOOKUP TABLES
131 : ! ************************************************
132 : !
133 : ! Soil LUT
134 13 : if(iFlag_soilDB .eq. 0) then
135 13 : fName = trim(adjustl(dirCommonFiles)) // trim(adjustl(file_soil_database))
136 0 : else if(iFlag_soilDB .eq. 1) then
137 0 : fName = trim(adjustl(dirCommonFiles)) // trim(adjustl(file_soil_database_1))
138 : end if
139 13 : call read_soil_LUT(trim(fName))
140 :
141 : ! Geological formation LUT
142 13 : fName = trim(adjustl(dirCommonFiles)) // trim(adjustl(file_geolut))
143 13 : call read_geoformation_lut(trim(fName), ugeolut, nGeoUnits, GeoUnitList, GeoUnitKar)
144 :
145 : ! LAI LUT
146 13 : if(timeStep_LAI_input .EQ. 0) then
147 12 : fName = trim(adjustl(dirCommonFiles)) // trim(adjustl(file_lailut))
148 12 : call read_lai_lut(trim(fName), ulailut, nLAIclass, LAIUnitList, LAILUT)
149 : end if
150 :
151 38 : domains: do iDomain = 1, domainMeta%nDomains
152 25 : domainID = domainMeta%indices(iDomain)
153 :
154 25 : level0_iDomain => level0(domainMeta%L0DataFrom(iDomain))
155 :
156 25 : call message(' Reading data for domain: ', trim(adjustl(num2str(domainID))), ' ...')
157 : ! check whether L0 data is shared
158 : ! ToDo: check change
159 25 : if (domainMeta%L0DataFrom(iDomain) < iDomain) then
160 : !
161 : call message(' Using data of domain ', &
162 0 : trim(adjustl(num2str(domainMeta%indices(domainMeta%L0DataFrom(iDomain))))), ' for domain: ',&
163 3 : trim(adjustl(num2str(domainID))), '...')
164 : ! DO NOT read L0 data
165 3 : cycle
166 :
167 : end if
168 :
169 22 : itimer = 2
170 22 : call timer_start(itimer)
171 :
172 : ! read slope and aspect - datatype real
173 66 : nVars_real : do iVar = 1, 2
174 22 : select case (iVar)
175 : case(1) ! slope
176 22 : call message(' Reading slope ...')
177 22 : fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_slope))
178 22 : nunit = uslope
179 : case(2) ! aspect
180 22 : call message(' Reading aspect ...')
181 22 : fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_aspect))
182 66 : nunit = uaspect
183 : end select
184 :
185 : ! reading
186 : call read_spatial_data_ascii(trim(fName), nunit, &
187 : level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
188 44 : level0_iDomain%yllcorner, level0_iDomain%cellsize, data_dp_2d, mask_2d)
189 : ! put global nodata value into array (probably not all grid cells have values)
190 5225368 : data_dp_2d = merge(data_dp_2d, nodata_dp, mask_2d)
191 : ! put data in variable
192 : select case (iVar)
193 : case(1) ! slope
194 22 : call append(L0_slope, pack(data_dp_2d, level0_iDomain%mask))
195 :
196 : case(2) ! aspect
197 44 : call append(L0_asp, pack(data_dp_2d, level0_iDomain%mask))
198 : end select
199 : ! deallocate arrays
200 66 : deallocate(data_dp_2d, mask_2d)
201 :
202 : end do nVars_real
203 :
204 : ! read datatype integer
205 : ! ***** CHANGE IS MADE HERE TO ACCOMODATE READ SEVERAL TYPES OF SOIL DATABASES
206 : ! change from the earlier code where everything was done in the DO loop (**see below **)
207 : ! here everything is more explicit and no do loop appears as was the case in previous version
208 :
209 : ! read soilID in both options
210 22 : call message(' Reading soil ids ...')
211 :
212 22 : nH = 1 !> by default; when iFlag_soilDB = 0
213 22 : if(iFlag_soilDB .eq. 1) nH = nSoilHorizons_mHM
214 : ! modified way to read multiple horizons specific soil class
215 44 : do iHorizon = 1, nH
216 22 : if(iFlag_soilDB .eq. 0) then
217 22 : fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_soilclass))
218 0 : else if(iFlag_soilDB .eq. 1) then
219 0 : write(fName, 172) iHorizon
220 : 172 format('soil_class_horizon_', i2.2, '.asc')
221 0 : fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(fName))
222 : end if
223 : call read_spatial_data_ascii(trim(fName), usoilclass, &
224 : level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
225 22 : level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
226 : ! put global nodata value into array (probably not all grid cells have values)
227 2612684 : data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
228 22 : call paste(dataMatrix_i4, pack(data_i4_2d, level0_iDomain%mask), nodata_i4)
229 44 : deallocate(data_i4_2d)
230 : end do
231 22 : call append(L0_soilId, dataMatrix_i4)
232 22 : deallocate(dataMatrix_i4)
233 :
234 : ! read geoUnit
235 22 : fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_hydrogeoclass))
236 : ! reading and transposing
237 : call read_spatial_data_ascii(trim(fName), uhydrogeoclass, &
238 : level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
239 22 : level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
240 : ! put global nodata value into array (probably not all grid cells have values)
241 2612684 : data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
242 22 : call append(L0_geoUnit, pack(data_i4_2d, level0_iDomain%mask))
243 22 : deallocate(data_i4_2d, mask_2d)
244 :
245 : ! LAI values
246 22 : call message(' Reading LAI ...')
247 : select case (timeStep_LAI_input)
248 : case(1) ! long term mean monthly gridded fields
249 0 : call prepare_gridded_mean_monthly_LAI_data(iDomain, level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%mask)
250 :
251 : case(0) ! long term mean monthly values per class with LUT
252 : ! only set if not yet allocated (e.g. domain 1)
253 21 : if (.not. allocated(LAIBoundaries)) then
254 12 : nLAI = int(YearMonths, i4)
255 12 : allocate(LAIBoundaries(nLAI+1))
256 492 : LAIBoundaries = [(iMon, iMon=1, nLAI+1)]
257 : end if
258 :
259 21 : fName = trim(adjustl(dirMorpho(iDomain))) // trim(adjustl(file_laiclass))
260 : ! reading and transposing
261 : call read_spatial_data_ascii(trim(fName), ulaiclass, &
262 : level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%xllcorner, &
263 21 : level0_iDomain%yllcorner, level0_iDomain%cellsize, data_i4_2d, mask_2d)
264 : ! put global nodata value into array (probably not all grid cells have values)
265 2487834 : data_i4_2d = merge(data_i4_2d, nodata_i4, mask_2d)
266 2487855 : allocate(dummy_i4(count(level0_iDomain%mask)))
267 21 : dummy_i4 = pack(data_i4_2d, level0_iDomain%mask)
268 21 : deallocate(data_i4_2d, mask_2d)
269 :
270 : call check_consistency_lut_map(dummy_i4, LAIUnitList, file_laiclass)
271 :
272 2487876 : allocate(data_dp_2d(count(level0_iDomain%mask), nLAI))
273 273 : do iMon = 1, nLAI
274 : ! determine LAIs per month
275 2793 : do ll = 1, size(LAILUT, dim = 1)
276 111171372 : data_dp_2d(:, iMon) = merge(LAILUT(ll, iMon), data_dp_2d(:, iMon), dummy_i4(:) .EQ. LAIUnitList(ll))
277 : end do
278 : end do
279 21 : call append(L0_gridded_LAI, data_dp_2d(:, :))
280 21 : deallocate(dummy_i4, data_dp_2d)
281 : ! correction for 0 LAI values to avoid numerical instabilities
282 24161913 : L0_gridded_LAI(:, :) = merge(1.00E-10_dp, L0_gridded_LAI(:, :), L0_gridded_LAI(:, :) .LT. 1.00E-10_dp)
283 24161934 : L0_gridded_LAI(:, :) = merge(30.0_dp, L0_gridded_LAI(:, :), L0_gridded_LAI(:, :) .GT. 30.0_dp)
284 : case(-3 : -1) ! daily, monthly or yearly gridded fields (time-series)
285 : call prepare_gridded_daily_LAI_data(iDomain, level0_iDomain%nrows, level0_iDomain%ncols, level0_iDomain%mask, &
286 22 : LAIPer(iDomain))
287 :
288 : end select
289 :
290 : ! read lat lon coordinates of each domain
291 22 : call message(' Reading latitude/logitude ...')
292 22 : call read_latlon(iDomain, "lon_l0", "lat_l0", "level0", level0_iDomain)
293 :
294 22 : call timer_stop(itimer)
295 57 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
296 :
297 : end do domains
298 13 : itimer = 1
299 13 : call timer_stop(itimer)
300 13 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
301 :
302 : !Soil
303 : ! determine name od soil class definition file based on input option iFlag_soilDB
304 13 : if(iFlag_soilDB .eq. 0) then
305 13 : fName = file_soil_database
306 0 : else if(iFlag_soilDB .eq. 1) then
307 0 : fName = file_soil_database_1
308 : end if
309 :
310 : ! If you are getting segmentation fault with intel on large file (e.g. here on soils), then
311 : ! make sure to include following into your Marefile: INTEL_EXCLUDE := mo_read_wrapper.f90
312 0 : call check_consistency_lut_map(reshape(L0_soilId, (/ size(L0_soilId, 1) * size(L0_soilId, 2) /)), &
313 26 : soilDB%id(:), fName)
314 : ! Geology
315 13 : call check_consistency_lut_map(L0_geoUnit, GeoUnitList, file_hydrogeoclass, dummy_i4)
316 :
317 : ! deactivate parameters of non existing geological classes in study domain for optimization
318 : ! loop over geological units in look up list
319 143 : do iVar = 1, size(GeoUnitList, 1)
320 : ! check if unit appears in geological map (dummy_i4 is unique number in L0_geoUnit)
321 689 : if (.not. ANY(dummy_i4 .EQ. GeoUnitList(iVar))) then
322 : ! deactivate optimization flag (dim=4 from global_parameters)
323 39 : global_parameters(processMatrix(9, 3) - processMatrix(9, 2) + iVar, 4) = 0
324 39 : call message('***WARNING: Geological unit ', trim(adjustl(num2str(GeoUnitList(iVar)))))
325 39 : call message(' is not appearing in study domain.')
326 : end if
327 : end do
328 :
329 13 : deallocate(dummy_i4) ! is allocated in subroutine check_consistency_lut_map - geology
330 :
331 : !----------------------------------------------------------------
332 : ! Correction for slope and aspect -- min value set above
333 : !----------------------------------------------------------------
334 : ! keep the colons (:) in the statements because of Intel's reallocation lhs problem
335 972963 : L0_slope(:) = merge(slope_minVal, L0_slope(:), (L0_slope(:) .lt. slope_minVal))
336 972963 : L0_asp(:) = merge(aspect_minVal, L0_asp(:), (L0_asp(:) .lt. aspect_minVal))
337 :
338 13 : end subroutine read_data
339 :
340 : ! ------------------------------------------------------------------
341 :
342 : ! NAME
343 : ! check_consistency_lut_map
344 :
345 : ! PURPOSE
346 : !> \brief Checks if classes in input maps appear in look up tables.
347 :
348 : !> \details Determines wether a class appearing in the morphological input
349 : !> is occuring in the respective look up table. mHM breaks if inconsistencies
350 : !> are discovered.
351 :
352 : ! INTENT(IN)
353 : !> \param[in] "integer(i4), dimension(:) :: data" map of study domain
354 : !> \param[in] "integer(i4), dimension(:) :: lookuptable" look up table corresponding to map
355 : !> \param[in] "character(*) :: filename" name of the lut file - ERORR warn
356 :
357 : ! INTENT(OUT), OPTIONAL
358 : !> \param[out] "integer(i4), dimension(:), optional :: unique_values" array of unique values in dataone
359 :
360 : ! HISTORY
361 : !> \authors Matthias Zink
362 :
363 : !> \date Nov 2016
364 :
365 : ! Modifications:
366 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
367 :
368 94 : subroutine check_consistency_lut_map(data, lookuptable, filename, unique_values)
369 :
370 13 : use mo_orderpack, only : unista
371 : use mo_string_utils, only : num2str
372 :
373 : implicit none
374 :
375 : ! map of study domain
376 : integer(i4), dimension(:), intent(in) :: data
377 :
378 : ! look up table corresponding to map
379 : integer(i4), dimension(:), intent(in) :: lookuptable
380 :
381 : ! name of the lut file - ERORR warn
382 : character(*), intent(in) :: filename
383 :
384 : ! array of unique values in dataone
385 : integer(i4), dimension(:), allocatable, intent(out), optional :: unique_values
386 :
387 : integer(i4) :: n_unique_elements
388 :
389 : integer(i4) :: ielement
390 :
391 47 : integer(i4), dimension(:), allocatable :: temp
392 :
393 :
394 141 : allocate(temp(size(data, 1)))
395 : ! copy L0_geoUnit because subroutine unista overwrites the nunit entries of the
396 : ! input array with the unique array values
397 2872399 : temp = data
398 : ! retrieve unique values of data
399 47 : call unista(temp, n_unique_elements)
400 : ! check if unit exists in look up table
401 824 : do ielement = 1, n_unique_elements
402 777 : if ( temp(ielement) == nodata_i4 ) call error_message( &
403 : '***ERROR: Class ', trim(adjustl(num2str(temp(ielement)))), &
404 : ' was searched in ', trim(adjustl(filename)), &
405 : ' which indicates a masking problem!' &
406 0 : )
407 310520 : if (.not. ANY(lookuptable .EQ. temp(ielement))) then
408 0 : call error_message('***ERROR: Class ', trim(adjustl(num2str(temp(ielement)))), ' is missing', raise=.false.)
409 0 : call error_message(' in input file ', trim(adjustl(filename)), ' ...')
410 : end if
411 : end do
412 :
413 : ! pass unique values if optional argument unique_values is given
414 47 : if (present(unique_values)) then
415 39 : allocate(unique_values(n_unique_elements))
416 117 : unique_values(:) = temp(1 : n_unique_elements)
417 : end if
418 47 : deallocate(temp)
419 :
420 94 : end subroutine check_consistency_lut_map
421 :
422 :
423 : END MODULE mo_read_wrapper
|