w3nogds.f

Go to the documentation of this file.

C> @file
C                .      .    .                                       .
C> SUBPROGRAM:  W3NOGDS        MAKE A COMPLETE GRIB MESSAGE
C>   PRGMMR: FARLEY           ORG: NMC421      DATE:94-11-22
C>
C> ABSTRACT: MAKES A COMPLETE GRIB MESSAGE FROM A USER SUPPLIED
C>   ARRAY OF FLOATING POINT OR INTEGER DATA.  THE USER HAS THE
C>   OPTION OF SUPPLYING THE PDS OR AN INTEGER ARRAY THAT WILL BE
C>   USED TO CREATE A PDS (WITH W3FI68).  THE USER MUST ALSO
C>   SUPPLY OTHER NECESSARY INFO; SEE USAGE SECTION BELOW.
C>
C> PROGRAM HISTORY LOG:
C>   97-02-24  M. FARLEY    MODIFIED W3FI72 - this routine allows for
C>                          NO GDS (errors in W3FI71 for GRIB grids
C>                          21-26, 61-64 forced the need for this routine).
C>   98-06-24  Gilbert      Added number of gridpoint values for grids
C>                          61-64, needed when igflag=2 ( no gds ).
C>   98-12-21  Gilbert    Replaced Function ICHAR with mova2i.
C>
C> USAGE:  CALL W3NOGDS(ITYPE,FLD,IFLD,IBITL,
C>        &            IPFLAG,ID,PDS,
C>        &            IGFLAG,IGRID,IGDS,ICOMP,
C>        &            IBFLAG,IBMAP,IBLEN,IBDSFL,
C>        &            IBDSFL,
C>        &            NPTS,KBUF,ITOT,JERR)
C>
C>   INPUT ARGUMENT LIST:
C>     ITYPE    - 0 = FLOATING POINT DATA SUPPLIED IN ARRAY 'FLD'
C>                1 = INTEGER DATA SUPPLIED IN ARRAY 'IFLD'
C>     FLD      - REAL ARRAY OF DATA (AT PROPER GRIDPOINTS) TO BE
C>                CONVERTED TO GRIB FORMAT IF ITYPE=0.
C>                SEE REMARKS #1 & 2.
C>     IFLD     - INTEGER ARRAY OF DATA (AT PROPER GRIDPOINTS) TO BE
C>                CONVERTED TO GRIB FORMAT IF ITYPE=1.
C>                SEE REMARKS #1 & 2.
C>     IBITL    - 0 = COMPUTER COMPUTES LENGTH FOR PACKING DATA FROM
C>                    POWER OF 2 (NUMBER OF BITS) BEST FIT OF DATA
C>                    USING 'VARIABLE' BIT PACKER W3FI58.
C>                8, 12, ETC. COMPUTER RESCALES DATA TO FIT INTO THAT
C>                    'FIXED' NUMBER OF BITS USING W3FI59.
C>                SEE REMARKS #3.
C>
C>     IPFLAG   - 0 = MAKE PDS FROM USER SUPPLIED ARRAY (ID)
C>                1 = USER SUPPLYING PDS
C>                NOTE: IF PDS IS GREATER THAN 30, USE IPLFAG=1.
C>                THE USER COULD CALL W3FI68 BEFORE HE CALLS
C>                W3NOGDS. THIS WOULD MAKE THE FIRST 30 BYTES OF
C>                THE PDS, USER THEN WOULD MAKE BYTES AFTER 30.
C>     ID       - INTEGER ARRAY OF  VALUES THAT W3FI68 WILL USE
C>                TO MAKE AN EDITION 1 PDS IF IPFLAG=0.  (SEE THE
C>                DOCBLOCK FOR W3FI68 FOR LAYOUT OF ARRAY)
C>     PDS      - CHARACTER ARRAY OF VALUES (VALID PDS SUPPLIED
C>                BY USER) IF IPFLAG=1. LENGTH MAY EXCEED 28 BYTES
C>                (CONTENTS OF BYTES BEYOND 28 ARE PASSED
C>                THROUGH UNCHANGED).
C>
C>     IGFLAG   - 0 = MAKE GDS BASED ON 'IGRID' VALUE.
C>                1 = MAKE GDS FROM USER SUPPLIED INFO IN 'IGDS'
C>                    AND 'IGRID' VALUE.
C>                SEE REMARKS #4.
C>                2 = NO GDS WILL BE INCLUDED...for international grids
C>                *** THIS IS AN EXCEPTION TO REMARKS #4!!!!
C>     IGRID    - #   = GRID IDENTIFICATION (TABLE B)
C>                255 = IF USER DEFINED GRID; IGDS MUST BE SUPPLIED
C>                      AND IGFLAG MUST =1.
C>     IGDS     - INTEGER ARRAY CONTAINING USER GDS INFO (SAME
C>                FORMAT AS SUPPLIED BY W3FI71 - SEE DOCKBLOCK FOR
C>                LAYOUT) IF IGFLAG=1.
C>     ICOMP    - RESOLUTION AND COMPONENT FLAG FOR BIT 5 OF GDS(17)
C>                0 = EARTH ORIENTED WINDS
C>                1 = GRID ORIENTED WINDS
C>
C>     IBFLAG   - 0 = MAKE BIT MAP FROM USER SUPPLIED DATA
C>                # = BIT MAP PREDEFINED BY CENTER
C>                SEE REMARKS #5.
C>     IBMAP    - INTEGER ARRAY CONTAINING BIT MAP
C>     IBLEN    - LENGTH OF BIT MAP WILL BE USED TO VERIFY LENGTH
C>                OF FIELD (ERROR IF IT DOESN'T MATCH).
C>
C>     IBDSFL   - INTEGER ARRAY CONTAINING TABLE 11 FLAG INFO
C>                BDS OCTET 4:
C>                (1) 0 = GRID POINT DATA
C>                    1 = SPHERICAL HARMONIC COEFFICIENTS
C>                (2) 0 = SIMPLE PACKING
C>                    1 = SECOND ORDER PACKING
C>                (3) ... SAME VALUE AS 'ITYPE'
C>                    0 = ORIGINAL DATA WERE FLOATING POINT VALUES
C>                    1 = ORIGINAL DATA WERE INTEGER VALUES
C>                (4) 0 = NO ADDITIONAL FLAGS AT OCTET 14
C>                    1 = OCTET 14 CONTAINS FLAG BITS 5-12
C>                (5) 0 = RESERVED - ALWAYS SET TO 0
C>         BYTE 6 OPTION 1 NOT AVAILABLE (AS OF 5-16-93)
C>                (6) 0 = SINGLE DATUM AT EACH GRID POINT
C>                    1 = MATRIX OF VALUES AT EACH GRID POINT
C>         BYTE 7 OPTION 0 WITH SECOND ORDER PACKING N/A (AS OF 5-16-93)
C>                (7) 0 = NO SECONDARY BIT MAPS
C>                    1 = SECONDARY BIT MAPS PRESENT
C>                (8) 0 = SECOND ORDER VALUES HAVE CONSTANT WIDTH
C>                    1 = SECOND ORDER VALUES HAVE DIFFERENT WIDTHS
C>
C>   OUTPUT ARGUMENT LIST:
C>     NPTS     - NUMBER OF GRIDPOINTS IN ARRAY FLD OR IFLD
C>     KBUF     - ENTIRE GRIB MESSAGE ('GRIB' TO '7777')
C>                EQUIVALENCE TO INTEGER ARRAY TO MAKE SURE IT
C>                IS ON WORD BOUNARY.
C>     ITOT     - TOTAL LENGTH OF GRIB MESSAGE IN BYTES
C>     JERR     - = 0, COMPLETED MAKING GRIB FIELD WITHOUT ERROR
C>                  1, IPFLAG NOT 0 OR 1
C>                  2, IGFLAG NOT 0 OR 1 OR 2
C>                  3, ERROR CONVERTING IEEE F.P. NUMBER TO IBM370 F.P.
C>                  4, W3FI71 ERROR/IGRID NOT DEFINED
C>                  5, W3FK74 ERROR/GRID REPRESENTATION TYPE NOT VALID
C>                  6, GRID TOO LARGE FOR PACKER DIMENSION ARRAYS
C>                     SEE AUTOMATION DIVISION FOR REVISION!
C>                  7, LENGTH OF BIT MAP NOT EQUAL TO SIZE OF FLD/IFLD
C>                  8, W3FI73 ERROR, ALL VALUES IN IBMAP ARE ZERO
C>
C>   OUTPUT FILES:
C>     FT06F001 - STANDARD FORTRAN OUTPUT PRINT FILE
C>
C>   SUBPROGRAMS CALLED:
C>     LIBRARY:
C>       W3LIB    - W3FI58, W3FI59, W3FI68, W3FI71, W3FI73, W3FI74
C>                  W3FI75, W3FI76, W3FI01
C>
C> REMARKS:
C>   1)  IF BIT MAP TO BE INCLUDED IN MESSAGE, NULL DATA SHOULD
C>       BE INCLUDED IN FLD OR IFLD.  THIS ROUTINE WILL TAKE CARE
C>       OF 'DISCARDING' ANY NULL DATA BASED ON THE BIT MAP.
C>   2)  UNITS MUST BE THOSE IN GRIB DOCUMENTATION:  NMC O.N. 388
C>       OR WMO PUBLICATION 306.
C>   3)  IN EITHER CASE, INPUT NUMBERS WILL BE MULTIPLIED BY
C>       '10 TO THE NTH' POWER FOUND IN ID(25) OR PDS(27-28),
C>       THE D-SCALING FACTOR, PRIOR TO BINARY PACKING.
C>   4)  ALL NMC PRODUCED GRIB FIELDS WILL HAVE A GRID DEFINITION
C>       SECTION INCLUDED IN THE GRIB MESSAGE.  ID(6) WILL BE
C>       SET TO '1'.
C>       - GDS WILL BE BUILT BASED ON GRID NUMBER (IGRID), UNLESS
C>         IGFLAG=1 (USER SUPPLYING IGDS).  USER MUST STILL SUPPLY
C>         IGRID EVEN IF IGDS PROVIDED.
C>   5)  IF BIT MAP USED THEN ID(7) OR PDS(8) MUST INDICATE THE
C>       PRESENCE OF A BIT MAP.
C>   6)  ARRAY KBUF SHOULD BE EQUIVALENCED TO AN INTEGER VALUE OR
C>       ARRAY TO MAKE SURE IT IS ON A WORD BOUNDARY.
C>   7)  SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C>
C> ATTRIBUTES:
C>   LANGUAGE: CRAY CFT77 FORTRAN
C>   MACHINE:  CRAY C916/256, Y-MP8/864, Y-MP EL92/256, J916/2048
C>
      SUBROUTINE w3nogds(ITYPE,FLD,IFLD,IBITL,
     &                  IPFLAG,ID,PDS,
     &                  IGFLAG,IGRID,IGDS,ICOMP,
     &                  IBFLAG,IBMAP,IBLEN,IBDSFL,
     &                  NPTS,KBUF,ITOT,JERR)
