Re: [gembud] GDDIAG - MASK FUNCTION BUG?

You're welcome, Evan!

I think the bug comes from the subroutine $GEMPAK/source/gemlib/dp/dppgrb.f. There needs to be a test for missing data at all grid points. Perhaps this has already been corrected @ NCEP after v5.11.4, but for now I have attached a new version which fixes the problem, with the previous version also attached.

--Kevin

______________________________________________________________________
Kevin Tyle, Systems Administrator               **********************
Dept. of Atmospheric & Environmental Sciences   ktyle@xxxxxxxxxxxxxxxx
University at Albany, ES-235                    518-442-4578 (voice)
1400 Washington Avenue                          518-442-5825 (fax)
Albany, NY 12222                                **********************
______________________________________________________________________

On 10/11/2010 04:17 PM, Evan Lowery wrote:

THANK YOU Kevin!

That did the trick. I was starting to realize that it wasn't a MASK or GTE grid function error, and it had to be occurring after the new grid was calculated and before/while it was being written to the GEMPAK grid file but couldn't exactly pin it down. So there must be an error when the data is being packed into the grid like you said. I'm including an email chain with Michael James which has some additional troubleshooting just in case anyone experiences a similar issue in the future.

Regards,
Evan



-----Original Message-----
From: elowery@xxxxxxxxxxxx
Sent: Monday, October 11, 2010 9:10am
To: "Michael James" <mjames@xxxxxxxxxxxxxxxx>
Subject: Re: [gembud] GDDIAG - MASK FUNCTION BUG?

Good morning Michael,

I've narrowed down the problem a little more this morning... There are no problems with the MASK and SGE grid functions. What seems to be the problem on all platforms (64bit / 32 bit CentOS / RHEL) is that when you have a grid of ALL RMISSD (in any projection), GEMPAK GDDIAG converts all of the values to 1E31. Another simple example is below.

In the same file I sent yesterday (test_mer.grd), I created another grid (EMISS). The output below shows that when all grid values = -9999.00, GEMPAK writes the values 1E31 (9999999848243207295109594873856.00). It seems as there must be some check before writing the new grid file which calculates this erroneous value. If there is a mixture of valid and RMISSD values within the newly calculated grid, this error does not occur... I've been trying and will continue trying to identify where this error is occuring in the source code. If I find the problem I'll let you know where it is.

Regards,
Evan

1) Create a grid within test_mer.grd wiht the value -9999.00
GEMPAK-GDDIAG>l
 GDFILE   = test_mer.grd
 GDOUTF   = test_mer.grd
 GFUNC    = -9999.00
 GDATTIM  = 921025/0000
 GLEVEL   = 0
 GVCORD   = none
 GRDNAM   = emiss
 GRDTYP   = S
 GPACK    =
 GRDHDR   =
 PROJ     =
 GRDAREA  =
 KXKY     =
 MAXGRD   =
 CPYFIL   =
 ANLYSS   =
 GEMPAK-GDDIAG>r

    TIME1             TIME2         LEVL1 LEVL2   VCORD PARM
921025/0000                             0          NONE EMISS
Enter a new grid parameter name, <cr> to accept or type EXIT:
 Parameters requested: GDFILE,GDOUTF,GFUNC,GDATTIM,GLEVEL,GVCORD,GRDNAM,
 GRDTYP,GPACK,GRDHDR,PROJ,GRDAREA,KXKY,MAXGRD,CPYFIL,ANLYSS.
 GEMPAK-GDDIAG>

2) Export grid (emiss) into a dat file using GDLIST
GEMPAK-GDLIST>l
 GDATTIM  = 921025/0000
 GLEVEL   = 0
 GVCORD   = none
 GFUNC    = emiss
 GDFILE   = test_mer.grd
 GAREA    = us
 PROJ     = mer
 SCALE    = 0
 OUTPUT   = f/emiss.dat
 GEMPAK-GDLIST>r
