From 7ceadb8259b99d21fc44933db5479e0c63c78675 Mon Sep 17 00:00:00 2001 From: Edward Hartnett <38856240+edwardhartnett@users.noreply.github.com> Date: Thu, 23 May 2024 07:30:50 +0200 Subject: [PATCH] introduce getgb2p2() (#697) * cleanup * more cleanup * more cleanup * more cleanup * more cleanup * more cleanup * more progress * more progress * adding test for getgb2rp2() * more testing * commenting out test * added test * fixed test interface * more test progress * test working! --- src/g2getgb2.F90 | 322 ++++++++++++++++++++++++---------------- src/g2index.F90 | 24 +-- tests/CMakeLists.txt | 1 + tests/test_getgb2p2.F90 | 104 +++++++++++++ 4 files changed, 316 insertions(+), 135 deletions(-) create mode 100644 tests/test_getgb2p2.F90 diff --git a/src/g2getgb2.F90 b/src/g2getgb2.F90 index 4c1869a8..5af377a9 100644 --- a/src/g2getgb2.F90 +++ b/src/g2getgb2.F90 @@ -455,6 +455,113 @@ end subroutine gf_unpack2 endif end subroutine getgb2l2 +!> Legacy subroutine to find and extract a GRIB2 message from a +!> file. Use getgb2p2() for new code. +!> +!> This subroutine reads a GRIB index file (or optionally the GRIB +!> file itself) to get the index buffer (i.e. table of contents) for +!> the GRIB file. It finds in the index buffer a reference to the +!> GRIB field requested. +!> +!> The GRIB field request specifies the number of fields to skip and +!> the unpacked identification section, grid definition template and +!> product defintion section parameters. (A requested parameter of +!> -9999 means to allow any value of this parameter to be found.) +!> +!> If the requested GRIB field is found, then it is read from the GRIB +!> file and unpacked. If the GRIB field is not found, then the return +!> code will be nonzero. +!> +!> The derived type @ref grib_mod::gribfield contains allocated memory +!> that must be freed by the caller with subroutine gf_free(). +!> +!> @note Specifing an index file may increase speed. +!> Do not engage the same logical unit from more than one processor. +!> +!> @param[in] lugb Unit of the unblocked GRIB data file. The +!> file must have been opened with [baopen() or baopenr()] +!> (https://noaa-emc.github.io/NCEPLIBS-bacio/) before calling this +!> routine. +!> @param[in] lugi Unit of the unblocked GRIB index file. If +!> nonzero, file must have been opened with [baopen() or baopenr()] +!> (https://noaa-emc.github.io/NCEPLIBS-bacio/) before calling this +!> subroutine. Set to 0 to get index buffer from the GRIB file. +!> @param[in] j Number of fields to skip (set to 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.) +!> @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[in] extract value indicating whether to return a +!> GRIB2 message with just the requested field, or the entire +!> GRIB2 message containing the requested field. +!> - .true. return GRIB2 message containing only the requested field. +!> - .false. return entire GRIB2 message containing the requested field. +!> @param[out] k field number unpacked. +!> @param[out] gribm returned GRIB message. +!> @param[out] leng length of returned GRIB message in bytes. +!> @param[out] iret integer return code +!> - 0 No error. +!> - 96 Error reading index. +!> - 97 Error reading GRIB file. +!> - 99 Request not found. +!> +!> @author Mark Iredell @date 1994-04-01 +subroutine getgb2p(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, k, gribm, leng, iret) + use grib_mod + implicit none + + integer, intent(in) :: lugb, lugi, j, jdisc + integer, dimension(:) :: jids(*) + integer, intent(in) :: jpdtn + integer, dimension(:) :: jpdt(*) + integer, intent(in) :: jgdtn + integer, dimension(:) :: jgdt(*) + logical, intent(in) :: extract + integer, intent(out) :: k + character(len = 1), pointer, dimension(:) :: gribm + integer, intent(out) :: leng, iret + integer :: idxver + + integer (kind = 8) :: leng8 + + interface + subroutine getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng8, iret) + integer, intent(in) :: lugb, lugi, j, jdisc + integer, dimension(:) :: jids(*) + integer, intent(in) :: jpdtn + integer, dimension(:) :: jpdt(*) + integer, intent(in) :: jgdtn + integer, dimension(:) :: jgdt(*) + logical, intent(in) :: extract + integer, intent(inout) :: idxver + integer, intent(out) :: k + character(len = 1), pointer, dimension(:) :: gribm + integer (kind = 8), intent(out) :: leng8 + integer, intent(out) :: iret + end subroutine getgb2p2 + end interface + + ! Call the new version of this subroutine, which handles messages > 2 GB. + idxver = 1 + call getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng8, iret) + leng = int(leng8, kind(4)) + +end subroutine getgb2p + !> Find and extract a GRIB2 message from a file. !> !> This subroutine reads a GRIB index file (or optionally the GRIB @@ -532,34 +639,42 @@ end subroutine getgb2l2 !> GRIB2 message containing the requested field. !> - .true. return GRIB2 message containing only the requested field. !> - .false. return entire GRIB2 message containing the requested field. +!> @param[in] idxver The index version, use 2 for new code, 1 for +!> legacy index files. !> @param[out] k field number unpacked. !> @param[out] gribm returned GRIB message. -!> @param[out] leng length of returned GRIB message in bytes. +!> @param[out] leng8 length of returned GRIB message in bytes. !> @param[out] iret integer return code !> - 0 No error. !> - 96 Error reading index. !> - 97 Error reading GRIB file. !> - 99 Request not found. !> -!> @author Mark Iredell @date 1994-04-01 -subroutine getgb2p(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & - extract, k, gribm, leng, iret) +!> @author Alex Richert, Edward Hartnett @date 2024-05-21 +subroutine getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng8, iret) use grib_mod implicit none - integer, intent(in) :: lugb, lugi, j, jdisc, jpdtn, jgdtn - integer, dimension(:) :: jids(*), jpdt(*), jgdt(*) + integer, intent(in) :: lugb, lugi, j, jdisc + integer, dimension(:) :: jids(*) + integer, intent(in) :: jpdtn + integer, dimension(:) :: jpdt(*) + integer, intent(in) :: jgdtn + integer, dimension(:) :: jgdt(*) logical, intent(in) :: extract - integer, intent(out) :: k, iret, leng + integer, intent(inout) :: idxver + integer, intent(out) :: k character(len = 1), pointer, dimension(:) :: gribm + integer (kind = 8), intent(out) :: leng8 + integer, intent(out) :: iret type(gribfield) :: gfld - integer :: msk1, irgi, irgs, jk, lpos, msk2, mskp, nlen, nmess, nnum - + integer :: irgi, irgs, jk, lpos, mskp, nlen, nmess, nnum character(len = 1), pointer, dimension(:) :: cbuf + integer (kind = 8) :: msk1, msk2 parameter(msk1 = 32000, msk2 = 4000) - ! Declare interfaces (required for cbuf pointer). interface subroutine getg2i(lugi, cbuf, nlen, nnum, iret) character(len = 1), pointer, dimension(:) :: cbuf @@ -572,22 +687,53 @@ subroutine getg2ir(lugb, msk1, msk2, mnum, cbuf, nlen, nnum, & integer, intent(in) :: lugb, msk1, msk2, mnum integer, intent(out) :: nlen, nnum, nmess, iret end subroutine getg2ir - subroutine getgb2rp(lugb, cindex, extract, gribm, leng, iret) + 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 + subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng8, iret) integer, intent(in) :: lugb + integer, intent(inout) :: idxver character(len = 1), intent(in) :: cindex(*) logical, intent(in) :: extract - integer, intent(out) :: leng, iret character(len = 1), pointer, dimension(:) :: gribm - end subroutine getgb2rp + integer (kind = 8), intent(out) :: leng8 + integer, intent(out) :: iret + end subroutine getgb2rp2 + 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 + nullify(gribm) + ! Initialize the index information in cbuf. irgi = 0 if (lugi .gt. 0) then - call getg2i(lugi, cbuf, nlen, nnum, irgi) + call getg2i2(lugi, cbuf, idxver, nlen, nnum, irgi) elseif (lugi .le. 0) then mskp = 0 - call getg2ir(lugb, msk1, msk2, mskp, cbuf, nlen, nnum, nmess, irgi) + call getg2i2r(lugb, msk1, msk2, mskp, idxver, cbuf, nlen, nnum, nmess, irgi) endif if (irgi .gt. 1) then iret = 96 @@ -595,7 +741,7 @@ end subroutine getgb2rp endif ! Find info from index and fill a grib_mod::gribfield variable. - call getgb2s(cbuf, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + call getgb2s2(cbuf, idxver, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & jk, gfld, lpos, irgs) if (irgs .ne. 0) then iret = 99 @@ -604,8 +750,7 @@ end subroutine getgb2rp endif ! Extract grib message from file. - nullify(gribm) - call getgb2rp(lugb, cbuf(lpos:), extract, gribm, leng, iret) + call getgb2rp2(lugb, idxver, cbuf(lpos:), extract, gribm, leng8, iret) k = jk @@ -613,98 +758,7 @@ end subroutine getgb2rp if (associated(cbuf)) deallocate(cbuf) call gf_free(gfld) -end subroutine getgb2p - -! subroutine getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & -! extract, k, gribm, leng, iret) -! use grib_mod -! implicit none - -! integer, intent(in) :: lugb, lugi, j, jdisc, jpdtn, jgdtn -! integer, dimension(:) :: jids(*), jpdt(*), jgdt(*) -! logical, intent(in) :: extract -! integer, intent(out) :: k, iret, leng -! character(len = 1), pointer, dimension(:) :: gribm - -! type(gribfield) :: gfld -! integer :: irgi, irgs, jk, lpos, mskp, nlen, nmess, nnum, idxver -! integer(kind = 8) :: msk1, msk2 - -! character(len = 1), pointer, dimension(:) :: cbuf -! parameter(msk1 = 32000, msk2 = 4000) - -! ! 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 -! 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 -! subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng, iret) -! integer, intent(in) :: lugb, idxver -! character(len = 1), intent(in) :: cindex(*) -! logical, intent(in) :: extract -! character(len = 1), pointer, dimension(:) :: gribm -! integer, intent(out) :: leng, iret -! end subroutine getgb2rp2 -! end interface - -! ! Initialize the index information in cbuf. -! irgi = 0 -! if (lugi .gt. 0) then -! call getg2i2(lugi, cbuf, idxver, nlen, nnum, irgi) -! elseif (lugi .le. 0) then -! mskp = 0 -! call getg2i2r(lugb, msk1, msk2, mskp, 2, cbuf, nlen, nnum, nmess, irgi) -! endif -! if (irgi .gt. 1) then -! iret = 96 -! return -! endif - -! ! Find info from index and fill a grib_mod::gribfield variable. -! call getgb2s2(cbuf, 2, nlen, nnum, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & -! jk, gfld, lpos, irgs) - -! if (irgs .ne. 0) then -! iret = 99 -! call gf_free(gfld) -! return -! endif - -! ! Extract grib message from file. -! nullify(gribm) -! call getgb2rp2(lugb, 2, cbuf(lpos:), .false., gribm, leng, iret) - -! k = jk - -! ! Free cbuf memory allocated in getg2i/getg2ir(). -! if (associated(cbuf)) deallocate(cbuf) - -! call gf_free(gfld) -! end subroutine getgb2p2 +end subroutine getgb2p2 !> Read and unpack sections 6 and 7 from a GRIB2 message using a !> version 1 index record. @@ -1021,11 +1075,14 @@ subroutine getgb2rp(lugb, cindex, extract, gribm, leng, iret) logical, intent(in) :: extract character(len = 1), pointer, dimension(:) :: gribm integer, intent(out) :: leng, iret + integer (kind = 8) :: leng8 + integer :: idxver interface subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng8, iret) - integer, intent(in) :: lugb, idxver + integer, intent(in) :: lugb + integer, intent(inout) :: idxver character(len = 1), intent(in) :: cindex(*) logical, intent(in) :: extract character(len = 1), pointer, dimension(:) :: gribm @@ -1036,7 +1093,8 @@ end subroutine getgb2rp2 ! Call the legacy version of this function. It will only work with ! GRIB messages < 2 GB. - call getgb2rp2(lugb, 1, cindex, extract, gribm, leng8, iret) + idxver = 1 + call getgb2rp2(lugb, idxver, cindex, extract, gribm, leng8, iret) leng = int(leng8, kind(4)) end subroutine getgb2rp @@ -1076,12 +1134,13 @@ subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng8, iret) use g2logging implicit none - integer, intent(in) :: lugb, idxver + integer, intent(in) :: lugb + integer, intent(inout) :: idxver character(len = 1), intent(in) :: cindex(*) logical, intent(in) :: extract + character(len = 1), pointer, dimension(:) :: gribm integer (kind = 8), intent(out) :: leng8 integer, intent(out) :: iret - character(len = 1), pointer, dimension(:) :: gribm integer, parameter :: zero = 0 character(len = 1), allocatable, dimension(:) :: csec2, csec6, csec7 @@ -1112,11 +1171,15 @@ subroutine g2_gbytec1(in, siout, iskip, nbits) 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 + end subroutine g2_gbytec81 end interface #ifdef LOGGING - write(g2_log_msg, '(a, i2, a, i1, a, l)') 'getgb2rp2: lugb ', lugb, ' idxver ', idxver, & - ' extract ', extract + write(g2_log_msg, *) 'getgb2rp2: lugb ', lugb, ' idxver ', idxver, ' extract ', extract call g2_log(1) #endif @@ -1142,7 +1205,7 @@ end subroutine g2_gbytec1 iskip = int(iskip8, kind(4)) call g2_gbytec81(cindex, iskp2_8, mypos, INT8_BITS) ! bytes to skip for section 2 mypos = mypos + INT8_BITS - mypos = mypos + 36 * INT1_BITS ! skip ahead in the cindex + mypos = mypos + 44 * INT1_BITS ! skip ahead in the cindex endif if (iskp2_8 .gt. 0) then call bareadl(lugb, iskip8 + iskp2_8, 4_8, lread8, ctemp) @@ -1154,19 +1217,14 @@ end subroutine g2_gbytec1 len2 = 0 endif call g2_gbytec1(cindex, len1, mypos, INT4_BITS) ! length of section 1 - ipos = 44 + len1 mypos = mypos + len1 * INT1_BITS ! skip ahead in the cindex call g2_gbytec1(cindex, len3, mypos, INT4_BITS) ! length of section 3 - ipos = ipos + len3 mypos = mypos + len3 * INT1_BITS ! skip ahead in the cindex call g2_gbytec1(cindex, len4, mypos, INT4_BITS) ! length of section 4 - ipos = ipos + len4 mypos = mypos + len4 * INT1_BITS ! skip ahead in the cindex call g2_gbytec1(cindex, len5, mypos, INT4_BITS) ! length of section 5 - ipos = ipos + len5 mypos = mypos + len5 * INT1_BITS ! skip ahead in the cindex call g2_gbytec1(cindex, len6, mypos, INT4_BITS) ! length of section 6 - ipos = ipos + 5 mypos = mypos + len6 * INT1_BITS ! skip ahead in the cindex call g2_gbytec1(cindex, ibmap, mypos, INT1_BITS) ! bitmap indicator if (ibmap .eq. 254) then @@ -1199,9 +1257,21 @@ end subroutine g2_gbytec1 len7_8 = len7 call bareadl(lugb, iskip8 + iskp7, len7_8, lread8, csec7) +#ifdef LOGGING + write(g2_log_msg, *) 'getgb2rp2: len0 ', len0, 'len1', len1, 'len2', len2 , 'len3', len3 + call g2_log(2) + write(g2_log_msg, *) 'getgb2rp2: len4', len4, 'len5', len5, 'len5', len5, 'len6', len6, 'len7', len7, 'len8', len8 + call g2_log(2) +#endif + ! Now we know the total length of the grib message. leng8 = len0 + len1 + len2 + len3 + len4 + len5 + len6 + len7 + len8 +#ifdef LOGGING + write(g2_log_msg, *) 'getgb2rp2: len7 ', len7, 'lread8', lread8, 'calculated leng8', leng8 + call g2_log(2) +#endif + ! Allocate storage for the message. if (.not. associated(gribm)) allocate(gribm(leng8)) @@ -1263,20 +1333,20 @@ end subroutine g2_gbytec1 if (allocated(csec7)) deallocate(csec7) else ! do not extract field from message : get entire message if (idxver .eq. 1) then - call g2_gbytec(cindex, iskip, mypos, INT4_BITS) ! bytes to skip in file + call g2_gbytec1(cindex, iskip, mypos, INT4_BITS) ! bytes to skip in file mypos = mypos + INT4_BITS mypos = mypos + 6 * INT4_BITS iskip8 = iskip else - call g2_gbytec8(cindex, iskip8, mypos, INT8_BITS) ! bytes to skip in file + call g2_gbytec81(cindex, iskip8, mypos, INT8_BITS) ! bytes to skip in file mypos = mypos + INT8_BITS - mypos = mypos + 2 * INT8_BITS + 4 * INT4_BITS + mypos = mypos + 4 * INT8_BITS + 2 * INT4_BITS endif ! Get the length of the GRIB2 message from the index. call g2_gbytec8(cindex, leng8, mypos, INT8_BITS) #ifdef LOGGING - write(g2_log_msg, *) ' iskip8 ', iskip8, ' mypos/8 ', mypos/8, ' leng8 ', leng8 + write(g2_log_msg, *) ' iskip8 ', iskip8, ' mypos/8 ', mypos/8, 'index leng8 ', leng8 call g2_log(2) #endif diff --git a/src/g2index.F90 b/src/g2index.F90 index 4c65764d..7cce6b8e 100644 --- a/src/g2index.F90 +++ b/src/g2index.F90 @@ -520,7 +520,6 @@ subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret) integer :: ios, istat, lbuf, lhead, nskp #ifdef LOGGING - ! Log results for debugging. write(g2_log_msg, *) 'getg2i2: lugi ', lugi call g2_log(1) #endif @@ -530,9 +529,17 @@ subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret) 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 - if (ios .eq. 0) then +#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 @@ -668,8 +675,7 @@ end subroutine ix2gb2 #ifdef LOGGING ! Log results for debugging. - write(g2_log_msg, '(a, i2, a, i7, a, i7, a, i5, a, i1)') 'getg2i2r: lugb ', lugb, ' msk1 ', msk1, ' msk2 ', msk2, & - ' mnum ', mnum, ' idxver ', idxver + write(g2_log_msg, *) 'getg2i2r: lugb ', lugb, ' msk1 ', msk1, ' msk2 ', msk2, 'idxver', idxver call g2_log(1) #endif @@ -984,7 +990,7 @@ end subroutine gf_unpack5 #ifdef LOGGING ! Log results for debugging. - write(g2_log_msg, '(a, i1, a, i5, a, i7, a, i3, a, i3)') 'getgb2s2: idxver ', idxver, ' nlen ', nlen, & + write(g2_log_msg, *) 'getgb2s2: idxver ', idxver, ' nlen ', nlen, & ' nnum ', nnum, ' j ', j, ' jdisc ', jdisc call g2_log(1) #endif @@ -1452,7 +1458,7 @@ end subroutine g2_gbytec1 #ifdef LOGGING write(g2_log_msg, *) ' writing pds location to index: mypos/8 ', mypos/8, & ' loc ', int(ibskip8 - lskip8, kind(4)) - call g2_log(2) + call g2_log(4) #endif !print '(i3, a8, i4)', mypos/8, ' locpds ', int(ibskip8 - lskip8, kind(4)) mypos = mypos + INT4_BITS @@ -1477,7 +1483,7 @@ end subroutine g2_gbytec1 mypos = mypos + INT4_BITS * 3 ! skip ahead in cbuf #ifdef LOGGING write(g2_log_msg, *) ' writing total len to index: mypos/8 ', mypos/8, lgrib8 - call g2_log(2) + call g2_log(4) #endif call g2_sbytec8(cindex, lgrib8, mypos, INT8_BITS) ! len of grib2 !print '(i3, a8, i4)', mypos/8, ' lgrib8 ', lgrib8 @@ -1522,14 +1528,14 @@ end subroutine g2_gbytec1 mypos = mypos + ilnpds #ifdef LOGGING write(g2_log_msg, *) ' after writing pds location to index: mypos/8 ', mypos/8 - call g2_log(3) + call g2_log(4) #endif elseif (numsec .eq. 5) then ! Write the byte offset to the DRS section into the cindex buffer. !mypos = (IXSDR + inc) * INT1_BITS #ifdef LOGGING write(g2_log_msg, *) ' before writing drs to index: ibskip8 - lskip8 ', ibskip8 - lskip8, IXDRS2 - call g2_log(3) + call g2_log(4) #endif ! Write the bytes to skip to the drs section into the index record. if (idxver .eq. 1) then diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 45def016..3aac0281 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -121,6 +121,7 @@ if(FTP_TEST_FILES) create_test(test_getgb2rp_2 ${kind}) create_test(test_getgb2s ${kind}) create_test(test_getgb2p ${kind}) + create_test(test_getgb2p2 ${kind}) create_test(test_getgb2r ${kind}) create_test(test_file_blend ${kind}) create_test(test_aqm ${kind}) diff --git a/tests/test_getgb2p2.F90 b/tests/test_getgb2p2.F90 new file mode 100644 index 00000000..2f470cd4 --- /dev/null +++ b/tests/test_getgb2p2.F90 @@ -0,0 +1,104 @@ +! This is a test program for NCEPLIBS-g2. +! +! This program tests getg2p2(). +! +! Edward Hartnett 10/21/24 +program test_getgb2p2 + use bacio_module + !use g2logging + implicit none + + integer :: lugi + integer :: lugb = 3 + integer :: iret + integer (kind = 8) :: leng8 + character(len=1), pointer, dimension(:) :: gribm + integer :: j, jdisc, jpdtn, jgdtn + integer :: jids(13), jpdt(100), jgdt(250) + logical :: extract + integer :: k + integer :: i + integer :: idxver, test_idx + + ! Interfaces are needed due to pointers in the parameter lists. + interface + subroutine getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng8, iret) + integer, intent(in) :: lugb, lugi, j, jdisc + integer, dimension(:) :: jids(*) + integer, intent(in) :: jpdtn + integer, dimension(:) :: jpdt(*) + integer, intent(in) :: jgdtn + integer, dimension(:) :: jgdt(*) + logical, intent(in) :: extract + integer, intent(inout) :: idxver + integer, intent(out) :: k + character(len = 1), pointer, dimension(:) :: gribm + integer (kind = 8), intent(out) :: leng8 + integer, intent(out) :: iret + end subroutine getgb2p2 + end interface + + print *, 'Testing the getgb2p2() subroutine...' + + ! Initialize search values to find the first message in the file. + lugi = 0 + j = 0 + jdisc = -1 + do i = 1, 13 + jids(i) = -9999 + end do + jpdtn = -1 + do i = 1, 100 + jpdt(i) = -9999 + end do + jgdtn = -1 + do i = 1, 250 + jgdt(i) = -9999 + end do + + ! Test with index version 1 and 2. + do test_idx = 1, 2 + ! Open a real GRIB2 file. + print *, 'Indexing a real GRIB2 file WW3_Regional_US_West_Coast_20220718_0000.grib2...' + call baopenr(lugb, "data/WW3_Regional_US_West_Coast_20220718_0000.grib2", iret) + if (iret .ne. 0) stop 100 + + !g2_log_level = 3 + extract = .true. + idxver = test_idx + print *, 'Try getgb2p2() with extract true, idxver:', idxver + + nullify(gribm) + call getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng8, iret) + if (iret .ne. 0) then + print *, iret + stop 101 + endif + if (k .ne. 1 .or. leng8 .ne. 11183) then + print *, k, leng8 + stop 110 + endif + + ! Deallocate buffer that got GRIB message. + deallocate(gribm) + print *, 'OK!' + + print *, 'Try getgb2p2() with extract false, idxver:', idxver + extract = .false. + nullify(gribm) + call getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng8, iret) + if (iret .ne. 0) stop 101 + if (k .ne. 1 .or. leng8 .ne. 11183) stop 110 + + deallocate(gribm) + print *, 'OK!' + + call baclose(lugb, iret) + if (iret .ne. 0) stop 199 + end do + print *, 'SUCCESS!...' + +end program test_getgb2p2