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