Creating process: gn for queue 860487684


 GDLIST PARAMETERS:

 Grid file: test_mer.grd

 GRID IDENTIFIER:
    TIME1             TIME2         LEVL1 LEVL2   VCORD PARM
921025/0000                             0          NONE EMISS

 GAREA:    us
 SCALE FACTOR : 10** 0
 OUTPUT:    FILE/


MINIMUM AND MAXIMUM VALUES9999999848243207295109594873856.009999999848243207295109594873856.00
Enter <cr> to accept parameters or type EXIT:
Parameters requested: GDATTIM,GLEVEL,GVCORD,GFUNC,GDFILE,GAREA,PROJ,SCALE,
 OUTPUT.
 GEMPAK-GDLIST>

-----Original Message-----
From: elowery@xxxxxxxxxxxx
Sent: Sunday, October 10, 2010 10:08am
To: "Michael James" <mjames@xxxxxxxxxxxxxxxx>
Subject: Re: [gembud] GDDIAG - MASK FUNCTION BUG?

Hi Michael,

Thanks for looking into this. I tried to perform a similar (simplified) calculation with a Mercator GEMPAK grid file and the same error occurred. I'm not sure if it's a projection problem or calculation within the MASK function when this scenario occurs. I'll dig into the GEMPAK code for the MASK grid function tomorrow morning to see if I can determine what's happening. The commands I used are below, maybe you'll notice that I'm doing something wrong within the commands?

Previous Scenario:
Projection: CED
When trying to mask a grid (TEMP) to identify temperatures >=80F, GEMPAK calculated an erroneous value because the TEMP grid had NO values >=80F.

New Scenario:
Projection: Mercator
When trying to mask a grid (ONE) to identify values >=2, GEMPAK calculated an erroneous value because the ONE grid had NO values >=2.

Regards,
Evan

1) Create new GEMPAK grid file using GDCFIL
 GEMPAK-GDCFIL>l
 GDOUTF   = test_mer.grd
 PROJ     = MER
 GRDAREA  = us
 KXKY     = 20;20
 MAXGRD   = 5000
 CPYFIL   =
 ANLYSS   = 1/1;1;1;1

2) Create a grid within test_mer.grd with the value 1.00.
GEMPAK-GDDIAG>l
 GDFILE   = test_mer.grd
 GDOUTF   = test_mer.grd
 GFUNC    = 1
 GDATTIM  = 921025/0000
 GLEVEL   = 0
 GVCORD   = none
 GRDNAM   = one
 GRDTYP   = S
 GPACK    =
 GRDHDR   =
 PROJ     =
 GRDAREA  =
 KXKY     =
 MAXGRD   =
 CPYFIL   =
 ANLYSS   =
 GEMPAK-GDDIAG>r

    TIME1             TIME2         LEVL1 LEVL2   VCORD PARM
921025/0000                             0          NONE ONE
Enter a new grid parameter name, <cr> to accept or type EXIT:
 Parameters requested: GDFILE,GDOUTF,GFUNC,GDATTIM,GLEVEL,GVCORD,GRDNAM,
 GRDTYP,GPACK,GRDHDR,PROJ,GRDAREA,KXKY,MAXGRD,CPYFIL,ANLYSS.

