forked from NOAA-EMC/NCEPLIBS-ip2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.iplib2
599 lines (481 loc) · 26.8 KB
/
README.iplib2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
Documentation of the general interpolation library 2 - ip2lib
--------------------------------------------------------------------------------
I. Introduction
The NCEP general interpolation library 2 (ip2lib) contains Fortran 90 subprograms
to be used for interpolating between nearly all grids used at NCEP.
The library is particularly efficient when interpolating many fields at one time.
The library has been extensively tested with AIX and Intel Fortran compilers.
There are currently six interpolation methods available in the library:
bilinear, bicubic, neighbor, budget, spectral and neighbor-budget.
Some of the methods have interpolation sub-options. A few methods
have restrictions on the type of input or output grids.
Also, several methods can perform interpolation on fields with bitmaps
(i.e. some points on the input grid may be undefined). In this case,
the bitmap is interpolated to the output grid. Only valid input points
are used to interpolate to valid output points. An output bitmap will also be
created to locate invalid data where the output grid extends outside the domain
of the input grid.
The driver routine for interpolating scalars is ipolates, while the routine
for interpolating vectors is ipolatev. The interpolation method is chosen
via the first argument of these routines (variable IP). Sub-options are
set via the IPOPT array.
Bilinear interpolation is chosen by setting IP=0. This method has two
sub-options. (1) The percent of valid input data required to make output
data (the default is 50%). (2) If valid input data is not found near an
a spiral search may be performed. The spiral search is only
an option for scalar data. The bilinear method also has no restrictions
and can interpolate with bitmaps.
Bicubic interpolation is chosen by setting IP=1. This method has two sub-options,
(1) A monotonic constraint option for straight bicubic or
for constraining the output value to be within the range of the four
surrounding input values. (2) The percent of valid input data
required to make output data, which defaults to 50%. Note: the bicubic method
cannot interpolate data with bitmaps.
Neighbor interpolation is chosen by setting IP=2. Neighbor interpolation
means that the output value is set to the nearest input value. It would be
appropriate for interpolating integer fields such as vegetation index.
This method has one sub-option: If valid input data is not found near an
an output point, a spiral search is optionally performed. The neighbor
method has no restrictions and can interpolate with bitmaps.
Budget interpolation is chosen by setting IP=3. Budget interpolation
means a low-order interpolation method that quasi-conserves area averages.
It would be appropriate for interpolating budget fields such as precipitation.
This method assumes that the field really represents box averages where each
box extends halfway to its neighboring grid point in each direction.
The method actually averages bilinearly interpolated values in a square array
of points distributed within each output grid box. This method can
interpolate data with bitmaps. There are several sub-options:
(1) The number of points in the radius of the square array may be set.
The default is 2, meaning that 25 sample points will be averaged
for each output value.
(2) The respective averaging weights for the radius points are adjustable.
The default is for all weights equal to 1, giving an unweighted average.
(3) Optionally, one may assume the boxes stretch nearly all the way to each of
the neighboring grid points and the weights are the adjoint of the bilinear
interpolation weights.
(4) The percent of valid input data required to make output data is
adjustable. The default is 50%.
(5) In cases where there is no or insufficient valid input data,
a spiral search may be invoked to search for the nearest valid data.
search square (scalar interpolation only).
Spectral interpolation is chosen by setting IP=4. This method has two sub-options,
to (1) set the spectral shape (triangular or rhomboidal) and (2) set the
spectral truncation. The input grid must be a global cylindrical grid
(either Gaussian or equidistant). This method cannot interpolate data
with bitmaps. Unless the output grid is a global cylindrical
grid, a polar stereographic grid centered at the pole, or a Mercator grid,
this method can be quite expensive.
Neighbor-budget interpolation is chosen by setting IP=6. This method
computes weighted averages of neighbor points arranged in a square box
centered around each output grid point and stretching nearly halfway
to each of the neighboring grid points. The main difference with
the budget interpolation (IP=3) is neighbor vs bilinear interpolation
of the square box of points. There are the following sub-options:
(1) The number of points in the radius of the square array may be set.
The default is 2, meaning that 25 sample points will be averaged
for each output value.
(2) The respective averaging weights for the radius points are adjustable.
The default is for all weights equal to 1, giving an unweighted average.
(3) The percent of valid input data required to make output data is
adjustable. The default is 50%.
The library can handle two-dimensional vector fields as well as scalar fields.
The input and output vectors are rotated if necessary so that they are
either resolved relative to their defined grid in the direction of
increasing x and y coordinates or resolved relative to eastward and northward
directions on the earth. The rotation is determined by the grid definitions.
Vectors are generally interpolated (by all methods but spectral interpolation)
by moving the relevant input vectors along a great circle to the output point,
keeping their orientations with respect to the great circle constant, before
independently interpolating the respective components. This ensures that vector
interpolation will be consistent over the whole globe including the poles.
The input and output grids are defined by their respective GRIB2 grid definition
template and template number as decoced by the NCEP G2 library. There
are six map projections recognized by the library:
Grid def. template # Map projection
-------------------- ---------------------------------
00 Equidistant cyclindrical
01 Rotated equidistant cylindrical
10 Mercator cyclindrical
20 Polar stereographic azimuthal
30 Lambert conformal conical
40 Gaussian equidistant cyclindrical
If the output grid definition template number is negative, then the
output data may be just a set of station points. In this case, the user must pass
the number of points to be output along with their latitudes and longitudes.
For vector interpolation, the vector rotations parameters must also be passed.
On the other hand, for non-negative output data representation types,
the number of output grid points and their latitudes and longitudes
(and the vector rotation parameters for vector interpolation) are all
returned by the interpolation subprograms.
If an output equidistant cylindrical grid contains multiple pole points, then
the pole points are forced to be self-consistent. That is, scalar fields
are obliged to be constant at the pole and vector components are obliged
to exhibit a wavenumber one variation at the pole.
Generally, only regular grids can be interpolated in this library. However,
the thinned WAFS grids may be expanded to a regular grid (or vice versa)
using subprograms ipxwafs/2/3. Eta data (with Arakawa "E" staggering)
on the "H" or "V" grid may be expanded to a filled regular grid (or vice versa)
using subprogram ipxetas.
The return code issued by an interpolation subprogram determines whether
it ran successfully or how it failed. Check nonzero return codes
against the docblock of the respective subprogram.
Developers are encouraged to create additional interpolation methods or
to create additional map projection "wizards" for ip2lib.
Questions may be directed to: [email protected]
II. Entry point list
Name Function
---- ------------------------------------------------------------------
Scalar field interpolation subprograms
IPOLATES IREDELL'S POLATE FOR SCALAR FIELDS
POLATES0 INTERPOLATE SCALAR FIELDS (BILINEAR)
POLATES1 INTERPOLATE SCALAR FIELDS (BICUBIC)
POLATES2 INTERPOLATE SCALAR FIELDS (NEIGHBOR)
POLATES3 INTERPOLATE SCALAR FIELDS (BUDGET)
POLATES4 INTERPOLATE SCALAR FIELDS (SPECTRAL)
POLATES6 INTERPOLATE SCALAR FIELDS (NEIGHBOR-BUDGET)
POLFIXS MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
Vector field interpolation subprograms
IPOLATEV IREDELL'S POLATE FOR VECTOR FIELDS
POLATEV0 INTERPOLATE VECTOR FIELDS (BILINEAR)
POLATEV1 INTERPOLATE VECTOR FIELDS (BICUBIC)
POLATEV2 INTERPOLATE VECTOR FIELDS (NEIGHBOR)
POLATEV3 INTERPOLATE VECTOR FIELDS (BUDGET)
POLATEV4 INTERPOLATE VECTOR FIELDS (SPECTRAL)
POLATEV6 INTERPOLATE VECTOR FIELDS (NEIGHBOR-BUDGET)
MOVECT MOVE A VECTOR ALONG A GREAT CIRCLE
POLFIXV MAKE MULTIPLE POLE VECTOR VALUES CONSISTENT
Grid description section decoders
GDSWZD GRID DESCRIPTION SECTION WIZARD
GDSWZD_C 'C' WRAPPER FOR CALLING GDSWZD
GDSWZD_EQUID_CYLIND GDS WIZARD FOR EQUIDISTANT CYCLINDRICAL
GDSWZD_MERCATOR GDS WIZARD FOR MERCATOR CYCLINDRICAL
GDSWZD_LAMBERT_CONF GDS WIZARD FOR LAMBERT CONFORMAL CONICAL
GDSWZD_GAUSSIAN GDS WIZARD FOR GAUSSIAN CYCLINDRICAL
GDSWZD_POLAR_STEREO GDS WIZARD FOR POLAR STEREOGRAPHIC
GDSWZD_ROT_EQUID_CYLIND_EGRID GDS WIZARD FOR ROTATED EQUIDISTANT CYCLINDRICAL
"E" STAGGER.
GDSWZD_ROT_EQUID_CYLIND GDS WIZARD FOR ROTATED EQUIDISTANT CYCLINDRICAL
NON "E" STAGGER.
IJKGDS0/1 RETURN FIELD POSITION FOR A GIVEN GRID POINT
Transform subprograms for special irregular grids
IPXWAFS/2/3 EXPAND OR CONTRACT WAFS GRIDS
III. How to invoke ip2lib: examples
***********************************************************************
Example 1. Read a grib 2 file of scalar data on a global regular
1-deg lat/lon grid and call ipolates to interpolate
it to NCEP standard grid 218, a lambert conformal grid.
Uses the NCEP G2 library to degrib the data.
***********************************************************************
program example_1
use grib_mod ! ncep grib 2 library
implicit none
character(len=100) :: input_file
integer :: iunit, iret, lugi
integer :: mi, mo, no
integer, allocatable :: ibi(:), ibo(:)
integer :: ip, ipopt(20)
integer :: j, jdisc, jpdtn, jgdtn, k, km
integer :: jids(200), jgdt(200), jpdt(200)
integer :: idim_input, jdim_input
integer :: idim_output, jdim_output
logical :: unpack
logical*1, allocatable :: input_bitmap(:,:), output_bitmap(:,:)
real, allocatable :: input_data(:,:)
real, allocatable :: output_rlat(:), output_rlon(:)
real, allocatable :: output_data(:,:)
type(gribfield) :: gfld_input
!---------------------------------------------------------------------------
! the output grid specs. this is ncep grid 218, a lambert conformal
! grid. the grid definition information is stored in section 3
! of a grib 2 message.
!---------------------------------------------------------------------------
integer, parameter :: igdtnum218 = 30 ! grid definition template number.
! "30" is lambert conformal.
integer, parameter :: igdtlen218 = 22 ! number of array elements needed
! for a lambert conf. grid definition
! template.
integer :: igdtmpl218(igdtlen218) ! the grid definition template.
! the entries are:
! 1 -shape of earth, oct 15
! 2 -scale factor, spherical earth, oct 16
! 3 -scaled value, spherical earth, octs 17-20
! 4 -scale factor, major axis of
! elliptical earth, oct 21
! 5 -scaled value of major axis of
! elliptical earth, octs 22-25
! 6 -scale factor, minor axis of
! elliptical earth, oct 26
! 7 -scaled value of minor axis of
! elliptical earth, octs 27-30
! 8 -number points along x-axis, octs 31-34
! 9 -number points along y-axis, octs 35-38
! 10-latitude of first point, octs 39-42
! 11-longitude of first point, octs 43-46
! 12-resolution and component flags, oct 47
! 13-latitude where grid lengths specified,
! octs 48-51
! 14-longitude parallel to y-axis, octs 52-55
! 15-x-direction grid length, octs 56-59
! 16-y-direction grid length, octs 60-63
! 17-projection center flag, oct 64
! 18-scanning mode, oct 65
! 19-first tangent latitude from pole, octs 66-69
! 20-second tangent latitude from pole, octs 70-73
! 21-latitude of south pole, octs 74-77
! 22-longitude of south pole, octs 78-81
integer, parameter :: missing=b'11111111111111111111111111111111'
data igdtmpl218 / 6, 255, missing, 255, missing, 255, missing, 614, 428, &
12190000, 226541000, 56, 25000000, 265000000, &
12191000, 12191000, 0, 64, 25000000, 25000000, -90000000, 0/
!---------------------------------------------------------------------------
! open the grib 2 file containing data to be interpolated. for this
! example, there are two data records.
!---------------------------------------------------------------------------
iunit=9
input_file="${path}/input.data.grib2"
call baopenr (iunit, input_file, iret)
!---------------------------------------------------------------------------
! prep for call to g2 library to degrib data. the data are on a regular
! lat/lon grid with i/j dimension of 360/181.
!---------------------------------------------------------------------------
idim_input = 360 ! the i/j dimensions of input grid
jdim_input = 181
mi = idim_input * jdim_input ! total number of pts, input grid
jdisc = -1 ! search for any discipline
jpdtn = -1 ! search for any product definition template number
jgdtn = 0 ! search for grid definition template number 0 - regular lat/lon grid
jids = -9999 ! array of values in identification section, set to wildcard
jgdt = -9999 ! array of values in grid definition template 3.m
jgdt(8) = idim_input ! search for grid with i/j of 360/181
jgdt(9) = jdim_input
jpdt = -9999 ! array of values in product definition template 4.n
unpack = .true. ! unpack data
lugi = 0 ! no index file
nullify(gfld_input%idsect)
nullify(gfld_input%local)
nullify(gfld_input%list_opt)
nullify(gfld_input%igdtmpl) ! holds the grid definition template information
nullify(gfld_input%ipdtmpl)
nullify(gfld_input%coord_list)
nullify(gfld_input%idrtmpl)
nullify(gfld_input%bmap) ! holds the bitmap
nullify(gfld_input%fld) ! holds the data
!---------------------------------------------------------------------------
! degrib the data. non-zero "iret" indicates a problem during degrib.
!---------------------------------------------------------------------------
km = 2 ! number of records to interpolate
allocate(ibi(km))
allocate(input_bitmap(mi,km))
allocate(input_data(mi,km))
do j = 0, (km-1) ! number of records to skip
call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
unpack, k, gfld_input, iret)
if (iret /= 0) stop
!---------------------------------------------------------------------------
! does input data have a bitmap?
!---------------------------------------------------------------------------
if (gfld_input%ibmap==0) then ! input data has bitmap
ibi(k) = 1 ! tell ipolates to use bitmap
input_bitmap(:,k) = gfld_input%bmap
else ! no bitmap, data everywhere
ibi(k) = 0 ! tell ipolates there is no bitmap
input_bitmap(:,k) = .true.
endif
input_data(:,k) = gfld_input%fld ! the input data field
enddo
call baclose (iunit, iret)
!---------------------------------------------------------------------------
! setup arguments for ipolates (scalar interpolation) call.
!---------------------------------------------------------------------------
ip = 0 ! bilinear interpolation
ipopt = 0 ! options for bilinear:
ipopt(1) = 75 ! set minimum mask to 75%
!---------------------------------------------------------------------------
! the i/j dimensions of the output grid.
!---------------------------------------------------------------------------
idim_output = igdtmpl218(8)
jdim_output = igdtmpl218(9)
mo = idim_output * jdim_output ! total number of output pts
!---------------------------------------------------------------------------
! will hold the latitude, longitude, data and bitmap on the output grid,
! which are computed in ipolates.
!---------------------------------------------------------------------------
allocate (ibo(km)) ! bitmap flags on output grid
allocate (output_rlat(mo))
allocate (output_rlon(mo))
allocate (output_data(mo,km))
allocate (output_bitmap(mo,km))
!---------------------------------------------------------------------------
! call ipolates to interpolate scalar data. non-zero "iret" indicates
! a problem.
!---------------------------------------------------------------------------
call ipolates(ip, ipopt, gfld_input%igdtnum, gfld_input%igdtmpl, &
gfld_input%igdtlen, igdtnum218, igdtmpl218, igdtlen218, &
mi, mo, km, ibi, input_bitmap, input_data, no, output_rlat, &
output_rlon, ibo, output_bitmap, output_data, iret)
if (iret /= 0) stop
!---------------------------------------------------------------------------
! write interpolated data to file. if ipolates computed a bitmap (ibo==1)
! for the output grid, one may mask out this data with a flag value.
!---------------------------------------------------------------------------
open (10, file="./output.bin", access='direct', recl=idim_output*jdim_output*4)
do k = 1, km
if(ibo(k)==1) where (.not. output_bitmap(:,k)) output_data(:,k) = -999.
write(10, rec=k) output_data(:,k)
enddo
write(10, rec=km+1) output_rlat
write(10, rec=km+2) output_rlon
close(10)
end program example_1
***********************************************************************
Example 2. Read a grib 2 file of u/v wind data on a global regular
1-deg lat/lon grid and call ipolatev to interpolate
it to four random station points. Uses the NCEP
G2 library to degrib the data.
***********************************************************************
program example_2
use grib_mod ! ncep grib 2 library
implicit none
character(len=100) :: input_file
integer :: iunit, iret, lugi
integer :: mi, mo, no
integer :: ibi, ibo
integer :: ip, ipopt(20)
integer :: j, jdisc, jpdtn, jgdtn, k, km
integer :: jids(200), jgdt(200), jpdt(200)
integer :: idim_input, jdim_input
logical :: unpack
logical*1, allocatable :: input_bitmap(:), output_bitmap(:)
real, allocatable :: input_u_data(:), input_v_data(:)
real, allocatable :: output_rlat(:), output_rlon(:)
real, allocatable :: output_crot(:), output_srot(:)
real, allocatable :: output_u_data(:), output_v_data(:)
type(gribfield) :: gfld_input
!---------------------------------------------------------------------------
! the output "grid" is a series of random station points. in this case,
! set the grid definition template number of a negative number.
! the grid definition template array information is not used, so set
! to a flag value.
!---------------------------------------------------------------------------
integer, parameter :: igdtnumo = -1
integer, parameter :: igdtleno = 1
integer :: igdtmplo(igdtleno)
data igdtmplo / -9999 /
!---------------------------------------------------------------------------
! open the grib 2 file containing data to be interpolated. for this
! example, there is one record of u-wind and v-wind.
!---------------------------------------------------------------------------
iunit=9
input_file="./reg_tests/copygb2/data/uv_wind.grb2"
call baopenr (iunit, input_file, iret)
!---------------------------------------------------------------------------
! prep for call to g2 library to degrib data. the data are on a regular
! lat/lon grid with i/j dimension of 360/181.
!---------------------------------------------------------------------------
idim_input = 360 ! the i/j dimensions of input grid
jdim_input = 181
mi = idim_input * jdim_input ! total number of pts, input grid
jdisc = -1 ! search for any discipline
jpdtn = -1 ! search for any product definition template number
jgdtn = 0 ! search for grid definition template number 0 - regular lat/lon grid
jids = -9999 ! array of values in identification section, set to wildcard
jgdt = -9999 ! array of values in grid definition template 3.m
jgdt(8) = idim_input ! search for grid with i/j of 360/181
jgdt(9) = jdim_input
jpdt = -9999 ! array of values in product definition template 4.n
unpack = .true. ! unpack data
lugi = 0 ! no index file
nullify(gfld_input%idsect)
nullify(gfld_input%local)
nullify(gfld_input%list_opt)
nullify(gfld_input%igdtmpl) ! holds the grid definition template information
nullify(gfld_input%ipdtmpl)
nullify(gfld_input%coord_list)
nullify(gfld_input%idrtmpl)
nullify(gfld_input%bmap) ! holds the bitmap
nullify(gfld_input%fld) ! holds the data
!---------------------------------------------------------------------------
! degrib the data. non-zero "iret" indicates a problem during degrib.
!---------------------------------------------------------------------------
allocate(input_bitmap(mi))
allocate(input_u_data(mi))
allocate(input_v_data(mi))
!---------------------------------------------------------------------------
! read u-wind record.
!---------------------------------------------------------------------------
j = 0
call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
unpack, k, gfld_input, iret)
if (iret /= 0) stop
!---------------------------------------------------------------------------
! does input data have a bitmap?
!---------------------------------------------------------------------------
if (gfld_input%ibmap==0) then ! input data has bitmap
ibi = 1 ! tell ipolates to use bitmap
input_bitmap = gfld_input%bmap
else ! no bitmap, data everywhere
ibi = 0 ! tell ipolates there is no bitmap
input_bitmap = .true.
endif
input_u_data = gfld_input%fld ! the input u-wind data
!---------------------------------------------------------------------------
! read v-wind record.
!---------------------------------------------------------------------------
j = 1
call getgb2(iunit, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
unpack, k, gfld_input, iret)
if (iret /= 0) stop
input_v_data = gfld_input%fld ! the input v-wind data
call baclose (iunit, iret)
!---------------------------------------------------------------------------
! setup arguments for ipolatev (vector interpolation) call.
!---------------------------------------------------------------------------
km = 1 ! number of records to interpolate
ip = 0 ! bilinear interpolation
ipopt = 0 ! options for bilinear:
ipopt(1) = 75 ! set minimum mask to 75%
!---------------------------------------------------------------------------
! interpolate to four random station points.
!---------------------------------------------------------------------------
mo = 4
no = mo
!---------------------------------------------------------------------------
! when interpolating to random station points, need to pass to ipolatev
! their latitude, longitude and the sines and cosines of the vector
! rotation angles. the vector rotation is defined:
!
! ugrid=crot*uearth-sort*vearth
! vgrid=srot*uearth+cort*vearth
!---------------------------------------------------------------------------
allocate (output_rlat(mo))
allocate (output_rlon(mo))
allocate (output_srot(mo))
allocate (output_crot(mo))
allocate (output_u_data(mo))
allocate (output_v_data(mo))
allocate (output_bitmap(mo))
output_rlat(1) = 45.0
output_rlon(1) = -100.0
output_rlat(2) = 35.0
output_rlon(2) = -100.0
output_rlat(3) = 40.0
output_rlon(3) = -90.0
output_rlat(4) = 35.0
output_rlon(4) = -120.0
output_srot = 0.0 ! no turning of wind
output_crot = 1.0
!---------------------------------------------------------------------------
! call ipolatev to interpolate vector data. non-zero "iret" indicates
! a problem.
!---------------------------------------------------------------------------
call ipolatev(ip, ipopt, gfld_input%igdtnum, gfld_input%igdtmpl, &
gfld_input%igdtlen, igdtnumo, igdtmplo, igdtleno, &
mi, mo, km, ibi, input_bitmap, input_u_data, input_v_data, &
no, output_rlat, output_rlon, output_crot, output_srot, &
ibo, output_bitmap, output_u_data, output_v_data, iret)
if (iret /= 0) stop
do k = 1, mo
print*,'station point ',k,' latitude ',output_rlat(k),' longitude ', &
output_rlon(k), ' u-wind ', output_u_data(k), ' v-wind ', output_v_data(k)
enddo
end program example_2