forked from erget/wgrib2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathiplib.doc
executable file
·2495 lines (2386 loc) · 110 KB
/
iplib.doc
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
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
Documentation of the general interpolation library iplib May 2, 1996
--------------------------------------------------------------------------------
I. Introduction
The general interpolation library iplib contains FORTRAN subprograms
to be used for interpolating between almost any grids used at NCEP.
The library has been optimized for the CRAY machines, taking full advantage
of both the vector and parallel capabilities. The library is particularly
efficient when interpolating many fields at one time. The library should be
transportable to other platforms with compilers allowing automatic arrays.
There are currently five interpolation methods available in the library.
They are respectively bilinear, bicubic, neighbor, budget and spectral
interpolation methods. Some of the methods have interpolation sub-options.
A few methods have restrictions on the type of input or output grids as well.
Also, several methods can perform interpolation on fields with bitmaps
(i.e. some points on the input grid may not have valid data). In this case,
the bitmap is interpolated to the output grid as well. 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. Thus an output bitmap must always be dimensioned.
Bilinear interpolation is chosen by setting IP=0. This method has no options.
This method also has no restrictions and can interpolate with bitmaps.
Bicubic interpolation is chosen by setting IP=1. This method has one option,
a monotonic constraint option IPOPT(1) which is 0 for straight bicubic or
is 1 for contraining the output value to be within the range of the four
surrounding input values. This method cannot interpolate 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 no options. This method also 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. The first option of this
method IPOPT(1) is the number of points in the radius of the square array.
If IPOPT(1)=-1, then the number defaults to 2 meaning that 25 sample points
will be averaged for each output value. Note that if IPOPT(1)=0, then
this method degenerates to simple bilinear interpolation. Further options
(IPOPT(2:2+IPOPT(1)) are the respective averaging weights for the radius points
starting at the center point. If IPOPT(2)=-1, then all the weights default
to 1 giving an unweighted average. This method can interpolate with bitmaps.
The output grid must be a well-defined grid and not a set of stations.
Spectral interpolation is chosen by setting IP=4. This method has two options,
the spectral shape IPOPT(1) which is 0 for triangular or 1 for rhomboidal
and the spectral truncation IPOPT(2). The input grid must be a global
cylindrical grid (either Gaussian or equidistant). This method cannot
interpolate 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.
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 GRIB grid description
sections (GDS) passed in integer form KGDS as decoded by subprogram w3fi63
in w3lib or by subprogram makgds in this library. That is, the interpolation
subprograms can readily interpolate from a GRIB field that is unpacked
by subprogram w3fi63; the interpolation subprograms can also readily
interpolate to an NCEP pre-defined grid that is expanded into KGDS form by
subprogram makgds (which in turn calls w3fi71). There are currently six grid
projections recognized in the library. The projections are respectively
equidistant cylindrical (KGDS(1)=000), Mercator cylindrical (KGDS(1)=001),
Lambert conformal conical (KGDS(1)=003), Gaussian cylindrical (KGDS(1)=004),
polar stereographic azimuthal (KGDS(1)=005), and filled rotated equidistant
cylindrical (KGDS(1)=202).
If the output data representation type is negative (KGDSO(1)<0), then the
output data may be just a set of station points. (This feature cannot be
used with the budget interpolation method.) 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 and the staggered eta grids can be interpolated by using
transform subprograms (ipxwafs and ipxetas, respectively) in this library that
will either expand the irregular grid to a regular grid or contract a regular
grid to an irregular grid as necessary.
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 iplib.
This documentation is divided into 4 chapters. Chapter I is this introduction.
Chapter II is a list of all entry points. Chapter III is a set of examples.
Chapter IV is a recapitulation of all the docblocks. The chapters all start
on a line number that is 1 modulo 60 in order to facilitate laser printing.
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)
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)
MOVECT MOVE A VECTOR ALONG A GREAT CIRCLE
POLFIXV MAKE MULTIPLE POLE VECTOR VALUES CONSISTENT
Grid description section decoders
GDSWIZ GRID DESCRIPTION SECTION WIZARD
GDSWIZ00 GDS WIZARD FOR EQUIDISTANT CYLINDRICAL
GDSWIZ01 GDS WIZARD FOR MERCATOR CYLINDRICAL
GDSWIZ03 GDS WIZARD FOR LAMBERT CONFORMAL CONICAL
GDSWIZ04 GDS WIZARD FOR GAUSSIAN CYLINDRICAL
GAUSSLAT COMPUTE GAUSSIAN LATITUDES
GDSWIZ05 GDS WIZARD FOR POLAR STEREOGRAPHIC AZIMUTHAL
GDSWIZCA GDS WIZARD FOR ROTATED EQUIDISTANT CYLINDRICAL
IJKGDS RETURN FIELD POSITION FOR A GIVEN GRID POINT
MAKGDS MAKE OR BREAK A GRID DESCRIPTION SECTION
Transform subprograms for special irregular grids
IPXWAFS EXPAND OR CONTRACT WAFS GRIDS
IPXETAS EXPAND OR CONTRACT ETA GRIDS
III. Examples
Example 1. Interpolate from an arbitrary input grid (probably 1x1)
to NCEP grid 27 (65x65 northern polar stereographic).
Interpolate heights bilinearly and winds bicubically.
Interpolate soil moisture and precipitation using bitmaps
with the budget method. Encode the soil moisture bitmap in GRIB.
Subprograms GETGB and PUTGB from w3lib are referenced.
c example of using ipolate package.
c see documentation of ipolates and ipolatev
c for further possible options.
integer ipopt(20)
integer jpds(25),jgds(22),kpdsi(25),kgdsi(22),kpdso(25),kgdso(22)
parameter(ji=360*181,ig=27,jo=65*65,km=4)
real rlat(jo),rlon(jo),crot(jo),srot(jo)
integer ibi(km),ibo(km)
logical li(ji,km),lo(jo,km)
real hi(ji,km),ri(ji),ui(ji),vi(ji)
real ho(jo,km),ro(jo),uo(jo),vo(jo)
character gdso(42)
integer lev(km)
data lev/1000,500,250,100/
c define 65x65 grid
call makgds(ig,kgdso,gdso,lengds,iret)
if(iret.ne.0) call exit(iret)
kgdso(4)=-20826! fix w3fi71 error
print *,'kgdso=',kgdso
ipopt=0
c interpolate 4 levels of height to 65x65 bilinearly
do k=1,km
jpds=-1
jpds(5)=7
jpds(6)=100
jpds(7)=lev(k)
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li(1,k),hi(1,k),iret)
if(iret.ne.0) call exit(iret)
call putgb(50,ki,kpdsi,kgdsi,li(1,k),hi(1,k),iret)
if(iret.ne.0) call exit(iret)
ibi(k)=mod(kpdsi(4)/64,2)
print *,'ibi(k)=',ibi(k)
enddo
call ipolates(0,ipopt,kgdsi,kgdso,ji,jo,km,ibi,li,hi,
& ko,rlat,rlon,ibo,lo,ho,iret)
if(iret.ne.0) call exit(iret)
kpdso=kpdsi
kpdso(3)=ig
do k=1,km
kpdso(7)=lev(k)
call putgb(51,ko,kpdso,kgdso,lo(1,k),ho(1,k),iret)
if(iret.ne.0) call exit(iret)
enddo
c interpolate precipitation to 65x65 with budget method
c (zero rain is masked out)
jpds=-1
jpds(5)=61
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li,ri,iret)
if(iret.ne.0) call exit(iret)
call putgb(50,ki,kpdsi,kgdsi,li,ri,iret)
if(iret.ne.0) call exit(iret)
ipopt(1)=-1
ipopt(2)=-1
li(1:ki,1)=ri(1:ki).gt.0.
call ipolates(3,ipopt,kgdsi,kgdso,ji,jo,1,1,li,ri,
& ko,rlat,rlon,ibo,lo,ro,iret)
if(iret.ne.0) call exit(iret)
kpdso=kpdsi
kpdso(3)=ig
call putgb(51,ko,kpdso,kgdso,lo,ro,iret)
if(iret.ne.0) call exit(iret)
c interpolate soil moisture to 65x65 with budget method
jpds=-1
jpds(5)=144
jpds(6)=112
jpds(7)=10
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li,ri,iret)
if(iret.ne.0) call exit(iret)
call putgb(50,ki,kpdsi,kgdsi,li,ri,iret)
if(iret.ne.0) call exit(iret)
ibi(1)=mod(kpdsi(4)/64,2)
ipopt(1)=-1
ipopt(2)=-1
call ipolates(3,ipopt,kgdsi,kgdso,ji,jo,1,ibi,li,ri,
& ko,rlat,rlon,ibo,lo,ro,iret)
if(iret.ne.0) call exit(iret)
kpdso=kpdsi
kpdso(3)=ig
kpdso(4)=128+64*ibo(1)
call putgb(51,ko,kpdso,kgdso,lo,ro,iret)
if(iret.ne.0) call exit(iret)
c interpolate 200 mb winds to 65x65 bicubically
jpds=-1
jpds(5)=33
jpds(6)=100
jpds(7)=200
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li,ui,iret)
if(iret.ne.0) call exit(iret)
call putgb(50,ki,kpdsi,kgdsi,li,ui,iret)
if(iret.ne.0) call exit(iret)
jpds(5)=34
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li,vi,iret)
if(iret.ne.0) call exit(iret)
call putgb(50,ki,kpdsi,kgdsi,li,vi,iret)
if(iret.ne.0) call exit(iret)
ipopt(1)=0
call ipolatev(1,ipopt,kgdsi,kgdso,ji,jo,1,0,li,ui,vi,
& ko,rlat,rlon,crot,srot,ibo,lo,uo,vo,iret)
if(iret.ne.0) call exit(iret)
kpdso=kpdsi
kpdso(3)=ig
kpdso(5)=33
call putgb(51,ko,kpdso,kgdso,lo,uo,iret)
if(iret.ne.0) call exit(iret)
kpdso(5)=34
call putgb(51,ko,kpdso,kgdso,lo,vo,iret)
if(iret.ne.0) call exit(iret)
stop
end
Example 2. Interpolate winds from an arbitrary input grid (probably 1x1)
to 4 station points while truncating spectrally to R30.
c read and unpack the 500 mb winds, truncate to R30,
c and interpolate to 4 corners of a box
integer ipopt(20)
integer jpds(25),jgds(22),kpdsi(25),kgdsi(22),kgdso(22)
parameter(jf=360*181,kp=4)
real rlat(kp),rlon(kp),crot(kp),srot(kp)
logical lgi(jf),lgo(kp)
real ui(jf),vi(jf),uo(kp),vo(kp)
jpds=-1
jpds(5)=33
jpds(6)=100
jpds(7)=500
call getgb(11,31,jf,0,jpds,jgds,kf,kr,kpdsi,kgdsi,
& lgi,ui,iret)
jpds(5)=34
call getgb(11,31,jf,0,jpds,jgds,kf,kr,kpdsi,kgdsi,
& lgi,vi,iret)
kgdso=-1
rlat(1)=20.
rlat(2)=20.
rlat(3)=10.
rlat(4)=10.
rlon(1)=-50.
rlon(2)=-40.
rlon(3)=-50.
rlon(4)=-40.
crot=1.
srot=0.
ipopt(1)=1
ipopt(2)=30
uo=-999
vo=-999
call ipolatev(4,ipopt,kgdsi,kgdso,jf,kp,1,0,lgi,ui,vi,
& kp,rlat,rlon,crot,srot,ibo,lgo,uo,vo,iret)
print '(2(2x,2f8.2))',(uo(k),vo(k),k=1,kp)
end
Example 3. Interpolate winds from an arbitrary input grid (probably 1x1)
bilinearly to 3 station points.
c read and unpack 4 levels of heights and winds
c and interpolate to 3 sonde sites.
integer ipopt(20)
integer jpds(25),jgds(22),kpdsi(25),kgdsi(22),kgdso(22)
parameter(ji=360*181,km=4,jo=3)
real rlat(jo),rlon(jo),crot(jo),srot(jo)
integer ibi(km),ibo(km)
logical li(ji,km),lo(jo,km)
real hi(ji,km),ui(ji,km),vi(ji,km),ho(jo,km),uo(jo,km),vo(jo,km)
integer lev(km)
data lev/1000,500,250,100/
c define output locations
kgdso=-1
rlat(1)=22.2
rlat(2)=33.3
rlat(3)=44.4
rlon(1)=-50.
rlon(2)=-40.
rlon(3)=-30.
crot=1.
srot=0.
c heights
do k=1,km
jpds=-1
jpds(5)=7
jpds(6)=100
jpds(7)=lev(k)
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li(1,k),hi(1,k),iret)
if(iret.ne.0) call exit(iret)
ibi(k)=mod(kpdsi(4)/64,2)
enddo
call ipolates(0,ipopt,kgdsi,kgdso,ji,jo,km,ibi,li,hi,
& jo,rlat,rlon,ibo,lo,ho,iret)
if(iret.ne.0) call exit(iret)
c winds
do k=1,km
jpds=-1
jpds(5)=33
jpds(6)=100
jpds(7)=lev(k)
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li(1,k),ui(1,k),iret)
if(iret.ne.0) call exit(iret)
jpds=-1
jpds(5)=34
jpds(6)=100
jpds(7)=lev(k)
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li(1,k),vi(1,k),iret)
if(iret.ne.0) call exit(iret)
ibi(k)=mod(kpdsi(4)/64,2)
enddo
call ipolatev(0,ipopt,kgdsi,kgdso,ji,jo,km,ibi,li,ui,vi,
& jo,rlat,rlon,crot,srot,ibo,lo,uo,vo,iret)
if(iret.ne.0) call exit(iret)
print '((i8,3(2x,3f8.2)))',
& (lev(k),(ho(j,k),uo(j,k),vo(j,k),j=1,jo),k=1,km)
end
Example 4. Interpolate 850 mb height and winds from the staggered meso-eta
to a regional 0.25 degree grid.
integer ipopt(20)
integer jpds(25),jgds(22),kpdsi(25),kgdsi(22),kpdso(25),kgdso(22)
integer kgdsi2(22)
parameter(ji=361*271,ig=255,jo=121*81)
real rlat(jo),rlon(jo),crot(jo),srot(jo)
logical li(ji),lo(jo)
real fi(ji),fi2(ji),fo(jo)
real vi(ji),vi2(ji),vo(jo)
character gdso(400)
kgdso=0
kgdso(1)=0
kgdso(2)=121
kgdso(3)=81
kgdso(4)=30000
kgdso(5)=-90000
kgdso(6)=128
kgdso(7)=50000
kgdso(8)=-60000
kgdso(9)=250
kgdso(10)=250
kgdso(11)=64
kgdso(19)=0
kgdso(20)=255
if(iret.ne.0) call exit(iret)
ipopt=0
jpds=-1
jpds(6)=100
jpds(7)=850
jpds(5)=7
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li,fi,iret)
if(iret.ne.0) call exit(iret)
call ipxetas(1,ji,ji,1,kgdsi,fi,kgdsi2,fi2,iret)
if(iret.ne.0) call exit(iret)
call ipolates(0,ipopt,kgdsi2,kgdso,ji,jo,1,0,li,fi2,
& ko,rlat,rlon,ibo,lo,fo,iret)
if(iret.ne.0) call exit(iret)
kpdso=kpdsi
kpdso(3)=ig
kpdso(4)=128+64*ibo
kpdso(22)=1
call putgb(51,ko,kpdso,kgdso,lo,fo,iret)
if(iret.ne.0) call exit(iret)
jpds(5)=33
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li,fi,iret)
if(iret.ne.0) call exit(iret)
call ipxetas(2,ji,ji,1,kgdsi,fi,kgdsi2,fi2,iret)
if(iret.ne.0) call exit(iret)
jpds(5)=34
call getgb(11,31,ji,0,jpds,jgds,ki,kr,kpdsi,kgdsi,
& li,vi,iret)
if(iret.ne.0) call exit(iret)
call ipxetas(2,ji,ji,1,kgdsi,vi,kgdsi2,vi2,iret)
if(iret.ne.0) call exit(iret)
call ipolatev(0,ipopt,kgdsi2,kgdso,ji,jo,1,0,li,fi2,vi2,
& ko,rlat,rlon,crot,srot,ibo,lo,fo,vo,iret)
if(iret.ne.0) call exit(iret)
kpdso=kpdsi
kpdso(3)=ig
kpdso(4)=128+64*ibo
kpdso(22)=1
kpdso(5)=33
call putgb(51,ko,kpdso,kgdso,lo,fo,iret)
if(iret.ne.0) call exit(iret)
kpdso(5)=34
call putgb(51,ko,kpdso,kgdso,lo,vo,iret)
if(iret.ne.0) call exit(iret)
end
IV. Docblocks
The primary documentation of iplib is via the docblocks in its subprograms.
The following recapitulation of docblocks is current as of May, 1996.
Docblock for ipolates.
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: IPOLATES IREDELL'S POLATE FOR SCALAR FIELDS
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM INTERPOLATES SCALAR FIELDS
C FROM ANY GRID TO ANY GRID (JOE IRWIN'S DREAM).
C ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
C THE FOLLOWING INTERPOLATION METHODS ARE POSSIBLE:
C (IP=0) BILINEAR
C (IP=1) BICUBIC
C (IP=2) NEIGHBOR
C (IP=3) BUDGET
C (IP=4) SPECTRAL
C SOME OF THESE METHODS HAVE INTERPOLATION OPTIONS AND/OR
C RESTRICTIONS ON THE INPUT OR OUTPUT GRIDS, BOTH OF WHICH
C ARE DOCUMENTED MORE FULLY IN THEIR RESPECTIVE SUBPROGRAMS.
C THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS
C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63).
C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS:
C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
C (KGDS(1)=001) MERCATOR CYLINDRICAL
C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
C (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)
C WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO.
C AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
C AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
C ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
C IF KGDSO(1)<0, IN WHICH CASE THE NUMBER OF POINTS
C AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT.
C INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS.
C OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID
C EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
C THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
C
C PROGRAM HISTORY LOG:
C 96-04-10 IREDELL
C
C USAGE: CALL IPOLATES(IP,IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
C & NO,RLAT,RLON,IBO,LO,GO,IRET)
C
C INPUT ARGUMENT LIST:
C IP - INTEGER INTERPOLATION METHOD
C (IP=0 FOR BILINEAR;
C IP=1 FOR BICUBIC;
C IP=2 FOR NEIGHBOR;
C IP=3 FOR BUDGET;
C IP=4 FOR SPECTRAL)
C IPOPT - INTEGER (20) INTERPOLATION OPTIONS
C (IP=0: (NO OPTIONS)
C IP=1: CONSTRAINT OPTION
C IP=2: (NO OPTIONS)
C IP=3: NUMBER IN RADIUS, RADIUS WEIGHTS ...
C IP=4: SPECTRAL SHAPE, SPECTRAL TRUNCATION)
C KGDSI - INTEGER (22) INPUT GDS PARAMETERS AS DECODED BY W3FI63
C KGDSO - INTEGER (22) OUTPUT GDS PARAMETERS AS DECODED BY W3FI63
C MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
C OR DIMENSION OF INPUT GRID FIELDS IF KM=1
C MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
C OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
C KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE
C IBI - INTEGER (KM) INPUT BITMAP FLAGS
C LI - LOGICAL (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
C GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
C RLAT - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
C RLON - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
C
C OUTPUT ARGUMENT LIST:
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
C RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
C RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
C IBO - INTEGER (KM) OUTPUT BITMAP FLAGS
C LO - LOGICAL (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
C GO - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
C IRET - INTEGER RETURN CODE
C 0 SUCCESSFUL INTERPOLATION
C 1 UNRECOGNIZED INTERPOLATION METHOD
C 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
C 3 UNRECOGNIZED OUTPUT GRID
C 1X INVALID BICUBIC METHOD PARAMETERS
C 3X INVALID BUDGET METHOD PARAMETERS
C 4X INVALID SPECTRAL METHOD PARAMETERS
C
C SUBPROGRAMS CALLED:
C POLATES0 INTERPOLATE SCALAR FIELDS (BILINEAR)
C POLATES1 INTERPOLATE SCALAR FIELDS (BICUBIC)
C POLATES2 INTERPOLATE SCALAR FIELDS (NEIGHBOR)
C POLATES3 INTERPOLATE SCALAR FIELDS (BUDGET)
C POLATES4 INTERPOLATE SCALAR FIELDS (SPECTRAL)
C
C REMARKS: EXAMPLES DEMONSTRATING RELATIVE CPU COSTS.
C THIS EXAMPLE IS INTERPOLATING 12 LEVELS OF TEMPERATURES
C FROM THE 360 X 181 GLOBAL GRID (NCEP GRID 3)
C TO THE 93 X 68 HAWAIIAN MERCATOR GRID (NCEP GRID 204).
C THE EXAMPLE TIMES ARE FOR THE C90. AS A REFERENCE, THE CP TIME
C FOR UNPACKING THE GLOBAL 12 TEMPERATURE FIELDS IS 0.04 SECONDS.
C
C METHOD IP IPOPT CP SECONDS
C -------- -- ------------- ----------
C BILINEAR 0 0.03
C BICUBIC 1 0 0.07
C BICUBIC 1 1 0.07
C NEIGHBOR 2 0.01
C BUDGET 3 -1,-1 0.46
C SPECTRAL 4 0,40 0.23
C SPECTRAL 4 1,40 0.25
C SPECTRAL 4 0,-1 0.43
C
C THE SPECTRAL INTERPOLATION IS FAST FOR THE MERCATOR GRID.
C HOWEVER, FOR SOME GRIDS THE SPECTRAL INTERPOLATION IS SLOW.
C THE FOLLOWING EXAMPLE IS INTERPOLATING 12 LEVELS OF TEMPERATURES
C FROM THE 360 X 181 GLOBAL GRID (NCEP GRID 3)
C TO THE 93 X 65 CONUS LAMBERT CONFORMAL GRID (NCEP GRID 211).
C
C METHOD IP IPOPT CP SECONDS
C -------- -- ------------- ----------
C BILINEAR 0 0.02
C BICUBIC 1 0 0.06
C BICUBIC 1 1 0.06
C NEIGHBOR 2 0.01
C BUDGET 3 -1,-1 0.45
C SPECTRAL 4 0,40 4.10
C SPECTRAL 4 1,40 5.23
C SPECTRAL 4 0,-1 11.75
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 77
C
C$$$
Docblock for polates0.
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: POLATES0 INTERPOLATE SCALAR FIELDS (BILINEAR)
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS BILINEAR INTERPOLATION
C FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
C NO OPTIONS ARE ALLOWED.
C ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
C THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS
C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63).
C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS:
C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
C (KGDS(1)=001) MERCATOR CYLINDRICAL
C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
C (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)
C WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO.
C AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
C AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
C ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
C IF KGDSO(1)<0, IN WHICH CASE THE NUMBER OF POINTS
C AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT.
C INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS.
C OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID
C EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
C THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
C
C PROGRAM HISTORY LOG:
C 96-04-10 IREDELL
C
C USAGE: CALL POLATES0(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
C & NO,RLAT,RLON,IBO,LO,GO,IRET)
C
C INPUT ARGUMENT LIST:
C IPOPT - INTEGER (20) INTERPOLATION OPTIONS (NO OPTIONS)
C KGDSI - INTEGER (22) INPUT GDS PARAMETERS AS DECODED BY W3FI63
C KGDSO - INTEGER (22) OUTPUT GDS PARAMETERS AS DECODED BY W3FI63
C (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
C MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
C OR DIMENSION OF INPUT GRID FIELDS IF KM=1
C MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
C OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
C KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE
C IBI - INTEGER (KM) INPUT BITMAP FLAGS
C LI - LOGICAL (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
C GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
C RLAT - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
C RLON - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
C
C OUTPUT ARGUMENT LIST:
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
C RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
C RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
C IBO - INTEGER (KM) OUTPUT BITMAP FLAGS
C LO - LOGICAL (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
C GO - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
C IRET - INTEGER RETURN CODE
C 0 SUCCESSFUL INTERPOLATION
C 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
C 3 UNRECOGNIZED OUTPUT GRID
C
C SUBPROGRAMS CALLED:
C GDSWIZ GRID DESCRIPTION SECTION WIZARD
C (IJKGDS) RETURN FIELD POSITION FOR A GIVEN GRID POINT
C POLFIXS MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 77
C
C$$$
Docblock for polates1.
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: POLATES1 INTERPOLATE SCALAR FIELDS (BICUBIC)
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS BICUBIC INTERPOLATION
C FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
C IT REQUIRES THAT NO INPUT FIELDS HAVE BITMAPS (IBI=0).
C OPTIONS ALLOW CHOICES BETWEEN STRAIGHT BICUBIC (IPOPT(1)=0)
C AND CONSTRAINED BICUBIC (IPOPT(1)=1) WHERE THE VALUE IS
C CONFINED WITHIN THE RANGE OF THE SURROUNDING 4 POINTS.
C BILINEAR USED WITHIN ONE GRID LENGTH OF BOUNDARIES.
C ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
C THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS
C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63).
C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS:
C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
C (KGDS(1)=001) MERCATOR CYLINDRICAL
C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
C (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)
C WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO.
C AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
C AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
C ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
C IF KGDSO(1)<0, IN WHICH CASE THE NUMBER OF POINTS
C AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT.
C OUTPUT BITMAPS WILL ONLY BE CREATED WHEN THE OUTPUT GRID
C EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
C THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
C
C PROGRAM HISTORY LOG:
C 96-04-10 IREDELL
C
C USAGE: CALL POLATES1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
C & NO,RLAT,RLON,IBO,LO,GO,IRET)
C
C INPUT ARGUMENT LIST:
C IPOPT - INTEGER (20) INTERPOLATION OPTIONS
C IPOPT(1)=0 FOR STRAIGHT BICUBIC;
C IPOPT(1)=1 FOR CONSTRAINED BICUBIC WHERE VALUE IS
C CONFINED WITHIN THE RANGE OF THE SURROUNDING 4 POINTS.
C KGDSI - INTEGER (22) INPUT GDS PARAMETERS AS DECODED BY W3FI63
C KGDSO - INTEGER (22) OUTPUT GDS PARAMETERS AS DECODED BY W3FI63
C (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
C MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
C OR DIMENSION OF INPUT GRID FIELDS IF KM=1
C MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
C OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
C KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE
C IBI - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
C LI - LOGICAL (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
C GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
C RLAT - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
C RLON - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
C
C OUTPUT ARGUMENT LIST:
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
C RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
C RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
C IBO - INTEGER (KM) OUTPUT BITMAP FLAGS
C LO - LOGICAL (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
C GO - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
C IRET - INTEGER RETURN CODE
C 0 SUCCESSFUL INTERPOLATION
C 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
C 3 UNRECOGNIZED OUTPUT GRID
C 11 INVALID INPUT BITMAPS
C
C SUBPROGRAMS CALLED:
C GDSWIZ GRID DESCRIPTION SECTION WIZARD
C (IJKGDS) RETURN FIELD POSITION FOR A GIVEN GRID POINT
C POLFIXS MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 77
C
C$$$
Docblock for polates2.
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: POLATES2 INTERPOLATE SCALAR FIELDS (NEIGHBOR)
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS NEIGHBOR INTERPOLATION
C FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
C NO OPTIONS ARE ALLOWED.
C ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
C THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS
C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63).
C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS:
C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
C (KGDS(1)=001) MERCATOR CYLINDRICAL
C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
C (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)
C WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO.
C AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
C AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
C ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
C IF KGDSO(1)<0, IN WHICH CASE THE NUMBER OF POINTS
C AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT.
C INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS.
C OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID
C EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
C THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
C
C PROGRAM HISTORY LOG:
C 96-04-10 IREDELL
C
C USAGE: CALL POLATES2(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
C & NO,RLAT,RLON,IBO,LO,GO,IRET)
C
C INPUT ARGUMENT LIST:
C IPOPT - INTEGER (20) INTERPOLATION OPTIONS (NO OPTIONS)
C KGDSI - INTEGER (22) INPUT GDS PARAMETERS AS DECODED BY W3FI63
C KGDSO - INTEGER (22) OUTPUT GDS PARAMETERS AS DECODED BY W3FI63
C (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
C MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
C OR DIMENSION OF INPUT GRID FIELDS IF KM=1
C MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
C OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
C KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE
C IBI - INTEGER (KM) INPUT BITMAP FLAGS
C LI - LOGICAL (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
C GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)
C RLAT - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0)
C RLON - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0)
C
C OUTPUT ARGUMENT LIST:
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0)
C RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0)
C RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0)
C IBO - INTEGER (KM) OUTPUT BITMAP FLAGS
C LO - LOGICAL (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
C GO - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
C IRET - INTEGER RETURN CODE
C 0 SUCCESSFUL INTERPOLATION
C 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
C 3 UNRECOGNIZED OUTPUT GRID
C
C SUBPROGRAMS CALLED:
C GDSWIZ GRID DESCRIPTION SECTION WIZARD
C (IJKGDS) RETURN FIELD POSITION FOR A GIVEN GRID POINT
C POLFIXS MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 77
C
C$$$
Docblock for polates3.
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: POLATES3 INTERPOLATE SCALAR FIELDS (BUDGET)
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS BUDGET INTERPOLATION
C FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
C IT REQUIRES A GRID FOR THE OUTPUT FIELDS (KGDSO(1)>=0).
C THE ALGORITHM SIMPLY COMPUTES (WEIGHTED) AVERAGES
C OF BILINEARLY INTERPOLATED POINTS ARRANGED IN A SQUARE BOX
C CENTERED AROUND EACH OUTPUT GRID POINT AND STRETCHING
C NEARLY HALFWAY TO EACH OF THE NEIGHBORING GRID POINTS.
C OPTIONS ALLOW CHOICES OF NUMBER OF POINTS IN EACH RADIUS
C FROM THE CENTER POINT (IPOPT(1)) WHICH DEFAULTS TO 2
C (IF IPOPT(1)=-1) MEANING THAT 25 POINTS WILL BE AVERAGED;
C FURTHER OPTIONS ARE THE RESPECTIVE WEIGHTS FOR THE RADIUS
C POINTS STARTING AT THE CENTER POINT (IPOPT(2:2+IPOPT(1))
C WHICH DEFAULTS TO ALL 1 (IF IPOPT(2)=-1.).
C ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
C THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS
C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63).
C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS:
C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
C (KGDS(1)=001) MERCATOR CYLINDRICAL
C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
C (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)
C WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO.
C AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
C AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
C INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS.
C OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID
C EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
C THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
C
C PROGRAM HISTORY LOG:
C 96-04-10 IREDELL
C
C USAGE: CALL POLATES3(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
C & NO,RLAT,RLON,IBO,LO,GO,IRET)
C
C INPUT ARGUMENT LIST:
C IPOPT - INTEGER (20) INTERPOLATION OPTIONS
C IPOPT(1) IS NUMBER OF RADIUS POINTS
C (DEFAULTS TO 2 IF IPOPT(1)=-1);
C IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS
C (DEFAULTS TO ALL 1 IF IPOPT(2)=-1).
C KGDSI - INTEGER (22) INPUT GDS PARAMETERS AS DECODED BY W3FI63
C KGDSO - INTEGER (22) OUTPUT GDS PARAMETERS AS DECODED BY W3FI63
C MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
C OR DIMENSION OF INPUT GRID FIELDS IF KM=1
C MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
C OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
C KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE
C IBI - INTEGER (KM) INPUT BITMAP FLAGS
C LI - LOGICAL (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
C GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
C
C OUTPUT ARGUMENT LIST:
C NO - INTEGER NUMBER OF OUTPUT POINTS
C RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES
C RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES
C IBO - INTEGER (KM) OUTPUT BITMAP FLAGS
C LO - LOGICAL (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
C GO - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
C IRET - INTEGER RETURN CODE
C 0 SUCCESSFUL INTERPOLATION
C 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
C 3 UNRECOGNIZED OUTPUT GRID
C 31 INVALID UNDEFINED OUTPUT GRID
C 32 INVALID BUDGET METHOD PARAMETERS
C
C SUBPROGRAMS CALLED:
C GDSWIZ GRID DESCRIPTION SECTION WIZARD
C (IJKGDS) RETURN FIELD POSITION FOR A GIVEN GRID POINT
C POLFIXS MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 77
C
C$$$
Docblock for polates4.
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: POLATES4 INTERPOLATE SCALAR FIELDS (SPECTRAL)
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
C
C ABSTRACT: THIS SUBPROGRAM PERFORMS SPECTRAL INTERPOLATION
C FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS.
C IT REQUIRES THAT THE INPUT FIELDS BE UNIFORMLY GLOBAL.
C OPTIONS ALLOW CHOICES BETWEEN TRIANGULAR SHAPE (IPOPT(1)=0)
C AND RHOMBOIDAL SHAPE (IPOPT(1)=1) WHICH HAS NO DEFAULT;
C A SECOND OPTION IS THE TRUNCATION (IPOPT(2)) WHICH DEFAULTS
C TO A SENSIBLE TRUNCATION FOR THE INPUT GRID (IF OPT(2)=-1).
C NOTE THAT IF THE OUTPUT GRID IS NOT FOUND IN A SPECIAL LIST,
C THEN THE TRANSFORM BACK TO GRID IS NOT VERY FAST.
C THIS SPECIAL LIST CONTAINS GLOBAL CYLINDRICAL GRIDS,
C POLAR STEREOGRAPHIC GRIDS CENTERED AT THE POLE
C AND MERCATOR GRIDS.
C ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
C THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS
C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63).
C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS:
C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL
C (KGDS(1)=001) MERCATOR CYLINDRICAL
C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL
C (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE)
C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL
C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE)
C WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO.
C AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
C AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
C ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
C IF KGDSO(1)<0, IN WHICH CASE THE NUMBER OF POINTS
C AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT.
C OUTPUT BITMAPS WILL NOT BE CREATED.
C
C PROGRAM HISTORY LOG:
C 96-04-10 IREDELL
C
C USAGE: CALL POLATES4(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI,
C & NO,RLAT,RLON,IBO,LO,GO,IRET)
C
C INPUT ARGUMENT LIST:
C IPOPT - INTEGER (20) INTERPOLATION OPTIONS
C IPOPT(1)=0 FOR TRIANGULAR, IPOPT(1)=1 FOR RHOMBOIDAL;
C IPOPT(2) IS TRUNCATION NUMBER
C (DEFAULTS TO SENSIBLE IF IPOPT(2)=-1).
C KGDSI - INTEGER (22) INPUT GDS PARAMETERS AS DECODED BY W3FI63
C KGDSO - INTEGER (22) OUTPUT GDS PARAMETERS AS DECODED BY W3FI63
C (KGDSO(1)<0 IMPLIES RANDOM STATION POINTS)
C MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
C OR DIMENSION OF INPUT GRID FIELDS IF KM=1
C MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
C OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
C KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE
C IBI - INTEGER (KM) INPUT BITMAP FLAGS (MUST BE ALL 0)
C LI - LOGICAL (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1)
C GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0)