3) Create a grid with the MASK function on "GRDNAM=one" where values >=2 (there are no values >=2, so all values in the new grid should be -9999.00/RMISSD
GEMPAK-GDDIAG>l
 GDFILE   = test_mer.grd
 GDOUTF   = test_mer.grd
 GFUNC    = mask(one,sge(one,2))
 GDATTIM  = 921025/0000
 GLEVEL   = 0
 GVCORD   = none
 GRDNAM   = onemsk
 GRDTYP   = S
 GPACK    =
 GRDHDR   =
 PROJ     = mer
 GRDAREA  =
 KXKY     =
 MAXGRD   =
 CPYFIL   =
 ANLYSS   =
 GEMPAK-GDDIAG>proj=
 GEMPAK-GDDIAG>r

    TIME1             TIME2         LEVL1 LEVL2   VCORD PARM
921025/0000                             0          NONE ONEMSK
Enter a new grid parameter name, <cr> to accept or type EXIT:
 Parameters requested: GDFILE,GDOUTF,GFUNC,GDATTIM,GLEVEL,GVCORD,GRDNAM,
 GRDTYP,GPACK,GRDHDR,PROJ,GRDAREA,KXKY,MAXGRD,CPYFIL,ANLYSS.

4) Export grids (one, onemsk) into a dat file using GDLIST to verify the values
GEMPAK-GDLIST>l
 GDATTIM  = 921025/0000
 GLEVEL   = 0
 GVCORD   = none
 GFUNC    = one
 GDFILE   = test_mer.grd
 GAREA    = us
 PROJ     = mer
 SCALE    = 0
 OUTPUT   = f/one.dat
 GEMPAK-GDLIST>r
Creating process: gn for queue 734035972


 GDLIST PARAMETERS:

 Grid file: test_mer.grd

 GRID IDENTIFIER:
    TIME1             TIME2         LEVL1 LEVL2   VCORD PARM
921025/0000                             0          NONE ONE

 GAREA:    us
 SCALE FACTOR : 10** 0
 OUTPUT:    FILE/


 MINIMUM AND MAXIMUM VALUES     1.00     1.00
Enter <cr> to accept parameters or type EXIT:
Parameters requested: GDATTIM,GLEVEL,GVCORD,GFUNC,GDFILE,GAREA,PROJ,SCALE,
 OUTPUT.

 GEMPAK-GDLIST>l
 GDATTIM  = 921025/0000
 GLEVEL   = 0
 GVCORD   = none
 GFUNC    = onemsk
 GDFILE   = test_mer.grd
 GAREA    = us
 PROJ     = mer
 SCALE    = 0
 OUTPUT   = f/onemsk.dat
 GEMPAK-GDLIST>r
Creating process: gn for queue 734101506


 GDLIST PARAMETERS:

 Grid file: test_mer.grd

 GRID IDENTIFIER:
    TIME1             TIME2         LEVL1 LEVL2   VCORD PARM
921025/0000                             0          NONE ONEMSK

 GAREA:    us
 SCALE FACTOR : 10** 0
 OUTPUT:    FILE/


MINIMUM AND MAXIMUM VALUES9999999848243207295109594873856.009999999848243207295109594873856.00
Enter <cr> to accept parameters or type EXIT:
Parameters requested: GDATTIM,GLEVEL,GVCORD,GFUNC,GDFILE,GAREA,PROJ,SCALE,
 OUTPUT.


-----Original Message-----
From: "Michael James" <mjames@xxxxxxxxxxxxxxxx>
Sent: Friday, October 8, 2010 5:11pm
To: "Evan Lowery" <elowery@xxxxxxxxxxxx>
Subject: Re: [gembud] GDDIAG - MASK FUNCTION BUG?

Evan,

I can confirm this on my end using your files.

However, I performed the same set of instructions using a GFS temperature grid and found min/max values were set to -9999/RMISSD as expected. This leads me to think the problem may be the specific projection in your file (but not really sure about anything at this point).

Michael James
Unidata



Sr. Business Meteorologist

Weather Trends International, Inc.

Phone 610-807-3197

http://www.wxtrends.com <http://www.wxtrends.com/>

This communication is privileged and may contain confidential information. It's intended only for the use of the person or entity named above. If you are not the intended recipient, do not distribute or copy this communication. If you have received this communication in error, please notify the sender immediately and return the original to the email address above. © Copyright 2010 Weather Trends International, Inc.

