NCEP Home > NCO Home > Systems Integration Branch > Decoders > BUFRLIB > BUFRLIB Table of Contents > Decoding NCEP PREPBUFR Files
Printer Friendly Version
Example Program: Decoding NCEP PREPBUFR files
In this example, we will demonstrate how to decode
(i.e. read) data from an NCEP PREPBUFR
file. These files are special BUFR files containing the entire set of data
(including quality control!)
that is input into the analysis step of a particular NCEP model run, and their
format is defined by the tables file
prepobs_prep.bufrtable.
The sample program below assumes that the PREPBUFR file to be read has already
been FORTRAN-blocked, via
the NCEP cwordsh utility or
otherwise, and that it is thus
directly readable as a result of using the FORM='UNFORMATTED' attribute within
the relevant FORTRAN
"OPEN" statement. However, aside from that, it will be
worthwhile to take the time now to discuss certain other aspects of the PREPBUFR
format in greater detail before actually delving into the sample program itself,
since such details will go a long way towards explaining and clarifying the
actual program design.
To begin with, a PREPBUFR file does not always contain,
within each single data subset, the data for an entire report!
Instead, for reports which contain
mass (i.e. temperature, moisture, etc.) as well as
wind (i.e. direction
and speed, U and V component, etc.) data values, such data values are stored
within two separate but adjacent (within the overall file) data subsets, where
each related subset, quite obviously,
contains the same report time, location, station
identification, etc. information as the other, but where the "mass" subset
contains the pressures and/or height levels at which "mass" data values
occur, while the corresponding "wind" subset contains the levels at which
"wind" data values occur. While it is true that this may, in some cases,
cause the same pressure and/or height level to appear in
both subsets, this separation is nonetheless
maintained for historical reasons peculiar to NCEP.
At any rate, the below program will actually merge all of the data
from both subsets into a single,
unified report in such cases, so that the final decoded output is clearer
and more intuitive.
In addition, and owing to the rich amount of quality control (QC) information
that
it contains, an NCEP PREPBUFR file contains an
events stack for each of the
basic data values such as pressure, temperature, wind, etc. in the file.
This "events stack" contains the history of QC activities
associated with that particular data value, including codes describing which
NCEP QC program modified that particular data value (and why?). The result is
a stack describing the entire QC history of that particular data value, from
the original reported value (at the bottom of the stack) to the "latest and
greatest" QC'ed value (at the top of the stack) that was used in the
respective NCEP model analysis. This "events stack" mechanism is implemented
within the BUFRLIB software via the use of subroutine
UFBEVN (and associated "[]"
notation surrounding the affected mnemonics within the
prepobs_prep.bufrtable file),
and it can perhaps be best pictured mentally
in the context of our usual 2-dimensional r8arr array,
containing data values for pressure, temperature, wind, etc. at some number of levels,
but where each
data value itself now has an "events stack" behind it, so that the
overall result is now a 3-dimensional "box" of data values, so to speak!
Of further note, since an NCEP PREPBUFR file contains the entire set of data
values that is input into a particular NCEP model analysis, it therefore
stands to reason that there would be many different types of data subsets
representing many different types of reports (e.g. rawinsonde, METAR, ship,
satellite, etc.) within such a file, and this is indeed reflected by the
large variety of Table A mnemonics that are listed in the tables file
prepobs_prep.bufrtable.
However, the method by which a PREPBUFR file is
constructed guarantees that all BUFR messages of a particular type (i.e.
Table A mnemonic) will occur adjacent to one another within the file, and this
feature, in tandem with the aforementioned realization that corresponding "mass"
and "wind" subsets for a particular report are always adjacent to one another
as well within the overall context of the file, is conveniently exploited
within the below sample program:
--------------------------- INCLUDE file "readpb.prm" ------------------------
REAL*8 R8BFMS
PARAMETER ( R8BFMS = 10.0E10 )
C* "Missing" value for BUFR data
C*
PARAMETER ( MXR8PM = 10 )
C* Maximum number of BUFR parameters
C*
PARAMETER ( MXR8LV = 255 )
C* Maximum number of BUFR levels
C*
PARAMETER ( MXR8VN = 10 )
C* Maximum number of BUFR event sequences
C*
PARAMETER ( MXR8VT = 6 )
C* Maximum number of BUFR variable types
C*
PARAMETER ( MXSTRL = 80 )
C* Maximum size of a string
C*
REAL*8 hdr ( MXR8PM ),
+ evns ( MXR8PM, MXR8LV, MXR8VN, MXR8VT )
C*
COMMON / PREPBC / hdr, evns, nlev
------------------------------------------------------------------------------
-------------------------- main program file "readpb.f" ----------------------
INCLUDE 'readpb.prm'
C*
CHARACTER outstg*(MXSTRL), subset*8
C*
CHARACTER var ( MXR8VT )
+ /'P','Q','T','Z','U','V'/
C*
PARAMETER ( NFILO = 15 )
INTEGER iunso ( NFILO )
+ / 51, 52, 53, 54, 55,
+ 56, 57, 58, 59, 60,
+ 61, 62, 63, 64, 65 /
CHARACTER*6 filo ( NFILO )
+ / 'ADPUPA', 'AIRCAR', 'AIRCFT', 'SATWND', 'PROFLR',
+ 'VADWND', 'SATBOG', 'SATEMP', 'ADPSFC', 'SFCSHP',
+ 'SFCBOG', 'SPSSMI', 'SYNDAT', 'ERS1DA', 'GOESND' /
C*
LOGICAL found
C-----------------------------------------------------------------------
C
C* Open the output files.
C
DO ii = 1, NFILO
OPEN ( UNIT = iunso ( ii ),
+ FILE = 'readpb.out.' // filo ( ii ) )
END DO
C
C* Open the input file.
C
OPEN ( UNIT = 11, FILE = 'prepbufr.in', FORM = 'UNFORMATTED' )
CALL OPENBF ( 11, 'IN', 11 )
CALL DATELEN ( 10 )
C
C* Get the next station report from the input file.
C
10 CALL READPB ( 11, subset, idate, ierrpb )
IF ( ierrpb .eq. -1 ) THEN
STOP
END IF
C
C* Set the appropriate output file unit number.
C
ii = 1
found = .false.
DO WHILE ( ( .not. found ) .and. ( ii .le. NFILO ) )
IF ( subset (1:6) .eq. filo ( ii ) ) THEN
found = .true.
iuno = iunso ( ii )
ELSE
ii = ii + 1
END IF
END DO
IF ( ( .not. found ) .and. ( ierrpb .eq. 0 ) ) THEN
GO TO 10
END IF
C
C* Print the HDR data for this station report.
C
WRITE ( UNIT = iuno,
+ FMT = '( /, A8, 1X, 2F7.2, 1X, F7.3, 1X, 2F8.1, ' //
+ '1X, F7.1, 1X, F6.1 )' )
+ ( hdr (ii), ii = 1, 8 )
C
C* Print the EVNS data for this station report.
C
DO lv = 1, nlev
WRITE ( UNIT = iuno, FMT = '( 80("-") )' )
WRITE ( UNIT = iuno, FMT = '( "lev ", I4, 7A9 )' )
+ lv, 'ob', 'qm', 'pc', 'rc', 'fc', 'an', 'cat'
WRITE ( UNIT = iuno, FMT = '( 80("-") )' )
DO kk = 1, MXR8VT
DO jj = 1, MXR8VN
WRITE ( UNIT = outstg, FMT = '( A8, 7(1X,F8.1) )' )
+ var (kk), ( evns ( ii, lv, jj, kk ), ii = 1, 7 )
DO mm = 1, MXSTRL
IF ( outstg (mm:mm) .eq. '*' ) THEN
outstg (mm:mm) = ' '
END IF
END DO
IF ( outstg (9:64) .ne. ' ' ) THEN
WRITE ( UNIT = iuno, FMT = '(A80)' ) outstg
ENDIF
END DO
END DO
END DO
C
IF ( ierrpb .eq. 0 ) THEN
GO TO 10
END IF
C*
STOP
END
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
SUBROUTINE READPB ( lunit, subset, idate, iret )
C
C* This subroutine will read and combine the mass and wind subsets
C* of the next station report in the prepbufr file. It is styled
C* after function IREADNS, and it only requires the prepbufr file
C* to be opened for reading with OPENBF. The combined station
C* report is returned to the caller in COMMON /PREPBC/.
C* This common area contains the number of levels in the report,
C* a one dimensional array with the header information, and a four
C* dimensional array containing all events from the variables POB,
C* QOB, TOB, ZOB, UOB, and VOB for the report.
C*
C* The header array contains the following list of mnemonics:
C*
C* SID XOB YOB DHR ELV TYP T29 ITP
C*
C* The 4-D array of data, EVNS ( ii, lv, jj, kk ), is indexed
C* as follows:
C*
C* "ii" indexes the event data types; these consist of:
C* 1) OBservation
C* 2) Quality Mark
C* 3) Program Code
C* 4) Reason Code
C* 5) ForeCast value
C* 6) ANalysed value
C* 7) office note CATegory
C* "lv" indexes the levels of the report
C* "jj" indexes the event stacks
C* "kk" indexes the variable types (p,q,t,z,u,v)
C*
C* Note that the structure of this array is identical to one
C* returned from UFBEVN, with an additional (4th) dimension to
C* include the six variable types into the same array.
C*
C* The return codes are as follows:
C* iret = 0 - normal return
C* = 1 - the station report within COMMON /PREPBC/ contains the
C* last available subset from within the prepbufr file
C* = -1 - there are no more subsets available from within the
C* prepbufr file
C*
INCLUDE 'readpb.prm'
C*
CHARACTER*(*) subset
C*
CHARACTER*(MXSTRL) head
+ / 'SID XOB YOB DHR ELV TYP T29 ITP' /
C*
CHARACTER*(MXSTRL) ostr ( MXR8VT )
+ / 'POB PQM PPC PRC PFC PAN CAT',
+ 'QOB QQM QPC QRC QFC QAN CAT',
+ 'TOB TQM TPC TRC TFC TAN CAT',
+ 'ZOB ZQM ZPC ZRC ZFC ZAN CAT',
+ 'UOB WQM WPC WRC UFC UAN CAT',
+ 'VOB WQM WPC WRC VFC VAN CAT' /
C*
REAL*8 hdr2 ( MXR8PM ),
+ evns2 ( MXR8PM, MXR8LV, MXR8VN, MXR8VT )
C*
REAL*8 r8sid, r8sid2, pob1, pob2
C*
CHARACTER*8 csid, csid2, subst2
C*
LOGICAL match / .true. /
C*
EQUIVALENCE ( r8sid, csid ), ( r8sid2, csid2 )
C*
SAVE match, subst2, idate2
C-----------------------------------------------------------------------
iret = 0
C*
C* If the previous call to this subroutine did not yield matching
C* mass and wind subsets, then IREADNS is already pointing at an
C* unmatched subset. Otherwise, call IREADNS to advance the subset
C* pointer to the next subset.
C*
IF ( match ) THEN
IF ( IREADNS ( lunit, subset, idate ) .ne. 0 ) THEN
iret = -1
RETURN
END IF
ELSE
subset = subst2
idate = idate2
END IF
C*
C* Read the HDR and EVNS data for the subset that is currently
C* being pointed to.
C*
CALL UFBINT ( lunit, hdr, MXR8PM, 1, jret, head )
DO ii = 1, MXR8VT
CALL UFBEVN ( lunit, evns ( 1, 1, 1, ii ), MXR8PM, MXR8LV,
+ MXR8VN, nlev, ostr (ii) )
END DO
C
C* Now, advance the subset pointer to the following subset and
C* read its HDR data.
C
IF ( IREADNS ( lunit, subset, idate ) .ne. 0 ) THEN
iret = 1
RETURN
END IF
CALL UFBINT ( lunit, hdr2, MXR8PM, 1, jret, head )
C
C* Check whether these two subsets have identical SID, YOB, XOB,
C* ELV, and DHR values. If so, then they are matching mass and
C* wind subsets for a single station report.
C
match = .true.
C
IF ( subset .ne. subst2 ) THEN
match = .false.
RETURN
END IF
C
r8sid = hdr (1)
r8sid2 = hdr2 (1)
IF ( csid .ne. csid2 ) THEN
match = .false.
RETURN
END IF
C
DO ii = 2, 5
IF ( hdr (ii) .ne. hdr2 (ii) ) THEN
match = .false.
RETURN
END IF
END DO
C
C* Read the EVNS data for the second of the two matching subsets.
C
DO ii = 1, MXR8VT
CALL UFBEVN ( lunit, evns2 ( 1, 1, 1, ii ), MXR8PM, MXR8LV,
+ MXR8VN, nlev2, ostr (ii) )
ENDDO
C
C* Combine the EVNS data for the two matching subsets into a
C* single 4-D array. Do this by merging the EVNS2 array into
C* the EVNS array.
C
DO 10 lv2 = 1, nlev2
DO lv = 1, nlev
pob1 = evns ( 1, lv, 1, 1 )
pob2 = evns2 ( 1, lv2, 1, 1 )
IF ( pob1 .eq. pob2 ) THEN
C
C* This pressure level from the second subset also exists
C* in the first subset, so overwrite any "missing" piece
C* of data for this pressure level in the first subset
C* with the corresponding piece of data from the second
C* subset (since this results in no net loss of data!).
C
DO kk = 1, MXR8VT
DO jj = 1, MXR8VN
DO ii = 1, MXR8PM
IF ( evns ( ii, lv, jj, kk ) .eq. R8BFMS ) THEN
evns ( ii, lv, jj, kk ) =
+ evns2 ( ii, lv2, jj, kk )
END IF
END DO
END DO
END DO
GO TO 10
ELSE IF ( ( pob2 .gt. pob1 ) .or.
+ ( lv .eq. nlev ) ) THEN
C
C* Either all remaining pressure levels within the first
C* subset are less than this pressure level from the
C* second subset (since levels within each subset are
C* guaranteed to be in descending order wrt pressure!)
C* *OR* there are more total levels within the second
C* subset than in the first subset. In either case, we
C* should now add this second subset level to the end of
C* the EVNS array.
C
nlev = nlev + 1
DO kk = 1, MXR8VT
DO jj = 1, MXR8VN
DO ii = 1, MXR8PM
evns ( ii, nlev, jj, kk ) =
+ evns2 ( ii, lv2, jj, kk )
END DO
END DO
END DO
GOTO 10
END IF
END DO
10 END DO
C*
RETURN
END
------------------------------------------------------------------------------
Unlike in some of our other examples, the above program can be directly compiled
and run "as is", without any additional code supplied by the user. The
program consists of a main routine and subroutine, where each of these
explicitly INCLUDEs the file "readpb.prm",
which itself contains an important
COMMON block as well as some other useful global definitions. The main routine
begins by opening a unique output file in which to store decoded data values
from each type of BUFR message that will be read. Then, the PREPBUFR file
itself, which presumably exists as filename "prepbufr.in" on the local system,
is assigned logical unit #11 and FORM='UNFORMATTED' since, again, we are
assuming that the file has already been FORTRAN-blocked for reading. Next, we
make the required call to OPENBF, where the astute reader
may have already
noticed that "11" was passed in as both the first and
third call argument,
meaning that the necessary BUFR tables information is embedded within the first
few BUFR messages of the file itself! Recall that this is the case for all
BUFR files (such as PREPBUFR) that are written out directly by the BUFRLIB
software; thus, when subsequently reading such files, we do not need to
go to the trouble of assigning an additional logical
unit to an ASCII-text version of the same BUFR tables file
(in this case, prepobs_prep.bufrtable)
for use by OPENBF.
Continuing on within the main program, a looping mechanism is now set up to
read each "complete" report, one at a time, from the PREPBUFR file.
By "complete", here we mean all of the basic data
(e.g. pressure,
temperature, moisture, wind, etc.) for a particular observing site
at a particular time and location, even when such data is stored within
the file across separate "mass" and "wind" data subsets! This is accomplished
via the included subroutine above, which uses function IREADNS to move
through the PREPBUFR file one data subset at a time, and which automatically
merges data from corresponding "mass" and "wind" subsets by taking advantage of
the fact that such corresponding subsets are guaranteed to be adjacent to each
other within the overall context of the file. The decoded data from each
so-called "complete" report is then written out to the appropriate output file
in an easy-to-read format, with each level of the report listed separately and
containing the entire "events stack" for every data value at that level.
Before leaving this discussion, we should point out that, for users who are
interested in working further with the NCEP PREPBUFR format, there is additional
documentation available, including explanations of the meanings of all of the
different QC program codes and reason codes used within PREPBUFR files. See the
NCEP/EMC PREPBUFR Processing Information Page
for all of the latest details.
|