Soil Water Balance (SWB2)
Loading...
Searching...
No Matches
data_catalog_entry.F90
Go to the documentation of this file.
1!
2! concept: for each significant gridded data input, keep track of the
3! native coordinates, the transformed base (project) coordinates,
4! and provide methods for extracting data from appropriate locations
5! as needed.
6!
7
9
11 use datetime
12 use exceptions
14 use logfiles
15 use fstring
16 use fstring_list
17 use grid
19 use parameters
20 use iso_c_binding
21 implicit none
22
23 private
24
25 integer (c_int), public, parameter :: netcdf_file_open = 27
26 integer (c_int), public, parameter :: netcdf_file_closed = 42
27
28 integer (c_int), parameter, public :: file_template_capitalized_monthname = 0
29 integer (c_int), parameter, public :: file_template_lowercase_monthname = 1
30 integer (c_int), parameter, public :: file_template_uppercase_monthname = 2
31
32 integer (c_int), parameter :: constant_grid = 0
33 integer (c_int), parameter :: static_grid = 1
34 integer (c_int), parameter :: static_netcdf_grid = 2
35 integer (c_int), public, parameter :: dynamic_grid = 3
36 integer (c_int), parameter :: dynamic_netcdf_grid = 4
37 integer (c_int), parameter :: no_grid = 5
38 integer (c_int), parameter :: table_lookup = 6
39
40 integer (c_int), parameter :: filetype_arc_ascii = 0
41 integer (c_int), parameter :: filetype_surfer = 1
42 integer (c_int), parameter :: filetype_netcdf = 2
43 integer (c_int), parameter :: filetype_ascii_table = 3
44 integer (c_int), parameter :: filetype_none = 4
45
46 type, public :: data_catalog_entry_t
47 character (len=:), allocatable :: skeyword
48 type (data_catalog_entry_t), pointer :: previous => null()
49 type (data_catalog_entry_t), pointer :: next => null()
50
51 integer (c_int) :: isourcedataform = no_grid ! constant, static grid, dynamic grid
52 integer (c_int) :: isourcedatatype = datatype_na ! real, short, integer, etc.
53 integer (c_int) :: isourcefiletype = filetype_none ! Arc ASCII, Surfer, NetCDF
54 integer (c_int) :: itargetdatatype = datatype_na ! Fortran real, integer, etc.
55
56 character (len=256) :: sdescription = ""
57 character (len=256) :: ssourceproj4_string = ""
58 character (len=256) :: stargetproj4_string = ""
59 character (len=256) :: ssourcefiletype = ""
60 character (len=512) :: ssourcefilename = ""
61 character (len=512) :: sfilenametemplate = ""
62 integer (c_int) :: ifilename_monthname_capitalization_rule = file_template_capitalized_monthname
63 character (len=512) :: soldfilename = ""
64 character (len=256) :: sdatecolumnname = ""
65 character (len=10) :: sdefaultdateformat = "YYYY-MM-DD"
66 character (len=256) :: svaluecolumnname = ""
67
68 real (c_float), allocatable :: table_values_real(:)
69 type (datetime_t), allocatable :: table_dates(:)
70 integer (c_int) :: table_indx
71
72 logical (c_bool) :: ltablevalueshavebeenretrieved = false
73
74 integer (c_int) :: ifilecount = -1
75 integer (c_int) :: ifilecountyear = -9999
76 real (c_float) :: rminallowedvalue = -rbigval ! default condition is to impose
77 real (c_float) :: rmaxallowedvalue = rbigval ! no bounds on data
78 integer (c_int) :: iminallowedvalue = -ibigval ! default condition is to impose
79 integer (c_int) :: imaxallowedvalue = ibigval ! no bounds on data
80 real (c_float) :: rmissingvaluescode = -rbigval
81 integer (c_int) :: imissingvaluescode = -ibigval
82 character (len=2) :: smissingvaluesoperator = "&&"
83 integer (c_int) :: imissingvaluesaction = 0
84
85 real (c_double) :: duserscalefactor = 1_c_double
86 real (c_double) :: duseraddoffset = 0_c_double
87 real (c_double) :: dusersuboffset = 0_c_double
88 real (c_double) :: rx_coord_addoffset = 0.0_c_double
89 real (c_double) :: ry_coord_addoffset = 0.0_c_double
90 real (c_double) :: rcoordinatetolerance = 0.0_c_double
91
92 real (c_double), allocatable :: rx_coordinate_subset(:)
93 real (c_double), allocatable :: ry_coordinate_subset(:)
94
95 logical (c_bool) :: lallowmissingfiles = false
96 logical (c_bool) :: lallowautomaticdataflipping = true
97 logical (c_bool) :: lfliphorizontal = false
98 logical (c_bool) :: lflipvertical = false
99 logical (c_bool) :: lusemajorityfilter = false
100 logical (c_bool) :: lrequirecompletespatialcoverage = true
101
102 integer (c_int) :: idaystopadatyearsend = 0
103 integer (c_int) :: idaystopadifleapyear = 1
104 integer (c_int) :: istartyear = -9999
105 integer (c_int) :: iendyear = -9999
106 logical (c_bool) :: lpadreplacewithzero = false
107 logical (c_bool) :: lpadvalues = false
108 logical (c_bool) :: lisannualgrid = false
109
110 ! the following are only used if data are being read from a NetCDF file
111 character (len=32) :: svariablename_x = "x"
112 character (len=32) :: svariablename_y = "y"
113 character (len=32) :: svariablename_z = ""
114 character (len=32) :: svariablename_time = "time"
115 character (len=8) :: svariableorder = "tyx"
116
117 type (grid_bounds_t) :: grid_bounds_native
118 type (grid_bounds_t) :: grid_bounds_base
119
120 integer (c_int) :: inc_file_status = netcdf_file_closed
121 type (t_netcdf4_file) :: ncfile
122
123 integer (c_int) :: inc_archive_status = netcdf_file_closed
124 type (t_netcdf4_file) :: ncfile_archive
125 integer (c_size_t) :: incfile_recnum = 0
126
127 integer (c_int) :: iconstantvalue = 0
128 real (c_float) :: rconstantvalue = 0.0
129
130 ! pGrdNative is a grid created to serve as an intermediary between
131 ! the native coordinate of the data source file and the project coordinates
132 ! in use by swb
133 type (general_grid_t), pointer :: pgrdnative => null()
134 logical (c_bool) :: lgridispersistent = false
135 logical (c_bool) :: lgridhaschanged = false
136 logical (c_bool) :: lperformfullinitialization = true
137 logical (c_bool) :: lcreatelocalnetcdfarchive = false
138
139 ! pGrdBase takes the coordinate system and dimensions as defined
140 ! for the overall SWB project (i.e. BASE_PROJECTION_DEFINITION )
141 type (general_grid_t), pointer :: pgrdbase => null()
142
143 contains
144
145 procedure :: setkey => set_keyword_sub
146
151 generic :: initialize => initialize_constant_int_data_object_sub, &
155
156 procedure :: initialize_netcdf => initialize_netcdf_data_object_sub
157
158 procedure :: set_scale => set_scale_sub
159 procedure :: set_sub_offset => set_sub_offset_sub
160 procedure :: set_add_offset => set_add_offset_sub
161
162 procedure :: set_x_offset => set_x_coord_offset_sub
163 procedure :: set_y_offset => set_y_coord_offset_sub
164 procedure :: set_coordinate_tolerance => set_coordinate_tolerance_sub
165
166 procedure :: set_majority_filter_flag => set_majority_filter_flag_sub
167
170 generic :: set_valid_minimum => set_minimum_allowable_value_int_sub, &
172
175 generic :: set_valid_maximum => set_maximum_allowable_value_int_sub, &
177
178 procedure :: set_grid_flip_horizontal => set_grid_flip_horizontal_sub
179 procedure :: set_grid_flip_vertical => set_grid_flip_vertical_sub
180 procedure :: allow_missing_files => set_allow_missing_files_flag_sub
181 procedure :: do_not_allow_netcdf_grid_data_flipping &
183
184 procedure :: getvalues_constant => getvalues_constant_sub
185 procedure :: getvalues_gridded => getvalues_gridded_sub
186
189 generic :: getvalues_netcdf => getvalues_dynamic_netcdf_sub, &
191
192 procedure :: getvalues => getvalues_sub
193
194 procedure :: get_value_int_sub
196 generic :: getvalue => get_value_int_sub, &
198
199 ! procedure :: update => update_data_object_sub
200 ! procedure :: destroy => create_data_object_sub
201 procedure :: get_filetype => get_source_filetype_fn
202
207
210 generic :: set_constant => set_constant_value_int, &
212
213 procedure :: make_filename => make_filename_from_template
214 procedure :: set_target_proj4 => set_target_proj4_string_sub
215 procedure :: set_source_proj4 => set_source_proj4_string_sub
216 procedure :: set_variable_order => set_variable_order_sub
217 procedure :: set_complete_spatial_coverage_flag => set_complete_spatial_coverage_flag_sub
218 procedure :: dump_data_structure => dump_data_structure_sub
219 procedure :: set_make_local_archive => set_archive_local_sub
220 procedure :: put_values_to_archive => put_values_to_local_netcdf_sub
221 procedure :: transform_native_to_base => transform_grid_to_grid_sub
222 procedure :: nullify_pointers => nullify_pointers_sub
223
226 generic :: enforce_limits => data_gridenforcelimits_real, &
228
231 generic :: handle_missing_values => data_gridhandlemissingdata_real, &
233
234 procedure :: calc_project_boundaries => calc_project_boundaries_sub
235 procedure :: test_for_need_to_pad_values => test_for_need_to_pad_values_fn
236
237 end type data_catalog_entry_t
238
239 integer (c_int), parameter, public :: missing_values_zero_out = 0
240 integer (c_int), parameter, public :: missing_values_replace_with_mean = 1
241
242 type (general_grid_t), public, pointer :: pgrd => null()
243
245 module procedure :: apply_scale_and_offset_float
246 module procedure :: apply_scale_and_offset_int
247 end interface apply_scale_and_offset
248
249contains
250
251 subroutine nullify_pointers_sub(this)
252
253 class(data_catalog_entry_t) :: this
254
255 if (associated(this%pGrdNative)) call grid_destroy(this%pGrdNative)
256 if (associated(this%pGrdBase)) call grid_destroy(this%pGrdBase)
257 if (associated(pgrd)) call grid_destroy(pgrd)
258
259 nullify(this%pGrdNative)
260 nullify(this%pGrdBase)
261 nullify( pgrd )
262
263 call netcdf_nullify_data_struct( ncfile=this%NCFILE )
264 call netcdf_nullify_data_struct( ncfile=this%NCFILE_ARCHIVE )
265
266 end subroutine nullify_pointers_sub
267
268!--------------------------------------------------------------------------------------------------
269
270 subroutine set_keyword_sub(this, sKeyword)
271
272 class(data_catalog_entry_t) :: this
273 character (len=*), intent(in) :: sKeyword
274
275 this%sKeyword = skeyword
276
277 end subroutine set_keyword_sub
278
279!--------------------------------------------------------------------------------------------------
280
281 subroutine get_value_int_sub(this, iCol, iRow, iValue)
282
283 class(data_catalog_entry_t), intent(in) :: this
284 integer (c_int), intent(in) :: iCol
285 integer (c_int), intent(in) :: iRow
286 integer (c_int), intent(out) :: iValue
287
288 if ( .not. associated(this%pGrdBase) ) &
289 call die("Internal programming error--attempt to use null pointer", __file__, __line__)
290
291 if (icol <= ubound(this%pGrdBase%iData,1) .and. irow <= ubound(this%pGrdBase%iData,2) ) then
292 ivalue = this%pGrdBase%iData(icol, irow)
293 else
294 call die ("Row/column indices out of bounds: ~row: "//ascharacter(irow)//"~ col:"//ascharacter(icol), &
295 __file__, __line__ )
296 endif
297
298 end subroutine get_value_int_sub
299
300!--------------------------------------------------------------------------------------------------
301
302 subroutine get_value_float_sub(this, iCol, iRow, fValue)
303
304 class(data_catalog_entry_t), intent(in) :: this
305 integer (c_int), intent(in) :: iCol
306 integer (c_int), intent(in) :: iRow
307 real (c_float), intent(out) :: fValue
308
309 if ( .not. associated(this%pGrdBase) ) &
310 call die("Internal programming error--attempt to use null pointer", __file__, __line__)
311
312 if (icol <= ubound(this%pGrdBase%iData,1) .and. irow <= ubound(this%pGrdBase%iData,2) ) then
313 fvalue = this%pGrdBase%rData(icol, irow)
314 else
315 call die ("Row/column indices out of bounds: ~row: "//ascharacter(irow)//"~ col:"//ascharacter(icol), &
316 __file__, __line__ )
317 endif
318
319 end subroutine get_value_float_sub
320
321!--------------------------------------------------------------------------------------------------
322
324 sDescription, &
325 iConstant )
326
327 class(data_catalog_entry_t) :: this
328 character (len=*), intent(in) :: sDescription
329 integer (c_int), intent(in) :: iConstant
330
331 this%iConstantValue = iconstant
332 this%sDescription = trim(sdescription)
333 this%iSourceDataForm = constant_grid
334 this%iSourceDataType = datatype_int
335 this%iTargetDataType = datatype_int
336 this%iSourceFileType = filetype_none
337
338 call this%nullify_pointers()
339
340 this%pGrdBase => grid_create(inx=bnds%iNumCols, iny=bnds%iNumRows, &
341 rx0=bnds%fX_ll, ry0=bnds%fY_ll, rgridcellsize=bnds%fGridCellSize, idatatype=datatype_int)
342
343 this%pGrdBase%sPROJ4_string = trim( bnds%sPROJ4_string )
344 this%pGrdBase%sFilename = "None: constant value entered from control file."
345
347
348!--------------------------------------------------------------------------------------------------
349
351 sDescription, &
352 rConstant )
353
354 class(data_catalog_entry_t) :: this
355 character (len=*), intent(in) :: sDescription
356 real (c_float), intent(in) :: rConstant
357
358 this%rConstantValue = rconstant
359 this%sDescription = trim(sdescription)
360 this%iSourceDataForm = constant_grid
361 this%iSourceDataType = datatype_real
362 this%iTargetDataType = datatype_real
363 this%iSourceFileType = filetype_none
364
365 call this%nullify_pointers()
366
367 this%pGrdBase => grid_create(inx=bnds%iNumCols, iny=bnds%iNumRows, &
368 rx0=bnds%fX_ll, ry0=bnds%fY_ll, rgridcellsize=bnds%fGridCellSize, idatatype=datatype_real)
369
370 this%pGrdBase%sPROJ4_string = trim( bnds%sPROJ4_string )
371 this%pGrdBase%sFilename = "None: constant value entered from control file."
372
374
375!--------------------------------------------------------------------------------------------------
376
377 subroutine initialize_table_sub( this, sDescription, sDateColumnName, sValueColumnName, sType)
378
379 class(data_catalog_entry_t) :: this
380 character (len=*), intent(in) :: sDescription
381 character (len=*), intent(in) :: sDateColumnName
382 character (len=*), intent(in) :: sValueColumnName
383 character (len=*), intent(in) :: sType
384
385 this%sDescription = trim(sdescription)
386 this%iSourceDataForm = table_lookup
387 this%iSourceDataType = datatype_real
388 this%iTargetDataType = datatype_real
389 this%iSourceFileType = filetype_ascii_table
390 this%sDateColumnName = sdatecolumnname
391 this%sValueColumnName = svaluecolumnname
392
393 call this%nullify_pointers()
394
395 select case (aslowercase(stype))
396
397 case ("float", "real" )
398
399 this%iSourceDataType = datatype_real
400 this%iTargetDataType = datatype_real
401 this%pGrdBase => grid_create(inx=bnds%iNumCols, iny=bnds%iNumRows, &
402 rx0=bnds%fX_ll, ry0=bnds%fY_ll, rgridcellsize=bnds%fGridCellSize, idatatype=datatype_real)
403
404 case ("int", "integer")
405
406 this%iSourceDataType = datatype_int
407 this%iTargetDataType = datatype_int
408 this%pGrdBase => grid_create(inx=bnds%iNumCols, iny=bnds%iNumRows, &
409 rx0=bnds%fX_ll, ry0=bnds%fY_ll, rgridcellsize=bnds%fGridCellSize, idatatype=datatype_int)
410
411 end select
412
413 this%pGrdBase%sPROJ4_string = bnds%sPROJ4_string
414 this%pGrdBase%sFilename = "None: daily value found in table of values."
415
416 end subroutine initialize_table_sub
417
418!--------------------------------------------------------------------------------------------------
419
421 sDescription, &
422 sFileType, &
423 iDataType, &
424 sFilename, &
425 sPROJ4_string )
426
427 class(data_catalog_entry_t) :: this
428 character (len=*), intent(in) :: sDescription
429 character (len=*), intent(in) :: sFileType
430 character (len=*), intent(in) :: sFilename
431 integer (c_int) :: iDataType
432 character (len=*), intent(in), optional :: sPROJ4_string
433
434 if (present(sproj4_string) ) then
435 this%sSourcePROJ4_string = trim(sproj4_string)
436 else
437 this%sSourcePROJ4_string = bnds%sPROJ4_string
438 endif
439
440 this%sSourceFilename = fully_qualified_filename( sfilename )
441
442 !> if either a '%' or '#' character is present in the filename
443 !! treat it as a template, not as a normal filename.
444 if ( scan(string=sfilename, set="%#") > 0 ) then
445
446 this%iSourceDataForm = dynamic_grid
447 this%lGridIsPersistent = true
448 this%sFilenameTemplate = trim( this%sSourceFilename )
449
450 else
451
452 this%iSourceDataForm = static_grid
453 this%lGridIsPersistent = false
454 this%sFilenameTemplate = ""
455
456 endif
457
458 this%sSourceFileType = sfiletype
459 this%iSourceFileType = this%get_filetype()
460
461 this%iSourceDataType = idatatype
462 this%iTargetDataType = idatatype
463
464 this%sDescription = trim(sdescription)
465
466 call assert(this%iSourceFileType == filetype_arc_ascii .or. &
467 this%iSourceFileType == filetype_surfer, "Only Arc ASCII or " &
468 //"Surfer grids are supported as static grid inputs (for now).", &
469 __file__, __line__)
470
471 call assert(this%iSourceDataType == datatype_int .or. &
472 this%iSourceDataType == datatype_real, "Only integer or " &
473 //"real data types are supported as static grid inputs.", &
474 __file__, __line__)
475
476 call this%nullify_pointers()
477
478 this%pGrdBase => grid_create(inx=bnds%iNumCols, iny=bnds%iNumRows, &
479 rx0=bnds%fX_ll, ry0=bnds%fY_ll, rgridcellsize=bnds%fGridCellSize, idatatype=idatatype)
480
481 this%pGrdBase%sPROJ4_string = bnds%sPROJ4_string
482 this%pGrdBase%sFilename = this%sSourceFilename
483
485
486!--------------------------------------------------------------------------------------------------
487
489 sDescription, &
490 iDataType, &
491 sFilename, &
492 sPROJ4_string )
493
494 class(data_catalog_entry_t) :: this
495 character (len=*), intent(in) :: sDescription
496 integer (c_int), intent(in) :: iDataType
497 character (len=*), intent(in) :: sFilename
498 character (len=*), intent(in), optional :: sPROJ4_string
499
500 ! [ LOCALS ]
501
502
503 if (present(sproj4_string) ) then
504 this%sSourcePROJ4_string = trim(sproj4_string)
505 else
506 this%sSourcePROJ4_string = bnds%sPROJ4_string
507 endif
508
509 this%sSourceFilename = fully_qualified_filename( sfilename )
510
511 !> if either a '%' or '#' character is present in the filename,
512 !! treat it as a template, not as a normal filename.
513 !! if there is a template, the implication is that there is
514 !! a series of files that will be read in successively, thus "dynamic" NetCDF
515! if ( scan(string=sFilename, set="%#") > 0 ) then
516
517
518!> @TODO Implement a way to read variables in via "static" NetCDF grid.
519!! In other words, a NetCDF grid having no "time" dimension or variable.
520
521
522 this%iSourceDataForm = dynamic_netcdf_grid
523 this%lGridIsPersistent = true
524 this%sFilenameTemplate = trim(sfilename)
525
526! else
527
528! !> intent of "static" NetCDF file is to house large non-changing
529! !! input grids (i.e. landuse, soils )
530! this%iSourceDataForm = STATIC_NETCDF_GRID
531! this%lGridIsPersistent = FALSE
532! this%sFilenameTemplate = ""
533
534! endif
535
536 call this%nullify_pointers()
537
538 this%pGrdBase => grid_create(inx=bnds%iNumCols, iny=bnds%iNumRows, &
539 rx0=bnds%fX_ll, ry0=bnds%fY_ll, rgridcellsize=bnds%fGridCellSize, idatatype=idatatype)
540
541 this%pGrdBase%sPROJ4_string = bnds%sPROJ4_string
542 this%pGrdBase%sFilename = this%sSourceFilename
543
544 this%sSourceFileType = "NETCDF"
545 this%iSourceFileType = this%get_filetype()
546
547 this%iTargetDataType = idatatype
548 this%iNC_FILE_STATUS = netcdf_file_closed
549
551
552!--------------------------------------------------------------------------------------------------
553
554 subroutine getvalues_sub( this, dt )
555
556 class(data_catalog_entry_t) :: this
557 type (DATETIME_T), optional :: dt
558
559 if(this%iSourceDataForm == dynamic_grid ) then
560
561 call getvalues_gridded_sub( this, dt )
562
563 elseif ( this%iSourceDataForm == dynamic_netcdf_grid ) then
564
565 call getvalues_dynamic_netcdf_sub( this, dt )
566
567 elseif ( this%iSourceDataForm == static_netcdf_grid ) then
568
569 call getvalues_static_netcdf_sub( this )
570
571 elseif( this%iSourceDataForm == table_lookup ) then
572
573 call getvalues_from_lookup_table( this, dt )
574
575 elseif(this%iSourceDataForm == static_grid ) then
576
577 call getvalues_gridded_sub( this )
578
579 elseif(this%iSourceDataForm == constant_grid ) then
580
581 call getvalues_constant_sub( this )
582
583 else
584
585 call assert(false, "Unsupported data source specified", &
586 __file__, __line__)
587
588 endif
589
590 ! if grid data hasn't changed this timestep, we do not want to *reapply* the
591 ! scale and offset values
592 if ( this%lGridHasChanged ) then
593
594 !> Now apply the user scale and offset amounts
595 if (this%iTargetDataType == datatype_real) then
596
597 call apply_scale_and_offset(fresult=this%pGrdBase%rData, fvalue=this%pGrdBase%rData, &
598 dusersuboffset=this%dUserSubOffset, &
599 duserscalefactor=this%dUserScaleFactor, &
600 duseraddoffset=this%dUserAddOffset )
601
602 elseif ( this%iTargetDataType == datatype_int ) then
603
604 call apply_scale_and_offset(iresult=this%pGrdBase%iData, ivalue=this%pGrdBase%iData, &
605 dusersuboffset=this%dUserSubOffset, &
606 duserscalefactor=this%dUserScaleFactor, &
607 duseraddoffset=this%dUserAddOffset )
608
609 else
610
611 call die("Unsupported data type specified", __file__, __line__)
612
613 endif
614
615 endif
616
617 end subroutine getvalues_sub
618
619!--------------------------------------------------------------------------------------------------
620
621elemental subroutine apply_scale_and_offset_float(fResult, fValue, dUserScaleFactor, dUserSubOffset, dUserAddOffset )
622
623 real (c_float), intent(inout) :: fresult
624 real (c_float), intent(in) :: fvalue
625 real (c_double), intent(in) :: duserscalefactor
626 real (c_double), intent(in) :: dusersuboffset
627 real (c_double), intent(in) :: duseraddoffset
628
629 fresult = real(( (fvalue - dusersuboffset) * duserscalefactor ) + duseraddoffset, c_float)
630
631end subroutine apply_scale_and_offset_float
632
633!--------------------------------------------------------------------------------------------------
634
635elemental subroutine apply_scale_and_offset_int(iResult, iValue, dUserScaleFactor, dUserSubOffset, dUserAddOffset )
636
637 integer (c_int), intent(inout) :: iresult
638 integer (c_int), intent(in) :: ivalue
639 real (c_double), intent(in) :: duserscalefactor
640 real (c_double), intent(in) :: dusersuboffset
641 real (c_double), intent(in) :: duseraddoffset
642
643 iresult = int(( ( real( ivalue, c_float) - dusersuboffset ) * duserscalefactor ) + duseraddoffset, c_int)
644
645end subroutine apply_scale_and_offset_int
646
647!--------------------------------------------------------------------------------------------------
648
649subroutine getvalues_constant_sub( this )
650
651 class(data_catalog_entry_t) :: this
652
653 if ( .not. associated(this%pGrdBase) ) &
654 call die("Internal programming error--attempt to use null pointer", __file__, __line__)
655
656 select case (this%iSourceDataType)
657
658 case ( datatype_real )
659
660 this%lGridHasChanged = true
661
662 this%pGrdBase%rData = this%rConstantValue
663
664 case ( datatype_int)
665
666 this%lGridHasChanged = true
667
668 this%pGrdBase%iData = this%iConstantValue
669
670 case default
671
672 call dump_data_structure_sub(this)
673
674 call assert(false, "INTERNAL PROGRAMMING ERROR - Unhandled data type: " &
675 //"name="//dquote(this%sDescription) &
676 //"; value="//trim(ascharacter(this%iSourceDataType)), &
677 __file__, __line__)
678
679 end select
680
681 end subroutine getvalues_constant_sub
682
683!--------------------------------------------------------------------------------------------------
684
685 subroutine getvalues_from_lookup_table( this, dt )
686
687 class(data_catalog_entry_t) :: this
688 type (DATETIME_T), intent(in) :: dt
689
690 integer (c_int) :: indx
691 integer (c_int) :: n
692 integer (c_int) :: status_code
693
694 ! [ LOCALS ]
695 type (FSTRING_LIST_T) :: slDateValues
696
697 if ( .not. associated(this%pGrdBase) ) &
698 call die("Internal programming error--attempt to use null pointer", __file__, __line__)
699
700 if ( .not. this%lTableValuesHaveBeenRetrieved ) then
701
702 this%lTableValuesHaveBeenRetrieved = true
703 this%table_indx = 1
704
705 select case (this%iSourceDataType)
706
707 case ( datatype_real )
708
709 call params%get_parameters(skey="date", slvalues=sldatevalues)
710
711 n = sldatevalues%count
712
713 !@TODO: more tests needed to ensure user can't feed in more than one column with same name,
714 !! dates out of order, dates missing, mismatched numbers of dates versus values, etc.
715 allocate(this%table_dates(n), stat=status_code)
716 allocate(this%table_values_real(n), stat=status_code)
717
718 call this%table_dates(1)%setDateFormat(this%sDefaultDateFormat)
719
720 call params%get_parameters(skey=this%sValueColumnName, fvalues=this%table_values_real)
721
722 if ( ( size( this%table_values_real,1) /= n) &
723 .or. ( size( this%table_values_real,1) == 1 ) ) &
724 call die("Did not find values associated with a required table entry ("//squote(this%sValueColumnName)//"). ", &
725 __file__, __line__)
726
727 do indx=1, n
728 call this%table_dates(indx)%parseDate(sldatevalues%get(indx),__file__, __line__)
729 enddo
730
731 case ( datatype_int)
732
733
734 case default
735
736 call dump_data_structure_sub(this)
737
738 call assert(false, "INTERNAL PROGRAMMING ERROR - Unhandled data type: " &
739 //"name="//dquote(this%sDescription) &
740 //"; value="//trim(ascharacter(this%iSourceDataType)), &
741 __file__, __line__)
742
743 end select
744
745 endif
746
747 do
748
749 if ( this%table_indx < lbound(this%table_dates,1)) exit
750 if ( this%table_indx > ubound(this%table_dates,1)) exit
751
752 if (this%table_dates(this%table_indx) < dt) then
753 this%table_indx = this%table_indx + 1
754 elseif (this%table_dates(this%table_indx) > dt) then
755 this%table_indx = this%table_indx - 1
756 else
757 exit
758 endif
759
760 enddo
761
762 if ( this%table_dates(this%table_indx) == dt ) then
763
764 select case (this%iSourceDataType)
765
766 case ( datatype_real )
767
768 this%lGridHasChanged = true
769
770 this%pGrdBase%rData = this%table_values_real(this%table_indx)
771
772 case ( datatype_int)
773
774 case default
775
776 call dump_data_structure_sub(this)
777
778 call assert(false, "INTERNAL PROGRAMMING ERROR - Unhandled data type: " &
779 //"name="//dquote(this%sDescription) &
780 //"; value="//trim(ascharacter(this%iSourceDataType)), &
781 __file__, __line__)
782
783 end select
784
785 else
786
787 call die("Missing or out-of-order value supplied for "//squote(this%sValueColumnName), &
788 __file__, __line__)
789
790 endif
791
792 end subroutine getvalues_from_lookup_table
793
794!--------------------------------------------------------------------------------------------------
795
796 subroutine dump_data_structure_sub(this)
797
798 class(data_catalog_entry_t) :: this
799
800 call logs%write("---------------------------------------------------")
801 call logs%write("DATA STRUCTURE DETAILS:")
802 call logs%write("---------------------------------------------------")
803
804 call logs%write(" catalog key word: "//dquote( this%sKeyWord ) )
805
806 call logs%write(" source data form: "//trim(ascharacter(this%iSourceDataForm)) )
807 call logs%write(" source data type: "//trim(ascharacter(this%iSourceDataType)) )
808 call logs%write(" source file type: "//trim(ascharacter(this%iSourceFileType)) )
809 call logs%write(" description: "//trim(this%sDescription) )
810 call logs%write(" source PROJ4 string: "//trim(this%sSourcePROJ4_string) )
811 call logs%write(" source file type: "//trim(this%sSourceFileType) )
812 call logs%write(" filename template: "//trim(this%sFilenameTemplate) )
813 call logs%write(" source filename: "//trim(this%sSourceFilename) )
814
815 if (associated(this%pGrdNative)) call grid_dumpgridextent(this%pGrdNative)
816 if (associated(this%pGrdBase)) call grid_dumpgridextent(this%pGrdBase)
817
818 end subroutine dump_data_structure_sub
819
820!--------------------------------------------------------------------------------------------------
821
822 subroutine getvalues_gridded_sub( this, dt )
823
824 class(data_catalog_entry_t) :: this
825 type (DATETIME_T), optional :: dt
826 logical :: lExist
827 logical :: lOpened
828
829 this%lGridHasChanged = false
830
831 do
832
833 call assert(this%iSourceFileType == filetype_arc_ascii .or. &
834 this%iSourceFileType == filetype_surfer, "INTERNAL PROGRAMMING ERROR -" &
835 //" improper file type in use for a call to this subroutine", &
836 __file__, __line__)
837
838 if(this%iSourceDataForm == dynamic_grid ) then
839
840 if(.not. present(dt) ) &
841 call assert(false, "INTERNAL PROGRAMMING ERROR - datetime object" &
842 //" must be supplied when calling this subroutine in a " &
843 //"dynamic mode.", __file__, __line__)
844
845
846 call this%make_filename(dt)
847
848 endif
849
850 ! if the source filename hasn't changed, we do not need to be here
851 if ( this%sOldFilename .strequal. this%sSourceFilename ) exit
852
853 this%sOldFilename = this%sSourceFilename
854
855 inquire(file=this%sSourceFilename, exist=lexist, opened=lopened)
856
857 ! if the file does not exist, EXIT
858 if (.not. lexist ) then
859 if ( this%lAllowMissingFiles ) then
860 exit
861 else
862 call assert( false, &
863 "Could not find input data file~filename:"//dquote(this%sSourceFilename) &
864 //"~data description: "//trim(this%sDescription))
865 endif
866 endif
867
868 call logs%write("Opening file "//dquote(this%sSourceFilename) &
869 //" for "//trim(this%sDescription)//" data.", iloglevel=log_all, lecho=true )
870
871 if ( this%lGridIsPersistent .and. associated(this%pGrdNative) ) then
872
873 call grid_readexisting ( sfilename=this%sSourceFilename, &
874 sfiletype=this%sSourceFileType, &
875 pgrd=this%pGrdNative )
876 else
877
878 ! create a grid in native coordinates of the source dataset.
879 this%pGrdNative => grid_read( sfilename=this%sSourceFilename, &
880 sfiletype=this%sSourceFileType, &
881 idatatype=this%iSourceDataType )
882
883 ! ensure that PROJ4 string is associated with the native grid
884 this%pGrdNative%sPROJ4_string = this%sSourcePROJ4_string
885 endif
886
887 this%lGridHasChanged = true
888
889 select case (this%iTargetDataType)
890
891 case ( grid_datatype_real )
892
893 call this%handle_missing_values(this%pGrdNative%rData)
894 call this%enforce_limits(this%pGrdNative%rData)
895
896 case ( grid_datatype_int )
897
898! call this%handle_missing_values(this%pGrdNative%iData)
899! call this%enforce_limits(this%pGrdNative%iData)
900
901 case default
902
903 call assert(false, "INTERNAL PROGRAMMING ERROR - Unhandled data type: value=" &
904 //trim(ascharacter(this%iSourceDataType)), &
905 __file__, __line__)
906
907 end select
908
909 call this%transform_native_to_base()
910
911 if ( .not. this%lGridIsPersistent ) call grid_destroy( this%pGrdNative )
912
913 exit
914
915 enddo
916
917 end subroutine getvalues_gridded_sub
918
919!--------------------------------------------------------------------------------------------------
920
921subroutine transform_grid_to_grid_sub(this, rX, rY)
922
923 class(data_catalog_entry_t) :: this
924 real (c_double), optional :: rX(:)
925 real (c_double), optional :: rY(:)
926
927 if (.not. associated(this%pGrdNative) ) &
928 call die("INTERNAL PROGRAMMING ERROR--Null pointer detected.", __file__, __line__)
929
930 if ( .not. associated(this%pGrdBase) ) &
931 this%pGrdBase => grid_create( inx=bnds%iNumCols, iny=bnds%iNumRows, rx0=bnds%fX_ll, ry0=bnds%fY_ll, &
932 rgridcellsize=bnds%fGridCellSize, idatatype=this%iTargetDataType )
933
934 ! only invoke the transform procedure if the PROJ4 strings are different
935 if (.not. ( this%pGrdNative%sPROJ4_string .strequal. this%pGrdBase%sPROJ4_string ) ) then
936
937 call logs%write("Transforming gridded data in file: "//dquote(this%sSourceFilename), ilinesbefore=1 )
938 call logs%write("FROM: "//squote(this%sSourcePROJ4_string), itab=2 )
939 call logs%write("TO: "//squote(this%pGrdBase%sPROJ4_string), itab=2 )
940
941 if (present(rx) .and. present(ry)) then
942 call grid_transform(pgrd=this%pGrdNative, &
943 sfromproj4=this%sSourcePROJ4_string, &
944 stoproj4=bnds%sPROJ4_string, &
945 rx=rx, &
946 ry=ry )
947 else
948 call grid_transform(pgrd=this%pGrdNative, &
949 sfromproj4=this%sSourcePROJ4_string, &
950 stoproj4=bnds%sPROJ4_string )
951
952 endif
953 !! following this call, the pGrdNative%rX and pGrdNative%rY values will be given in the
954 !! base SWB project projection
955
956 endif
957
958 if ( this%lRequireCompleteSpatialCoverage ) &
959 call assert( grid_completelycover( this%pGrdBase, this%pGrdNative ), &
960 "Transformed grid read from file "//dquote(this%sSourceFilename) &
961 //" does not completely cover your model domain.")
962
963 select case (this%iTargetDataType)
964
965 case ( grid_datatype_real )
966
967 call grid_gridtogrid_sgl(pgrdfrom=this%pGrdNative,&
968 pgrdto=this%pGrdBase )
969
970 case ( grid_datatype_int )
971
972 call grid_gridtogrid_int(pgrdfrom=this%pGrdNative, &
973 pgrdto=this%pGrdBase, &
974 lusemajorityfilter=this%lUseMajorityFilter )
975
976 case default
977
978 call assert(false, "INTERNAL PROGRAMMING ERROR - Unhandled data type: value=" &
979 //trim(ascharacter(this%iSourceDataType)), &
980 __file__, __line__)
981
982 end select
983
984
985end subroutine transform_grid_to_grid_sub
986
987!--------------------------------------------------------------------------------------------------
988
989subroutine set_constant_value_int( this, iValue )
990
991 class(data_catalog_entry_t) :: this
992 integer (c_int) :: iValue
993
994 this%iConstantValue = ivalue
995
996end subroutine set_constant_value_int
997
998!--------------------------------------------------------------------------------------------------
999
1000subroutine set_constant_value_real( this, rValue )
1001
1002 class(data_catalog_entry_t) :: this
1003 real (c_float) :: rValue
1004
1005 this%rConstantValue = rvalue
1006
1007end subroutine set_constant_value_real
1008
1009!--------------------------------------------------------------------------------------------------
1010
1011 subroutine set_filecount( this, iValue, iYear)
1012
1013 class(data_catalog_entry_t) :: this
1014 integer (c_int) :: iValue
1015 integer (c_int), optional :: iYear
1016
1017 this%iFileCount = ivalue
1018
1019 if (present(iyear) ) this%iFileCountYear = iyear
1020
1021 end subroutine set_filecount
1022
1023!--------------------------------------------------------------------------------------------------
1024
1025 subroutine increment_filecount( this )
1026
1027 class(data_catalog_entry_t) :: this
1028
1029 this%iFileCount = this%iFileCount + 1
1030
1031 end subroutine increment_filecount
1032
1033!--------------------------------------------------------------------------------------------------
1034
1035 subroutine reset_filecount( this )
1036
1037 class(data_catalog_entry_t) :: this
1038
1039 this%iFileCount = 0
1040
1041 end subroutine reset_filecount
1042
1043!--------------------------------------------------------------------------------------------------
1044
1045 subroutine reset_at_yearend_filecount( this, iYear )
1046
1047 class(data_catalog_entry_t) :: this
1048 integer (c_int) :: iYear
1049
1050 if (iyear /= this%iFileCountYear ) then
1051 this%iFileCount = 0
1052 this%iFileCountYear = iyear
1053 endif
1054
1055 end subroutine reset_at_yearend_filecount
1056
1057!--------------------------------------------------------------------------------------------------
1058
1059 subroutine make_filename_from_template( this, dt )
1060
1061 class(data_catalog_entry_t) :: this
1062 type (DATETIME_T), intent(in), optional :: dt
1063
1064 ! [ LOCALS ]
1065 character (len=512) :: sNewFilename
1066 character (len=256) :: sCWD
1067 character (len=256) :: sBuf2
1068 integer (c_int) :: iPos_Y, iPos_D, iPos_M, iPos_0D, iPos_0M, iPos_B, &
1069 iPos_BF, iPos_j, iPos, iPos2, iLen, iCount
1070 integer (c_int) :: iNumZeros, iNumZerosToPrint
1071 logical (c_bool) :: lMatch
1072 character (len=16) :: sBuf
1073 character (len=12) :: sNumber
1074 character (len=1) :: sDelimiter
1075 logical (c_bool) :: lAnnual
1076
1077 ipos_y = 0; ipos_m = 0; ipos_d = 0; ipos = 0; ipos_b = 0; ipos_bf = 0; snumber = ""
1078 ipos_j = 0
1079 lannual = false
1080
1081 ! EXAMPLES of the kinds of templates that we need to be able to understand:
1082 ! tars1980\prcp.nc template => "tars%Y\prcp.nc"
1083 ! prcp_1980_00.nc template => "prcp_%Y_%m.nc"
1084
1085! iStatus = getcwd(sCWD )
1086
1087! call assert(iStatus==0, "Problem detemining what the current working" &
1088! //" directory is", __FILE__, __LINE__)
1089
1090 snewfilename = this%sFilenameTemplate
1091
1092 icount = 0
1093
1094 do
1095
1096 lmatch = false
1097
1098 ipos_y = max(index(snewfilename, "%Y"), index(snewfilename, "%y") )
1099
1100 if (ipos_y > 0) then
1101 lmatch = true
1102 ilen=len_trim(snewfilename)
1103 snewfilename = snewfilename(1:ipos_y - 1)//trim(ascharacter(dt%iYear)) &
1104 //snewfilename(ipos_y + 2:ilen)
1105
1106 lannual = true
1107
1108 endif
1109
1110 ! evaluate template string for "#" characters
1111 ipos = index(snewfilename, "#")
1112
1113 if (ipos > 0) then
1114
1115 ! example: %000#
1116 ! trying to determine how many zero values have been inserted between % and # characters
1117 ipos2 = index(snewfilename(1:ipos),"%", back=true)
1118 sbuf2 = trim(ascharacter(this%iFileCount))
1119 inumzeros = max(0, ipos - ipos2 - 1)
1120
1121 if (inumzeros > 0) then
1122 inumzerostoprint = max(0,inumzeros - len_trim(sbuf2) + 1)
1123 snumber = repeat("0", inumzerostoprint )//trim(sbuf2)
1124 else
1125 snumber = trim(sbuf2)
1126 endif
1127
1128 lmatch = true
1129 lannual = false
1130 ilen=len_trim(snewfilename)
1131 snewfilename = snewfilename(1:ipos-2-inumzeros)//trim(snumber) &
1132 //snewfilename(ipos+1:ilen)
1133 endif
1134
1135 ! evaluate template string for "%m": month number
1136
1137 ipos_m = index(snewfilename, "%m")
1138 ipos_0m = index(snewfilename, "%0m")
1139 ipos_b = index(snewfilename, "%b")
1140 ipos_bf = index(snewfilename, "%B")
1141
1142 if ( ipos_0m > 0 ) then
1143
1144 lmatch = true
1145 lannual = false
1146 write (unit=sbuf, fmt="(i2.2)") dt%iMonth
1147
1148 ilen=len_trim(snewfilename)
1149 snewfilename = snewfilename(1:ipos_0m - 1)//trim(sbuf) &
1150 //snewfilename(ipos_0m + 3:ilen)
1151
1152 elseif ( ipos_m > 0 ) then
1153
1154 lmatch = true
1155 lannual = false
1156 sbuf = ascharacter( dt%iMonth )
1157
1158 ilen=len_trim(snewfilename)
1159 snewfilename = snewfilename(1:ipos_m - 1)//trim(sbuf) &
1160 //snewfilename(ipos_m + 2:ilen)
1161
1162 elseif ( ipos_b > 0 ) then
1163
1164 lmatch = true
1165 lannual = false
1166
1167 select case ( this% iFilename_Monthname_Capitalization_Rule )
1168
1170
1171 sbuf = months( dt%iMonth )%sName
1172 call touppercase( sbuf )
1173
1175
1176 sbuf = months( dt%iMonth )%sName
1177 call tolowercase ( sbuf )
1178
1179 case default
1180
1181 sbuf = months( dt%iMonth )%sName
1182
1183 end select
1184
1185 ilen=len_trim(snewfilename)
1186 snewfilename = snewfilename(1:ipos_b - 1)//trim(sbuf) &
1187 //snewfilename(ipos_b + 2:ilen)
1188
1189 elseif ( ipos_bf > 0 ) then
1190
1191 lmatch = true
1192 lannual = false
1193
1194 select case ( this% iFilename_Monthname_Capitalization_Rule )
1195
1197
1198 sbuf = months( dt%iMonth )%sFullName
1199 call touppercase( sbuf )
1200
1202
1203 sbuf = months( dt%iMonth )%sFullName
1204 call tolowercase( sbuf )
1205
1206 case default
1207
1208 sbuf = months( dt%iMonth )%sFullName
1209
1210 end select
1211
1212 ilen=len_trim(snewfilename)
1213 snewfilename = snewfilename(1:ipos_bf - 1)//trim(sbuf) &
1214 //snewfilename( ( ipos_bf + len_trim(sbuf) - 1):ilen)
1215
1216 endif
1217
1218 ! evaluate template string for DOY number (%j)
1219 ipos_j = max(index(snewfilename, "%J"),index(snewfilename, "%j") )
1220
1221 if (ipos_j > 0) then
1222 lmatch = true
1223 lannual = false
1224 write (unit=sbuf, fmt="(i3.3)") dt%getDayOfYear()
1225 ilen=len_trim(snewfilename)
1226 snewfilename = snewfilename(1:ipos_j - 1)//trim(sbuf) &
1227 //snewfilename(ipos_j + 2:ilen)
1228
1229 endif
1230
1231 ! evaluate template string for "%d": day number
1232
1233 ipos_d = max(index(snewfilename, "%D"),index(snewfilename, "%d") )
1234 ipos_0d = max(index(snewfilename, "%0D"), index(snewfilename, "%0d") )
1235
1236 if (ipos_0d > 0) then
1237 lmatch = true
1238 lannual = false
1239 write (unit=sbuf, fmt="(i2.2)") dt%iDay
1240 ilen=len_trim(snewfilename)
1241 snewfilename = snewfilename(1:ipos_0d - 1)//trim(sbuf) &
1242 //snewfilename(ipos_0d + 3:ilen)
1243
1244 elseif ( ipos_d > 0 ) then
1245
1246 lmatch = true
1247 lannual = false
1248 sbuf = ascharacter( dt%iDay )
1249
1250 ilen=len_trim(snewfilename)
1251 snewfilename = snewfilename(1:ipos_d - 1)//trim(sbuf) &
1252 //snewfilename(ipos_d + 2:ilen)
1253
1254 endif
1255
1256 if (.not. lmatch) exit
1257
1258 icount = icount + 1
1259
1260 ! failsafe
1261 if (icount > 4) exit
1262
1263 enddo
1264
1265 if( index(string=scwd, substring=forwardslash) > 0 ) then
1266 sdelimiter = forwardslash
1267 else
1268 sdelimiter = backslash
1269 endif
1270
1271! this%sSourceFilename = trim(sCWD)//trim(sDelimiter)//trim(sNewFilename)
1272 this%sSourceFilename = trim(snewfilename)
1273
1274 this%lIsAnnualGrid = lannual
1275
1276 end subroutine make_filename_from_template
1277
1278!--------------------------------------------------------------------------------------------------
1279
1280 function test_for_need_to_pad_values_fn(this, dt ) result(lNeedToPadData)
1281
1282 class(data_catalog_entry_t) :: this
1283 type (datetime_t), intent(in) :: dt
1284
1285 ! [ LOCALS ]
1286 logical :: lexist
1287 integer (c_int) :: idaysleftinmonth
1288 integer (c_int) :: ipos
1289 logical (c_bool) :: lneedtopaddata
1290
1291 do
1292
1293 lneedtopaddata = false
1294
1295 ipos = scan(string=trim(this%sSourceFilename), set="http://")
1296
1297 ! if this is a URL, we do not want to test for file existence using
1298 ! the Fortran "inquire" function
1299 if (this%sSourceFilename(ipos:ipos+6) == "http://") then
1300
1301 exit
1302
1303 else
1304
1305 ! does this file actually exist?
1306 inquire( file=this%sSourceFilename, exist=lexist )
1307
1308 ! if the file exists, do not bother with padding any values, unless
1309 ! we are dealing with a file that does exist but drops the last day
1310 ! or two of data (looking at you, DayMet)
1311 if ( lexist .and. ( .not. this%lIsAnnualGrid ) ) exit
1312
1313 ! if file does not exist, and we are close to the end of the year,
1314 ! assume that we should pad values at the end of the year
1315 if (dt%iMonth == 12 ) then
1316
1317 idaysleftinmonth = 31 - dt%iDay
1318
1319 if ( dt%isLeapYear() ) then
1320
1321 if ( idaysleftinmonth <= this%iDaysToPadIfLeapYear ) then
1322
1323 lneedtopaddata = true
1324 exit
1325
1326 endif
1327
1328 else ! it's not leap year
1329
1330 if ( idaysleftinmonth <= this%iDaysToPadAtYearsEnd ) then
1331
1332 lneedtopaddata = true
1333 exit
1334
1335 endif
1336
1337 endif
1338
1339 endif
1340
1341 ! if we have reached this point, we cannot locate the proper file and
1342 ! we are not within the proper range of dates to allow for padding.
1343 call assert(lexist, "The filename created from your template refers to " &
1344 //"a nonexistent file. ~ Attempted to open filename "&
1345 //dquote(this%sSourceFilename), __file__, __line__)
1346
1347 exit
1348
1349 endif
1350
1351 enddo
1352
1354
1355!--------------------------------------------------------------------------------------------------
1356
1357 subroutine getvalues_dynamic_netcdf_sub( this, dt )
1358
1359 class(data_catalog_entry_t) :: this
1360 type (DATETIME_T), intent(in) :: dt
1361
1362 ! [ LOCALS ]
1363 logical (c_bool) :: lDateTimeFound
1364 real (c_double) :: dAddOffset
1365 real (c_double) :: dScaleFactor
1366 integer (c_int) :: iMonth
1367 integer (c_int) :: iDay
1368 integer (c_int) :: iYear
1369
1370 if ( .not. associated(this%pGrdBase) ) &
1371 call die("Internal programming error--attempt to use null pointer", __file__, __line__)
1372
1373 this%lPadValues = false
1374
1375 ! call once at start of run...
1376 if ( this%iFileCountYear < 0 ) call this%set_filecount(-1, dt%iYear)
1377
1378 do
1379
1380 if (this%iNC_FILE_STATUS == netcdf_file_open) then
1381
1382 ! check to see whether currently opened file is within date range
1383 ! if past date range, close file
1384
1385 if ( netcdf_date_within_range(ncfile=this%NCFILE, &
1386 ijulianday=int(dt%iJulianDay, c_int) ) ) then
1387 exit
1388 else
1389 call netcdf_close_file( ncfile=this%NCFILE )
1390 this%iNC_FILE_STATUS = netcdf_file_closed
1391 endif
1392
1393 endif ! end of block to process if NETCDF_FILE_OPEN
1394
1395 if ( this%iNC_FILE_STATUS == netcdf_file_closed ) then
1396
1397 ! increment or reset file counter based on current year value
1398 call this%increment_filecount()
1399
1400 ! the numerical counter used in creating filenames is reset at the end of each year
1401 call this%reset_at_yearend_filecount(dt%iYear)
1402
1403 ! based on the template information, create the filename that SWB
1404 ! is to look for
1405 call this%make_filename( dt )
1406
1407 this%lPadValues = this%test_for_need_to_pad_values(dt)
1408
1409 ! call to test_for_need_to_pad_values return value of "TRUE" if
1410 ! if attempts to open a nonexistent file within the last few days of a year.
1411 ! The assumption is that values missing at the end of a calendar year
1412 ! translates into a missing file at the year's end
1413
1414 if (.not. this%lPadValues) then
1415
1416 if (this%lPerformFullInitialization ) then
1417
1418 if( ( len_trim( this%sSourcePROJ4_string ) > 0 ) &
1419 .and. ( .not. ( this%sSourcePROJ4_string .strequal. "<NA>") ) ) then
1420
1421 ! calculate the project boundaries in the coordinate system of
1422 ! the native data file
1423 call this%calc_project_boundaries(pgrdbase=this%pGrdBase)
1424
1425 if ( this%lRequireCompleteSpatialCoverage ) then
1426
1427 call netcdf_open_and_prepare_as_input(ncfile=this%NCFILE, &
1428 sfilename=this%sSourceFilename, &
1429 lfliphorizontal=this%lFlipHorizontal, &
1430 lflipvertical=this%lFlipVertical, &
1431 lallowautomaticdataflipping=this%lAllowAutomaticDataFlipping, &
1432 rx_coord_addoffset = this%rX_Coord_AddOffset, &
1433 ry_coord_addoffset = this%rY_Coord_AddOffset, &
1434 svariableorder=this%sVariableOrder, &
1435 svarname_x=this%sVariableName_x, &
1436 svarname_y=this%sVariableName_y, &
1437 svarname_z=this%sVariableName_z, &
1438 svarname_time=this%sVariableName_time, &
1439 rcoordinatetolerance=this%rCoordinateTolerance, &
1440 tgridbounds=this%GRID_BOUNDS_NATIVE )
1441
1442 else
1443
1444 call netcdf_open_and_prepare_as_input(ncfile=this%NCFILE, &
1445 sfilename=this%sSourceFilename, &
1446 lfliphorizontal=this%lFlipHorizontal, &
1447 lflipvertical=this%lFlipVertical, &
1448 lallowautomaticdataflipping=this%lAllowAutomaticDataFlipping, &
1449 rx_coord_addoffset = this%rX_Coord_AddOffset, &
1450 ry_coord_addoffset = this%rY_Coord_AddOffset, &
1451 svariableorder=this%sVariableOrder, &
1452 svarname_x=this%sVariableName_x, &
1453 svarname_y=this%sVariableName_y, &
1454 svarname_z=this%sVariableName_z, &
1455 svarname_time=this%sVariableName_time, &
1456 rcoordinatetolerance=this%rCoordinateTolerance )
1457
1458 endif
1459
1460 else ! PROJ4 string is blank
1461
1462 ! assume source NetCDF file is in same projection and
1463 ! of same dimensions as base grid
1464 call netcdf_open_and_prepare_as_input(ncfile=this%NCFILE, &
1465 sfilename=this%sSourceFilename, &
1466 lfliphorizontal=this%lFlipHorizontal, &
1467 lflipvertical=this%lFlipVertical, &
1468 lallowautomaticdataflipping=this%lAllowAutomaticDataFlipping, &
1469 svariableorder=this%sVariableOrder, &
1470 svarname_x=this%sVariableName_x, &
1471 svarname_y=this%sVariableName_y, &
1472 svarname_z=this%sVariableName_z, &
1473 rcoordinatetolerance=this%rCoordinateTolerance, &
1474 svarname_time=this%sVariableName_time )
1475
1476 this%NCFILE%iNX = this%pGrdBase%iNX
1477 this%NCFILE%iNY = this%pGrdBase%iNY
1478 this%NCFILE%rX(nc_left) = this%pGrdBase%rX0
1479 this%NCFILE%rY(nc_bottom) = this%pGrdBase%rY0
1480 this%NCFILE%rX(nc_right) = this%pGrdBase%rX1
1481 this%NCFILE%rY(nc_top) = this%pGrdBase%rY1
1482
1483 endif
1484
1485 this%rX_coordinate_subset = this%NCFILE%rX_Coords(this%NCFILE%iColBounds(nc_left):this%NCFILE%iColBounds(nc_right))
1486 this%rY_coordinate_subset = this%NCFILE%rY_Coords(this%NCFILE%iRowBounds(nc_top):this%NCFILE%iRowBounds(nc_bottom))
1487
1488 ! print *, 'getvalues_dynamic_netcdf'
1489 ! print *, trim(__FILE__), ': ', __LINE__
1490 ! print *, 'bounds (x): ', this%NCFILE%iColBounds(NC_LEFT), this%NCFILE%iColBounds(NC_RIGHT)
1491 ! print *, 'bounds (y): ', this%NCFILE%iRowBounds(NC_TOP), this%NCFILE%iRowBounds(NC_BOTTOM)
1492 ! print *, 'y-coords: ', this%rY_coordinate_subset
1493
1494 this%iNC_FILE_STATUS = netcdf_file_open
1495
1496 this%iSourceDataType = this%NCFILE%iVarType(nc_z)
1497
1498 ! Amongst other things, the call to netcdf_open_and_prepare
1499 ! finds the nearest column and row that correspond to the
1500 ! project bounds, then back-calculates the coordinate values
1501 ! of the column and row numbers in the *NATIVE* coordinate system
1502 if ( associated(this%pGrdNative) ) call grid_destroy (this%pGrdNative)
1503
1504 this%pGrdNative => grid_create( inx=this%NCFILE%iNX, &
1505 iny=this%NCFILE%iNY, &
1506 rx0=this%NCFILE%rX(nc_left), &
1507 ry0=this%NCFILE%rY(nc_bottom), &
1508 rx1=this%NCFILE%rX(nc_right), &
1509 ry1=this%NCFILE%rY(nc_top), &
1510 idatatype=this%iTargetDataType )
1511
1512 ! print *, ''
1513 ! print *, 'routine getvalues_netcdf_dynamic'
1514 ! print *, 'processing netCDF file: ', trim(this%sSourceFilename)
1515 ! print *, 'pGrdNative created with the following values:'
1516 ! print *, 'iNX=', this%NCFILE%iNX, &
1517 ! 'iNY=',this%NCFILE%iNY, &
1518 ! 'rX0=',this%NCFILE%rX(NC_LEFT), &
1519 ! 'rY0=',this%NCFILE%rY(NC_BOTTOM), &
1520 ! 'rX1=',this%NCFILE%rX(NC_RIGHT), &
1521 ! 'rY1=',this%NCFILE%rY(NC_TOP)
1522
1523 if( len_trim( this%sSourcePROJ4_string ) > 0 ) then
1524 ! ensure that PROJ4 string is associated with the native grid
1525 this%pGrdNative%sPROJ4_string = this%sSourcePROJ4_string
1526 endif
1527
1528 this%pGrdNative%sFilename = this%sSourceFilename
1529
1530 ! we do not need to perform all these steps for the next file; we are
1531 ! assuming, of course, that all of the subsequent files cover the same
1532 ! extents and are in the same projection as this first file
1533 this%lPerformFullInitialization = false
1534
1535 else
1536 ! Projection settings can be left alone; read values from new
1537 ! NetCDF file with same grid boundaries, projection, etc.
1538 call netcdf_open_file(ncfile=this%NCFILE, sfilename=this%sSourceFilename)
1539
1540 this%iNC_FILE_STATUS = netcdf_file_open
1541
1542 endif
1543
1544 if ( netcdf_date_within_range(ncfile=this%NCFILE, &
1545 ijulianday=int(dt%iJulianDay, c_int) ) ) then
1546
1547 exit
1548
1549 elseif ( scan(this%sSourceFilename, "#") /= 0 ) then
1550
1551 call netcdf_close_file( ncfile=this%NCFILE )
1552 this%iNC_FILE_STATUS = netcdf_file_closed
1553 call logs%write("Did not find the current date in the file "//dquote(this%sSourceFilename)//"~" &
1554 //"JD range: "//ascharacter(this%NCFILE%iFirstDayJD)//" to "//ascharacter(this%NCFILE%iLastDayJD) &
1555 //"~current JD: "//ascharacter(dt%iJulianDay)//"~ Will increment sequential file number and try again.", &
1556 ilinesbefore=1, ilinesafter=1 )
1557
1558 else
1559
1560 call gregorian_date(this%NCFILE%iFirstDayJD,iyear, imonth, iday )
1561 call logs%write("NetCDF start date: "//trim(ascharacter(imonth, "i2.2")) &
1562 //"/"//trim(ascharacter(iday, "i2.2"))//"/"//trim(ascharacter(iyear, "i4.4")))
1563 !//" to "//trim(asCharacter(this%NCFILE%iLastDayJD)) )
1564
1565 call gregorian_date(this%NCFILE%iLastDayJD,iyear, imonth, iday )
1566 call logs%write("NetCDF end date: "//trim(ascharacter(imonth, "i2.2")) &
1567 //"/"//trim(ascharacter(iday, "i2.2"))//"/"//trim(ascharacter(iyear, "i4.4")))
1568
1569 call gregorian_date(int(dt%iJulianDay, c_int),iyear, imonth, iday )
1570 call logs%write("Current SWB simulation date: "//trim(ascharacter(imonth, "i2.2")) &
1571 //"/"//trim(ascharacter(iday, "i2.2"))//"/"//trim(ascharacter(iyear, "i4.4")))
1572
1573 call assert (false, "Date range for currently open NetCDF file" &
1574 //" does not include the current simulation date.", &
1575 __file__, __line__)
1576
1577 endif
1578
1579 endif ! if(lPadValues)
1580
1581 exit
1582
1583 endif ! If (NC_FILE_STATUS == NETCDF_CLOSED)
1584
1585 enddo
1586
1587 if (.not. this%lPadValues) then
1588
1589 do
1590 ldatetimefound = netcdf_update_time_starting_index(ncfile=this%NCFILE, &
1591 ijulianday=int(dt%iJulianDay, c_int) )
1592
1593 if (.not. ldatetimefound) then
1594 this%lPadValues = true
1595 exit
1596 endif
1597
1598 call netcdf_get_variable_slice(ncfile=this%NCFILE, rvalues=this%pGrdNative%rData)
1599 this%lGridHasChanged = true
1600
1601 ! this initialization must take place here so that initialization may
1602 ! occur *after* the netCDF file has been opened. previously initialization
1603 ! took place tens of lines above, which resulted in an 'add_offset' of 0.0
1604 ! and a 'scale_factor' of 1.0 being applied for the first time step.
1605 daddoffset = this%NCFILE%rAddOffset(nc_z)
1606 dscalefactor = this%NCFILE%rScaleFactor(nc_z)
1607
1608 this%pGrdNative%rData = real(this%pGrdNative%rData * dscalefactor + daddoffset, c_float)
1609
1610 call this%handle_missing_values(this%pGrdNative%rData)
1611 call this%enforce_limits(this%pGrdNative%rData)
1612 exit
1613 enddo
1614
1615 endif
1616
1617 if (this%lPadValues) then
1618
1619 if (this%lPadReplaceWithZero) then
1620
1621 this%pGrdNative%rData = 0_c_float
1622 this%pGrdNative%iData = 0_c_int
1623
1624 endif
1625
1626 call logs%write( repeat("=", 60) )
1627 call logs%write( "Missing day found in NetCDF file - padding values" )
1628 call logs%write( repeat("=", 60) )
1629
1630 endif
1631
1632 if (this%lCreateLocalNetCDFArchive) &
1633 call this%put_values_to_archive(int(dt%iMonth,c_int), &
1634 int(dt%iDay,c_int), dt%iYear)
1635
1636 call this%transform_native_to_base( rx=this%rX_coordinate_subset, &
1637 ry=this%rY_coordinate_subset)
1638
1639 end subroutine getvalues_dynamic_netcdf_sub
1640
1641
1642
1643!--------------------------------------------------------------------------------------------------
1644
1646
1647 class(data_catalog_entry_t) :: this
1648
1649 ! [ LOCALS ]
1650 real (c_double) :: dAddOffset
1651 real (c_double) :: dScaleFactor
1652
1653 if ( .not. associated(this%pGrdBase) ) &
1654 call die("Internal programming error--attempt to use null pointer", __file__, __line__)
1655
1656 if ( this%iNC_FILE_STATUS == netcdf_file_closed ) then
1657
1658 if (this%lPerformFullInitialization ) then
1659
1660 if( ( len_trim( this%sSourcePROJ4_string ) > 0 ) &
1661 .and. ( .not. ( this%sSourcePROJ4_string .strequal. "<NA>") ) ) then
1662
1663 ! calculate the project boundaries in the coordinate system of
1664 ! the native data file
1665 call this%calc_project_boundaries(pgrdbase=this%pGrdBase)
1666
1667 if ( this%lRequireCompleteSpatialCoverage ) then
1668 call netcdf_open_and_prepare_as_input(ncfile=this%NCFILE, &
1669 sfilename=this%sSourceFilename, &
1670 lfliphorizontal=this%lFlipHorizontal, &
1671 lflipvertical=this%lFlipVertical, &
1672 svariableorder=this%sVariableOrder, &
1673 svarname_x=this%sVariableName_x, &
1674 svarname_y=this%sVariableName_y, &
1675 svarname_z=this%sVariableName_z, &
1676 svarname_time=this%sVariableName_time, &
1677 tgridbounds=this%GRID_BOUNDS_NATIVE )
1678 else
1679 call netcdf_open_and_prepare_as_input(ncfile=this%NCFILE, &
1680 sfilename=this%sSourceFilename, &
1681 lfliphorizontal=this%lFlipHorizontal, &
1682 lflipvertical=this%lFlipVertical, &
1683 svariableorder=this%sVariableOrder, &
1684 svarname_x=this%sVariableName_x, &
1685 svarname_y=this%sVariableName_y, &
1686 svarname_z=this%sVariableName_z, &
1687 svarname_time=this%sVariableName_time )
1688 endif
1689
1690 else ! PROJ4 string is blank
1691
1692 ! assume source NetCDF file is in same projection and
1693 ! of same dimensions as base grid
1694 call netcdf_open_and_prepare_as_input(ncfile=this%NCFILE, &
1695 sfilename=this%sSourceFilename, &
1696 lfliphorizontal=this%lFlipHorizontal, &
1697 lflipvertical=this%lFlipVertical, &
1698 svariableorder=this%sVariableOrder, &
1699 svarname_x=this%sVariableName_x, &
1700 svarname_y=this%sVariableName_y, &
1701 svarname_z=this%sVariableName_z, &
1702 svarname_time=this%sVariableName_time )
1703
1704 this%NCFILE%iNX = this%pGrdBase%iNX
1705 this%NCFILE%iNY = this%pGrdBase%iNY
1706 this%NCFILE%rX(nc_left) = this%pGrdBase%rX0
1707 this%NCFILE%rY(nc_bottom) = this%pGrdBase%rY0
1708 this%NCFILE%rX(nc_right) = this%pGrdBase%rX1
1709 this%NCFILE%rY(nc_top) = this%pGrdBase%rY1
1710
1711 endif
1712
1713 this%iNC_FILE_STATUS = netcdf_file_open
1714
1715 this%iSourceDataType = this%NCFILE%iVarType(nc_z)
1716
1717 ! Amongst other things, the call to netcdf_open_and_prepare
1718 ! finds the nearest column and row that correspond to the
1719 ! project bounds, then back-calculates the coordinate values
1720 ! of the column and row numbers in the *NATIVE* coordinate system
1721 if ( associated(this%pGrdNative) ) call grid_destroy (this%pGrdNative)
1722
1723 this%pGrdNative => grid_create( inx=this%NCFILE%iNX, &
1724 iny=this%NCFILE%iNY, &
1725 rx0=this%NCFILE%rX(nc_left), &
1726 ry0=this%NCFILE%rY(nc_bottom), &
1727 rx1=this%NCFILE%rX(nc_right), &
1728 ry1=this%NCFILE%rY(nc_top), &
1729 idatatype=this%iTargetDataType )
1730
1731 ! print *, ''
1732 ! print *, 'routine getvalues_netcdf_static'
1733 ! print *, 'processing netCDF file: ', trim(this%sSourceFilename)
1734 ! print *, 'pGrdNative created with the following values:'
1735 ! print *, 'iNX=', this%NCFILE%iNX, &
1736 ! 'iNY=',this%NCFILE%iNY, &
1737 ! 'rX0=',this%NCFILE%rX(NC_LEFT), &
1738 ! 'rY0=',this%NCFILE%rY(NC_BOTTOM), &
1739 ! 'rX1=',this%NCFILE%rX(NC_RIGHT), &
1740 ! 'rY1=',this%NCFILE%rY(NC_TOP)
1741
1742
1743 if( len_trim( this%sSourcePROJ4_string ) > 0 ) then
1744 ! ensure that PROJ4 string is associated with the native grid
1745 this%pGrdNative%sPROJ4_string = this%sSourcePROJ4_string
1746 endif
1747
1748 this%pGrdNative%sFilename = this%sSourceFilename
1749
1750 ! we do not need to perform all these steps for the next file; we are
1751 ! assuming, of course, that all of the subsequent files cover the same
1752 ! extents and are in the same projection as this first file
1753 this%lPerformFullInitialization = false
1754
1755 else
1756 ! Projection settings can be left alone; read values from new
1757 ! NetCDF file with same grid boundaries, projection, etc.
1758
1759! call netcdf_open_file(NCFILE=this%NCFILE, sFilename=this%sSourceFilename, iLU=LU_LOG)
1760 call netcdf_open_file(ncfile=this%NCFILE, sfilename=this%sSourceFilename)
1761
1762 this%iNC_FILE_STATUS = netcdf_file_open
1763
1764 endif
1765
1766
1767 endif ! If (NC_FILE_STATUS == NETCDF_CLOSED)
1768
1769 call netcdf_get_variable_slice(ncfile=this%NCFILE, rvalues=this%pGrdNative%rData)
1770
1771 daddoffset = this%NCFILE%rAddOffset(nc_z)
1772 dscalefactor = this%NCFILE%rScaleFactor(nc_z)
1773 this%pGrdNative%rData = real(this%pGrdNative%rData * dscalefactor + daddoffset, c_float)
1774
1775 call this%handle_missing_values(this%pGrdNative%rData)
1776
1777 call this%enforce_limits(this%pGrdNative%rData)
1778
1779 call this%transform_native_to_base( )
1780
1781 end subroutine getvalues_static_netcdf_sub
1782
1783!--------------------------------------------------------------------------------------------------
1784
1785 subroutine put_values_to_local_netcdf_sub(this, iMonth, iDay, iYear)
1786
1787 class(data_catalog_entry_t) :: this
1788 integer (c_int) :: iMonth
1789 integer (c_int) :: iDay
1790 integer (c_int) :: iYear
1791
1792 ! [ LOCALS ]
1793 integer (c_size_t) :: iNumRows, iNumCols, iNumRecs
1794
1795 if (this%iNC_ARCHIVE_STATUS == netcdf_file_closed) then
1796
1797 call netcdf_open_and_prepare_as_output_archive(ncfile=this%NCFILE, &
1798 ncfile_archive=this%NCFILE_ARCHIVE, &
1799 ioriginmonth=imonth, ioriginday=iday, ioriginyear=iyear, &
1800 istartyear=this%iStartYear, iendyear=this%iEndYear)
1801
1802 this%iNC_ARCHIVE_STATUS = netcdf_file_open
1803
1804 endif
1805
1806 inumrows = int(size(this%pGrdNative%rData, 2), c_size_t)
1807 inumcols = int(size(this%pGrdNative%rData, 1), c_size_t)
1808 inumrecs = this%iNCFILE_RECNUM
1809
1810 call netcdf_put_variable_array(ncfile=this%NCFILE_ARCHIVE, &
1811 ivarid=this%NCFILE_ARCHIVE%iVarID(nc_z), &
1812 istart=[inumrecs, 0_c_size_t, 0_c_size_t], &
1813 icount=[1_c_size_t, inumrows, inumcols], &
1814 istride=[1_c_size_t,1_c_size_t,1_c_size_t], &
1815 rvalues=this%pGrdNative%rData)
1816
1817 call netcdf_put_variable_vector(ncfile=this%NCFILE_ARCHIVE, &
1818 ivarid=this%NCFILE_ARCHIVE%iVarID(nc_time), &
1819 istart=[this%iNCFILE_RECNUM], &
1820 icount=[1_c_size_t], &
1821 istride=[1_c_size_t], &
1822 rvalues=[real(this%iNCFILE_RECNUM, c_float)])
1823
1824 this%iNCFILE_RECNUM = this%iNCFILE_RECNUM + 1
1825
1826 end subroutine put_values_to_local_netcdf_sub
1827
1828!--------------------------------------------------------------------------------------------------
1829
1830 function get_source_filetype_fn(this) result(iFileType)
1831
1832 class(data_catalog_entry_t) :: this
1833 integer (c_int) :: ifiletype
1834
1835 ifiletype = -1
1836
1837 if ( (this%sSourceFileType .strequal. "ARC_GRID") &
1838 .or. (this%sSourceFileType .strequal. "ARC_ASCII") ) then
1839
1840 ifiletype = filetype_arc_ascii
1841
1842 elseif ( this%sSourceFileType .strequal. "SURFER" ) then
1843
1844 ifiletype = filetype_surfer
1845
1846 elseif ( this%sSourceFileType .strequal. "NETCDF" ) then
1847
1848 ifiletype = filetype_netcdf
1849
1850 else
1851
1852 call assert(false, "Unknown input file type specified. ~"&
1853 //" filename: "//dquote(this%sSourceFilename) &
1854 //"~ file type specified as: "//dquote(this%sSourceFileType), &
1855 __file__, __line__)
1856
1857 endif
1858
1859 end function get_source_filetype_fn
1860
1861!--------------------------------------------------------------------------------------------------
1862
1863 subroutine set_source_proj4_string_sub(this, sPROJ4_string)
1864
1865 class(data_catalog_entry_t) :: this
1866 character (len=*), optional :: sPROJ4_string
1867
1868 this%sSourcePROJ4_string = sproj4_string
1869
1870 end subroutine set_source_proj4_string_sub
1871
1872!--------------------------------------------------------------------------------------------------
1873
1874 subroutine set_target_proj4_string_sub(this, sPROJ4_string)
1875
1876 class(data_catalog_entry_t) :: this
1877 character (len=*), optional :: sPROJ4_string
1878
1879 this%sTargetPROJ4_string = sproj4_string
1880
1881 end subroutine set_target_proj4_string_sub
1882
1883!--------------------------------------------------------------------------------------------------
1884
1886
1887 class(data_catalog_entry_t) :: this
1888
1889 this%lAllowAutomaticDataFlipping = false
1890
1892
1893!--------------------------------------------------------------------------------------------------
1894
1896
1897 class(data_catalog_entry_t) :: this
1898
1899 this%lFlipHorizontal = true
1900
1901 end subroutine set_grid_flip_horizontal_sub
1902
1903!--------------------------------------------------------------------------------------------------
1904
1906
1907 class(data_catalog_entry_t) :: this
1908
1909 this%lFlipVertical = true
1910
1911 end subroutine set_grid_flip_vertical_sub
1912
1913 !--------------------------------------------------------------------------------------------------
1914
1916
1917 class(data_catalog_entry_t) :: this
1918
1919 this%lAllowMissingFiles = true
1920
1922
1923!--------------------------------------------------------------------------------------------------
1924
1925 subroutine set_variable_order_sub(this, sVariableOrder)
1926
1927 class(data_catalog_entry_t) :: this
1928 character (len=*) :: sVariableOrder
1929
1930 this%sVariableOrder = svariableorder
1931
1932 end subroutine set_variable_order_sub
1933
1934!--------------------------------------------------------------------------------------------------
1935
1936subroutine set_scale_sub(this, dScaleFactor)
1937
1938 class(data_catalog_entry_t) :: this
1939 real (c_double) :: dScaleFactor
1940
1941 this%dUserScaleFactor = dscalefactor
1942
1943end subroutine set_scale_sub
1944
1945!--------------------------------------------------------------------------------------------------
1946
1947subroutine set_archive_local_sub(this, lValue)
1948
1949 class(data_catalog_entry_t) :: this
1950 logical (c_bool) :: lValue
1951
1952 this%lCreateLocalNetCDFArchive = lvalue
1953
1954end subroutine set_archive_local_sub
1955
1956!--------------------------------------------------------------------------------------------------
1957
1958subroutine set_x_coord_offset_sub(this, rXOffset)
1959
1960 class(data_catalog_entry_t) :: this
1961 real (c_double) :: rXOffset
1962
1963 this%rX_Coord_AddOffset = rxoffset
1964
1965end subroutine set_x_coord_offset_sub
1966
1967!----------------------------------------------------------------------
1968
1969subroutine set_y_coord_offset_sub(this, rYOffset)
1970
1971 class(data_catalog_entry_t) :: this
1972 real (c_double) :: rYOffset
1973
1974 this%rY_Coord_AddOffset = ryoffset
1975
1976end subroutine set_y_coord_offset_sub
1977
1978!----------------------------------------------------------------------
1979
1980subroutine set_coordinate_tolerance_sub(this, rCoordinateTolerance)
1981
1982 class(data_catalog_entry_t) :: this
1983 real (c_double) :: rCoordinateTolerance
1984
1985 this%rCoordinateTolerance = rcoordinatetolerance
1986
1987end subroutine set_coordinate_tolerance_sub
1988
1989!----------------------------------------------------------------------
1990
1991subroutine set_add_offset_sub(this, dAddOffset)
1992
1993 class(data_catalog_entry_t) :: this
1994 real (c_double) :: dAddOffset
1995
1996 this%dUserAddOffset = daddoffset
1997
1998end subroutine set_add_offset_sub
1999
2000!----------------------------------------------------------------------
2001
2002subroutine set_sub_offset_sub(this, dSubOffset)
2003
2004 class(data_catalog_entry_t) :: this
2005 real (c_double) :: dSubOffset
2006
2007 this%dUserSubOffset = dsuboffset
2008
2009end subroutine set_sub_offset_sub
2010
2011!--------------------------------------------------------------------------------------------------
2012
2013subroutine set_majority_filter_flag_sub(this, lUseMajorityFilter)
2014
2015 class(data_catalog_entry_t) :: this
2016 logical (c_bool) :: lUseMajorityFilter
2017
2018 this%lUseMajorityFilter = lusemajorityfilter
2019
2020end subroutine set_majority_filter_flag_sub
2021
2022!--------------------------------------------------------------------------------------------------
2023
2024
2025!--------------------------------------------------------------------------------------------------
2026
2027subroutine set_complete_spatial_coverage_flag_sub(this, lRequireCompleteSpatialCoverage )
2028
2029 class(data_catalog_entry_t) :: this
2030 logical (c_bool), intent(in) :: lRequireCompleteSpatialCoverage
2031
2032 this%lRequireCompleteSpatialCoverage = lrequirecompletespatialcoverage
2033
2035
2036!--------------------------------------------------------------------------------------------------
2037
2039
2040 class(data_catalog_entry_t) :: this
2041 integer (c_int) :: iMinVal
2042
2043 this%iMinAllowedValue = iminval
2044
2046
2047!--------------------------------------------------------------------------------------------------
2048
2050
2051 class(data_catalog_entry_t) :: this
2052 integer (c_int) :: iMaxVal
2053
2054 this%iMaxAllowedValue = imaxval
2055
2057
2058!--------------------------------------------------------------------------------------------------
2059
2061
2062 class(data_catalog_entry_t) :: this
2063 real (c_float) :: rMinVal
2064
2065 this%rMinAllowedValue = rminval
2066
2068
2069!--------------------------------------------------------------------------------------------------
2070
2072
2073 class(data_catalog_entry_t) :: this
2074 real (c_float) :: rMaxVal
2075
2076 this%rMaxAllowedValue = rmaxval
2077
2079
2080!--------------------------------------------------------------------------------------------------
2081
2082 subroutine calc_project_boundaries_sub(this, pGrdBase)
2083
2084 class(data_catalog_entry_t) :: this
2085 type ( GENERAL_GRID_T ), pointer :: pGrdBase
2086
2087 ! [ LOCALS ]
2088 integer (c_int) :: iRetVal
2089 real (c_double), dimension(4) :: rX, rY
2090
2091 ! ensure that there is sufficient coverage on all sides of grid
2092 rx(1) = pgrdbase%rX0 ! - pGrdBase%rGridCellSize * rMultiplier ! Xll
2093 ry(1) = pgrdbase%rY0 ! - pGrdBase%rGridCellSize * rMultiplier ! Yll
2094 rx(2) = pgrdbase%rX1 ! + pGrdBase%rGridCellSize * rMultiplier ! Xlr
2095 ry(2) = pgrdbase%rY0 ! - pGrdBase%rGridCellSize * rMultiplier ! Ylr
2096 rx(3) = pgrdbase%rX0 ! - pGrdBase%rGridCellSize * rMultiplier ! Xul
2097 ry(3) = pgrdbase%rY1 ! + pGrdBase%rGridCellSize * rMultiplier ! Yul
2098 rx(4) = pgrdbase%rX1 ! + pGrdBase%rGridCellSize * rMultiplier ! Xur
2099 ry(4) = pgrdbase%rY1 ! + pGrdBase%rGridCellSize * rMultiplier ! Yur
2100
2101 ! don't invoke PROJ4 unless projections are at least superficially different
2102 if ( .not. trim( pgrdbase%sPROJ4_string) == trim(this%sSourcePROJ4_string)) then
2103
2104 ! now transform the project coordinates to native coordinates so we can
2105 ! use the native coordinate boundaries to "cookie-cut" only the data
2106 ! pertinent to our project area.
2107 iretval = pj_init_and_transform(trim(pgrdbase%sPROJ4_string)//c_null_char, &
2108 trim(this%sSourcePROJ4_string)//c_null_char, &
2109 __file__//c_null_char, &
2110 __line__, &
2111 4_c_long, &
2112 rx, ry )
2113
2114 call grid_checkforproj4error(iretval=iretval, &
2115 sfromproj4=trim(pgrdbase%sPROJ4_string), &
2116 stoproj4=trim(this%sSourcePROJ4_string))
2117
2118 endif
2119
2120 ! because PROJ4 works in RADIANS if data are unprojected (i.e. GEOGRAPHIC),
2121 ! we need to convert back to degrees on the assumption that the coordinates
2122 ! referenced in the file will also be i degrees
2123 ! if( index(string=trim(this%sSourcePROJ4_string), substring="latlon") > 0 &
2124 ! .or. index(string=trim(this%sSourcePROJ4_string), substring="lonlat") > 0 ) then
2125
2126 if ( ( this%sSourcePROJ4_string .containssimilar. "latlon" ) &
2127 .or. ( this%sSourcePROJ4_string .containssimilar. "latlong" ) &
2128 .or. ( this%sSourcePROJ4_string .containssimilar. "lonlat" ) &
2129 .or. ( this%sSourcePROJ4_string .containssimilar. "longlat" ) ) then
2130
2131 rx = rad_to_deg(rx)
2132 ry = rad_to_deg(ry)
2133
2134 endif
2135
2136 ! GRID_BOUNDS_NATIVE will contain the SWB model bounds as defined in the
2137 ! projection of the particular dataset
2138 this%GRID_BOUNDS_NATIVE%rXll = rx(1); this%GRID_BOUNDS_NATIVE%rXlr = rx(2)
2139 this%GRID_BOUNDS_NATIVE%rYll = ry(1); this%GRID_BOUNDS_NATIVE%rYlr = ry(2)
2140 this%GRID_BOUNDS_NATIVE%rXul = rx(3); this%GRID_BOUNDS_NATIVE%rXur = rx(4)
2141 this%GRID_BOUNDS_NATIVE%rYul = ry(3); this%GRID_BOUNDS_NATIVE%rYur = ry(4)
2142
2143#ifdef DEBUG_PRINT
2144 print *, " "
2145 print *, " routine 'calc_project_boundaries'"
2146 print *, trim(__file__), ": ", __line__
2147 print *, "-- BASE GRID BOUNDS projected to DATA NATIVE COORDS"
2148 print *, "FROM: ", dquote(pgrdbase%sPROJ4_string)
2149 print *, "TO: ", dquote(this%sSourcePROJ4_string)
2150 print *, "file: ", dquote(this%sSourceFileName)
2151 print *, " X Y"
2152 print *, "LL: ", this%GRID_BOUNDS_NATIVE%rXll, this%GRID_BOUNDS_NATIVE%rYll
2153 print *, "LR: ", this%GRID_BOUNDS_NATIVE%rXlr, this%GRID_BOUNDS_NATIVE%rYlr
2154 print *, "UL: ", this%GRID_BOUNDS_NATIVE%rXul, this%GRID_BOUNDS_NATIVE%rYul
2155 print *, "UR: ", this%GRID_BOUNDS_NATIVE%rXur, this%GRID_BOUNDS_NATIVE%rYur
2156#endif
2157
2158 end subroutine calc_project_boundaries_sub
2159
2160!--------------------------------------------------------------------------------------------------
2161
2162 subroutine data_gridenforcelimits_int(this, iValues)
2163
2164 class(data_catalog_entry_t) :: this
2165 integer (c_int), dimension(:,:) :: iValues
2166
2167 ! [ LOCALS ]
2168 integer (c_int) :: iMin, iMax
2169
2170 imin = this%iMinAllowedValue
2171 imax = this%iMaxAllowedValue
2172
2173 where ( ivalues < imin ) ivalues = imin
2174 where ( ivalues > imax ) ivalues = imax
2175
2176 end subroutine data_gridenforcelimits_int
2177
2178!--------------------------------------------------------------------------------------------------
2179
2180 subroutine data_gridenforcelimits_real(this, rValues)
2181
2182 class(data_catalog_entry_t) :: this
2183 real (c_float), dimension(:,:) :: rValues
2184
2185 ! [ LOCALS ]
2186 real (c_float) :: rMin, rMax
2187
2188 rmin = real(this%rMinAllowedValue, c_float)
2189 rmax = real(this%rMaxAllowedValue, c_float)
2190
2191 where ( rvalues < rmin ) rvalues = rmin
2192 where ( rvalues > rmax ) rvalues = rmax
2193
2194 end subroutine data_gridenforcelimits_real
2195
2196!--------------------------------------------------------------------------------------------------
2197
2198 subroutine data_gridhandlemissingdata_real(this, rValues)
2199
2200 class(data_catalog_entry_t) :: this
2201 real (c_float), dimension(:,:), intent(inout) :: rValues
2202
2203 ! [ LOCALS ]
2204 real (c_float) :: rMissing, rMean
2205
2206 rmissing = real(this%rMissingValuesCode, c_float)
2207
2208 ! changing the default operation to "do nothing"
2209 ! user must actively choose what to do with missing values
2210 ! by specifying a valid operator
2211 if ( trim(this%sMissingValuesOperator) .ne. "&&" ) then
2212
2213 select case (this%iMissingValuesAction)
2214
2216
2217 select case (trim(this%sMissingValuesOperator))
2218
2219 case ("<=")
2220
2221 where (rvalues <= rmissing) rvalues = rzero
2222
2223 case ("<")
2224
2225 where (rvalues < rmissing) rvalues = rzero
2226
2227 case (">=")
2228
2229 where (rvalues >= rmissing) rvalues = rzero
2230
2231 case (">")
2232
2233 where (rvalues > rmissing) rvalues = rzero
2234
2235 case default
2236
2237 call assert(false, "Unknown missing values code was supplied " &
2238 //"for processing data "//squote(this%sDescription)//": " &
2239 //dquote(this%sMissingValuesOperator) )
2240
2241 end select
2242
2244
2245 select case (this%sMissingValuesOperator)
2246
2247 case ("<=")
2248
2249 rmean = sum(rvalues, rvalues > rmissing ) / count(rvalues > rmissing )
2250
2251 where (rvalues <= rmissing) rvalues = rmean
2252
2253 case ("<")
2254
2255 rmean = sum(rvalues, rvalues >= rmissing ) / count(rvalues >= rmissing )
2256
2257 where (rvalues < rmissing) rvalues = rmean
2258
2259 case (">=")
2260
2261 rmean = sum(rvalues, rvalues < rmissing ) / count(rvalues < rmissing )
2262
2263 where (rvalues >= rmissing) rvalues = rmean
2264
2265 case (">")
2266
2267 rmean = sum(rvalues, rvalues <= rmissing ) / count(rvalues <= rmissing )
2268
2269 where (rvalues > rmissing) rvalues = rmean
2270
2271 case default
2272
2273 call assert(false, "Unknown missing values code was supplied " &
2274 //"for processing data "//squote(this%sDescription)//": " &
2275 //dquote(this%sMissingValuesOperator) )
2276
2277 end select
2278
2279 case default
2280
2281 call assert(false, "INTERNAL PROGRAMMING ERROR - unhandled iMissingValuesAction", &
2282 __file__, __line__)
2283
2284 end select
2285
2286 endif
2287 end subroutine data_gridhandlemissingdata_real
2288
2289!--------------------------------------------------------------------------------------------------
2290
2291 subroutine data_gridhandlemissingdata_int(this, iValues)
2292
2293 class(data_catalog_entry_t) :: this
2294 integer (c_int), dimension(:,:), intent(inout) :: iValues
2295
2296 ! [ LOCALS ]
2297 integer (c_int) :: iMissing, iMean
2298
2299 imissing = this%iMissingValuesCode
2300 if ( trim(this%sMissingValuesOperator) .ne. "&&" ) then
2301
2302 select case (this%iMissingValuesAction)
2303
2305
2306 select case (trim(this%sMissingValuesOperator))
2307
2308 case ("<=")
2309
2310 where (ivalues <= imissing) ivalues = izero
2311
2312 case ("<")
2313
2314 where (ivalues < imissing) ivalues = izero
2315
2316 case (">=")
2317
2318 where (ivalues >= imissing) ivalues = izero
2319
2320 case (">")
2321
2322 where (ivalues > imissing) ivalues = izero
2323
2324 case default
2325
2326 call assert(false, "Unknown missing values code was supplied " &
2327 //"for processing data "//squote(this%sDescription)//": " &
2328 //dquote(this%sMissingValuesOperator) )
2329
2330 end select
2331
2333
2334 select case (this%sMissingValuesOperator)
2335
2336 case ("<=")
2337
2338 imean = sum(ivalues, ivalues > imissing ) &
2339 / count(ivalues > imissing )
2340
2341 where (ivalues <= imissing) ivalues = imean
2342
2343 case ("<")
2344
2345 imean = sum(ivalues, ivalues >= imissing ) &
2346 / count(ivalues >= imissing )
2347
2348 where (ivalues < imissing) ivalues = imean
2349
2350 case (">=")
2351
2352 imean = sum(ivalues, ivalues < imissing ) &
2353 / count(ivalues < imissing )
2354
2355 where (ivalues >= imissing) ivalues = imean
2356
2357 case (">")
2358
2359 imean = sum(ivalues, ivalues <= imissing ) &
2360 / count(ivalues <= imissing )
2361
2362 where (ivalues > imissing) ivalues = imean
2363
2364 case default
2365
2366 call assert(false, "Unknown missing values code was supplied " &
2367 //"for processing data "//squote(this%sDescription)//": " &
2368 //dquote(this%sMissingValuesOperator) )
2369
2370 end select
2371
2372 case default
2373
2374 call assert(false, "INTERNAL PROGRAMMING ERROR - unhandled iMissingValuesAction", &
2375 __file__, __line__)
2376
2377 end select
2378
2379 endif
2380 end subroutine data_gridhandlemissingdata_int
2381
2382end module data_catalog_entry
This module contains physical constants and convenience functions aimed at performing unit conversion...
logical(c_bool), parameter, public true
real(c_float), parameter, public rzero
integer(c_int), parameter datatype_real
logical(c_bool), parameter, public false
real(c_float), parameter, public rbigval
integer(c_int), parameter datatype_int
integer(c_int), parameter datatype_na
integer(c_int), parameter, public ibigval
subroutine initialize_gridded_data_object_sub(this, sdescription, sfiletype, idatatype, sfilename, sproj4_string)
subroutine set_coordinate_tolerance_sub(this, rcoordinatetolerance)
subroutine data_gridhandlemissingdata_real(this, rvalues)
subroutine set_minimum_allowable_value_real_sub(this, rminval)
integer(c_int), parameter filetype_ascii_table
integer(c_int), parameter static_grid
subroutine set_majority_filter_flag_sub(this, lusemajorityfilter)
subroutine getvalues_dynamic_netcdf_sub(this, dt)
subroutine data_gridhandlemissingdata_int(this, ivalues)
subroutine make_filename_from_template(this, dt)
integer(c_int), parameter table_lookup
type(general_grid_t), pointer, public pgrd
subroutine transform_grid_to_grid_sub(this, rx, ry)
subroutine set_variable_order_sub(this, svariableorder)
subroutine nullify_pointers_sub(this)
subroutine get_value_float_sub(this, icol, irow, fvalue)
elemental subroutine apply_scale_and_offset_int(iresult, ivalue, duserscalefactor, dusersuboffset, duseraddoffset)
subroutine set_maximum_allowable_value_real_sub(this, rmaxval)
subroutine set_do_not_allow_netcdf_grid_data_flipping_sub(this)
subroutine set_y_coord_offset_sub(this, ryoffset)
subroutine get_value_int_sub(this, icol, irow, ivalue)
subroutine set_sub_offset_sub(this, dsuboffset)
subroutine set_maximum_allowable_value_int_sub(this, imaxval)
integer(c_int), parameter, public missing_values_replace_with_mean
subroutine set_x_coord_offset_sub(this, rxoffset)
subroutine initialize_netcdf_data_object_sub(this, sdescription, idatatype, sfilename, sproj4_string)
subroutine set_filecount(this, ivalue, iyear)
subroutine calc_project_boundaries_sub(this, pgrdbase)
subroutine set_keyword_sub(this, skeyword)
subroutine set_target_proj4_string_sub(this, sproj4_string)
integer(c_int), parameter, public file_template_uppercase_monthname
subroutine data_gridenforcelimits_real(this, rvalues)
integer(c_int), parameter, public dynamic_grid
subroutine set_complete_spatial_coverage_flag_sub(this, lrequirecompletespatialcoverage)
subroutine getvalues_gridded_sub(this, dt)
subroutine getvalues_from_lookup_table(this, dt)
subroutine data_gridenforcelimits_int(this, ivalues)
integer(c_int), parameter filetype_arc_ascii
integer(c_int), parameter, public missing_values_zero_out
integer(c_int), parameter static_netcdf_grid
integer(c_int), parameter, public netcdf_file_open
subroutine set_minimum_allowable_value_int_sub(this, iminval)
integer(c_int), parameter filetype_surfer
subroutine set_grid_flip_vertical_sub(this)
subroutine increment_filecount(this)
integer(c_int), parameter no_grid
integer(c_int), parameter filetype_none
integer(c_int), parameter filetype_netcdf
subroutine set_constant_value_real(this, rvalue)
subroutine put_values_to_local_netcdf_sub(this, imonth, iday, iyear)
integer(c_int), parameter, public netcdf_file_closed
elemental subroutine apply_scale_and_offset_float(fresult, fvalue, duserscalefactor, dusersuboffset, duseraddoffset)
integer(c_int) function get_source_filetype_fn(this)
integer(c_int), parameter, public file_template_capitalized_monthname
subroutine getvalues_static_netcdf_sub(this)
subroutine initialize_table_sub(this, sdescription, sdatecolumnname, svaluecolumnname, stype)
subroutine reset_at_yearend_filecount(this, iyear)
subroutine set_constant_value_int(this, ivalue)
subroutine set_grid_flip_horizontal_sub(this)
subroutine getvalues_constant_sub(this)
logical(c_bool) function test_for_need_to_pad_values_fn(this, dt)
subroutine initialize_constant_int_data_object_sub(this, sdescription, iconstant)
subroutine set_archive_local_sub(this, lvalue)
subroutine set_add_offset_sub(this, daddoffset)
subroutine dump_data_structure_sub(this)
integer(c_int), parameter, public file_template_lowercase_monthname
subroutine set_allow_missing_files_flag_sub(this)
subroutine set_scale_sub(this, dscalefactor)
subroutine reset_filecount(this)
integer(c_int), parameter dynamic_netcdf_grid
integer(c_int), parameter constant_grid
subroutine initialize_constant_real_data_object_sub(this, sdescription, rconstant)
subroutine getvalues_sub(this, dt)
subroutine set_source_proj4_string_sub(this, sproj4_string)
This module contains the DATETIME_T class and associated time and date-related routines,...
Definition datetime.F90:9
type(month_t), dimension(12), target, public months
Month information.
Definition datetime.F90:148
subroutine, public gregorian_date(ijd, iyear, imonth, iday, iorigin)
Definition datetime.F90:571
subroutine, public die(smessage, smodule, iline, shints, scalledby, icalledbyline)
character(len=:) function, allocatable, public fully_qualified_filename(filename, pathname)
character(len=1), parameter, public forwardslash
Definition fstring.F90:170
character(len=1), parameter, public backslash
Definition fstring.F90:169
Provides support for input and output of gridded ASCII data, as well as for creation and destruction ...
Definition grid.F90:8
logical(c_bool) function, public grid_completelycover(pbasegrd, pothergrd, rtolerance)
Definition grid.F90:1328
type(general_grid_t) function, pointer, public grid_read(sfilename, sfiletype, idatatype)
Definition grid.F90:447
subroutine, public grid_dumpgridextent(pgrd)
Definition grid.F90:2264
integer(c_int), parameter, public grid_datatype_int
Definition grid.F90:25
subroutine, public grid_destroy(pgrd)
Definition grid.F90:366
subroutine, public grid_gridtogrid_int(pgrdfrom, pgrdto, lusemajorityfilter)
Definition grid.F90:2342
integer(c_int), parameter, public grid_datatype_real
Definition grid.F90:26
subroutine, public grid_readexisting(sfilename, sfiletype, pgrd)
Definition grid.F90:472
subroutine, public grid_gridtogrid_sgl(pgrdfrom, pgrdto)
Definition grid.F90:2420
subroutine, public grid_checkforproj4error(iretval, sfromproj4, stoproj4)
Definition grid.F90:1666
subroutine, public grid_transform(pgrd, sfromproj4, stoproj4, rx, ry)
Call PROJ4 to transform coordinates.
Definition grid.F90:1538
type(logfile_t), public logs
Definition logfiles.F90:62
Provide support for use of netCDF files as input for time-varying, gridded meteorlogic data,...
integer(c_int), parameter, public nc_bottom
subroutine, public netcdf_nullify_data_struct(ncfile)
integer(c_int), parameter, public nc_top
integer(c_int), parameter, public nc_right
subroutine, public netcdf_get_variable_slice(ncfile, rvalues, dpvalues, ivalues)
subroutine, public netcdf_close_file(ncfile)
integer(c_int), parameter, public nc_time
@TODO: implement a more flexible way of tracking variable IDs; presently the code can break if lat an...
subroutine, public netcdf_open_file(ncfile, sfilename, ilu)
subroutine, public netcdf_open_and_prepare_as_input(ncfile, sfilename, lfliphorizontal, lflipvertical, lallowautomaticdataflipping, rx_coord_addoffset, ry_coord_addoffset, svariableorder, svarname_x, svarname_y, svarname_z, svarname_time, rcoordinatetolerance, tgridbounds, ilu)
logical(c_bool) function, public netcdf_date_within_range(ncfile, ijulianday)
integer(c_int), parameter, public nc_left
subroutine, public netcdf_put_variable_array(ncfile, ivarid, istart, icount, istride, ivalues, i2values, rvalues, dpvalues)
subroutine, public netcdf_put_variable_vector(ncfile, ivarid, istart, icount, istride, ivalues, i2values, rvalues, dpvalues)
subroutine, public netcdf_open_and_prepare_as_output_archive(ncfile, ncfile_archive, ioriginmonth, ioriginday, ioriginyear, istartyear, iendyear)
logical(c_bool) function, public netcdf_update_time_starting_index(ncfile, ijulianday)
integer(c_int), parameter, public nc_z
type(parameters_t), public params
interface to C code that provides a simplified entry point to PROJ4 capabilities: it has been modifie...
Definition grid.F90:42