C
      parameter(mxsize=260000)
C     ALLOW UP TO 24 BITS PER POINT
      parameter(mxsiz3=mxsize*3)
      parameter(mxsizb=mxsize/8+6)
C     FOR 64 BIT CRAY
      parameter(mxsizi=mxsiz3/8)
C     FOR 32 BIT WORKSTATIONS AND HDS
C     PARAMETER       (MXSIZI=MXSIZ3/4)
C
      REAL            FLD(*)
C
      INTEGER         IBDSFL(*)
      INTEGER         IBMAP(*)
      INTEGER         ID(*)
      INTEGER         IFLD(*)
      INTEGER         IGDS(*)
      INTEGER         IPFLD(MXSIZI)
      INTEGER         IB(4)
C
      CHARACTER * 1   BDS11(11)
      CHARACTER * 1   KBUF(*)
      CHARACTER * 1   PDS(*)
      CHARACTER * 1   GDS(200)
      CHARACTER * 1   BMS(MXSIZB)
      CHARACTER * 1   PFLD(MXSIZ3)
      CHARACTER * 1   SEVEN
      CHARACTER * 1   ZERO
C
      equivalence(ipfld(1),pfld(1))
      equivalence(bds11(1),idummy)
C
C   ASCII REP OF  /'G', 'R', 'I', 'B'/
C
      DATA  ib    / 71,  82,  73,  66/