*From:* gembud-bounces@xxxxxxxxxxxxxxxx [mailto:gembud-bounces@xxxxxxxxxxxxxxxx] *On Behalf Of *Kevin R. Tyle
*Sent:* Monday, October 11, 2010 12:03 PM
*To:* gembud@xxxxxxxxxxxxxxxx
*Subject:* Re: [gembud] GDDIAG - MASK FUNCTION BUG?

Hi Evan,

In GDDIAG, try setting the packing variable "GPACK" to "NONE". The default, if "GPACK" is set to <blank> is to use "GRIB/16" packing. I think you may have unconvered a bug in the packing process which is not treating a grid whose points are all set equal to "RMISSD" properly.

--Kevin

______________________________________________________________________
Kevin Tyle, Systems Administrator               **********************
Dept. of Atmospheric & Environmental Sciences ktyle@xxxxxxxxxxxxxxxx <mailto:ktyle@xxxxxxxxxxxxxxxx>
University at Albany, ES-235                    518-442-4578 (voice)
1400 Washington Avenue                          518-442-5825 (fax)
Albany, NY 12222                                **********************
______________________________________________________________________

On 10/08/2010 12:20 PM, Evan Lowery wrote:

Hello all,

I've found a possible bug (or user error) with GDDIAG and the MASK grid function, but wanted to check with everyone here before sending out an official support request.

Within a GEMPAK grid file (test.grd), I have a temperature field (TEMPA) which needs to be masked to only show values between a certain temperature range (>=80F). I run this process daily, and have never had a problem up until this point. When NO temperatures in TEMPA are >=80F, the MASK function generates an erroneously large number rather than -9999.00 (RMISSD).

Dataset:                                               test.grd

gdlist

GDATTIM=101007/0000F001

GLEVEL=0

GVCORD=NONE

GFUNC=TEMPA

Using GDLIST (tempa.dat) I see that TEMPA:

MINIMUM AND MAXIMUM VALUES     1.03    60.00

Goal: only keep temperatures >=80F

gddiag

GDATTIM=101007/0000F001

GLEVEL=0

GVCORD=NONE

GFUNC=MASK(TEMPA,SGE(TEMPA,80))

GRDNAM=TEMPB

GRDTYP=S

GPACK=

GRDHDR=

PROJ=

GRDAREA=

KXKY=

MAXGRD=

CPYFIL=

ANLYSS=

Using GDLIST (tempb.dat) I see that TEMPB:

MINIMUM AND MAXIMUM VALUES9999999848243207295109594873856.009999999848243207295109594873856.00

I would expect all TEMPB values to be -9999.00 (RMISSD) since no temperatures are greater than 80F, but instead it blows up and returns a very large value.

http://www.unidata.ucar.edu/cgi-bin/gempak/manual/apxB_index

MASK  Masking function        MASK (S1, S2) = RMISSD IF S2 = RMISSD

                                            = S1 otherwise

In this example TEMPA had no values >=80F, but my csh scripts are constantly mining through temperature grids, and "usually" there are values >=80F.

Has anyone ever experienced this type of result? If yes, do you know a work around to get the grid (TEMPB) with all values = -9999.00 (RMISSD) rather than erroneously large values?

Regards,

Evan Lowery

_______________________________________________
gembud mailing list
gembud@xxxxxxxxxxxxxxxx  <mailto:gembud@xxxxxxxxxxxxxxxx>
For list information or to unsubscribe, visit:http://www.unidata.ucar.edu/mailing_lists/
        SUBROUTINE DP_PGRB  ( grid, igx, igy, nbits, idata, lendat,
     +                        qmin, scale, iret )
