Line data Source code
1 : !> \file mo_smi_read.f90
2 : !> \brief \copybrief mo_smi_read
3 : !> \details \copydetails mo_smi_read
4 :
5 : !> \brief Reading routines for SMI
6 : !> \author Luis Samaniego
7 : !> \date 09.02.2011
8 : !> \author Stephan Thober
9 : !> \date 07.11.2017
10 : !! - switched to mo_netcdf
11 : ! - added invert_SMI
12 : !> \copyright Copyright 2005-\today, the CHS Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
13 : !! SMI is released under the LGPLv3+ license \license_note
14 : module mo_smi_read
15 :
16 : USE mo_kind, only : i4, sp, dp
17 :
18 : IMPLICIT NONE
19 :
20 : PRIVATE
21 :
22 : PUBLIC :: ReadDataMain ! read everything
23 :
24 : ! ------------------------------------------------------------------
25 :
26 : CONTAINS
27 :
28 : ! ------------------------------------------------------------------
29 :
30 :
31 : !> \brief Reads main file
32 : !> \author Luis E. Samaniego-Eguiguren, UFZ
33 : !> \date 23.05.2007
34 : !> \note packed fields are stored in dim1->dim2 sequence
35 6 : subroutine ReadDataMain( SMI_in, do_cluster, ext_smi, invert_SMI, read_opt_h, silverman_h, opt_h, lats_1d, lons_1d,&
36 : lats_2d, lons_2d,easting, northing, basin_flag, mask, SM_kde, SM_eval, &
37 : Basin_Id, SMI_thld, outpath, cellsize, thCellClus, nCellInter, do_sad, deltaArea, &
38 : nCalendarStepsYear, per_kde, per_eval, per_smi )
39 :
40 : use mo_kind, only: i4
41 : use mo_message, only: message
42 : use mo_utils, only: notequal, equal
43 : use mo_netcdf, only: NcDataset, NcVariable
44 : use mo_constants, only: nodata_dp
45 : use mo_smi_global_variables, only: period
46 : use mo_smi_info, only: version, version_date, file_namelist
47 : use mo_os, only: check_path_isfile, path_isfile
48 :
49 : implicit none
50 :
51 : ! input / output Variables
52 : real(sp), dimension(:,:), allocatable, intent(out) :: SMI_in !< SMI only read if ext_smi is TRUE
53 : logical, intent(out) :: do_sad !< do SAD analysis
54 : logical, intent(out) :: do_cluster !< do cluster calculation
55 : logical, intent(out) :: ext_smi !< clsutering external data
56 : logical, intent(out) :: invert_SMI !< do inversion of SMI
57 : logical, intent(out) :: read_opt_h !< read kernel width
58 : logical, intent(out) :: silverman_h !< optimize kernel width
59 : logical, intent(out) :: basin_flag !< basin flag
60 : logical, dimension(:,:), allocatable, intent(out) :: mask !< grid mask
61 : integer(i4), intent(out) :: thCellClus !< treshold for cluster formation in space ~ 640 km2
62 : integer(i4), intent(out) :: nCellInter !< number cells for joining clusters in time ~ 6400 km2
63 : integer(i4), intent(out) :: deltaArea !< number of cells per area interval
64 : !> number of calendar time steps per year (month=12, day=365)
65 : integer(i4), intent(out) :: nCalendarStepsYear
66 : integer(i4), dimension(:,:), allocatable, intent(out) :: Basin_Id !< IDs for basinwise drought analysis
67 : real(sp), intent(out) :: cellsize !< cell edge lenght of input data
68 : real(sp), dimension(:,:), allocatable, intent(out) :: SM_kde !< daily / monthly fields packed for estimation
69 : real(sp), dimension(:,:), allocatable, intent(out) :: SM_eval !< monthly fields packed for evaluation
70 : real(dp), dimension(:,:), allocatable, intent(out) :: opt_h !< optimized kernel width
71 : real(dp), dimension(:), allocatable, intent(out) :: lats_1d, lons_1d !< latitude and longitude vectors of input
72 : real(dp), dimension(:,:), allocatable, intent(out) :: lats_2d !< latitude vectors of input
73 : real(dp), dimension(:,:), allocatable, intent(out) :: lons_2d !< longitude vectors of input
74 : real(dp), dimension(:), allocatable, intent(out) :: easting !< easting coordinates of input
75 : real(dp), dimension(:), allocatable, intent(out) :: northing !< easting coordinates of input
76 : real(sp), intent(out) :: SMI_thld !< SMI threshold for clustering
77 : character(len=256), intent(out) :: outpath !< ouutput path for results
78 : type(period), intent(out) :: per_kde !< period contain information during estimation
79 : type(period), intent(out) :: per_eval !< period during evaluation
80 : type(period), intent(out) :: per_smi !< period of smi file
81 :
82 : ! local Variables
83 : logical :: read_sm
84 : integer(i4) :: ii
85 : integer(i4) :: nCells ! number of effective cells
86 : integer(i4) :: lag! number of lag days for daily files (0,7,31)
87 : integer(i4) :: nTimeSteps ! number of time steps excluding leap days
88 :
89 : ! directories, filenames, attributes
90 : character(256) :: maskfName
91 : character(256) :: opt_h_file
92 : character(256) :: basinfName
93 :
94 : ! variable names in netcdf input files
95 : character(256) :: soilmoist_file
96 : character(256) :: mask_vname
97 : character(256) :: SM_vname
98 : character(256) :: ext_smi_file
99 : character(256) :: basin_vname
100 : character(256) :: opt_h_vname
101 : character(256) :: opt_h_sm_vname
102 : real(sp) :: nodata_value ! local data nodata value in nc for mask creation
103 :
104 6 : real(sp), dimension(:,:), allocatable :: dummy_D2_sp
105 6 : real(sp), dimension(:,:,:), allocatable :: dummy_D3_sp
106 6 : real(dp), dimension(:,:,:), allocatable :: dummy_D3_dp
107 :
108 : type(NcDataset) :: nc_in
109 : type(NcVariable) :: nc_var
110 :
111 : ! read main config
112 : namelist/mainconfig/basin_flag, basinfName, basin_vname, maskfName, mask_vname, &
113 : soilmoist_file, SM_vname, outpath, nCalendarStepsYear, lag, &
114 : silverman_h, read_opt_h, opt_h_vname, opt_h_sm_vname, opt_h_file, &
115 : SMI_thld, do_cluster, ext_smi, ext_smi_file, cellsize, thCellClus, nCellInter, &
116 : do_sad, deltaArea, invert_SMI
117 :
118 :
119 6 : call message("")
120 6 : call message("====================================================")
121 6 : call message("||")
122 6 : call message("|| SOIL MOISTURE INDEX")
123 6 : call message("|| VERSION ", version, " (", version_date, ")")
124 6 : call message("||")
125 6 : call message("====================================================")
126 6 : call message("")
127 :
128 : ! dummy init to avoid gnu compiler complaint
129 6 : outpath = '-999'; cellsize=-999.0_sp; thCellClus=-999; nCellInter=-999; deltaArea=-999
130 :
131 6 : do_cluster = .FALSE.
132 6 : do_sad = .FALSE.
133 6 : SMI_thld = 0.0_dp
134 6 : silverman_h = .FALSE.
135 :
136 : ! read namelist
137 6 : call check_path_isfile(path=file_namelist, raise=.true.)
138 6 : open (unit=10, file=file_namelist, status='old')
139 6 : read(10, nml=mainconfig)
140 6 : close (10)
141 6 : call message(file_namelist, ' read ...ok')
142 :
143 : ! consistency check
144 6 : if ( .not.(( nCalendarStepsYear == 12) .or. ( nCalendarStepsYear == 365) ) ) &
145 0 : stop '***ERROR: number of time steps per year different from 12 or 365'
146 :
147 6 : if ( ( lag .gt. 0 ) .and. ( nCalendarStepsYear == 12) ) &
148 0 : stop '***ERROR: lag = 0 for montly SM values'
149 :
150 6 : if ( do_sad .and. (.not. do_cluster) ) &
151 0 : stop '***ERROR: flag do_cluster must be set to .TRUE. to perform SAD analisys'
152 :
153 6 : if (invert_SMI .and. (.not. ext_smi)) &
154 0 : stop '***ERROR: if invert_SMI is .TRUE., then ext_smi must be .TRUE. and ext_smi_file must be given.'
155 :
156 : ! read sm if not external smi or invert is given
157 6 : read_sm = (.not. ext_smi)
158 :
159 : ! read main mask
160 9 : if (path_isfile(maskfName)) then
161 3 : nc_in = NcDataset(maskfName, "r")
162 3 : nc_var = nc_in%getVariable(trim(mask_vname))
163 3 : call nc_var%getData(dummy_D2_sp)
164 3 : call nc_var%getAttribute('missing_value', nodata_value)
165 3 : call read_latlon( nc_in, lats_1d, lons_1d, lats_2d, lons_2d,easting, northing)
166 3 : call nc_in%close()
167 :
168 : ! create mask
169 126 : mask = merge( .true., .false., notequal(dummy_D2_sp, nodata_value ) )
170 3 : deallocate( dummy_D2_sp )
171 3 : nCells = n_cells_consistency(mask)
172 3 : call message('mask read ...ok')
173 : else
174 3 : call message('Mask file ' // trim(maskfName) // ' does not exist!')
175 : end if
176 :
177 : ! read basin mask
178 6 : if ( basin_flag .AND. (.NOT. ext_smi)) then
179 0 : nc_in = NcDataset(basinfName, "r")
180 0 : nc_var = nc_in%getVariable(trim(basin_vname))
181 0 : call nc_var%getData(Basin_id)
182 0 : call nc_var%getAttribute('missing_value', nodata_value)
183 0 : call nc_in%close()
184 :
185 : ! consistency check
186 0 : if ( ( size( Basin_Id, 1) .ne. size( mask, 1 ) ) .or. &
187 : ( size( Basin_Id, 2) .ne. size( mask, 2 ) ) ) then
188 0 : call message('***ERROR: size mismatch between basin field and given mask file')
189 0 : stop 1
190 : end if
191 : ! intersect mask and basin_id
192 0 : mask = merge( .true., .false., mask .and. ( Basin_Id .ne. int( nodata_value, i4) ) )
193 0 : Basin_Id = merge( Basin_Id, int( nodata_value, i4), mask )
194 0 : call message('basin ID read ...ok')
195 : end if
196 :
197 : ! *************************************************************************
198 : ! >>> read SM field
199 : ! *************************************************************************
200 6 : if ( read_sm ) then
201 4 : nc_in = NcDataset(soilmoist_file, "r")
202 4 : nc_var = nc_in%getVariable(trim(SM_vname))
203 4 : call nc_var%getData(dummy_D3_sp)
204 4 : if (.not. allocated(mask)) then
205 1 : call nc_var%getAttribute('missing_value', nodata_value)
206 32 : mask = (dummy_D3_sp(:, :, 1) .ne. nodata_value)
207 1 : nCells = n_cells_consistency(mask)
208 1 : call message('mask read ...ok')
209 : end if
210 :
211 : ! consistency check
212 4 : if ( ( size( dummy_D3_sp, 1) .ne. size( mask, 1 ) ) .or. &
213 : ( size( dummy_D3_sp, 2) .ne. size( mask, 2 ) ) ) then
214 0 : call message('***ERROR: size mismatch between SM field and given mask file')
215 0 : call nc_in%close()
216 0 : stop 1
217 : end if
218 :
219 : ! find number of leap years in SM data set
220 : ! get timepoints in days and masks for climatologies at calendar day or month
221 4 : call get_time(nc_in, size( dummy_D3_sp, 3 ), nCalendarStepsYear, per_eval)
222 :
223 : ! read lats and lon from file
224 4 : call read_latlon( nc_in, lats_1d, lons_1d, lats_2d, lons_2d,easting, northing)
225 4 : call nc_in%close()
226 :
227 :
228 : ! if ( nCalendarStepsYear .eq. YearMonths ) then
229 16 : allocate( SM_eval( nCells, size( dummy_D3_sp, 3 ) ) )
230 : ! no lag average, use monthly data as read
231 1029 : do ii = 1, size( dummy_D3_sp, 3 )
232 39029 : SM_eval(:,ii) = pack( real(dummy_D3_sp(:,:,ii),sp), mask )
233 : end do
234 :
235 4 : deallocate( dummy_D3_sp )
236 :
237 : ! calculate moving-lag average
238 8 : if (lag .gt. 0) then
239 1 : call message('calculate moving average')
240 1 : nTimeSteps = size(SM_eval, dim=2)
241 4 : allocate(dummy_d2_sp(size(SM_eval, dim=1), size(SM_eval, dim=2)))
242 3902 : dummy_d2_sp = SM_eval
243 301 : do ii = 1, nTimeSteps
244 301 : if ( ii .lt. lag ) then
245 2327 : SM_eval(:,ii) = dummy_d2_sp(:,ii)
246 : else
247 262933 : SM_eval(:,ii) = sum(dummy_D2_sp(:,ii-lag+1:ii), DIM=2) / real(lag,sp)
248 : end if
249 : end do
250 1 : deallocate(dummy_d2_sp)
251 : end if
252 :
253 : end if
254 :
255 : ! *************************************************************************
256 : ! >>> read_opt_h
257 : ! *************************************************************************
258 6 : if (allocated(mask)) then
259 16 : allocate( opt_h(ncells, nCalendarStepsYear))
260 5577 : opt_h = nodata_dp
261 : end if
262 : !
263 6 : if ( read_opt_h ) then
264 :
265 4 : nc_in = NcDataset(trim( opt_h_file ), 'r')
266 :
267 4 : nc_var = nc_in%getVariable(trim(opt_h_sm_vname))
268 4 : call nc_var%getData(dummy_D3_sp)
269 : ! read mask if not done already
270 4 : if (.not. allocated(mask)) then
271 1 : call nc_var%getAttribute('missing_value', nodata_value)
272 42 : mask = (dummy_D3_sp(:, :, 1) .ne. nodata_value)
273 1 : nCells = n_cells_consistency(mask)
274 4 : allocate( opt_h(ncells, nCalendarStepsYear))
275 277 : opt_h = nodata_dp
276 1 : call message('mask read ...ok')
277 : end if
278 : ! unpack values
279 16 : allocate(SM_kde(nCells, size(dummy_D3_sp, 3)))
280 35417 : do ii = 1, size( dummy_D3_sp, 3 )
281 35417 : SM_kde(:, ii) = pack( dummy_D3_sp(:,:,ii), mask)
282 : end do
283 4 : deallocate( dummy_D3_sp )
284 4 : call message('read kde soil moisture values from file... ok')
285 :
286 : ! read optimized kernel width from file
287 4 : nc_var = nc_in%getVariable(trim(opt_h_vname))
288 4 : call nc_var%getData(dummy_D3_dp)
289 405 : do ii = 1, size( dummy_D3_dp, 3 )
290 17967 : opt_h( :, ii ) = pack( real( dummy_D3_dp(:,:,ii),sp ), mask )
291 : end do
292 4 : deallocate( dummy_D3_dp )
293 5577 : if ( any( equal( opt_h, nodata_dp ) ) ) &
294 : ! if ( any( opt_h .eq. nodata_dp ) ) &
295 0 : stop '***ERROR kernel width contains nodata values'
296 4 : call message('read kde kernel width from file... ok')
297 :
298 :
299 : ! get timepoints in days and mask of months
300 4 : call get_time(nc_in, size(SM_kde, 2), nCalendarStepsYear, per_kde)
301 : ! close file
302 4 : call nc_in%close()
303 :
304 : else
305 : ! SM_eval is SM_kde too
306 2 : if (allocated(SM_eval)) then
307 8284 : SM_kde = SM_eval
308 1 : per_kde = per_eval
309 : end if
310 : end if
311 :
312 6 : if (allocated(SM_kde)) then
313 5 : if (((nCalendarStepsYear .eq. 12) .and. (modulo(size( SM_kde, 2 ), nCalendarStepsYear) .ne. 0_i4)) .or. &
314 5 : ((nCalendarStepsYear .ne. 12) .and. (modulo(size( SM_kde, 2 ), nCalendarStepsYear) .ne. per_kde%n_leap_days))) then
315 : #ifdef SMIDEBUG
316 : print *, modulo(size( SM_kde, 2 ), nCalendarStepsYear), per_kde%n_leap_days
317 : #endif
318 0 : print *, '***ERROR: timesteps in provided SM_kde field must be multiple of nCalendarStepsYear (', &
319 0 : nCalendarStepsYear, '; leap days are accounted for)'
320 0 : stop
321 : end if
322 : end if
323 :
324 :
325 : ! *************************************************************************
326 : ! >>> read external SMI field
327 : ! *************************************************************************
328 6 : if (ext_smi) then
329 2 : nc_in = NcDataset(trim(ext_smi_file), 'r')
330 2 : nc_var = nc_in%getVariable('SMI')
331 2 : call nc_var%getData(dummy_D3_sp)
332 2 : if (.not. allocated(mask)) then
333 1 : call nc_var%getAttribute('missing_value', nodata_value)
334 42 : mask = (dummy_D3_sp(:, :, 1) .ne. nodata_value)
335 1 : nCells = n_cells_consistency(mask)
336 1 : call message('mask read ...ok')
337 : end if
338 : ! consistency check
339 2 : if ( ( size( dummy_D3_sp, 1) .ne. size( mask, 1 ) ) .or. &
340 : ( size( dummy_D3_sp, 2) .ne. size( mask, 2 ) ) ) then
341 0 : call message('***ERROR: size mismatch between SMI field and given mask file')
342 0 : stop
343 : end if
344 :
345 8 : allocate(SMI_in( nCells, size( dummy_D3_sp, 3 ) ) )
346 367 : do ii = 1, size( dummy_D3_sp, 3 )
347 367 : SMI_in(:,ii) = pack(dummy_D3_sp(:,:,ii), mask)
348 : end do
349 2 : deallocate( dummy_D3_sp )
350 :
351 : ! get timepoints in days and mask of months
352 2 : call get_time( nc_in, size(SMI_in, dim=2), nCalendarStepsYear, per_smi)
353 :
354 2 : call read_latlon( nc_in, lats_1d, lons_1d, lats_2d, lons_2d,easting, northing)
355 2 : call nc_in%close()
356 : else
357 4 : per_smi = per_kde
358 : end if
359 :
360 6 : end subroutine ReadDataMain
361 :
362 :
363 : !> \brief Convert input time into months & check timesteps & determine mask for months
364 : !> \author Matthias Zink, Oct 2012
365 : !> \author Stephan Thober
366 : !> \date Nov 2017
367 : !! - convert to mo_netcdf
368 : !> \author Stephan Thober
369 : !> \date Aug 2019
370 : !! - added type period
371 0 : subroutine get_time(nc_in, in_time_steps, nCalendarStepsYear, per_out) !, mask)
372 : !
373 6 : use mo_julian, only: date2dec, dec2date
374 : use mo_message, only: message
375 : use mo_string_utils, only: DIVIDE_STRING
376 : use mo_netcdf, only: NcDataset, NcVariable
377 : use mo_constants, only: YearMonths, DayHours
378 : use mo_smi_global_variables, only: period, period_init
379 :
380 : implicit none
381 :
382 : type(NcDataset), intent(in) :: nc_in ! NetCDF dataset
383 : integer(i4), intent(in) :: in_time_steps
384 : integer(i4), intent(in) :: nCalendarStepsYear
385 : type(period), intent(out) :: per_out
386 : integer(i4) :: yStart ! start year of the dataser
387 : integer(i4) :: mStart ! start month of the dataser
388 : integer(i4) :: dStart ! start day of the dataser
389 : integer(i4) :: yEnd ! end year of the dataset
390 : integer(i4) :: mEnd ! end month of the dataset
391 : integer(i4) :: dEnd ! end day of the dataset
392 : integer(i4) :: nTimeSteps ! size of the time dimension
393 10 : integer(i4), dimension(:), allocatable :: timepoints
394 : ! timestep in days or months
395 :
396 : type(ncVariable) :: nc_var
397 : integer(i4) :: i ! loop variable
398 : integer(i4) :: yRef, dRef, mRef ! reference time of NetCDF (unit attribute of
399 : integer(i4) :: month, year, d !
400 : !
401 10 : integer(i4), dimension(:), allocatable :: timesteps ! time variable of NetCDF in input units
402 : !
403 : character(256) :: AttValues ! netcdf attribute values
404 10 : character(256), dimension(:), allocatable :: strArr ! dummy for netcdf attribute handling
405 10 : character(256), dimension(:), allocatable :: date ! dummy for netcdf attribute handling
406 : real(dp) :: ref_jday ! refernce date of dateset as julian day
407 : !
408 : ! get unit attribute of variable 'time'
409 10 : nc_var = nc_in%getVariable('time')
410 10 : call nc_var%getAttribute('units', AttValues)
411 : ! AttValues looks like "<unit> since YYYY-MM-DD HH:MM:SS"
412 10 : call DIVIDE_STRING(trim(AttValues), ' ', strArr)
413 : !
414 10 : call DIVIDE_STRING(trim(strArr(3)), '-', date)
415 10 : read(date(1),*) yRef
416 10 : read(date(2),*) mRef
417 10 : read(date(3),*) dRef
418 :
419 : ! allocate and initialize
420 10 : call nc_var%getData(timesteps)
421 10 : nTimeSteps = size(timesteps, dim=1)
422 30 : allocate(timepoints ( nTimeSteps ))
423 36813 : timepoints = 0
424 :
425 :
426 : ! strArr(1) is <unit>
427 10 : ref_jday = date2dec(dd=dRef, mm=mRef, yy=yRef)
428 36813 : do i = 1, nTimeSteps
429 38908 : select case (strArr(1))
430 : case('hours')
431 2105 : call dec2date(real(timesteps(i), dp)/DayHours + ref_jday, dd=d, mm=month, yy=year)
432 2105 : timepoints(i) = timepoints(i) + nint( real(timesteps(i) - timesteps(1), dp) / DayHours , i4)
433 : case('days')
434 34698 : call dec2date(real(timesteps(i), dp) + ref_jday, dd=d, mm=month, yy=year)
435 34698 : timepoints(i) = timepoints(i) + timesteps(i) - timesteps(1)
436 : case('months')
437 0 : d = dRef ! set to dref as default
438 0 : month = mod( (timesteps(i) + mRef ), int(YearMonths, i4) )
439 0 : year = yRef + floor(real( timesteps(i) + mRef, dp) / YearMonths )
440 : ! correct month and year for december
441 0 : if (month .eq. 0 ) then
442 0 : month = 12
443 0 : year = year - 1
444 : end if
445 : !
446 0 : if (i .EQ. 1) then
447 0 : timepoints(i) = 0
448 : else
449 0 : timepoints(i) = timepoints(i) + nint(date2dec(dd=15 ,mm=month, yy=year) - date2dec(dd=15 ,mm=mStart, yy=yStart), i4)
450 : end if
451 : case DEFAULT
452 0 : call message('***ERROR: Time slicing in ' , trim(strArr(1)), ' not implemented!')
453 36803 : stop
454 : end select
455 :
456 : ! save start year and month output timestamp
457 36803 : if (i .EQ. 1) then
458 10 : yStart = year
459 10 : mStart = month
460 10 : dStart = d
461 : end if
462 36813 : if (i .EQ. nTimeSteps) then
463 10 : yEnd = year
464 10 : mEnd = month
465 10 : dEnd = d
466 : end if
467 : end do
468 :
469 10 : per_out = period_init(yStart, mStart, dStart, yEnd, mEnd, dEnd, timepoints, strArr(1), nCalendarStepsYear)
470 :
471 : ! consistency check
472 10 : if (in_time_steps .ne. size(per_out%time_points, dim=1)) then
473 0 : call message('***ERROR: time steps of array and time axis are different by more than leap days')
474 0 : stop 1
475 : end if
476 :
477 10 : end subroutine get_time
478 :
479 6 : integer(i4) function n_cells_consistency(mask)
480 :
481 10 : use mo_message, only: message
482 :
483 : implicit none
484 :
485 : logical, intent(in) :: mask(:, :)
486 :
487 236 : n_cells_consistency = count( mask )
488 6 : if ( n_cells_consistency .eq. 0 ) then
489 0 : call message('***ERROR: no cell selected in mask')
490 0 : stop 1
491 : end if
492 6 : end function n_cells_consistency
493 :
494 9 : subroutine read_latlon(nc_in, lats_1d, lons_1d, lats_2d, lons_2d,easting, northing)
495 :
496 : use mo_kind, only: dp,i4
497 : use mo_netcdf, only: NcDataset, NcVariable
498 :
499 : implicit none
500 :
501 : type(NcDataset), intent(in) :: nc_in ! NetCDF dataset
502 : real(dp), dimension(:), allocatable, intent(out) :: lats_1d, lons_1d ! latitude and longitude vectors of input
503 : real(dp), dimension(:,:), allocatable, intent(out) :: lats_2d, lons_2d ! latitude and longitude vectors of input
504 : real(dp), dimension(:), allocatable, intent(out) :: easting ! easting coordinates of input
505 : real(dp), dimension(:), allocatable, intent(out) :: northing ! easting coordinates of input
506 : integer(i4) :: lat_lon_Ndim
507 : type(NcVariable) :: nc_var
508 :
509 : ! read lats and lon from file
510 9 : if (nc_in%hasVariable('lat')) then
511 8 : nc_var = nc_in%getVariable('lat')
512 8 : lat_lon_Ndim=nc_var%getNoDimensions()
513 8 : if (lat_lon_Ndim == 1) then
514 0 : call nc_var%getData(lats_1d)
515 8 : else if (lat_lon_Ndim == 2) then
516 8 : call nc_var%getData(lats_2d)
517 : ! nc_var = nc_in%getVariable('northing')
518 : ! call nc_var%getData(northing)
519 : end if
520 : end if
521 9 : if (nc_in%hasVariable('lon')) then
522 8 : nc_var = nc_in%getVariable('lon')
523 8 : lat_lon_Ndim=nc_var%getNoDimensions()
524 8 : if (lat_lon_Ndim == 1) then
525 0 : call nc_var%getData(lons_1d)
526 8 : else if (lat_lon_Ndim == 2) then
527 8 : call nc_var%getData(lons_2d)
528 : ! nc_var = nc_in%getVariable('easting')
529 : ! call nc_var%getData(easting)
530 : end if
531 : end if
532 :
533 9 : end subroutine read_latlon
534 :
535 : end module mo_smi_read
|