C
      ier    = 0
      iberr  = 0
      jerr   = 0
      igribl = 8
      ipdsl  = 0
      lengds = 0
      lenbms = 0
      lenbds = 0
      itoss  = 0
C
C$           1.0   PRODUCT DEFINITION SECTION(PDS).
C
C   SET ID(6) TO 1 ...OR... MODIFY PDS(8) ...
C      REGARDLESS OF USER SPECIFICATION...
C   NMC GRIB FIELDS WILL ALWAYS HAVE A GDS
C ***** exception for international GRIB GRIDS 21-26, 61-64
C ***** which will NOT contain a GDS until subroutine W3FI71 is fixed!
C
      IF (ipflag .EQ.0) THEN
        id(6) = 1
      if (igflag .eq. 2) then
        id(6) = 0
      endif
        CALL w3fi68(id,pds)
      ELSE IF (ipflag .EQ. 1) THEN
        IF (iand(mova2i(pds(8)),64) .EQ. 64) THEN
C         BOTH GDS AND BMS
          pds(8) = char(192)
        ELSE IF (mova2i(pds(8)) .EQ. 0) THEN
C         GDS ONLY
          pds(8) = char(128)
        END IF
        CONTINUE
      ELSE
C       PRINT *,' W3NOGDS ERROR, IPFLAG IS NOT 0 OR 1 IPFLAG = ',IPFLAG
        jerr = 1
        GO TO 900
      END IF
