w3fi75.f

Go to the documentation of this file.

C> @file
C> @brief GRIB pack data and form bds octets(1-11)
C> @author M. Farley @date 1992-07-10
 
C> This routine packs a grib field and forms octets(1-11)
C> of the binary data section (bds).
C>
C> Program history log:
C> - M. Farley 1992-07-10 Original author
C> - Ralph Jones 1992-10-01 Correction for field of constant data
C> - Ralph Jones 1992-10-16 Get rid of arrays fp and int
C> - Bill Cavanaugh 1993-08-06 Added routines fi7501, fi7502, fi7503
C> To allow second order packing in pds.
C> - John Stackpole 1993-07-21 Assorted repairs to get 2nd diff pack in
C> - Bill Cavanaugh 1993-10-28 Commented out nonoperational prints and
C> Write statements
C> - Bill Cavanaugh 1993-12-15 Corrected location of start of first order
C> Values and start of second order values to
C> Reflect a byte location in the bds instead
C> Of an offset in subroutine fi7501().
C> - Bill Cavanaugh 1994-01-27 Added igds as input argument to this routine
C> And added pds and igds arrays to the call to
C> W3fi82 to provide information needed for
C> Boustrophedonic processing.
C> - Bill Cavanaugh 1994-05-25 Subroutine fi7503 has been added to provide
C> For row by row or column by column second
C> Order packing.  this feature can be activated
C> By setting ibdsfl(7) to zero.
C> - Bill Cavanaugh 1994-07-08 Commented out print statements used for debug
C> - M. Farley 1994-11-22 Enlarged work arrays to handle .5degree grids
C> - Ralph Jones 1995-06-01 Correction for number of unused bits at end
C> Of section 4, in bds byte 4, bits 5-8.
C> - Mark Iredell 1995-10-31 Removed saves and prints
C> - Stephen Gilbert 2001-06-06 Changed gbyte/sbyte calls to refer to
C> Wesley ebisuzaki's endian independent
C> versions gbytec/sbytec.
C> Use f90 standard routine bit_size to get
C> number of bits in an integer instead of w3fi01.
C>
C> @param[in] IBITL
C> - 0, computer computes packing length from power of 2 that best fits the data.
C> - 8, 12, etc. computer rescales data to fit into set number of bits.
C> @param[in] ITYPE
C> - 0 = if input data is floating point (fld)
C> - 1 = If input data is integer (ifld)
C> @param[in] ITOSS
C> - 0 = no bit map is included (don't toss data)
C> - 1 = Toss null data according to ibmap
C> @param[in] FLD Real array of data to be packed if itype=0
C> @param[in] IFLD Integer array to be packed if itype=1
C> @param[in] IBMAP Bit map supplied from user
C> @param[in] IBDSFL Integer array containing table 11 flag info
C> BDS octet 4:
C> - (1)
C>  - 0 = grid point data
C>  - 1 = spherical harmonic coefficients
C> - (2)
C>  - 0 = simple packing
C>  - 1 = second order packing
C> - (3)
C>  - 0 = original data were floating point values
C>  - 1 = original data were integer values
C> - (4)
C>  - 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> - (6)
C>  - 0 = single datum at each grid point
C>  - 1 = matrix of values at each grid point
C> - (7)
C>  - 0 = no secondary bit maps
C>  - 1 = secondary bit maps present
C> - (8)
C>  - 0 = second order values have constant width
C>  - 1 = second order values have different widths
C> @param[in] NPTS Number of gridpoints in array to be packed
C> @param[in] IGDS Array of gds information
C> @param[out] BDS11 First 11 octets of bds
C> @param[out] PFLD Packed grib field
C> @param[out] LEN Length of pfld
C> @param[out] LENBDS Length of bds
C> @param[out] IBERR 1, error converting ieee f.p. number to ibm370 f.p.
C> @param IPFLD
C> @param PDS
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
      SUBROUTINE w3fi75 (IBITL,ITYPE,ITOSS,FLD,IFLD,IBMAP,IBDSFL,
     &  NPTS,BDS11,IPFLD,PFLD,LEN,LENBDS,IBERR,PDS,IGDS)
C
      REAL            FLD(*)
C     REAL            FWORK(260000)
C
C     FWORK CAN USE DYNAMIC ALLOCATION OF MEMORY ON CRAY
C
      REAL            FWORK(NPTS)
      REAL            RMIN,REFNCE
C
      character(len=1)         IPFLD(*)
      INTEGER         IBDSFL(*)
      INTEGER         IBMAP(*)
      INTEGER         IFLD(*),IGDS(*)
C     INTEGER         IWORK(260000)
C
C     IWORK CAN USE DYNAMIC ALLOCATION OF MEMORY ON CRAY
C
      INTEGER         IWORK(NPTS)
C
      LOGICAL         CONST
C
      CHARACTER * 1   BDS11(11),PDS(*)
      CHARACTER * 1   PFLD(*)
C
C            1.0   PACK THE FIELD.
C
C            1.1   TOSS DATA IF BITMAP BEING USED,
C                  MOVING 'DATA' TO WORK AREA...
C
      const = .false.
      iberr = 0
      iw    = 0
C
      IF (itoss .EQ. 1) THEN
        IF (itype .EQ. 0) THEN
          DO 110 it=1,npts
            IF (ibmap(it) .EQ. 1) THEN
              iw = iw + 1
              fwork(iw) = fld(it)
            ENDIF
  110     CONTINUE
          npts = iw
        ELSE IF (itype .EQ. 1) THEN
          DO 111 it=1,npts
            IF (ibmap(it) .EQ. 1) THEN
              iw = iw + 1
              iwork(iw) = ifld(it)
            ENDIF
  111     CONTINUE
          npts = iw
        ENDIF
C
C             ELSE, JUST MOVE DATA TO WORK ARRAY
C
      ELSE IF (itoss .EQ. 0) THEN
        IF (itype .EQ. 0) THEN
          DO 112 it=1,npts
            fwork(it) = fld(it)
  112     CONTINUE
        ELSE IF (itype .EQ. 1) THEN
          DO 113 it=1,npts
            iwork(it) = ifld(it)
  113     CONTINUE
        ENDIF
      ENDIF
C
C            1.2   CONVERT DATA IF NEEDED PRIOR TO PACKING.
C                  (INTEGER TO F.P. OR F.P. TO INTEGER)
C     ITYPE = 0...FLOATING POINT DATA
C       IBITL = 0...PACK IN LEAST # BITS...CONVERT TO INTEGER
C     ITYPE = 1...INTEGER DATA
C       IBITL > 0...PACK IN FIXED # BITS...CONVERT TO FLOATING POINT
C
      IF (itype .EQ. 0 .AND. ibitl .EQ. 0) THEN
        DO 120 if=1,npts
          iwork(if) = nint(fwork(if))
  120   CONTINUE
      ELSE IF (itype .EQ. 1 .AND. ibitl .NE. 0) THEN
        DO 123 if=1,npts
          fwork(if) = float(iwork(if))
  123   CONTINUE
      ENDIF
C
C            1.3   PACK THE DATA.
C
      IF (ibdsfl(2).NE.0) THEN
C                                    SECOND ORDER PACKING
C
C            PRINT*,'  DOING SECOND ORDER PACKING...'
          IF (ibitl.EQ.0) THEN
C
C             PRINT*,'    AND VARIABLE BIT PACKING'
C
C                           WORKING WITH INTEGER VALUES
C                           SINCE DOING VARIABLE BIT PACKING
C
              max  = iwork(1)
              min  = iwork(1)
              DO 300 i = 2, npts
                  IF (iwork(i).LT.min) THEN
                      min  = iwork(i)
                  ELSE IF (iwork(i).GT.max) THEN
                      max  = iwork(i)
                  END IF
  300         CONTINUE
C                           EXTRACT MINIMA
              DO 400 i = 1, npts
C                 IF (IWORK(I).LT.0) THEN
C                     PRINT *,'MINIMA 400',I,IWORK(I),NPTS
C                 END IF
                  iwork(i)  = iwork(i) - min
  400         CONTINUE
              refnce  = min
              idiff   = max - min
C             PRINT *,'REFERENCE VALUE',REFNCE
C
C             WRITE (6,FMT='(''  MINIMA REMOVED      = '',/,
C    &              10(3X,10I10,/))') (IWORK(I),I=1,6)
C             WRITE (6,FMT='(''  END OF ARRAY  = '',/,
C    &              10(3X,10I10,/))') (IWORK(I),I=NPTS-5,NPTS)
C
C                      FIND BIT WIDTH OF IDIFF
C
              CALL fi7505 (idiff,kwide)
C             PRINT*,'  BIT WIDTH FOR ORIGINAL DATA', KWIDE
              iscal2 = 0
C
C             MULTIPLICATIVE SCALE FACTOR SET TO 1
C             IN ANTICIPATION OF POSSIBLE USE IN GLAHN 2DN DIFF
C
              scal2 = 1.
C
          ELSE
C
C             PRINT*,'   AND FIXED BIT PACKING, IBITL = ', IBITL
C                               FIXED BIT PACKING
C                               - LENGTH OF FIELD IN IBITL
C                               - MUST BE REAL DATA
C                            FLOATING POINT INPUT
C
              rmax  = fwork(1)
              rmin  = fwork(1)
              DO 100 i = 2, npts
                  IF (fwork(i).LT.rmin) THEN
                      rmin  = fwork(i)
                  ELSE IF (fwork(i).GT.rmax) THEN
                      rmax  = fwork(i)
                  END IF
  100         CONTINUE
              refnce  = rmin
C             PRINT *,'100 REFERENCE',REFNCE
C                             EXTRACT MINIMA
              DO 200 i = 1, npts
                  fwork(i)  = fwork(i) - rmin
  200         CONTINUE
C             PRINT *,'REFERENCE VALUE',REFNCE
C             WRITE (6,FMT='(''  MINIMA REMOVED      = '',/,
C    &              10(3X,10F8.2,/))') (FWORK(I),I=1,6)
C             WRITE (6,FMT='(''  END OF ARRAY  = '',/,
C    &              10(3X,10F8.2,/))') (FWORK(I),I=NPTS-5,NPTS)
C                                FIND LARGEST DELTA
              idelt  = nint(rmax - rmin)
C                                DO BINARY SCALING
C                                   FIND OUT WHAT BINARY SCALE FACTOR
C                                       PERMITS CONTAINMENT OF
C                                       LARGEST DELTA
              CALL fi7505 (idelt,iwide)
C
C                                   BINARY SCALING
C
              iscal2  = iwide - ibitl
C             PRINT *,'SCALING NEEDED TO FIT =',ISCAL2
C             PRINT*,'  RANGE OF  = ',IDELT
C
C                                EXPAND DATA WITH BINARY SCALING
C                                CONVERT TO INTEGER
              scal2  = 2.0**iscal2
              scal2  = 1./ scal2
              DO 600 i = 1, npts
                  iwork(i)  = nint(fwork(i) * scal2)
  600         CONTINUE
              kwide = ibitl
          END IF
C
C  *****************************************************************
C
C           FOLLOWING IS FOR GLAHN SECOND DIFFERENCING
C           NOT STANDARD GRIB
C
C            TEST FOR SECOND DIFFERENCE PACKING
C            BASED OF SIZE OF PDS - SIZE IN FIRST 3 BYTES
C
          CALL gbytec(pds,ipdsiz,0,24)
          IF (ipdsiz.EQ.50) THEN
C             PRINT*,'  DO SECOND DIFFERENCE PACKING '
C
C                   GLAHN PACKING TO 2ND DIFFS
C
C             WRITE (6,FMT='(''  CALL TO W3FI82 WITH = '',/,
C    &                  10(3X,10I6,/))') (IWORK(I),I=1,NPTS)
C
               CALL w3fi82 (iwork,fval1,fdiff1,npts,pds,igds)
C
C             PRINT *,'GLAHN',FVAL1,FDIFF1
C             WRITE (6,FMT='(''  OUT FROM W3FI82 WITH = '',/,
C    &                  10(3X,10I6,/))') (IWORK(I),I=1,NPTS)
C
C             MUST NOW RE-REMOVE THE MINIMUM VALUE
C             OF THE SECOND DIFFERENCES TO ASSURE
C             ALL POSITIVE NUMBERS FOR SECOND ORDER GRIB PACKING
C
C             ORIGINAL REFERENCE VALUE ADDED TO FIRST POINT
C             VALUE FROM THE 2ND DIFF PACKER TO BE ADDED
C             BACK IN WHEN THE 2ND DIFF VALUES ARE
C             RECONSTRUCTED BACK TO THE BASIC VALUES
C
C             ALSO, THE REFERENCE VALUE IS
C             POWER-OF-TWO SCALED TO MATCH
C             FVAL1.  ALL OF THIS SCALING
C             WILL BE REMOVED AFTER THE
C             GLAHN SECOND DIFFERENCING IS UNDONE.
C             THE SCALING FACTOR NEEDED TO DO THAT
C             IS SAVED IN THE PDS AS A SIGNED POSITIVE
C             TWO BYTE INTEGER
C
C             THE SCALING FOR THE 2ND DIF PACKED
C             VALUES IS PROPERLY SET TO ZERO
C
              fval1 = fval1 + refnce*scal2
C                                          FIRST TEST TO SEE IF
C                                          ON 32 OR 64 BIT COMPUTER
C              CALL W3FI01(LW)
              IF (bit_size(lw).EQ.32) THEN
                  CALL w3fi76 (fval1,iexp,imant,32)
              ELSE
                  CALL w3fi76 (fval1,iexp,imant,64)
              END IF
              CALL sbytec(pds,iexp,320,8)
              CALL sbytec(pds,imant,328,24)
C
              IF (bit_size(lw).EQ.32) THEN
                  CALL w3fi76 (fdiff1,iexp,imant,32)
              ELSE
                  CALL w3fi76 (fdiff1,iexp,imant,64)
              END IF
              CALL sbytec(pds,iexp,352,8)
              CALL sbytec(pds,imant,360,24)
C
C             TURN ISCAL2 INTO SIGNED POSITIVE INTEGER
C             AND STORE IN TWO BYTES
C
              IF(iscal2.GE.0)  THEN
                CALL sbytec(pds,iscal2,384,16)
              ELSE
                CALL sbytec(pds,1,384,1)
                iscal2 = - iscal2
                CALL sbytec( pds,iscal2,385,15)
              ENDIF
C
              max  = iwork(1)
              min  = iwork(1)
              DO 700 i = 2, npts
                  IF (iwork(i).LT.min) THEN
                      min  = iwork(i)
                  ELSE IF (iwork(i).GT.max) THEN
                      max  = iwork(i)
                  END IF
  700         CONTINUE
C                           EXTRACT MINIMA
              DO 710 i = 1, npts
                  iwork(i)  = iwork(i) - min
  710         CONTINUE
              refnce  = min
C             PRINT *,'710 REFERENCE',REFNCE
              iscal2 = 0
C
C             AND RESET VALUE OF KWIDE - THE BIT WIDTH
C             FOR THE RANGE OF THE VALUES
C
              idiff = max - min
              CALL fi7505 (idiff,kwide)
C
C             PRINT*,'BIT WIDTH (KWIDE) OF 2ND DIFFS', KWIDE
C
C  **************************** END OF GLAHN PACKING  ************
          ELSE IF (ibdsfl(2).EQ.1.AND.ibdsfl(7).EQ.0) THEN
C                        HAVE SECOND ORDER PACKING WITH NO SECOND ORDER
C                        BIT MAP. ERGO ROW BY ROW - COL BY COL
              CALL fi7503 (iwork,ipfld,npts,ibdsfl,bds11,
     *              len,lenbds,pds,refnce,iscal2,kwide,igds)
              RETURN
          END IF
C         WRITE (6,FMT='(''  CALL TO FI7501 WITH = '',/,
C    &                  10(3X,10I6,/))') (IWORK(I),I=1,NPTS)
C         WRITE (6,FMT='(''  END OF ARRAY = '',/,
C    &                  10(3X,10I6,/))') (IWORK(I),I=NPTS-5,NPTS)
C         PRINT*,' REFNCE,ISCAL2, KWIDE AT CALL TO FI7501',
C    &             REFNCE, ISCAL2,KWIDE
C
C                         SECOND ORDER PACKING
C
          CALL fi7501 (iwork,ipfld,npts,ibdsfl,bds11,
     *             len,lenbds,pds,refnce,iscal2,kwide)
C
C              BDS COMPLETELY ASSEMBLED IN FI7501 FOR SECOND ORDER
C              PACKING.
C
      ELSE
C                                      SIMPLE PACKING
C
C                PRINT*,'  SIMPLE FIRST ORDER PACKING...'
          IF (ibitl.EQ.0) THEN
C                PRINT*,' WITH VARIABLE BIT LENGTH'
C
C                  WITH VARIABLE BIT LENGTH, ADJUSTED
C                  TO ACCOMMODATE LARGEST VALUE
C                  BINARY SCALING ALWAYS = 0
C
              CALL w3fi58(iwork,npts,iwork,pfld,nbits,len,kmin)
              rmin   = kmin
              refnce  = rmin
              iscale = 0
C             PRINT*,'  BIT LENGTH CAME OUT AT ...',NBITS
C
C           SET CONST .TRUE. IF ALL VALUES ARE THE SAME
C
              IF (len.EQ.0.AND.nbits.EQ.0) const = .true.
C
          ELSE
C           PRINT*,' FIXED BIT LENGTH, IBITL = ', IBITL
C
C             FIXED BIT LENGTH PACKING (VARIABLE PRECISION)
C             VALUES SCALED BY POWER OF 2 (ISCALE) TO
C             FIT LARGEST VALUE INTO GIVEN BIT LENGTH (IBITL)
C
              CALL w3fi59(fwork,npts,ibitl,iwork,pfld,iscale,len,rmin)
              refnce = rmin
C             PRINT *,' SCALING NEEDED TO FIT IS ...', ISCALE
              nbits = ibitl
C
C           SET CONST .TRUE. IF ALL VALUES ARE THE SAME
C
              IF (len.EQ.0) THEN
                  const = .true.
                  nbits = 0
              END IF
          END IF
C
C$        COMPUTE LENGTH OF BDS IN OCTETS
C
          inum  = npts * nbits + 88
C         PRINT *,'NUMBER OF BITS BEFORE FILL ADDED',INUM
C
C                  NUMBER OF FILL BITS
          nfill  = 0
          nleft  = mod(inum,16)
          IF (nleft.NE.0) THEN
              inum  = inum + 16 - nleft
              nfill = 16 - nleft
          END IF
C         PRINT *,'NUMBER OF BITS AFTER FILL ADDED',INUM
C                  LENGTH OF BDS IN BYTES
          lenbds = inum / 8
C
C                2.0   FORM THE BINARY DATA SECTION (BDS).
C
C                 CONCANTENATE ALL FIELDS FOR BDS
C
C                               BYTES 1-3
          CALL sbytec (bds11,lenbds,0,24)
C
C                               BYTE  4
C                                       FLAGS
          CALL sbytec (bds11,ibdsfl(1),24,1)
          CALL sbytec (bds11,ibdsfl(2),25,1)
          CALL sbytec (bds11,ibdsfl(3),26,1)
          CALL sbytec (bds11,ibdsfl(4),27,1)
C                                        NR OF FILL BITS
          CALL sbytec (bds11,nfill,28,4)
C
C$      FILL OCTETS 5-6 WITH THE SCALE FACTOR.
C
C                               BYTE  5-6
          IF (iscale.LT.0) THEN
              CALL sbytec (bds11,1,32,1)
              iscale  = - iscale
              CALL sbytec (bds11,iscale,33,15)
          ELSE
              CALL sbytec (bds11,iscale,32,16)
          END IF
C
C$  FILL OCTET 7-10 WITH THE REFERENCE VALUE
C   CONVERT THE FLOATING POINT OF YOUR MACHINE TO IBM370 32 BIT
C   FLOATING POINT NUMBER
C
C                               BYTE  7-10
C                                        REFERENCE VALUE
C                                          FIRST TEST TO SEE IF
C                                          ON 32 OR 64 BIT COMPUTER
C          CALL W3FI01(LW)
          IF (bit_size(lw).EQ.32) THEN
              CALL w3fi76 (refnce,iexp,imant,32)
          ELSE
              CALL w3fi76 (refnce,iexp,imant,64)
          END IF
          CALL sbytec (bds11,iexp,48,8)
          CALL sbytec (bds11,imant,56,24)
C
C
C$                        FILL OCTET 11 WITH THE NUMBER OF BITS.
C
C                               BYTE  11
          CALL sbytec (bds11,nbits,80,8)
      END IF
C
      RETURN
      END
C
C> @brief BDS second order packing.
C> @author Bill Cavanaugh @date 1993-08-06
 
C> Perform secondary packing on grid point data, generating all BDS information.
C>
C> Program history log:
C> - Bill Cavanaugh 1993-08-06
C> - Bill Cavanaugh 1993-12-15 Corrected location of start of first order
C> values and start of second order values to reflect a byte location in the
C> BDS instead of an offset.
C> - Mark Iredell 1995-10-31 Removed saves and prints
C>
C> @param[in] IWORK Integer source array
C> @param[in] NPTS Number of points in iwork
C> @param[in] IBDSFL Flags
C> @param[out] IPFLD Contains bds from byte 12 on
C> @param[out] BDS11 Contains first 11 bytes for bds
C> @param[out] LEN Number of bytes from 12 on
C> @param[out] LENBDS Total length of bds
C> @param PDS
C> @param REFNCE
C> @param ISCAL2
C> @param KWIDE
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
C> @author Bill Cavanaugh @date 1993-08-06
      SUBROUTINE fi7501 (IWORK,IPFLD,NPTS,IBDSFL,BDS11,
     *           LEN,LENBDS,PDS,REFNCE,ISCAL2,KWIDE)
 
      CHARACTER*1     BDS11(*),PDS(*)
C
      REAL            REFNCE
C
      INTEGER         ISCAL2,KWIDE
      INTEGER         LENBDS
      CHARACTER(len=1)         IPFLD(*)
      INTEGER         LEN,KBDS(22)
      INTEGER         IWORK(*)
C                        OCTET NUMBER IN SECTION, FIRST ORDER PACKING
C     INTEGER         KBDS(12)
C                        FLAGS
      INTEGER         IBDSFL(*)
C                        EXTENDED FLAGS
C     INTEGER         KBDS(14)
C                        OCTET NUMBER FOR SECOND ORDER PACKING
C     INTEGER         KBDS(15)
C                        NUMBER OF FIRST ORDER VALUES
C     INTEGER         KBDS(17)
C                        NUMBER OF SECOND ORDER PACKED VALUES
C     INTEGER         KBDS(19)
C                        WIDTH OF SECOND ORDER PACKING
      character(len=1)         ISOWID(400000)
C                        SECONDARY BIT MAP
      character(len=1)         ISOBMP(65600)
C                        FIRST ORDER PACKED VALUES
      character(len=1)         IFOVAL(400000)
C                        SECOND ORDER PACKED VALUES
      character(len=1)         ISOVAL(800000)
C
C     INTEGER         KBDS(11)
C                        BIT WIDTH TABLE
      INTEGER         IBITS(31)
C
      DATA            ibits/1,3,7,15,31,63,127,255,511,1023,
     *                      2047,4095,8191,16383,32767,65535,131072,
     *                      262143,524287,1048575,2097151,4194303,
     *                      8388607,16777215,33554431,67108863,
     *                      134217727,268435455,536870911,
     *                      1073741823,2147483647/
C  ----------------------------------
C                       INITIALIZE ARRAYS
 
      DO i = 1, 400000
          ifoval(i)  = char(0)
          isowid(i)  = char(0)
      ENDDO
C
      DO 101 i = 1, 65600
          isobmp(i)  = char(0)
  101 CONTINUE
      DO 102 i = 1, 800000
          isoval(i)  = char(0)
  102 CONTINUE
C                      INITIALIZE POINTERS
C                            SECONDARY BIT WIDTH POINTER
      iwdptr  = 0
C                            SECONDARY BIT MAP POINTER
      ibmp2p  = 0
C                            FIRST ORDER VALUE POINTER
      ifoptr  = 0
C                            BYTE POINTER TO START OF 1ST ORDER VALUES
      kbds(12)  = 0
C                            BYTE POINTER TO START OF 2ND ORDER VALUES
      kbds(15)  = 0
C                            TO CONTAIN NUMBER OF FIRST ORDER VALUES
      kbds(17)  = 0
C                            TO CONTAIN NUMBER OF SECOND ORDER VALUES
      kbds(19)  = 0
C                            SECOND ORDER PACKED VALUE POINTER
      isoptr  = 0
C  =======================================================
C
C                         DATA IS IN IWORK
C
      kbds(11)  = kwide
C
C       DATA PACKING
C
      iter    = 0
      inext   = 1
      istart  = 1
C  -----------------------------------------------------------
      kount = 0
C     DO 1 I = 1, NPTS, 10
C         PRINT *,I,(IWORK(K),K=I, I+9)
C   1 CONTINUE
 2000 CONTINUE
      iter  = iter + 1
C     PRINT *,'NEXT ITERATION STARTS AT',ISTART
       IF (istart.GT.npts) THEN
           GO TO 4000
       ELSE IF (istart.EQ.npts) THEN
           kpts    = 1
           mxdiff  = 0
           GO TO 2200
       END IF
C
C                     LOOK FOR REPITITIONS OF A SINGLE VALUE
       CALL fi7502 (iwork,istart,npts,isame)
       IF (isame.GE.15) THEN
           kount = kount + 1
C          PRINT *,'FI7501 - FOUND IDENTICAL SET OF ',ISAME
           mxdiff  = 0
           kpts    = isame
       ELSE
C
C                     LOOK FOR SETS OF VALUES IN TREND SELECTED RANGE
           CALL fi7513 (iwork,istart,npts,nmax,nmin,inrnge)
C          PRINT *,'ISTART  ',ISTART,' INRNGE',INRNGE,NMAX,NMIN
           iend  = istart + inrnge - 1
C          DO 2199 NM = ISTART, IEND, 10
C              PRINT *,'  ',(IWORK(NM+JK),JK=0,9)
C2199      CONTINUE
           mxdiff  = nmax - nmin
           kpts    = inrnge
       END IF
 2200 CONTINUE
C     PRINT *,'                 RANGE ',MXDIFF,' MAX',NMAX,' MIN',NMIN
C                 INCREMENT NUMBER OF FIRST ORDER VALUES
      kbds(17)  = kbds(17) + 1
C                 ENTER FIRST ORDER VALUE
      IF (mxdiff.GT.0) THEN
          DO 2220 lk = 0, kpts-1
              iwork(istart+lk)  = iwork(istart+lk) - nmin
 2220     CONTINUE
          CALL sbytec (ifoval,nmin,ifoptr,kbds(11))
      ELSE
          CALL sbytec (ifoval,iwork(istart),ifoptr,kbds(11))
      END IF
      ifoptr  = ifoptr + kbds(11)
C                  PROCESS SECOND ORDER BIT WIDTH
      IF (mxdiff.GT.0) THEN
          DO 2330 kwide = 1, 31
              IF (mxdiff.LE.ibits(kwide)) THEN
                  GO TO 2331
              END IF
 2330     CONTINUE
 2331     CONTINUE
      ELSE
          kwide  = 0
      END IF
      CALL sbytec (isowid,kwide,iwdptr,8)
      iwdptr  = iwdptr + 8
C         PRINT *,KWIDE,' IFOVAL=',NMIN,IWORK(ISTART),KPTS
C               IF KWIDE NE 0, SAVE SECOND ORDER VALUE
      IF (kwide.GT.0) THEN
          CALL sbytesc (isoval,iwork(istart),isoptr,kwide,0,kpts)
          isoptr  = isoptr + kpts * kwide
          kbds(19)  = kbds(19) + kpts
C         PRINT *,'            SECOND ORDER VALUES'
C         PRINT *,(IWORK(ISTART+I),I=0,KPTS-1)
      END IF
C                 ADD TO SECOND ORDER BITMAP
      CALL sbytec (isobmp,1,ibmp2p,1)
      ibmp2p  = ibmp2p + kpts
      istart  = istart + kpts
      GO TO 2000
C  --------------------------------------------------------------
 4000 CONTINUE
C     PRINT *,'THERE WERE ',ITER,' SECOND ORDER GROUPS'
C     PRINT *,'THERE WERE ',KOUNT,' STRINGS OF CONSTANTS'
C                 CONCANTENATE ALL FIELDS FOR BDS
C
C                   REMAINDER GOES INTO IPFLD
      iptr  = 0
C                               BYTES 12-13
C                                          VALUE FOR N1
C                                          LEAVE SPACE FOR THIS
      iptr   = iptr + 16
C                               BYTE 14
C                                          EXTENDED FLAGS
      CALL sbytec (ipfld,ibdsfl(5),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(6),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(7),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(8),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(9),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(10),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(11),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(12),iptr,1)
      iptr  = iptr + 1
C                               BYTES 15-16
C                 SKIP OVER VALUE  FOR N2
      iptr  = iptr + 16
C                               BYTES 17-18
C                                     P1
      CALL sbytec (ipfld,kbds(17),iptr,16)
      iptr  = iptr + 16
C                               BYTES 19-20
C                                   P2
      CALL sbytec (ipfld,kbds(19),iptr,16)
      iptr  = iptr + 16
C                               BYTE 21 - RESERVED LOCATION
      CALL sbytec (ipfld,0,iptr,8)
      iptr  = iptr + 8
C                               BYTES 22 - ?
C                                      WIDTHS OF SECOND ORDER PACKING
      ix    = (iwdptr + 32) / 32
C      CALL SBYTESC (IPFLD,ISOWID,IPTR,32,0,IX)
      ijk=iwdptr/8
      jst=(iptr/8)+1
      ipfld(jst:jst+ijk)=isowid(1:ijk)
      iptr  = iptr + iwdptr
C                                      SECONDARY BIT MAP
      ij    = (ibmp2p + 32) / 32
C      CALL SBYTESC (IPFLD,ISOBMP,IPTR,32,0,IJ)
      ijk=(ibmp2p/8)+1
      jst=(iptr/8)+1
      ipfld(jst:jst+ijk)=isobmp(1:ijk)
      iptr  = iptr + ibmp2p
      IF (mod(iptr,8).NE.0) THEN
          iptr  = iptr + 8 - mod(iptr,8)
      END IF
C                                         DETERMINE LOCATION FOR START
C                                         OF FIRST ORDER PACKED VALUES
      kbds(12)  = iptr / 8 + 12
C                                        STORE LOCATION
      CALL sbytec (ipfld,kbds(12),0,16)
C                                     MOVE IN FIRST ORDER PACKED VALUES
      ipass   = (ifoptr + 32) / 32
C      CALL SBYTESC (IPFLD,IFOVAL,IPTR,32,0,IPASS)
      ijk=(ifoptr/8)+1
      jst=(iptr/8)+1
      ipfld(jst:jst+ijk)=ifoval(1:ijk)
      iptr  = iptr + ifoptr
      IF (mod(iptr,8).NE.0) THEN
          iptr  = iptr + 8 - mod(iptr,8)
      END IF
C     PRINT *,'IFOPTR =',IFOPTR,' ISOPTR =',ISOPTR
C                DETERMINE LOCATION FOR START
C                     OF SECOND ORDER VALUES
      kbds(15)  = iptr / 8 + 12
C                                   SAVE LOCATION OF SECOND ORDER VALUES
      CALL sbytec (ipfld,kbds(15),24,16)
C                  MOVE IN SECOND ORDER PACKED VALUES
      ix    = (isoptr + 32) / 32
c      CALL SBYTESC (IPFLD,ISOVAL,IPTR,32,0,IX)
      ijk=(isoptr/8)+1
      jst=(iptr/8)+1
      ipfld(jst:jst+ijk)=isoval(1:ijk)
      iptr  = iptr + isoptr
      nleft  = mod(iptr+88,16)
      IF (nleft.NE.0) THEN
          nleft  = 16 - nleft
          iptr   = iptr + nleft
      END IF
C                                COMPUTE LENGTH OF DATA PORTION
      len     = iptr / 8
C                                    COMPUTE LENGTH OF BDS
      lenbds  = len + 11
C  -----------------------------------
C                               BYTES 1-3
C                                   THIS FUNCTION COMPLETED BELOW
C                                   WHEN LENGTH OF BDS IS KNOWN
      CALL sbytec (bds11,lenbds,0,24)
C                               BYTE  4
      CALL sbytec (bds11,ibdsfl(1),24,1)
      CALL sbytec (bds11,ibdsfl(2),25,1)
      CALL sbytec (bds11,ibdsfl(3),26,1)
      CALL sbytec (bds11,ibdsfl(4),27,1)
C                              ENTER NUMBER OF FILL BITS
      CALL sbytec (bds11,nleft,28,4)
C                               BYTE  5-6
      IF (iscal2.LT.0) THEN
          CALL sbytec (bds11,1,32,1)
          iscal2 = - iscal2
      ELSE
          CALL sbytec (bds11,0,32,1)
      END IF
      CALL sbytec (bds11,iscal2,33,15)
C
C$  FILL OCTET 7-10 WITH THE REFERENCE VALUE
C   CONVERT THE FLOATING POINT OF YOUR MACHINE TO IBM370 32 BIT
C   FLOATING POINT NUMBER
C                                        REFERENCE VALUE
C                                          FIRST TEST TO SEE IF
C                                          ON 32 OR 64 BIT COMPUTER
C          CALL W3FI01(LW)
          IF (bit_size(lw).EQ.32) THEN
              CALL w3fi76 (refnce,iexp,imant,32)
          ELSE
              CALL w3fi76 (refnce,iexp,imant,64)
          END IF
          CALL sbytec (bds11,iexp,48,8)
          CALL sbytec (bds11,imant,56,24)
C
C                               BYTE  11
C
      CALL sbytec (bds11,kbds(11),80,8)
C
      RETURN
      END
C
C> @brief Second order same value collection.
C> @author Bill Cavanaugh @date 1993-06-23
 
C> Collect sequential same values for processing
C> as second order value for grib messages.
C>
C> Program history log:
C> - Bill Cavanaugh 1993-06-23
C> - Mark Iredell 1995-10-31 Removed saves and prints
C>
C> @param[in] IWORK Array containing source data
C> @param[in] ISTART Starting location for this test
C> @param[in] NPTS Number of points in iwork
C> @param[out] ISAME Number of sequential points having the same value
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
C> @author Bill Cavanaugh @date 1993-06-23
      SUBROUTINE fi7502 (IWORK,ISTART,NPTS,ISAME)
 
      INTEGER        IWORK(*)
      INTEGER        ISTART
      INTEGER        ISAME
      INTEGER        K
      INTEGER        NPTS
C  -------------------------------------------------------------
      isame  = 0
      DO 100 k = istart, npts
          IF (iwork(k).NE.iwork(istart)) THEN
              RETURN
          END IF
          isame  = isame + 1
  100 CONTINUE
      RETURN
      END
C
C> @brief Row by row, col by col packing.
C> @author Bill Cavanaugh @date 1993-08-06
 
C> Perform row by row or column by column packing
C> generating all bds information.
C>
C> Program history log:
C> - Bill Cavanaugh 1993-08-06
C> - Mark Iredell 1995-10-31 Removed saves and prints
C>
C> @param[in] IWORK Integer source array
C> @param[in] NPTS Number of points in iwork
C> @param[in] IBDSFL Flags
C> @param[out] IPFLD Contains bds from byte 12 on
C> @param[out] BDS11 Contains first 11 bytes for bds
C> @param[out] LEN Number of bytes from 12 on
C> @param[out] LENBDS Total length of bds
C> @param PDS
C> @param REFNCE
C> @param ISCAL2
C> @param KWIDE
C> @param IGDS
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
C> @author Bill Cavanaugh @date 1993-08-06
      SUBROUTINE fi7503 (IWORK,IPFLD,NPTS,IBDSFL,BDS11,
     *           LEN,LENBDS,PDS,REFNCE,ISCAL2,KWIDE,IGDS)
 
      CHARACTER*1     BDS11(*),PDS(*),IPFLD(*)
C
      REAL            REFNCE
C
      INTEGER         ISCAL2,KWIDE
      INTEGER         LENBDS
      INTEGER         IGDS(*)
      INTEGER         LEN,KBDS(22)
      INTEGER         IWORK(*)
C                        OCTET NUMBER IN SECTION, FIRST ORDER PACKING
C     INTEGER         KBDS(12)
C                        FLAGS
      INTEGER         IBDSFL(*)
C                        EXTENDED FLAGS
C     INTEGER         KBDS(14)
C                        OCTET NUMBER FOR SECOND ORDER PACKING
C     INTEGER         KBDS(15)
C                        NUMBER OF FIRST ORDER VALUES
C     INTEGER         KBDS(17)
C                        NUMBER OF SECOND ORDER PACKED VALUES
C     INTEGER         KBDS(19)
C                        WIDTH OF SECOND ORDER PACKING
      character(len=1)         ISOWID(400000)
C                        SECONDARY BIT MAP
      character(len=1)         ISOBMP(65600)
C                        FIRST ORDER PACKED VALUES
      character(len=1)         IFOVAL(400000)
C                        SECOND ORDER PACKED VALUES
      character(len=1)         ISOVAL(800000)
C
C     INTEGER         KBDS(11)
C  ----------------------------------
C                       INITIALIZE ARRAYS
C
      DO i = 1, 400000
          ifoval(i)  = char(0)
          isowid(i)  = char(0)
      ENDDO
C
      DO 101 i = 1, 65600
          isobmp(i)  = char(0)
  101 CONTINUE
      DO 102 i = 1, 800000
          isoval(i)  = char(0)
  102 CONTINUE
C                      INITIALIZE POINTERS
C                            SECONDARY BIT WIDTH POINTER
      iwdptr  = 0
C                            SECONDARY BIT MAP POINTER
      ibmp2p  = 0
C                            FIRST ORDER VALUE POINTER
      ifoptr  = 0
C                            BYTE POINTER TO START OF 1ST ORDER VALUES
      kbds(12)  = 0
C                            BYTE POINTER TO START OF 2ND ORDER VALUES
      kbds(15)  = 0
C                            TO CONTAIN NUMBER OF FIRST ORDER VALUES
      kbds(17)  = 0
C                            TO CONTAIN NUMBER OF SECOND ORDER VALUES
      kbds(19)  = 0
C                            SECOND ORDER PACKED VALUE POINTER
      isoptr  = 0
C  =======================================================
C                         BUILD SECOND ORDER BIT MAP IN EITHER
C                         ROW BY ROW OR COL BY COL FORMAT
      IF (iand(igds(13),32).NE.0) THEN
C                              COLUMN BY COLUMN
          kout  = igds(4)
          kin   = igds(5)
C         PRINT *,'COLUMN BY COLUMN',KOUT,KIN
      ELSE
C                              ROW BY ROW
          kout  = igds(5)
          kin   = igds(4)
C         PRINT *,'ROW BY ROW',KOUT,KIN
      END IF
      kbds(17)  = kout
      kbds(19)  = npts
C
C     DO 4100 J = 1, NPTS, 53
C         WRITE (6,4101) (IWORK(K),K=J,J+52)
 4101     FORMAT (1x,25i4)
C         PRINT *,' '
C4100 CONTINUE
C
C                             INITIALIZE BIT MAP POINTER
      ibmp2p = 0
C                             CONSTRUCT WORKING BIT MAP
      DO 2000 i = 1, kout
          DO 1000 j = 1, kin
              IF (j.EQ.1) THEN
                  CALL sbytec (isobmp,1,ibmp2p,1)
              ELSE
                  CALL sbytec (isobmp,0,ibmp2p,1)
              END IF
              ibmp2p  = ibmp2p + 1
 1000     CONTINUE
 2000 CONTINUE
      len  = ibmp2p / 32 + 1
C     CALL BINARY(ISOBMP,LEN)
C
C                       PROCESS OUTER LOOP OF ROW BY ROW OR COL BY COL
C
      kptr  = 1
      kbds(11)  = kwide
      DO 6000 i = 1, kout
C                       IN CURRENT ROW OR COL
C                              FIND FIRST ORDER VALUE
          jptr  = kptr
          lowest  = iwork(jptr)
          DO 4000 j = 1, kin
              IF (iwork(jptr).LT.lowest) THEN
                  lowest = iwork(jptr)
              END IF
              jptr  = jptr + 1
 4000     CONTINUE
C                            SAVE FIRST ORDER VALUE
          CALL sbytec (ifoval,lowest,ifoptr,kwide)
          ifoptr  = ifoptr + kwide
C         PRINT *,'FOVAL',I,LOWEST,KWIDE
C                            SUBTRACT FIRST ORDER VALUE FROM OTHER VALS
C                                         GETTING SECOND ORDER VALUES
          jptr  = kptr
          ibig  = iwork(jptr) - lowest
          DO 4200 j = 1, kin
              iwork(jptr)  = iwork(jptr) - lowest
              IF (iwork(jptr).GT.ibig) THEN
                  ibig  = iwork(jptr)
              END IF
              jptr  = jptr + 1
 4200     CONTINUE
C                            HOW MANY BITS TO CONTAIN LARGEST SECOND
C                                         ORDER VALUE IN SEGMENT
          CALL fi7505 (ibig,nwide)
C                            SAVE BIT WIDTH
          CALL sbytec (isowid,nwide,iwdptr,8)
          iwdptr  = iwdptr + 8
C         PRINT *,I,'SOVAL',IBIG,' IN',NWIDE,' BITS'
C         WRITE (6,4101) (IWORK(K),K=KPTR,KPTR+52)
C                            SAVE SECOND ORDER VALUES OF THIS SEGMENT
          DO 5000 j = 0, kin-1
              CALL sbytec (isoval,iwork(kptr+j),isoptr,nwide)
              isoptr  = isoptr + nwide
 5000     CONTINUE
          kptr    = kptr + kin
 6000 CONTINUE
C  =======================================================
C                 CONCANTENATE ALL FIELDS FOR BDS
C
C                   REMAINDER GOES INTO IPFLD
      iptr  = 0
C                               BYTES 12-13
C                                          VALUE FOR N1
C                                          LEAVE SPACE FOR THIS
      iptr   = iptr + 16
C                               BYTE 14
C                                          EXTENDED FLAGS
      CALL sbytec (ipfld,ibdsfl(5),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(6),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(7),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(8),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(9),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(10),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(11),iptr,1)
      iptr  = iptr + 1
      CALL sbytec (ipfld,ibdsfl(12),iptr,1)
      iptr  = iptr + 1
C                               BYTES 15-16
C                 SKIP OVER VALUE  FOR N2
      iptr  = iptr + 16
C                               BYTES 17-18
C                                     P1
      CALL sbytec (ipfld,kbds(17),iptr,16)
      iptr  = iptr + 16
C                               BYTES 19-20
C                                   P2
      CALL sbytec (ipfld,kbds(19),iptr,16)
      iptr  = iptr + 16
C                               BYTE 21 - RESERVED LOCATION
      CALL sbytec (ipfld,0,iptr,8)
      iptr  = iptr + 8
C                               BYTES 22 - ?
C                                      WIDTHS OF SECOND ORDER PACKING
      ix    = (iwdptr + 32) / 32
C      CALL SBYTESC (IPFLD,ISOWID,IPTR,32,0,IX)
      ijk=iwdptr/8
      jst=(iptr/8)+1
      ipfld(jst:jst+ijk)=isowid(1:ijk)
      iptr  = iptr + iwdptr
C     PRINT *,'ISOWID',IWDPTR,IX
C     CALL BINARY (ISOWID,IX)
C
C                     NO SECONDARY BIT MAP
 
C                                         DETERMINE LOCATION FOR START
C                                         OF FIRST ORDER PACKED VALUES
      kbds(12)  = iptr / 8 + 12
C                                        STORE LOCATION
      CALL sbytec (ipfld,kbds(12),0,16)
C                                     MOVE IN FIRST ORDER PACKED VALUES
      ipass   = (ifoptr + 32) / 32
c      CALL SBYTESC (IPFLD,IFOVAL,IPTR,32,0,IPASS)
      ijk=(ifoptr/8)+1
      jst=(iptr/8)+1
      ipfld(jst:jst+ijk)=ifoval(1:ijk)
      iptr  = iptr + ifoptr
C     PRINT *,'IFOVAL',IFOPTR,IPASS,KWIDE
C     CALL BINARY (IFOVAL,IPASS)
      IF (mod(iptr,8).NE.0) THEN
          iptr  = iptr + 8 - mod(iptr,8)
      END IF
C     PRINT *,'IFOPTR =',IFOPTR,' ISOPTR =',ISOPTR
C                DETERMINE LOCATION FOR START
C                     OF SECOND ORDER VALUES
      kbds(15)  = iptr / 8 + 12
C                                   SAVE LOCATION OF SECOND ORDER VALUES
      CALL sbytec (ipfld,kbds(15),24,16)
C                  MOVE IN SECOND ORDER PACKED VALUES
      ix    = (isoptr + 32) / 32
C      CALL SBYTESC (IPFLD,ISOVAL,IPTR,32,0,IX)
      ijk=(isoptr/8)+1
      jst=(iptr/8)+1
      ipfld(jst:jst+ijk)=isoval(1:ijk)
      iptr  = iptr + isoptr
C     PRINT *,'ISOVAL',ISOPTR,IX
C     CALL BINARY (ISOVAL,IX)
      nleft  = mod(iptr+88,16)
      IF (nleft.NE.0) THEN
          nleft  = 16 - nleft
          iptr   = iptr + nleft
      END IF
C                                COMPUTE LENGTH OF DATA PORTION
      len     = iptr / 8
C                                    COMPUTE LENGTH OF BDS
      lenbds  = len + 11
C  -----------------------------------
C                               BYTES 1-3
C                                   THIS FUNCTION COMPLETED BELOW
C                                   WHEN LENGTH OF BDS IS KNOWN
      CALL sbytec (bds11,lenbds,0,24)
C                               BYTE  4
      CALL sbytec (bds11,ibdsfl(1),24,1)
      CALL sbytec (bds11,ibdsfl(2),25,1)
      CALL sbytec (bds11,ibdsfl(3),26,1)
      CALL sbytec (bds11,ibdsfl(4),27,1)
C                              ENTER NUMBER OF FILL BITS
      CALL sbytec (bds11,nleft,28,4)
C                               BYTE  5-6
      IF (iscal2.LT.0) THEN
          CALL sbytec (bds11,1,32,1)
          iscal2 = - iscal2
      ELSE
          CALL sbytec (bds11,0,32,1)
      END IF
      CALL sbytec (bds11,iscal2,33,15)
C
C$  FILL OCTET 7-10 WITH THE REFERENCE VALUE
C   CONVERT THE FLOATING POINT OF YOUR MACHINE TO IBM370 32 BIT
C   FLOATING POINT NUMBER
C                                        REFERENCE VALUE
C                                          FIRST TEST TO SEE IF
C                                          ON 32 OR 64 BIT COMPUTER
C      CALL W3FI01(LW)
      IF (bit_size(lw).EQ.32) THEN
          CALL w3fi76 (refnce,iexp,imant,32)
      ELSE
          CALL w3fi76 (refnce,iexp,imant,64)
      END IF
      CALL sbytec (bds11,iexp,48,8)
      CALL sbytec (bds11,imant,56,24)
C
C                               BYTE  11
C
      CALL sbytec (bds11,kbds(11),80,8)
C
      klen  = lenbds / 4 + 1
C     PRINT *,'BDS11 LISTING',4,LENBDS
C     CALL BINARY (BDS11,4)
C     PRINT *,'IPFLD LISTING'
C     CALL BINARY (IPFLD,KLEN)
      RETURN
      END
C
C> @brief Determine number of bits to contain value.
C> @author Bill Cavanaugh @date 1993-06-23
 
C> Calculate number of bits to contain value n, with a maximum of 32 bits.
C>
C> Program history log:
C> - Bill Cavanaugh 1993-06-23
C> - Mark Iredell 1995-10-31 Removed saves and prints
C>
C> @param[in] N Integer value
C> @param[out] NBITS Number of bits to contain n
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
C> @author Bill Cavanaugh @date 1993-06-23
      SUBROUTINE fi7505 (N,NBITS)
 
      INTEGER        N,NBITS
      INTEGER        IBITS(31)
C
      DATA           ibits/1,3,7,15,31,63,127,255,511,1023,2047,
     *               4095,8191,16383,32767,65535,131071,262143,
     *               524287,1048575,2097151,4194303,8388607,
     *               16777215,33554431,67108863,134217727,268435455,
     *               536870911,1073741823,2147483647/
C  ----------------------------------------------------------------
C
      DO 1000 nbits = 1, 31
          IF (n.LE.ibits(nbits)) THEN
              RETURN
          END IF
 1000 CONTINUE
      RETURN
      END
C
C> @brief Select block of data for packing.
C> @author Bill Cavanaugh @date 1994-01-21
 
C> Select a block of data for packing
C>
C> Program history log:
C> - Bill Cavanaugh 1994-01-21
C> - Mark Iredell 1995-10-31 Removed saves and prints
C>
C> - Return address if encounter set of same values
C> @param[in] IWORK
C> @param[in] ISTART
C> @param[in] NPTS
C> @param[out] MAX
C> @param[out] MIN
C> @param[out] INRNGE
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
C> @author Bill Cavanaugh @date 1994-01-21
      SUBROUTINE fi7513 (IWORK,ISTART,NPTS,MAX,MIN,INRNGE)
 
      INTEGER        IWORK(*),NPTS,ISTART,INRNGE,INRNGA,INRNGB
      INTEGER        MAX,MIN,MXVAL,MAXB,MINB,MXVALB
      INTEGER        IBITS(31)
C
      DATA           ibits/1,3,7,15,31,63,127,255,511,1023,2047,
     *               4095,8191,16383,32767,65535,131071,262143,
     *               524287,1048575,2097151,4194303,8388607,
     *               16777215,33554431,67108863,134217727,268435455,
     *               536870911,1073741823,2147483647/
C  ----------------------------------------------------------------
C                        IDENTIFY NEXT BLOCK OF DATA FOR PACKING AND
C                           RETURN TO CALLER
C  ********************************************************************
      istrta  = istart
C
C                     GET BLOCK A
      CALL fi7516 (iwork,npts,inrnga,istrta,
     *                                  max,min,mxval,lwide)
C  ********************************************************************
C
      istrtb  = istrta + inrnga
 2000 CONTINUE
C                         IF HAVE PROCESSED ALL DATA, RETURN
      IF (istrtb.GT.npts) THEN
C                         NO MORE DATA TO LOOK AT
          inrnge  = inrnga
          RETURN
      END IF
C                     GET BLOCK B
      CALL fi7502 (iwork,istrtb,npts,isame)
      IF (isame.GE.15) THEN
C         PRINT *,'BLOCK B HAS ALL IDENTICAL VALUES'
C         PRINT *,'BLOCK A HAS INRNGE =',INRNGA
C                     BLOCK B CONTAINS ALL IDENTICAL VALUES
          inrnge  = inrnga
C                     EXIT WITH BLOCK A
          RETURN
      END IF
C                     GET BLOCK B
C
      istrtb  = istrta + inrnga
      CALL fi7516 (iwork,npts,inrngb,istrtb,
     *                                  maxb,minb,mxvalb,lwideb)
C     PRINT *,'BLOCK A',INRNGA,' BLOCK B',INRNGB
C  ********************************************************************
C                     PERFORM TREND ANALYSIS TO DETERMINE
C                     IF DATA COLLECTION CAN BE IMPROVED
C
      ktrnd  = lwide - lwideb
C     PRINT *,'TREND',LWIDE,LWIDEB
      IF (ktrnd.LE.0) THEN
C         PRINT *,'BLOCK A - SMALLER, SHOULD EXTEND INTO BLOCK B'
          mxval   = ibits(lwide)
C
C                     IF BLOCK A REQUIRES THE SAME OR FEWER BITS
C                             LOOK AHEAD
C                        AND GATHER THOSE DATA POINTS THAT CAN
C                        BE RETAINED IN BLOCK A
C                        BECAUSE THIS BLOCK OF DATA
C                            USES FEWER BITS
C
          CALL fi7518 (iret,iwork,npts,istrta,inrnga,inrngb,
     *                          max,min,lwide,mxval)
          IF(iret.EQ.1) GO TO 8000
C         PRINT *,'18 INRNGA IS NOW ',INRNGA
          IF (inrngb.LT.20) THEN
              RETURN
          ELSE
              GO TO 2000
          END IF
      ELSE
C         PRINT *,'BLOCK A - LARGER, B SHOULD EXTEND BACK INTO A'
          mxvalb  = ibits(lwideb)
C
C                     IF BLOCK B REQUIRES FEWER BITS
C                             LOOK BACK
C                            SHORTEN BLOCK A BECAUSE NEXT BLOCK OF DATA
C                            USES FEWER BITS
C
          CALL fi7517 (iret,iwork,npts,istrtb,inrnga,
     *                               maxb,minb,lwideb,mxvalb)
          IF(iret.EQ.1) GO TO 8000
C         PRINT *,'17 INRNGA IS NOW ',INRNGA
      END IF
C
C                           PACK UP BLOCK A
C                           UPDATA POINTERS
 8000 CONTINUE
      inrnge  = inrnga
C                           GET NEXT BLOCK A
 9000 CONTINUE
      RETURN
      END
C
C> @brief Scan number of points.
C> @author Bill Cavanaugh @date 1994-01-21
 
C> Scan forward from current position. collect points and
C> determine maximum and minimum values and the number
C> of points that are included. Forward search is terminated
C> by encountering a set of identical values, by reaching
C> the number of points selected or by reaching the end
C> of data.
C>
C> Program history log:
C> - Bill Cavavnaugh 1994-01-21
C> - Mark Iredell 1995-10-31 Removed saves and prints
C>
C> - Return address if encounter set of same values
C> @param[in] IWORK Data array
C> @param[in] NPTS Number of points in data array
C> @param[in] ISTART Starting location in data
C> @param[out] INRNG Number of points selected
C> @param[out] MAX Maximum value of points
C> @param[out] MIN Minimum value of points
C> @param[out] MXVAL Maximum value that can be contained in lwidth bits
C> @param[out] LWIDTH Number of bits to contain max diff
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
C> @author Bill Cavanaugh @date 1994-01-21
      SUBROUTINE fi7516 (IWORK,NPTS,INRNG,ISTART,MAX,MIN,MXVAL,LWIDTH)
 
      INTEGER        IWORK(*),NPTS,ISTART,INRNG,MAX,MIN,LWIDTH,MXVAL
      INTEGER        IBITS(31)
C
      DATA           ibits/1,3,7,15,31,63,127,255,511,1023,2047,
     *               4095,8191,16383,32767,65535,131071,262143,
     *               524287,1048575,2097151,4194303,8388607,
     *               16777215,33554431,67108863,134217727,268435455,
     *               536870911,1073741823,2147483647/
C  ----------------------------------------------------------------
C
      inrng  = 1
      jq        = istart + 19
      max       = iwork(istart)
      min       = iwork(istart)
      DO 1000 i = istart+1, jq
          CALL fi7502 (iwork,i,npts,isame)
          IF (isame.GE.15) THEN
              GO TO 5000
          END IF
          inrng  = inrng + 1
          IF (iwork(i).GT.max) THEN
              max  = iwork(i)
          ELSE IF (iwork(i).LT.min) THEN
              min  = iwork(i)
          END IF
 1000 CONTINUE
 5000 CONTINUE
      krng   = max - min
C
      DO 9000 lwidth = 1, 31
          IF (krng.LE.ibits(lwidth)) THEN
C             PRINT *,'RETURNED',INRNG,' VALUES'
              RETURN
          END IF
 9000 CONTINUE
      RETURN
      END
C
C> @brief Scan backward.
C> @author Bill Cavanaugh @date 1994-01-21
 
C> Scan backwards until a value exceeds range of group b this may shorten group a
C>
C> Program history log:
C> - Bill Cavanaugh 1994-01-21
C> - Mark Iredell 1995-10-31 Removed saves and prints
C> - Mark Iredell 1998-06-17 Removed alternate return
C>
C> @param[in] IWORK
C> @param[in] ISTRTB
C> @param[in] NPTS
C> @param[in] INRNGA
C> @param[out] IRET
C> @param[out] MAXB
C> @param[out] MINB
C> @param MXVALB
C> @param LWIDEB
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
C> @author Bill Cavanaugh @date 1994-01-21
      SUBROUTINE fi7517 (IRET,IWORK,NPTS,ISTRTB,INRNGA,
     *                           MAXB,MINB,MXVALB,LWIDEB)
 
      INTEGER        IWORK(*),NPTS,ISTRTB,INRNGA
      INTEGER        MAXB,MINB,LWIDEB,MXVALB
      INTEGER        IBITS(31)
C
      DATA           ibits/1,3,7,15,31,63,127,255,511,1023,2047,
     *               4095,8191,16383,32767,65535,131071,262143,
     *               524287,1048575,2097151,4194303,8388607,
     *               16777215,33554431,67108863,134217727,268435455,
     *               536870911,1073741823,2147483647/
C  ----------------------------------------------------------------
      iret=0
C     PRINT *,'          FI7517'
      npos  = istrtb - 1
      itst  = 0
      kset  = inrnga
C
 1000 CONTINUE
C     PRINT *,'TRY NPOS',NPOS,IWORK(NPOS),MAXB,MINB
      itst  = itst + 1
      IF (itst.LE.kset) THEN
          IF (iwork(npos).GT.maxb) THEN
              IF ((iwork(npos)-minb).GT.mxvalb) THEN
C                 PRINT *,'WENT OUT OF RANGE AT',NPOS
                  iret=1
                  RETURN
              ELSE
                  maxb    = iwork(npos)
              END IF
          ELSE IF (iwork(npos).LT.minb) THEN
              IF ((maxb-iwork(npos)).GT.mxvalb) THEN
C                 PRINT *,'WENT OUT OF RANGE AT',NPOS
                  iret=1
                  RETURN
              ELSE
                  minb    = iwork(npos)
              END IF
          END IF
          inrnga  = inrnga - 1
          npos  = npos - 1
          GO TO 1000
      END IF
C  ----------------------------------------------------------------
C
 9000 CONTINUE
      RETURN
      END
C
C> @brief Scan forward.
C> @author Bill Cavanaugh @date 1994-01-21
 
C> Scan forward from start of block b towards end of block b
C> if next point under test forces a larger maxvala then
C> terminate indicating last point tested for inclusion
C> into block a.
C>
C> Program history log:
C> - Bill Cavanaugh 1994-01-21
C> - Mark Iredell 1995-10-31 Removed saves and prints
C> - Mark Iredell 1998-06-17 Removed alternate return
C>
C> @param IWORK
C> @param ISTRTA
C> @param INRNGA
C> @param INRNGB
C> @param MAXA
C> @param MINA
C> @param LWIDEA
C> @param MXVALA
C> @param[in] NPTS
C> @param[out] IRET
C>
C> @note Subprogram can be called from a multiprocessing environment.
C>
C> @author Bill Cavanaugh @date 1994-01-21
      SUBROUTINE fi7518 (IRET,IWORK,NPTS,ISTRTA,INRNGA,INRNGB,
     *                          MAXA,MINA,LWIDEA,MXVALA)
 
      INTEGER        IWORK(*),NPTS,ISTRTA,INRNGA
      INTEGER        MAXA,MINA,LWIDEA,MXVALA
      INTEGER        IBITS(31)
C
      DATA           IBITS/1,3,7,15,31,63,127,255,511,1023,2047,
     *               4095,8191,16383,32767,65535,131071,262143,
     *               524287,1048575,2097151,4194303,8388607,
     *               16777215,33554431,67108863,134217727,268435455,
     *               536870911,1073741823,2147483647/
C  ----------------------------------------------------------------
      iret=0
C     PRINT *,'          FI7518'
      npos  = istrta + inrnga
      itst  = 0
C
 1000 CONTINUE
      itst  = itst + 1
      IF (itst.LE.inrngb) THEN
C         PRINT *,'TRY NPOS',NPOS,IWORK(NPOS),MAXA,MINA
          IF (iwork(npos).GT.maxa) THEN
              IF ((iwork(npos)-mina).GT.mxvala) THEN
C                 PRINT *,'FI7518A -',ITST,' RANGE EXCEEDS MAX'
                  iret=1
                  RETURN
              ELSE
                  maxa    = iwork(npos)
              END IF
          ELSE IF (iwork(npos).LT.mina) THEN
              IF ((maxa-iwork(npos)).GT.mxvala) THEN
C                 PRINT *,'FI7518B -',ITST,' RANGE EXCEEDS MAX'
                  iret=1
                  RETURN
              ELSE
                  mina    = iwork(npos)
              END IF
          END IF
          inrnga  = inrnga + 1
C         PRINT *,'               ',ITST,INRNGA
          npos  = npos +1
          GO TO 1000
      END IF
C  ----------------------------------------------------------------
 9000 CONTINUE
      RETURN
      END