adding skgb8() to handle GRIB2 files > 2 GB (#493)

* adding skgb8() * new attempt * new attempt * better attempt * converted to bareadl() * fixed memory problem * more progress * eliminating unneeded vars * more cleanup * more cleanup * more cleanup * more cleanup * more cleanup * more cleanup * more cleanup * more cleanup * more cleanup * removed print statement
NOAA-EMC · Aug 23, 2023 · dcd372d · dcd372d
1 parent 2dc352f
commit dcd372d
Show file tree

Hide file tree

Showing 2 changed files with 76 additions and 20 deletions.
diff --git a/src/skgb.F90 b/src/skgb.F90
@@ -20,43 +20,79 @@
 subroutine skgb(lugb, iseek, mseek, lskip, lgrib)
   implicit none
 
-  integer lseek, lugb, iseek, mseek, lskip, lgrib, i1, i4, k, k4, kg, km, ks, kz, kn
+  integer :: lugb, iseek, mseek, lskip, lgrib
+  integer*8 :: iseek8, mseek8, lskip8, lgrib8
+
+  iseek8 = iseek
+  mseek8 = mseek
+  call skgb8(lugb, iseek8, mseek8, lskip8, lgrib8)
+  lskip = int(lskip8, kind = 4)
+  lgrib = int(lgrib8, kind = 4)
+
+end subroutine skgb
+
+!> Search a file for the next GRIB1 or GRIB2 message.
+!>
+!> This function works for GRIB2 files larger than 2 GB, but still
+!> requires that each GRIB2 message be less than 2 GB.
+!>
+!> A grib message is identified by its indicator section, which is an
+!> 8-byte sequence with 'GRIB' in bytes 1-4 and a '1' or '2' in byte
+!> 8. If found, the length of the message is decoded from bytes
+!> 5-7. The search is done over a given section of the file. The
+!> search is terminated if an eof or i/o error is encountered.
+!>
+!> @param[in] lugb Unit of the unblocked grib file.
+!> @param[in] iseek8 Number of bytes to skip before search.
+!> @param[in] mseek8 Maximum number of bytes to search.
+!> @param[out] lskip8 Number of bytes to skip before message.
+!> @param[out] lgrib8 Number of bytes in message (0 if not found).
+!>
+!> @author Edward Hartnett @date 2023-08-21
+subroutine skgb8(lugb, iseek8, mseek8, lskip8, lgrib8)
+  implicit none
+
+  integer*8 iseek8, mseek8, lskip8, lgrib8
+  integer*8 ks8, kn8, kz8, k48, km8
+  integer lseek, lugb, i1, i4, k, kg
   parameter(lseek = 512)
   character z(lseek)
   character z4(4)
 
-  lgrib = 0
-  ks = iseek
-  kn = min(lseek, mseek)
-  kz = lseek
+  lgrib8 = 0
+  ks8 = iseek8
+  kn8 = min(int(lseek, kind = 8), mseek8)
+  kz8 = lseek
 
   !  loop until grib message is found
-  do while (lgrib .eq. 0 .and. kn .ge. 8 .and. kz .eq. lseek)
+  do while (lgrib8 .eq. 0 .and. kn8 .ge. 8 .and. kz8 .eq. lseek)
      !  read partial section
-     call baread(lugb, ks, kn, kz, z)
-     km = kz - 8 + 1
+     call bareadl(lugb, ks8, kn8, kz8, z)
+
+     km8 = kz8 - 8 + 1
      k = 0
      !  look for 'grib...1' in partial section
-     do while (lgrib .eq. 0 .and. k .lt. km)
+     do while (lgrib8 .eq. 0 .and. k .lt. km8)
         call g2_gbytec(z, i4, (k + 0) * 8, 4 * 8)
         call g2_gbytec(z, i1, (k + 7) * 8, 1 * 8)
         if (i4 .eq. 1196575042 .and. (i1 .eq. 1 .or. i1 .eq. 2)) then
            !  look for '7777' at end of grib message
            if (i1 .eq. 1) call g2_gbytec(z, kg, (k + 4) * 8, 3 * 8)
            if (i1 .eq. 2) call g2_gbytec(z, kg, (k + 12) * 8, 4 * 8)
-           call baread(lugb, ks + k + kg-4, 4, k4, z4)
-           if (k4 .eq. 4) then
+
+           call bareadl(lugb, ks8 + k + kg - 4, 4_8, k48, z4)
+           if (k48 .eq. 4) then
               call g2_gbytec(z4, i4, 0, 4 * 8)
               if (i4 .eq. 926365495) then
                  !  grib message found
-                 lskip = ks + k
-                 lgrib = kg
+                 lskip8 = ks8 + k
+                 lgrib8 = kg
               endif
            endif
         endif
         k = k + 1
      enddo
-     ks = ks + km
-     kn = min(lseek, iseek + mseek - ks)
+     ks8 = ks8 + km8
+     kn8 = min(int(lseek, kind = 8), iseek8 + mseek8 - ks8)
   enddo
-end subroutine skgb
+end subroutine skgb8
diff --git a/tests/test_skgb.F90 b/tests/test_skgb.F90
@@ -17,6 +17,7 @@ program test_skgb
 
   integer :: lugb = 5
   integer :: lskip, lgrib
+  integer*8 :: lskip8, lgrib8
   integer :: NUM_MSG
   parameter (NUM_MSG = 26)
   integer :: expected_lskip(NUM_MSG) = (/ &
@@ -26,20 +27,21 @@ program test_skgb
        11183, 18917, 12644, 15067, 15042, 10362, 11326, 21168, 15123, 19821, 22210, 15849, 12631, 21895, &
        12785, 15053, 11342, 12957, 19567, 15006, 14962, 12467, 19677, 10438, 22082, 12497 /)
   integer :: i, start
+  integer*8 :: start8
   integer :: iret
 
   print *, 'Testing skgb()...'
 
   ! Open a real GRIB2 file.
   call baopenr(lugb, TEST_FILE_WW3_WEST, iret)
-  if (iret .ne. 0) stop 100
+  if (iret .ne. 0) stop 1
 
   ! Loop through the file, checking location of each message.
   start = 0
   do i = 1, NUM_MSG
      call skgb(lugb, start, 10000, lskip, lgrib)
-     print *, i, lskip, lgrib
-     if (lskip .ne. expected_lskip(i) .or. lgrib .ne. expected_lgrib(i)) stop 101
+     print *, 'i', i, ' lskip ', lskip, ' lgrib ',lgrib
+     if (lskip .ne. expected_lskip(i) .or. lgrib .ne. expected_lgrib(i)) stop 10
      start = start + lgrib
   end do
 
@@ -50,8 +52,26 @@ program test_skgb
 
   ! Close the file.
   call baclose(lugb, iret)
-  if (iret .ne. 0) stop 199
+  if (iret .ne. 0) stop 20
+
+  print *, 'OK!'
+  print *, 'Testing skgb8()...'
+  ! Open a real GRIB2 file.
+  call baopenr(lugb, TEST_FILE_WW3_WEST, iret)
+  if (iret .ne. 0) stop 100
 
+  ! Loop through the file, checking location of each message.
+  start8 = 0
+  do i = 1, NUM_MSG
+     call skgb8(lugb, start8, 10000_8, lskip8, lgrib8)
+     print *, 'i', i, ' lskip8 ', lskip8, ' lgrib8 ',lgrib8
+     if (lskip8 .ne. expected_lskip(i) .or. lgrib8 .ne. expected_lgrib(i)) stop 110
+     start8 = start8 + lgrib8
+  end do
+
+  ! Close the file.
+  call baclose(lugb, iret)
+  if (iret .ne. 0) stop 120
   print *, 'SUCCESS!...'
 
 end program test_skgb