C
C     GET LENGTH OF PDS
C
      ipdsl = mova2i(pds(1)) * 65536 + mova2i(pds(2)) * 256 +
     &        mova2i(pds(3))
C
C$           2.0   GRID DEFINITION SECTION (GDS).
C
C     IF IGFLAG=1 THEN USER IS SUPPLYING THE IGDS INFORMATION
C     IF IGFLAG=2 THEN USER doesn't want a GDS and this section
C                will be skipped...LENGDS=0
C
      IF (igflag .EQ. 0) THEN
        CALL w3fi71(igrid,igds,igerr)
        IF (igerr .EQ. 1) THEN
C         PRINT *,' W3FI71 ERROR, GRID TYPE NOT DEFINED...',IGRID
          jerr = 4
          GO TO 900
        END IF
      END IF
      IF (igflag .EQ. 0  .OR.  igflag .EQ.1) THEN
        CALL w3fi74(igds,icomp,gds,lengds,npts,igerr)
        IF (igerr .EQ. 1) THEN
C         PRINT *,' W3FI74 ERROR, GRID REP TYPE NOT VALID...',IGDS(3)
          jerr = 5
          GO TO 900
        ELSE
        END IF
        IF (npts .GT. mxsize) THEN
C         PRINT *,' W3NOGDS ERROR, GRID TOO LARGE FOR PACKER ARRAY',
C    &            ' DIMENSIONS'
          jerr = 6
          GO TO 900
        END IF
      else if (igflag .eq. 2) then
    lengds = 0
        if (igrid.eq.21) then
      npts=1333
        else if (igrid.eq.22) then
      npts=1333
        else if (igrid.eq.23) then
      npts=1333
        else if (igrid.eq.24) then
      npts=1333
    else if (igrid.eq.25) then
      npts=1297
    else if (igrid.eq.26) then
      npts=1297
    else if ((igrid.ge.61).and.(igrid.le.64)) then
      npts=4096
    end if 
      ELSE
C       PRINT *,' W3NOGDS ERROR, IGFLAG IS NOT 0-2 IGFLAG = ',IGFLAG
        GO TO 900
      END IF
C
C$           3.0   BIT MAP SECTION (BMS).
C
C     SET ITOSS=1 IF BITMAP BEING USED.  W3FI75 WILL TOSS DATA
C     PRIOR TO PACKING.  LATER CODING WILL BE NEEDED WHEN THE
C     'PREDEFINED' GRIDS ARE FINALLY 'DEFINED'.
C
      IF (mova2i(pds(8)) .EQ. 64 .OR.
     &    mova2i(pds(8)) .EQ. 192)   THEN
        itoss = 1
        IF (ibflag .EQ. 0) THEN
          IF (iblen .NE. npts) THEN
C           PRINT *,' W3NOGDS ERROR, IBLEN .NE. NPTS = ',IBLEN,NPTS
            jerr = 7
            GO TO 900
          END IF
          CALL w3fi73(ibflag,ibmap,iblen,bms,lenbms,ier)
          IF (ier .NE. 0) THEN
C           PRINT *,' W3FI73 ERROR, IBMAP VALUES ARE ALL ZERO'
            jerr = 8
            GO TO 900
          END IF
        ELSE
C         PRINT *,'   BIT MAP PREDEFINED BY CENTER, IBFLAG = ',IBFLAG
        END IF
      END IF
C
C$           4.0   BINARY DATA SECTION (BDS).
C
C$           4.1   SCALE THE DATA WITH D-SCALE FROM PDS(27-28)
C
      jscale = mova2i(pds(27)) * 256 + mova2i(pds(28))
      IF (iand(jscale,32768).NE.0) THEN
        jscale = - iand(jscale,32767)
      END IF
      scale  = 10.0 ** jscale
      IF (itype .EQ. 0) THEN
        DO 410 i = 1,npts
          fld(i) = fld(i) * scale
  410   CONTINUE
      ELSE
        DO 411 i = 1,npts
          ifld(i) = nint(float(ifld(i)) * scale)
  411   CONTINUE
      END IF
