196 lines
7.1 KiB
Fortran
196 lines
7.1 KiB
Fortran
C$PRAGMA C (GBF_READ)
|
|
c =====================================================================
|
|
c pgm: ungribg (iupr,ibug,mbuf,kbuf,mfld,fld,
|
|
c mifld,ifld,kptr,kpds,kgds,kbms,
|
|
c lwcol,lsrow,lncol,lnrow,ngrib,gridf,istat)
|
|
c
|
|
c in: iupr .... unit number of output
|
|
c in: ibug .... debug control
|
|
c in: mbuf .... number of words in kbuf
|
|
c in: kbuf .... array to hold packed message
|
|
c in: mfld .... number of words in fld
|
|
c out: fld .... real array of unpacked values
|
|
c in: mifld .... number of words in ifld
|
|
c out: ifld .... integer array of unpacked values
|
|
c out: kptr .... decoded GRIB pointer array
|
|
c out: kpds .... decoded Product Definition Section parameters
|
|
c out: kgds .... decoded Grid Definition Section parameters
|
|
c out: kbms .... decoded bitmap array
|
|
c out: lwcol .... decoded western most column
|
|
c out: lsrow .... decoded southern most row
|
|
c out: lncol .... decoded number of HRAP columns
|
|
c out: lnrow .... decoded number of HRAP rows
|
|
c in: gridf .... grid factor 1=HRAP 4= 1/4 HRAP
|
|
c out: istat .... completion code
|
|
c =====================================================================
|
|
subroutine ungribg(iupr,ibug,mbuf,kbuf,mfld,fld,
|
|
+ mifld,ifld, kptr,kpds,kgds,kbms,
|
|
+ lwcol,lsrow,lncol,lnrow,ngrib,gridf,istat)
|
|
c.......................................................................
|
|
c This routine unpacks a GRIB file.
|
|
c
|
|
c.......................................................................
|
|
c ungrib Initially written by
|
|
c Tim Sweeney, HRL Jun 1999
|
|
c
|
|
c Added units labels
|
|
c Tim Sweeney, HRL Oct 1999
|
|
c
|
|
c Changed to read length of file from within GRIB, then read
|
|
c entire GRIB array
|
|
c Tim Sweeney, HRL Apr 2000
|
|
c
|
|
c Moved parameter print statements into new routine prgrib
|
|
c Tim Sweeney, HL Jan 2001
|
|
c
|
|
c Adapted from ungrib to use a grid factor to account for
|
|
c 1/4 HRAP grid used by EMPE for use by gribit routine
|
|
c
|
|
c David T. Miller, RSIS, OHD/HSEB Nov 2007
|
|
c.......................................................................
|
|
c
|
|
logical*1 kbms(*)
|
|
character*1 kbuf(mbuf)
|
|
c
|
|
dimension fld(mfld)
|
|
dimension ifld(mifld)
|
|
dimension kptr(*)
|
|
dimension kpds(*)
|
|
dimension kgds(*)
|
|
real*8 gridf
|
|
|
|
C
|
|
C ================================= RCS keyword statements ==========
|
|
CHARACTER*68 RCSKW1,RCSKW2
|
|
DATA RCSKW1,RCSKW2 / '
|
|
.$Source: /fs/hseb/ob82/ohd/pproc/src/gribit/RCS/ungrib.f,v $
|
|
. $', '
|
|
.$Id: ungrib.f,v 1.1 2006/10/19 16:06:04 dsa Exp $
|
|
. $' /
|
|
C ===================================================================
|
|
C
|
|
c
|
|
if (ibug.gt.0) write (iupr,*) 'enter ungribg - istat=',istat
|
|
c
|
|
if (istat.lt.0) istatt = 0
|
|
c
|
|
c progess GRIB
|
|
10 istat = 0
|
|
ipos = 1
|
|
nb = 1
|
|
c
|
|
c read to find 'GRIB' and to get length of file
|
|
20 call gbf_read (nb,kbuf(ipos),istat)
|
|
if (ibug.gt.0) write (iupr,30) istat
|
|
30 format (' in ungribg - after gbf_read: istat=',i4)
|
|
if (istat.ne.0) then
|
|
if (ipos.eq.1.and.ngrib.eq.0) then
|
|
write (iupr,40)
|
|
40 format (' ERROR: cannot read first byte of GRIB file. ',
|
|
+ 'File may not be a GRIB file.')
|
|
istat = -20
|
|
else
|
|
if (ibug.gt.0) write (iupr,50) nb
|
|
50 format (' NOTE: in ungribg - file not read for ',i4,' bytes. ',
|
|
+ 'Probably end-of-file.')
|
|
istat = 20
|
|
endif
|
|
go to 160
|
|
endif
|
|
if (nb.eq.1) istatt = istatt + 1
|
|
if (istatt.gt.100.and.ibug.gt.0) then
|
|
write (iupr,60)
|
|
60 format (' NOTE: in ungribg - another GRIB indicator not found ',
|
|
+ 'Possibly end-of-file')
|
|
istat = 20
|
|
go to 160
|
|
endif
|
|
if (kbuf(1).eq.'G'.and.nb.eq.1) then
|
|
ipos = 2
|
|
nb = 3
|
|
else if (nb.eq.3) then
|
|
istatt = 0
|
|
if (kbuf(2).ne.'R') go to 10
|
|
if (kbuf(3).ne.'I') go to 10
|
|
if (kbuf(4).ne.'B') go to 10
|
|
nb = 4
|
|
ipos = 5
|
|
else if (nb.eq.4) then
|
|
len=ichar(kbuf(5))*256*256+ichar(kbuf(6))*256+ichar(kbuf(7))
|
|
if (ibug.gt.2) then
|
|
ned = ichar(kbuf(8))
|
|
write (iupr,70) (kbuf(i),i=1,8),len,ned
|
|
70 format (' in ungrib - kbuf(1...8)=',8a1,' len=',i7,' ned=',i4)
|
|
write (iupr,80) ichar(kbuf(5)),ichar(kbuf(6)),
|
|
+ ichar(kbuf(7))
|
|
80 format(' in ungrib - kbuf(5)=',i4,' kbuf(6)=',i4,' kbuf(7)=',i4)
|
|
endif
|
|
nb = len - 8
|
|
ipos = 9
|
|
else if (nb.gt.8) then
|
|
if (ibug.gt.2) then
|
|
write (iupr,90) (kbuf(i),i=1,60)
|
|
90 format (' in ungrib - kbuf(1...60)=',60a1)
|
|
k = len - 60
|
|
write (iupr,100) k,len,(kbuf(i),i=k,len)
|
|
100 format (' in ungrib - k=',i2,' len=',i6,
|
|
+ ' kbuf(k...len)=',(1x,60a1 /))
|
|
endif
|
|
go to 110
|
|
endif
|
|
go to 20
|
|
c
|
|
c initialize arrays
|
|
110 do 120 i=1,25
|
|
kpds(i) = 0
|
|
120 continue
|
|
do 130 i=1,20
|
|
kptr(i) = 0
|
|
kgds(i) = 0
|
|
130 continue
|
|
c
|
|
c get GRIB field
|
|
call w3fi63 (kbuf,kpds,kgds,kbms,fld,kptr,istat)
|
|
if (istat.ne.0) then
|
|
write (iupr,140) istat
|
|
140 format (' ERROR: in GRIB decoding - w3fi63 istat=',i2 /
|
|
+ ' Values for istat:' /
|
|
+ 5x,' 1 = ''GRIB'' not found in first 100 characters' /
|
|
+ 5x,' 2 = ''7777'' not in correct location' /
|
|
+ 5x,' 3 = Unpacked field is larger than 65160' /
|
|
+ 5x,' 4 = GDS/ grid not one of currently accepted values' /
|
|
+ 5x,' 5 = Grid not currently available for center' /
|
|
+ 5x,' 7 = Edition indicated not included in decoder' /
|
|
+ 5x,' 8 = Temp GDS indicated, but GDS flag is off' /
|
|
+ 5x,' 9 = GDS indicates size mismatch with std grid' /
|
|
+ 5x,'10 = Incorrect center indicator' /
|
|
+ 5x,'11 = Binary data section (BDS) not completely processed' /
|
|
+ 5x,' Program is not set to process flag combinations' /
|
|
+ 5x,' shown in octets 4 and 14' /
|
|
+ 5x,'12 = Binary data section (BDS) not completely procesed.' /
|
|
+ 5x,' Program is not set to process flag combinations')
|
|
ccc go to 9999
|
|
endif
|
|
if (ibug.gt.0) write (iupr,150) (kbms(j),j=1,20)
|
|
150 format (' in ungrib - (kbms(j),j=1,20)=',20i2)
|
|
c
|
|
c get grid coordinates of southwest corner of field
|
|
rlat = kgds(4)/1000.
|
|
rlon = -kgds(5)/1000.
|
|
lncol = kgds(2)
|
|
lnrow = kgds(3)
|
|
if (kgds(1).eq.5) then
|
|
illgd = 1
|
|
npair = 1
|
|
call cvllgdg (rlon,rlat,npair,x,y,illgd,gridf,istat)
|
|
lwcol = int(x + 0.05)
|
|
lsrow = int(y + 0.05)
|
|
if (ibug.gt.0) write (iupr,*) 'in ungrib - lwcol=',lwcol,
|
|
+ ' lsrow=',lsrow,' lncol=',lncol,' lnrow=',lnrow
|
|
endif
|
|
c
|
|
ngrib=ngrib+1
|
|
c
|
|
160 return
|
|
c
|
|
end
|