diff --git a/.github/workflows/developer.yml b/.github/workflows/developer.yml index 4df37c1d..66d264c0 100644 --- a/.github/workflows/developer.yml +++ b/.github/workflows/developer.yml @@ -103,7 +103,7 @@ jobs: uses: actions/cache@v4 with: path: ~/data - key: data-developer-2 + key: data-developer-3 - name: asan if: matrix.config == 'asan/code coverage' diff --git a/src/g2getgb2.F90 b/src/g2getgb2.F90 index b1fbf05e..5ce3d5db 100644 --- a/src/g2getgb2.F90 +++ b/src/g2getgb2.F90 @@ -430,6 +430,110 @@ end subroutine gf_unpack2 endif end subroutine getgb2l2 +!> Find and extract a GRIB2 message from a file. +!> +!> 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 + + interface + subroutine getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng, iret) + use grib_mod + 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, intent(out) :: leng, iret + end subroutine getgb2p2 + end interface + + ! Call the updated version of this subroutine, which handles index + ! version 2. But if this subroutine is called, then use version 1. + idxver = 1 + call getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng, iret) +end subroutine getgb2p + !> Find and extract a GRIB2 message from a file. !> !> This subroutine reads a GRIB index file (or optionally the GRIB @@ -507,6 +611,8 @@ 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 Index version. Use version 2. Version 1 is for +!> legacy purposes. !> @param[out] k field number unpacked. !> @param[out] gribm returned GRIB message. !> @param[out] leng length of returned GRIB message in bytes. @@ -516,37 +622,68 @@ end subroutine getgb2l2 !> - 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 Ed Hartnett @date 05-14-2024 +subroutine getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & + extract, idxver, k, gribm, leng, 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, intent(out) :: leng, iret + type(gribfield) :: gfld integer :: msk1, irgi, irgs, jk, lpos, msk2, mskp, nlen, nmess, nnum + integer (kind = 8) :: msk1_8, msk2_8 + parameter(msk1_8 = 32000_8, msk2_8 = 4000_8) character(len = 1), pointer, dimension(:) :: cbuf 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 + subroutine getg2i2(lugi, cbuf, idxver, nlen, nnum, iret) integer, intent(in) :: lugi - integer, intent(out) :: nlen, nnum, iret - end subroutine getg2i - subroutine getg2ir(lugb, msk1, msk2, mnum, cbuf, nlen, nnum, & - nmess, iret) + 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(in) :: lugb, msk1, msk2, mnum integer, intent(out) :: nlen, nnum, nmess, iret - end subroutine getg2ir + 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 subroutine getgb2rp(lugb, cindex, extract, gribm, leng, iret) integer, intent(in) :: lugb character(len = 1), intent(in) :: cindex(*) @@ -559,10 +696,10 @@ end subroutine getgb2rp ! 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_8, msk2_8, mskp, idxver, cbuf, nlen, nnum, nmess, irgi) endif if (irgi .gt. 1) then iret = 96 @@ -570,7 +707,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 @@ -588,7 +725,7 @@ end subroutine getgb2rp if (associated(cbuf)) deallocate(cbuf) call gf_free(gfld) -end subroutine getgb2p +end subroutine getgb2p2 ! subroutine getgb2p2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, & ! extract, k, gribm, leng, iret) @@ -949,17 +1086,22 @@ subroutine getgb2rp(lugb, cindex, extract, gribm, leng, iret) character(len = 1), pointer, dimension(:) :: gribm integer, intent(out) :: leng, iret + integer (kind = 8) :: leng8 + interface - subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng, iret) + subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng8, 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 + integer(kind = 8), intent(out) :: leng8 + integer, intent(out) :: iret end subroutine getgb2rp2 end interface - call getgb2rp2(lugb, 1, cindex, extract, gribm, leng, iret) + ! Call with idxver = 1 for legacy purposes. + call getgb2rp2(lugb, 1, cindex, extract, gribm, leng8, iret) + leng = int(leng8, kind(4)) end subroutine getgb2rp @@ -994,13 +1136,14 @@ end subroutine getgb2rp !> - 97 Error reading grib file. !> !> @author Edward Hartnett, Stephen Gilbert @date Feb 13, 2024 -subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng, iret) +subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng8, iret) implicit none integer, intent(in) :: lugb, idxver character(len = 1), intent(in) :: cindex(*) logical, intent(in) :: extract - integer, intent(out) :: leng, iret + integer(kind = 8), intent(out) :: leng8 + integer, intent(out) :: iret character(len = 1), pointer, dimension(:) :: gribm integer, parameter :: zero = 0 @@ -1013,7 +1156,7 @@ subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng, iret) integer :: INT1_BITS, INT2_BITS, INT4_BITS, INT8_BITS parameter(INT1_BITS = 8, INT2_BITS = 16, INT4_BITS = 32, INT8_BITS = 64) integer :: mypos, inc = 0 - integer (kind = 8) :: lread8, iskip8, leng8, len2_8, len7_8, len6_8 + integer (kind = 8) :: lread8, iskip8, len2_8, len7_8, len6_8 iret = 0 @@ -1080,8 +1223,8 @@ subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng, iret) len7_8 = len7 call bareadl(lugb, iskip8 + iskp7, len7_8, lread8, csec7) - leng = len0 + len1 + len2 + len3 + len4 + len5 + len6 + len7 + len8 - if (.not. associated(gribm)) allocate(gribm(leng)) + leng8 = len0 + len1 + len2 + len3 + len4 + len5 + len6 + len7 + len8 + if (.not. associated(gribm)) allocate(gribm(leng8)) ! Create Section 0 gribm(1) = 'G' @@ -1096,7 +1239,7 @@ subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng, iret) gribm(10) = char(0) gribm(11) = char(0) gribm(12) = char(0) - call g2_sbytec(gribm, leng, 12*8, INT4_BITS) + call g2_sbytec8(gribm, leng8, 12*8, INT8_BITS) ! Copy Section 1 gribm(17:16 + len1) = cindex(45 + inc:44 + inc + len1) @@ -1153,9 +1296,8 @@ subroutine getgb2rp2(lugb, idxver, cindex, extract, gribm, leng, iret) mypos = mypos + INT8_BITS endif mypos = mypos + 7 * INT4_BITS - call g2_gbytec(cindex, leng, mypos, INT4_BITS) ! length of grib message - if (.not. associated(gribm)) allocate(gribm(leng)) - leng8 = leng + call g2_gbytec8(cindex, leng8, mypos, INT8_BITS) ! length of grib message + if (.not. associated(gribm)) allocate(gribm(leng8)) call bareadl(lugb, iskip8, leng8, lread8, gribm) if (leng8 .ne. lread8) then deallocate(gribm) diff --git a/src/g2index.F90 b/src/g2index.F90 index 3234a90f..9b6e3c00 100644 --- a/src/g2index.F90 +++ b/src/g2index.F90 @@ -1232,7 +1232,7 @@ subroutine ix2gb2(lugb, lskip8, idxver, lgrib8, cbuf, numfld, mlen, iret) integer (kind = 8) :: ibread8, lbread8, ibskip8, lengds8 integer (kind = 8) :: ilnpds8, ilndrs8 integer :: lensec, lensec1 - integer :: mypos, inc, i + integer :: mypos, inc ! Parameters. ! Size of the internal char buffers used in this subroutine. diff --git a/tests/test_getgb2p_2.F90 b/tests/test_getgb2p_2.F90 index d98bef21..93ac7650 100755 --- a/tests/test_getgb2p_2.F90 +++ b/tests/test_getgb2p_2.F90 @@ -53,9 +53,7 @@ END SUBROUTINE GETGB2P fileo='test_tocgrib2.output.grib2' call baopenw(lugo,fileo,iret1) if (iret1 .ne. 0) then - write(6,fmt='(" Error opening output transmission file: ", & - A200)') fileo - write(6,fmt='(" baopenw error = ",I5)') iret1 + print *, " Error opening output transmission file: ", fileo stop 20 endif @@ -77,15 +75,12 @@ END SUBROUTINE GETGB2P READ (12,GRIBIDS,iostat=ios) if (ios .ne. 0) then - write(6,fmt='(" Error reading PDS from input file. iostat = " & - ,i5)') ios + print *, " Error reading PDS from input file. iostat = ", ios cycle endif nrec = nrec + 1 ! Echo input record - WRITE(6,FMT='(/,''***********************************'', & - ''********************************************'')') write(6,'(A,I0)') ' Start new record no. = ',nrec write(6,'(73A)') ' DESC=',DESC(1:73) write(6,'(11A)') ' WMOHEAD=',WMOHEAD(1:11)