C************************************************************************
C* DP_PGRB                                                              *
C*                                                                      *
C* This subroutine packs a grid into the GEMPAK GRIB format using the   *
C* number of bits specified.  The packing and unpacking equations are:  *
C*                                                                      *
C*      IDATA  =  NINT ( ( GRID - QMIN ) / SCALE )                      *
C*      GRID   =  QMIN  +  IDATA * SCALE                                *
C*                                                                      *
C* DP_PGRB  ( GRID, IGX, IGY, NBITS, IDATA, LENDAT, QMIN, SCALE,        *
C*            IRET )                                                    *
C*                                                                      *
C* Input parameters:                                                    *
C*      GRID (IGX,IGY)  REAL            Grid data                       *
C*      IGX             INTEGER         Number of points in x dir       *
C*      IGY             INTEGER         Number of points in y dir       *
C*      NBITS           INTEGER         Number of bits                  *
C*                                                                      *
C* Output parameters:                                                   *
C*      IDATA (LENDAT)  INTEGER         Packed data                     *
C*      LENDAT          INTEGER         Length of packed data array     *
C*      QMIN            REAL            Minimum value of grid           *
C*      SCALE           REAL            Scaling factor                  *
C*      IRET            INTEGER         Return code                     *
C*                                        0 = normal return             *
C*                                      -10 = NBITS invalid             *
C*                                      -11 = invalid data range        *
C**                                                                     *
C* Log:                                                                 *
C* M. desJardins/GSFC    3/89                                           *
C* K. Brill/GSC          2/90   Fix to find negative max on grid        *
C* K. Brill/NMC         03/92   Fix for constant grid                   *
C* S. Chiswell/Unidata  10/02   Check for qmax-qmin > 2**32             *
C* H. Zeng/SAIC         09/07   Fixed constant grib with missing data   *
C* K. Tyle/UAlbany      10/10   Fixed for case with all missing data    *
C************************************************************************
        INCLUDE         'GEMPRM.PRM'
C*
        REAL            grid  (*)
        INTEGER         idata (*)
        LOGICAL         hvmis, allmis
C*
        INCLUDE         'ERMISS.FNC'
C------------------------------------------------------------------------
        iret = 0
        kxky = igx * igy
C
C*      Check for valid input.
C
        IF  ( ( nbits .le. 1 ) .or. ( nbits .gt. 31 ) )  THEN
            iret = -10
            RETURN
        END IF
C
C*      Compute the number of output words and initialize the output
C*      buffer and scaling parameter.
C
        lendat = FLOAT ( nbits * kxky ) / 32.0
        IF  ( lendat * 32 .ne. nbits * kxky ) lendat = lendat + 1
        DO  ii = 1, lendat
            idata (ii) = 0
        END DO
        scale = 1.
C
C*      Read through the grid finding the minimum and maximum values.
C
        hvmis = .false.
        allmis = .true.
        qmin = 1.0E+31
        qmax = -1.0E+31
        DO  ii = 1, kxky
            IF  ( .not. ERMISS ( grid (ii) ) )  THEN
                IF  ( grid (ii) .lt. qmin )  qmin = grid (ii)
                IF  ( grid (ii) .gt. qmax )  qmax = grid (ii)
                allmis = .false.
            ELSE
                hvmis = .true.
            END IF
        END DO
C
C*      If all points in the grid are set to missing, set the min and max 
values to missing.
C
        IF (allmis) THEN
            qmin = RMISSD
            qmax = RMISSD
        END IF
C
C*      Find the data range and check that it is valid.
C
        qdif = qmax - qmin
        IF  ( qdif .lt. 0.0 )  THEN
            iret = -11
            RETURN
        ELSE IF ( qdif .eq. 0.0 .and. .not. hvmis ) THEN
            RETURN
        END IF
C
C*      Find the scaling factor.  The scaling factor is set to a 
C*      power of 2.
C
        nnnn = 0
        idat = NINT (qdif)
        imax = 2 ** nbits - 1
        rtest = qdif * 2. ** nnnn

        IF  ( rtest .gt. ( 2. ** 31 - 1 ) ) THEN
            idat = 214748367
            write (*,*) 'Warning: range too big ', qdif, qmax, qmin,
     +                  idat, nbits
        ELSE
            idat = NINT ( qdif * 2. ** nnnn )
        END IF

        IF  ( qdif .ne. 0.0 )  THEN
