From 02f633c91cac291924ec437bb73abb12972e88a1 Mon Sep 17 00:00:00 2001 From: Edward Hartnett <38856240+edwardhartnett@users.noreply.github.com> Date: Tue, 5 Sep 2023 12:50:47 +0200 Subject: [PATCH] converted gf_unpack4, 5, and 7 to F90 (#548) * converted gf_unpack4, 5, and 7 to F90 * converted ixbg2() to F90 --- src/CMakeLists.txt | 4 +- src/gf_unpack4.F90 | 141 +++++++++++++++++++++++++++++++ src/gf_unpack4.f | 142 -------------------------------- src/gf_unpack5.F90 | 120 +++++++++++++++++++++++++++ src/gf_unpack5.f | 119 --------------------------- src/gf_unpack7.F90 | 102 +++++++++++++++++++++++ src/gf_unpack7.f | 110 ------------------------- src/ixgb2.F90 | 197 ++++++++++++++++++++++++++++++++++++++++++++ src/ixgb2.f | 201 --------------------------------------------- 9 files changed, 562 insertions(+), 574 deletions(-) create mode 100644 src/gf_unpack4.F90 delete mode 100644 src/gf_unpack4.f create mode 100644 src/gf_unpack5.F90 delete mode 100644 src/gf_unpack5.f create mode 100644 src/gf_unpack7.F90 delete mode 100644 src/gf_unpack7.f create mode 100644 src/ixgb2.F90 delete mode 100644 src/ixgb2.f diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1fb50a26..cb61a8ef 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -9,10 +9,10 @@ comunpack.f drstemplates.F90 g2_gbytesc.F90 g2grids.F90 gb_info.F90 getdim.f getfield.F90 getg2i.F90 getg2ir.F90 getgb2.F90 getgb2l.F90 getgb2p.F90 getgb2r.F90 getgb2rp.F90 getgb2s.F90 getidx.F90 getlocal.F90 getpoly.f gettemplates.F90 gf_free.F90 gf_getfld.F90 gf_unpack1.F90 gf_unpack2.F90 -gf_unpack3.F90 gf_unpack4.f gf_unpack5.f gf_unpack6.F90 gf_unpack7.f +gf_unpack3.F90 gf_unpack4.F90 gf_unpack5.F90 gf_unpack6.F90 gf_unpack7.F90 gribcreate.F90 gribend.F90 gribinfo.F90 ${CMAKE_CURRENT_BINARY_DIR}/gribmod.F90 gridtemplates.F90 intmath.f -ixgb2.f jpcpack.F90 jpcunpack.F90 misspack.f mkieee.F90 pack_gp.f +ixgb2.F90 jpcpack.F90 jpcunpack.F90 misspack.f mkieee.F90 pack_gp.f params_ecmwf.F90 params.F90 pdstemplates.F90 pngpack.F90 pngunpack.F90 putgb2.F90 rdieee.F90 realloc.f reduce.f simpack.f simunpack.F90 skgb.F90 specpack.F90 specunpack.F90) diff --git a/src/gf_unpack4.F90 b/src/gf_unpack4.F90 new file mode 100644 index 00000000..f07eda5a --- /dev/null +++ b/src/gf_unpack4.F90 @@ -0,0 +1,141 @@ +!> @file +!> @brief Unpack Section 4 ([Product Definition +!> Section] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect4.shtml)) +!> of a GRIB2 message. +!> @author Stephen Gilbert @date 2000-05-26 + +!> Unpack Section 4 ([Product Definition Section] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect4.shtml)) +!> of a GRIB2 message, starting at octet 6 of that Section. +!> +!> @param[in] cgrib Character array that contains the GRIB2 message. +!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib. +!> @param[inout] iofst Bit offset of the beginning/end(returned) of +!> Section 4. +!> @param[out] ipdsnum Product Definition Template Number ([Code Table 4.0] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-0.shtml)). +!> @param[out] ipdstmpl Contains the data values for the +!> Product Definition Template specified by ipdsnum. A safe dimension for this array +!> can be obtained in advance from maxvals(4), which is returned from +!> subroutine gribinfo. +!> @param[out] mappdslen Number of elements in ipdstmpl. i.e. number +!> of entries in Product Defintion Template specified by ipdsnum. +!> @param[out] coordlist Pointer to real array containing floating +!> point values intended to document the vertical discretisation +!> associated to model data on hybrid coordinate vertical +!> levels (part of Section 4). Must be deallocated by caller. +!> @param[out] numcoord Number of values in array coordlist. +!> @param[out] ierr Error return code. +!> - 0 no error. +!> - 5 "GRIB" message contains an undefined Grid Definition Template. +!> - 6 memory allocation error. +!> +!> @author Stephen Gilbert @date 2000-05-26 +subroutine gf_unpack4(cgrib, lcgrib, iofst, ipdsnum, ipdstmpl, & + mappdslen, coordlist, numcoord, ierr) + use pdstemplates + use re_alloc ! needed for subroutine realloc + implicit none + + 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 + + real(4), allocatable :: coordieee(:) + integer, allocatable :: mappds(:) + integer :: mappdslen + logical needext + integer :: lensec, nbits, newmappdslen + integer :: istat1, istat, isign, iret, i + + ierr = 0 + nullify(ipdstmpl, coordlist) + + ! Get Length of Section. + call g2_gbytec(cgrib, lensec, iofst, 32) + iofst = iofst + 32 + iofst = iofst + 8 ! skip section number + allocate(mappds(lensec)) + + ! Get num of coordinate values. + call g2_gbytec(cgrib, numcoord, iofst, 16) + iofst = iofst + 16 + ! Get Prod. Def Template num. + call g2_gbytec(cgrib, ipdsnum, iofst, 16) + iofst = iofst + 16 + ! Get Product Definition Template. + call getpdstemplate(ipdsnum, mappdslen, mappds, needext, iret) + if (iret.ne.0) then + ierr = 5 + if (allocated(mappds)) deallocate(mappds) + return + endif + + ! Unpack each value into array ipdstmpl from the the appropriate + ! number of octets, which are specified in corresponding entries in + ! array mappds. + istat = 0 + if (mappdslen.gt.0) allocate(ipdstmpl(mappdslen), stat = istat) + if (istat.ne.0) then + ierr = 6 + nullify(ipdstmpl) + if (allocated(mappds)) deallocate(mappds) + return + endif + do i = 1, mappdslen + nbits = iabs(mappds(i))*8 + if (mappds(i).ge.0) then + call g2_gbytec(cgrib, ipdstmpl(i), iofst, nbits) + else + call g2_gbytec(cgrib, isign, iofst, 1) + call g2_gbytec(cgrib, ipdstmpl(i), iofst + 1, nbits-1) + if (isign.eq.1) ipdstmpl(i) = -ipdstmpl(i) + endif + iofst = iofst + nbits + enddo + + ! Check to see if the Product Definition Template needs to be + ! extended. The number of values in a specific template may vary + ! depending on data specified in the "static" part of the template. + if (needext) then + call extpdstemplate(ipdsnum, ipdstmpl, newmappdslen, mappds) + call realloc(ipdstmpl, mappdslen, newmappdslen, istat) + ! Unpack the rest of the Product Definition Template. + do i = mappdslen + 1, newmappdslen + nbits = iabs(mappds(i))*8 + if (mappds(i).ge.0) then + call g2_gbytec(cgrib, ipdstmpl(i), iofst, nbits) + else + call g2_gbytec(cgrib, isign, iofst, 1) + call g2_gbytec(cgrib, ipdstmpl(i), iofst + 1, nbits-1) + if (isign.eq.1) ipdstmpl(i) = -ipdstmpl(i) + endif + iofst = iofst + nbits + enddo + mappdslen = newmappdslen + endif + if (allocated(mappds)) deallocate(mappds) + + ! Get Optional list of vertical coordinate values + ! after the Product Definition Template, if necessary. + nullify(coordlist) + if (numcoord .ne. 0) then + allocate (coordieee(numcoord), stat = istat1) + allocate(coordlist(numcoord), stat = istat) + if ((istat1 + istat).ne.0) then + ierr = 6 + nullify(coordlist) + if (allocated(coordieee)) deallocate(coordieee) + return + endif + call g2_gbytesc(cgrib, coordieee, iofst, 32, 0, numcoord) + call rdieee(coordieee, coordlist, numcoord) + deallocate (coordieee) + iofst = iofst + (32 * numcoord) + endif +end subroutine gf_unpack4 diff --git a/src/gf_unpack4.f b/src/gf_unpack4.f deleted file mode 100644 index dbb2a050..00000000 --- a/src/gf_unpack4.f +++ /dev/null @@ -1,142 +0,0 @@ -!> @file -!> @brief Unpack Section 4 ([Product Definition -!> Section] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect4.shtml)) -!> of a GRIB2 message. -!> @author Stephen Gilbert @date 2000-05-26 - -!> Unpack Section 4 ([Product Definition Section] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect4.shtml)) -!> of a GRIB2 message, starting at octet 6 of that Section. -!> -!> @param[in] cgrib Character array that contains the GRIB2 message. -!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib. -!> @param[inout] iofst Bit offset of the beginning/end(returned) of -!> Section 4. -!> @param[out] ipdsnum Product Definition Template Number ([Code Table 4.0] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-0.shtml)). -!> @param[out] ipdstmpl Contains the data values for the -!> Product Definition Template specified by ipdsnum. A safe dimension for this array -!> can be obtained in advance from maxvals(4), which is returned from -!> subroutine gribinfo. -!> @param[out] mappdslen Number of elements in ipdstmpl. i.e. number -!> of entries in Product Defintion Template specified by ipdsnum. -!> @param[out] coordlist Pointer to real array containing floating -!> point values intended to document the vertical discretisation -!> associated to model data on hybrid coordinate vertical -!> levels (part of Section 4). Must be deallocated by caller. -!> @param[out] numcoord Number of values in array coordlist. -!> @param[out] ierr Error return code. -!> - 0 no error. -!> - 5 "GRIB" message contains an undefined Grid Definition Template. -!> - 6 memory allocation error. -!> -!> @author Stephen Gilbert @date 2000-05-26 - subroutine gf_unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl, - & mappdslen,coordlist,numcoord,ierr) - - use pdstemplates - use re_alloc ! needed for subroutine realloc - implicit none - - 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 - - real(4),allocatable :: coordieee(:) - integer,allocatable :: mappds(:) - integer :: mappdslen - logical needext - integer :: lensec, nbits, newmappdslen - integer :: istat1, istat, isign, iret, i - - ierr=0 - nullify(ipdstmpl,coordlist) - - call g2_gbytec(cgrib,lensec,iofst,32) ! Get Length of Section - iofst=iofst+32 - iofst=iofst+8 ! skip section number - allocate(mappds(lensec)) - - call g2_gbytec(cgrib,numcoord,iofst,16) ! Get num of coordinate values - iofst=iofst+16 - call g2_gbytec(cgrib,ipdsnum,iofst,16) ! Get Prod. Def Template num. - iofst=iofst+16 - ! Get Product Definition Template - call getpdstemplate(ipdsnum,mappdslen,mappds,needext,iret) - if (iret.ne.0) then - ierr=5 - if( allocated(mappds) ) deallocate(mappds) - return - endif - -! Unpack each value into array ipdstmpl from the -! the appropriate number of octets, which are specified in -! corresponding entries in array mappds. - istat=0 - if (mappdslen.gt.0) allocate(ipdstmpl(mappdslen),stat=istat) - if (istat.ne.0) then - ierr=6 - nullify(ipdstmpl) - if( allocated(mappds) ) deallocate(mappds) - return - endif - do i=1,mappdslen - nbits=iabs(mappds(i))*8 - if ( mappds(i).ge.0 ) then - call g2_gbytec(cgrib,ipdstmpl(i),iofst,nbits) - else - call g2_gbytec(cgrib,isign,iofst,1) - call g2_gbytec(cgrib,ipdstmpl(i),iofst+1,nbits-1) - if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i) - endif - iofst=iofst+nbits - enddo - - ! Check to see if the Product Definition Template needs to be - ! extended. - ! The number of values in a specific template may vary - ! depending on data specified in the "static" part of the - ! template. - if ( needext ) then - call extpdstemplate(ipdsnum,ipdstmpl,newmappdslen,mappds) - call realloc(ipdstmpl,mappdslen,newmappdslen,istat) - ! Unpack the rest of the Product Definition Template - do i=mappdslen+1,newmappdslen - nbits=iabs(mappds(i))*8 - if ( mappds(i).ge.0 ) then - call g2_gbytec(cgrib,ipdstmpl(i),iofst,nbits) - else - call g2_gbytec(cgrib,isign,iofst,1) - call g2_gbytec(cgrib,ipdstmpl(i),iofst+1,nbits-1) - if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i) - endif - iofst=iofst+nbits - enddo - mappdslen=newmappdslen - endif - if( allocated(mappds) ) deallocate(mappds) - - ! Get Optional list of vertical coordinate values - ! after the Product Definition Template, if necessary. - nullify(coordlist) - if ( numcoord .ne. 0 ) then - allocate (coordieee(numcoord),stat=istat1) - allocate(coordlist(numcoord),stat=istat) - if ((istat1+istat).ne.0) then - ierr=6 - nullify(coordlist) - if( allocated(coordieee) ) deallocate(coordieee) - return - endif - call g2_gbytesc(cgrib,coordieee,iofst,32,0,numcoord) - call rdieee(coordieee,coordlist,numcoord) - deallocate (coordieee) - iofst=iofst+(32*numcoord) - endif - end - diff --git a/src/gf_unpack5.F90 b/src/gf_unpack5.F90 new file mode 100644 index 00000000..f818834e --- /dev/null +++ b/src/gf_unpack5.F90 @@ -0,0 +1,120 @@ +!> @file +!> @brief Unpack Section 5 ([Data Representation +!> Section] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect5.shtml)) +!> of a GRIB2 message. +!> @author Stephen Gilbert @date 2000-05-26 + +!> Unpack Section 5 ([Data Representation Section] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect5.shtml)) +!> of a GRIB2 message, starting at octet 6 of that Section. +!> +!> @param[in] cgrib Character array that contains the GRIB2 message. +!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib. +!> @param[inout] iofst Bit offset of the beginning/end(returned) of +!> Section 5. +!> @param[out] ndpts Number of data points unpacked and returned. +!> @param[out] idrsnum Data Representation Template Number ([Code Table 5.0] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-0.shtml)). +!> @param[out] idrstmpl Contains the data values for the specified +!> Data Representation Template (N = idrsnum). Each element of this +!> integer array contains an entry (in the order specified) of +!> Product Defintion Template 5.N A safe dimension for this array can +!> be obtained in advance from maxvals(6), which is returned from +!> subroutine gribinfo. +!> @param[out] mapdrslen Number of elements in idrstmpl. i.e. number +!> of entries in Data Representation Template 5.N (N = idrsnum). +!> @param[out] ierr Error return code. +!> - 0 no error. +!> - 6 memory allocation error. +!> - 7 "GRIB" message contains an undefined Grid Definition Template. +!> +!> @author Stephen Gilbert @date 2000-05-26 +subroutine gf_unpack5(cgrib, lcgrib, iofst, ndpts, idrsnum, idrstmpl, & + mapdrslen, ierr) + use drstemplates + use re_alloc ! needed for subroutine realloc + implicit none + + 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 + + integer, allocatable :: mapdrs(:) + integer :: mapdrslen + logical needext + integer :: newmapdrslen, nbits, istat, isign, lensec, iret, i + + ierr = 0 + nullify(idrstmpl) + + call g2_gbytec(cgrib, lensec, iofst, 32) ! Get Length of Section + iofst = iofst + 32 + iofst = iofst + 8 ! skip section number + allocate(mapdrs(lensec)) + + ! Get num of data points. + call g2_gbytec(cgrib, ndpts, iofst, 32) + iofst = iofst + 32 + ! Get Data Rep Template Num. + call g2_gbytec(cgrib, idrsnum, iofst, 16) + iofst = iofst + 16 + ! Gen Data Representation Template. + call getdrstemplate(idrsnum, mapdrslen, mapdrs, needext, iret) + if (iret .ne. 0) then + ierr = 7 + if (allocated(mapdrs)) deallocate(mapdrs) + return + endif + + ! Unpack each value into array ipdstmpl from the the appropriate + ! number of octets, which are specified in corresponding entries in + ! array mappds. + istat = 0 + if (mapdrslen .gt. 0) allocate(idrstmpl(mapdrslen), stat = istat) + if (istat .ne. 0) then + ierr = 6 + nullify(idrstmpl) + if (allocated(mapdrs)) deallocate(mapdrs) + return + endif + do i = 1, mapdrslen + nbits = iabs(mapdrs(i)) * 8 + if (mapdrs(i) .ge. 0) then + call g2_gbytec(cgrib, idrstmpl(i), iofst, nbits) + else + call g2_gbytec(cgrib, isign, iofst, 1) + call g2_gbytec(cgrib, idrstmpl(i), iofst + 1, nbits-1) + if (isign .eq. 1) idrstmpl(i) = -idrstmpl(i) + endif + iofst = iofst + nbits + enddo + + ! Check to see if the Data Representation Template needs to be + ! extended. The number of values in a specific template may vary + ! depending on data specified in the "static" part of the template. + if (needext) then + call extdrstemplate(idrsnum, idrstmpl, newmapdrslen, mapdrs) + call realloc(idrstmpl, mapdrslen, newmapdrslen, istat) + + ! Unpack the rest of the Data Representation Template. + do i = mapdrslen + 1, newmapdrslen + nbits = iabs(mapdrs(i)) * 8 + if (mapdrs(i) .ge. 0) then + call g2_gbytec(cgrib, idrstmpl(i), iofst, nbits) + else + call g2_gbytec(cgrib, isign, iofst, 1) + call g2_gbytec(cgrib, idrstmpl(i), iofst + 1, nbits - 1) + if (isign.eq.1) idrstmpl(i) = -idrstmpl(i) + endif + iofst = iofst + nbits + enddo + mapdrslen = newmapdrslen + endif + if (allocated(mapdrs)) deallocate(mapdrs) + +end subroutine gf_unpack5 + diff --git a/src/gf_unpack5.f b/src/gf_unpack5.f deleted file mode 100644 index d547a80d..00000000 --- a/src/gf_unpack5.f +++ /dev/null @@ -1,119 +0,0 @@ -!> @file -!> @brief Unpack Section 5 ([Data Representation -!> Section] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect5.shtml)) -!> of a GRIB2 message. -!> @author Stephen Gilbert @date 2000-05-26 - -!> Unpack Section 5 ([Data Representation Section] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect5.shtml)) -!> of a GRIB2 message, starting at octet 6 of that Section. -!> -!> @param[in] cgrib Character array that contains the GRIB2 message. -!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib. -!> @param[inout] iofst Bit offset of the beginning/end(returned) of -!> Section 5. -!> @param[out] ndpts Number of data points unpacked and returned. -!> @param[out] idrsnum Data Representation Template Number ([Code Table 5.0] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-0.shtml)). -!> @param[out] idrstmpl Contains the data values for the specified -!> Data Representation Template (N=idrsnum). Each element of this -!> integer array contains an entry (in the order specified) of -!> Product Defintion Template 5.N A safe dimension for this array can -!> be obtained in advance from maxvals(6), which is returned from -!> subroutine gribinfo. -!> @param[out] mapdrslen Number of elements in idrstmpl. i.e. number -!> of entries in Data Representation Template 5.N (N=idrsnum). -!> @param[out] ierr Error return code. -!> - 0 no error. -!> - 6 memory allocation error. -!> - 7 "GRIB" message contains an undefined Grid Definition Template. -!> -!> @author Stephen Gilbert @date 2000-05-26 - subroutine gf_unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,idrstmpl, - & mapdrslen,ierr) - use drstemplates - use re_alloc ! needed for subroutine realloc - implicit none - - 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 - - integer,allocatable :: mapdrs(:) - integer :: mapdrslen - logical needext - integer :: newmapdrslen, nbits, istat, isign, lensec, iret, i - - ierr=0 - nullify(idrstmpl) - - call g2_gbytec(cgrib,lensec,iofst,32) ! Get Length of Section - iofst=iofst+32 - iofst=iofst+8 ! skip section number - allocate(mapdrs(lensec)) - - call g2_gbytec(cgrib,ndpts,iofst,32) ! Get num of data points - iofst=iofst+32 - call g2_gbytec(cgrib,idrsnum,iofst,16) ! Get Data Rep Template Num. - iofst=iofst+16 - ! Gen Data Representation Template - call getdrstemplate(idrsnum,mapdrslen,mapdrs,needext,iret) - if (iret.ne.0) then - ierr=7 - if( allocated(mapdrs) ) deallocate(mapdrs) - return - endif - - ! Unpack each value into array ipdstmpl from the - ! the appropriate number of octets, which are specified in - ! corresponding entries in array mappds. - istat=0 - if (mapdrslen.gt.0) allocate(idrstmpl(mapdrslen),stat=istat) - if (istat.ne.0) then - ierr=6 - nullify(idrstmpl) - if( allocated(mapdrs) ) deallocate(mapdrs) - return - endif - do i=1,mapdrslen - nbits=iabs(mapdrs(i))*8 - if ( mapdrs(i).ge.0 ) then - call g2_gbytec(cgrib,idrstmpl(i),iofst,nbits) - else - call g2_gbytec(cgrib,isign,iofst,1) - call g2_gbytec(cgrib,idrstmpl(i),iofst+1,nbits-1) - if (isign.eq.1) idrstmpl(i)=-idrstmpl(i) - endif - iofst=iofst+nbits - enddo - - ! Check to see if the Data Representation Template needs to be - ! extended. - ! The number of values in a specific template may vary - ! depending on data specified in the "static" part of the - ! template. - if ( needext ) then - call extdrstemplate(idrsnum,idrstmpl,newmapdrslen,mapdrs) - call realloc(idrstmpl,mapdrslen,newmapdrslen,istat) - ! Unpack the rest of the Data Representation Template - do i=mapdrslen+1,newmapdrslen - nbits=iabs(mapdrs(i))*8 - if ( mapdrs(i).ge.0 ) then - call g2_gbytec(cgrib,idrstmpl(i),iofst,nbits) - else - call g2_gbytec(cgrib,isign,iofst,1) - call g2_gbytec(cgrib,idrstmpl(i),iofst+1,nbits-1) - if (isign.eq.1) idrstmpl(i)=-idrstmpl(i) - endif - iofst=iofst+nbits - enddo - mapdrslen=newmapdrslen - endif - if( allocated(mapdrs) ) deallocate(mapdrs) - - end - diff --git a/src/gf_unpack7.F90 b/src/gf_unpack7.F90 new file mode 100644 index 00000000..4b04c34b --- /dev/null +++ b/src/gf_unpack7.F90 @@ -0,0 +1,102 @@ +!> @file +!> @brief Unpack Section 7 ([Data Section] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect7.shtml)) +!> of a GRIB2 message. +!> @author Stephen Gilbert @date 2002-01-24 + +!> Unpack Section 7 ([Data Section] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect7.shtml)) +!> of a GRIB2 message. +!> +!> @param[in] cgrib Character array that contains the GRIB2 message. +!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib. +!> @param[inout] iofst Bit offset of the beginning/end(returned) of +!> Section 7. +!> @param[in] igdsnum Grid Definition Template Number (Code Table +!> 3.0) (Only required to unpack DRT 5.51). +!> @param[in] igdstmpl Pointer to an integer array containing the +!> data values for the specified Grid Definition. Template +!> (N = igdsnum). Each element of this integer array contains an entry +!> (in the order specified) of Grid Definition Template 3.N (Only +!> required to unpack DRT 5.51). +!> @param[in] idrsnum Data Representation Template Number ([Code Table +!> 5.0] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-0.shtml)). +!> @param[in] idrstmpl Pointer to an integer array containing the data +!> values for the Data Representation Template specified by idrsnum. A +!> safe dimension for this array can be obtained in advance from +!> maxvals(6), which is returned from subroutine gribinfo. +!> @param[in] ndpts Number of data points unpacked and returned. +!> @param[out] fld Pointer to a real array containing the unpacked +!> data field. +!> @param[out] ierr Error return code. +!> - 0 no error. +!> - 4 Unrecognized Data Representation Template. +!> - 5 One of GDT 3.50 through 3.53 required to unpack DRT 5.51. +!> - 6 memory allocation error. +!> - 7 corrupt section 7. +!> +!> @author Stephen Gilbert @date 2002-01-24 +subroutine gf_unpack7(cgrib, lcgrib, iofst, igdsnum, igdstmpl, & + idrsnum, idrstmpl, ndpts, fld, ierr) + implicit none + + character(len = 1), intent(in) :: cgrib(lcgrib) + integer, intent(in) :: lcgrib, ndpts, igdsnum, idrsnum + integer, intent(inout) :: iofst + integer, pointer, dimension(:) :: igdstmpl, idrstmpl + integer, intent(out) :: ierr + real, pointer, dimension(:) :: fld + integer :: ier, ipos, istat, lensec, ieee + + ierr = 0 + nullify(fld) + + call g2_gbytec(cgrib, lensec, iofst, 32) ! Get Length of Section + iofst = iofst + 32 + iofst = iofst + 8 ! skip section number + + ipos = (iofst/8) + 1 + istat = 0 + allocate(fld(ndpts), stat = istat) + if (istat.ne.0) then + ierr = 6 + return + endif + + if (idrsnum .eq. 0) then + call simunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld) + elseif (idrsnum.eq.2.or.idrsnum.eq.3) then + call comunpack(cgrib(ipos), lensec-5, lensec, idrsnum, idrstmpl, ndpts, fld, ier) + if (ier .ne. 0) then + ierr = 7 + return + endif + elseif (idrsnum .eq. 50) then ! Spectral simple + call simunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts-1, fld(2)) + ieee = idrstmpl(5) + call rdieee(ieee, fld(1), 1) + elseif (idrsnum .eq. 51) then ! Spectral complex + if (igdsnum.ge.50.AND.igdsnum.le.53) then + call specunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, & + igdstmpl(1), igdstmpl(2), igdstmpl(3), fld) + else + print *, 'gf_unpack7: Cannot use GDT 3.', igdsnum, ' to unpack Data Section 5.51.' + ierr = 5 + nullify(fld) + return + endif + elseif (idrsnum .eq. 40 .OR. idrsnum .eq. 40000) then + call jpcunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld) + elseif (idrsnum .eq. 41 .OR. idrsnum .eq. 40010) then + call pngunpack(cgrib(ipos), lensec-5, idrstmpl, ndpts, fld) + else + print *, 'gf_unpack7: Data Representation Template ', idrsnum, ' not yet implemented.' + ierr = 4 + nullify(fld) + return + endif + + iofst = iofst + (8 * lensec) + +end subroutine gf_unpack7 diff --git a/src/gf_unpack7.f b/src/gf_unpack7.f deleted file mode 100644 index 646950c5..00000000 --- a/src/gf_unpack7.f +++ /dev/null @@ -1,110 +0,0 @@ -!> @file -!> @brief Unpack Section 7 ([Data Section] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect7.shtml)) -!> of a GRIB2 message. -!> @author Stephen Gilbert @date 2002-01-24 - -!> Unpack Section 7 ([Data Section] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect7.shtml)) -!> of a GRIB2 message. -!> -!> @param[in] cgrib Character array that contains the GRIB2 message. -!> @param[in] lcgrib Length (in bytes) of GRIB message array cgrib. -!> @param[inout] iofst Bit offset of the beginning/end(returned) of -!> Section 7. -!> @param[in] igdsnum Grid Definition Template Number (Code Table -!> 3.0) (Only required to unpack DRT 5.51). -!> @param[in] igdstmpl Pointer to an integer array containing the -!> data values for the specified Grid Definition. Template -!> (N=igdsnum). Each element of this integer array contains an entry -!> (in the order specified) of Grid Definition Template 3.N (Only -!> required to unpack DRT 5.51). -!> @param[in] idrsnum Data Representation Template Number ([Code Table -!> 5.0] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-0.shtml)). -!> @param[in] idrstmpl Pointer to an integer array containing the data -!> values for the Data Representation Template specified by idrsnum. A -!> safe dimension for this array can be obtained in advance from -!> maxvals(6), which is returned from subroutine gribinfo. -!> @param[in] ndpts Number of data points unpacked and returned. -!> @param[out] fld Pointer to a real array containing the unpacked -!> data field. -!> @param[out] ierr Error return code. -!> - 0 no error. -!> - 4 Unrecognized Data Representation Template. -!> - 5 One of GDT 3.50 through 3.53 required to unpack DRT 5.51. -!> - 6 memory allocation error. -!> - 7 corrupt section 7. -!> -!> @author Stephen Gilbert @date 2002-01-24 - subroutine gf_unpack7(cgrib,lcgrib,iofst,igdsnum,igdstmpl, - & idrsnum,idrstmpl,ndpts,fld,ierr) - implicit none - - character(len=1),intent(in) :: cgrib(lcgrib) - integer,intent(in) :: lcgrib,ndpts,igdsnum,idrsnum - integer,intent(inout) :: iofst - integer,pointer,dimension(:) :: igdstmpl,idrstmpl - integer,intent(out) :: ierr - real,pointer,dimension(:) :: fld - integer :: ier, ipos, istat, lensec, ieee - - ierr=0 - nullify(fld) - - call g2_gbytec(cgrib,lensec,iofst,32) ! Get Length of Section - iofst=iofst+32 - iofst=iofst+8 ! skip section number - - ipos=(iofst/8)+1 - istat=0 - allocate(fld(ndpts),stat=istat) - if (istat.ne.0) then - ierr=6 - return - endif - - if (idrsnum.eq.0) then - call simunpack(cgrib(ipos),lensec-5,idrstmpl,ndpts,fld) - elseif (idrsnum.eq.2.or.idrsnum.eq.3) then - call comunpack(cgrib(ipos),lensec-5,lensec,idrsnum,idrstmpl, - & ndpts,fld,ier) - if ( ier .NE. 0 ) then - ierr=7 - return - endif - elseif (idrsnum.eq.50) then ! Spectral simple - call simunpack(cgrib(ipos),lensec-5,idrstmpl,ndpts-1, - & fld(2)) - ieee=idrstmpl(5) - call rdieee(ieee,fld(1),1) - elseif (idrsnum.eq.51) then ! Spectral complex - if (igdsnum.ge.50.AND.igdsnum.le.53) then - call specunpack(cgrib(ipos),lensec-5,idrstmpl,ndpts, - & igdstmpl(1),igdstmpl(2),igdstmpl(3),fld) - else - print *,'gf_unpack7: Cannot use GDT 3.',igdsnum, - & ' to unpack Data Section 5.51.' - ierr=5 - nullify(fld) - return - endif - - elseif (idrsnum.eq.40 .OR. idrsnum.eq.40000) then - call jpcunpack(cgrib(ipos),lensec-5,idrstmpl,ndpts,fld) - - - elseif (idrsnum.eq.41 .OR. idrsnum.eq.40010) then - call pngunpack(cgrib(ipos),lensec-5,idrstmpl,ndpts,fld) - - else - print *,'gf_unpack7: Data Representation Template ',idrsnum, - & ' not yet implemented.' - ierr=4 - nullify(fld) - return - endif - - iofst=iofst+(8*lensec) - - end diff --git a/src/ixgb2.F90 b/src/ixgb2.F90 new file mode 100644 index 00000000..06c6d72e --- /dev/null +++ b/src/ixgb2.F90 @@ -0,0 +1,197 @@ +!> @file +!> @brief Generate an index record for each field in a GRIB2 message. +!> @author Mark Iredell @date 1995-10-31 + +!> Generate an index record for each field in a GRIB2 message. The index +!> records are written to index buffer pointed to by cbuf. All integers +!> in the index are in big-endian format. +!> +!> This subroutine is called by getg2ir(), which packages the index +!> records into an index file. +!> +!> The index buffer returned contains index records with the +!> format: +!> - byte 001 - 004 length of index record +!> - byte 005 - 008 bytes to skip in data file before GRIB message +!> - byte 009 - 012 bytes to skip in message before lus (local use) set = 0, if no local section. +!> - byte 013 - 016 bytes to skip in message before gds +!> - byte 017 - 020 bytes to skip in message before pds +!> - byte 021 - 024 bytes to skip in message before drs +!> - byte 025 - 028 bytes to skip in message before bms +!> - byte 029 - 032 bytes to skip in message before data section +!> - byte 033 - 040 bytes total in the message +!> - byte 041 - 041 GRIB version number (2) +!> - byte 042 - 042 message discipline +!> - byte 043 - 044 field number within GRIB2 message +!> - byte 045 - ii identification section (ids) +!> - byte ii + 1- jj grid definition section (gds) +!> - byte jj + 1- kk product definition section (pds) +!> - byte kk + 1- ll the data representation section (drs) +!> - byte ll + 1-ll + 6 first 6 bytes of the bit map section (bms) +!> +!> @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] lskip Number of bytes to skip before GRIB message. +!> @param[in] lgrib Number of bytes in GRIB message. When subroutine is +!> called, this must be set to the size of the cbuf buffer. +!> @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] numfld Number of index records created. +!> @param[out] mlen Total length of all index records. +!> @param[out] iret Return code +!> - 0 No error +!> - 1 Not enough memory available to hold full index buffer. +!> - 2 I/O error in read. +!> - 3 GRIB message is not edition 2. +!> - 4 Not enough memory to allocate extent to index buffer. +!> - 5 Unidentified GRIB section encountered. +!> +!> @author Mark Iredell @date 1995-10-31 +SUBROUTINE IXGB2(LUGB, LSKIP, LGRIB, CBUF, NUMFLD, MLEN, IRET) + USE RE_ALLOC ! NEEDED FOR SUBROUTINE REALLOC + implicit none + + CHARACTER(LEN = 1), POINTER, DIMENSION(:) :: CBUF + CHARACTER CVER, CDISC + CHARACTER(LEN = 4) :: CTEMP + INTEGER LOCLUS, LOCGDS, LENGDS, LOCBMS + integer :: indbmp, numsec, next, newsize, mova2i, mbuf, lindex + integer :: linmax, ixskp + integer :: mxspd, mxskp, mxsgd, mxsdr, mxsbm, mxlus + integer :: mxlen, mxds, mxfld, mxbms + integer :: init, ixlus, lugb, lskip, lgrib, numfld, mlen, iret + integer :: ixsgd, ibread, ibskip, ilndrs, ilnpds, istat, ixds + integer :: ixspd, ixfld, ixids, ixlen, ixsbm, ixsdr + integer :: lbread, lensec, lensec1 + PARAMETER(LINMAX = 5000, INIT = 50000, NEXT = 10000) + PARAMETER(IXSKP = 4, IXLUS = 8, IXSGD = 12, IXSPD = 16, IXSDR = 20, IXSBM = 24, & + IXDS = 28, IXLEN = 36, IXFLD = 42, IXIDS = 44) + PARAMETER(MXSKP = 4, MXLUS = 4, MXSGD = 4, MXSPD = 4, MXSDR = 4, MXSBM = 4, & + MXDS = 4, MXLEN = 4, MXFLD = 2, MXBMS = 6) + CHARACTER CBREAD(LINMAX), CINDEX(LINMAX) + CHARACTER CIDS(LINMAX), CGDS(LINMAX) + + LOCLUS = 0 + IRET = 0 + MLEN = 0 + NUMFLD = 0 + NULLIFY(CBUF) + MBUF = INIT + ALLOCATE(CBUF(MBUF), STAT = ISTAT) ! ALLOCATE INITIAL SPACE FOR CBUF + IF (ISTAT .NE. 0) THEN + IRET = 1 + RETURN + ENDIF + + ! Read sections 0 and 1 for versin number and discipline. + IBREAD = MIN(LGRIB, LINMAX) + CALL BAREAD(LUGB, LSKIP, IBREAD, LBREAD, CBREAD) + IF(LBREAD .NE. IBREAD) THEN + IRET = 2 + RETURN + ENDIF + IF(CBREAD(8) .NE. CHAR(2)) THEN ! NOT GRIB EDITION 2 + IRET = 3 + RETURN + ENDIF + CVER = CBREAD(8) + CDISC = CBREAD(7) + CALL G2_GBYTEC(CBREAD, LENSEC1, 16 * 8, 4 * 8) + LENSEC1 = MIN(LENSEC1, IBREAD) + CIDS(1:LENSEC1) = CBREAD(17:16 + LENSEC1) + IBSKIP = LSKIP + 16 + LENSEC1 + + ! Loop through remaining sections creating an index for each field. + IBREAD = MAX(5, MXBMS) + DO + CALL BAREAD(LUGB, IBSKIP, IBREAD, LBREAD, CBREAD) + CTEMP = CBREAD(1)//CBREAD(2)//CBREAD(3)//CBREAD(4) + IF (CTEMP .EQ. '7777') RETURN ! END OF MESSAGE FOUND + IF(LBREAD .NE. IBREAD) THEN + IRET = 2 + RETURN + ENDIF + CALL G2_GBYTEC(CBREAD, LENSEC, 0 * 8, 4 * 8) + CALL G2_GBYTEC(CBREAD, NUMSEC, 4 * 8, 1 * 8) + + IF (NUMSEC .EQ. 2) THEN ! SAVE LOCAL USE LOCATION + LOCLUS = IBSKIP-LSKIP + ELSEIF (NUMSEC .EQ. 3) THEN ! SAVE GDS INFO + LENGDS = LENSEC + CGDS = CHAR(0) + CALL BAREAD(LUGB, IBSKIP, LENGDS, LBREAD, CGDS) + IF (LBREAD .NE. LENGDS) THEN + IRET = 2 + RETURN + ENDIF + LOCGDS = IBSKIP-LSKIP + ELSEIF (NUMSEC .EQ. 4) THEN ! FOUND PDS + CINDEX = CHAR(0) + CALL G2_SBYTEC(CINDEX, LSKIP, 8 * IXSKP, 8 * MXSKP) ! BYTES TO SKIP + CALL G2_SBYTEC(CINDEX, LOCLUS, 8 * IXLUS, 8 * MXLUS) ! LOCATION OF LOCAL USE + CALL G2_SBYTEC(CINDEX, LOCGDS, 8 * IXSGD, 8 * MXSGD) ! LOCATION OF GDS + CALL G2_SBYTEC(CINDEX, IBSKIP-LSKIP, 8 * IXSPD, 8 * MXSPD) ! LOCATION OF PDS + CALL G2_SBYTEC(CINDEX, LGRIB, 8 * IXLEN, 8 * MXLEN) ! LEN OF GRIB2 + CINDEX(41) = CVER + CINDEX(42) = CDISC + CALL G2_SBYTEC(CINDEX, NUMFLD + 1, 8 * IXFLD, 8 * MXFLD) ! FIELD NUM + CINDEX(IXIDS + 1:IXIDS + LENSEC1) = CIDS(1:LENSEC1) + LINDEX = IXIDS + LENSEC1 + CINDEX(LINDEX + 1:LINDEX + LENGDS) = CGDS(1:LENGDS) + LINDEX = LINDEX + LENGDS + ILNPDS = LENSEC + CALL BAREAD(LUGB, IBSKIP, ILNPDS, LBREAD, CINDEX(LINDEX + 1)) + IF (LBREAD .NE. ILNPDS) THEN + IRET = 2 + RETURN + ENDIF + LINDEX = LINDEX + ILNPDS + ELSEIF (NUMSEC .EQ. 5) THEN ! FOUND DRS + CALL G2_SBYTEC(CINDEX, IBSKIP-LSKIP, 8 * IXSDR, 8 * MXSDR) ! LOCATION OF DRS + ILNDRS = LENSEC + CALL BAREAD(LUGB, IBSKIP, ILNDRS, LBREAD, CINDEX(LINDEX + 1)) + IF (LBREAD .NE. ILNDRS) THEN + IRET = 2 + RETURN + ENDIF + LINDEX = LINDEX + ILNDRS + ELSEIF (NUMSEC .EQ. 6) THEN ! FOUND BMS + INDBMP = MOVA2I(CBREAD(6)) + IF (INDBMP.LT.254) THEN + LOCBMS = IBSKIP-LSKIP + CALL G2_SBYTEC(CINDEX, LOCBMS, 8 * IXSBM, 8 * MXSBM) ! LOC. OF BMS + ELSEIF (INDBMP.EQ.254) THEN + CALL G2_SBYTEC(CINDEX, LOCBMS, 8 * IXSBM, 8 * MXSBM) ! LOC. OF BMS + ELSEIF (INDBMP.EQ.255) THEN + CALL G2_SBYTEC(CINDEX, IBSKIP-LSKIP, 8 * IXSBM, 8 * MXSBM) ! LOC. OF BMS + ENDIF + CINDEX(LINDEX + 1:LINDEX + MXBMS) = CBREAD(1:MXBMS) + LINDEX = LINDEX + MXBMS + CALL G2_SBYTEC(CINDEX, LINDEX, 0, 8 * 4) ! NUM BYTES IN INDEX RECORD + ELSEIF (NUMSEC .EQ. 7) THEN ! FOUND DATA SECTION + CALL G2_SBYTEC(CINDEX, IBSKIP-LSKIP, 8 * IXDS, 8 * MXDS) ! LOC. OF DATA SEC. + NUMFLD = NUMFLD + 1 + IF ((LINDEX + MLEN) .GT. MBUF) THEN ! ALLOCATE MORE SPACE IF NECESSARY + NEWSIZE = MAX(MBUF + NEXT, MBUF + LINDEX) + CALL REALLOC(CBUF, MLEN, NEWSIZE, ISTAT) + IF (ISTAT .NE. 0) THEN + NUMFLD = NUMFLD-1 + IRET = 4 + RETURN + ENDIF + MBUF = NEWSIZE + ENDIF + CBUF(MLEN + 1:MLEN + LINDEX) = CINDEX(1:LINDEX) + MLEN = MLEN + LINDEX + ELSE ! UNRECOGNIZED SECTION + IRET = 5 + RETURN + ENDIF + IBSKIP = IBSKIP + LENSEC + ENDDO +END SUBROUTINE IXGB2 diff --git a/src/ixgb2.f b/src/ixgb2.f deleted file mode 100644 index 1a780bb0..00000000 --- a/src/ixgb2.f +++ /dev/null @@ -1,201 +0,0 @@ -!> @file -!> @brief Generate an index record for each field in a GRIB2 message. -!> @author Mark Iredell @date 1995-10-31 - -!> Generate an index record for each field in a GRIB2 message. The index -!> records are written to index buffer pointed to by cbuf. All integers -!> in the index are in big-endian format. -!> -!> This subroutine is called by getg2ir(), which packages the index -!> records into an index file. -!> -!> The index buffer returned contains index records with the -!> format: -!> - byte 001 - 004 length of index record -!> - byte 005 - 008 bytes to skip in data file before GRIB message -!> - byte 009 - 012 bytes to skip in message before lus (local use) set = 0, if no local section. -!> - byte 013 - 016 bytes to skip in message before gds -!> - byte 017 - 020 bytes to skip in message before pds -!> - byte 021 - 024 bytes to skip in message before drs -!> - byte 025 - 028 bytes to skip in message before bms -!> - byte 029 - 032 bytes to skip in message before data section -!> - byte 033 - 040 bytes total in the message -!> - byte 041 - 041 GRIB version number (2) -!> - byte 042 - 042 message discipline -!> - byte 043 - 044 field number within GRIB2 message -!> - byte 045 - ii identification section (ids) -!> - byte ii+1- jj grid definition section (gds) -!> - byte jj+1- kk product definition section (pds) -!> - byte kk+1- ll the data representation section (drs) -!> - byte ll+1-ll+6 first 6 bytes of the bit map section (bms) -!> -!> @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] lskip Number of bytes to skip before GRIB message. -!> @param[in] lgrib Number of bytes in GRIB message. When subroutine is -!> called, this must be set to the size of the cbuf buffer. -!> @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] numfld Number of index records created. -!> @param[out] mlen Total length of all index records. -!> @param[out] iret Return code -!> - 0 No error -!> - 1 Not enough memory available to hold full index buffer. -!> - 2 I/O error in read. -!> - 3 GRIB message is not edition 2. -!> - 4 Not enough memory to allocate extent to index buffer. -!> - 5 Unidentified GRIB section encountered. -!> -!> @author Mark Iredell @date 1995-10-31 - SUBROUTINE IXGB2(LUGB,LSKIP,LGRIB,CBUF,NUMFLD,MLEN,IRET) - USE RE_ALLOC ! NEEDED FOR SUBROUTINE REALLOC - implicit none - - CHARACTER(LEN=1),POINTER,DIMENSION(:) :: CBUF - - CHARACTER CVER,CDISC - CHARACTER(LEN=4) :: CTEMP - INTEGER LOCLUS,LOCGDS,LENGDS,LOCBMS - integer :: indbmp, numsec, next, newsize, mova2i, mbuf, lindex - integer :: linmax, ixskp - integer :: mxspd, mxskp, mxsgd, mxsdr, mxsbm, mxlus - integer :: mxlen, mxds, mxfld, mxbms - integer :: init, ixlus, lugb, lskip, lgrib, numfld, mlen, iret - integer :: ixsgd, ibread, ibskip, ilndrs, ilnpds, istat, ixds - integer :: ixspd, ixfld, ixids, ixlen, ixsbm, ixsdr - integer :: lbread, lensec, lensec1 - PARAMETER(LINMAX=5000,INIT=50000,NEXT=10000) - PARAMETER(IXSKP=4,IXLUS=8,IXSGD=12,IXSPD=16,IXSDR=20,IXSBM=24, - & IXDS=28,IXLEN=36,IXFLD=42,IXIDS=44) - PARAMETER(MXSKP=4,MXLUS=4,MXSGD=4,MXSPD=4,MXSDR=4,MXSBM=4, - & MXDS=4,MXLEN=4,MXFLD=2,MXBMS=6) - - CHARACTER CBREAD(LINMAX),CINDEX(LINMAX) - CHARACTER CIDS(LINMAX),CGDS(LINMAX) - - LOCLUS=0 - IRET=0 - MLEN=0 - NUMFLD=0 - NULLIFY(CBUF) - MBUF=INIT - ALLOCATE(CBUF(MBUF),STAT=ISTAT) ! ALLOCATE INITIAL SPACE FOR CBUF - IF (ISTAT.NE.0) THEN - IRET=1 - RETURN - ENDIF - -C READ SECTIONS 0 AND 1 FOR VERSIN NUMBER AND DISCIPLINE - IBREAD=MIN(LGRIB,LINMAX) - CALL BAREAD(LUGB,LSKIP,IBREAD,LBREAD,CBREAD) - IF(LBREAD.NE.IBREAD) THEN - IRET=2 - RETURN - ENDIF - IF(CBREAD(8).NE.CHAR(2)) THEN ! NOT GRIB EDITION 2 - IRET=3 - RETURN - ENDIF - CVER=CBREAD(8) - CDISC=CBREAD(7) - CALL G2_GBYTEC(CBREAD,LENSEC1,16*8,4*8) - LENSEC1=MIN(LENSEC1,IBREAD) - CIDS(1:LENSEC1)=CBREAD(17:16+LENSEC1) - IBSKIP=LSKIP+16+LENSEC1 - -C LOOP THROUGH REMAINING SECTIONS CREATING AN INDEX FOR EACH FIELD - IBREAD=MAX(5,MXBMS) - DO - CALL BAREAD(LUGB,IBSKIP,IBREAD,LBREAD,CBREAD) - CTEMP=CBREAD(1)//CBREAD(2)//CBREAD(3)//CBREAD(4) - IF (CTEMP.EQ.'7777') RETURN ! END OF MESSAGE FOUND - IF(LBREAD.NE.IBREAD) THEN - IRET=2 - RETURN - ENDIF - CALL G2_GBYTEC(CBREAD,LENSEC,0*8,4*8) - CALL G2_GBYTEC(CBREAD,NUMSEC,4*8,1*8) - - IF (NUMSEC.EQ.2) THEN ! SAVE LOCAL USE LOCATION - LOCLUS=IBSKIP-LSKIP - ELSEIF (NUMSEC.EQ.3) THEN ! SAVE GDS INFO - LENGDS=LENSEC - CGDS=CHAR(0) - CALL BAREAD(LUGB,IBSKIP,LENGDS,LBREAD,CGDS) - IF(LBREAD.NE.LENGDS) THEN - IRET=2 - RETURN - ENDIF - LOCGDS=IBSKIP-LSKIP - ELSEIF (NUMSEC.EQ.4) THEN ! FOUND PDS - CINDEX=CHAR(0) - CALL G2_SBYTEC(CINDEX,LSKIP,8*IXSKP,8*MXSKP) ! BYTES TO SKIP - CALL G2_SBYTEC(CINDEX,LOCLUS,8*IXLUS,8*MXLUS) ! LOCATION OF LOCAL USE - CALL G2_SBYTEC(CINDEX,LOCGDS,8*IXSGD,8*MXSGD) ! LOCATION OF GDS - CALL G2_SBYTEC(CINDEX,IBSKIP-LSKIP,8*IXSPD,8*MXSPD) ! LOCATION OF PDS - CALL G2_SBYTEC(CINDEX,LGRIB,8*IXLEN,8*MXLEN) ! LEN OF GRIB2 - CINDEX(41)=CVER - CINDEX(42)=CDISC - CALL G2_SBYTEC(CINDEX,NUMFLD+1,8*IXFLD,8*MXFLD) ! FIELD NUM - CINDEX(IXIDS+1:IXIDS+LENSEC1)=CIDS(1:LENSEC1) - LINDEX=IXIDS+LENSEC1 - CINDEX(LINDEX+1:LINDEX+LENGDS)=CGDS(1:LENGDS) - LINDEX=LINDEX+LENGDS - ILNPDS=LENSEC - CALL BAREAD(LUGB,IBSKIP,ILNPDS,LBREAD,CINDEX(LINDEX+1)) - IF(LBREAD.NE.ILNPDS) THEN - IRET=2 - RETURN - ENDIF - ! CINDEX(LINDEX+1:LINDEX+ILNPDS)=CBREAD(1:ILNPDS) - LINDEX=LINDEX+ILNPDS - ELSEIF (NUMSEC.EQ.5) THEN ! FOUND DRS - CALL G2_SBYTEC(CINDEX,IBSKIP-LSKIP,8*IXSDR,8*MXSDR) ! LOCATION OF DRS - ILNDRS=LENSEC - CALL BAREAD(LUGB,IBSKIP,ILNDRS,LBREAD,CINDEX(LINDEX+1)) - IF(LBREAD.NE.ILNDRS) THEN - IRET=2 - RETURN - ENDIF - ! CINDEX(LINDEX+1:LINDEX+ILNDRS)=CBREAD(1:ILNDRS) - LINDEX=LINDEX+ILNDRS - ELSEIF (NUMSEC.EQ.6) THEN ! FOUND BMS - INDBMP=MOVA2I(CBREAD(6)) - IF ( INDBMP.LT.254 ) THEN - LOCBMS=IBSKIP-LSKIP - CALL G2_SBYTEC(CINDEX,LOCBMS,8*IXSBM,8*MXSBM) ! LOC. OF BMS - ELSEIF ( INDBMP.EQ.254 ) THEN - CALL G2_SBYTEC(CINDEX,LOCBMS,8*IXSBM,8*MXSBM) ! LOC. OF BMS - ELSEIF ( INDBMP.EQ.255 ) THEN - CALL G2_SBYTEC(CINDEX,IBSKIP-LSKIP,8*IXSBM,8*MXSBM) ! LOC. OF BMS - ENDIF - CINDEX(LINDEX+1:LINDEX+MXBMS)=CBREAD(1:MXBMS) - LINDEX=LINDEX+MXBMS - CALL G2_SBYTEC(CINDEX,LINDEX,0,8*4) ! NUM BYTES IN INDEX RECORD - ELSEIF (NUMSEC.EQ.7) THEN ! FOUND DATA SECTION - CALL G2_SBYTEC(CINDEX,IBSKIP-LSKIP,8*IXDS,8*MXDS) ! LOC. OF DATA SEC. - NUMFLD=NUMFLD+1 - IF ((LINDEX+MLEN).GT.MBUF) THEN ! ALLOCATE MORE SPACE IF NECESSARY - NEWSIZE=MAX(MBUF+NEXT,MBUF+LINDEX) - CALL REALLOC(CBUF,MLEN,NEWSIZE,ISTAT) - IF ( ISTAT .NE. 0 ) THEN - NUMFLD=NUMFLD-1 - IRET=4 - RETURN - ENDIF - MBUF=NEWSIZE - ENDIF - CBUF(MLEN+1:MLEN+LINDEX)=CINDEX(1:LINDEX) - MLEN=MLEN+LINDEX - ELSE ! UNRECOGNIZED SECTION - IRET=5 - RETURN - ENDIF - IBSKIP=IBSKIP+LENSEC - ENDDO - END