Soil Water Balance (SWB2)
Loading...
Searching...
No Matches
grid.F90
Go to the documentation of this file.
1!> @file
2!> Contains a single module, grid, which
3!> provides support for gridded ASCII data file and data structure operations
4!> @ingroup grid
5
6!> Provides support for input and output of gridded ASCII data,
7!> as well as for creation and destruction of grid data structures (defined types).
8module grid
9
10 use iso_c_binding
11 use constants_and_conversions
12 use exceptions
13 use fstring
14 use logfiles
15
16 implicit none
17
18 private
19
20 integer (c_int), public, parameter :: grid_nodata_int = -9999_c_int
21 real (c_float), public, parameter :: grid_nodata_real = -9999._c_float
22 real (c_double), public, parameter :: grid_nodata_double = -9999._c_double
23
24 ! @TODO: why redefine these when types have been defined in 'constants_and_conversions'?
25 integer (c_int), public, parameter :: grid_datatype_int = datatype_int
26 integer (c_int), public, parameter :: grid_datatype_real = datatype_real
27 integer (c_int), public, parameter :: grid_datatype_double = datatype_double
28
29 integer (c_int), public, parameter :: output_surfer = 0
30 integer (c_int), public, parameter :: output_arc = 1
31
32 integer (c_int), public, parameter :: grid_active_cell = 1
33 integer(c_int), parameter :: nc_fill_int = grid_nodata_int
34 real(c_float), parameter :: nc_fill_float = grid_nodata_real
35 real(c_double), parameter :: nc_fill_double = grid_nodata_double
36
37 !> interface to C code that provides a simplified entry point to PROJ4
38 !> capabilities: it has been modified so that all C pointers are kept within the
39 !> C code; no pointers are returned to fortran
40 public :: pj_init_and_transform
41 interface
42 function pj_init_and_transform(from_projection, to_projection, caller_name, &
43 caller_linenum, point_count, x, y) bind(c,name='pj_init_and_transform')
44 import
45 character(c_char) :: from_projection(*)
46 character(c_char) :: to_projection(*)
47 character(c_char) :: caller_name(*)
48 integer(c_int),value :: caller_linenum
49 integer(c_long),value :: point_count
50 real(c_double) :: x(*)
51 real(c_double) :: y(*)
52 integer(c_int) :: pj_init_and_transform
53 end function pj_init_and_transform
54 end interface
55
56
57 type, public :: general_grid_t
58 integer (c_int) :: inx ! Number of cells in the x-direction
59 integer (c_int) :: iny ! Number of cells in the y-direction
60 integer (c_int) :: inumgridcells ! Total number of grid cells
61 integer (c_int) :: idatatype ! Data type contained in the grid (integer, real, SWB cell)
62 character (len=:), allocatable :: sproj4_string ! proj4 string defining coordinate system of grid
63 character (len=:), allocatable :: sfilename ! original file name that the data was read from
64 real (c_double) :: rgridcellsize ! size of one side of a grid cell
65 integer (c_int) :: ilengthunits= -99999 ! length units code
66 real (c_double) :: rx0, rx1 ! World-coordinate range in X
67 real (c_double) :: ry0, ry1 ! World-coordinate range in Y
68
69 integer (c_int), dimension(:,:), allocatable :: idata ! Integer data
70 real (c_float), dimension(:,:), allocatable :: rdata ! Real data
71 real (c_float), dimension(:,:), allocatable :: fdata ! Float data
72 real (c_double), dimension(:,:), allocatable :: dpdata ! Douple-precision data
73 real (c_double), dimension(:,:), allocatable :: rx ! x coordinate associated with data
74 real (c_double), dimension(:,:), allocatable :: ry ! y coordinate associated with data
75 integer (c_int), dimension(:,:), allocatable :: imask ! Mask for processing results
76
77 integer (c_int) :: inodatavalue = nc_fill_int
78 real (c_float) :: rnodatavalue = nc_fill_float
79 real (c_double) :: dpnodatavalue = nc_fill_double
80! type (T_CELL), dimension(:,:), pointer :: Cells ! T_CELL objects
81 end type general_grid_t
82
83
84 type, public :: grid_bounds_t
85 real (c_double) :: rxll, ryll
86 real (c_double) :: rxul, ryul
87 real (c_double) :: rxlr, rylr
88 real (c_double) :: rxur, ryur
89 character (len=256) :: sproj4_string
90 end type grid_bounds_t
91
92
93 type error_message_t
94 character(len=40) :: cerrormessagetext
95 integer :: ierrornumber
96 end type error_message_t
97
98 ! error messages adapted from file "pj_strerrno.c" in PROJ.4.8.0
99 type(error_message_t), dimension(49) :: terrormessage = (/ &
100 error_message_t("no arguments in initialization list ", -1), &
101 error_message_t("no options found in 'init' file ", -2), &
102 error_message_t("no colon in init= string ", -3), &
103 error_message_t("projection not named ", -4), &
104 error_message_t("unknown projection id ", -5), &
105 error_message_t("effective eccentricity = 1. ", -6), &
106 error_message_t("unknown unit conversion id ", -7), &
107 error_message_t("invalid boolean param argument ", -8), &
108 error_message_t("unknown elliptical parameter name ", -9), &
109 error_message_t("reciprocal flattening (1/f) = 0 ", -10), &
110 error_message_t("|radius reference latitude| > 90 ", -11), &
111 error_message_t("squared eccentricity < 0 ", -12), &
112 error_message_t("major axis or radius = 0 or not given ", -13), &
113 error_message_t("latitude or longitude exceeded limits ", -14), &
114 error_message_t("invalid x or y ", -15), &
115 error_message_t("improperly formed DMS value ", -16), &
116 error_message_t("non-convergent inverse meridional dist ", -17), &
117 error_message_t("non-convergent inverse phi2 ", -18), &
118 error_message_t("acos/asin: |arg| >1.+1e-14 ", -19), &
119 error_message_t("tolerance condition error ", -20), &
120 error_message_t("conic lat_1 = -lat_2 ", -21), &
121 error_message_t("lat_1 >= 90 ", -22), &
122 error_message_t("lat_1 = 0 ", -23), &
123 error_message_t("lat_ts >= 90 ", -24), &
124 error_message_t("no distance between control points ", -25), &
125 error_message_t("projection not selected to be rotated ", -26), &
126 error_message_t("W <= 0 or M <= 0 ", -27), &
127 error_message_t("lsat not in 1-5 range ", -28), &
128 error_message_t("path not in range ", -29), &
129 error_message_t("h <= 0 ", -30), &
130 error_message_t("k <= 0 ", -31), &
131 error_message_t("lat_0 = 0 or 90 or alpha = 90 ", -32), &
132 error_message_t("lat_1=lat_2 or lat_1=0 or lat_2=90 ", -33), &
133 error_message_t("elliptical usage required ", -34), &
134 error_message_t("invalid UTM zone number ", -35), &
135 error_message_t("arg(s) out of range for Tcheby eval ", -36), &
136 error_message_t("failed to find projection to be rotated ", -37), &
137 error_message_t("failed to load datum shift file ", -38), &
138 error_message_t("both n ), & m must be spec'd and > 0 ", -39), &
139 error_message_t("n <= 0, n > 1 or not specified ", -40), &
140 error_message_t("lat_1 or lat_2 not specified ", -41), &
141 error_message_t("|lat_1| == |lat_2| ", -42), &
142 error_message_t("lat_0 is pi/2 from mean lat ", -43), &
143 error_message_t("unparseable coordinate system definition ", -44), &
144 error_message_t("geocentric transformation missing z or ellps ", -45), &
145 error_message_t("unknown prime meridian conversion id ", -46), &
146 error_message_t("illegal axis orientation combination ", -47), &
147 error_message_t("point not within available datum shift grids ", -48), &
148 error_message_t("invalid sweep axis, choose x or y ", -49) /)
149
150 public :: grid_gridtogrid_int, grid_gridtogrid_sgl, grid_create, grid_destroy, grid_read, &
151 grid_checkforproj4error, grid_completelycover, &
152 grid_readexisting, grid_dumpgridextent, &
153 grid_conform, grid_transform, grid_interpolate, grid_populatexy, &
154 grid_getgridx, grid_getgridy, &
155 grid_set_nodata_value, &
156 grid_createcomplete, grid_createsimple, grid_writegrid, &
157 grid_writearcgrid, grid_writesurfergrid, grid_set_output_directory_name, &
158 grid_getgridrownum, grid_getgridcolnum, &
159 grid_rowcolfallsinsidegrid, grid_coordinatesfallinsidegrid
160
161 interface grid_create
162 module procedure grid_createsimple
163 module procedure grid_createcomplete
164 end interface grid_create
165
166 ! global parameters for use with the majority filter
167 integer (c_int), parameter :: four_cells = 1
168 integer (c_int), parameter :: eight_cells = 2
169
170 integer (c_int), parameter, private :: column = 1
171 integer (c_int), parameter, private :: row = 2
172
173 !> @todo change these global (module) variables to local variables
174 integer (c_int) :: lu_temp, lu_grid
175
176 character (len=256) :: output_grid_directory_name = ""
177
178contains
179
180
181 subroutine grid_set_output_directory_name( sDirName )
182
183 character (len=*), intent(in) :: sdirname
184
185 output_grid_directory_name = trim(sdirname)
186
187 call logs%write("ASCII grids will be written to subdirectory " &
188 //dquote( output_grid_directory_name ), iloglevel=log_all, &
189 ilinesbefore=1, lecho=true )
190
191 end subroutine grid_set_output_directory_name
192
193
194!--------------------------------------------------------------------------
195
196!> Creates a grid of a specified type.
197!>
198!> Creates a grid pointer object and allocates memory for the data
199!> associated with the grid (REAL, INTEGER, or T_CELL).
200!
201!> @param iNX Number of grid cells in the x direction
202!> @param iNY Number of grid cells in the y direction
203!> @param rX0 X coordinate for the lower left corner of the grid
204!> @param rY0 Y coordinate for the lower left corner of the grid
205!> @param rX1 X coordinate for the upper right corner of the grid
206!> @param rY1 Y coordinate for the upper right corner of the grid
207!> @param iDataType Integer value corresponding to the type of data contained in the grid
208!>
209!> @return pGrd Pointer to a grid object
210function grid_createcomplete ( iNX, iNY, rX0, rY0, rX1, rY1, iDataType ) result ( pGrd )
211
212 ! ARGUMENTS
213 integer (c_int), intent(in) :: inx, iny ! Grid dimensions
214 real (c_double), intent(in) :: rx0, ry0 ! Lower-left corner (world coords)
215 real (c_double), intent(in) :: rx1, ry1 ! Upper-right corner (world coords)
216 integer (c_int), intent(in) :: idatatype ! Data type (DATATYPE_INT, etc.)
217 ! RETURN VALUE
218 type (general_grid_t), pointer :: pgrd
219 ! LOCALS
220 integer (c_int) :: istat
221
222 allocate ( pgrd, stat=istat )
223 call assert ( istat == 0, &
224 "Could not allocate pointer to T_GRID object", __file__,__line__ )
225
226 if (inx <= 0 .or. iny <= 0) then
227 call logs%write("Illegal grid dimensions: ")
228 call logs%write("iNX: "//ascharacter(inx) )
229 call logs%write("iNY: "//ascharacter(iny) )
230 call logs%write("rX0: "//ascharacter(rx0) )
231 call logs%write("rY0: "//ascharacter(ry0) )
232 call logs%write("rX1: "//ascharacter(rx1) )
233 call logs%write("rY1: "//ascharacter(ry1) )
234 call assert ( false, &
235 "INTERNAL PROGRAMMING ERROR? - Illegal grid dimensions specified", __file__,__line__)
236 endif
237
238 select case (idatatype)
239 case ( grid_datatype_int )
240
241 allocate ( pgrd%iData( inx, iny ), stat=istat )
242 call assert (istat == 0, &
243 "Could not allocate integer data", &
244 __file__,__line__)
245 pgrd%iData = 0
246
247 case ( grid_datatype_real )
248 allocate ( pgrd%rData( inx, iny ), stat=istat )
249 call assert (istat == 0, &
250 "Could not allocate real data", &
251 __file__,__line__)
252 pgrd%rData = rzero
253
254 case ( grid_datatype_double )
255 allocate ( pgrd%dpData( inx, iny ), stat=istat )
256 call assert (istat == 0, &
257 "Could not allocate double-precision data", &
258 __file__,__line__)
259 pgrd%dpData = 0.0_c_double
260
261 case default
262 call assert ( false, 'Internal error -- illegal grid data type' )
263 end select
264
265 pgrd%iDataType = idatatype
266 pgrd%iNX = inx
267 pgrd%iNY = iny
268 pgrd%rX0 = rx0
269 pgrd%rX1 = rx1
270 pgrd%rY0 = ry0
271 pgrd%rY1 = ry1
272 pgrd%rGridCellSize = (pgrd%rX1 - pgrd%rX0) / real(pgrd%iNX, c_float)
273 pgrd%iNumGridCells = inx * iny
274
275 allocate(pgrd%iMask(inx, iny))
276 pgrd%iMask = grid_active_cell
277
278end function grid_createcomplete
279
280function grid_createsimple ( iNX, iNY, rX0, rY0, rGridCellSize, iDataType ) result ( pGrd )
281 !! Creates a new iNX-by-iNY T_GRID of data type iDataType, over the extent
282 !! (rX0,rY0)-(rX1,rY1), and returns a pointer.
283 ! ARGUMENTS
284 integer (c_int), intent(in) :: inx, iny ! Grid dimensions
285 real (c_double), intent(in) :: rx0, ry0 ! Lower-left corner (world coords)
286 real (c_double), intent(in) :: rgridcellsize
287 integer (c_int), intent(in) :: idatatype ! Data type (DATATYPE_INT, etc.)
288
289 ! RETURN VALUE
290 type (general_grid_t), pointer :: pgrd
291 ! LOCALS
292 integer (c_int) :: istat
293
294 allocate ( pgrd, stat=istat )
295 call assert ( istat == 0, &
296 "Could not allocate pointer to T_GRID object", __file__,__line__ )
297 call assert ( inx > 0 .and. iny > 0, &
298 "Illegal grid dimensions specified: NX="//ascharacter(inx) &
299 //" NY="//ascharacter(iny), __file__,__line__)
300
301 select case (idatatype)
302 case ( grid_datatype_int )
303
304 allocate ( pgrd%iData( inx, iny ), stat=istat )
305 call assert (istat == 0, &
306 "Could not allocate integer data", &
307 __file__,__line__)
308 pgrd%iData = pgrd%iNoDataValue
309
310 case ( grid_datatype_real )
311 allocate ( pgrd%rData( inx, iny ), stat=istat )
312 call assert (istat == 0, &
313 "Could not allocate real data", &
314 __file__,__line__)
315 pgrd%rData = pgrd%rNoDataValue
316
317 case ( grid_datatype_double )
318 allocate ( pgrd%dpData( inx, iny ), stat=istat )
319 call assert (istat == 0, &
320 "Could not allocate double-precision data", &
321 __file__,__line__)
322 pgrd%dpData = pgrd%dpNoDataValue
323
324 case default
325 call assert ( false, 'Internal error -- illegal grid data type' )
326 end select
327
328 pgrd%iDataType = idatatype
329 pgrd%iNX = inx
330 pgrd%iNY = iny
331 pgrd%rX0 = rx0
332 pgrd%rX1 = rx0 + real(inx, c_double) * rgridcellsize
333 pgrd%rY0 = ry0
334 pgrd%rY1 = ry0 + real(iny, c_double) * rgridcellsize
335 pgrd%rGridCellSize = rgridcellsize
336 pgrd%iNumGridCells = inx * iny
337
338 allocate(pgrd%iMask(inx, iny))
339 pgrd%iMask = grid_active_cell
340
341end function grid_createsimple
342
343!!***
344
345!--------------------------------------------------------------------------
346!!****f* grid/grid_Destroy
347! NAME
348! grid_Destroy - Destroys a grid of a specified type.
349!
350! SYNOPSIS
351! Destroys a grid pointer object and deallocates memory for the data
352! associated with the grid (REAL, INTEGER, or T_CELL).
353!
354! INPUTS
355! pGrd - Pointer to a grid object
356!
357! OUTPUTS
358! None
359!
360! NOTES
361! Code refers to parameters that are set within types.f95.
362!
363! SOURCE
364
365subroutine grid_destroy ( pGrd )
366 !! Destroys the data in the GENERAL_GRID_T pGrd
367 ! ARGUMENTS
368 type ( general_grid_t ), pointer :: pgrd
369 ! LOCALS
370 integer (c_int) :: istat
371
372 if(associated(pgrd) )then
373
374 if ( pgrd%iDataType == grid_datatype_int ) then
375 deallocate ( pgrd%iData, stat=istat )
376 call assert ( istat == 0, "Failed to deallocate integer grid" )
377 else if ( pgrd%iDataType == grid_datatype_real ) then
378 deallocate ( pgrd%rData, stat=istat )
379 call assert ( istat == 0, "Failed to deallocate real grid" )
380 else if ( pgrd%iDataType == grid_datatype_double ) then
381 deallocate ( pgrd%dpData, stat=istat )
382 call assert ( istat == 0, "Failed to deallocate double-precision grid" )
383
384! else if ( pGrd%iDataType == GRID_DATATYPE_CELL_GRID ) then
385! deallocate ( pGrd%Cells, stat=iStat )
386! call assert ( iStat == 0, "Failed to deallocate cell grid" )
387! else if ( pGrd%iDataType == GRID_DATATYPE_ALL ) then
388! deallocate ( pGrd%iData, stat=iStat )
389! call assert ( iStat == 0, "Failed to deallocate integer grid" )
390! deallocate ( pGrd%rData, stat=iStat )
391! call assert ( iStat == 0, "Failed to deallocate real grid" )
392! deallocate ( pGrd%Cells, stat=iStat )
393! call assert ( iStat == 0, "Failed to deallocate cell grid" )
394 else
395 call assert ( false, "Internal error -- unknown grid type", &
396 __file__, __line__)
397 end if
398
399 if( allocated(pgrd%rX) ) then
400 deallocate( pgrd%rX, stat=istat)
401 call assert ( istat == 0, "Failed to deallocate X-coordinate data structure associated with grid", &
402 __file__, __line__ )
403 endif
404
405 if( allocated(pgrd%rY) ) then
406 deallocate( pgrd%rY, stat=istat)
407 call assert ( istat == 0, "Failed to deallocate Y-coordinate data structure associated with grid", &
408 __file__, __line__ )
409 endif
410
411 endif
412
413 pgrd => null()
414
415end subroutine grid_destroy
416!!***
417
418!--------------------------------------------------------------------------
419!!****f* grid/grid_Read
420! NAME
421! grid_Read - Reads a grid of a specified type.
422!
423! SYNOPSIS
424! Reads a grid of the
425! INPUTS
426! sFileName - Character string containing the name of the file to be read
427! sFileType - Character string indicating the type of file to be read
428! iDataType - Integer value corresponding to the type of data contained
429! in the input data file
430!
431! OUTPUTS
432! pGrd - Pointer to a grid object
433!
434! NOTES
435!
436! Current legal file types are
437! "ARC_GRID" ESRI ARC Grid (ASCII)
438! "SURFER" Golden Software SURFER grid (ASCII)
439!
440! Valid data types are (see module 'types'):
441! DATATYPE_INT
442! DATATYPE_REAL
443!
444! SOURCE
445
446function grid_read ( sFilename, sFileType, iDataType ) result ( pGrd )
447
448 ! ARGUMENTS
449 character (len=*), intent(in) :: sfilename ! Name of the grid input file
450 character (len=*), intent(in) :: sfiletype ! File type (see above)
451 integer (c_int), intent(in) :: idatatype ! T_GRID_INT or T_GRID_REAL
452 ! RETURN VALUE
453 type (general_grid_t), pointer :: pgrd
454
455 if ( trim(sfiletype) == "ARC_GRID" ) then
456 pgrd => grid_readarcgrid_fn( sfilename, idatatype )
457 else if ( trim(sfiletype) == "SURFER" ) then
458 pgrd => grid_readsurfergrid_fn( sfilename, idatatype )
459 else
460 call assert( false, "Illegal grid file type requested" )
461 end if
462
463 pgrd%sFilename = trim(sfilename)
464
465end function grid_read
466
467!--------------------------------------------------------------------------
468
469subroutine grid_readexisting ( sFileName, sFileType, pGrd )
470
471 ! ARGUMENTS
472 character (len=*), intent(in) :: sfilename ! Name of the grid input file
473 character (len=*), intent(in) :: sfiletype ! File type (see above)
474 type (general_grid_t), pointer :: pgrd
475
476 if ( trim(sfiletype) == "ARC_GRID" ) then
477 call grid_readarcgrid_sub( sfilename, pgrd )
478 else if ( trim(sfiletype) == "SURFER" ) then
479 call grid_readsurfergrid_sub( sfilename, pgrd )
480 else
481 call assert( false, "Illegal grid file type requested" )
482 end if
483
484 pgrd%sFilename = trim(sfilename)
485
486end subroutine grid_readexisting
487
488!!***
489!--------------------------------------------------------------------------
490!!****f* grid/grid_ReadArcGrid_fn
491! NAME
492! grid_ReadArcGrid_fn - Reads an ARC ASCII grid of a specified type.
493!
494! SYNOPSIS
495! Reads an ARC ASCII grid of specified data type
496!
497! INPUTS
498! sFileName - Character string containing the name of the file to be read
499! iDataType - Integer value corresponding to the type of data contained
500! in the input data file
501!
502! OUTPUTS
503! pGrd - Pointer to a grid object
504!
505! NOTES
506!
507! Valid data types are (see module 'types'):
508! DATATYPE_INT
509! DATATYPE_REAL
510!
511! SOURCE
512
513function grid_readarcgrid_fn ( sFileName, iDataType ) result ( pGrd )
514
515 ! ARGUMENTS
516 character (len=*), intent(in) :: sfilename ! Name of the grid input file
517 integer (c_int), intent(in) :: idatatype ! T_GRID_INT or T_GRID_REAL
518 ! RETURN VALUE
519 type (general_grid_t), pointer :: pgrd
520 ! LOCALS
521 character (len=256) :: sinputrecord ! Text record for input
522 character (len=256) :: sdirective ! Directive for input
523 character (len=256) :: sargument ! Argument for keyword directives
524 character (len=256) :: snodatavalue ! String to hold nodata value
525 character (len=8192) :: sbuf
526 integer (c_int) :: istat ! For "iostat="
527 integer (c_int) :: inx,iny ! Grid dimensions
528 integer (c_int) :: ihdrrecs ! Number of records in header
529 integer (c_int) :: icol,irow,k ! Loop indices for grid reading
530 real (c_double) :: rx0,rx1 ! Limits in X
531 real (c_double) :: ry0,ry1 ! Limits in Y
532 real (c_double) :: rcellsize ! Cell size
533 integer (c_int) :: icount,icumlcount
534 logical (c_bool) :: lxllcenter, lyllcenter ! Flags XLLCENTER / XLLCORNER
535 logical (c_bool) :: lfileexists
536 logical (c_bool) :: lisopen
537
538 ! Pre-scan for the number of header records and read the header
539 inquire(file=trim(sfilename), exist=lfileexists)
540 call assert( lfileexists, "The Arc ASCII grid file "//dquote(sfilename)// &
541 " could not be found.",__file__,__line__)
542
543 inquire(unit=lu_grid, opened=lisopen )
544 if (lisopen ) close( unit=lu_grid )
545
546 open ( newunit=lu_grid, iostat=istat, file=trim(sfilename) )
547 call assert( istat == 0, &
548 "Could not open input file " // trim(sfilename) )
549
550 ihdrrecs = 0
551 lxllcenter = false
552 lyllcenter = false
553 snodatavalue = ""
554
555 do
556 read ( unit=lu_grid, fmt="(a256)", iostat=istat ) sinputrecord
557 call assert ( istat == 0, &
558 "Could not read input record - file:"//trim(sfilename) )
559 call chomp ( sinputrecord, sdirective )
560 call touppercase ( sdirective )
561 call chomp ( sinputrecord, sargument )
562 if ( sdirective == "NCOLS" ) then
563 read ( unit=sargument, fmt=*, iostat=istat ) inx
564 call assert ( istat == 0, "Could not read NCOLS" )
565 ihdrrecs = ihdrrecs + 1
566 else if ( sdirective == "NROWS" ) then
567 read ( unit=sargument, fmt=*, iostat=istat ) iny
568 call assert ( istat == 0, "Could not read NROWS" )
569 ihdrrecs = ihdrrecs + 1
570 else if ( sdirective == "XLLCENTER" ) then
571 read ( unit=sargument, fmt=*, iostat=istat ) rx0
572 call assert ( istat == 0, "Could not read XLLCENTER" )
573 lxllcenter = true
574 ihdrrecs = ihdrrecs + 1
575 else if ( sdirective == "YLLCENTER" ) then
576 read ( unit=sargument, fmt=*, iostat=istat ) ry0
577 call assert ( istat == 0, "Could not read YLLCENTER" )
578 lxllcenter = true
579 ihdrrecs = ihdrrecs + 1
580 else if ( sdirective == "XLLCORNER" ) then
581 read ( unit=sargument, fmt=*, iostat=istat ) rx0
582 call assert ( istat == 0, "Could not read XLLCORNER" )
583 ihdrrecs = ihdrrecs + 1
584 else if ( sdirective == "YLLCORNER" ) then
585 read ( unit=sargument, fmt=*, iostat=istat ) ry0
586 call assert ( istat == 0, "Could not read YLLCORNER" )
587 ihdrrecs = ihdrrecs + 1
588 else if ( sdirective == "CELLSIZE" ) then
589 read ( unit=sargument, fmt=*, iostat=istat ) rcellsize
590 call assert ( istat == 0, "Could not read CELLSIZE" )
591 ihdrrecs = ihdrrecs + 1
592 else if ( sdirective == "NODATA_VALUE" ) then
593 snodatavalue = trim(sargument)
594 ihdrrecs = ihdrrecs + 1
595 else
596 ! Found the data -- construct the grid
597 if ( lxllcenter ) rx0 = rx0 - 0.5_c_float*rcellsize
598 if ( lyllcenter ) ry0 = ry0 - 0.5_c_float*rcellsize
599 rx1 = rx0 + real(inx, c_double) * rcellsize
600 ry1 = ry0 + real(iny, c_double) * rcellsize
601
602 pgrd => grid_create( inx, iny, rx0, ry0, rx1, ry1, idatatype )
603
604 pgrd%rGridCellSize = rcellsize
605 ! Go back to the top, skip the header...
606 rewind( unit=lu_grid, iostat=istat )
607 call assert ( istat == 0, "Failed to rewind grid file" )
608 do icol=1,ihdrrecs
609 read ( unit=lu_grid, fmt="(a256)", iostat=istat ) sinputrecord
610 call assert ( istat == 0, &
611 "Could not read input record - file: "//trim(sfilename))
612 end do
613 ! ... and read the data.
614 select case ( idatatype )
615
616 case ( datatype_int )
617
618 do irow=1,pgrd%iNY
619 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%iData(:,irow)
620 call assert ( istat == 0, &
621 "Failed to read integer grid data - file: " &
622 //trim(sfilename)//" row num: "//trim( ascharacter(irow)), &
623 __file__,__line__ )
624 end do
625 if(len_trim(snodatavalue) > 0) then
626 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%iNoDataValue
627 call assert ( istat == 0, &
628 "Failed to read NODATA value in grid data - file: " &
629 //dquote(sfilename), __file__,__line__ )
630 endif
631
632 case ( datatype_real )
633
634 do irow=1,pgrd%iNY
635 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%rData(:,irow)
636 call assert ( istat == 0, &
637 "Failed to read real grid data - file: " &
638 //trim(sfilename)//" row num: "//trim( ascharacter(irow)), &
639 __file__,__line__ )
640 end do
641 if(len_trim(snodatavalue) > 0) then
642 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%rNoDataValue
643 call assert ( istat == 0, &
644 "Failed to read NODATA value in grid data - file: " &
645 //trim(sfilename), __file__,__line__ )
646 endif
647
648 case ( datatype_double )
649
650 do irow=1,pgrd%iNY
651 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%dpData(:,irow)
652 call assert ( istat == 0, &
653 "Failed to read double-precision grid data - file: " &
654 //trim(sfilename)//" row num: "//trim( ascharacter(irow)), &
655 __file__,__line__ )
656 end do
657 if(len_trim(snodatavalue) > 0) then
658 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%dpNoDataValue
659 call assert ( istat == 0, &
660 "Failed to read NODATA value in grid data - file: " &
661 //trim(sfilename), __file__,__line__ )
662 endif
663
664 case default
665
666 call assert ( false, &
667 "Internal error -- illegal ARC GRID data type", &
668 __file__,__line__)
669
670 end select
671 exit
672 end if
673 end do
674
675 if ( idatatype == datatype_int ) call grid_checkintegergridvalues(pgrd, sfilename)
676
677 open ( unit=lu_grid, iostat=istat )
678 call assert ( istat == 0, "Failed to close grid file" )
679
680end function grid_readarcgrid_fn
681
682!--------------------------------------------------------------------------
683
684subroutine grid_readarcgrid_sub ( sFileName, pGrd )
685
686 ! ARGUMENTS
687 character (len=*), intent(in) :: sFileName ! Name of the grid input file
688 type (GENERAL_GRID_T), pointer :: pGrd
689
690 ! LOCALS
691 character (len=256) :: sInputRecord ! Text record for input
692 character (len=256) :: sDirective ! Directive for input
693 character (len=256) :: sArgument ! Argument for keyword directives
694 character (len=256) :: sNoDataValue ! String to hold nodata value
695 character (len=8192) :: sBuf
696 integer (c_int) :: iStat ! For "iostat="
697 integer (c_int) :: iNX,iNY ! Grid dimensions
698 integer (c_int) :: iHdrRecs ! Number of records in header
699 integer (c_int) :: iCol,iRow,k ! Loop indices for grid reading
700 real (c_double) :: rX0,rX1 ! Limits in X
701 real (c_double) :: rY0,rY1 ! Limits in Y
702 real (c_double) :: rCellSize ! Cell size
703 integer (c_int) :: iCount,iCumlCount
704 logical (c_bool) :: lXLLCenter, lYLLCenter ! Flags XLLCENTER / XLLCORNER
705 logical (c_bool) :: lFileExists
706 logical (c_bool) :: lIsOpen
707
708 ! Pre-scan for the number of header records and read the header
709
710 inquire(file=trim(sfilename), exist=lfileexists, opened=lisopen )
711 call assert( lfileexists, "The Arc ASCII grid file "//dquote(sfilename)// &
712 " could not be found.",__file__,__line__)
713
714 inquire(unit=lu_grid, opened=lisopen )
715 if (lisopen ) close( unit=lu_grid )
716
717 open ( newunit=lu_grid, iostat=istat, file=trim(sfilename) )
718 call assert( istat == 0, &
719 "Could not open input file " // trim(sfilename) )
720
721 ihdrrecs = 0
722 lxllcenter = false
723 lyllcenter = false
724 snodatavalue = ""
725
726 do
727 read ( unit=lu_grid, fmt="(a256)", iostat=istat ) sinputrecord
728 call assert ( istat == 0, &
729 "Could not read input record - file:"//trim(sfilename) )
730 call chomp ( sinputrecord, sdirective )
731 call touppercase ( sdirective )
732 call chomp ( sinputrecord, sargument )
733 if ( sdirective == "NCOLS" ) then
734 read ( unit=sargument, fmt=*, iostat=istat ) inx
735 call assert ( istat == 0, "Could not read NCOLS" )
736 ihdrrecs = ihdrrecs + 1
737 else if ( sdirective == "NROWS" ) then
738 read ( unit=sargument, fmt=*, iostat=istat ) iny
739 call assert ( istat == 0, "Could not read NROWS" )
740 ihdrrecs = ihdrrecs + 1
741 else if ( sdirective == "XLLCENTER" ) then
742 read ( unit=sargument, fmt=*, iostat=istat ) rx0
743 call assert ( istat == 0, "Could not read XLLCENTER" )
744 lxllcenter = true
745 ihdrrecs = ihdrrecs + 1
746 else if ( sdirective == "YLLCENTER" ) then
747 read ( unit=sargument, fmt=*, iostat=istat ) ry0
748 call assert ( istat == 0, "Could not read YLLCENTER" )
749 lxllcenter = true
750 ihdrrecs = ihdrrecs + 1
751 else if ( sdirective == "XLLCORNER" ) then
752 read ( unit=sargument, fmt=*, iostat=istat ) rx0
753 call assert ( istat == 0, "Could not read XLLCORNER" )
754 ihdrrecs = ihdrrecs + 1
755 else if ( sdirective == "YLLCORNER" ) then
756 read ( unit=sargument, fmt=*, iostat=istat ) ry0
757 call assert ( istat == 0, "Could not read YLLCORNER" )
758 ihdrrecs = ihdrrecs + 1
759 else if ( sdirective == "CELLSIZE" ) then
760 read ( unit=sargument, fmt=*, iostat=istat ) rcellsize
761 call assert ( istat == 0, "Could not read CELLSIZE" )
762 ihdrrecs = ihdrrecs + 1
763 else if ( sdirective == "NODATA_VALUE" ) then
764 snodatavalue = trim(sargument)
765 ihdrrecs = ihdrrecs + 1
766 else
767 ! Found the data -- construct the grid
768
769 ! Go back to the top, skip the header...
770 rewind( unit=lu_grid, iostat=istat )
771 call assert ( istat == 0, "Failed to rewind grid file" )
772 do icol=1,ihdrrecs
773 read ( unit=lu_grid, fmt="(a256)", iostat=istat ) sinputrecord
774 call assert ( istat == 0, &
775 "Could not read input record - file: "//trim(sfilename))
776 end do
777 ! ... and read the data.
778 select case ( pgrd%iDataType )
779
780 case ( datatype_int )
781 do irow=1,pgrd%iNY
782 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%iData(:,irow)
783 call assert ( istat == 0, &
784 "Failed to read integer grid data - file: " &
785 //trim(sfilename)//"; row num: "//trim( ascharacter(irow)) &
786 //"; line num: "//trim( ascharacter(irow + ihdrrecs)) &
787 //'; error code: '//as_character(istat), &
788 __file__,__line__ )
789 end do
790 if(len_trim(snodatavalue) > 0) then
791 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%iNoDataValue
792 call assert ( istat == 0, &
793 "Failed to read NODATA value in grid data - file: " &
794 //trim(sfilename), __file__,__line__ )
795 endif
796
797 case ( datatype_real )
798
799 do irow=1,pgrd%iNY
800 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%rData(:,irow)
801
802 call assert ( istat == 0, &
803 "Failed to read real grid data - file: " &
804 //trim(sfilename)//"; row num: "//trim( ascharacter(irow)) &
805 //"; line num: "//trim( ascharacter(irow + ihdrrecs)) &
806 //'; error code: '//as_character(istat), &
807 __file__,__line__ )
808 end do
809 if(len_trim(snodatavalue) > 0) then
810 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%rNoDataValue
811 call assert ( istat == 0, &
812 "Failed to read NODATA value in grid data - file: " &
813 //trim(sfilename), __file__,__line__ )
814 endif
815
816 case ( datatype_double )
817 do irow=1,pgrd%iNY
818 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%dpData(:,irow)
819 call assert ( istat == 0, &
820 "Failed to read double-precision grid data - file: " &
821 //trim(sfilename)//"; row num: "//trim( ascharacter(irow)) &
822 //"; line num: "//trim( ascharacter(irow + ihdrrecs)) &
823 //'; error code: '//as_character(istat), &
824 __file__,__line__ )
825 end do
826 if(len_trim(snodatavalue) > 0) then
827 read(unit=snodatavalue, fmt=*, iostat=istat) pgrd%dpNoDataValue
828 call assert ( istat == 0, &
829 "Failed to read NODATA value in grid data - file: " &
830 //trim(sfilename), __file__,__line__ )
831 endif
832
833 case default
834 call assert ( false, &
835 "Internal error -- illegal ARC GRID data type", &
836 __file__,__line__)
837 end select
838 exit
839 end if
840 end do
841
842 if ( pgrd%iDataType == datatype_int ) call grid_checkintegergridvalues(pgrd, sfilename)
843
844 open ( unit=lu_grid, iostat=istat )
845 call assert ( istat == 0, "Failed to close grid file" )
846
847end subroutine grid_readarcgrid_sub
848
849!!***
850!--------------------------------------------------------------------------
851!!****f* grid/grid_ReadSurferGrid_fn
852! NAME
853! grid_ReadSurferGrid_fn - Reads an Surfer ASCII grid of a specified type.
854!
855! SYNOPSIS
856! Reads an Surfer ASCII grid of specified data type
857!
858! INPUTS
859! sFileName - Character string containing the name of the file to be read
860! iDataType - Integer value corresponding to the type of data contained
861! in the input data file
862!
863! OUTPUTS
864! pGrd - Pointer to a grid object
865!
866! NOTES
867!
868! Valid data types are (see module 'types'):
869! DATATYPE_INT
870! DATATYPE_REAL
871!
872! SOURCE
873
874function grid_readsurfergrid_fn ( sFileName, iDataType ) result ( pGrd )
875
876 ! ARGUMENTS
877 character (len=*), intent(in) :: sfilename ! Name of the grid input file
878 integer (c_int), intent(in) :: idatatype ! T_GRID_INT or T_GRID_REAL
879 ! RETURN VALUE
880 type (general_grid_t), pointer :: pgrd
881 ! LOCALS
882 character (len=4) :: ssentinel ! Better be "DSAA" for SURFER!
883 integer (c_int) :: istat ! For "iostat="
884 integer (c_int) :: inx,iny ! Grid dimensions
885 integer (c_int) :: icol,irow ! Loop indices for grid reading
886 real (c_double) :: rx0,rx1 ! Limits in X
887 real (c_double) :: ry0,ry1 ! Limits in Y
888 real (c_float) :: rz0,rz1 ! Limits in Z (not used)
889 logical (c_bool) :: lfileexists
890 logical (c_bool) :: lisopen
891
892 inquire(file=trim(sfilename), exist=lfileexists)
893 call assert( lfileexists, "The Surfer ASCII grid file "//dquote(sfilename)// &
894 " could not be found.",__file__,__line__)
895
896 inquire(unit=lu_grid, opened=lisopen )
897 if (lisopen ) close( unit=lu_grid )
898
899 open ( newunit=lu_grid, iostat=istat, file=trim(sfilename) )
900 call assert( istat == 0, &
901 "Could not open input file " // trim(sfilename) )
902
903 read ( unit=lu_grid, fmt=*, iostat=istat ) ssentinel
904 call assert ( istat == 0, &
905 "Could not read first record of SURFER grid" )
906 call assert ( LOGICAL(trim(sSentinel) == "DSAA",c_bool), &
907 "This is not a SURFER grid" )
908 read ( unit=lu_grid, fmt=*, iostat=istat ) inx, iny
909 call assert ( istat == 0, &
910 "Error reading SURFER grid dimensions" )
911 read ( unit=lu_grid, fmt=*, iostat=istat ) rx0, rx1
912 call assert ( istat == 0, &
913 "Error reading SURFER X limits" )
914 read ( unit=lu_grid, fmt=*, iostat=istat ) ry0, ry1
915 call assert ( istat == 0, &
916 "Error reading SURFER y limits" )
917 read ( unit=lu_grid, fmt=*, iostat=istat ) rz0, rz1
918 call assert ( istat == 0, &
919 "Error reading SURFER Z limits" )
920
921 pgrd => grid_create( inx, iny, rx0, ry0, rx1, ry1, idatatype )
922 select case ( idatatype )
923 case ( datatype_int )
924 do irow=1, iny
925 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%iData(:,irow)
926 call assert ( istat == 0, "Failed to read integer grid data" )
927 end do
928 case ( datatype_real )
929 do irow=1, iny
930 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%rData(:,irow)
931 call assert ( istat == 0, "Failed to read real grid data" )
932 end do
933 case ( datatype_double )
934 do irow=1, iny
935 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%dpData(:,irow)
936 call assert ( istat == 0, "Failed to read double-precision grid data" )
937 end do
938 case default
939 call assert ( false, "Internal error -- illegal SURFER grid data type" )
940 end select
941
942 call assert(inx>0, "Must have a non-zero number of columns surfer grid file...")
943 call assert(iny>0, "Must have a non-zero number of rows in a surfer grid file...")
944
945 pgrd%rGridCellSize = (rx1-rx0)/inx
946
947 if ( idatatype == datatype_int ) call grid_checkintegergridvalues(pgrd, sfilename)
948
949 open ( unit=lu_grid, iostat=istat )
950 call assert ( istat == 0, "Failed to close grid file" )
951
952end function grid_readsurfergrid_fn
953
954!------------------------------------------------------------------------------
955
956subroutine grid_readsurfergrid_sub ( sFileName, pGrd )
957
958 ! ARGUMENTS
959 character (len=*), intent(in) :: sFileName ! Name of the grid input file
960 type (GENERAL_GRID_T), pointer :: pGrd
961
962 ! LOCALS
963 character (len=4) :: sSentinel ! Better be "DSAA" for SURFER!
964 integer (c_int) :: iStat ! For "iostat="
965 integer (c_int) :: iNX,iNY ! Grid dimensions
966 integer (c_int) :: iCol,iRow ! Loop indices for grid reading
967 real (c_double) :: rX0,rX1 ! Limits in X
968 real (c_double) :: rY0,rY1 ! Limits in Y
969 real (c_float) :: rZ0,rZ1 ! Limits in Z (not used)
970 logical (c_bool) :: lFileExists
971 logical (c_bool) :: lIsOpen
972
973 inquire(file=trim(sfilename), exist=lfileexists)
974 call assert( lfileexists, "The Surfer ASCII grid file "//dquote(sfilename)// &
975 " could not be found.",__file__,__line__)
976
977 inquire(unit=lu_grid, opened=lisopen )
978 if (lisopen ) close( unit=lu_grid )
979
980 open ( newunit=lu_grid, iostat=istat, file=trim(sfilename) )
981 call assert( istat == 0, &
982 "Could not open input file " // trim(sfilename) )
983
984 read ( unit=lu_grid, fmt=*, iostat=istat ) ssentinel
985 call assert ( istat == 0, &
986 "Could not read first record of SURFER grid" )
987 call assert ( LOGICAL(trim(sSentinel) == "DSAA",c_bool), &
988 "This is not a SURFER grid" )
989 read ( unit=lu_grid, fmt=*, iostat=istat ) inx, iny
990 call assert ( istat == 0, &
991 "Error reading SURFER grid dimensions" )
992 read ( unit=lu_grid, fmt=*, iostat=istat ) rx0, rx1
993 call assert ( istat == 0, &
994 "Error reading SURFER X limits" )
995 read ( unit=lu_grid, fmt=*, iostat=istat ) ry0, ry1
996 call assert ( istat == 0, &
997 "Error reading SURFER y limits" )
998 read ( unit=lu_grid, fmt=*, iostat=istat ) rz0, rz1
999 call assert ( istat == 0, &
1000 "Error reading SURFER Z limits" )
1001
1002 select case ( pgrd%iDataType )
1003 case ( datatype_int )
1004 do irow=1, iny
1005 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%iData(:,irow)
1006 call assert ( istat == 0, "Failed to read integer grid data" )
1007 end do
1008 case ( datatype_real )
1009 do irow=1, iny
1010 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%rData(:,irow)
1011 call assert ( istat == 0, "Failed to read real grid data" )
1012 end do
1013 case ( datatype_double )
1014 do irow=1, iny
1015 read ( unit=lu_grid, fmt=*, iostat=istat ) pgrd%dpData(:,irow)
1016 call assert ( istat == 0, "Failed to read double-precision grid data" )
1017 end do
1018 case default
1019 call assert ( false, "Internal error -- illegal SURFER grid data type" )
1020 end select
1021
1022 call assert(inx > 0,"Must have a non-zero number of grid cells in a surfer grid file...")
1023 call assert(iny > 0,"Must have a non-zero number of grid cells in a surfer grid file...")
1024
1025 if ( pgrd%iDataType == datatype_int ) call grid_checkintegergridvalues(pgrd, sfilename)
1026
1027 open ( unit=lu_grid, iostat=istat )
1028 call assert ( istat == 0, "Failed to close grid file" )
1029
1030end subroutine grid_readsurfergrid_sub
1031
1032!------------------------------------------------------------------------------
1033
1034subroutine grid_writegrid(sFilename, pGrd, iOutputFormat)
1035
1036 ! [ ARGUMENTS ]
1037 character (len=*),intent(in) :: sfilename
1038 type (general_grid_t), pointer :: pgrd
1039 integer (c_int) :: ioutputformat
1040
1041 if ( ioutputformat == output_arc ) then
1042
1043 call grid_writearcgrid(trim(output_grid_directory_name)//sfilename, pgrd)
1044
1045 elseif ( ioutputformat == output_surfer ) then
1046
1047 call grid_writesurfergrid(trim(output_grid_directory_name)//sfilename, pgrd)
1048
1049 endif
1050
1051end subroutine grid_writegrid
1052
1053!------------------------------------------------------------------------------
1054
1055subroutine grid_writearcgrid(sFilename, pGrd)
1056
1057 ! [ ARGUMENTS ]
1058 character (len=*),intent(in) :: sfilename
1059 type (general_grid_t), pointer :: pgrd
1060
1061
1062 ! [ LOCALS ]
1063 integer (c_int) :: icol,irow
1064 integer (c_int) :: inumcols, inumrows
1065 integer (c_int) :: istat
1066 character(len=256) :: sbuf
1067
1068 if ( pgrd%iDataType == datatype_int ) then
1069 inumcols = size(pgrd%iData,1)
1070 inumrows = size(pgrd%iData,2)
1071 elseif ( pgrd%iDataType == datatype_real ) then
1072 inumcols = size(pgrd%rData,1)
1073 inumrows = size(pgrd%rData,2)
1074 elseif ( pgrd%iDataType == datatype_double ) then
1075 inumcols = size(pgrd%dpData,1)
1076 inumrows = size(pgrd%dpData,2)
1077 else
1078 call assert(false, "Internal programming error - Unsupported grid type", &
1079 __file__, __line__)
1080 endif
1081
1082 ! dynamically create the Fortran output format
1083 write(sbuf,fmt="(a,a,a)") '(',trim( ascharacter(inumcols)),'(a,1x))'
1084
1085 open ( lu_temp, file=trim(output_grid_directory_name)//sfilename, iostat=istat, status="REPLACE" )
1086 call assert( istat==0, "Could not open output file "//dquote(sfilename), &
1087 __file__,__line__)
1088
1089 write ( unit=lu_temp, fmt="('NCOLS ',i10)", iostat=istat ) inumcols
1090 call assert( istat==0, "Error writing grid file header", __file__, __line__)
1091
1092 write ( unit=lu_temp, fmt="('NROWS ',i10)", iostat=istat ) inumrows
1093 call assert( istat==0, "Error writing grid file header", __file__, __line__)
1094
1095 write ( unit=lu_temp, fmt="('XLLCORNER ',f14.3)", iostat=istat ) pgrd%rX0
1096 call assert( istat==0, "Error writing X limits", __file__, __line__)
1097
1098 write ( unit=lu_temp, fmt="('YLLCORNER ',f14.3)", iostat=istat ) pgrd%rY0
1099 call assert( istat==0, "Error writing Y limits", __file__, __line__)
1100
1101 write ( unit=lu_temp, fmt="('CELLSIZE ',f14.3)", iostat=istat ) pgrd%rGridCellSize
1102 call assert( istat==0, "Error writing cell size", __file__, __line__)
1103
1104 if ( pgrd%iDataType == datatype_int ) then
1105
1106 write ( unit=lu_temp, fmt="('NODATA_VALUE ',i14)", iostat=istat ) pgrd%iNoDataValue
1107 call assert( istat==0, "Error writing NODATA value", __file__, __line__)
1108 do irow=1,inumrows
1109 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1110 (trim( ascharacter(pgrd%iData(icol,irow) ) ),icol=1,inumcols)
1111 call assert( istat==0, "Error writing Arc ASCII INTEGER grid data", &
1112 __file__, __line__)
1113 end do
1114
1115 elseif ( pgrd%iDataType == datatype_real ) then
1116
1117 write ( unit=lu_temp, fmt="('NODATA_VALUE ',g0.4)", iostat=istat ) pgrd%rNoDataValue
1118 call assert( istat==0, "Error writing NODATA value", __file__, __line__)
1119 do irow=1,inumrows
1120 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1121 (trim(ascharacter( pgrd%rData(icol,irow), fmt_string="g0.5" )),icol=1,inumcols)
1122 call assert( istat==0, "Error writing Arc ASCII REAL grid data", &
1123 __file__, __line__)
1124 end do
1125
1126 elseif ( pgrd%iDataType == datatype_double ) then
1127
1128 write ( unit=lu_temp, fmt="('NODATA_VALUE ',g0.4)", iostat=istat ) pgrd%dpNoDataValue
1129 call assert( istat==0, "Error writing NODATA value", __file__, __line__)
1130 do irow=1,inumrows
1131 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1132 (trim(ascharacter( pgrd%dpData(icol,irow), fmt_string="g0.5" )),icol=1,inumcols)
1133 call assert( istat==0, "Error writing Arc ASCII REAL grid data", &
1134 __file__, __line__)
1135 end do
1136
1137 endif
1138
1139 close (unit=lu_temp)
1140
1141end subroutine grid_writearcgrid
1142
1143!------------------------------------------------------------------------------
1144
1145subroutine grid_writesurfergrid(sFilename, pGrd)
1146
1147 ! [ ARGUMENTS ]
1148 character (len=*),intent(in) :: sfilename
1149 type (general_grid_t), pointer :: pgrd
1150
1151 ! [ LOCALS ]
1152 integer (c_int) :: icol,irow
1153 integer (c_int) :: inumcols, inumrows
1154 integer (c_int) :: istat
1155 character(len=256) :: sbuf
1156 real (c_float) :: fhalfcell
1157
1158 if ( pgrd%iDataType == datatype_int ) then
1159 inumcols = size(pgrd%iData,1)
1160 inumrows = size(pgrd%iData,2)
1161 elseif ( pgrd%iDataType == datatype_real ) then
1162 inumcols = size(pgrd%rData,1)
1163 inumrows = size(pgrd%rData,2)
1164 elseif ( pgrd%iDataType == datatype_double ) then
1165 inumcols = size(pgrd%dpData,1)
1166 inumrows = size(pgrd%dpData,2)
1167 else
1168 call assert(false, "Internal programming error - Unsupported grid type", &
1169 __file__, __line__)
1170 endif
1171
1172 fhalfcell = pgrd%rGridCellSize * 0.5_c_float
1173
1174 ! dynamically create the Fortran output format
1175 write(sbuf,fmt="(a,a,a)") '(',trim( ascharacter(inumcols)),'(a,1x))'
1176
1177 open ( lu_temp, file=trim(output_grid_directory_name)//sfilename, iostat=istat, status="REPLACE" )
1178 call assert( istat==0, "Could not open output file "//dquote(sfilename), &
1179 __file__,__line__)
1180
1181 write ( unit=lu_temp, fmt="('DSAA')", iostat=istat )
1182 call assert( istat==0, "Error writing SURFER header", &
1183 __file__, __line__)
1184
1185 write ( unit=lu_temp, fmt="(2i8)", iostat=istat ) inumcols, inumrows
1186 call assert( istat==0, "Error writing SURFER dimensions", &
1187 __file__, __line__)
1188
1189 write ( unit=lu_temp, fmt="(2f14.3)", iostat=istat ) &
1190 pgrd%rX0 + fhalfcell , pgrd%rX1 - fhalfcell
1191 call assert( istat==0, "Error writing SURFER X limits", &
1192 __file__, __line__)
1193
1194 write ( unit=lu_temp, fmt="(2f14.3)", iostat=istat ) &
1195 pgrd%rY0 + fhalfcell, pgrd%rY1 - fhalfcell
1196 call assert( istat==0, "Error writing SURFER Y limits", &
1197 __file__, __line__)
1198
1199 if ( pgrd%iDataType == datatype_int ) then
1200
1201 write ( unit=lu_temp, fmt="(2i14)", iostat=istat ) minval(pgrd%iData),maxval(pgrd%iData)
1202 call assert( istat==0, "Error writing SURFER Z limits", &
1203 __file__, __line__)
1204
1205
1206 do irow=inumrows,1,-1
1207 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1208 (trim(ascharacter(pgrd%iData(icol,irow) ) ),icol=1,inumcols)
1209 call assert( istat==0, "Error writing SURFER grid data" , &
1210 __file__, __line__)
1211 end do
1212
1213 elseif ( pgrd%iDataType == datatype_real ) then
1214
1215 write ( unit=lu_temp, fmt="(2f14.3)", iostat=istat ) minval(pgrd%rData),maxval(pgrd%rData)
1216 call assert( istat==0, "Error writing SURFER Z limits", &
1217 __file__, __line__)
1218
1219 do irow=inumrows,1,-1
1220 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1221 (trim( ascharacter(pgrd%rData(icol,irow) ) ),icol=1,inumcols)
1222 call assert( istat==0, "Error writing SURFER grid data" , &
1223 __file__, __line__)
1224 end do
1225
1226 elseif ( pgrd%iDataType == datatype_double ) then
1227
1228 write ( unit=lu_temp, fmt="(2f14.3)", iostat=istat ) minval(pgrd%dpData),maxval(pgrd%dpData)
1229 call assert( istat==0, "Error writing SURFER Z limits", &
1230 __file__, __line__)
1231
1232 do irow=inumrows,1,-1
1233 write( unit=lu_temp, fmt=trim(sbuf), iostat=istat ) &
1234 (trim( ascharacter(pgrd%dpData(icol,irow) ) ),icol=1,inumcols)
1235 call assert( istat==0, "Error writing SURFER grid data" , &
1236 __file__, __line__)
1237 end do
1238
1239 endif
1240
1241 close (unit=lu_temp)
1242
1243end subroutine grid_writesurfergrid
1244!!***
1245!--------------------------------------------------------------------------
1246!!****f* grid/grid_Conform
1247! NAME
1248! grid_WriteArcGrid - Reads an Surfer ASCII grid of a specified type.
1249!
1250! SYNOPSIS
1251! Write an ARC ASCII grid of specified data type
1252!
1253! INPUTS
1254!
1255! pGrd1 - Pointer to first grid
1256! pGrd2 - Pointer to second grid
1257! rTolerance - OPTIONAL - Allowable tolerance between the two grids
1258!
1259! OUTPUTS
1260! lConform - Logical, TRUE if the grids are compatible within the
1261! specified tolerance, FALSE if otherwise
1262!
1263! SOURCE
1264
1265function grid_conform ( pGrd1, pGrd2, rTolerance ) result ( lConform )
1266 !! Returns .true. if the T_GRID objects conform (in terms of cell sizes and extents)
1267 !! The optional argument rTolerance is the precision for checking the (floating-point)
1268 !! extent coordinates (this defaults to rDEFAULT_TOLERANCE, below). The tolerance
1269 !! is set to abs(rTolerance * ( rX1-rX1 ) )
1270 ! ARGUMENTS
1271 type (general_grid_t), pointer :: pgrd1,pgrd2
1272 real (c_float), intent(in), optional :: rtolerance
1273 ! RETURN VALUE
1274 logical (c_bool) :: lconform
1275 ! LOCALS
1276 real (c_float) :: rtol
1277 real (c_float), parameter :: rdefault_tolerance = 0.5_c_float
1278
1279 if ( present ( rtolerance ) ) then
1280 rtol = rtolerance * ( pgrd1%rX1 - pgrd1%rX0 )
1281 else
1282 rtol = rdefault_tolerance * ( pgrd1%rX1 - pgrd1%rX0 )
1283 end if
1284
1285 lconform = true
1286
1287 if ( pgrd1%iNX /= pgrd2%iNX ) then
1288 lconform = false
1289 call logs%write("Unequal number of columns between grids", iloglevel=log_all)
1290 endif
1291
1292 if ( pgrd1%iNY /= pgrd2%iNY ) then
1293 lconform = false
1294 call logs%write("Unequal number of rows between grids", iloglevel=log_all)
1295 endif
1296
1297 if( abs( pgrd1%rX0 - pgrd2%rX0 ) > rtol ) then
1298 call logs%write("Lower left-hand side X coordinates do not match:", iloglevel=log_all)
1299 call logs%write("Grid 1 value: "//ascharacter(pgrd1%rX0) &
1300 //"; grid 2 value: "//ascharacter(pgrd2%rX0), iloglevel=log_all)
1301 lconform = false
1302 endif
1303
1304 if( abs( pgrd1%rY0 - pgrd2%rY0 ) > rtol ) then
1305 call logs%write("Lower left-hand side Y coordinates do not match:", iloglevel=log_all)
1306 call logs%write("Grid 1 value: "//ascharacter(pgrd1%rY0) &
1307 //"; grid 2 value: "//ascharacter(pgrd2%rY0), iloglevel=log_all)
1308 lconform = false
1309 endif
1310
1311 if( abs( pgrd1%rX1 - pgrd2%rX1 ) > rtol ) then
1312 call logs%write("Upper right-hand side X coordinates do not match:", iloglevel=log_all)
1313 call logs%write("Grid 1 value: "//ascharacter(pgrd1%rX1) &
1314 //"; grid 2 value: "//ascharacter(pgrd2%rX1), iloglevel=log_all)
1315 lconform = false
1316 endif
1317
1318 if( abs( pgrd1%rY1 - pgrd2%rY1 ) > rtol ) then
1319 call logs%write("Upper right-hand side Y coordinates do not match:", iloglevel=log_all)
1320 call logs%write("Grid 1 value: "//ascharacter(pgrd1%rY1) &
1321 //"; grid 2 value: "//ascharacter(pgrd2%rY1), iloglevel=log_all)
1322 lconform = false
1323 endif
1324
1325end function grid_conform
1326
1327!------------------------------------------------------------------------------
1328
1329function grid_completelycover( pBaseGrd, pOtherGrd, rTolerance ) result ( lCompletelyCover )
1330 !! Returns .true. if the T_GRID objects conform (in terms of cell sizes and extents)
1331 !! The optional argument rTolerance is the precision for checking the (floating-point)
1332 !! extent coordinates (this defaults to rDEFAULT_TOLERANCE, below). The tolerance
1333 !! is set to abs(rTolerance * ( rX1-rX1 ) )
1334 ! ARGUMENTS
1335 type (general_grid_t), pointer :: pbasegrd,pothergrd
1336 real (c_float), intent(in), optional :: rtolerance
1337 ! RETURN VALUE
1338 logical (c_bool) :: lcompletelycover
1339 ! LOCALS
1340 real (c_float) :: rtol
1341 real (c_float) :: rdefault_tolerance
1342
1343 rdefault_tolerance = -pbasegrd%rGridCellSize / 100.
1344 !rDEFAULT_TOLERANCE = rZERO
1345
1346 if ( present ( rtolerance ) ) then
1347 rtol = rtolerance
1348 else
1349 rtol = rdefault_tolerance
1350 end if
1351
1352 if ( (pbasegrd%rX0 - pothergrd%rX0 >= rtol) .and. &
1353 (pbasegrd%rY0 - pothergrd%rY0 >= rtol) .and. &
1354 (pothergrd%rX1 - pbasegrd%rX1 >= rtol) .and. &
1355 (pothergrd%rY1 - pbasegrd%rY1 >= rtol) ) then
1356 lcompletelycover = true
1357 else
1358 lcompletelycover = false
1359
1360 call logs%write("Extents of the data grid file "//dquote(pothergrd%sFilename) &
1361 //" do not cover the base grid extents.")
1362
1363 call logs%write("BASE GRID EXTENTS DATA GRID EXTENTS DIFFERENCE", &
1364 itab=18, ilinesbefore=1)
1365 call logs%write("X (lower-left): "//trim(ascharacter(pbasegrd%rX0)) &
1366 //" "//trim(ascharacter(pothergrd%rX0)) &
1367 //" "//trim(ascharacter(pbasegrd%rX0 - pothergrd%rX0)) )
1368 call logs%write("Y (lower-left): "//trim(ascharacter(pbasegrd%rY0)) &
1369 //" "//trim(ascharacter(pothergrd%rY0)) &
1370 //" "//trim(ascharacter(pbasegrd%rY0 - pothergrd%rY0)) )
1371 call logs%write("X (upper-right): "//trim(ascharacter(pbasegrd%rX1)) &
1372 //" "//trim(ascharacter(pothergrd%rX1)) &
1373 //" "//trim(ascharacter(pothergrd%rX1 - pbasegrd%rX1)) )
1374 call logs%write("Y (upper-right): "//trim(ascharacter(pbasegrd%rY1)) &
1375 //" "//trim(ascharacter(pothergrd%rY1)) &
1376 //" "//trim(ascharacter(pothergrd%rY1 - pbasegrd%rY1)), ilinesafter=1 )
1377
1378 end if
1379
1380end function grid_completelycover
1381
1382!!***
1383!--------------------------------------------------------------------------
1384!!****s* grid/grid_LookupColumn
1385! NAME
1386! grid_LookupColumn - Finds the column position of the value rXval
1387! in the grid
1388!
1389! SYNOPSIS
1390! Finds the column position of the value rXval in the grid and returns:
1391! iBefore = column before the value rYval
1392! iAfter = column after the value rYval
1393! rFrac = fraction of the distance between iBefore and iAfter
1394! If iBefore or iAfter are outside the grid, they're set to -1
1395!
1396! INPUTS
1397! pGrd - Pointer to a data grid
1398! rXval - Real value to be interpolated from the gridded data.
1399! rTolerance - OPTIONAL - Allowable tolerance between the two grids
1400!
1401! OUTPUTS
1402! iBefore - column before the value rXval
1403! iAfter - column after the value rXval
1404! rFrac - fraction of the distance between iBefore and iAfter
1405!
1406! If iBefore or iAfter are outside the grid, they're set to -1
1407!
1408! NOTES
1409! The grid lookup functions were designed specifically to work with the
1410! Thornthwaite-Mather soil moisture retention table. The functions will
1411! work with any data that are similarly arranged.
1412!
1413! If X=0.65, for example, and the following table were queried, the subroutine
1414! will return: iBefore=1, iAfter=2, rFrac=0.3
1415!
1416! X=Max SM Capacity
1417!
1418! 0.50 1.00 1.50
1419! _________________l
1420! 0.0 | 0.50 1.00 1.50
1421! Y=APWL 0.1 | 0.45 0.90 1.40
1422! 0.2 | 0.40 0.80 1.30
1423!
1424! SOURCE
1425
1426subroutine grid_lookupcolumn(pGrd,rXval,iBefore,iAfter,rFrac)
1427
1428 ! [ ARGUMENTS ]
1429 type ( GENERAL_GRID_T ),pointer :: pGrd
1430 real (c_float),intent(in) :: rXval
1431 integer (c_int),intent(out) :: iBefore,iAfter
1432 real (c_float),intent(out) :: rFrac
1433 ! [ LOCALS ]
1434 real (c_float) :: rColPosition
1435
1436 rcolposition = (pgrd%iNX - 1) * (rxval - pgrd%rX0) / (pgrd%rX1 - pgrd%rX0)
1437 rfrac = rzero
1438 ibefore = floor(rcolposition) + 1
1439 iafter = ibefore + 1
1440 if ( ibefore > pgrd%iNX .or. ibefore < 1) then
1441 ibefore = -1
1442 iafter = -1
1443 else if ( iafter > pgrd%iNX ) then
1444 iafter = -1
1445 else
1446 rfrac = mod(rcolposition, 1.0_c_float )
1447 end if
1448
1449end subroutine grid_lookupcolumn
1450!!***
1451!--------------------------------------------------------------------------
1452!!****s* grid/grid_LookupRow
1453! NAME
1454! grid_LookupRow - Finds the row position of the value rYval
1455! in the grid
1456!
1457! SYNOPSIS
1458! Finds the column position of the value rYval in the grid and returns:
1459! iBefore = column before the value rYval
1460! iAfter = column after the value rYval
1461! rFrac = fraction of the distance between iBefore and iAfter
1462! If iBefore or iAfter are outside the grid, they're set to -1
1463!
1464! INPUTS
1465! pGrd - Pointer to a data grid
1466! rYval - Real value to be interpolated from the gridded data.
1467! rTolerance - OPTIONAL - Allowable tolerance between the two grids
1468!
1469! OUTPUTS
1470! iBefore - column before the value rXval
1471! iAfter - column after the value rXval
1472! rFrac - fraction of the distance between iBefore and iAfter
1473!
1474! If iBefore or iAfter are outside the grid, they're set to -1
1475!
1476! NOTES
1477! The grid lookup functions were designed specifically to work with the
1478! Thornthwaite-Mather soil moisture retention table. The functions will
1479! work with any data that are similarly arranged.
1480!
1481! If Y=0.15, for example, and the following table were queried, the subroutine
1482! will return: iBefore=2, iAfter=3, rFrac=0.5
1483!
1484! X=Max SM Capacity
1485!
1486! 0.50 1.00 1.50
1487! _________________l
1488! 0.0 | 0.50 1.00 1.50
1489! Y=APWL 0.1 | 0.45 0.90 1.40
1490! 0.2 | 0.40 0.80 1.30
1491!
1492! SOURCE
1493
1494subroutine grid_lookuprow(pGrd,rYval,iBefore,iAfter,rFrac)
1495
1496 ! [ ARGUMENTS ]
1497 type ( GENERAL_GRID_T ),pointer :: pGrd
1498 real (c_float),intent(in) :: rYval
1499 integer (c_int),intent(out) :: iBefore,iAfter
1500 real (c_float),intent(out) :: rFrac
1501 ! [ LOCALS ]
1502 real (c_float) :: rColPosition
1503
1504 rcolposition = (pgrd%iNY - 1) * (pgrd%rY1 - ryval) / (pgrd%rY1 - pgrd%rY0)
1505 rfrac = rzero
1506 ibefore = floor(rcolposition) + 1
1507 iafter = ibefore + 1
1508 if ( ibefore > pgrd%iNY ) then
1509 ibefore = -1
1510 iafter = -1
1511 else if ( ibefore < 1) then
1512 ibefore = 1
1513 iafter = 1
1514 else if ( iafter > pgrd%iNY ) then
1515 iafter = ibefore
1516 else
1517 rfrac = mod(rcolposition, 1.0_c_float )
1518 end if
1519
1520end subroutine grid_lookuprow
1521
1522!------------------------------------------------------------------------------
1523
1524!> Call PROJ4 to transform coordinates.
1525!! @details This subroutine calls a Fortran wrapper to the C library
1526!! PROJ4. A set of input coordinates is transformed to the base SWB
1527!! coordinate system.
1528!!
1529!! The idea is to first create a 2D array of the X and Y coordinates in
1530!! the projection system of the input data set. These coordinate values are
1531!! then fed to PROJ4, which modifies the values and returns the coordinate
1532!! values in the projection system of the base grid. Thereafter, mapping
1533!! source to target is a matter of finding the value of the cell closest to
1534!! a SWB grid coordinate pair.
1535!! @param[inout] pGrd
1536!! @param[in] sFromPROJ4
1537!! @param[in] sToPROJ4
1538
1539subroutine grid_transform(pGrd, sFromPROJ4, sToPROJ4, rX, rY )
1540
1541 type ( general_grid_t ),pointer :: pgrd
1542 character (len=*) :: sfromproj4, stoproj4
1543 character (kind=c_char, len=len_trim(sFromPROJ4)) :: csfromproj4
1544 character (kind=c_char, len=len_trim(sToPROJ4)) :: cstoproj4
1545 real (c_double), optional :: rx(:)
1546 real (c_double), optional :: ry(:)
1547
1548 ! [ LOCALS ]
1549 integer (c_int) :: iretval
1550 integer (c_int) :: i
1551 logical (c_bool), dimension(pGrd%iNY, pGrd%iNX) :: lmask
1552
1553 if (present(rx) .and. present(ry)) then
1554
1555 !> supply X and Y for the untransformed data
1556 call grid_populatexy(pgrd, rx, ry)
1557
1558 else
1559
1560 !> calculate X and Y for the untransformed data
1561 call grid_populatexy(pgrd)
1562
1563 endif
1564
1565 csfromproj4 = trim(sfromproj4)
1566 cstoproj4 = trim(stoproj4)
1567
1568 !> PROJ4 expects unprojected coordinates (i.e. lat lon) to be provided
1569 !> in RADIANS. Therefore, we convert to radians prior to the call...
1570
1571 if ( ( csfromproj4 .containssimilar. "latlon" ) &
1572 .or. ( csfromproj4 .containssimilar. "latlong" ) &
1573 .or. ( csfromproj4 .containssimilar. "lonlat" ) &
1574 .or. ( csfromproj4 .containssimilar. "longlat" ) ) then
1575
1576 pgrd%rX = pgrd%rX * degrees_to_radians
1577 pgrd%rY = pgrd%rY * degrees_to_radians
1578
1579 endif
1580
1581 iretval = pj_init_and_transform(csfromproj4//c_null_char, &
1582 cstoproj4//c_null_char, &
1583 __file__//c_null_char, &
1584 __line__, &
1585 int(pgrd%iNumGridCells, c_long), pgrd%rX, pgrd%rY)
1586
1587 call grid_checkforproj4error(iretval, sfromproj4, stoproj4)
1588
1589
1590 !> If the coordinates have been converted TO latlon, convert back to degrees
1591
1592 if( index(string=cstoproj4, substring="latlon") > 0 &
1593 .or. index(string=cstoproj4, substring="lonlat") > 0 ) then
1594
1595 pgrd%rX = pgrd%rX * radians_to_degrees
1596 pgrd%rY = pgrd%rY * radians_to_degrees
1597
1598 endif
1599
1600 ! now update the grid boundaries based on the transformed coordinate values
1601 pgrd%rGridCellSize = ( maxval(pgrd%rX) - minval(pgrd%rX) ) &
1602 / real(pgrd%iNX - 1, c_double)
1603
1604 pgrd%rX0 = minval(pgrd%rX) - pgrd%rGridCellSize / 2_c_float
1605 pgrd%rX1 = maxval(pgrd%rX) + pgrd%rGridCellSize / 2_c_float
1606 pgrd%rY0 = minval(pgrd%rY) - pgrd%rGridCellSize / 2_c_float
1607 pgrd%rY1 = maxval(pgrd%rY) + pgrd%rGridCellSize / 2_c_float
1608
1609 ! finally, change the projection string to reflect the new coordinate system
1610 pgrd%sPROJ4_string = trim(stoproj4)
1611
1612 ! at this point, what we have is the same old values that were read into the
1613 ! input grid, but a new set of coordinate and boundary definitions that reference
1614 ! the SWB grid projection.
1615
1616end subroutine grid_transform
1617
1618!--------------------------------------------------------------------------
1619
1620subroutine grid_checkintegergridvalues(pGrd, sFilename)
1621
1622 type ( GENERAL_GRID_T ),pointer :: pGrd
1623 character (len=*), intent(in) :: sFilename
1624
1625 ! [ LOCALS ]
1626 integer (c_int) :: iRunningSum
1627 integer (c_int) :: iIndex
1628 integer (c_int) :: iCount
1629 integer (c_int) :: iRecord
1630 character (len=10) :: sBuf0
1631 character (len=14) :: sBuf1
1632 character (len=10) :: sBuf2
1633 character (len=40) :: sBuf3
1634
1635 irunningsum = 0
1636 irecord = 0
1637
1638! call LOGS%set_echo(FALSE)
1639! call LOGS%write("### Summary of integer grid data values for file "//dquote(sFilename)//" ###", &
1640! iLogLevel=LOG_DEBUG, iLinesBefore=1, iLinesAfter=1 )
1641
1642! call LOGS%write("number | count | value ", iLogLevel=LOG_DEBUG )
1643! call LOGS%write("---------- | -------------- | ----------", iLogLevel=LOG_DEBUG )
1644! do iIndex=0,maxval(pGrd%iData)
1645! iCount=COUNT( pGrd%iData==iIndex )
1646! if ( iCount > 0 ) then
1647! iRecord = iRecord + 1
1648! iRunningSum = iRunningSum + iCount
1649! write (sBuf0, fmt="(i10)") iRecord
1650! write (sBuf1, fmt="(i14)") iCount
1651! write (sBuf2, fmt="(i10)") iIndex
1652! write (sBuf3, fmt="(a10,' | ', a14,' | ',a10)") adjustl(sBuf0), adjustl(sBuf1), adjustl(sBuf2)
1653! call LOGS%write( sBuf3, iLogLevel=LOG_DEBUG )
1654! end if
1655! end do
1656
1657! call LOGS%write(" Total number of grid cells with value NODATA: " &
1658! //asCharacter( COUNT(pGrd%iData == pGrd%iNoDataValue ) ), iLinesBefore=1, iLogLevel=LOG_ALL )
1659
1660! call LOGS%write(" Total number of grid cells: "//asCharacter( size(pGrd%iData) ), iLogLevel=LOG_ALL )
1661
1662! call LOGS%write(" Total number of grid cells with value >= 0: "//asCharacter(iRunningSum), iLogLevel=LOG_ALL )
1663
1664
1665! if (size(pGrd%iData) /= iRunningSum) then
1666! call LOGS%write(repeat("*",80), iLogLevel=LOG_ALL)
1667! call LOGS%write("Possible illegal or missing values in integer grid file: "//trim(sFileName), iLogLevel=LOG_ALL)
1668! call LOGS%write(repeat("*",80), iLogLevel=LOG_ALL)
1669! endif
1670
1671end subroutine grid_checkintegergridvalues
1672
1673!--------------------------------------------------------------------------
1674
1675subroutine grid_checkforproj4error(iRetVal, sFromPROJ4, sToPROJ4)
1676
1677 ! [ ARGUMENTS ]
1678 integer (c_int) :: iretval
1679 character (len=*) :: sfromproj4
1680 character (len=*) :: stoproj4
1681
1682 ! [ LOCALS ]
1683 integer (c_int) :: i
1684 logical (c_bool) :: lfound
1685 character (len=256) :: serrormessage
1686
1687 serrormessage = ""
1688
1689 if (iretval /= 0) then
1690
1691 write(serrormessage,fmt="(a,a,a,a,a)") &
1692 "There was an error transforming a grid from this projection:~", &
1693 dquote(sfromproj4), " to:~", &
1694 dquote(stoproj4)
1695
1696 do i=1,49
1697 if(iretval == terrormessage(i)%iErrorNumber) then
1698 serrormessage = serrormessage &
1699 //"~ PROJ4 Error number "// ascharacter(iretval)//" reported.~" &
1700 //" ==> error description: "//trim(terrormessage(i)%cErrorMessageText)
1701 lfound = true
1702 exit
1703 endif
1704 enddo
1705
1706 call assert(false, trim(serrormessage), __file__, __line__)
1707
1708 endif
1709
1710end subroutine grid_checkforproj4error
1711
1712!!***
1713!--------------------------------------------------------------------------
1714!!****f* grid/grid_Interpolate
1715! NAME
1716! grid_Interpolate - Returns an interpolated table value given a row
1717! and column position
1718!
1719
1720! SYNOPSIS
1721! Returns an interpolated table value given a row and column position
1722!
1723! INPUTS
1724! pGrd - Pointer to a data grid
1725! rXval - X value for which we want to interpolate table data.
1726! rYval - Y value for which we want to interpolate table data.
1727!
1728! OUTPUTS
1729! rValue - Value interpolated from input grid.
1730!
1731! NOTES
1732! The grid lookup functions were designed specifically to work with the
1733! Thornthwaite-Mather soil moisture retention table. The functions will
1734! work with any data that are similarly arranged.
1735!
1736! SOURCE
1737
1738function grid_interpolate(pGrd,rXval,rYval) result ( rValue )
1739 !! Interpolates values from the grid 'grd' for the row position 'yval' and
1740 !! the column position 'xval'. Assumes that the row and column spacing
1741 !! are constant. Applicable only to DATATYPE_REAL grids.
1742 ! [ ARGUMENTS ]
1743 type ( general_grid_t ),pointer :: pgrd
1744 real (c_float),intent(in) :: rxval,ryval
1745 ! [ RETURN VALUE ]
1746 real (c_float) :: rvalue
1747 ! [ LOCALS ]
1748 integer (c_int) :: ib,jb,ia,ja
1749 real (c_float) :: ylocal,u,v
1750
1751 call grid_lookupcolumn(pgrd,rxval,ib,ia,u)
1752
1753 call assert (ib>0 .and. ia>0 .and. ib <= ubound(pgrd%rData,1) &
1754 .and. ia <= ubound(pgrd%rData,1), &
1755 "Internal programming error: illegal bounds caught~requested column value " &
1756 //trim( ascharacter(rxval, fmt_string="F0.3")) &
1757 //" out of range", __file__, __line__)
1758
1759 ! In some cases, when things really dry out, the y value
1760 ! goes out of range - enforce bounds.
1761 if ( ryval < pgrd%rY0 ) then
1762 ylocal = pgrd%rY0
1763 else if ( ryval > pgrd%rY1 ) then
1764 ylocal = pgrd%rY1
1765 else
1766 ylocal = ryval
1767 end if
1768
1769 call grid_lookuprow(pgrd=pgrd, &
1770 ryval=ylocal, &
1771 ibefore=jb, &
1772 iafter=ja, &
1773 rfrac=v)
1774
1775 call assert (jb>0 .and. ja>0 .and. jb <= ubound(pgrd%rData,2) &
1776 .and. ja <= ubound(pgrd%rData,2), &
1777 "Internal programming error: illegal bounds caught~requested row value " &
1778 //trim( ascharacter(rxval, fmt_string="F0.3")) &
1779 //" out of range", __file__, __line__)
1780
1781 rvalue = ( 1.0_c_float -u) * ( 1.0_c_float -v) * pgrd%rData(ib,jb) + &
1782 u * ( 1.0_c_float -v) * pgrd%rData(ib,ja) + &
1783 ( 1.0_c_float -u) * v * pgrd%rData(ia,jb) + &
1784 u * v * pgrd%rData(ia,ja)
1785
1786end function grid_interpolate
1787!!***
1788!--------------------------------------------------------------------------
1789!!****f* grid/grid_SearchColumn
1790! NAME
1791! grid_SearchColumn - Returns the value of table Y value given the X value
1792! and the table Z value
1793!
1794! SYNOPSIS
1795! Searches in the y-direction, given the x value and the table value
1796! z, and returns the value of y. Assumes that the row and
1797! column spacing are constant. The search begins at the _top_l of the
1798! table (y=pGrd%rY1). Applicable only to DATATYPE_REAL grids.
1799! Parameter 'rNoData' is the missing value code for the grid.
1800
1801! INPUTS
1802! pGrd - Pointer to a data grid
1803! rXval - X value corresponding to one or more rows of the table.
1804! rZval - Table value to scan for.
1805! rNoData - Real value assigned to "No Data" cells.
1806!
1807! OUTPUTS
1808! rYval - Returned value associated with the input x and z values.
1809!
1810! NOTES
1811! The grid lookup functions were designed specifically to work with the
1812! Thornthwaite-Mather soil moisture retention table. The functions will
1813! work with any data that are similarly arranged.
1814!
1815! If X=1.00 and Z=0.90, for example, and the following table were queried,
1816! the subroutine will return: rValue=0.1
1817!
1818! X=Max SM Capacity
1819!
1820! 0.50 1.00 1.50
1821! _________________l
1822! 0.0 | 0.50 1.00 1.50
1823! Y=APWL 0.1 | 0.45 0.90 1.40
1824! 0.2 | 0.40 0.80 1.30
1825!
1826! SOURCE
1827
1828function grid_searchcolumn(pGrd,rXval,rZval,rNoData) result ( rValue )
1829 !! Searches in the y-direction, given the x value 'xval', for the value
1830 !! 'zval', and returns the value of y in 'ry'. Assumes that the row and
1831 !! column spacing are constant. The search begins at the _top_l of the
1832 !! table (y=pGrd%rY1). Applicable only to DATATYPE_REAL grids.
1833 !! Parameter 'rmv' is the missing value code for the grid.
1834 ! [ ARGUMENTS ]
1835 type ( general_grid_t ),pointer :: pgrd
1836 real (c_float),intent(in) :: rxval,rzval,rnodata
1837 ! [ RETURN VALUE ]
1838 real (c_float) :: rvalue
1839 ! [ LOCALS ]
1840 integer (c_int) :: ib,ia,istat,irow
1841 real (c_float) :: u,v,frac,rprev
1842 real (c_float),dimension(:),allocatable :: rcol
1843 character (len=128) :: buf
1844
1845 ! allocate space for all ROWS of data that come from a single COLUMN
1846 allocate ( rcol(pgrd%iNY), stat=istat )
1847 call assert( LOGICAL(istat==0,c_bool), &
1848 "Couldn't allocate temporary array" )
1849
1850 call grid_lookupcolumn(pgrd=pgrd, &
1851 rxval=rxval, &
1852 ibefore=ib, &
1853 iafter=ia, &
1854 rfrac=u)
1855
1856 call assert (ib>0 .and. ia>0 .and. ib <= ubound(pgrd%rData,1) &
1857 .and. ia <= ubound(pgrd%rData,1), &
1858 "Internal programming error: requested X value " &
1859 //trim( ascharacter(rxval, fmt_string="F0.3")) &
1860 //" out of range", __file__, __line__)
1861
1862 call assert(ubound(rcol,1) == ubound(pgrd%rData,2), &
1863 "Internal programming error: upper bound of rCol /= upper " &
1864 //"bound of first array element of pGrd%rData~" &
1865 //"ubound(rCol)="//trim( ascharacter(ubound(rcol,1))) &
1866 //"~ubound(pGrd%rData)="//trim( ascharacter(ubound(pgrd%rData,1))), &
1867 __file__, __line__)
1868
1869 ! interpolate the column of values based on the columns of values
1870 ! that bracket the real value rXval
1871! rCol = u*pGrd%rData(:,ia) + ( 1.0_c_float -u)*pGrd%rData(:,ib)
1872 rcol = u * pgrd%rData(ia,:) + ( 1.0_c_float -u)*pgrd%rData(ib,:)
1873 ! Fix missing values
1874
1875 do irow=1,pgrd%iNY
1876 if ( pgrd%rData(ia, irow) == rnodata .and. pgrd%rData(ib, irow) == rnodata ) then
1877 rcol(irow) = rnodata
1878 else if ( pgrd%rData(ib,irow) == rnodata .and. u>0.9_c_float ) then
1879 rcol(irow) = pgrd%rData(ia,irow)
1880 else if ( pgrd%rData(ia,irow) == rnodata .and. u<0.1_c_float ) then
1881 rcol(irow) = pgrd%rData(ib,irow)
1882 else if ( pgrd%rData(ia,irow) == rnodata .or. pgrd%rData(ib,irow) == rnodata ) then
1883 rcol(irow) = rnodata
1884 end if
1885 end do
1886
1887 ! Search for the specified value, assuming that the data vary monotonically
1888 ! Skip MVs
1889 rprev = rnodata
1890 rvalue = rnodata
1891 do irow=1,pgrd%iNY
1892 if ( rcol(irow) == rnodata ) then
1893 if ( rprev == rnodata ) then
1894 ! keep skipping missing values...
1895 continue
1896 else
1897 ! end of the line -- we didn't find the value, so return the limit
1898 v = real(irow-1) / real(pgrd%iNY-1)
1899 rvalue = ( 1.0_c_float -v)*pgrd%rY0 + v*pgrd%rY1
1900 exit
1901 endif
1902 else
1903 if ( rprev == rnodata ) then
1904 rprev = rcol(irow)
1905 else
1906 if ( sign( 1.0_c_float ,rcol(irow)-rzval) /= sign( 1.0_c_float ,rprev-rzval) ) then
1907 ! Found it!
1908 frac = (rzval-rcol(irow-1)) / (rcol(irow)-rcol(irow-1))
1909 ! Note: it's 'iRow-2' in the next expr to account for one-based indexing
1910 v = ( real(irow-2)+frac ) / real(pgrd%iNY-1)
1911 rvalue = v*pgrd%rY0 + ( 1.0_c_float -v)*pgrd%rY1
1912 exit
1913 else if ( rzval < rprev .and. rcol(irow) > rprev ) then
1914 ! does not exist, choose the limit
1915 v = real(irow-2) / real(pgrd%iNY-1)
1916 rvalue = v*pgrd%rY0 + ( 1.0_c_float -v)*pgrd%rY1
1917 exit
1918 else if ( rzval > rprev .and. rcol(irow) < rprev ) then
1919 v = real(irow-2) / real(pgrd%iNY-1)
1920 rvalue = v*pgrd%rY0 + ( 1.0_c_float -v)*pgrd%rY1
1921 exit
1922 else if ( irow == pgrd%iNY ) then
1923 v = 1.0_c_float
1924 rvalue = v*pgrd%rY0 + ( 1.0_c_float -v)*pgrd%rY1
1925 exit
1926 end if
1927 end if
1928 end if
1929 end do
1930
1931 deallocate ( rcol )
1932
1933end function grid_searchcolumn
1934!!***
1935!--------------------------------------------------------------------------
1936!!****f* grid/grid_LookupReal
1937! NAME
1938! grid_LookupReal - Returns the grid cell value nearest to that for a
1939! given a row and column position
1940!
1941! SYNOPSIS
1942! Returns the grid cell value nearest to that for a given a row and
1943! column position. No interpolation is performed.
1944!
1945! INPUTS
1946! pGrd - Pointer to a data grid
1947! rXval - X value for which we want to interpolate table data.
1948! rYval - Y value for which we want to interpolate table data.
1949!
1950! OUTPUTS
1951! rValue - Value of the grid cell nearest the given row, column combination.
1952!
1953! NOTES
1954! The grid lookup functions were designed specifically to work with the
1955! Thornthwaite-Mather soil moisture retention table. The functions will
1956! work with any data that are similarly arranged.
1957!
1958! SOURCE
1959
1960function grid_lookupreal(pGrd,rXval,rYval) result(rValue)
1961 !! Returns the grid value for the cell containing (rXval,rYval)
1962 type ( general_grid_t ),pointer :: pgrd
1963 real (c_float),intent(in) :: rxval,ryval
1964 real (c_float) :: rvalue
1965 integer (c_int) :: icol,irow
1966
1967 irow = int( pgrd%iNY * (pgrd%rY0 - ryval) / (pgrd%rY0 - pgrd%rY1) ) + 1
1968 if ( irow > pgrd%iNY ) irow = pgrd%iNY
1969 icol = int( pgrd%iNX * (rxval - pgrd%rX0) / (pgrd%rX1 - pgrd%rX0) ) + 1
1970 if ( icol > pgrd%iNX ) icol = pgrd%iNX
1971 rvalue = pgrd%rData(icol,irow)
1972
1973end function grid_lookupreal
1974!!***
1975
1976!----------------------------------------------------------------------
1977
1978function grid_coordinatesfallinsidegrid(pGrd, rX, rY) result(coordinates_fall_inside_grid)
1979
1980 type ( general_grid_t ),pointer :: pgrd
1981 real (c_double) :: rx
1982 real (c_double) :: ry
1983 logical (c_bool) :: coordinates_fall_inside_grid
1984
1985 coordinates_fall_inside_grid = ( ( rx >= pgrd%rX0 ) .and. ( rx <= pgrd%rX1 ) &
1986 .and. ( ry >= pgrd%rY0 ) .and. ( ry <= pgrd%rY1 ) )
1987
1988end function grid_coordinatesfallinsidegrid
1989
1990!----------------------------------------------------------------------
1991
1992function grid_rowcolfallsinsidegrid(pGrd, iRow, iCol) result(row_col_falls_inside_grid)
1993
1994 type ( general_grid_t ),pointer :: pgrd
1995 integer (c_int) :: irow
1996 integer (c_int) :: icol
1997 logical (c_bool) :: row_col_falls_inside_grid
1998
1999 row_col_falls_inside_grid = ( ( irow >= 1 ) .and. ( irow <= pgrd%iNY ) &
2000 .and. ( icol >= 1 ) .and. ( icol <= pgrd%iNX ) )
2001
2002end function grid_rowcolfallsinsidegrid
2003
2004!----------------------------------------------------------------------
2005
2006function grid_getgridcolnum(pGrd,rX) result(iColumnNumber)
2007
2008 type ( general_grid_t ),pointer :: pgrd
2009 real (c_double) :: rx
2010 integer (c_int) :: icolumnnumber
2011
2012 !! this only works if the data are in the originally supplied projection. If the coordinates have been
2013 !! transformed, there is no guarantee that the assumption used in this calculation will hold.
2014 icolumnnumber = nint(real(pgrd%iNX, c_double) &
2015 * ( rx - pgrd%rX0 ) / (pgrd%rX1 - pgrd%rX0) + 0.5_c_double, c_int)
2016
2017! * ( rX - pGrd%rX0 ) / (pGrd%rX1 - pGrd%rX0) + 0.5_c_double, c_int)
2018
2019 ! print *, "iColumnNumber = ", iColumnNumber
2020 ! print *, "calc: ", real(pGrd%iNX, c_double) &
2021 ! * ( rX - pGrd%rX0 ) / (pGrd%rX1 - pGrd%rX0) + 0.5_c_double
2022 ! print *, "numerator: ", real(pGrd%iNX, c_double) * ( rX - pGrd%rX0 )
2023 ! print *, "numerator (LHS): ", real(pGrd%iNX, c_double)
2024 ! print *, "numerator (RHS): ", rX, pGrd%rX0, ( rX - pGrd%rX0 )
2025 ! print *, "denominator: ", (pGrd%rX1 - pGrd%rX0)
2026 !
2027! if ( iColumnNumber < 1 .or. iColumnNumber > pGrd%iNX ) then
2028! call grid_DumpGridExtent(pGrd)
2029! write(*, fmt="(a)") "was attempting to find column associated with X: "//trim(asCharacter(rX))
2030! call assert(FALSE, "INTERNAL PROGRAMMING ERROR: Column number out of bounds (value: " &
2031! //trim( asCharacter(iColumnNumber))//")", &
2032! __FILE__, __LINE__)
2033! endif
2034
2035end function grid_getgridcolnum
2036
2037!----------------------------------------------------------------------
2038
2039function grid_getgridrownum(pGrd,rY) result(iRowNumber)
2040
2041 type ( general_grid_t ),pointer :: pgrd
2042 real (c_double) :: ry
2043 integer (c_int) :: irownumber
2044
2045 irownumber = pgrd%iNY - nint(real(pgrd%iNY, c_double) &
2046 * ( ry - pgrd%rY0 ) / (pgrd%rY1 - pgrd%rY0) - 0.5_c_double, c_int)
2047
2048! if ( iRowNumber < 1 .or. iRowNumber > pGrd%iNY ) then
2049! call grid_DumpGridExtent(pGrd)
2050! write(*, fmt="(a)") "was attempting to find row associated with Y: "//trim(asCharacter(rY))
2051! call assert(FALSE, "INTERNAL PROGRAMMING ERROR: Row number out of bounds (value: " &
2052! //trim( asCharacter(iRowNumber))//")", &
2053! __FILE__, __LINE__)
2054! endif
2055
2056end function grid_getgridrownum
2057
2058!----------------------------------------------------------------------
2059
2060function grid_getgridcolrownum(pGrd, rX, rY) result(iColRow)
2061
2062 type ( general_grid_t ),pointer :: pgrd
2063 real (c_double) :: rx, ry
2064 integer (c_int), dimension(2) :: icolrow
2065
2066 ! [ LOCALS ]
2067 real (c_float) :: rdist, rmindistance, rdist2
2068
2069 integer (c_int), save :: ilastcolnum, ilastrownum
2070 integer (c_int) :: istartingcolnum, istartingrownum
2071 integer (c_int) :: irowboundlower, irowboundupper
2072 integer (c_int) :: icolboundlower, icolboundupper
2073 integer (c_int) :: icol, irow
2074 integer (c_int) :: icandidaterow, icandidatecol
2075 integer (c_int) :: istat
2076 logical (c_bool) :: lchanged
2077
2078 ! this will get us close to where we need to be; however, since the
2079 ! transformed grid may contain significant nonlinearity between
2080 ! adjacent grid coordinates, we need to do a more thorough search
2081 ! in the neighborhood of these points
2082
2083 istartingcolnum = grid_getgridcolnum(pgrd,rx)
2084 istartingrownum = grid_getgridrownum(pgrd,ry)
2085
2086 if (istartingcolnum > 0 .and. istartingcolnum <= pgrd%iNX &
2087 .and. istartingrownum > 0 .and. istartingrownum <= pgrd%iNY) then
2088 rdist = hypot(pgrd%rX(istartingcolnum,istartingrownum) - rx, &
2089 pgrd%rY(istartingcolnum,istartingrownum) - ry)
2090
2091 ! need to ensure that the leftover column and row numbers from
2092 ! perhaps an entirely different grid are not used as indices
2093 if (ilastcolnum > 0 .and. ilastcolnum <= pgrd%iNX &
2094 .and. ilastrownum > 0 .and. ilastrownum <= pgrd%iNY) then
2095 rdist2 = hypot(pgrd%rX(ilastcolnum,ilastrownum) - rx, &
2096 pgrd%rY(ilastcolnum,ilastrownum) - ry)
2097 else
2098 rdist2 = rbigval
2099 endif
2100
2101 if (rdist > rdist2) then
2102 icandidaterow = ilastrownum
2103 icandidatecol = ilastcolnum
2104 else
2105 icandidaterow = istartingrownum
2106 icandidatecol = istartingcolnum
2107 endif
2108
2109! print *, 'filename: ', trim(pGrd%sFilename)
2110! print *, 'start of iteration: ', iStartingColNum, iStartingRowNum, rDist, rDist2, iCandidateCol, iCandidateRow
2111
2112 rmindistance = rbigval
2113
2114 do
2115
2116 !> need to ensure that whatever bound is calculated
2117 !> is within the declared array bounds or we get a segfault
2118 !iRowBoundLower = min(max( 1, iCandidateRow - 1), pGrd%iNY)
2119 irowboundlower = clip(value=icandidaterow - 1, minval=1, maxval=pgrd%iNY)
2120 irowboundupper = clip(value=icandidaterow + 1, minval=1, maxval=pgrd%iNY)
2121
2122 !iColBoundLower = min(max( 1, iCandidateCol - 1), pGrd%iNX)
2123 icolboundlower = clip(value=icandidatecol - 1, minval=1, maxval=pgrd%iNX)
2124 icolboundupper = clip(value=icandidatecol + 1, minval=1, maxval=pgrd%iNX)
2125
2126 lchanged = false
2127
2128 do irow=irowboundlower,irowboundupper
2129 do icol=icolboundlower,icolboundupper
2130
2131 rdist = hypot(pgrd%rX(icol,irow) - rx, pgrd%rY(icol,irow) - ry)
2132
2133 if (rdist < rmindistance ) then
2134
2135 rmindistance = rdist
2136 icandidatecol = icol
2137 icandidaterow = irow
2138 lchanged = true
2139
2140 endif
2141
2142 enddo
2143 enddo
2144
2145! print *, 'iterating: rdist: ', rDist, ' cand X: ', pGrd%rX(iCandidateCol,iCandidateRow), ' X: ', rX, ' cand Y: ', pGrd%rY(iCandidateCol,iCandidateRow), ' Y: ', rY
2146
2147 if (.not. lchanged ) exit
2148
2149 enddo
2150
2151 ilastcolnum = icandidatecol
2152 ilastrownum = icandidaterow
2153
2154 icolrow(column) = icandidatecol
2155 icolrow(row) = icandidaterow
2156
2157 else
2158
2159 icolrow(column) = istartingcolnum
2160 icolrow(row) = istartingrownum
2161 endif
2162
2163end function grid_getgridcolrownum
2164
2165!----------------------------------------------------------------------
2166subroutine grid_set_nodata_value( pGrd, iValue, fValue )
2167 type ( general_grid_t ),pointer :: pgrd
2168 integer (c_int), intent(in), optional :: ivalue
2169 real (c_float), intent(in), optional :: fvalue
2170
2171 if ( present( ivalue ) ) pgrd%iNoDataValue = ivalue
2172 if ( present( fvalue ) ) pgrd%rNoDataValue = fvalue
2173end subroutine grid_set_nodata_value
2174!--------------------------------------------------------------------------------------------------
2175
2176function grid_getgridx(pGrd,iColumnNumber) result(rX)
2177
2178 type ( general_grid_t ),pointer :: pgrd
2179 real (c_double) :: rx
2180 integer (c_int) :: icolumnnumber
2181
2182 rx = pgrd%rX0 + pgrd%rGridCellSize * (real(icolumnnumber, c_double) - 0.5_c_double )
2183
2184end function grid_getgridx
2185
2186!----------------------------------------------------------------------
2187
2188function grid_getgridy(pGrd,iRowNumber) result(rY)
2189
2190 type ( general_grid_t ),pointer :: pgrd
2191 real (c_double) :: ry
2192 integer (c_int) :: irownumber
2193
2194 ry = pgrd%rY1 &
2195 - pgrd%rGridCellSize * (real(irownumber, c_double) - 0.5_c_double )
2196
2197end function grid_getgridy
2198
2199!----------------------------------------------------------------------
2200
2201subroutine grid_populatexy(pGrd, rX, rY)
2202
2203 ! [ ARGUMENTS ]
2204 type ( general_grid_t ),pointer :: pgrd
2205 real (c_double), optional :: rx(:)
2206 real (c_double), optional :: ry(:)
2207
2208 ! [ LOCALS ]
2209 integer (c_int) :: icol,irow
2210 integer (c_int) :: istat
2211 integer (c_int) :: nx, ny
2212
2213 if ( .not. allocated(pgrd%rX) ) then
2214
2215 ALLOCATE (pgrd%rX(pgrd%iNX, pgrd%iNY), stat=istat)
2216 call assert( istat == 0, &
2217 "Could not allocate memory for x-coordinates within grid data structure", &
2218 __file__, __line__)
2219 endif
2220
2221 if ( .not. allocated(pgrd%rY) ) then
2222 ALLOCATE (pgrd%rY(pgrd%iNX, pgrd%iNY), stat=istat)
2223 call assert( istat == 0, &
2224 "Could not allocate memory for y-coordinates within grid data structure", &
2225 __file__, __line__)
2226 endif
2227
2228 if (present(rx) .and. present(ry)) then
2229
2230 nx = size(rx)
2231 ny = size(ry)
2232
2233 call assert( nx == pgrd%iNX, &
2234 "Internal programming error: number of 'X' values does not match number of x values in grid", &
2235 __file__, __line__)
2236
2237 call assert( ny == pgrd%iNY, &
2238 "Internal programming error: number of 'X' values does not match number of x values in grid", &
2239 __file__, __line__)
2240
2241 do irow=1,pgrd%iNY
2242 do icol=1,pgrd%iNX
2243 pgrd%rX(icol, irow) = rx(icol)
2244 pgrd%rY(icol, irow) = ry(irow)
2245 enddo
2246 enddo
2247
2248 else
2249
2250 do irow=1,pgrd%iNY
2251 do icol=1,pgrd%iNX
2252 pgrd%rX(icol, irow) = grid_getgridx(pgrd, icol)
2253 pgrd%rY(icol, irow) = grid_getgridy(pgrd, irow)
2254 enddo
2255 enddo
2256
2257 endif
2258
2259 ! print *, 'grid_PopulateXY: obtained grid coordinates, routine ',__FILE__,' line ', __LINE__
2260 ! print *, '(col=1, row=1): ', pGrd%rX(1,1), pGrd%rY(1,1)
2261 ! print *, '(col=1, row=', pGrd%iNY,'): ', pGrd%rX(1, pGrd%iNY),pGrd%rY(1, pGrd%iNY)
2262 ! print *, '(col=',pGrd%iNX,', row=1): ', pGrd%rX(pGrd%iNX, 1),pGrd%rY(pGrd%iNX, 1)
2263 ! print *, '(col=',pGrd%iNX,', row=', pGrd%iNY,'): ', pGrd%rX(pGrd%iNX, pGrd%iNY),pGrd%rY(pGrd%iNX, pGrd%iNY)
2264
2265
2266end subroutine grid_populatexy
2267
2268!----------------------------------------------------------------------
2269
2270subroutine grid_dumpgridextent(pGrd)
2271
2272 type ( general_grid_t ),pointer :: pgrd
2273
2274 if ( associated( pgrd) ) then
2275
2276 call logs%write("---------------------------------------------------")
2277 call logs%write("GRID DETAILS:")
2278 call logs%write("---------------------------------------------------")
2279 call logs%write("file: "//dquote(pgrd%sFilename) )
2280 call logs%write("nx: "//trim( ascharacter(pgrd%iNX) ) )
2281 call logs%write("ny: "//trim( ascharacter(pgrd%iNY) ) )
2282 call logs%write("cellsize: "//trim(ascharacter(pgrd%rGridCellSize) ) )
2283 call logs%write("X0: "//trim(ascharacter(pgrd%rX0) ) )
2284 call logs%write("Y0: "//trim(ascharacter(pgrd%rY0) ) )
2285 call logs%write("X1: "//trim(ascharacter(pgrd%rX1) ) )
2286 call logs%write("Y1: "//trim(ascharacter(pgrd%rY1) ) )
2287 call logs%write("Type: "//trim(ascharacter(pgrd%iDataType) ) )
2288 call logs%write("PROJ4 string: "//dquote(pgrd%sPROJ4_string) )
2289 call logs%write("---------------------------------------------------")
2290
2291 else
2292
2293 call logs%write("---------------------------------------------------")
2294 call logs%write(" ** grid not allocated ** ")
2295 call logs%write("---------------------------------------------------")
2296
2297 endif
2298
2299
2300end subroutine grid_dumpgridextent
2301
2302!----------------------------------------------------------------------
2303function grid_gridtopoint_int(pGrdFrom, pGrdTo, iCol, iRow) result(iValue)
2304
2305 type ( general_grid_t ), pointer :: pgrdfrom
2306 type ( general_grid_t ), pointer :: pgrdto
2307 integer (c_int), intent(in) :: icol
2308 integer (c_int), intent(in) :: irow
2309 integer (c_int) :: ivalue
2310
2311 ! [ LOCALS ]
2312 integer (c_int), dimension(2) :: icolrow
2313 integer (c_int) :: isrccol, isrcrow
2314 integer (c_int) :: ispread
2315
2316 ! must ensure that there are coordinates associated with the "to" grid...
2317 ! by default, these are left unpopulated during a "normal" swb run
2318 if( .not. allocated(pgrdto%rX) ) call grid_populatexy(pgrdto)
2319 if( .not. allocated(pgrdfrom%rX) ) call grid_populatexy(pgrdfrom)
2320
2321 icolrow = grid_getgridcolrownum(pgrd=pgrdfrom, &
2322 rx=real(pgrdto%rX(icol, irow), c_double), &
2323 ry=real(pgrdto%rY(icol, irow), c_double))
2324
2325 if ( icolrow(column) < 1 .or. icolrow(column) > pgrdfrom%iNX ) &
2326 call die( "Illegal column number supplied: "//trim(ascharacter(icolrow(column))), &
2327 __file__, __line__)
2328
2329 if ( icolrow(row) < 1 .or. icolrow(row) > pgrdfrom%iNY ) &
2330 call die( "Illegal row number supplied: "//trim(ascharacter(icolrow(row))), &
2331 __file__, __line__)
2332
2333 !> @todo check the logic here...intent is to ensure that the majority filter is
2334 !! searching an area that adequately covers corresponding gridcell areas
2335 ispread = max(1, nint(pgrdto%rGridCellSize / pgrdfrom%rGridCellSize / 2.))
2336
2337 ivalue = grid_majorityfilter_int(pgrdfrom=pgrdfrom, &
2338 itargetcol=icolrow(column), &
2339 itargetrow=icolrow(row), &
2340 inodatavalue=pgrdfrom%iNoDataValue, &
2341 ispread=ispread)
2342
2343
2344end function grid_gridtopoint_int
2345
2346
2347!----------------------------------------------------------------------
2348
2349subroutine grid_gridtogrid_int( pGrdFrom, pGrdTo, lUseMajorityFilter )
2350
2351 type ( general_grid_t ), pointer :: pgrdfrom ! pointer to source grid
2352 type ( general_grid_t ), pointer :: pgrdto ! pointer to destination grid
2353 logical (c_bool), intent(in) :: lusemajorityfilter
2354
2355 ! [ LOCALS ]
2356 integer (c_int) :: icol, irow
2357 integer (c_int), dimension(2) :: icolrow
2358 integer (c_int) :: isrccol, isrcrow
2359 integer (c_int) :: ispread
2360 real (c_float) :: fgridcellratio
2361
2362 ! must ensure that there are coordinates associated with the "to" grid...
2363 ! by default, these are left unpopulated during a "normal" swb run
2364 if(.not. allocated(pgrdto%rX) ) call grid_populatexy(pgrdto)
2365 if(.not. allocated(pgrdfrom%rX) ) call grid_populatexy(pgrdfrom)
2366
2367 fgridcellratio = pgrdto%rGridCellSize / pgrdfrom%rGridCellSize
2368
2369 ! Allow USER to trigger whether the majority filter is employed or not
2370 if ( lusemajorityfilter ) then
2371
2372 call logs%write( "** Majority filter in use for data from grid file "//dquote(pgrdfrom%sFilename), lecho=true )
2373
2374 ispread = max( 1, nint( fgridcellratio / 2.0_c_float ) )
2375
2376 do irow=1,pgrdto%iNY
2377 do icol=1,pgrdto%iNX
2378
2379 icolrow = grid_getgridcolrownum(pgrd=pgrdfrom, &
2380 rx=real(pgrdto%rX(icol, irow), c_double), &
2381 ry=real(pgrdto%rY(icol, irow), c_double) )
2382
2383 call assert(icolrow(column) > 0 .and. icolrow(column) <= pgrdfrom%iNX, &
2384 "Illegal column number supplied: "//trim(ascharacter(icolrow(column))), &
2385 __file__, __line__)
2386
2387 call assert(icolrow(row) > 0 .and. icolrow(row) <= pgrdfrom%iNY, &
2388 "Illegal row number supplied: "//trim(ascharacter(icolrow(row))), &
2389 __file__, __line__)
2390
2391 pgrdto%iData(icol,irow) = grid_majorityfilter_int( pgrdfrom=pgrdfrom, &
2392 itargetcol=icolrow(column), &
2393 itargetrow=icolrow(row), &
2394 inodatavalue=pgrdfrom%iNoDataValue, &
2395 ispread=ispread)
2396
2397 enddo
2398 enddo
2399
2400 else ! target grid resolution similar to or more dense than source grid resolution: NEAREST NEIGHBOR [DEFAULT]
2401
2402 do irow=1,pgrdto%iNY
2403 do icol=1,pgrdto%iNX
2404
2405 icolrow = grid_getgridcolrownum(pgrd=pgrdfrom, &
2406 rx=real(pgrdto%rX(icol, irow), c_double), &
2407 ry=real(pgrdto%rY(icol, irow), c_double))
2408
2409 call assert(icolrow(column) > 0 .and. icolrow(column) <= pgrdfrom%iNX, &
2410 "Illegal column number supplied: "//trim(ascharacter(icolrow(column))), &
2411 __file__, __line__)
2412
2413 call assert(icolrow(row) > 0 .and. icolrow(row) <= pgrdfrom%iNY, &
2414 "Illegal row number supplied: "//trim(ascharacter(icolrow(row))), &
2415 __file__, __line__)
2416
2417 pgrdto%iData(icol,irow) = pgrdfrom%iData( icolrow(column), icolrow(row) )
2418
2419 enddo
2420 enddo
2421
2422 endif
2423
2424end subroutine grid_gridtogrid_int
2425
2426!----------------------------------------------------------------------
2427
2428subroutine grid_gridtogrid_sgl(pGrdFrom, pGrdTo )
2429
2430 ! [ ARGUMENTS ]
2431 type ( general_grid_t ),pointer :: pgrdfrom ! pointer to source grid
2432 type ( general_grid_t ),pointer :: pgrdto ! pointer to destination grid
2433
2434 ! [ LOCALS ]
2435 integer (c_int), dimension(2) :: icolrow
2436 integer (c_int) :: icol, irow
2437 integer (c_int) :: isrccol, isrcrow
2438 real (c_float), dimension(3,3) :: rkernel
2439 logical :: printout
2440
2441 rkernel = 1.
2442 rkernel(2,2) = 8.
2443
2444 ! must ensure that there are coordinates associated with the "to" grid...
2445 ! by default, these are left unpopulated during a "normal" swb run
2446 if(.not. allocated(pgrdto%rX) ) call grid_populatexy(pgrdto)
2447 if(.not. allocated(pgrdfrom%rX) ) call grid_populatexy(pgrdfrom)
2448
2449!!! *$OMP PARALLEL DO ORDERED PRIVATE(iRow, iCol, iColRow)
2450
2451 do irow=1,pgrdto%iNY
2452 do icol=1,pgrdto%iNX
2453
2454 icolrow = grid_getgridcolrownum(pgrd=pgrdfrom, &
2455 rx=real(pgrdto%rX(icol, irow), c_double), &
2456 ry=real(pgrdto%rY(icol, irow), c_double))
2457
2458 if ( icolrow(column) < 1 .or. icolrow(column) > pgrdfrom%iNX) cycle
2459 if ( icolrow(row) < 1 .or. icolrow(row) > pgrdfrom%iNY) cycle
2460
2461 pgrdto%rData(icol,irow) = pgrdfrom%rData( icolrow(column), icolrow(row) )
2462
2463
2464 ! printout = (iRow==1 .and. iCol==1) .or. (iRow==1 .and. iCol==pGrdTo%iNX) &
2465 ! .or. (iRow==pGrdTo%iNY .and. iCol==1) .or. (iRow==pGrdTo%iNY .and. iCol==pGrdTo%iNX)
2466
2467 ! if (printout) then
2468 ! print *, trim(__FILE__), ': ', __LINE__
2469 ! print *, 'getting data from the file: ', trim(pGrdFrom%sFilename)
2470 ! print *, ' coord-to : ', pGrdTo%rX(iCol, iRow), pGrdTo%rY(iCol, iRow)
2471 ! print *, ' to: (iCol, iRow): ', iCol, iRow
2472 ! print *, ' from: (col, row): ', iColRow(COLUMN), iColRow(ROW)
2473 ! endif
2474
2475 enddo
2476 enddo
2477
2478!!! *$OMP END PARALLEL DO
2479
2480end subroutine grid_gridtogrid_sgl
2481
2482!----------------------------------------------------------------------
2483
2484function grid_majorityfilter_int(pGrdFrom, iTargetCol, &
2485 iTargetRow, iNoDataValue, iSpread) result(iRetVal)
2486
2487 type ( general_grid_t ),pointer :: pgrdfrom ! pointer to source grid
2488 integer (c_int) :: itargetcol ! column number of target cell
2489 integer (c_int) :: itargetrow ! row number of target cell
2490 integer (c_int) :: inodatavalue
2491 integer (c_int) :: iretval
2492 integer (c_int) :: ispread ! integer representing how many
2493 ! cells should be searched
2494
2495 ! [ LOCALS ]
2496 integer (c_int), dimension(625) :: ivalue
2497 integer (c_int), dimension(625) :: icount
2498 logical (c_bool) :: lmatch
2499 integer (c_int) :: irow
2500 integer (c_int) :: icol
2501 integer (c_int) :: i, isize, ilast, iindex, icellnum
2502
2503 ivalue = 0
2504 icount = 0
2505 ilast = 0
2506 icellnum = 0
2507
2508 ! need to set a reasonable upper bound on the spread to consider
2509 ispread = min(ispread, 12)
2510
2511 do irow=max(1,itargetrow - ispread), min(pgrdfrom%iNY, itargetrow + ispread)
2512 do icol=max(1,itargetcol - ispread), min(pgrdfrom%iNX, itargetcol + ispread)
2513!
2514 icellnum = icellnum + 1
2515!
2516 lmatch = false
2517!
2518 do i=1, ilast
2519 if (ivalue(i) == pgrdfrom%iData(icol, irow) ) then
2520 icount(i) = icount(i) + 1
2521 lmatch = true
2522 exit
2523 endif
2524 enddo
2525
2526 ! this integer value has not been found previously; add the value to the next position
2527 ! and increment counter
2528 if (.not. lmatch) then
2529 ilast = ilast + 1
2530 ivalue(ilast) = pgrdfrom%iData(icol, irow)
2531 icount(ilast) = icount(ilast) + 1
2532 lmatch = true
2533 endif
2534
2535 enddo
2536 enddo
2537
2538 iindex = maxloc(icount, dim=1)
2539 iretval = ivalue(iindex)
2540
2541end function grid_majorityfilter_int
2542
2543function grid_convolve_sgl( pGrdFrom, iTargetCol, &
2544 iTargetRow, rKernel, rNoDataValue) result(rRetVal)
2545
2546 type ( general_grid_t ),pointer :: pgrdfrom ! pointer to source grid
2547 integer (c_int) :: inx
2548 integer (c_int) :: iny
2549 integer (c_int) :: itargetcol ! column number of target cell
2550 integer (c_int) :: itargetrow ! row number of target cell
2551 real (c_float), dimension(:,:) :: rkernel
2552 real (c_float), optional :: rnodatavalue
2553 real (c_float) :: rretval
2554
2555 ! [ LOCALS ]
2556 integer (c_int) :: irowmin, irowmax
2557 integer (c_int) :: icolmin, icolmax
2558 integer (c_int) :: ikernelsize ! i.e. 3, 5, 7, 9, 11, etc.
2559 integer (c_int) :: iincvalue
2560 integer (c_int) :: icol, irow, irownum, icolnum
2561 real (c_float) :: rkernelsum
2562
2563 rkernelsum = rzero
2564 rretval = rzero
2565
2566 ikernelsize = size(rkernel, dim=1)
2567 iincvalue = (ikernelsize - 1) / 2
2568
2569 inx = ubound(pgrdfrom%rData,1)
2570 iny = ubound(pgrdfrom%rData,2)
2571
2572 irowmin = max(1,itargetrow - iincvalue)
2573 irowmax = min(iny, itargetrow + iincvalue)
2574
2575 icolmin = max(1,itargetcol - iincvalue)
2576 icolmax = min(iny, itargetcol + iincvalue)
2577
2578 if( (irowmax - irowmin + 1) /= ikernelsize &
2579 .or. (icolmax - icolmin + 1) /= ikernelsize ) then
2580
2581 ! This is a simple, but less than desirable way to treat cells that
2582 ! fall near the edge of the source grid. Ideally, the source grid will
2583 ! be greater than the target grid by a wide enough margin that this
2584 ! code will rarely be used.
2585 rretval = pgrdfrom%rData(itargetcol, itargetrow)
2586
2587 else
2588
2589 do icol=0,ikernelsize-1
2590 do irow=0,ikernelsize-1
2591
2592 irownum = irow + irowmin
2593 icolnum = icol + icolmin
2594
2595 if (irownum > iny .or. icolnum > inx) then
2596 cycle ! our calculated row or column number is outside of the
2597 ! bounds of the rValues array
2598 else
2599 rretval = rretval + pgrdfrom%rData(icolnum,irownum) * rkernel(icol+1, irow+1)
2600 rkernelsum = rkernelsum + rkernel(icol+1, irow+1)
2601 endif
2602 enddo
2603 enddo
2604
2605 if (rkernelsum > 0.) rretval = rretval / rkernelsum
2606
2607 endif
2608
2609end function grid_convolve_sgl
2610
2611end module grid
constants_and_conversions::clip
Definition constants_and_conversions.F90:180
fstring::as_character
Definition fstring.F90:55
grid::grid_create
Definition grid.F90:161
grid::pj_init_and_transform
interface to C code that provides a simplified entry point to PROJ4 capabilities: it has been modifie...
Definition grid.F90:42
constants_and_conversions
This module contains physical constants and convenience functions aimed at performing unit conversion...
Definition constants_and_conversions.F90:8
constants_and_conversions::true
logical(c_bool), parameter, public true
Definition constants_and_conversions.F90:38
constants_and_conversions::rzero
real(c_float), parameter, public rzero
Definition constants_and_conversions.F90:53
constants_and_conversions::datatype_real
integer(c_int), parameter datatype_real
Definition constants_and_conversions.F90:17
constants_and_conversions::false
logical(c_bool), parameter, public false
Definition constants_and_conversions.F90:39
constants_and_conversions::degrees_to_radians
real(c_double), parameter, public degrees_to_radians
Definition constants_and_conversions.F90:30
constants_and_conversions::rbigval
real(c_float), parameter, public rbigval
Definition constants_and_conversions.F90:41
constants_and_conversions::datatype_int
integer(c_int), parameter datatype_int
Definition constants_and_conversions.F90:15
constants_and_conversions::radians_to_degrees
real(c_double), parameter, public radians_to_degrees
Definition constants_and_conversions.F90:31
constants_and_conversions::datatype_double
integer(c_int), parameter datatype_double
Definition constants_and_conversions.F90:18
exceptions::die
subroutine, public die(smessage, smodule, iline, shints, scalledby, icalledbyline)
Definition exceptions.F90:139
grid
Provides support for input and output of gridded ASCII data, as well as for creation and destruction ...
Definition grid.F90:8
grid::grid_createcomplete
type(general_grid_t) function, pointer, public grid_createcomplete(inx, iny, rx0, ry0, rx1, ry1, idatatype)
Creates a grid of a specified type.
Definition grid.F90:211
grid::grid_getgridcolnum
integer(c_int) function, public grid_getgridcolnum(pgrd, rx)
Definition grid.F90:2007
grid::grid_coordinatesfallinsidegrid
logical(c_bool) function, public grid_coordinatesfallinsidegrid(pgrd, rx, ry)
Definition grid.F90:1979
grid::terrormessage
type(error_message_t), dimension(49) terrormessage
Definition grid.F90:99
grid::grid_populatexy
subroutine, public grid_populatexy(pgrd, rx, ry)
Definition grid.F90:2202
grid::grid_nodata_real
real(c_float), parameter, public grid_nodata_real
Definition grid.F90:21
grid::grid_searchcolumn
real(c_float) function grid_searchcolumn(pgrd, rxval, rzval, rnodata)
Definition grid.F90:1829
grid::output_grid_directory_name
character(len=256) output_grid_directory_name
Definition grid.F90:176
grid::grid_completelycover
logical(c_bool) function, public grid_completelycover(pbasegrd, pothergrd, rtolerance)
Definition grid.F90:1330
grid::eight_cells
integer(c_int), parameter eight_cells
Definition grid.F90:168
grid::grid_active_cell
integer(c_int), parameter, public grid_active_cell
Definition grid.F90:32
grid::grid_read
type(general_grid_t) function, pointer, public grid_read(sfilename, sfiletype, idatatype)
Definition grid.F90:447
grid::grid_lookupreal
real(c_float) function grid_lookupreal(pgrd, rxval, ryval)
Definition grid.F90:1961
grid::four_cells
integer(c_int), parameter four_cells
Definition grid.F90:167
grid::grid_convolve_sgl
real(c_float) function grid_convolve_sgl(pgrdfrom, itargetcol, itargetrow, rkernel, rnodatavalue)
Definition grid.F90:2545
grid::grid_writegrid
subroutine, public grid_writegrid(sfilename, pgrd, ioutputformat)
Definition grid.F90:1035
grid::nc_fill_float
real(c_float), parameter nc_fill_float
Definition grid.F90:34
grid::grid_dumpgridextent
subroutine, public grid_dumpgridextent(pgrd)
Definition grid.F90:2271
grid::grid_getgridrownum
integer(c_int) function, public grid_getgridrownum(pgrd, ry)
Definition grid.F90:2040
grid::lu_grid
integer(c_int) lu_grid
Definition grid.F90:174
grid::grid_readarcgrid_fn
type(general_grid_t) function, pointer grid_readarcgrid_fn(sfilename, idatatype)
Definition grid.F90:514
grid::grid_createsimple
type(general_grid_t) function, pointer, public grid_createsimple(inx, iny, rx0, ry0, rgridcellsize, idatatype)
Definition grid.F90:281
grid::grid_datatype_int
integer(c_int), parameter, public grid_datatype_int
Definition grid.F90:25
grid::grid_destroy
subroutine, public grid_destroy(pgrd)
Definition grid.F90:366
grid::grid_writearcgrid
subroutine, public grid_writearcgrid(sfilename, pgrd)
Definition grid.F90:1056
grid::row
integer(c_int), parameter, private row
Definition grid.F90:171
grid::grid_lookuprow
subroutine grid_lookuprow(pgrd, ryval, ibefore, iafter, rfrac)
Definition grid.F90:1495
grid::lu_temp
integer(c_int) lu_temp
Definition grid.F90:174
grid::grid_nodata_double
real(c_double), parameter, public grid_nodata_double
Definition grid.F90:22
grid::grid_set_nodata_value
subroutine, public grid_set_nodata_value(pgrd, ivalue, fvalue)
Definition grid.F90:2167
grid::grid_lookupcolumn
subroutine grid_lookupcolumn(pgrd, rxval, ibefore, iafter, rfrac)
Definition grid.F90:1427
grid::grid_gridtopoint_int
integer(c_int) function grid_gridtopoint_int(pgrdfrom, pgrdto, icol, irow)
Definition grid.F90:2304
grid::grid_gridtogrid_int
subroutine, public grid_gridtogrid_int(pgrdfrom, pgrdto, lusemajorityfilter)
Definition grid.F90:2350
grid::grid_getgridcolrownum
integer(c_int) function, dimension(2) grid_getgridcolrownum(pgrd, rx, ry)
Definition grid.F90:2061
grid::grid_rowcolfallsinsidegrid
logical(c_bool) function, public grid_rowcolfallsinsidegrid(pgrd, irow, icol)
Definition grid.F90:1993
grid::grid_nodata_int
integer(c_int), parameter, public grid_nodata_int
Definition grid.F90:20
grid::grid_datatype_real
integer(c_int), parameter, public grid_datatype_real
Definition grid.F90:26
grid::grid_readexisting
subroutine, public grid_readexisting(sfilename, sfiletype, pgrd)
Definition grid.F90:470
grid::output_arc
integer(c_int), parameter, public output_arc
Definition grid.F90:30
grid::grid_interpolate
real(c_float) function, public grid_interpolate(pgrd, rxval, ryval)
Definition grid.F90:1739
grid::column
integer(c_int), parameter, private column
Definition grid.F90:170
grid::output_surfer
integer(c_int), parameter, public output_surfer
Definition grid.F90:29
grid::grid_readsurfergrid_sub
subroutine grid_readsurfergrid_sub(sfilename, pgrd)
Definition grid.F90:957
grid::nc_fill_int
integer(c_int), parameter nc_fill_int
Definition grid.F90:33
grid::nc_fill_double
real(c_double), parameter nc_fill_double
Definition grid.F90:35
grid::grid_gridtogrid_sgl
subroutine, public grid_gridtogrid_sgl(pgrdfrom, pgrdto)
Definition grid.F90:2429
grid::grid_getgridy
real(c_double) function, public grid_getgridy(pgrd, irownumber)
Definition grid.F90:2189
grid::grid_writesurfergrid
subroutine, public grid_writesurfergrid(sfilename, pgrd)
Definition grid.F90:1146
grid::grid_checkintegergridvalues
subroutine grid_checkintegergridvalues(pgrd, sfilename)
Definition grid.F90:1621
grid::grid_conform
logical(c_bool) function, public grid_conform(pgrd1, pgrd2, rtolerance)
Definition grid.F90:1266
grid::grid_checkforproj4error
subroutine, public grid_checkforproj4error(iretval, sfromproj4, stoproj4)
Definition grid.F90:1676
grid::grid_set_output_directory_name
subroutine, public grid_set_output_directory_name(sdirname)
Definition grid.F90:182
grid::grid_transform
subroutine, public grid_transform(pgrd, sfromproj4, stoproj4, rx, ry)
Call PROJ4 to transform coordinates.
Definition grid.F90:1540
grid::grid_getgridx
real(c_double) function, public grid_getgridx(pgrd, icolumnnumber)
Definition grid.F90:2177
grid::grid_readsurfergrid_fn
type(general_grid_t) function, pointer grid_readsurfergrid_fn(sfilename, idatatype)
Definition grid.F90:875
grid::grid_datatype_double
integer(c_int), parameter, public grid_datatype_double
Definition grid.F90:27
grid::grid_readarcgrid_sub
subroutine grid_readarcgrid_sub(sfilename, pgrd)
Definition grid.F90:685
grid::grid_majorityfilter_int
integer(c_int) function grid_majorityfilter_int(pgrdfrom, itargetcol, itargetrow, inodatavalue, ispread)
Definition grid.F90:2486
logfiles::logs
type(logfile_t), public logs
Definition logfiles.F90:62
grid::error_message_t
Definition grid.F90:93
grid::general_grid_t
Definition grid.F90:57
grid::grid_bounds_t
Definition grid.F90:84