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.