-
Notifications
You must be signed in to change notification settings - Fork 17
/
Copy pathg2get.F90
1679 lines (1552 loc) · 61.6 KB
/
g2get.F90
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
!> @file
!> @brief Get information from a GRIB2 message.
!> @author Edward Hartnett @date Mar 6, 2024
!> Find the number of gridded fields and Local Use Sections in a GRIB2
!> message. Also various checks are performed to see if the message is a
!> valid GRIB2 message.
!>
!> This function is similar to gribinfo(), but returns less information.
!>
!> @param[in] cgrib Character array that contains the GRIB2 message.
!> @param[in] lcgrib Length (in bytes) of GRIB message in array
!> cgrib.
!> @param[out] listsec0 Contains information decoded from GRIB
!> Indicator Section 0. Must be dimensioned >= 2.
!> - listsec0(1) Discipline-GRIB Master Table Number (see [Code Table
!> - 0.0](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table0-0.shtml))
!> - listsec0(2) GRIB Edition Number (currently 2)
!> - listsec0(3) Length of GRIB message
!> @param[out] listsec1 Contains information read from GRIB
!> Identification Section 1. Must be dimensioned >= 13.
!> - listsec1(1) Id of orginating centre (Common Code Table C-1)
!> - listsec1(2) Id of orginating sub-centre (local table)
!> - listsec1(3) GRIB Master Tables Version Number (Code Table 1.0)
!> - listsec1(4) GRIB Local Tables Version Number
!> - listsec1(5) Significance of Reference Time (Code Table 1.1)
!> - listsec1(6) Reference Time - Year (4 digits)
!> - listsec1(7) Reference Time - Month
!> - listsec1(8) Reference Time - Day
!> - listsec1(9) Reference Time - Hour
!> - listsec1(10) Reference Time - Minute
!> - listsec1(11) Reference Time - Second
!> - listsec1(12) Production status of data (Code Table 1.2)
!> - listsec1(13) Type of processed data (Code Table 1.3)
!> @param[out] numfields The number of gridded fieldse found in the
!> GRIB message.
!> @param[out] numlocal The number of Local Use Sections (Section 2)
!> found in the GRIB message.
!> @param[out] maxlocal The size of the largest Local Use Section
!> (Section 2). Can be used to ensure that the return array passed
!> to subroutine getlocal is dimensioned large enough.
!> @param[out] ierr Error return code.
!> - 0 no error.
!> - 1 Beginning characters "GRIB" not found.
!> - 2 GRIB message is not Edition 2.
!> - 3 Could not find Section 1, where expected.
!> - 4 End string "7777" found, but not where expected.
!> - 5 End string "7777" not found at end of message.
!> - 6 Invalid section number found.
!>
!> @author Stephen Gilbert @date 2000-05-25
subroutine gb_info(cgrib, lcgrib, listsec0, listsec1, &
numfields, numlocal, maxlocal, ierr)
use g2logging
implicit none
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib
integer, intent(out) :: listsec0(3), listsec1(13)
integer, intent(out) :: numlocal, numfields, maxlocal, ierr
character(len = 4), parameter :: grib = 'GRIB', c7777 = '7777'
character(len = 4) :: ctemp
integer, parameter :: zero = 0, one = 1
integer, parameter :: mapsec1len = 13
integer, parameter :: mapsec1(mapsec1len) = (/ 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1 /)
integer :: iofst, istart
integer :: nbits, lensec1, lensec0, lensec, lenposs, lengrib, j
integer (kind = 8) :: lengrib8
integer :: i, ipos, isecnum
interface
subroutine g2_gbytec(in, iout, iskip, nbits)
character*1, intent(in) :: in(*)
integer, intent(inout) :: iout(*)
integer, intent(in) :: iskip, nbits
end subroutine g2_gbytec
subroutine g2_gbytec1(in, siout, iskip, nbits)
character*1, intent(in) :: in(*)
integer, intent(inout) :: siout
integer, intent(in) :: iskip, nbits
end subroutine g2_gbytec1
subroutine g2_gbytec81(in, siout, iskip, nbits)
character*1, intent(in) :: in(*)
integer (kind = 8), intent(inout) :: siout
integer, intent(in) :: iskip, nbits
integer (kind = 8) :: iout(1)
end subroutine g2_gbytec81
end interface
#ifdef LOGGING
write(g2_log_msg, *) 'gb_info: lcgrib ', lcgrib
call g2_log(1)
#endif
ierr = 0
numlocal = 0
numfields = 0
maxlocal = 0
! Check for beginning of GRIB message in the first 100 bytes.
istart = 0
do j = 1, 100
ctemp = cgrib(j) // cgrib(j + 1) // cgrib(j + 2) // cgrib(j+3)
if (ctemp .eq. grib ) then
istart = j
exit
endif
enddo
if (istart .eq. 0) then
print *, 'gb_info: Beginning characters GRIB not found.'
ierr = 1
return
endif
! Unpack Section 0 - Indicator Section.
iofst = 8 * (istart + 5)
call g2_gbytec(cgrib, listsec0(1), iofst, 8) ! Discipline
iofst = iofst + 8
call g2_gbytec(cgrib, listsec0(2), iofst, 8) ! GRIB edition number
iofst = iofst + 8
#ifdef LOGGING
! Log results for debugging.
write(g2_log_msg, *) 'before getting len: iofst/8 ', iofst/8
call g2_log(2)
#endif
! Bytes 9-16 contain the length of GRIB message. This is an 8-byte
! value, but unfortunately we only have a 4-byte subroutine
! parameter for it.
call g2_gbytec81(cgrib, lengrib8, iofst, 64)
lengrib = int(lengrib8, kind(4))
iofst = iofst + 32
iofst = iofst + 32
listsec0(3) = lengrib
lensec0 = 16
ipos = istart + lensec0
! The g2 library only handles GRIB2.
if (listsec0(2) .ne. 2) then
print *, 'gb_info: can only decode GRIB edition 2.'
ierr = 2
return
endif
! Unpack Section 1 - Identification Section.
call g2_gbytec1(cgrib, lensec1, iofst, 32) ! Length of Section 1
iofst = iofst + 32
call g2_gbytec1(cgrib, isecnum, iofst, 8) ! Section number ( 1 )
iofst = iofst + 8
if (isecnum .ne. 1) then
print *, 'gb_info: Could not find section 1.'
ierr = 3
return
endif
! Unpack each input value in array listsec1 into the the appropriate
! number of octets, which are specified in corresponding entries in
! array mapsec1.
do i = 1, mapsec1len
nbits = mapsec1(i) * 8
call g2_gbytec(cgrib, listsec1(i), iofst, nbits)
iofst = iofst + nbits
enddo
ipos = ipos + lensec1
! Loop through the remaining sections to see if they are valid.
! Also count the number of times Section 2 and Section 4 appear.
do
ctemp = cgrib(ipos) // cgrib(ipos + 1) // cgrib(ipos + 2) // cgrib(ipos + 3)
if (ctemp .eq. c7777 ) then
ipos = ipos + 4
if (ipos .ne. (istart + lengrib)) then
print *, 'gb_info: "7777" found, but not where expected.'
ierr = 4
return
endif
exit
endif
iofst = (ipos - 1) * 8
call g2_gbytec1(cgrib, lensec, iofst, 32) ! Get Length of Section
iofst = iofst + 32
call g2_gbytec1(cgrib, isecnum, iofst, 8) ! Get Section number
iofst = iofst + 8
ipos = ipos+lensec ! Update beginning of section pointer
if (ipos .gt. (istart + lengrib)) then
print *, 'gb_info: "7777" not found at end of GRIB message.'
ierr = 5
return
endif
if (isecnum .ge. 2 .AND. isecnum .le. 7) then
if (isecnum .eq. 2) then ! Local Section 2
! Increment counter for total number of local sections found.
numlocal = numlocal + 1
lenposs = lensec - 5
if (lenposs .gt. maxlocal) maxlocal = lenposs
elseif (isecnum .eq. 4) then
! Increment counter for total number of fields found.
numfields = numfields + 1
endif
else
print *, 'gb_info: Invalid section number found in GRIB message: ', isecnum
ierr = 6
return
endif
enddo
end subroutine gb_info
!> Find the number of Local Use Sections and gridded fields in
!> a GRIB2 message, and the maximum sizes of template arrays.
!>
!> This subroutine also performs various checks to see if the message
!> is a valid GRIB2 message. Also, a list of safe array dimensions is
!> returned for use in allocating return arrays from routines
!> getlocal(), gettemplates(), and getfields().
!>
!> Array maxvals contains the maximum possible number of values
!> that will be returned in argument arrays for routines getlocal(),
!> gettemplates() and getfields(). Users can use this info to determine
!> if their arrays are dimensioned large enough for the data that may
!> be returned from the above routines, or to dynamically allocate
!> arrays with a reasonable size.
!>
!> This function is similar to gb_info(), but returns more information.
!>
!> @note The actual number of values in these arrays will likely be
!> less than the values calculated by this routine.
!>
!> @param[in] cgrib Character that contains the GRIB2 message.
!> @param[in] lcgrib Length (in bytes) of array cgrib.
!> @param[out] listsec0 Contains information needed for GRIB
!> Indicator Section 0. Must be dimensioned >= 2.
!> - listsec0(1) Discipline-GRIB Master Table Number ([Code Table 0.0]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table0-0.shtml)).
!> - listsec0(2) GRIB Edition Number (currently 2)
!> - listsec0(3) Length of GRIB message.
!> @param[out] listsec1 Contains information needed for GRIB
!> Identification Section 1. Must be dimensioned >= 13.
!> - listsec1(1) Id of orginating centre ([Table 0]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/on388/table0.html)).
!> - listsec1(2) Id of orginating sub-centre ([Table C]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/on388/tablec.html)).
!> - listsec1(3) GRIB Master Tables Version Number ([Table 1.0]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-0.shtml)).
!> - listsec1(4) GRIB Local Tables Version Number ([Table 1.1]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-1.shtml)).
!> - listsec1(5) Significance of Reference Time ([Table 1.2]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-1.shtml))
!> - listsec1(6) Reference Time - Year (4 digits)
!> - listsec1(7) Reference Time - Month
!> - listsec1(8) Reference Time - Day
!> - listsec1(9) Reference Time - Hour
!> - listsec1(10) Reference Time - Minute
!> - listsec1(11) Reference Time - Second
!> - listsec1(12) Production status of data ([Table 1.3]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-3.shtml)).
!> - listsec1(13) Type of processed data ([Table 1.4]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-4.shtml)).
!> @param[out] numlocal The number of Local Use Sections (Section 2)
!> found in the GRIB message.
!> @param[out] numfields The number of gridded fieldse found in the
!> GRIB message.
!> @param[out] maxvals The maximum number of elements that could be
!> returned in various arrays from this GRIB2 message.
!> - maxvals(1) max length of local section 2 (for getlocal()).
!> - maxvals(2) max length of GDS Template (for gettemplates() and getfield()).
!> - maxvals(3) max length of GDS Optional list (for getfield()).
!> - maxvals(4) max length of PDS Template (for gettemplates() and getfield()).
!> - maxvals(5) max length of PDS Optional list (for getfield()).
!> - maxvals(6) max length of DRS Template (for gettemplates() and getfield()).
!> - maxvals(7) max number of gridpoints (for getfield()).
!> @param[out] ierr Error return code.
!> - 0 no error.
!> - 1 Beginning characters "GRIB" not found.
!> - 2 GRIB message is not Edition 2.
!> - 3 Could not find Section 1, where expected.
!> - 4 End string "7777" found, but not where expected.
!> - 5 End string "7777" not found at end of message.
!>
!> @author Stephen Gilbert @date 2000-05-25
subroutine gribinfo(cgrib, lcgrib, listsec0, listsec1, &
numlocal, numfields, maxvals, ierr)
implicit none
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib
integer, intent(out) :: listsec0(3), listsec1(13), maxvals(7)
integer, intent(out) :: numlocal, numfields, ierr
integer :: i, ipos, isecnum, j, lengrib, lenposs, lensec, lensec0, lensec1
integer :: maxcoordlist, maxdeflist, maxdrstmpl, maxgridpts, maxpdstmpl, maxgdstmpl
integer :: maxsec2len, nbits, nbyte, ngdpts, numcoord
character(len = 4), parameter :: grib = 'GRIB', c7777 = '7777'
character(len = 4) :: ctemp
integer, parameter :: zero = 0, one = 1
integer, parameter :: mapsec1len = 13
integer, parameter :: mapsec1(mapsec1len) = (/ 2, 2, 1, 1, 1, 2, 1, 1, &
1, 1, 1, 1, 1 /)
integer iofst, istart
integer (kind = 8) :: lengrib8
ierr = 0
numlocal = 0
numfields = 0
maxsec2len = 1
maxgdstmpl = 1
maxdeflist = 1
maxpdstmpl = 1
maxcoordlist = 1
maxdrstmpl = 1
maxgridpts = 0
! Check for beginning of GRIB message in the first 100 bytes.
istart = 0
do j = 1, 100
ctemp = cgrib(j) // cgrib(j + 1) // cgrib(j + 2) // cgrib(j + 3)
if (ctemp .eq. grib) then
istart = j
exit
endif
enddo
if (istart .eq. 0) then
print *, 'gribinfo: Beginning characters GRIB not found.'
ierr = 1
return
endif
! Unpack Section 0 - Indicator Section.
iofst = 8 * (istart + 5)
call g2_gbytec(cgrib, listsec0(1), iofst, 8) ! Discipline
iofst = iofst + 8
call g2_gbytec(cgrib, listsec0(2), iofst, 8) ! GRIB edition number
iofst = iofst + 8
! Bytes 9-16 contain the length of GRIB message. This is an 8-byte
! value, but unfortunately we only have a 4-byte subroutine
! parameter for it.
call g2_gbytec81(cgrib, lengrib8, iofst, 64) ! Length of GRIB message
lengrib = int(lengrib8, kind(4))
iofst = iofst + 32
iofst = iofst + 32
listsec0(3) = lengrib
lensec0 = 16
ipos = istart + lensec0
! Currently handles only GRIB Edition 2.
if (listsec0(2) .ne. 2) then
print *, 'gribinfo: can only decode GRIB edition 2.'
ierr = 2
return
endif
! Unpack Section 1 - Identification Section.
call g2_gbytec1(cgrib, lensec1, iofst, 32) ! Length of Section 1
iofst = iofst + 32
call g2_gbytec1(cgrib, isecnum, iofst, 8) ! Section number (1)
iofst = iofst + 8
if (isecnum .ne. 1) then
print *, 'gribinfo: Could not find section 1.'
ierr = 3
return
endif
! Unpack each input value in array listsec1 into the the appropriate
! number of octets, which are specified in corresponding entries in
! array mapsec1.
do i = 1, mapsec1len
nbits = mapsec1(i) * 8
call g2_gbytec(cgrib, listsec1(i), iofst, nbits)
iofst = iofst + nbits
enddo
ipos = ipos + lensec1
! Loop through the remaining sections keeping track of the length of
! each. Also count the number of times Section 2 and Section 4
! appear.
do
ctemp = cgrib(ipos) // cgrib(ipos + 1) // cgrib(ipos + 2) // &
cgrib(ipos + 3)
if (ctemp .eq. c7777) then
ipos = ipos + 4
if (ipos .ne. (istart + lengrib)) then
print *, 'gribinfo: "7777" found, but not where expected.'
ierr = 4
return
endif
exit
endif
iofst = (ipos - 1) * 8
call g2_gbytec1(cgrib, lensec, iofst, 32) ! Get Length of Section
iofst = iofst + 32
call g2_gbytec1(cgrib, isecnum, iofst, 8) ! Get Section number
iofst = iofst + 8
ipos = ipos + lensec ! Update beginning of section pointer
if (ipos .gt. (istart + lengrib)) then
print *, 'gribinfo: "7777" not found at end of GRIB message.'
ierr = 5
return
endif
if (isecnum .eq. 2) then ! Local Section 2
! Increment counter for total number of local sections found
! and determine largest Section 2 in message.
numlocal = numlocal + 1
lenposs = lensec-5
if (lenposs .gt. maxsec2len) maxsec2len = lenposs
elseif (isecnum .eq. 3) then
iofst = iofst + 8 ! skip source of grid def.
call g2_gbytec1(cgrib, ngdpts, iofst, 32) ! Get Num of Grid Points
iofst = iofst + 32
call g2_gbytec1(cgrib, nbyte, iofst, 8) ! Get Num octets for opt. list
iofst = iofst + 8
if (ngdpts .gt. maxgridpts) maxgridpts = ngdpts
lenposs = lensec - 14
if (lenposs .gt. maxgdstmpl) maxgdstmpl = lenposs
if (nbyte .ne. 0) then
lenposs = lenposs / nbyte
if (lenposs .gt. maxdeflist) maxdeflist = lenposs
endif
elseif (isecnum .eq. 4) then
numfields = numfields + 1
call g2_gbytec1(cgrib, numcoord, iofst, 16) ! Get Num of Coord Values
iofst = iofst + 16
if (numcoord .ne. 0) then
if (numcoord .gt. maxcoordlist) maxcoordlist = numcoord
endif
lenposs = lensec - 9
if (lenposs.gt.maxpdstmpl) maxpdstmpl = lenposs
elseif (isecnum .eq. 5) then
lenposs = lensec-11
if (lenposs .gt. maxdrstmpl) maxdrstmpl = lenposs
endif
enddo
maxvals(1) = maxsec2len
maxvals(2) = maxgdstmpl
maxvals(3) = maxdeflist
maxvals(4) = maxpdstmpl
maxvals(5) = maxcoordlist
maxvals(6) = maxdrstmpl
maxvals(7) = maxgridpts
end subroutine gribinfo
!> Return the Grid Definition, Product Definition, Bit-map (if
!> applicable), and the unpacked data for a data field. Since there
!> can be multiple data fields packed into a GRIB2 message, the
!> calling routine indicates which field is being requested with the
!> ifldnum argument.
!>
!> @note Note that subroutine gribinfo() can be used to first
!> determine how many data fields exist in a given GRIB message.
!>
!> @param[in] cgrib Character array that contains the GRIB2 message.
!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib.
!> @param[in] ifldnum Specifies which field in the GRIB2 message to
!> return.
!> @param[out] igds Contains information read from the appropriate
!> GRIB Grid Definition Section 3 for the field being returned. Must
!> be dimensioned >= 5.
!> - igds(1) Source of grid definition (see [Code Table
!> 3.0](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-0.shtml)).
!> - igds(2) Number of grid points in the defined grid.
!> - igds(3) Number of octets needed for each additional grid points
!> definition. Used to define number of points in each row (or
!> column) for non-regular grids. = 0, if using regular grid.
!> - igds(4) Interpretation of list for optional points definition
!> ([Code Table
!> 3.11](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-11.shtml)).
!> - igds(5) Grid Definition Template Number ([Code Table
!> 3.1](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-1.shtml)).
!> @param[out] igdstmpl Contains the data values for the specified
!> Grid Definition Template (NN=igds(5)). Each element of this
!> integer array contains an entry (in the order specified) of Grid
!> Defintion Template 3.NN. A safe dimension for this array can be
!> obtained in advance from maxvals(2), which is returned from
!> subroutine gribinfo().
!> @param[out] igdslen Number of elements in igdstmpl. i.e. number of
!> entries in Grid Defintion Template 3.NN (NN=igds(5)).
!> @param[out] ideflist (Used if igds(3) .ne. 0) This array contains
!> the number of grid points contained in each row (or column). (part
!> of Section 3) A safe dimension for this array can be obtained in
!> advance from maxvals(3), which is returned from subroutine
!> gribinfo().
!> @param[out] idefnum (Used if igds(3) .ne. 0) The number of entries
!> in array ideflist - i.e. number of rows (or columns) for which
!> optional grid points are defined.
!> @param[out] ipdsnum Product Definition Template Number (see [Code
!> Table 4.0]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-0.shtml)).
!> @param[out] ipdstmpl Contains the data values for the specified
!> Product Definition Template (N=ipdsnum). Each element of this
!> integer array contains an entry (in the order specified) of
!> Product Defintion Template 4.N. A safe dimension for this array
!> can be obtained in advance from maxvals(4), which is returned from
!> subroutine gribinfo().
!> @param[out] ipdslen Number of elements in ipdstmpl - i.e. number
!> of entries in Product Defintion Template 4.N (N=ipdsnum).
!> @param[out] coordlist Array containg floating point values
!> intended to document the vertical discretisation associated to
!> model data on hybrid coordinate vertical levels (part of Section
!> 4). The dimension of this array can be obtained in advance from
!> maxvals(5), which is returned from subroutine gribinfo().
!> @param[out] numcoord number of values in array coordlist.
!> @param[out] ndpts Number of data points unpacked and returned.
!> @param[out] idrsnum Data Representation Template Number (see [Code
!> Table 5.0]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-0.shtml)).
!> @param[out] idrstmpl Contains the data values for the specified
!> Data Representation Template (N=idrsnum). Each element of this
!> integer array contains an entry (in the order specified) of
!> Product Defintion Template 5.N A safe dimension for this array can
!> be obtained in advance from maxvals(6), which is returned from
!> subroutine gribinfo().
!> @param[out] idrslen Number of elements in idrstmpl. i.e. number of
!> entries in Data Representation Template specified by idrsnum.
!> @param[out] ibmap Bitmap indicator (see [Code Table
!> 6.0](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table6-0.shtml)).
!> - 0 bitmap applies and is included in Section 6.
!> - 1-253 Predefined bitmap applies.
!> - 254 Previously defined bitmap applies to this field.
!> - 255 Bit map does not apply to this product.
!> @param[out] bmap Logical*1 array containing decoded bitmap (if
!> ibmap=0). The dimension of this array can be obtained in advance
!> from maxvals(7), which is returned from subroutine gribinfo().
!> @param[out] fld Array of ndpts unpacked data points. A safe
!> dimension for this array can be obtained in advance from
!> maxvals(7), which is returned from subroutine gribinfo().
!> @param[out] ierr Error return code.
!> - 0 no error.
!> - 1 Beginning characters "GRIB" not found.
!> - 2 GRIB message is not Edition 2.
!> - 3 The data field request number was not positive.
!> - 4 End string "7777" found, but not where expected.
!> - 6 GRIB message did not contain the requested number of data fields.
!> - 7 End string "7777" not found at end of message.
!> - 9 Data Representation Template 5.NN not yet implemented.
!> - 10 Error unpacking Section 3.
!> - 11 Error unpacking Section 4.
!> - 12 Error unpacking Section 5.
!> - 13 Error unpacking Section 6.
!> - 14 Error unpacking Section 7.
!>
!> @author Stephen Gilbert @date 2000-05-26
subroutine getfield(cgrib, lcgrib, ifldnum, igds, igdstmpl, &
igdslen, ideflist, idefnum, ipdsnum, ipdstmpl, ipdslen, &
coordlist, numcoord, ndpts, idrsnum, idrstmpl, idrslen, &
ibmap, bmap, fld, ierr)
implicit none
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib, ifldnum
integer, intent(out) :: igds(*), igdstmpl(*), ideflist(*)
integer, intent(out) :: ipdsnum, ipdstmpl(*)
integer, intent(out) :: idrsnum, idrstmpl(*)
integer, intent(out) :: ndpts, ibmap, idefnum, numcoord
integer, intent(out) :: igdslen, ipdslen, idrslen
real, intent(out) :: fld(*), coordlist(*)
logical*1, intent(out) :: bmap(*)
integer, intent(out) :: ierr
character(len = 4), parameter :: grib = 'GRIB', c7777 = '7777'
character(len = 4) :: ctemp
integer :: listsec0(2)
integer :: iofst, istart
real (kind = 4) :: ieee(1)
logical :: have3, have4, have5, have6, have7
integer (kind = 8) :: lengrib8
integer :: numfld, j, lengrib, lensec0, ipos
integer :: lensec, isecnum, jerr, ier, numlocal
have3 = .false.
have4 = .false.
have5 = .false.
have6 = .false.
have7 = .false.
ierr = 0
numfld = 0
numlocal = 0
! Check for valid request number.
if (ifldnum .le. 0) then
print *, 'getfield: Request for field number ' &
,'must be positive.'
ierr = 3
return
endif
! Check for beginning of GRIB message in the first 100 bytes.
istart = 0
do j = 1, 100
ctemp = cgrib(j) // cgrib(j + 1) // cgrib(j + 2) // cgrib(j + 3)
if (ctemp .eq. grib) then
istart = j
exit
endif
enddo
if (istart .eq. 0) then
print *, 'getfield: Beginning characters GRIB not found.'
ierr = 1
return
endif
! Unpack Section 0 - Indicator Section.
iofst = 8 * (istart + 5)
call g2_gbytec(cgrib, listsec0(1), iofst, 8) ! Discipline
iofst = iofst + 8
call g2_gbytec(cgrib, listsec0(2), iofst, 8) ! GRIB edition number
iofst = iofst + 8
! Bytes 9-16 contain the length of GRIB message. This is an 8-byte
! value, but unfortunately we only have a 4-byte subroutine
! parameter for it.
call g2_gbytec81(cgrib, lengrib8, iofst, 64) ! Length of GRIB message
lengrib = int(lengrib8, kind(4))
iofst = iofst + 32
iofst = iofst + 32
lensec0 = 16
ipos = istart + lensec0
! The g2 library only handles GRIB2.
if (listsec0(2) .ne. 2) then
print *, 'getfield: can only decode GRIB edition 2.'
ierr = 2
return
endif
! Loop through the remaining sections keeping track of the length of
! each. Also keep the latest Grid Definition Section info. Unpack
! the requested field number.
do
! Check to see if we are at end of GRIB message
ctemp = cgrib(ipos) // cgrib(ipos + 1) // cgrib(ipos + 2) // &
cgrib(ipos + 3)
if (ctemp .eq. c7777) then
ipos = ipos + 4
! If end of GRIB message not where expected, issue error
if (ipos.ne.(istart + lengrib)) then
print *, 'getfield: "7777" found, but not ' &
,'where expected.'
ierr = 4
return
endif
exit
endif
! Get length of Section and Section number
iofst = (ipos - 1) * 8
call g2_gbytec1(cgrib, lensec, iofst, 32) ! Get Length of Section
iofst = iofst + 32
call g2_gbytec1(cgrib, isecnum, iofst, 8) ! Get Section number
iofst = iofst + 8
! If found Section 3, unpack the GDS info using the appropriate
! template. Save in case this is the latest grid before the
! requested field.
if (isecnum .eq. 3) then
iofst = iofst - 40 ! reset offset to beginning of section
call unpack3(cgrib, lcgrib, iofst, igds, igdstmpl, &
igdslen, ideflist, idefnum, jerr)
if (jerr .eq. 0) then
have3 = .true.
else
ierr = 10
return
endif
endif
! If found Section 4, check to see if this field is the one
! requested.
if (isecnum .eq. 4) then
numfld = numfld + 1
if (numfld .eq. ifldnum) then
iofst = iofst - 40 ! reset offset to beginning of section
call unpack4(cgrib, lcgrib, iofst, ipdsnum, ipdstmpl, &
ipdslen, coordlist, numcoord, jerr)
if (jerr .eq. 0) then
have4 = .true.
else
ierr = 11
return
endif
endif
endif
! If found Section 5, check to see if this field is the one
! requested.
if ((isecnum .eq. 5) .and. (numfld .eq. ifldnum)) then
iofst = iofst - 40 ! reset offset to beginning of section
call unpack5(cgrib, lcgrib, iofst, ndpts, idrsnum, &
idrstmpl, idrslen, jerr)
if (jerr .eq. 0) then
have5 = .true.
else
ierr = 12
return
endif
endif
! If found Section 6, Unpack bitmap. Save in case this is the
! latest bitmap before the requested field.
if (isecnum .eq. 6) then
iofst = iofst - 40 ! reset offset to beginning of section
call unpack6(cgrib, lcgrib, iofst, igds(2), ibmap, bmap, &
jerr)
if (jerr .eq. 0) then
have6 = .true.
else
ierr = 13
return
endif
endif
! If found Section 7, check to see if this field is the one
! requested.
if ((isecnum .eq. 7) .and. (numfld .eq. ifldnum)) then
if (idrsnum .eq. 0) then
call simunpack(cgrib(ipos + 5), lensec - 6, idrstmpl, &
ndpts, fld)
have7 = .true.
elseif (idrsnum .eq. 2 .or. idrsnum .eq. 3) then
call comunpack(cgrib(ipos + 5), lensec - 6, lensec, &
idrsnum,idrstmpl, ndpts, fld, ier)
if (ier .ne. 0) then
ierr = 14
return
endif
have7 = .true.
elseif (idrsnum .eq. 50) then
call simunpack(cgrib(ipos + 5), lensec - 6, idrstmpl, &
ndpts - 1, fld(2))
ieee = transfer(idrstmpl(5), ieee, 1)
call rdieee(ieee, fld(1), 1)
have7 = .true.
elseif (idrsnum .eq. 40 .or. idrsnum .eq. 40000) then
call jpcunpack(cgrib(ipos + 5), lensec - 5, idrstmpl, &
ndpts, fld)
have7 = .true.
elseif (idrsnum .eq. 41 .or. idrsnum .eq. 40010) then
call pngunpack(cgrib(ipos + 5), lensec - 5, idrstmpl, &
ndpts, fld)
have7 = .true.
#ifdef USE_AEC
elseif (idrsnum .eq. 42) then
call aecunpack(cgrib(ipos + 5), lensec - 5, idrstmpl, &
ndpts, fld)
have7 = .true.
#endif /* USE_AEC */
else
print *, 'getfield: Data Representation Template ', &
idrsnum, ' not yet implemented.'
ierr = 9
return
endif
endif
! Check to see if we read pass the end of the GRIB message and
! missed the terminator string '7777'.
ipos = ipos + lensec ! Update beginning of section pointer
if (ipos .gt. (istart + lengrib)) then
print *, 'getfield: "7777" not found at end' &
,' of GRIB message.'
ierr = 7
return
endif
if (have3 .and. have4 .and. have5 .and. have6 .and. have7) &
return
enddo
! If exited from above loop, the end of the GRIB message was reached
! before the requested field was found.
print *, 'getfield: GRIB message contained ', numlocal, &
' different fields.'
print *, 'getfield: The request was for the ', ifldnum, &
' field.'
ierr = 6
end subroutine getfield
!> This subroutine unpacks Section 3 (Grid Definition Section)
!> starting at octet 6 of that Section.
!>
!> @param[in] cgrib Character array that contains the GRIB2 message.
!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib.
!> @param[inout] iofst Bit offset of the beginning (in) or the end
!> (out) of Section 3.
!> @param[out] igds Contains information read from the appropriate
!> GRIB Grid Definition Section 3 for the field being returned. Must
!> be dimensioned >= 5.
!> - igds(1) Source of grid definition (see [Code Table - 3.0]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-0.shtml))
!> - igds(2) Number of grid points in the defined grid.
!> - igds(3) Number of octets needed for each additional grid points
!> definition. Used to define number of points in each row (or
!> column) for non-regular grids. = 0, if using regular grid.
!> - igds(4) Interpretation of list for optional points
!> definition. ([Code Table 3.11]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-11.shtml)).
!> - igds(5) Grid Definition Template Number ([Code Table 3.1]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-1.shtml)).
!> @param[out] igdstmpl Contains the data values for the specified
!> Grid Definition Template (NN=igds(5)). Each element of this
!> integer array contains an entry (in the order specified) of Grid
!> Defintion Template 3.NN.
!> @param[out] mapgridlen Number of elements in igdstmpl -
!> i.e. number of entries in Grid Defintion Template 3.NN
!> (NN=igds(5)).
!> @param[out] ideflist (Used if igds(3) .ne. 0). This array contains
!> the number of grid points contained in each row (or column) (part
!> of Section 3).
!> @param[out] idefnum (Used if igds(3) .ne. 0). The number of
!> entries in array ideflist - i.e. number of rows (or columns) for
!> which optional grid points are defined.
!> @param[out] ierr Error return code.
!> - 0 no error.
!> - 5 "GRIB" message contains an undefined Grid Definition Template.
!>
!> @author Stephen Gilbert @date 2000-05-26
subroutine unpack3(cgrib, lcgrib, iofst, igds, igdstmpl, &
mapgridlen, ideflist, idefnum, ierr)
use gridtemplates
implicit none
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib
integer, intent(inout) :: iofst
integer, intent(out) :: igds(*), igdstmpl(*), ideflist(*)
integer, intent(out) :: ierr, idefnum
integer, allocatable :: mapgrid(:)
integer :: mapgridlen, ibyttem
logical needext
integer :: lensec, iret, i, nbits, isign, newmapgridlen
ierr = 0
call g2_gbytec1(cgrib, lensec, iofst, 32) ! Get Length of Section
iofst = iofst + 32
iofst = iofst + 8 ! skip section number
call g2_gbytec1(cgrib, igds(1), iofst, 8) ! Get source of Grid def.
iofst = iofst + 8
call g2_gbytec1(cgrib, igds(2), iofst, 32) ! Get number of grid pts.
iofst = iofst + 32
call g2_gbytec1(cgrib, igds(3), iofst, 8) ! Get num octets for opt. list
iofst = iofst + 8
call g2_gbytec1(cgrib, igds(4), iofst, 8) ! Get interpret. for opt. list
iofst = iofst + 8
call g2_gbytec1(cgrib, igds(5), iofst, 16) ! Get Grid Def Template num.
iofst = iofst + 16
if (igds(1) .eq. 0) then
allocate(mapgrid(lensec))
! Get Grid Definition Template
call getgridtemplate(igds(5), mapgridlen, mapgrid, needext, &
iret)
if (iret .ne. 0) then
ierr = 5
return
endif
else
mapgridlen = 0
needext = .false.
endif
! Unpack each value into array igdstmpl from the the appropriate
! number of octets, which are specified in corresponding entries in
! array mapgrid.
ibyttem = 0
do i = 1, mapgridlen
nbits = iabs(mapgrid(i)) * 8
if (mapgrid(i) .ge. 0) then
call g2_gbytec1(cgrib, igdstmpl(i), iofst, nbits)
else
call g2_gbytec1(cgrib, isign, iofst, 1)
call g2_gbytec1(cgrib, igdstmpl(i), iofst + 1, nbits-1)
if (isign .eq. 1) igdstmpl(i) = -igdstmpl(i)
endif
iofst = iofst + nbits
ibyttem = ibyttem + iabs(mapgrid(i))
enddo
! Check to see if the Grid Definition Template needs to be
! extended. The number of values in a specific template may vary
! depending on data specified in the "static" part of the template.
if (needext) then
call extgridtemplate(igds(5), igdstmpl, newmapgridlen, &
mapgrid)
! Unpack the rest of the Grid Definition Template
do i = mapgridlen + 1, newmapgridlen
nbits = iabs(mapgrid(i)) * 8
if (mapgrid(i) .ge. 0) then
call g2_gbytec1(cgrib, igdstmpl(i), iofst, nbits)
else
call g2_gbytec1(cgrib, isign, iofst, 1)
call g2_gbytec1(cgrib, igdstmpl(i), iofst + 1, nbits - &
1)
if (isign .eq. 1) igdstmpl(i) = -igdstmpl(i)
endif
iofst = iofst + nbits
ibyttem = ibyttem + iabs(mapgrid(i))
enddo
mapgridlen = newmapgridlen
endif
! Unpack optional list of numbers defining number of points in each
! row or column, if included. This is used for non regular grids.
if (igds(3) .ne. 0) then
nbits = igds(3) * 8
idefnum = (lensec - 14 - ibyttem) / igds(3)
call g2_gbytesc(cgrib, ideflist, iofst, nbits, 0, idefnum)
iofst = iofst + (nbits * idefnum)
else
idefnum = 0
endif
if (allocated(mapgrid)) deallocate(mapgrid)
end subroutine unpack3
!> This subroutine unpacks Section 4 (Product Definition Section)
!> starting at octet 6 of that Section.
!>
!> @param[in] cgrib Character array that contains the GRIB2 message.
!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib.
!> @param[inout] iofst Bit offset of the beginning (in) or the end
!> (out) of Section 4.
!> @param[out] ipdsnum Product Definition Template Number (see [Code
!> Table 4.0]
!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-0.shtml)).
!> @param[out] ipdstmpl Contains the data values for the specified
!> Product Definition Template (N=ipdsnum). Each element of this
!> integer array contains an entry (in the order specified) of
!> Product Defintion Template 4.N.
!> @param[out] mappdslen Number of elements in ipdstmpl. i.e. number
!> of entries in Product Defintion Template 4.N (N=ipdsnum).
!> @param[out] coordlist- Array containg floating point values
!> intended to document the vertical discretisation associated to
!> model data on hybrid coordinate vertical levels (part of Section
!> 4).
!> @param[out] numcoord number of values in array coordlist.
!> @param[out] ierr Error return code.
!> - 0 no error.
!> - 5 GRIB message contains an undefined Product Definition Template.
!>
!> @author Stephen Gilbert @date 2000-05-26
subroutine unpack4(cgrib, lcgrib, iofst, ipdsnum, ipdstmpl, &
mappdslen, coordlist, numcoord, ierr)
use pdstemplates
implicit none
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib
integer, intent(inout) :: iofst
real, intent(out) :: coordlist(*)
integer, intent(out) :: ipdsnum, ipdstmpl(*)
integer, intent(out) :: ierr, numcoord
real(4), allocatable :: coordieee(:)
integer, allocatable :: mappds(:)
integer :: mappdslen
logical needext
integer :: lensec, iret, i, nbits, isign, newmappdslen
ierr = 0
call g2_gbytec1(cgrib, lensec, iofst, 32) ! Get Length of Section
iofst = iofst + 32
iofst = iofst + 8 ! skip section number
allocate(mappds(lensec))
call g2_gbytec1(cgrib, numcoord, iofst, 16) ! Get num of coordinate values
iofst = iofst + 16
call g2_gbytec1(cgrib, ipdsnum, iofst, 16) ! Get Prod. Def Template num.
iofst = iofst + 16
! Get Product Definition Template.
call getpdstemplate(ipdsnum, mappdslen, mappds, needext, iret)
if (iret.ne.0) then
ierr = 5
return
endif
! Unpack each value into array ipdstmpl from the the appropriate
! number of octets, which are specified in corresponding entries in
! array mappds.
do i = 1, mappdslen
nbits = iabs(mappds(i)) * 8
if (mappds(i) .ge. 0) then
call g2_gbytec1(cgrib, ipdstmpl(i), iofst, nbits)
else
call g2_gbytec1(cgrib, isign, iofst, 1)
call g2_gbytec1(cgrib, ipdstmpl(i), iofst + 1, nbits-1)
if (isign .eq. 1) ipdstmpl(i) = -ipdstmpl(i)
endif
iofst = iofst + nbits
enddo
! Check to see if the Product Definition Template needs to be
! extended. The number of values in a specific template may vary
! depending on data specified in the "static" part of the template.
if (needext) then
call extpdstemplate(ipdsnum, ipdstmpl, newmappdslen, mappds)
! Unpack the rest of the Product Definition Template.