-
Notifications
You must be signed in to change notification settings - Fork 17
/
Copy pathg2index.F90
1657 lines (1529 loc) · 60.9 KB
/
g2index.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 Subroutines for dealing with indexes.
!> @author Edward Hartnett @date Jan 31, 2024
!> Create a version 1 or 2 index file for a GRIB2 file.
!>
!> @param[in] lugb Logical unit of opened GRIB2 file.
!> @param[in] lugi Logical unit file opened to write index to.
!> @param[in] idxver Index version.
!> @param[in] filename Name of GRIB2 file. This file must already be
!> open, with logical unit passed to the lugb parameter. The filename
!> is also encoded in one of the index header records.
!> @param[out] iret Return code:
!> - 0 success
!> - 90 problem opening GRIB2 file.
!> - 91 problem opening index file.
!> - 92 no messages found in GRIB2 file.
!>
!> @author Ed Hartnett, Mark Iredell @date Feb 15, 2024
subroutine g2_create_index(lugb, lugi, idxver, filename, iret)
implicit none
integer, intent(in) :: lugb, lugi, idxver
character*(*) :: filename
integer, intent(out) :: iret
integer (kind = 8) :: msk1, msk2
parameter(msk1 = 32000_8, msk2 = 4000_8)
character(len=1), pointer, dimension(:) :: cbuf
integer :: numtot, nnum, nlen, mnum, kw
integer :: irgi, iw, nmess
interface
subroutine getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, &
nlen, nnum, nmess, iret)
integer, intent(in) :: lugb
integer (kind = 8), intent(in) :: msk1, msk2
integer, intent(in) :: mnum, idxver
character(len = 1), pointer, dimension(:) :: cbuf
integer, intent(out) :: nlen, nnum, nmess, iret
end subroutine getg2i2r
end interface
! Assume success.
iret = 0
numtot = 0
nlen = 0
! Generate index records for all messages in file, or until memory
! runs out.
mnum = 0
call getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, &
nlen, nnum, nmess, irgi)
if (irgi .gt. 1 .or. nnum .eq. 0 .or. nlen .eq. 0) then
iret = 92
if (associated(cbuf)) then
deallocate(cbuf)
nullify(cbuf)
endif
return
endif
numtot = numtot + nnum
mnum = mnum + nmess
! Write headers.
call g2_write_index_headers(lugi, nlen, numtot, idxver, filename)
iw = 162
! Write the index data we have so far.
call bawrite(lugi, iw, nlen, kw, cbuf)
iw = iw + nlen
! Extend index file if index buffer length too large to hold in memory.
if (irgi .eq. 1) then
do while (irgi .eq. 1 .and. nnum .gt. 0)
if (associated(cbuf)) then
deallocate(cbuf)
nullify(cbuf)
endif
call getg2i2r(11, msk1, msk2, mnum, idxver, cbuf, &
nlen, nnum, nmess, irgi)
if (irgi .le. 1 .and. nnum .gt. 0) then
numtot = numtot + nnum
mnum = mnum + nmess
call bawrite(lugi, iw, nlen, kw, cbuf)
iw = iw + nlen
endif
enddo
! Go back and overwrite headers with new info.
call g2_write_index_headers(lugi, iw, numtot, idxver, filename)
endif
deallocate(cbuf)
end subroutine g2_create_index
!> Write index headers.
!>
!> @param[in] lugi integer logical unit of output index file
!> @param[in] nlen integer total length of index records
!> @param[in] nnum integer number of index records
!> @param[in] idxver Index version, 1 for legacy, 2 for files that may
!> be > 2 GB.
!> @param[in] filename character name of GRIB file
!>
!> @author Iredell, Ed Hartnett @date 93-11-22
subroutine g2_write_index_headers(lugi, nlen, nnum, idxver, filename)
implicit none
integer, intent(in) :: lugi, nlen, nnum, idxver
character, intent(in) :: filename*(*)
character cd8*8, ct10*10, hostname*15
#ifdef __GFORTRAN__
integer istat
#else
character hostnam*15
integer hostnm
#endif
character chead(2)*81
integer :: kw
! fill first 81-byte header
call date_and_time(cd8, ct10)
chead(1) = '!GFHDR!'
chead(1)(9:10) = ' 1'
chead(1)(12:14) = ' 1'
write(chead(1)(16:20),'(i5)') 162
chead(1)(22:31) = cd8(1:4) // '-' // cd8(5:6) // '-' // cd8(7:8)
chead(1)(33:40) = ct10(1:2) // ':' // ct10(3:4) // ':' // ct10(5:6)
chead(1)(42:47) = 'GB2IX1'
chead(1)(49:54) = ' '
#ifdef __GFORTRAN__
istat = hostnm(hostname)
if (istat .eq. 0) then
chead(1)(56:70) = '0000'
else
chead(1)(56:70) = '0001'
endif
#else
chead(1)(56:70) = hostnam(hostname)
#endif
chead(1)(72:80) = 'grb2index'
chead(1)(81:81) = char(10)
! fill second 81-byte header
if (idxver .eq. 1) then
chead(2) = 'IX1FORM:'
else
chead(2) = 'IX2FORM:'
endif
write(chead(2)(9:38),'(3i10)') 162, nlen, nnum
chead(2)(41:80) = filename
chead(2)(81:81) = char(10)
! write headers at beginning of index file
call bawrite(lugi, 0, 162, kw, chead)
return
end subroutine g2_write_index_headers
!> Find, read or generate a version 1 GRIB2 index for a GRIB2 file
!> (which must be < 2 GB).
!>
!> If the index already exists in library memory, it is returned,
!> otherwise, the index is read from an existing index file associated
!> with unit lugi or generated from the GRIB2 file lugb.
!>
!> Users can force a regeneration of an index: if lugi equals lugb,
!> the index will be regenerated from the data in file lugb. If lugi
!> is less than zero, then the index is re-read from index file
!> abs(lugi).
!>
!> This subroutine allocates memory and stores the resulting pointers
!> in an array that is a Fortran "save" variable. The result is that
!> the memory will not be freed by the library and cannot be reached
!> by the caller. To free this memory call gf_finalize() after all
!> library operations are complete.
!>
!> @note The file unit numbers must be in range 1 - 9999.
!>
!> @param[in] lugb integer unit of the GRIB2 data file. File must
!> have been opened with [baopen() or baopenr()]
!> (https://noaa-emc.github.io/NCEPLIBS-bacio/) before calling this
!> routine. If 0, then all saved memory will be released (necessary
!> for g2_finalize()).
!> @param[in] lugi integer unit of the GRIB2 index file.
!> If nonzero, file must have been opened with [baopen() or baopenr()]
!> (https://noaa-emc.github.io/NCEPLIBS-bacio/) before
!> calling this routine. Set to 0 to get index information from the GRIB2 file.
!> @param[inout] cindex character*1 Pointer to a buffer that will get
!> index records.
!> @param[out] nlen integer Total length of all index records.
!> @param[out] nnum integer Number of index records.
!> @param[out] iret integer Return code:
!> - 0 No error.
!> - 90 Unit number out of range.
!> - 96 Error reading/creating index file.
!> - 97 Index version 2 detected.
!>
!> @author Stephen Gilbert, Ed Hartnett @date 2005-03-15
subroutine getidx(lugb, lugi, cindex, nlen, nnum, iret)
implicit none
integer, intent(in) :: lugb, lugi
character(len = 1), pointer, dimension(:) :: cindex
integer, intent(out) :: nlen, nnum, iret
integer :: idxver
interface
subroutine getidx2(lugb, lugi, idxver, cindex, nlen, nnum, iret)
integer, intent(in) :: lugb, lugi
integer, intent(inout) :: idxver
character(len = 1), pointer, dimension(:) :: cindex
integer, intent(out) :: nlen, nnum, iret
end subroutine getidx2
end interface
! When getidx() is called, always use index version 1. Call
! getidx2() for a chance to set the index version.
idxver = 1
call getidx2(lugb, lugi, idxver, cindex, nlen, nnum, iret)
! If index version 2 is being used, return error.
if (iret .eq. 0 .and. idxver .eq. 2) iret = 97
end subroutine getidx
!> Find, read or generate a version 1 or 2 GRIB2 index for a GRIB2
!> file (which may be > 2 GB).
!>
!> If the index already exists in library memory, it is returned,
!> otherwise, the index is read from an existing index file associated
!> with unit lugi or generated from the GRIB2 file lugb.
!>
!> Users can force a regeneration of an index: if lugi equals lugb,
!> the index will be regenerated from the data in file lugb. If lugi
!> is less than zero, then the index is re-read from index file
!> abs(lugi).
!>
!> This subroutine allocates memory and stores the resulting pointers
!> in an array that is a Fortran "save" variable. The result is that
!> the memory will not be freed by the library and cannot be reached
!> by the caller. To free this memory call gf_finalize() after all
!> library operations are complete.
!>
!> @note The file unit numbers must be in range 1 - 9999.
!>
!> @param[in] lugb integer unit of the GRIB2 data file. File must
!> have been opened with [baopen() or baopenr()]
!> (https://noaa-emc.github.io/NCEPLIBS-bacio/) before calling this
!> routine. If 0, then all saved memory will be released (necessary
!> for g2_finalize()).
!> @param[in] lugi integer unit of the GRIB2 index file.
!> If nonzero, file must have been opened with [baopen() or baopenr()]
!> (https://noaa-emc.github.io/NCEPLIBS-bacio/) before
!> calling this routine. Set to 0 to get index information from the GRIB2 file.
!> @param[in] idxver Index version, 1 for legacy, 2 if files may be > 2 GB.
!> index records.
!> @param[inout] cindex character*1 Pointer to a buffer that will get
!> index records.
!> @param[out] nlen integer Total length of all index records.
!> @param[out] nnum integer Number of index records.
!> @param[out] iret integer Return code:
!> - 0 No error.
!> - 90 Unit number out of range.
!> - 96 Error reading/creating index file.
!>
!> @author Stephen Gilbert, Ed Hartnett @date Feb 9, 2024
subroutine getidx2(lugb, lugi, idxver, cindex, nlen, nnum, iret)
use g2logging
implicit none
integer, intent(in) :: lugb, lugi
integer, intent(inout) :: idxver
character(len = 1), pointer, dimension(:) :: cindex
integer, intent(out) :: nlen, nnum, iret
integer, parameter :: maxidx = 10000
integer (kind = 8), parameter :: msk1 = 32000_8, msk2 = 4000_8
integer :: lux
integer :: irgi, mskp, nmess, i
type gindex
integer :: nlen
integer :: nnum
integer :: idxver
character(len = 1), pointer, dimension(:) :: cbuf
end type gindex
type(gindex), save :: idxlist(10000)
data lux/0/
! declare interfaces (required for cbuf pointer)
interface
subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret)
integer, intent(in) :: lugi
character(len=1), pointer, dimension(:) :: cbuf
integer, intent(out) :: idxver, nlen, nnum, iret
end subroutine getg2i2
subroutine getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, &
nlen, nnum, nmess, iret)
integer, intent(in) :: lugb
integer (kind = 8), intent(in) :: msk1, msk2
integer, intent(in) :: mnum, idxver
character(len = 1), pointer, dimension(:) :: cbuf
integer, intent(out) :: nlen, nnum, nmess, iret
end subroutine getg2i2r
end interface
#ifdef LOGGING
! Log results for debugging.
write(g2_log_msg, '(a, i2, a, i2, a, i1)') 'getidx2: lugb ', lugb, ' lugi ', lugi, &
' idxver ', idxver
call g2_log(1)
#endif
! Free all associated memory and exit.
if (lugb .eq. 0) then
!print *, 'getidx: Freeing all memory'
do i = 1, 10000
if (associated(idxlist(i)%cbuf)) then
!print *, 'deallocating ', loc(idxlist(i)%cbuf)
deallocate(idxlist(i)%cbuf)
nullify(idxlist(i)%cbuf)
endif
end do
iret = 0
return
endif
! Determine whether index buffer needs to be initialized.
lux = 0
iret = 0
if (lugb .le. 0 .or. lugb .gt. 9999) then
print *, ' file unit number out of range'
print *, ' use unit numbers in range: 0 - 9999 '
iret = 90
return
endif
! Force regeneration of index from GRIB2 file.
if (lugi .eq. lugb) then
if (associated(idxlist(lugb)%cbuf)) &
deallocate(idxlist(lugb)%cbuf)
!print *, 'Force regeneration'
nullify(idxlist(lugb)%cbuf)
idxlist(lugb)%nlen = 0
idxlist(lugb)%nnum = 0
lux = 0
endif
! Force re-read of index from indexfile.
if (lugi .lt. 0) then
! associated with unit abs(lugi)
if (associated(idxlist(lugb)%cbuf)) &
deallocate(idxlist(lugb)%cbuf)
!print *, 'Force re-read'
nullify(idxlist(lugb)%cbuf)
idxlist(lugb)%nlen = 0
idxlist(lugb)%nnum = 0
lux = abs(lugi)
endif
! Check if index already exists in memory.
if (associated(idxlist(lugb)%cbuf)) then
!print *, 'Index exists in memory!'
cindex => idxlist(lugb)%cbuf
nlen = idxlist(lugb)%nlen
nnum = idxlist(lugb)%nnum
idxver = idxlist(lugb)%idxver
return
endif
! Either read index from index file, or generate it from the GRIB2
! file. This is an index for all messages in the file, each message
! gets an index record, all stuffed into idxlist(lugbb)%cbuf.
irgi = 0
if (lux .gt. 0) then
call getg2i2(lux, idxlist(lugb)%cbuf, idxver, nlen, nnum, irgi)
elseif (lux .le. 0) then
mskp = 0
call getg2i2r(lugb, msk1, msk2, mskp, idxver, idxlist(lugb)%cbuf, &
nlen, nnum, nmess, irgi)
endif
! Handle errors.
if (irgi .ne. 0) then
nlen = 0
nnum = 0
print *, ' error reading index file '
iret = 96
return
endif
! Fill these values.
cindex => idxlist(lugb)%cbuf
idxlist(lugb)%nlen = nlen
idxlist(lugb)%nnum = nnum
idxlist(lugb)%idxver = idxver
end subroutine getidx2
!> Read a version 1 index file and return its contents.
!>
!> The index file may be generated by the grb2index utility of the
!> [NCEPLIBS-grib_util](https://github.com/NOAA-EMC/NCEPLIBS-grib_util)
!> project.
!>
!> For the contents of the index file, see getgi2().
!>
!> This subroutine is maintained for backward compatibility, and may
!> only be used with version 1 of the index files. To handle both
!> version 1 and 2 index files, use getg2i2().
!>
!> @param[in] lugi Integer unit of the GRIB index file. Must be opened
!> by [baopen() or baopenr()]
!> (https://noaa-emc.github.io/NCEPLIBS-bacio/).
!> @param[out] cbuf Pointer to a buffer that will get the index
!> records. Memory will be allocated within this function, so callers
!> must free the memory that cbuf points to, using deallocate(cbuf)
!> when cbuf is no longer needed.
!> @param[out] nlen Total length of all index records.
!> @param[out] nnum Number of index records.
!> @param[out] iret Return code.
!> - 0 No error.
!> - 2 not enough memory to hold index buffer
!> - 3 error reading index file buffer
!> - 4 error reading index file header
!> - 5 index format 2 detected
!>
!> @author Mark Iredell, Ed Hartnett @date 2000-05-26
subroutine getg2i(lugi, cbuf, nlen, nnum, iret)
implicit none
integer, intent(in) :: lugi
character(len=1), pointer, dimension(:) :: cbuf
integer, intent(out) :: nlen, nnum, iret
integer :: idxver
interface
subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret)
integer, intent(in) :: lugi
character(len=1), pointer, dimension(:) :: cbuf
integer, intent(out) :: idxver, nlen, nnum, iret
end subroutine getg2i2
end interface
! Call the version of this subroutine that handles version 1 and
! version 2. If getg2i() is called on a version 2 index, return
! error.
call getg2i2(lugi, cbuf, idxver, nlen, nnum, iret)
if (idxver .eq. 2) iret = 5
end subroutine getg2i
!> Read a version 1 or 2 index file and return its contents.
!>
!> The index file may be generated by the grb2index utility of the
!> [NCEPLIBS-grib_util](https://github.com/NOAA-EMC/NCEPLIBS-grib_util)
!> project.
!>
!> The index file has two header records:
!> 1. an 81-byte header with 'GB2IX1' in columns 42-47
!> 2. an 81-byte header with the index version number, the number of
!> bytes to skip before index records, total length in bytes of the
!> index records, number of index records, and the GRIB file basename.
!>
!> Each record in the index table contains the following fields. All
!> integers are in big-endian format in the file. The only difference
!> between index version 1 and index version 2 is the size of the
!> field containing the number of bytes to skip in file before
!> message. To accomodate files > 2 GB, this must be a 64-bit int.
!>
!> Index Version 1 | Index Version 2 | Contents
!> ----------------|-----------------|---------
!> 001 - 004 | 001 - 004 | length of index record (4 bytes)
!> 005 - 008 | 005 - 012 | bytes to skip in data file before grib message (4/8 bytes)
!> 009 - 012 | 013 - 020 | bytes to skip in message before lus (local use) set = 0, if no local section. (4/8 bytes)
!> 013 - 016 | 021 - 028 | bytes to skip in message before gds (4/8 bytes)
!> 017 - 020 | 029 - 036 | bytes to skip in message before pds (4/8 bytes)
!> 021 - 024 | 037 - 044 | bytes to skip in message before drs (4/8 bytes)
!> 025 - 028 | 045 - 052 | bytes to skip in message before bms (4/8 bytes)
!> 029 - 032 | 053 - 060 | bytes to skip in message before data section (4/8 bytes)
!> 033 - 040 | 061 - 068 | bytes total in the message (8 bytes)
!> 041 - 041 | 069 - 069 | grib version number (always 2) (1 byte)
!> 042 - 042 | 070 - 070 | message discipline (1 byte)
!> 043 - 044 | 071 - 072 | field number within grib2 message (2 bytes)
!> 045 - ii | 073 - ii | identification section (ids) (character)
!> ii+1- jj | ii+1- jj | grid definition section (gds) (character)
!> jj+1- kk | jj+1- kk | product definition section (pds) (character)
!> kk+1- ll | kk+1- ll | the data representation section (drs) (character)
!> ll+1-ll+6 | ll+1-ll+6 | first 6 bytes of the bit map section (bms) (character)
!>
!> @note Subprogram can be called from a multiprocessing environment.
!> Do not engage the same logical unit from more than one processor.
!>
!> @param[in] lugi Integer unit of the unblocked GRIB index file. Must
!> be opened by [baopen() or baopenr()]
!> (https://noaa-emc.github.io/NCEPLIBS-bacio/).
!> @param[out] cbuf Pointer to a buffer that will get the index
!> records. Memory will be allocated within this function, so callers
!> must free the memory that cbuf points to, using deallocate(cbuf)
!> when cbuf is no longer needed.
!> @param[out] idxver Index version of this index, will be 1 or 2.
!> @param[out] nlen Total length of all index records.
!> @param[out] nnum Number of index records.
!> @param[out] iret Return code.
!> - 0 No error.
!> - 2 not enough memory to hold index buffer
!> - 3 error reading index file buffer
!> - 4 error reading index file header
!>
!> @author Ed Hartnett, Mark Iredell @date Feb 9, 2024
subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret)
use g2logging
implicit none
integer, intent(in) :: lugi
character(len=1), pointer, dimension(:) :: cbuf
integer, intent(out) :: idxver, nlen, nnum, iret
character chead*162
integer :: ios, istat, lbuf, lhead, nskp
#ifdef LOGGING
write(g2_log_msg, *) 'getg2i2: lugi ', lugi
call g2_log(1)
#endif
nullify(cbuf)
nlen = 0
nnum = 0
iret = 4
call baread(lugi, 0, 162, lhead, chead)
#ifdef LOGGING
write(g2_log_msg, *) 'lhead', lhead
call g2_log(2)
#endif
if (lhead .eq. 162 .and. chead(42:47) .eq. 'GB2IX1') then
read(chead(82:162), '(2x, i1, 5x, 3i10, 2x, a40)', iostat = ios) idxver, nskp, nlen, nnum
#ifdef LOGGING
write(g2_log_msg, *) 'ios', ios, 'idxver', idxver, 'nskp', nskp, 'nlen', nlen, 'nnum', nnum
call g2_log(2)
#endif
if (ios .eq. 0) then
allocate(cbuf(nlen), stat = istat) ! Allocate space for cbuf.
if (istat .ne. 0) then
iret = 2
return
endif
iret = 0
call baread(lugi, nskp, nlen, lbuf, cbuf)
if (lbuf .ne. nlen) iret = 3
endif
endif
end subroutine getg2i2
!> Generate a version 1 index record for each message in a GRIB2 file.
!>
!> The index record contains byte offsets to the message, it's length,
!> and byte offsets within the message to each section. The index file
!> record format is documented in subroutine ixgb2().
!>
!> @note Subprogram can be called from a multiprocessing environment.
!> Do not engage the same logical unit from more than one processor.
!>
!> @param[in] lugb Unit of the unblocked GRIB file. Must
!> be opened by [baopen() or baopenr()]
!> (https://noaa-emc.github.io/NCEPLIBS-bacio/).
!> @param[in] msk1 Number of bytes to search for first message.
!> @param[in] msk2 Number of bytes to search for other messages.
!> @param[in] mnum Number of GRIB messages to skip (usually 0).
!> @param[out] cbuf Pointer to a buffer that will get the index
!> records. If any memory is associated with cbuf when this subroutine
!> is called, cbuf will be nullified in the subroutine. Initially cbuf
!> will get an allocation of 5000 bytes. realloc() will be used to
!> increase the size if necessary. Users must free memory that cbuf
!> points to when cbuf is no longer needed.
!> @param[out] nlen Total length of index record buffer in bytes.
!> @param[out] nnum Number of index records, zero if no GRIB
!> messages are found.
!> @param[out] nmess Last GRIB message in file successfully processed
!> @param[out] iret Return code.
!> - 0 No error.
!> - 1 Not enough memory available to hold full index buffer.
!> - 2 Not enough memory to allocate initial index buffer.
!> - 3 Error deallocating memory.
!>
!> @author Mark Iredell, Ed Hartnett @date 1995-10-31
subroutine getg2ir(lugb, msk1, msk2, mnum, cbuf, nlen, nnum, nmess, iret)
use re_alloc ! needed for subroutine realloc
implicit none
integer, intent(in) :: lugb
integer, intent(in) :: msk1, msk2
integer, intent(in) :: mnum
character(len = 1), pointer, dimension(:) :: cbuf
integer, intent(out) :: nlen, nnum, nmess, iret
integer (kind = 8) :: msk1_8, msk2_8
interface
subroutine getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, &
nlen, nnum, nmess, iret)
integer, intent(in) :: lugb
integer (kind = 8), intent(in) :: msk1, msk2
integer, intent(in) :: mnum, idxver
character(len = 1), pointer, dimension(:) :: cbuf
integer, intent(out) :: nlen, nnum, nmess, iret
end subroutine getg2i2r
end interface
msk1_8 = msk1
msk2_8 = msk2
call getg2i2r(lugb, msk1_8, msk2_8, mnum, 1, cbuf, nlen, nnum, nmess, iret)
end subroutine getg2ir
!> Generate a version 1 or 2 index record for each message in a GRIB2
!> file.
!>
!> The index record contains byte offsets to the message, it's length,
!> and byte offsets within the message to each section. The index file
!> record format is documented in subroutine ixgb2().
!>
!> @note Subprogram can be called from a multiprocessing environment.
!> Do not engage the same logical unit from more than one processor.
!>
!> @param[in] lugb Unit of the GRIB file. Must be opened by [baopen()
!> or baopenr()] (https://noaa-emc.github.io/NCEPLIBS-bacio/).
!> @param[in] msk1 Number of bytes to search for first message.
!> @param[in] msk2 Number of bytes to search for other messages.
!> @param[in] mnum Number of GRIB messages to skip (usually 0).
!> @param[in] idxver Index version number. Use 1 for legacy, 2 for files > 2 GB.
!> @param[out] cbuf Pointer to a buffer that will get the index
!> records. If any memory is associated with cbuf when this subroutine
!> is called, cbuf will be nullified in the subroutine. Initially cbuf
!> will get an allocation of 5000 bytes. realloc() will be used to
!> increase the size if necessary. Users must free memory that cbuf
!> points to when cbuf is no longer needed.
!> @param[out] nlen Total length of index record buffer in bytes.
!> @param[out] nnum Number of index records, zero if no GRIB
!> messages are found.
!> @param[out] nmess Last GRIB message in file successfully processed
!> @param[out] iret Return code.
!> - 0 No error.
!> - 1 Not enough memory available to hold full index buffer.
!> - 2 Not enough memory to allocate initial index buffer.
!> - 3 Error deallocating memory.
!>
!> @author Mark Iredell, Ed Hartnett @date 1995-10-31
subroutine getg2i2r(lugb, msk1, msk2, mnum, idxver, cbuf, nlen, nnum, nmess, iret)
use g2logging
use re_alloc ! needed for subroutine realloc
implicit none
integer, intent(in) :: lugb
integer (kind = 8), intent(in) :: msk1, msk2
integer, intent(in) :: mnum, idxver
character(len = 1), pointer, dimension(:) :: cbuf
integer, intent(out) :: nlen, nnum, nmess, iret
character(len = 1), pointer, dimension(:) :: cbuftmp
integer :: nbytes, newsize, next, numfld, m, mbuf
integer (kind = 8) :: iseek, lskip, lgrib
integer :: istat, init, iret1, lgrib4
parameter(init = 50000, next = 10000)
interface ! required for cbuf pointer
subroutine ix2gb2(lugb, lskip8, idxver, lgrib8, cbuf, numfld, mlen, iret)
integer, intent(in) :: lugb
integer (kind = 8), intent(in) :: lskip8
integer, intent(in) :: idxver
integer (kind = 8), intent(in) :: lgrib8
character(len = 1), pointer, intent(inout), dimension(:) :: cbuf
integer, intent(out) :: numfld, mlen, iret
end subroutine ix2gb2
end interface
#ifdef LOGGING
! Log results for debugging.
write(g2_log_msg, *) 'getg2i2r: lugb ', lugb, ' msk1 ', msk1, ' msk2 ', msk2, 'idxver', idxver
call g2_log(1)
#endif
! Initialize.
iret = 0
nullify(cbuf)
mbuf = init
allocate(cbuf(mbuf), stat = istat) ! allocate initial space for cbuf.
if (istat .ne. 0) then
iret = 2
return
endif
! Search for first grib message.
iseek = 0_8
call skgb8(lugb, iseek, msk1, lskip, lgrib)
do m = 1, mnum
if (lgrib .gt. 0) then
iseek = lskip + lgrib
call skgb8(lugb, iseek, msk2, lskip, lgrib)
endif
enddo
! Get index records for every grib message found.
nlen = 0
nnum = 0
nmess = mnum
do while (iret .eq. 0 .and. lgrib .gt. 0)
lgrib4 = int(lgrib, kind(4))
call ix2gb2(lugb, lskip, idxver, lgrib, cbuftmp, numfld, nbytes, iret1)
if (iret1 .ne. 0) print *, ' SAGT ', numfld, nbytes, iret1
if (nbytes + nlen .gt. mbuf) then ! Allocate more space, if necessary.
newsize = max(mbuf + next, mbuf + nbytes)
call realloc(cbuf, nlen, newsize, istat)
if (istat .ne. 0) then
iret = 1
return
endif
mbuf = newsize
endif
! If index records were returned in cbuftmp from ixgb2,
! copy cbuftmp into cbuf, then deallocate cbuftmp when done.
if (associated(cbuftmp)) then
cbuf(nlen + 1 : nlen + nbytes) = cbuftmp(1 : nbytes)
deallocate(cbuftmp, stat = istat)
if (istat .ne. 0) then
print *, ' deallocating cbuftmp ... ', istat
iret = 3
return
endif
nullify(cbuftmp)
nnum = nnum + numfld
nlen = nlen + nbytes
nmess = nmess + 1
endif
! Look for next grib message.
iseek = lskip + lgrib
call skgb8(lugb, iseek, msk2, lskip, lgrib)
enddo
end subroutine getg2i2r
!> Find information about a GRIB field from the index and fill a @ref
!> grib_mod::gribfield.
!>
!> For a description of the index record see getg2i2().
!>
!> Users of this routine will need to include the line "use grib_mod"
!> in their calling routine.
!>
!> The unpacked bitmap and bitmap data field are the only components
!> of the @ref grib_mod::gribfield type not set by this routine.
!>
!> @note This subprogram is intended for private use by getgb2()
!> routines only. Note that derived type @ref grib_mod::gribfield contains
!> pointers to many arrays of data. The memory for these arrays is
!> allocated when the values in the arrays are set. Users must free this
!> memory, when it is no longer needed, by a call to subroutine
!> gf_free().
!>
!> @param[in] cbuf Buffer (of size nlen bytes) containing index data.
!> @param[in] nlen Total length of all index records.
!> @param[in] nnum Number of index records.
!> @param[in] j Number of fields to skip (0 to search from beginning).
!> @param[in] jdisc GRIB2 discipline number of requested field. See
!> [GRIB2 - TABLE 0.0 -
!> DISCIPLINE](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table0-0.shtml).
!> Use -1 to accept any discipline.
!> @param[in] jids Array of values in the identification
!> section. (Set to -9999 for wildcard.)
!> - jids(1) Identification of originating centre. See [TABLE 0 -
!> NATIONAL/INTERNATIONAL ORIGINATING
!> CENTERS](https://www.nco.ncep.noaa.gov/pmb/docs/on388/table0.html).
!> - jids(2) Identification of originating sub-centre. See [TABLE C -
!> NATIONAL
!> SUB-CENTERS](https://www.nco.ncep.noaa.gov/pmb/docs/on388/tablec.html).
!> - jids(3) GRIB master tables version number. See [GRIB2 - TABLE 1.0
!> - GRIB Master Tables Version
!> Number](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-0.shtml).
!> - jids(4) GRIB local tables version number. See [GRIB2 - TABLE 1.1
!> - GRIB Local Tables Version
!> Number](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-1.shtml).
!> - jids(5) Significance of reference time. See [GRIB2 - TABLE 1.2 -
!> Significance of Reference
!> Time](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-2.shtml).
!> - jids(6) year (4 digits)
!> - jids(7) month
!> - jids(8) day
!> - jids(9) hour
!> - jids(10) minute
!> - jids(11) second
!> - jids(12) Production status of processed data. See [GRIB2 - TABLE
!> 1.3 - Production Status of
!> Data](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-3.shtml).
!> - jids(13) Type of processed data. See [GRIB2 - TABLE 1.4 - TYPE OF
!> DATA](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-4.shtml).
!> @param[in] jpdtn Product Definition Template (PDT) number (n)
!> (if = -1, don't bother matching PDT - accept any).
!> @param[in] jpdt Array of values defining the Product Definition
!> Template of the field for which to search (=-9999 for wildcard).
!> @param[in] jgdtn Grid Definition Template (GDT) number (if = -1,
!> don't bother matching GDT - accept any).
!> @param[in] jgdt array of values defining the Grid Definition
!> Template of the field for which to search (=-9999 for wildcard).
!> @param[out] k Field number unpacked.
!> @param[out] gfld Derived type @ref grib_mod::gribfield.
!> @param[out] lpos Starting position of the found index record
!> within the complete index buffer, CBUF. = 0, if request not found.
!> @param[out] iret integer return code:
!> - 0 No error.
!> - 97 Error reading GRIB file.
!> - other gf_getfld GRIB2 unpacker return code.
!>
!> @author Stephen Gilbert, Ed Hartnett @date 2002-01-15
subroutine getgb2s(cbuf, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, &
jgdt, k, gfld, lpos, iret)
use grib_mod
implicit none
character(len = 1), intent(in) :: cbuf(nlen)
integer, intent(in) :: nlen, nnum, j, jdisc
integer, dimension(:) :: jids(*)
integer, intent(in) :: jpdtn
integer, dimension(:) :: jpdt(*)
integer, intent(in) :: jgdtn
integer, dimension(:) :: jgdt(*)
integer, intent(out) :: k
type(gribfield), intent(out) :: gfld
integer, intent(out) :: lpos, iret
interface
subroutine getgb2s2(cbuf, idxver, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, &
jgdt, k, gfld, lpos, iret)
import gribfield
character(len = 1), intent(in) :: cbuf(nlen)
integer, intent(in) :: idxver, nlen, nnum, j, jdisc
integer, dimension(:) :: jids(*)
integer, intent(in) :: jpdtn
integer, dimension(:) :: jpdt(*)
integer, intent(in) :: jgdtn
integer, dimension(:) :: jgdt(*)
integer, intent(out) :: k
type(gribfield), intent(out) :: gfld
integer, intent(out) :: lpos, iret
end subroutine getgb2s2
end interface
! When getgb2s() is called, always use index version 1. Call
! getgb2s2() to handle version 1 or 2.
call getgb2s2(cbuf, 1, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, &
jgdt, k, gfld, lpos, iret)
end subroutine getgb2s
!> Find information about a GRIB field from the index and fill a @ref
!> grib_mod::gribfield.
!>
!> For a description of the index record see getg2i2().
!>
!> Users of this routine will need to include the line "use grib_mod"
!> in their calling routine.
!>
!> The unpacked bitmap and bitmap data field are the only components
!> of the @ref grib_mod::gribfield type not set by this routine.
!>
!> @note This subprogram is intended for private use by getgb2()
!> routines only. Note that derived type @ref grib_mod::gribfield contains
!> pointers to many arrays of data. The memory for these arrays is
!> allocated when the values in the arrays are set. Users must free this
!> memory, when it is no longer needed, by a call to subroutine
!> gf_free().
!>
!> @param[in] cbuf Buffer (of size nlen bytes) containing index data.
!> @param[in] idxver Index version, 1 for legacy, 2 if files may be > 2 GB.
!> @param[in] nlen Total length of all index records.
!> @param[in] nnum Number of index records.
!> @param[in] j Number of fields to skip (0 to search from beginning).
!> @param[in] jdisc GRIB2 discipline number of requested field. See
!> [GRIB2 - TABLE 0.0 -
!> DISCIPLINE](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table0-0.shtml).
!> Use -1 to accept any discipline.
!> @param[in] jids Array of values in the identification
!> section. (Set to -9999 for wildcard.)
!> - jids(1) Identification of originating centre. See [TABLE 0 -
!> NATIONAL/INTERNATIONAL ORIGINATING
!> CENTERS](https://www.nco.ncep.noaa.gov/pmb/docs/on388/table0.html).
!> - jids(2) Identification of originating sub-centre. See [TABLE C -
!> NATIONAL
!> SUB-CENTERS](https://www.nco.ncep.noaa.gov/pmb/docs/on388/tablec.html).
!> - jids(3) GRIB master tables version number. See [GRIB2 - TABLE 1.0
!> - GRIB Master Tables Version
!> Number](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-0.shtml).
!> - jids(4) GRIB local tables version number. See [GRIB2 - TABLE 1.1
!> - GRIB Local Tables Version
!> Number](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-1.shtml).
!> - jids(5) Significance of reference time. See [GRIB2 - TABLE 1.2 -
!> Significance of Reference
!> Time](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-2.shtml).
!> - jids(6) year (4 digits)
!> - jids(7) month
!> - jids(8) day
!> - jids(9) hour
!> - jids(10) minute
!> - jids(11) second
!> - jids(12) Production status of processed data. See [GRIB2 - TABLE
!> 1.3 - Production Status of
!> Data](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-3.shtml).
!> - jids(13) Type of processed data. See [GRIB2 - TABLE 1.4 - TYPE OF
!> DATA](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table1-4.shtml).
!> @param[in] jpdtn Product Definition Template (PDT) number (n)
!> (if = -1, don't bother matching PDT - accept any).
!> @param[in] jpdt Array of values defining the Product Definition
!> Template of the field for which to search (=-9999 for wildcard).
!> @param[in] jgdtn Grid Definition Template (GDT) number (if = -1,
!> don't bother matching GDT - accept any).
!> @param[in] jgdt array of values defining the Grid Definition
!> Template of the field for which to search (=-9999 for wildcard).
!> @param[out] k Field number unpacked.
!> @param[out] gfld Derived type @ref grib_mod::gribfield.
!> @param[out] lpos Starting position of the found index record
!> within the complete index buffer, CBUF. = 0, if request not found.
!> @param[out] iret integer return code:
!> - 0 No error.
!> - 97 Error reading GRIB file.
!> - other gf_getfld GRIB2 unpacker return code.
!>
!> @author Stephen Gilbert, Ed Hartnett @date Feb 9 2024
subroutine getgb2s2(cbuf, idxver, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, &
jgdt, k, gfld, lpos, iret)
use grib_mod
use g2logging
implicit none
character(len = 1), intent(in) :: cbuf(nlen)
integer, intent(in) :: idxver, nlen, nnum, j, jdisc
integer, dimension(:) :: jids(*)
integer, intent(in) :: jpdtn
integer, dimension(:) :: jpdt(*)
integer, intent(in) :: jgdtn
integer, dimension(:) :: jgdt(*)
integer, intent(out) :: k
type(gribfield), intent(out) :: gfld
integer, intent(out) :: lpos, iret
integer :: kgds(5)
logical :: match1, match3, match4
integer :: i, icnd, inlen, iof, ipos, jpos, lsec1, lsec3, lsec4, lsec5, numgdt, numpdt, inc
interface
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 gf_unpack1(cgrib, lcgrib, iofst, ids, idslen, ierr)
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib
integer, intent(inout) :: iofst
integer, pointer, dimension(:) :: ids
integer, intent(out) :: ierr, idslen
end subroutine gf_unpack1
subroutine gf_unpack3(cgrib, lcgrib, iofst, igds, igdstmpl, &
mapgridlen, ideflist, idefnum, ierr)
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib
integer, intent(inout) :: iofst
integer, pointer, dimension(:) :: igdstmpl, ideflist
integer, intent(out) :: igds(5)
integer, intent(out) :: ierr, idefnum
end subroutine gf_unpack3
subroutine gf_unpack4(cgrib, lcgrib, iofst, ipdsnum, ipdstmpl, &
mappdslen, coordlist, numcoord, ierr)
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib
integer, intent(inout) :: iofst
real, pointer, dimension(:) :: coordlist
integer, pointer, dimension(:) :: ipdstmpl
integer, intent(out) :: ipdsnum
integer, intent(out) :: ierr, numcoord
end subroutine gf_unpack4
subroutine gf_unpack5(cgrib, lcgrib, iofst, ndpts, idrsnum, &
idrstmpl, mapdrslen, ierr)
character(len = 1), intent(in) :: cgrib(lcgrib)
integer, intent(in) :: lcgrib
integer, intent(inout) :: iofst
integer, intent(out) :: ndpts, idrsnum
integer, pointer, dimension(:) :: idrstmpl
integer, intent(out) :: ierr
end subroutine gf_unpack5
end interface
#ifdef LOGGING
! Log results for debugging.
write(g2_log_msg, *) 'getgb2s2: idxver ', idxver, ' nlen ', nlen, &
' nnum ', nnum, ' j ', j, ' jdisc ', jdisc
call g2_log(1)
#endif