C
            IF  ( idat .ge. imax )  THEN
                DO WHILE  ( idat .ge. imax )
                    nnnn = nnnn - 1
                    idat = qdif * 2.0 ** nnnn
                END DO
             ELSE
                DO WHILE  ( NINT ( qdif * 2.0 ** (nnnn+1) ) .lt. imax )
                    nnnn = nnnn + 1
                END DO
             END IF
C
        END IF
        scale = 2. ** nnnn
C
C*      Add data points to output buffer.
C
        iword = 1
        ibit  = 1
        DO  ii = 1, kxky
C
C*          Turn grid value into an integer.
C
            IF  ( ERMISS ( grid (ii) ) )  THEN
                idat = imax
              ELSE
                gggg = grid (ii) - qmin
                IF  ( gggg .lt. 0.0 ) gggg = 0.0
                idat = NINT ( gggg * scale )
            END IF
C
C*          Shift bits to correct location in word.
C
            jshft = 33 - nbits - ibit 
            idat2 = ISHFT ( idat, jshft )
            idata (iword) = IOR ( idata (iword), idat2 )
C
C*          Check to see if packed integer overflows into next word.
C
            IF  ( jshft .lt. 0 )  THEN
                jshft = 32 + jshft 
                idata (iword+1) = ISHFT ( idat, jshft )
            END IF
C
C*          Set location for next word.
C
            ibit = ibit + nbits
            IF  ( ibit .gt. 32 )  THEN
                ibit  = ibit - 32
                iword = iword + 1
            END IF
        END DO
C
C*      The scale factor to be saved is really the reciprocal of the
C*      scale which was computed.
C
        scale = 1.0 / scale
C*
        RETURN
        END
        SUBROUTINE DP_PGRB  ( grid, igx, igy, nbits, idata, lendat,
     +                        qmin, scale, iret )
C************************************************************************
C* DP_PGRB                                                              *
C*                                                                      *
C* This subroutine packs a grid into the GEMPAK GRIB format using the   *
C* number of bits specified.  The packing and unpacking equations are:  *
C*                                                                      *
C*      IDATA  =  NINT ( ( GRID - QMIN ) / SCALE )                      *
C*      GRID   =  QMIN  +  IDATA * SCALE                                *
C*                                                                      *
C* DP_PGRB  ( GRID, IGX, IGY, NBITS, IDATA, LENDAT, QMIN, SCALE,        *
C*            IRET )                                                    *
C*                                                                      *
C* Input parameters:                                                    *
C*      GRID (IGX,IGY)  REAL            Grid data                       *
C*      IGX             INTEGER         Number of points in x dir       *
C*      IGY             INTEGER         Number of points in y dir       *
C*      NBITS           INTEGER         Number of bits                  *
C*                                                                      *
C* Output parameters:                                                   *
C*      IDATA (LENDAT)  INTEGER         Packed data                     *
C*      LENDAT          INTEGER         Length of packed data array     *
C*      QMIN            REAL            Minimum value of grid           *
C*      SCALE           REAL            Scaling factor                  *
C*      IRET            INTEGER         Return code                     *
C*                                        0 = normal return             *
C*                                      -10 = NBITS invalid             *
C*                                      -11 = invalid data range        *
C**                                                                     *
C* Log:                                                                 *
C* M. desJardins/GSFC    3/89                                           *
C* K. Brill/GSC          2/90   Fix to find negative max on grid        *
C* K. Brill/NMC         03/92   Fix for constant grid                   *
C* S. Chiswell/Unidata  10/02   Check for qmax-qmin > 2**32             *
C* H. Zeng/SAIC         09/07   Fixed constant grib with missing data   *
C************************************************************************
        INCLUDE         'GEMPRM.PRM'