C
C$           4.2   CALL W3FI75 TO PACK DATA AND MAKE BDS.
C
      CALL w3fi75(ibitl,itype,itoss,fld,ifld,ibmap,ibdsfl,
     &         npts,bds11,ipfld,pfld,len,lenbds,iberr,pds,igds)
        IF (iberr .EQ. 1) THEN
          jerr = 3
          GO TO 900
        END IF
C            4.3   IF D-SCALE NOT 0, RESCALE INPUT FIELD TO
C                   ORIGINAL VALUE
C
      IF (jscale.NE.0) THEN
          dscale = 1.0 / scale
          IF (itype.EQ.0) THEN
              DO 412 i = 1, npts
                  fld(i)  = fld(i) * dscale
  412         CONTINUE
          ELSE
              DO 413 i = 1, npts
                  fld(i)  = nint(float(ifld(i)) * dscale)
  413         CONTINUE
          END IF
      END IF
C
C$           5.0   OUTPUT SECTION.
C
C$           5.1   ZERO OUT THE OUTPUT ARRAY KBUF.
C
      zero    = char(00)
      itot    = igribl + ipdsl + lengds + lenbms + lenbds + 4
C     PRINT *,'IGRIBL  =',IGRIBL
C     PRINT *,'IPDSL   =',IPDSL
C     PRINT *,'LENGDS  =',LENGDS
C     PRINT *,'LENBMS  =',LENBMS
C     PRINT *,'LENBDS  =',LENBDS
C     PRINT *,'ITOT    =',ITOT
C
C     KBUF MUST BE ON A WORD BOUNDRY, EQUIVALENCE TO AN
C     INTEGER ARRAY IN THE MAIN PROGRAM TO MAKE SURE IT IS.
C     THIS IS BOTH COMPUTER AND COMPILER DEPENDENT, W3FI01
C     IS USED TO FILL OUT IF THE COMPUTER IS A 64 BIT OR
C     32 BIT WORD SIZE COMPUTER. LW IS SET TO 4 FOR 32 BIT
C     COMPUTER, 8 FOR 64 BIT COMPUTER.
C
      CALL w3fi01(lw)
      iwords = itot / lw
      CALL xstore(kbuf,0,iwords)
      IF (mod(itot,lw).NE.0) THEN
        ibytes = itot - iwords * lw
        DO 510 i = 1,ibytes
          kbuf(iwords * lw + i) = zero
  510   CONTINUE
      END IF
C
C$           5.2   MOVE SECTION 0 - 'IS' INTO KBUF (8 BYTES).
C
      istart  = 0
      DO 520 i = 1,4
        kbuf(i) = char(ib(i))
  520 CONTINUE
C
      kbuf(5) = char(mod(itot / 65536,256))
      kbuf(6) = char(mod(itot /   256,256))
      kbuf(7) = char(mod(itot        ,256))
      kbuf(8) = char(1)
C
C$           5.3   MOVE SECTION 1 - 'PDS' INTO KBUF (28 BYTES).
C
      istart  = istart + igribl
      IF (ipdsl.GT.0) THEN
        CALL xmovex(kbuf(istart+1),pds,ipdsl)
      ELSE
C       PRINT *,'LENGTH OF PDS LESS OR EQUAL 0, IPDSL = ',IPDSL
      END IF
C
C$           5.4   MOVE SECTION 2 - 'GDS' INTO KBUF.
C
      istart  = istart + ipdsl
      IF (lengds .GT. 0) THEN
        CALL xmovex(kbuf(istart+1),gds,lengds)
      END IF
C
C$           5.5   MOVE SECTION 3 - 'BMS' INTO KBUF.
C
      istart  = istart + lengds
      IF (lenbms .GT. 0) THEN
        CALL xmovex(kbuf(istart+1),bms,lenbms)
      END IF
C
C$           5.6   MOVE SECTION 4 - 'BDS' INTO KBUF.
C
C$                 MOVE THE FIRST 11 OCTETS OF THE BDS INTO KBUF.
C
      istart  = istart + lenbms
      CALL xmovex(kbuf(istart+1),bds11,11)
C
C$                 MOVE THE PACKED DATA INTO THE KBUF
C
      istart  = istart + 11
      IF (len.GT.0) THEN
        CALL xmovex(kbuf(istart+1),pfld,len)
      END IF
C
C$                 ADD '7777' TO END OFF KBUF
C   NOTE THAT THESE 4 OCTETS NOT INCLUDED IN ACTUAL SIZE OF BDS.
C
      seven  = char(55)
      istart = itot - 4
      DO 562 i = 1,4
        kbuf(istart+i) = seven
 562  CONTINUE
C
 900  CONTINUE
      RETURN
      END