From dcd372d4814124b64132b53b4826b1052e0acda6 Mon Sep 17 00:00:00 2001 From: Edward Hartnett <38856240+edwardhartnett@users.noreply.github.com> Date: Wed, 23 Aug 2023 08:40:04 +0200 Subject: [PATCH] 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 --- src/skgb.F90 | 68 ++++++++++++++++++++++++++++++++++----------- tests/test_skgb.F90 | 28 ++++++++++++++++--- 2 files changed, 76 insertions(+), 20 deletions(-) diff --git a/src/skgb.F90 b/src/skgb.F90 index 8a3f1b6f..c1433821 100644 --- 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 index 2dfc765b..e816e7fb 100644 --- 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