C*
        REAL            grid  (*)
        INTEGER         idata (*)
        LOGICAL         hvmis
C*
        INCLUDE         'ERMISS.FNC'
C------------------------------------------------------------------------
        iret = 0
        kxky = igx * igy
C
C*      Check for valid input.
C
        IF  ( ( nbits .le. 1 ) .or. ( nbits .gt. 31 ) )  THEN
            iret = -10
            RETURN
        END IF
C
C*      Compute the number of output words and initialize the output
C*      buffer and scaling parameter.
C
        lendat = FLOAT ( nbits * kxky ) / 32.0
        IF  ( lendat * 32 .ne. nbits * kxky ) lendat = lendat + 1
        DO  ii = 1, lendat
            idata (ii) = 0
        END DO
        scale = 1.
C
C*      Read through the grid finding the minimum and maximum values.
C
        hvmis = .false.
        qmin = 1.0E+31
        qmax = -1.0E+31
        DO  ii = 1, kxky
            IF  ( .not. ERMISS ( grid (ii) ) )  THEN
                IF  ( grid (ii) .lt. qmin )  qmin = grid (ii)
                IF  ( grid (ii) .gt. qmax )  qmax = grid (ii)
            ELSE
                hvmis = .true.
            END IF
        END DO
C
C*      Find the data range and check that it is valid.
C
        qdif = qmax - qmin
        IF  ( qdif .lt. 0.0 )  THEN
            iret = -11
            RETURN
        ELSE IF ( qdif .eq. 0.0 .and. .not. hvmis ) THEN
            RETURN
        END IF
C
C*      Find the scaling factor.  The scaling factor is set to a 
C*      power of 2.
C
        nnnn = 0
        idat = NINT (qdif)
        imax = 2 ** nbits - 1
        rtest = qdif * 2. ** nnnn

        IF  ( rtest .gt. ( 2. ** 31 - 1 ) ) THEN
            idat = 214748367
            write (*,*) 'Warning: range too big ', qdif, qmax, qmin,
     +                  idat, nbits
        ELSE
            idat = NINT ( qdif * 2. ** nnnn )
        END IF

        IF  ( qdif .ne. 0.0 )  THEN
C
            IF  ( idat .ge. imax )  THEN
                DO WHILE  ( idat .ge. imax )
                    nnnn = nnnn - 1
                    idat = qdif * 2.0 ** nnnn
                END DO
             ELSE
                DO WHILE  ( NINT ( qdif * 2.0 ** (nnnn+1) ) .lt. imax )
                    nnnn = nnnn + 1
                END DO
             END IF
C
        END IF
        scale = 2. ** nnnn
C
C*      Add data points to output buffer.
C
        iword = 1
        ibit  = 1
        DO  ii = 1, kxky
C
C*          Turn grid value into an integer.
C
            IF  ( ERMISS ( grid (ii) ) )  THEN
                idat = imax
              ELSE
                gggg = grid (ii) - qmin
                IF  ( gggg .lt. 0.0 ) gggg = 0.0
                idat = NINT ( gggg * scale )
            END IF
C
C*          Shift bits to correct location in word.
C
            jshft = 33 - nbits - ibit 
            idat2 = ISHFT ( idat, jshft )
            idata (iword) = IOR ( idata (iword), idat2 )
C
C*          Check to see if packed integer overflows into next word.
C
            IF  ( jshft .lt. 0 )  THEN
                jshft = 32 + jshft 
                idata (iword+1) = ISHFT ( idat, jshft )
            END IF
C
C*          Set location for next word.
C
            ibit = ibit + nbits
            IF  ( ibit .gt. 32 )  THEN
                ibit  = ibit - 32
                iword = iword + 1
            END IF
        END DO
C
C*      The scale factor to be saved is really the reciprocal of the
C*      scale which was computed.
C
        scale = 1.0 / scale
C*
        RETURN
        END