Line data Source code
1 : !> \file mo_edk_read_data.f90
2 : !> \brief \copybrief mo_edk_read_data
3 : !> \details \copydetails mo_edk_read_data
4 :
5 : !> \brief Module containing routines to read data.
6 : !> \author Luis Samaniego
7 : !> \date 21.03.2006
8 : !> \date 11.06.2010
9 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
10 : !! EDK is released under the LGPLv3+ license \license_note
11 : module mo_edk_read_data
12 :
13 : use mo_kind, only: dp
14 :
15 : implicit none
16 :
17 : private
18 :
19 : public :: ReadData
20 :
21 : !> \class extend
22 : !> \brief domain extend description
23 : type extend
24 : real(dp) :: left !< left border
25 : real(dp) :: top !< top border
26 : real(dp) :: right !< right border
27 : real(dp) :: bottom !< bottom border
28 : end type extend
29 :
30 : contains
31 :
32 : !> \brief Routine to read data for EDK.
33 : !> \details Read Database Precipitation
34 : !! Read grid DEM
35 : !! Reads parameters variogram
36 : !> \author Luis Samaniego
37 : !> \date 21.03.2006
38 : !> \date 11.06.2010
39 2 : subroutine ReadData()
40 : use mainVar, only: noDataValue, cellFactor, DataConvertFactor, OffSet, &
41 : yStart, mStart, dStart, yEnd, mEnd, dEnd, tBuffer,DEMNcFlag
42 : use mo_kind, only: i4, dp
43 : use mo_julian , only: NDAYS
44 :
45 : use kriging, only: maxDist
46 : use VarFit, only: vType, nParam, dh, hMax, beta
47 : use runControl, only: interMth, fnameDEM, DataPathOut, DataPathIn, fNameSta, correctNeg, &
48 : distZero, flagVario, fNameVario, flagEDK
49 : use NetCDFVar, only: FileOut, originator, crs, source, institution, title, contact,&
50 : variable_name, variable_unit, &
51 : variable_standard_name, variable_calendar_type, &
52 : variable_long_name, ncIn_variable_name, ncIn_dem_variable_name, &
53 : ncIn_yCoord_name, ncIn_xCoord_name, invert_y, &
54 : ncOut_dem_variable_name, ncOut_dem_yCoord_name, ncOut_dem_xCoord_name, &
55 : ncOut_dem_Latitude, ncOut_dem_Longitude
56 : use mo_message, only: message
57 : use mo_string_utils, only: divide_string
58 : use mo_edk_info, only: file_namelist
59 :
60 : implicit none
61 : !
62 : logical :: isNcFile
63 : logical :: isDEMNcFile
64 : integer(i4) :: i, ios
65 : character(256) :: dummy, fileName
66 2 : character(256), allocatable :: DataPathIn_parts(:)
67 2 : character(256), allocatable :: fNameDEM_parts(:)
68 : !
69 : !===============================================================
70 : ! namelist definition
71 : !===============================================================
72 : !
73 : namelist/mainVars/noDataValue, DataPathIn, fNameDEM, &
74 : DataPathOut, FileOut, fNameSTA, cellFactor, DataConvertFactor, OffSet, InterMth, correctNeg, &
75 : distZero, originator, crs, source, institution, title, contact,&
76 : variable_name, variable_unit, variable_long_name, &
77 : variable_standard_name, variable_calendar_type, &
78 : yStart, mStart, dStart, yEnd, mEnd, dEnd,tBuffer, maxDist, flagVario, vType, nParam, &
79 : fNameVario, dh, hMax, ncIn_variable_name, ncIn_dem_variable_name, ncIn_yCoord_name, &
80 : ncIn_xCoord_name, invert_y, &
81 : ncOut_dem_variable_name, ncOut_dem_yCoord_name, ncOut_dem_xCoord_name, ncOut_dem_Latitude,&
82 : ncOut_dem_Longitude
83 : !
84 : ! -----------------------------------------------------------------------
85 : ! MAIN.DAT
86 : ! -----------------------------------------------------------------------
87 2 : open(unit=10, file=file_namelist, STATUS='OLD', ACTION='read')
88 2 : read(10, nml=mainVars)
89 2 : close(10)
90 :
91 : ! -----------------------------------------------------------------------
92 : ! consistency checks
93 : ! -----------------------------------------------------------------------
94 2 : if ((interMth .gt. 2) .or. (interMth .lt. 0)) then
95 0 : call message(' >>> Interpolated Method Value not valid')
96 0 : call message(' Please choose one of the following:')
97 0 : call message(' 0 - no interpolation')
98 0 : call message(' 1 - ordinary kriging')
99 0 : call message(' 2 - external drift kriging')
100 0 : call message('')
101 0 : stop 1
102 : end if
103 :
104 2 : FileOut = trim(DataPathOut) // trim(FileOut)
105 : !
106 : ! *************************************
107 : ! Read initial control parameters
108 : ! *************************************
109 :
110 2 : if (flagVario .AND. flagEDK) then
111 0 : print*, '***Warning: Both flags flagVario and flagEDK should not be activated at the same time!'
112 : ! stop
113 : end if
114 :
115 2 : ios = 0
116 6 : allocate (beta(nParam))
117 :
118 2 : if ( .not. flagVario ) then
119 1 : filename = trim(fNameVARIO)
120 1 : open (unit=20, file=filename, STATUS='OLD', ACTION='READ')
121 1 : read (20,1) dummy
122 4 : read (20, *,iostat=ios) (beta(i), i=1,nParam), i
123 1 : close (20)
124 : end if
125 :
126 : ! if vario is read from file take read in vario type
127 2 : if ( .not. flagVario ) then
128 1 : vType = i
129 : end if
130 2 : if (.NOT. ((vType .EQ. 1) .OR. (vType .EQ. 2))) then
131 0 : print*, "***ERROR: Variogram type not known: ", vType
132 : end if
133 :
134 2 : write(*,*) 'variogram parameters: '
135 2 : write(*,'(3(a12), a18)') 'nugget', 'sill', 'range', 'Variogramtype'
136 8 : write(*,'(3(f12.2), i18)') (beta(i), i=1,nParam), vType
137 :
138 : ! determine whether DEM is an nc file or not
139 2 : call divide_string(fNameDEM, '.', fNameDEM_parts)
140 2 : isDEMNcFile = (fNameDEM_parts(size(fNameDEM_parts,1)) .eq. 'nc')
141 2 : deallocate(fNameDEM_parts)
142 2 : DEMNcFlag = 0
143 : ! read DEM
144 2 : if (isDEMNcFile) then
145 : ! read the netcdf DEM file
146 0 : call ReadDEMNc
147 0 : DEMNcFlag = 1
148 : else
149 2 : call ReadDEM
150 : end if
151 :
152 : ! determine whether input path is actually a nc file
153 2 : call divide_string(DataPathIn, '.', DataPathIn_parts)
154 2 : isNcFile = (DataPathIn_parts(size(DataPathIn_parts,1)) .eq. 'nc')
155 2 : deallocate(DataPathIn_parts)
156 :
157 2 : if (isNcFile) then
158 : ! read 3d meteo file
159 0 : call initMetStaNc
160 :
161 : else
162 : ! read look up table for meteorological stations
163 2 : call ReadStationLut
164 :
165 : ! read whole METEO data
166 2 : call ReadDataMeteo
167 : end if
168 :
169 : 1 format (a256)
170 :
171 6 : end subroutine ReadData
172 :
173 :
174 2 : subroutine ReadStationLut
175 2 : use mo_kind, only: i4, dp
176 :
177 : use mainVar, only: nSta, MetSta, grid, noDataValue
178 : use kriging, only: maxDist, edk_dist
179 : use VarFit, only: hmax
180 : use runControl
181 :
182 : implicit none
183 : !
184 : integer(i4) :: iSta
185 : integer(i4) :: n_stations_all ! number of stations in file
186 : character(256) :: dummy, filename
187 2 : real(dp), dimension(:), allocatable :: temp_id ! ID of read in station
188 2 : real(dp), dimension(:), allocatable :: temp_x ! easting coordinate of station
189 2 : real(dp), dimension(:), allocatable :: temp_y ! northing coordinate of station
190 2 : real(dp), dimension(:), allocatable :: temp_h ! heigth of station
191 : type(extend) :: domain
192 :
193 : !-----------------------------------------------------------------------
194 : ! STATIONS COORDINATES AND ELEVATION
195 : !-----------------------------------------------------------------------
196 :
197 2 : call getSearchDomain(domain)
198 :
199 2 : filename = trim(fNameSTA)
200 2 : open (unit=20, file=filename, status='old', action='read')
201 2 : read (20,*) dummy, n_stations_all
202 2 : read (20, 1) dummy
203 :
204 : ! allocate
205 8 : allocate(temp_id(n_stations_all)); allocate(temp_x(n_stations_all))
206 6 : allocate( temp_y(n_stations_all)); allocate(temp_h(n_stations_all))
207 :
208 : ! initialize
209 82 : temp_id = grid%nodata_value; temp_h = grid%nodata_value
210 82 : temp_x = grid%nodata_value; temp_y = grid%nodata_value
211 : ! init
212 2 : nSta = 1
213 :
214 : ! new format (reads only part of the line)
215 40 : do iSta = 1 , n_stations_all
216 38 : read (20, *) temp_id(nSta), temp_x(nSta), temp_y(nSta), temp_h(nSta)
217 :
218 : ! check if station is within search distance
219 : ! if yes increase station counting by one
220 : ! otherwise overwrite station in next iteration
221 38 : if ((temp_x(nSta) .GE. domain%left) .AND. (temp_x(nSta) .LE. domain%right) .AND. &
222 78 : (temp_y(nSta) .GE. domain%bottom) .AND. (temp_y(nSta) .LE. domain%top) ) then
223 38 : nSta = nSta + 1
224 : end if
225 : end do
226 :
227 2 : nSta = nSta - 1
228 :
229 2 : close(20)
230 :
231 : ! save relevant stations in MetSta
232 44 : allocate(MetSta(nSta))
233 2 : edk_dist%nstat = nSta
234 2 : edk_dist%noDataValue = noDataValue
235 6 : allocate(edk_dist%stat_pos(nSta, 2))
236 :
237 40 : do iSta = 1, nSta
238 38 : MetSta(iSta)%Id = temp_id(iSta)
239 38 : MetSta(iSta)%x = temp_x(iSta)
240 38 : MetSta(iSta)%y = temp_y(iSta)
241 38 : MetSta(iSta)%h = temp_h(iSta)
242 116 : edk_dist%stat_pos(iSta, :) = [MetSta(iSta)%x, MetSta(iSta)%y]
243 : end do
244 :
245 2 : deallocate(temp_id, temp_x, temp_y, temp_h)
246 :
247 2 : print*, "Only " , nSta , ' of the given ' , n_stations_all, ' are considered since they are within the search distance.'
248 :
249 : ! ! test filtering
250 : ! open(20, file='station_test.txt', status='unknown', action='write')
251 : ! write(20,*) dummy
252 : ! do iSta = 1 , nSta
253 : ! write(20,"(i10,' ',f12.2,' ',f12.2,' ',f12.2)") MetSta(iSta)%Id, MetSta(iSta)%x, MetSta(iSta)%y, MetSta(iSta)%h
254 : ! end do
255 : ! close(20)
256 : ! !pause
257 :
258 : ! formats
259 : 1 format (a256)
260 :
261 2 : end subroutine ReadStationLut
262 :
263 : ! -----------------------------------------------------------------------
264 : ! DEM Blocks (all must be EQUAL in size!)
265 : ! -----------------------------------------------------------------------
266 2 : subroutine ReadDEM
267 2 : use mainVar
268 : use mo_kind, only : i4
269 : use runControl
270 :
271 : implicit none
272 :
273 : integer(i4) :: i, j
274 : character(256) :: dummy, fileName
275 : !
276 :
277 2 : fileName = trim(fNameDEM)
278 2 : call readHeader(15, fileName)
279 8 : if (.not. allocated(G) ) allocate (G(grid%nrows , grid%ncols))
280 802 : do i=1,grid%nrows
281 192802 : read (15, *) (G(i,j)%h, j=1,grid%ncols)
282 : end do
283 :
284 2 : close (15)
285 : ! Assumed all blocks equal
286 2 : nCell = gridMeteo%ncols * gridMeteo%nrows
287 : !
288 2 : print *, 'DEM read ...'
289 : ! save gridMeteo*
290 2 : fileName= trim(dataPathOut)//trim('header.txt')
291 2 : call writeHeader(20, fileName)
292 : !
293 : ! formats
294 : 100 format (i1.1, a4)
295 : 101 format (i2.2, a4)
296 :
297 2 : end subroutine ReadDEM
298 :
299 : ! -----------------------------------------------------------------------
300 : ! DEM Blocks (for Irregular grids in NetCDF format )
301 : ! -----------------------------------------------------------------------
302 0 : subroutine ReadDEMNc
303 2 : use mainVar
304 : use mo_kind, only : i4
305 : use runControl
306 : use mo_netCDF, only:NcDataset, NcVariable
307 : use NetCDFVar, only:ncOut_dem_variable_name, ncOut_dem_yCoord_name, &
308 : ncOut_dem_xCoord_name, ncOut_dem_Latitude, ncOut_dem_Longitude
309 :
310 : implicit none
311 :
312 : integer(i4) :: i,j
313 : character(256) :: dummy, fileName
314 0 : real(dp), allocatable :: dem_x(:,:) ! easting
315 0 : real(dp), allocatable :: dem_y(:,:) ! northing
316 0 : real(dp), allocatable :: dem_data(:,:) ! dem data
317 : type(NcDataset) :: ncin
318 : type(NcVariable) :: ncvar
319 :
320 0 : fileName = trim(fNameDEM)
321 :
322 : ! get northing and easting information
323 0 : ncin = NcDataset(fileName, 'r')
324 0 : ncvar = ncin%getVariable(ncOut_dem_yCoord_name)
325 0 : call ncvar%getData(dem_y)
326 0 : ncvar = ncin%getVariable(ncOut_dem_xCoord_name)
327 0 : call ncvar%getData(dem_x)
328 :
329 : ! get number of rows and columns
330 0 : grid%nrows = size(dem_x, dim=1)
331 0 : grid%ncols = size(dem_x, dim=2)
332 :
333 : !write(*,*),"Nrows: ",grid%nrows
334 : !write(*,*),"Ncols: ",grid%ncols
335 :
336 : ! get latitude and longitude
337 0 : if(.not. allocated(gridMeteo%latitude)) allocate(gridMeteo%latitude(grid%ncols))
338 0 : if(.not. allocated(gridMeteo%longitude)) allocate(gridMeteo%longitude(grid%nrows))
339 :
340 0 : ncvar = ncin%getVariable(ncOut_dem_Latitude)
341 0 : call ncvar%getData(gridMeteo%latitude)
342 0 : ncvar = ncin%getVariable(ncOut_dem_Longitude)
343 0 : call ncvar%getData(gridMeteo%longitude)
344 :
345 : ! read in the dem data
346 0 : if (.not. allocated(G)) allocate(G(grid%nrows, grid%ncols))
347 0 : if (.not. allocated(gridMeteo%easting)) allocate(gridMeteo%easting(grid%nrows, grid%ncols))
348 0 : if (.not. allocated(gridMeteo%northing)) allocate(gridMeteo%northing(grid%nrows, grid%ncols))
349 : !if (.not. allocated(grid%easting)) allocate(grid%easting(grid%nrows, grid%ncols))
350 : !if (.not. allocated(grid%northing)) allocate(grid%northing(grid%nrows, grid%ncols))
351 :
352 :
353 0 : ncvar = ncin%getVariable(ncOut_dem_variable_name)
354 0 : call ncvar%getData(dem_data)
355 0 : call ncvar%getAttribute('_FillValue',grid%nodata_value)
356 :
357 : ! store the dem data in G%h
358 0 : do i=1,grid%ncols
359 0 : do j=1,grid%nrows
360 0 : G(j,i)%h = dem_data(j,i)
361 0 : gridMeteo%northing(j,i) = dem_y(j,i)
362 0 : gridMeteo%easting(j,i) = dem_x(j,i)
363 : !grid%northing(i,j) = dem_y(i,j)
364 : !grid%easting(i,j) = dem_x(i,j)
365 : end do
366 : end do
367 :
368 :
369 : !write(*,*),"DEM data: ",shape(G%h)
370 : !write(*,*),"Northing: ",shape(gridMeteo%northing)
371 : !write(*,*),"Easting: ",shape(gridMeteo%easting)
372 : !write(*,*),"Northing,Easting First Row, First Column: ",gridMeteo%northing(1,1),",",gridMeteo%easting(1,1)
373 : !write(*,*),"Northing,Easting Last Row, Last Column: ",gridMeteo%northing(192,64),",",gridMeteo%easting(192,64)
374 :
375 0 : gridMeteo%nodata_value = grid%nodata_value
376 0 : gridMeteo%nrows = grid%nrows
377 0 : gridMeteo%ncols = grid%ncols
378 0 : nCell = grid%ncols * grid%nrows
379 0 : deallocate(dem_data, dem_y, dem_x)
380 : ! close netcdf file
381 0 : call ncin%close()
382 :
383 0 : end subroutine ReadDEMNc
384 :
385 : !
386 : !
387 : ! *************************************************************************
388 : !
389 : ! SUBROUTINE READ HEADER GRIDs
390 : !
391 : ! *************************************************************************
392 2 : subroutine readHeader(ic, fName)
393 0 : use mainVar
394 : use mo_kind, only : i4, dp
395 : implicit none
396 : !
397 : character(50) :: dummy
398 : integer(i4), intent(in) :: ic ! input channel
399 : character(256), intent(in) :: fName
400 : !
401 2 : open (unit=ic, file=fName, status='old', action='read')
402 2 : read (ic, *) dummy, grid%ncols ! number of columns
403 2 : read (ic, *) dummy, grid%nrows ! number of rows
404 2 : read (ic, *) dummy, grid%xllcorner ! x coordinate of the lowerleft corner
405 2 : read (ic, *) dummy, grid%yllcorner ! y coordinate of the lowerleft corner
406 2 : read (ic, *) dummy, grid%cellsize ! cellsize x = cellsize y
407 2 : read (ic, *) dummy, grid%nodata_value ! code to define the mask
408 : !
409 2 : if (mod(float(grid%ncols),float(cellFactor)) /= 0.0_dp .or. &
410 : mod(float(grid%nrows),float(cellFactor)) /= 0.0_dp ) then
411 0 : print *, 'nrows or ncols are not exactly divisible by cellFactor ...'
412 0 : stop
413 : end if
414 : !
415 2 : gridMeteo%ncols = ceiling(float(grid%ncols)/float(cellFactor))
416 2 : gridMeteo%nrows = ceiling(float(grid%nrows)/float(cellFactor))
417 : !
418 2 : gridMeteo%cellsize = grid%cellsize * cellFactor
419 2 : gridMeteo%xllcorner = grid%xllcorner
420 : ! to correct the increment that may be produced by ceiling (ensure consistent output in surfer)
421 : gridMeteo%yllcorner = grid%yllcorner + float(grid%nrows) * float(grid%cellsize) &
422 2 : - float(gridMeteo%nrows) * float(gridMeteo%cellsize)
423 2 : gridMeteo%nodata_value = grid%nodata_value
424 2 : thresholdDist = 7.07d-1 * gridMeteo%cellsize
425 2 : end subroutine readHeader
426 :
427 :
428 : ! **************************************************************************
429 : !
430 : ! SUBROUTINE READ METEOROLOGICAL DATABASE
431 : ! :UPDATES Sa.2010-07-06 reading formats / conversion units
432 : !
433 : ! **************************************************************************
434 2 : subroutine ReadDataMeteo
435 2 : use mo_kind, only : i4, dp
436 : use mo_julian, only : NDAYS
437 : use kriging, only : edk_dist
438 : use mainVar
439 : use runControl
440 : implicit none
441 : integer(i4) :: i, ios
442 : integer(i4) :: jDay, jS, jE
443 : integer(i4) :: ds, ms, ys
444 : integer(i4) :: de, me, ye
445 : character(256) :: dummy
446 : character(12) :: dateStart
447 : character(12) :: dateEnd
448 : character(256) :: fileName
449 : integer(i4) :: nStaEfec
450 2 : real(dp) :: p
451 : !
452 : ! allocate vectors for interpolation period
453 2 : jStart = NDAYS (dStart, mStart, yStart)
454 2 : jEnd = NDAYS (dEnd, mEnd, yEnd)
455 :
456 8 : allocate(edk_dist%stat_z(nSta, jStart:jEnd))
457 :
458 2 : nStaEfec = 0
459 40 : do i = 1, nSta
460 114 : allocate (MetSta(i)%z(jStart:jEnd))
461 13908 : MetSta(i)%z = noDataValue
462 : !
463 : ! read yearly data file
464 38 : write (dummy, 2) MetSta(i)%Id
465 38 : fileName = trim(dataPathIn)//trim(dummy)
466 : ! print *, 'read file: '//trim(fileName)
467 38 : open (60, file=fileName, status='old', action='read', iostat=ios)
468 38 : read (60, *) dummy
469 : !
470 : !>>> -------------------------------------------------------------
471 : ! this is temporary / better, robust free format MUST be used
472 : !read(60,1) dummy
473 : !k = scan(dummy, '/', back=.true.)
474 38 : read(60,*) dummy, dateStart , dummy, dateEnd
475 38 : dateStart = dateStart(3:12)
476 38 : dateEnd = dateEnd(3:12)
477 38 : read(dateStart, 3) ys, ms, ds
478 38 : read(dateEnd, 3) ye, me, de
479 : !<<< -------------------------------------------------------------
480 : !
481 : ! set limits
482 38 : jS = NDAYS (ds,ms,ys)
483 38 : jE = NDAYS (de,me,ye)
484 38 : jDay = jS - 1
485 : ! check whether the time series is within the period of interest
486 38 : if ( (jE >= jStart ) .and. (jS <= jEnd ) ) then
487 571186 : do while (.NOT. ios < 0)
488 571154 : read (60, *,iostat=ios) p
489 571154 : jDay = jDay + 1
490 : ! check limits
491 571154 : if (jDay < jStart ) cycle
492 182468 : if (jDay > jEnd ) cycle
493 571154 : MetSta(i)%z(jDay) = p
494 : end do
495 32 : jDay = jDay - 1 ! because iostat is reading an additional line at the end
496 32 : close (60)
497 32 : nStaEfec = nStaEfec + 1
498 : ! unit conversion
499 : !if ( DataConvertFactor /= 1.0_dp ) then
500 11712 : where ( MetSta(i)%z(:) /= noDataValue ) MetSta(i)%z(:) = MetSta(i)%z(:) * DataConvertFactor + OffSet
501 : !end if
502 : ! check consistency
503 32 : if (jDay /= jE ) then
504 0 : print *, 'File ', trim(fileName), ' has missing values.', jE - jDay
505 0 : stop
506 : end if
507 : else
508 6 : close (60)
509 : end if
510 13910 : edk_dist%stat_z(i,:) = MetSta(i)%z ! fill dist container
511 : end do
512 2 : print *, nStaEfec, ' of ', nSta, ' have data from ', yStart , ' to ', yEnd
513 : !
514 : ! formats
515 : 1 format (a256)
516 : 2 format (i5.5,'.dat')
517 : 3 format (i4, 1x, i2, 1x, i2)
518 : !
519 2 : end subroutine ReadDataMeteo
520 :
521 0 : subroutine initMetStaNc
522 :
523 2 : use mo_kind, only: i4, dp
524 :
525 : use mainVar, only: nSta, MetSta, grid, period, noDataValue, &
526 : yStart, mStart, dStart, yEnd, mEnd, dEnd, jStart, jEnd, DataConvertFactor, OffSet
527 : use kriging, only: maxDist, edk_dist
528 : use VarFit, only: hmax
529 : use runControl, only: interMth, DataPathIn
530 : use mo_netCDF, only: NcDataset, NcVariable
531 : use NetCDFVar, only: ncIn_variable_name, ncIn_yCoord_name, ncIn_xCoord_name, ncIn_dem_variable_name
532 : use mo_edk_get_nc_time, only: get_time_vector_and_select
533 :
534 : implicit none
535 : !
536 0 : logical, allocatable :: mask(:, :)
537 : integer(i4) :: iSta, i, j, nrows, ncols
538 : integer(i4) :: n_stations_all ! number of stations in file
539 : integer(i4) :: time_start, time_count
540 : character(256) :: dummy, filename
541 0 : real(dp) :: missing_value
542 0 : real(dp), allocatable :: temp_x(:, :) ! easting coordinate of station
543 0 : real(dp), allocatable :: temp_y(:, :) ! northing coordinate of station
544 0 : real(dp), allocatable :: temp_h(:, :) ! heigth of station
545 0 : real(dp), allocatable :: met_data(:, :, :) ! meteorological data at station
546 : type(period) :: target_period
547 : type(extend) :: domain
548 : type(NcDataset) :: ncin
549 : type(NcVariable) :: ncvar
550 : type(NcVariable) :: ncdem
551 :
552 0 : call getSearchDomain(domain)
553 :
554 : ! setup target period
555 0 : call target_period%init(dStart, mStart, yStart, dEnd, mEnd, yEnd)
556 0 : jStart = target_period%julStart
557 0 : jEnd = target_period%julEnd
558 0 : print *, jStart, jEnd
559 :
560 : ! open netcdf file
561 0 : ncin = NcDataset(dataPathIn, 'r')
562 :
563 0 : ncvar = ncin%getVariable(ncIn_yCoord_name)
564 0 : call ncvar%getData(temp_y)
565 0 : ncvar = ncin%getVariable(ncIn_xCoord_name)
566 0 : call ncvar%getData(temp_x)
567 :
568 : !print *,"Meteo Easting Min: ",minval(temp_x)
569 : !print *,"Meteo Easting Max: ",maxval(temp_x)
570 : !print *,"Meteo Northing Min: ",minval(temp_y)
571 : !print *,"Meteo Northing Max: ",maxval(temp_y)
572 :
573 0 : nrows = size(temp_x, dim=1)
574 0 : ncols = size(temp_x, dim=2)
575 :
576 : ! get time variable
577 0 : ncvar = ncin%getVariable('time')
578 : ! read the time vector and get start index and count of selection
579 0 : call get_time_vector_and_select(ncvar, dataPathIn, -1, time_start, time_count, target_period)
580 : ! ! extract data and select time slice
581 : ! call ncvar%getData(data, start = (/1, 1, time_start/), cnt = (/nRows, nCols, time_count/))
582 :
583 : ! read dem from nc
584 0 : if (interMth .eq. 2) then
585 0 : ncdem = ncin%getVariable(ncIn_dem_variable_name)
586 0 : call ncdem%getData(temp_h, start = (/1, 1/), cnt = (/nRows, nCols/))
587 : end if
588 :
589 : !write(*,*),"DEM Source: ",shape(temp_h)
590 :
591 : ! read meteo data from nc
592 0 : ncvar = ncin%getVariable(ncIn_variable_name)
593 0 : call ncvar%getData(met_data, start = (/1, 1, time_start/), cnt = (/nRows, nCols, time_count/))
594 : !call ncvar%getAttribute('missing_value', missing_value)
595 0 : call ncvar%getAttribute('_FillValue',missing_value)
596 :
597 : !write(*,*),"Data Source: ",shape(met_data)
598 :
599 : ! close netcdf file
600 0 : call ncin%close()
601 :
602 : ! calculate initial mask based on meteorological data
603 0 : mask = (met_data(:, :, 1) .ne. missing_value)
604 :
605 0 : n_stations_all = count(mask)
606 0 : nSta = 0
607 0 : do i = 1, nrows
608 0 : do j = 1, ncols
609 0 : if (.not. mask(i, j)) cycle
610 :
611 : ! check if station is within search distance
612 : ! if yes increase station counting by one
613 : ! otherwise overwrite station in next iteration
614 : !print *,"Meteo Grid Easting: ",temp_x(i,j)
615 : !print *,"DEM Grid Easting LL: ",domain%left
616 : !print *,"DEM Grid Easting TR: ", domain%right
617 0 : if ((temp_x(i, j) .GE. domain%left) .AND. (temp_x(i, j) .LE. domain%right) .AND. &
618 0 : (temp_y(i, j) .GE. domain%bottom) .AND. (temp_y(i, j) .LE. domain%top) ) then
619 0 : nSta = nSta + 1
620 : else
621 0 : mask(i, j) = .False.
622 : end if
623 : end do
624 : end do
625 :
626 0 : print*, ""
627 0 : print*, "Only " , nSta , ' of the given ' , n_stations_all, ' are considered since they are within the search distance.'
628 :
629 : ! >>> populate MetSta object
630 0 : allocate(MetSta(nSta))
631 : ! >>> populate the dynamic distance object
632 0 : edk_dist%nstat = nSta
633 0 : edk_dist%noDataValue = noDataValue
634 0 : allocate(edk_dist%stat_pos(nSta, 2))
635 0 : allocate(edk_dist%stat_z(nSta, target_period%julStart:target_period%julEnd))
636 :
637 0 : iSta = 0
638 0 : do i = 1, nrows
639 0 : do j = 1, ncols
640 0 : if (.not. mask(i,j)) cycle
641 0 : iSta = iSta + 1
642 0 : MetSta(iSta)%Id = iSta
643 0 : MetSta(iSta)%x = temp_x(i, j)
644 0 : MetSta(iSta)%y = temp_y(i, j)
645 0 : MetSta(iSta)%h = temp_h(i, j)
646 : !
647 : ! add meteorological data
648 0 : allocate(MetSta(iSta)%z(target_period%julStart:target_period%julEnd))
649 0 : MetSta(iSta)%z = met_data(i, j, :) * DataConvertFactor + OffSet
650 : !
651 0 : edk_dist%stat_pos(iSta, :) = [MetSta(iSta)%x, MetSta(iSta)%y]
652 0 : edk_dist%stat_z(iSta, :) = MetSta(iSta)%z
653 : end do
654 : end do
655 :
656 0 : deallocate(mask, temp_x, temp_y, met_data) !, temp_h)
657 :
658 0 : end subroutine initMetStaNc
659 :
660 2 : subroutine getSearchDomain(domain)
661 :
662 0 : use mainVar, only: grid,gridMeteo,DEMNcFlag
663 : use kriging, only: maxDist
664 : use VarFit, only: hMax
665 : use runControl, only: flagVario
666 :
667 : implicit none
668 :
669 : type(extend), intent(out) :: domain
670 :
671 2 : real(dp) :: search_distance
672 :
673 : ! take only stations which are in search distance
674 : ! two possibilties for search distance for vario - hMax for EDK - maxDist
675 2 : if (flagVario) then
676 1 : search_distance = hMax
677 : else
678 1 : search_distance = maxDist
679 : end if
680 :
681 : !print *,"DEMNcFlag: ",DEMNcFlag
682 2 : if (DEMNcFlag == 0) then
683 : ! derive borders for station search
684 2 : domain%bottom = -search_distance + grid%yllcorner
685 2 : domain%top = search_distance + grid%yllcorner + grid%nrows * grid%cellsize
686 2 : domain%left = -search_distance + grid%xllcorner
687 2 : domain%right = search_distance + grid%xllcorner + grid%ncols * grid%cellsize
688 : else
689 0 : domain%bottom = -search_distance + minval(gridMeteo%northing)
690 0 : domain%top = search_distance + maxval(gridMeteo%northing)
691 0 : domain%left = -search_distance + minval(gridMeteo%easting)
692 0 : domain%right = search_distance + maxval(gridMeteo%easting)
693 : end if
694 :
695 : !print *,"DEM Min Northing: ",minval(gridMeteo%northing)
696 : !print *,"DEM Max Northing: ",maxval(gridMeteo%northing)
697 : !print *,"DEM Min Easting: ",minval(gridMeteo%easting)
698 : !print *,"DEM Max Easting: ",maxval(gridMeteo%easting)
699 :
700 2 : end subroutine getSearchDomain
701 :
702 :
703 : ! *************************************************************************
704 : !
705 : ! SUBROUTINE WRITE HEADER GRIDs
706 : !
707 : ! *************************************************************************
708 2 : subroutine writeHeader(ic, fName)
709 2 : use mainVar
710 : use mo_kind, only : i4
711 : implicit none
712 : integer(i4), intent(in) :: ic ! input channel
713 : character(256), intent(in) :: fName
714 : !
715 2 : open (unit=ic, file=fName, status='unknown')
716 2 : write (ic, 1) 'ncols ', gridMeteo%ncols
717 2 : write (ic, 1) 'nrows ', gridMeteo%nrows
718 2 : write (ic, 2) 'xllcorner ', gridMeteo%xllcorner
719 2 : write (ic, 2) 'yllcorner ', gridMeteo%yllcorner
720 2 : write (ic, 1) 'cellsize ', gridMeteo%cellsize
721 2 : write (ic, 1) 'NODATA_value', gridMeteo%nodata_value
722 2 : close (ic)
723 : ! formats
724 : 1 format (a12, 2x, i10)
725 : 2 format (a12, 2x, f10.1)
726 : !
727 2 : end subroutine writeHeader
728 2 : end module mo_edk_read_data
|
