getgb1.f

Go to the documentation of this file.

C> @file
C> @brief Find and unpacks a grib message.
C> @author Mark Iredell @date 1994-04-01
 
C> Find and unpack a grib message.
C> Read an associated grib index file (unless it already was read).
C> Find in the index file a reference to the grib message requested.
C> The grib message request specifies the number of messages to skip
C> and the unpacked pds and gds parameters. (A requested parameter
C> of -1 means to allow any value of this parameter to be found.)
C> If the requested grib message is found, then it is read from the
C> grib file and unpacked. Its message number is returned along with
C> the unpacked pds and gds parameters, the unpacked bitmap (if any),
C> and the unpacked data. If the grib message is not found, then the
C> return code will be nonzero.
C>
C> Program history log:
C>  - Mark Iredell 1994-04-01
C>  - Ralph Jones 1995-05-10  Add one more parameter to getgb and
C> change name to getgb1.
C>
C> @param[in] lugb logical unit of the unblocked grib data file.
C> @param[in] lugi logical unit of the unblocked grib index file.
C> @param[in] jf integer maximum number of data points to unpack.
C> @param[in] j integer number of messages to skip (=0 to search from beginning)
C> (<0 to reopen index file and search from beginning).
C> @param[in] jpds integer (25) pds parameters for which to search
C> (=-1 for wildcard) look in doc block of w3fi63 for array kpds
C> for list of order of unpacked pds values.
C> In most cases you only need to set 4 or 5 values to pick up record.
C> @param[in] jgds integer (22) gds parameters for which to search
C> (only searched if jpds(3)=255) (=-1 for wildcard).
C> @param[out] grib Grib data array before it is unpacked.
C> @param[out] kf Integer number of data points unpacked.
C> @param[out] k Integer message number unpacked
C> (can be same as j in calling program
C> in order to facilitate multiple searches).
C> @param[out] kpds Integer (25) unpacked pds parameters.
C> @param[out] kgds Integer (22) unpacked gds parameters.
C> @param[out] lb Logical (kf) unpacked bitmap if present.
C> @param[out] f Real (kf) unpacked data.
C> @param[out] iret Integer return code.
C>  - 0 All ok.
C>  - 96 Error reading index file.
C>  - 97 Error reading grib file.
C>  - 98 Number of data points greater than jf.
C>  - 99 Request not found.
C>  - other w3fi63 grib unpacker return code.
C>
C> @author Mark Iredell @date 1994-04-01
      SUBROUTINE getgb1(LUGB,LUGI,JF,J,JPDS,JGDS,
     &                       GRIB,KF,K,KPDS,KGDS,LB,F,IRET)
C
      parameter(mbuf=8192*128)
      parameter(lpds=23,lgds=22)
C
      INTEGER      JPDS(25),JGDS(*),KPDS(25),KGDS(*)
      INTEGER      IPDSP(LPDS),JPDSP(LPDS),IGDSP(LGDS)
      INTEGER      JGDSP(LGDS)
      INTEGER      KPTR(20)
C
      LOGICAL      LB(*)
C
      REAL         F(*)
C
      CHARACTER    CBUF(MBUF)
      CHARACTER*81 CHEAD(2)
      CHARACTER*1  CPDS(28)
      CHARACTER*1  CGDS(42)
      CHARACTER*1  GRIB(*)
C
C     SAVE LUX,NSKP,NLEN,NNUM,CBUF
      SAVE
C
      DATA lux/0/
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  READ INDEX FILE
      IF(j.LT.0.OR.lugi.NE.lux) THEN
C        REWIND LUGI
C        READ(LUGI,fmt='(2A81)',IOSTAT=IOS) CHEAD
        CALL baread(lugi,0,162,ios,chead)
        IF(ios.EQ.162.AND.chead(1)(42:47).EQ.'GB1IX1') THEN
          lux=0
          READ(chead(2),'(8X,3I10,2X,A40)',iostat=ios) nskp,nlen,nnum
          IF(ios.EQ.0) THEN
            nbuf=nnum*nlen
            IF(nbuf.GT.mbuf) THEN
              print *,'GETGB1: INCREASE BUFFER FROM ',mbuf,' TO ',nbuf
              nnum=mbuf/nlen
              nbuf=nnum*nlen
            ENDIF
            CALL baread(lugi,nskp,nbuf,lbuf,cbuf)
            IF(lbuf.EQ.nbuf) THEN
              lux=lugi
              j=max(j,0)
            ENDIF
          ENDIF
        ENDIF
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  SEARCH FOR REQUEST
      lgrib=0
      kj=j
      k=j
      kf=0
      IF(j.GE.0.AND.lugi.EQ.lux) THEN
        lpdsp=0
        DO i=1,lpds
          IF(jpds(i).NE.-1) THEN
            lpdsp=lpdsp+1
            ipdsp(lpdsp)=i
            jpdsp(lpdsp)=jpds(i)
          ENDIF
        ENDDO
        lgdsp=0
        IF(jpds(3).EQ.255) THEN
          DO i=1,lgds
            IF(jgds(i).NE.-1) THEN
              lgdsp=lgdsp+1
              igdsp(lgdsp)=i
              jgdsp(lgdsp)=jgds(i)
            ENDIF
          ENDDO
        ENDIF
        iret=99
        dowhile(lgrib.EQ.0.AND.kj.LT.nnum)
          kj=kj+1
          lt=0
          IF(lpdsp.GT.0) THEN
            cpds=cbuf((kj-1)*nlen+26:(kj-1)*nlen+53)
            kptr=0
            CALL gbyte(cbuf,kptr(3),(kj-1)*nlen*8+25*8,3*8)
            CALL fi632(cpds,kptr,kpds,iret)
            DO i=1,lpdsp
              ip=ipdsp(i)
              lt=lt+abs(jpds(ip)-kpds(ip))
            ENDDO
          ENDIF
          IF(lt.EQ.0.AND.lgdsp.GT.0) THEN
            cgds=cbuf((kj-1)*nlen+54:(kj-1)*nlen+95)
            kptr=0
            CALL fi633(cgds,kptr,kgds,iret)
            DO i=1,lgdsp
              ip=igdsp(i)
              lt=lt+abs(jgds(ip)-kgds(ip))
            ENDDO
          ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  READ AND UNPACK GRIB DATA
          IF(lt.EQ.0) THEN
            CALL gbyte(cbuf,lskip,(kj-1)*nlen*8,4*8)
            CALL gbyte(cbuf,lgrib,(kj-1)*nlen*8+20*8,4*8)
            cgds=cbuf((kj-1)*nlen+54:(kj-1)*nlen+95)
            kptr=0
            CALL fi633(cgds,kptr,kgds,iret)
C  BSM      IF(LGRIB.LE.200+17*JF/8.AND.KGDS(2)*KGDS(3).LE.JF) THEN
C  Change number of bits that can be handled to 25
            IF(lgrib.LE.200+25*jf/8.AND.kgds(2)*kgds(3).LE.jf) THEN
              CALL baread(lugb,lskip,lgrib,lread,grib)
              IF(lread.EQ.lgrib) THEN
                CALL w3fi63(grib,kpds,kgds,lb,f,kptr,iret)
                IF(iret.EQ.0) THEN
                  k=kj
                  kf=kptr(10)
                ENDIF
              ELSE
                iret=97
              ENDIF
            ELSE
              iret=98
            ENDIF
          ENDIF
        ENDDO
      ELSE
        iret=96
      ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      RETURN
      END