C$$$  MAIN PROGRAM DOCUMENTATION BLOCK
C  
C MAIN PROGRAM:  READ_PREPBUFR
C   PRGMMR: KEYSER           ORG: NP22        DATE: 2002-01-28
C
C ABSTRACT: READS SUBSETS (REPORTS) FROM A PREPBUFR FILE, MERGING
C   MASS AND WIND SUBSETS FROM THE SAME ORIGINAL REPORT.  MERGED
C   REPORTS ARE PARSED ACCORDING TO THEIR BUFR MESSAGE TYPE AND
C   ARE LISTED IN OUTPUT TEXT FILES.
C
C PROGRAM HISTORY LOG:
C ????-??-??  WOOLLEN/GSC ORIGINAL AUTHOR
C ????-??-??  ATOR        STREAMLINED AND DOCUMENTED
C 2001-10-26  KEYSER      FURTHER DOCUMENTATION FOR WEB, UPDATED
C                         TO REPLACE OBSOLETE MESSAGE TYPE "SATBOG"
C                         WITH "QKSWND"; ADDED MORE VARAIBLES AND
C                         MANY OTHER UPDATES
C 2002-01-28  KEYSER      NOW THAT "PCAT" (PRECISION OF TEMPERATURE
C                         OBS.) HAS BEEN ADDED TO PREPBUFR FILE FOR
C                         "AIRCFT" AND "AIRCAR" MESSAGE TYPES, ADDED
C                         THIS TO LISTING
C
C USAGE:
C   INPUT FILES:
C     UNIT 11  - PREPBUFR FILE
C
C   OUTPUT FILES: 
C     UNIT 06  - UNIT 6 (STANDARD PRINTFILE)
C     UNIT 51  - LISTING OF REPORTS IN MESSAGE TYPE "ADPUPA"
C     UNIT 52  - LISTING OF REPORTS IN MESSAGE TYPE "AIRCAR"
C     UNIT 53  - LISTING OF REPORTS IN MESSAGE TYPE "AIRCFT"
C     UNIT 54  - LISTING OF REPORTS IN MESSAGE TYPE "SATWND"
C     UNIT 55  - LISTING OF REPORTS IN MESSAGE TYPE "PROFLR"
C     UNIT 56  - LISTING OF REPORTS IN MESSAGE TYPE "VADWND"
C     UNIT 57  - LISTING OF REPORTS IN MESSAGE TYPE "SATEMP"
C     UNIT 58  - LISTING OF REPORTS IN MESSAGE TYPE "ADPSFC"
C     UNIT 59  - LISTING OF REPORTS IN MESSAGE TYPE "SFCSHP"
C     UNIT 60  - LISTING OF REPORTS IN MESSAGE TYPE "SFCBOG"
C     UNIT 61  - LISTING OF REPORTS IN MESSAGE TYPE "SPSSMI"
C     UNIT 62  - LISTING OF REPORTS IN MESSAGE TYPE "SYNDAT"
C     UNIT 63  - LISTING OF REPORTS IN MESSAGE TYPE "ERS1DA"
C     UNIT 64  - LISTING OF REPORTS IN MESSAGE TYPE "GOESND"
C     UNIT 65  - LISTING OF REPORTS IN MESSAGE TYPE "QKSWND"
C
C   SUBPROGRAMS CALLED:
C     UNIQUE:    - READPB   INDEXF
C     LIBRARY:
C       W3LIB    - ERREXIT
C       BUFR     - OPENBF   DATELEN  READNS   UFBINT   UFBEVN
C                  UFBQCP
C
C   EXIT STATES:
C     COND =   0 - SUCCESSFUL RUN
C             22 - ABORT, NUMBER OF BALLOON DRIFT LEVELS RETURNED
C                  DOES NOT EQUAL NUMBER OF P/T/Z/Q/U/V LEVELS
C                  RETURNED (ADPUPA MESSAGE TYPE ONLY)
C             33 - ABORT, NUMBER OF VIRTUAL TEMP./DEWPOINT TEMP.
C                  LEVELS RETURNED DOES NOT EQUAL NUMBER OF
C                  P/T/Z/Q/U/V LEVELS RETURNED (ADPUPA MESSAGE
C                  TYPE ONLY)
C
C REMARKS: IN THIS PROGRAM, THE TERM "SUBSET" REFERS TO THE
C   COMPONENTS IN A BUFR MESSAGE.  THESE COULD BE CONSIDERED TO
C   BE "REPORTS" IN THE PREPBUFR FILE.  HERE, "REPORT" IS USED TO
C   REFER TO THE FINAL PRODUCT WHICH CONTAINS MERGED MASS AND WIND
C   SUBSETS FROM THE SAME ORIGINAL OBSERVATION.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE:  IBM SP, SGI
C
C$$$

      PROGRAM   READ_PREPBUFR

      REAL*8    R8BFMS

      PARAMETER ( R8BFMS = 9.0E08 )! Missing value threshhold for BUFR
      PARAMETER ( NHR8PM = 14 )    ! Actual number of BUFR parameters
                                   !  in header
      PARAMETER ( MXR8PM =  8 )    ! Maximum number of BUFR parameters
                                   !  in level data (non-radiance) or
                                   !  in channel data (radiance)
      PARAMETER ( MXR8LV = 255 )   ! Maximum number of BUFR levels/
                                   !                        channels
      PARAMETER ( MXR8VN = 10 )    ! Max. number of BUFR event sequences
                                   !  (non-radiance reports)
      PARAMETER ( MXR8VT = 17)     ! Max. number of BUFR variable types
                                   !  (non-radiance reports)
      PARAMETER ( NFILO = 15 )

      REAL*8    hdr ( NHR8PM ), qkswnd_hdr(3), aircar_hdr(6),
     +          aircft_hdr(4),  adpupa_hdr(1), goesnd1_hdr(6),
     +          goesnd2_hdr(MXR8LV), adpsfc_hdr(5), sfcshp_hdr(2),
     +          satwnd_hdr(1), evns ( MXR8PM, MXR8LV, MXR8VN, MXR8VT ),
     +          drft ( 3, MXR8LV ), brts ( MXR8PM, MXR8LV ),
     +          toth ( 2, MXR8LV )

      INTEGER   ifirst ( NFILO ), iunso ( NFILO ), indxvr1 ( 100:299 ),
     +          indxvr2 ( 100:299 )

      CHARACTER    msgtyp*8, qc_step*8
      CHARACTER*18 var ( MXR8VT )
      CHARACTER*6  filo ( NFILO )

      LOGICAL   found

      COMMON  / PREPBC /      hdr, evns, drft, toth, brts, qkswnd_hdr,
     +  aircar_hdr, aircft_hdr, adpupa_hdr, goesnd1_hdr, goesnd2_hdr,
     +  adpsfc_hdr, sfcshp_hdr, satwnd_hdr, indxvr1, indxvr2, nlev, nchn

      DATA      iunso
     +        /   51,   52,   53,   54,   55,
     +            56,   57,   58,   59,   60,
     +            61,   62,   63,   64,   65  /
      DATA         var /
     +  'PRESSURE (MB)     ','SP HUMIDITY(MG/KG)','TEMPERATURE (C)   ',
     +  'HEIGHT (METERS)   ','U-COMP WIND (M/S) ','V-COMP WIND (M/S) ',
     +  'WIND DIR (DEG)    ','WIND SPEED (KNOTS)','TOTAL PWATER (MM) ',
     +  'RAIN RATE (MM/HR) ','1.-.9 sig PWAT(MM)','.9-.7 sig PWAT(MM)',
     +  '.7-.3 sig PWAT(MM)','.3-0. sig PWAT(MM)','CLOUD-TOP PRES(MB)',
     +  'CLOUD-TOP TEMP (K)','TOT. CLD. COVER(%)' /
      DATA         filo
     +        / 'ADPUPA', 'AIRCAR', 'AIRCFT', 'SATWND', 'PROFLR',
     +          'VADWND', 'SATEMP', 'ADPSFC', 'SFCSHP', 'SFCBOG',
     +          'SPSSMI', 'SYNDAT', 'ERS1DA', 'GOESND', 'QKSWND'/

C-----------------------------------------------------------------------

      IFIRST = 0

C       Open the output files.
C       ----------------------

      DO ifile = 1, NFILO
          OPEN  ( UNIT = iunso ( ifile ),
     +            FILE = 'readpb.out.' // filo ( ifile ) )
      END DO

C       Open the input PREPBUFR file.
C       -----------------------------
      OPEN  ( UNIT = 11, FILE = 'prepbufr.in', FORM = 'UNFORMATTED' )
      CALL OPENBF(11,'IN',11 )  ! BUFRLIB routine to open file
      CALL DATELEN(10)  ! BUFRLIB routine to use 10-digit date

C       Get the next subset from the input file.  Merge mass and
C        wind subsets for the same "true" report.
C       --------------------------------------------------------

  10  CALL READPB(11,msgtyp,idate,ierrpb )
      IF ( ierrpb .eq. -1 )  THEN
          print *, '1-All subsets read in and processed'
          STOP      ! All subsets read in and processed
      END IF

C       Set the appropriate output file unit number based on the
C        BUFR message type the subset is in.
C       --------------------------------------------------------

      ifile = 1
      found = .false.
      DO WHILE  ( ( .not. found ) .and. ( ifile .le. NFILO ) )
          IF  ( msgtyp (1:6) .eq. filo ( ifile ) )  THEN
              found = .true.
              iuno = iunso ( ifile )
          ELSE
              ifile = ifile + 1
          END IF
      END DO
      IF  ( ( .not. found ) .and. ( ierrpb .eq. 0 ) )  THEN
          GO TO 10
      END IF

      IF ( ifirst(ifile) .eq. 0 )  THEN
         WRITE  ( UNIT = iuno,
     +    FMT = '("LISTING FOR MESSAGE TYPE ",A6//' //
     +    '"Key for header:"/' //
     +    '"SID  -- Station identification"/' //
     +    '"YOB  -- Latitude (N+, S-)"/' //
     +    '"XOB  -- Longitude (E)"/' //
     +    '"ELV  -- Elevation (meters)"/' //
     +    '"DHR  -- Observation time minus cycle time (hours)"/' //
     +    '"RPT  -- Reported observation time (hours)"/' //
     +    '"TCOR -- Indicator whether observation time used to ",' //
     +     '"generate DHR was corrected (0- Observation time ",' //
     +     '"is reported time, no"/9X,"correction, 1- Observation ",' //
     +     '"time corrected based on reported launch time, 2- ",' //
     +     '"Observation time corrected even though"/9X,"launch ", ' //
     +     '"time is missing)"/' //
     +    '"TYP  -- PREPBUFR report type"/' //
     +    '"TSB  -- PREPBUFR report subtype"/' //
     +    '"T29  -- Input report type"/' //
     +    '"ITP  -- Instrument type (BUFR code table 0-02-001)"/' //
     +    '"SQN  -- Report sequence number"/' //
     +    '"PROCN - Process number for this MPI run (obtained from ",'//
     +     '"script)"/' //
     +    '"SAID -- Satellite identifier (BUFR code table 0-01-007)")')
     +    filo(ifile)
         IF (filo(ifile) .eq. 'ADPUPA' ) then
            WRITE  ( UNIT = iuno, FMT =
     +       '("SIRC -- Rawinsonde solar & infrared radiation ",' //
     +        '"correction indicator (BUFR code table 0-02-013)")')
         ELSE IF (filo(ifile) .eq. 'GOESND' .or.
     +            filo(ifile) .eq. 'SATEMP' ) then
            WRITE  ( UNIT = iuno, FMT =
     +       '("For reports with non-radiance data:"/' //
     +       '"  ACAV -- Total number with respect to accumulation ",'//
     +        '"or average"/' //
     +       '"  PRSS -- Surface pressure observation (mb)"/' //
     +       '"For reports with radiance data:"/' //
     +       '"  ELEV -- Satellite elevation (zenith angle - deg.)"/' //
     +       '"  SOEL -- Solar elevation (zenith angle - deg.)"/' //
     +       '"  OZON -- Total ozone (Dobson units)"/' //
     +       '"  TMSK -- Skin temperature (Kelvin)"/' //
     +       '"  CLAM -- Cloud Amount (BUFR code table 0-20-011)")')
         ELSE IF (filo(ifile) .eq. 'QKSWND' ) then
            WRITE  ( UNIT = iuno, FMT =
     +       '("ATRN -- Along track row number"/' //
     +       '"CTCN -- Cross track cell number"/' //
     +       '"SPRR -- Seawinds probability of rain (%)")')
         ELSE IF (filo(ifile) .eq. 'AIRCAR' ) then
            WRITE  ( UNIT = iuno, FMT =
     +       '("PCAT -- Precision of temperature obs. (Kelvin)"/' //
     +       '"POAF -- Phase of aircraft flight (BUFR code table ",' //
     +        '"0-08-004)"/' //
     +       '"TRBX10 -- Turbulence index for period TOB-1 min to ",' //
     +        '"TOB"/' //
     +       '"TRBX21 -- Turbulence index for period TOB-2 min to ",' //
     +        '"TOB-1 min"/' //
     +       '"TRBX32 -- Turbulence index for period TOB-3 min to ",' //
     +        '"TOB-2 min"/' //
     +       '"TRBX43 -- Turbulence index for period TOB-4 min to ",' //
     +        '"TOB-3 min")')
         ELSE IF (filo(ifile) .eq. 'AIRCFT' ) then
            WRITE  ( UNIT = iuno, FMT =
     +       '("RCT  -- Receipt time (hours)"/' //
     +       '"PCAT -- Precision of temperature obs. (Kelvin)"/' //
     +       '"POAF -- Phase of aircraft flight (BUFR code table ",' //
     +        '"0-08-004)"/' //
     +       '"DGOT -- Degree of turbulence (BUFR code table ",' //
     +        '"0-11-031)")')
         ELSE IF (filo(ifile) .eq. 'ADPSFC' ) then
            WRITE  ( UNIT = iuno, FMT =
     +       '("PMO  -- Mean sea-level pressure observation (mb)"/' //
     +       '"PMQ  -- Mean sea-level pressure quality marker"/' //
     +       '"ALSE -- Altimeter setting observation (mb)"/' //
     +       '"SOB  -- Wind speed observation (stored when ",' //
     +        '"direction is missing) (m/s)"/' //
     +       '"SQM  -- Wind speed quality marker (stored when ",' //
     +        '"direction is missing)")')
         ELSE IF (filo(ifile) .eq. 'SFCSHP' ) then
            WRITE  ( UNIT = iuno, FMT =
     +       '("PMO  -- Mean sea-level pressure observation (mb)"/' //
     +       '"PMQ  -- Mean sea-level pressure quality marker")')
         ELSE IF (filo(ifile) .eq. 'SATWND' ) then
            WRITE  ( UNIT = iuno, FMT =
     +       '("RFFL -- NESDIS recursive filter flag")')
         END IF
         WRITE  ( UNIT = iuno, FMT = '(//"PREPBUFR FILE DATE IS: ",
     +    I10)') idate
         WRITE  ( UNIT = iuno, FMT = '(/)')
         ifirst(ifile)=1
      END IF

      WRITE  ( UNIT = iuno,
     + FMT = '(/"SID=",A8,", YOB=",F6.2,", XOB=",F6.2,", ELV=",I5,' //
     + '", DHR=",F6.3,", RPT=",F8.3,", TCOR=",I1,", TYP=",I3,' //
     + '", TSB=",I1,", T29=",I3,", ITP=",I2,", SQN=",I6/,"PROCN=",I2,'//
     + '", SAID=",I3)')
     + (hdr(ii),ii=1,3),nint(hdr(4)),hdr(5),hdr(6),(nint(hdr(ii)),
     + ii=7,14)

      IF (filo(ifile) .eq. 'ADPUPA' ) then
         WRITE  ( UNIT = iuno, FMT = '("SIRC=",I2)') nint(adpupa_hdr)
      ELSE IF (filo(ifile) .eq. 'GOESND' .or.
     +         filo(ifile) .eq. 'SATEMP' ) then
         IF (nlev.gt.0)  THEN
            WRITE  ( UNIT = iuno, FMT = '("ACAV=",I3,", PRSS=",F6.1)')
     +       nint(goesnd1_hdr(1)),goesnd2_hdr(1)*0.01
         END IF
         IF (nchn.gt.0)  THEN
            WRITE  ( UNIT = iuno, FMT = '("ELEV=",F6.2,", SOEL=",' //
     +       'F6.2,", OZON=",I4,", TMSK=",F5.1,", CLAM=",I2)')
     +       goesnd1_hdr(2),goesnd1_hdr(3),nint(goesnd1_hdr(4)),
     +       goesnd1_hdr(5),nint(goesnd1_hdr(6))
         END IF
      ELSE IF (filo(ifile) .eq. 'QKSWND' ) then
         WRITE  ( UNIT = iuno, FMT = '("ATRN=",I3,", CTCN=",I3,' //
     +   '", SPRR=",I3)') (nint(qkswnd_hdr(ii)),ii=1,3)
      ELSE IF (filo(ifile) .eq. 'AIRCAR' ) then
         WRITE  ( UNIT = iuno, FMT = '("PCAT=",F5.2,", POAF=",I3,' //
     + '", TRBX10=",I3,", TRBX21=",I3,", TRBX32=",I3,", TRBX43=",I3)')
     +    aircar_hdr(1),(nint(aircar_hdr(ii)),ii=2,6)
      ELSE IF (filo(ifile) .eq. 'AIRCFT' ) then
         WRITE  ( UNIT = iuno, FMT = '("RCT=",F7.2,", PCAT=",F5.2,' //
     +    '", POAF=",I3,", DGOT=",I3)') aircft_hdr(1),aircft_hdr(2),
     +    (nint(aircft_hdr(ii)),ii=3,4)
      ELSE IF (filo(ifile) .eq. 'ADPSFC' ) then
         WRITE  ( UNIT = iuno, FMT = '("PMO=",F6.1,", PMQ=",I2,' //
     +    '", ALSE=",F6.1,", SOB=",F5.1,", SQM=",I2)')
     +   adpsfc_hdr(1),nint(adpsfc_hdr(2)),adpsfc_hdr(3)*0.01,
     +   adpsfc_hdr(4),nint(adpsfc_hdr(5))
      ELSE IF (filo(ifile) .eq. 'SFCSHP' ) then
         WRITE  ( UNIT = iuno, FMT = '("PMO=",F6.1,", PMQ=",I2)')
     +     sfcshp_hdr(1),nint(sfcshp_hdr(2))
      ELSE IF (filo(ifile) .eq. 'SATWND' ) then
         WRITE  ( UNIT = iuno, FMT = '("RFFL=",I3)') nint(satwnd_hdr(1))
      END IF

      IF (nlev .gt. 0 )  THEN

C       Print the level (non-radiance) data for report (if separate
C        mass & wind pieces in PREPBUFR file, these have been merged).
C       --------------------------------------------------------------

         DO lv = 1, nlev  ! loop through the levels
            WRITE  ( UNIT = iuno, FMT = '( 94("-") )' )
            WRITE  ( UNIT = iuno,
     +               FMT = '( 10X,"level ", I4, 2X, A9,A8,A10,5A9 )' )
     +       lv, 'obs      ', 'qmark   ', 'qc_step   ', 'rcode    ',
     +       'fcst     ', 'anal     ', 'oberr    ', 'category '
            WRITE  ( UNIT = iuno, FMT = '( 94("-") )' )
            DO kk = indxvr1(nint(hdr(8))),indxvr2(nint(hdr(8)))
                             ! loop through the variables
               IF( kk .eq. 10 )  ! convert rain rate from mm/s to mm/hr
     +                 evns( 1, lv, :,kk) = evns( 1, lv, :,kk)*3600.0_8
               IF( kk .eq. 15 ) ! convert cloud-top press. from Pa to mb
     +                 evns( 1, lv, :,kk) = evns( 1, lv, :,kk)*0.01
               DO jj = 1, MXR8VN  ! loop through the events
                  IF(min(evns(1,lv,jj,kk),evns(2,lv,jj,kk),
     +                   evns(3,lv,jj,kk),evns(4,lv,jj,kk),
     +                   evns(5,lv,jj,kk),evns(6,lv,jj,kk),
     +                   evns(7,lv,jj,kk))
     +                   .lt. R8BFMS) THEN
                     IF(evns(3,lv,jj,kk) .le. R8BFMS) THEN

C       Copy forecast value, analyzed value, observation error and
C        PREPBUFR category from latest event sequence to all earlier
C        valid event sequences
C       --------------------------------------------------------------

                        evns(5:8,lv,jj,kk) = evns(5:8,lv,1,kk)
                        pcode = evns(3,lv,jj,kk)
                        CALL UFBQCP(11,pcode,QC_STEP)
                     ELSE
                        QC_STEP = '        '
                     END IF
                     WRITE  ( UNIT = iuno,
     + FMT = '( A18,1X,F8.1,2X,F6.0,4X,A8,2X,F5.0,3(1X,F8.1),2X,F7.0)' )
     +               var (kk), ( evns ( ii, lv, jj, kk ), ii = 1, 2 ),
     +               qc_step,  ( evns ( ii, lv, jj, kk ), ii = 4, 8 )
                  END IF
               END DO
            END DO
            IF ( msgtyp .eq. 'ADPUPA' )  THEN
               WRITE  ( UNIT = iuno,
     +    FMT = '( 10X,"BALLOON DRIFT: D-TIME(HRS)=",F7.3,", LAT(N)=",
     +             F6.2,", LON(E)=",F6.2)' ) ( drft ( ii, lv), ii = 1,3)
               IF(min(toth(1,lv),toth(2,lv)) .lt. R8BFMS)
     +          WRITE  ( UNIT = iuno,
     + FMT = '( 10X,"Through PREPRO STEP only: VIRTUAL TEMP (C)=",F8.1,
     +       ", DEWPOINT TEMP (C)=",F8.1)' ) toth (1,lv), toth (2,lv)
            END IF
         END DO
      END IF

      IF (nchn .gt. 0 )  THEN

C       Print the channel (radiance) data for report
C       --------------------------------------------

         WRITE  ( UNIT = iuno, FMT = '( 42("-") )' )
         WRITE  ( UNIT = iuno,
     +            FMT = '("BRIGHTNESS TEMP  channel  observation (K)")')
         WRITE  ( UNIT = iuno, FMT = '( 42("-") )' )
         DO lv = 1, nchn  ! loop through the channels
            IF(min(brts(1,lv),brts(2,lv)) .lt. R8BFMS) THEN
               WRITE  ( UNIT = iuno, FMT = '(19X,I3,7X,F6.2)' )
     +          nint(brts(1,lv)), brts(2,lv)
            END IF
         END DO
      END IF

      IF  ( ierrpb .eq. 0 )  GO TO 10

      print *, '2-All subsets read in and processed'

      STOP      ! All subsets read in and processed

      END

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    READPB
C   PRGMMR: KEYSER           ORG: NP12        DATE: 2002-01-28
C
C ABSTRACT: THIS SUBROUTINE READS IN SUBSETS FROM THE PREPBUFR
C   FILE.  SINCE ACTUAL OBSERVATIONS ARE SPLIT INTO MASS AND
C   WIND SUBSETS, THIS ROUTINE COMBINES THEM INTO A SINGLE
C   REPORT FOR OUTPUT.  IT IS STYLED AFTER BUFRLIB ENTRY POINT
C   READNS, AND IT ONLY REQUIRES THE PREPBUFR FILE TO BE OPENED
C   FOR READING WITH OPENBF.  THE COMBINED REPORT IS RETURNED TO
C   THE CALLER IN COMMON /PREPBC/.  THIS COMMON AREA CONTAINS
C   THE NUMBER OF LEVELS (NON-RADIANCES) IN THE REPORT, THE
C   NUMBER OF CHANNELS (RADIANCES) IN THE REPORT, A ONE
C   DIMENSIONAL ARRAY WITH THE HEADER INFORMATION, A FOUR
C   DIMENSIONAL ARRAY CONTAINING ALL EVENTS FROM ALL
C   NON-RADIANCE OBSERVATIONAL VARIABLE TYPES, AND A TWO
C   DIMENSIONAL ARRAY CONTAINING ALL CHANNELS OF BRIGHTNESS
C   TEMPERATURE FOR REPORTS CONTAINING RADIANCE DATA.
C
C PROGRAM HISTORY LOG:
C ????-??-??  WOOLLEN/GSC ORIGINAL AUTHOR
C ????-??-??  ATOR        STREAMLINED AND DOCUMENTED
C 2001-10-26  KEYSER      FURTHER DOCUMENTATION FOR WEB, UPDATED
C                         TO REPLACE OBSOLETE MESSAGE TYPE "SATBOG"
C                         WITH "QKSWND"; ADDED MORE VARAIBLES AND
C                         MANY OTHER UPDATES
C 2002-01-28  KEYSER      NOW THAT "PCAT" (PRECISION OF TEMPERATURE
C                         OBS.) HAS BEEN ADDED TO PREPBUFR FILE FOR
C                         "AIRCFT" AND "AIRCAR" MESSAGE TYPES, ADDED
C                         THIS TO LISTING
C
C USAGE:    CALL READPB  (LUNIT, MSGTYP, IDATE, IRET)
C   INPUT ARGUMENT LIST:
C     LUNIT    - UNIT NUMBER OF INPUT PREPBUFR FILE
C
C   OUTPUT ARGUMENT LIST:
C     MSGTYP   - BUFR MESSAGE TYPE (CHARACTER)
C     IDATE    - BUFR MESSAGE DATE IN FORM YYYYMMDDHH
C     IRET     - ERROR RETURN CODE (0 - normal return; 1 - the report
C                within COMMON /PREPBC/ contains the last available
C                report from within the prepbufr file; -1 - there are
C                no more reports available from within the PREPBUFR
C                file, all subsets have been processed)
C
C REMARKS: 
C 
C   The header array HDR contains the following list of mnemonics:
C         HDR(1)  Station identification (SID)
C         HDR(2)  Latitude  (YOB)
C         HDR(3)  Longitude (XOB)
C         HDR(4)  Elevation (ELV)
C         HDR(5)  Observation time minus cycle time (DHR)
C         HDR(6)  Reported observation time (RPT)
C         HDR(7)  Indicator whether observation time used to generate
C                  "DHR" was corrected (TCOR)
C         HDR(8)  PREPBUFR report type (TYP)
C         HDR(9) PREPBUFR report subtype (TSB)
C         HDR(10) Input report type (T29)
C         HDR(11) Instrument type (ITP)
C         HDR(12) Report sequence number (SQN)
C         HDR(13) Process number for this MPI run (obtained from
C                  script (PROCN)
C         HDR(14) Satellite identifier (SAID)
C 
C   The 4-D array of non-radiance data, EVNS (ii, lv, jj, kk), is
C    indexed as follows:
C     "ii" indexes the event data types; these consist of:
C         1) Observation       (e.g., POB, ZOB, UOB, VOB, TOB, QOB, PWO)
C         2) Quality mark      (e.g., PQM, ZRM, WQM, TQM, QQM, PWQ)
C         3) Program code      (e.g., PPC, ZPC, WPC, TPC, QPC, PWP)
C         4) Reason code       (e.g., PRC, ZRC, WRC, TRC, QRC, PWR)
C         5) Forecast value    (e.g., PFC, ZFC, UFC, VFC, TFC, QFC, PWF)
C         6) Analyzed value    (e.g., PAN, ZAN, UAN, VAN, TAN, QAN, PWA)
C         7) Observation error (e.g., POE, ZOE, WOE, TOE, QOE, PWO)
C         8) PREPBUFR data level category (CAT)
C     "lv" indexes the levels of the report
C         1) Lowest level
C     "jj" indexes the event stacks
C         1) N'th event
C         2) (N-1)'th event (if present)
C         3) (N-2)'th event (if present)
C                ...
C        10) (N-9)'th event (if present)
C     "kk" indexes the variable types
C         1) Pressure
C         2) Specific humidity
C         3) Temperature
C         4) Height
C         5) U-component wind
C         6) V-component wind
C         7) Wind direction
C         8) Wind speed
C         9) Total precipitable water
C        10) Rain rate
C        11) 1.0 to 0.9 sigma layer precipitable water
C        12) 0.9 to 0.7 sigma layer precipitable water
C        13) 0.7 to 0.3 sigma layer precipitable water
C        14) 0.3 to 0.0 sigma layer precipitable water
C        15) Cloud-top pressure
C        16) Cloud-top temperature
C        17) Total cloud cover
C 
C     Note that the structure of this array is identical to one
C      returned from BUFRLIB routine UFBEVN, with an additional (4'th)
C      dimension to include the 14 variable types into the same
C      array.
C 
C   The 2-D array of radiance data, BRTS (ii, lv), is indexed
C    as follows:
C     "ii" indexes the event data types; these consist of:
C         1) Observation       (CHNM, TMBR)
C     "lv" indexes the channels (CHNM gives channel number)
C 
C     Note that this array is directly returned from BUFRLIB
C      routine UFBINT.
C
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE:  IBM SP
C
C$$$
      SUBROUTINE READPB  ( lunit, msgtyp, idate, iret )

      REAL*8    R8BFMS

      PARAMETER ( R8BFMS = 9.0E08 )! Missing value threshhold for BUFR
      PARAMETER ( NHR8PM = 14 )    ! Actual number of BUFR parameters
                                   !  in header
      PARAMETER ( MXR8PM =  8 )    ! Maximum number of BUFR parameters
                                   !  in level data (non-radiance) or
                                   !  in channel data (radiance)
      PARAMETER ( MXR8LV = 255 )   ! Maximum number of BUFR levels/
                                   !                        channels
      PARAMETER ( MXR8VN = 10 )    ! Max. number of BUFR event sequences
                                   !  (non-radiance reports)
      PARAMETER ( MXR8VT = 17)     ! Max. number of BUFR variable types
                                   !  (non-radiance reports)
      PARAMETER ( MXSTRL = 80 )    ! Maximum size of a string

      REAL*8    hdr ( NHR8PM ), qkswnd_hdr(3), aircar_hdr(6),
     +          aircft_hdr(4),  adpupa_hdr(1), goesnd1_hdr(6),
     +          goesnd2_hdr(MXR8LV), adpsfc_hdr(5), sfcshp_hdr(2),
     +          satwnd_hdr(1), evns ( MXR8PM, MXR8LV, MXR8VN, MXR8VT ),
     +          drft ( 3, MXR8LV ), brts ( MXR8PM, MXR8LV ),
     +          toth ( 2, MXR8LV )
      REAL*8    hdr2 ( NHR8PM ), qkswnd_hdr2(3), aircar_hdr2(6),
     +          aircft_hdr2(4),  adpupa_hdr2(1), goesnd1_hdr2(6),
     +          goesnd2_hdr2(MXR8LV), adpsfc_hdr2(5), sfcshp_hdr2(2),
     +          satwnd_hdr2(1), evns2 ( MXR8PM, MXR8LV, MXR8VN, MXR8VT),
     +          drft2 ( 3, MXR8LV ), toth2 ( 2, MXR8LV )
      REAL*8    hdr_save  ( NHR8PM ), adpsfc_hdr_save(5),
     +          sfcshp_hdr_save(2),
     +          evns_save ( MXR8PM, MXR8LV, MXR8VN, MXR8VT ),
     +          drft_save ( 3, MXR8LV ), toth_save ( 2, MXR8LV )
      REAL*8    r8sid, r8sid2, pob1, pob2

      REAL      p ( MXR8LV )

      INTEGER  indr ( MXR8LV ), indxvr1 ( 100:299 ), indxvr2 ( 100:299 )

      CHARACTER    msgtyp*8, msgtyp_process*8
      CHARACTER*8  csid, csid2, msgtp2
      CHARACTER*(MXSTRL)      head1, head2, ostr ( MXR8VT )

      LOGICAL         match, seq_match, single_msgtyp, radiance,
     +                not_radiance

      EQUIVALENCE     ( r8sid, csid ), ( r8sid2, csid2 )

      COMMON  / PREPBC /      hdr, evns, drft, toth, brts, qkswnd_hdr,
     +  aircar_hdr, aircft_hdr, adpupa_hdr, goesnd1_hdr, goesnd2_hdr,
     +  adpsfc_hdr, sfcshp_hdr, satwnd_hdr, indxvr1, indxvr2, nlev, nchn

      DATA   head1
     +        / 'SID YOB XOB ELV DHR RPT TCOR TYP TSB T29 ITP SQN' /
      DATA   head2
     +        / 'PROCN SAID' /
      DATA   ostr
     +        / 'POB PQM PPC PRC PFC PAN POE CAT',
     +          'QOB QQM QPC QRC QFC QAN QOE CAT',
     +          'TOB TQM TPC TRC TFC TAN TOE CAT',
     +          'ZOB ZQM ZPC ZRC ZFC ZAN ZOE CAT',
     +          'UOB WQM WPC WRC UFC UAN WOE CAT',
     +          'VOB WQM WPC WRC VFC VAN WOE CAT',
     +          'DDO DFQ DFP DFR NULL NULL NULL CAT',
     +          'FFO DFQ DFP DFR NULL NULL NULL CAT',
     +          'PWO PWQ PWP PWR PWF PWA PWE CAT',
     + 'REQ6 REQ6_QM REQ6_PC REQ6_RC REQ6_FC REQ6_AN REQ6_OE CAT',
     +          'PW1O PW1Q PW1P PW1R PW1F PW1A PW1E CAT',
     +          'PW2O PW2Q PW2P PW2R PW2F PW2A PW2E CAT',
     +          'PW3O PW3Q PW3P PW3R PW3F PW3A PW3E CAT',
     +          'PW4O PW4Q PW4P PW4R PW4F PW4A PW4E CAT',
     + 'CDTP CDTP_QM CDTP_PC CDTP_RC CDTP_TC CDTP_AN CDTP_OE CAT',
     +          'GCDTT NULL NULL NULL NULL NULL NULL CAT' ,
     +          'TOCC  NULL NULL NULL NULL NULL NULL CAT' /
      DATA            match / .true. /, seq_match / .false. /

      SAVE            match, msgtp2, idate2
C-----------------------------------------------------------------------


      single_msgtyp = .false.   ! set to true if you want to process
                                !  reports from only ONE message type


      msgtyp_process = 'ADPUPA' ! if single_msgtyp=T, then only
                                !  reports in this msgtyp will be
                                !  processed

 1000 continue

      iret = 0

C       If the previous call to this subroutine did not yield matching
C        mass and wind subsets, then BUFRLIB routine READNS is already
C        pointing at an unmatched subset.  Otherwise, call READNS to
C        advance the subset pointer to the next subset.
C       --------------------------------------------------------------
      IF  ( match )  THEN
          CALL READNS(lunit,msgtyp,idate,jret) ! BUFRLIB routine to
                                               !  read in next subset
          IF  ( jret .ne. 0 )  THEN
             iret = -1 ! there are no more subsets in the PREPBUFR file
                       !  and all subsets have been processed
             RETURN
          END IF
          if(single_msgtyp .and. msgtyp.ne.msgtyp_process) go to 1000
      ELSE
          msgtyp = msgtp2
          idate = idate2
      END IF

C       Read the report header based on mnemonic string "head" and
C        transfer 1 for 1 into array HDR for the subset that is
C        currently being pointed to.
C       -----------------------------------------------------------

                       !BUFRLIB routine to read specific info from
                       ! report (no event info though)
      CALL UFBINT(lunit,hdr,           12,1,jret,head1)
      CALL UFBINT(lunit,hdr(13),NHR8PM-12,1,jret,head2)

C       Read header information that is specific to the various
C        message types
C       -------------------------------------------------------

      qkswnd_hdr=10E10
      aircar_hdr=10E10
      aircft_hdr=10E10
      adpupa_hdr=10E10
      goesnd1_hdr=10E10
      goesnd2_hdr=10E10
      adpsfc_hdr=10E10
      sfcshp_hdr=10E10
      satwnd_hdr=10E10

      IF( msgtyp .eq. 'QKSWND')  THEN
         CALL UFBINT(lunit,qkswnd_hdr,3,1,jret,'CTCN ATRN SPRR ')
      ELSE IF( msgtyp .eq. 'AIRCAR')  THEN
         CALL UFBINT(lunit,aircar_hdr,6,1,jret,
     +    'PCAT POAF TRBX10 TRBX21 TRBX32 TRBX43 ')
      ELSE IF( msgtyp .eq. 'AIRCFT')  THEN
         CALL UFBINT(lunit,aircft_hdr,4,1,jret,'RCT PCAT POAF DGOT ')
      ELSE IF( msgtyp .eq. 'ADPUPA')  THEN
         CALL UFBINT(lunit,adpupa_hdr,1,1,jret,'SIRC ')
      ELSE IF( msgtyp .eq. 'GOESND' .or.  msgtyp .eq. 'SATEMP' )  THEN
         CALL UFBINT(lunit,goesnd1_hdr,6,1,jret,
     +               'ACAV ELEV SOEL OZON TMSK CLAM')
         IF( msgtyp .eq. 'GOESND' ) CALL UFBINT(lunit,goesnd2_hdr,1,
     +                                          MXR8LV,jret,'PRSS ')
      ELSE IF( msgtyp .eq. 'ADPSFC')  THEN
         CALL UFBINT(lunit,adpsfc_hdr,5,1,jret,'PMO PMQ ALSE SOB SQM ')
      ELSE IF( msgtyp .eq. 'SFCSHP')  THEN
         CALL UFBINT(lunit,sfcshp_hdr,2,1,jret,'PMO PMQ ')
      ELSE IF( msgtyp .eq. 'SATWND')  THEN
         CALL UFBINT(lunit,satwnd_hdr,1,1,jret,'RFFL ')
      END IF

C       From PREPBUFR report type, determine if this report contains
C        only non-radiance data, only radiance data, or both
C       ------------------------------------------------------------

      radiance = (  nint(hdr(8)) .eq. 102 .or.
     +            ( nint(hdr(8)) .ge. 160 .and. nint(hdr(8)) .le. 179 ))
      not_radiance = ( .not.radiance .or. 
     +          ( nint(hdr(8)) .ge. 160 .and. nint(hdr(8)) .le. 163 )
     +     .or. ( nint(hdr(8)) .ge. 170 .and. nint(hdr(8)) .le. 173 ))

      EVNS = 10E10
      BRTS = 10E10
      nlev = 0
      nchn = 0

      IF ( not_radiance )  THEN

C       For reports with non-radiance data, read the report level data
C        for the variable types based on the variable-specific mnemonic
C        string "ostr(kk)" and transfer all levels, all events 1 for 1
C        into array EVNS for the subset that is currently being pointed
C        to.
C       ---------------------------------------------------------------


         DO kk = indxvr1(nint(hdr(8))),indxvr2(nint(hdr(8)))
                                 ! loop through the variables
             CALL UFBEVN(lunit,evns(1,1,1,kk),MXR8PM,MXR8LV,MXR8VN,nlev,
     +        ostr(kk)) ! BUFRLIB routine to read specific info from
                        ! report, accounting for all possible events
         END DO

         drft=10E10
         toth=10E10
         IF( msgtyp .eq. 'ADPUPA')  THEN

C       Read balloon drift level data for message type ADPUPA
C       -----------------------------------------------------

            CALL UFBINT(lunit,drft,3,MXR8LV,nlevd,'HRDR YDR XDR ')
            if(nlevd .ne. nlev)  CALL ERREXIT(22)

C       Read virtual temperature and dewpoint temperature (through
C        "PREPRO" step only) data for message type ADPUPA
C       ----------------------------------------------------------

            CALL UFBINT(lunit,toth,2,MXR8LV,nlevt,'TVO TDO ')
            if(nlevt .ne. nlev)  CALL ERREXIT(33)

         END IF

      END IF

      IF ( radiance)  THEN

C       For reports with radiance data, read the report channel data
C        for the variable types based on the mnemonic string
C        "CHNM TMBR" and transfer all channels 1 for 1 into array BRTS
C        for the subset that is currently being pointed to.
C       --------------------------------------------------------------

          CALL UFBINT(lunit,brts,MXR8PM,MXR8LV,nchn,'CHNM TMBR')

      END IF

 2000 CONTINUE

C      Now, advance the subset pointer to the next subset and read
C       its report header.
C      -----------------------------------------------------------

      CALL READNS(lunit,msgtp2,idate2,jret)
      IF  ( jret .ne. 0 )  THEN
         iret = 1 ! there are no more subsets in the PREPBUFR file
                  !  but we must still return and process the
                  !  previous report in common /PREPBC/
         RETURN
      END IF
      if(single_msgtyp .and. msgtp2.ne.msgtyp_process) go to 2000
      CALL UFBINT(lunit,hdr2,           12,1,jret,head1)
      CALL UFBINT(lunit,hdr2(13),NHR8PM-12,1,jret,head2)

C       Next, check whether this subset and the previous one are
C        matching mass and wind components for a single "true" report.
C       --------------------------------------------------------------

      match = .true.
      seq_match = .false.

      IF ( max( hdr (12), hdr2 (12) ) .le. R8BFMS )  THEN
         IF ( hdr (13)  .gt. R8BFMS )  hdr (13) = 0.
         IF ( hdr2(13)  .gt. R8BFMS )  hdr2(13) = 0.
         seq  = (hdr (13) * 1000000.) + hdr (12) ! combine seq. nunmber
         seq2 = (hdr2(13) * 1000000.) + hdr2(12) !  and MPI proc. number
         IF  ( seq .ne. seq2 )  THEN
            match = .false.
            RETURN  ! the two process/sequence numbers do not match,
                    !  return and process information in common /PREPBC/
         ELSE
            seq_match = .true.  ! the two process/sequence numbers match
         END IF
      END IF

      IF  ( msgtyp .ne. msgtp2 )  THEN
          match = .false.
          RETURN   ! the two message types do not match, return and
                   !  process information in common /PREPBC/
      END IF

      IF ( .not. seq_match )  THEN

C       The two process/sequence numbers are missing, so as a last
C        resort must check whether this subset and the previous one
C        have identical values for SID, YOB, XOB, ELV, and DHR in
C        which case they match.
C       -----------------------------------------------------------

         r8sid  = hdr (1)
         r8sid2 = hdr2(1)
         IF  ( csid .ne. csid2 )  THEN
             match = .false.
             RETURN   ! the two report id's do not match, return and
                      !  process information in common /PREPBC/
         END IF

         DO ii = 5, 2, -1
            IF  ( hdr (ii) .ne. hdr2 (ii) )  THEN
               match = .false.
               RETURN   ! the two headers do not match, return and
                        !  process information in common /PREPBC/
            END IF
         END DO
      END IF

C       The two headers match, read level data for the second subset.
C        (Note: This can only happen for non-radiance reports)
C       -------------------------------------------------------------

      DO kk = indxvr1(nint(hdr(8))),indxvr2(nint(hdr(8)))
                                 ! loop through the variables
          CALL UFBEVN(lunit,evns2 (1,1,1,kk),MXR8PM,MXR8LV,MXR8VN,
     +     nlev2,ostr(kk))
      ENDDO

      drft2=10E10
      toth2=10E10
      IF( msgtp2 .eq. 'ADPUPA')  THEN

C       Read balloon drift level data for message type ADPUPA
C        for the second subset
C       -----------------------------------------------------

         CALL UFBINT(lunit,drft2,3,MXR8LV,nlev2d,'HRDR YDR XDR ')
         if(nlev2d .ne. nlev2)  CALL ERREXIT(22)

C       Read virtual temperature and dewpoint temperature (through
C        "PREPRO" step only) data for message type ADPUPA for the
C        second subset
C       ----------------------------------------------------------

         CALL UFBINT(lunit,toth2,2,MXR8LV,nlev2t,'TVO TDO ')
         if(nlev2t .ne. nlev2)  CALL ERREXIT(33)

      END IF

C       Read header information that is specific to the various
C        message types for the second supset
C       -------------------------------------------------------

      qkswnd_hdr2=10E10
      aircar_hdr2=10E10
      aircft_hdr2=10E10
      adpupa_hdr2=10E10
      goesnd1_hdr2=10E10
      goesnd2_hdr2=10E10
      adpsfc_hdr2=10E10
      sfcshp_hdr2=10E10
      satwnd_hdr2=10E10
      IF( msgtp2 .eq. 'QKSWND')  THEN
         CALL UFBINT(lunit,qkswnd_hdr2,3,1,jret,'CTCN ATRN SPRR ')
      ELSE IF( msgtp2 .eq. 'AIRCAR')  THEN
         CALL UFBINT(lunit,aircar_hdr2,6,1,jret,
     +    'PCAT POAF TRBX10 TRBX21 TRBX32 TRBX43 ')
      ELSE IF( msgtp2 .eq. 'AIRCFT')  THEN
         CALL UFBINT(lunit,aircft_hdr2,4,1,jret,'RCT PCAT POAF DGOT ')
      ELSE IF( msgtp2 .eq. 'ADPUPA')  THEN
         CALL UFBINT(lunit,adpupa_hdr2,1,1,jret,'SIRC ')
      ELSE IF( msgtp2 .eq. 'GOESND' .or.  msgtp2 .eq. 'SATEMP' )  THEN
         CALL UFBINT(lunit,goesnd1_hdr2,6,1,jret,
     +               'ACAV ELEV SOEL OZON TMSK CLAM')
         IF( msgtp2 .eq. 'GOESND' ) CALL UFBINT(lunit,goesnd2_hdr2,1,
     +                                          MXR8LV,jret,'PRSS ')
      ELSE IF( msgtp2 .eq. 'ADPSFC')  THEN
         CALL UFBINT(lunit,adpsfc_hdr2,5,1,jret,'PMO PMQ ALSE SOB SQM ')
      ELSE IF( msgtp2 .eq. 'SFCSHP')  THEN
         CALL UFBINT(lunit,sfcshp_hdr2,2,1,jret,'PMO PMQ ')
      ELSE IF( msgtp2 .eq. 'SATWND')  THEN
         CALL UFBINT(lunit,satwnd_hdr2,1,1,jret,'RFFL ')
      END IF

      IF (nint(hdr(8)) .ge. 280)  THEN

C       If this is a surface report, the wind subset precedes the
C        mass subset - switch the subsets around in order to combine
C        the surface pressure properly
C       ------------------------------------------------------------

         evns_save = evns2
         evns2 = evns
         evns = evns_save
         hdr_save = hdr2
         hdr2 = hdr
         hdr = hdr_save
         adpsfc_hdr_save = adpsfc_hdr2
         adpsfc_hdr2 = adpsfc_hdr
         adpsfc_hdr = adpsfc_hdr_save
         sfcshp_hdr_save = sfcshp_hdr2
         sfcshp_hdr2 = sfcshp_hdr
         sfcshp_hdr = sfcshp_hdr_save
      END IF

C       Combine the message type-specific header data for the two
C        matching subsets into a single array. 
C       ---------------------------------------------------------

      qkswnd_hdr(:) = min(qkswnd_hdr(:),qkswnd_hdr2(:))
      aircar_hdr(:) = min(aircar_hdr(:),aircar_hdr2(:))
      aircft_hdr(:) = min(aircft_hdr(:),aircft_hdr2(:))
      adpupa_hdr(:) = min(adpupa_hdr(:),adpupa_hdr2(:))
      goesnd1_hdr(:) = min(goesnd1_hdr(:),goesnd1_hdr2(:))
      goesnd2_hdr(:) = min(goesnd2_hdr(:),goesnd2_hdr2(:))
      adpsfc_hdr(:) = min(adpsfc_hdr(:),adpsfc_hdr2(:))
      sfcshp_hdr(:) = min(sfcshp_hdr(:),sfcshp_hdr2(:))
      satwnd_hdr(:) = min(satwnd_hdr(:),satwnd_hdr2(:))

C       Combine the data for the two matching subsets into a single 4-D
C        array.  Do this by merging the EVNS2 array into the EVNS array.
C       ----------------------------------------------------------------

      LOOP1: DO lv2 = 1, nlev2  ! loop through the levels of subset 2
         DO lv1 = 1, nlev  ! loop through the levels of subset 1
            pob1 = evns  ( 1, lv1, 1, 1 )
            pob2 = evns2 ( 1, lv2, 1, 1 )
            IF  ( pob1 .eq. pob2 )  THEN

C       This pressure level from the second subset also exists in the
C        first subset, so overwrite any "missing" piece of data for
C        this pressure level in the first subset with the corresponding
C        piece of data from the second subset (since this results in no
C        net loss of data!).
C       ---------------------------------------------------------------

               DO kk = indxvr1(nint(hdr(8))),indxvr2(nint(hdr(8)))
                                 ! loop through the variables
                  DO jj = 1,MXR8VN  ! loop through the events
                     DO ii = 1,MXR8PM-1 ! loop through BUFR parameters
                                        ! skip the CAT parameter
                        IF ( evns (ii,lv1,jj,kk ) .gt. R8BFMS ) THEN
                           evns  ( ii, lv1, jj, kk ) =
     +                     evns2 ( ii, lv2, jj, kk )
                           IF ( evns2 (ii,lv2,jj,kk).le.R8BFMS) THEN

C       If any piece of data in the first subset is missing and the
C        corresponding piece of data is NOT missing in the second
C        subset, also overwrite the PREPBUFR category parameter in the
C        first subset with that from the second subset regardless of
C        its value in either subset
C       --------------------------------------------------------------

                              evns  (  8, lv1, jj, kk ) =
     +                        evns2 (  8, lv2, jj, kk )
                           END IF
                        END IF
                     END DO
                  END DO
               END DO
               CYCLE LOOP1
            ELSE IF  (  ( pob2 .gt. pob1 )  .or.
     +                  ( lv1  .eq. nlev )  )  THEN

C       Either all remaining pressure levels within the first subset
C        are less than this pressure level from the second subset
C        (since levels within each subset are guaranteed to be in
C        descending order wrt pressure!) *OR* there are more total
C        levels within the second subset than in the first subset.
C        In either case, we should now add this second subset level
C        to the end of the EVNS array.
C       ------------------------------------------------------------

               nlev = nlev + 1
               evns  ( :, nlev, :, : ) = evns2 ( :, lv2,  :, : )
               drft  ( :, nlev ) = drft2 ( :, lv2 )
               toth  ( :, nlev ) = toth2 ( :, lv2 )
               CYCLE LOOP 1
            END IF
         END DO
      END DO  LOOP1

C       Sort combined report according to decreasing pressure.
C       ------------------------------------------------------

      p(:) = evns  ( 1, :, 1, 1)
      if(nlev .gt. 0)  call indexf(nlev,p,indr)
      evns_save = evns
      drft_save = drft
      toth_save = toth
      do i = 1,nlev
         j = indr(nlev-i+1)
         evns ( :, i, :, :) = evns_save ( :, j, :, :)
         drft ( :, i) = drft_save ( :, j)
         toth ( :, i) = toth_save ( :, j)
      end do

C       Return to calling program with a complete report in common
C        /PREPBC/.
C       ----------------------------------------------------------

      RETURN

      END

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    INDEXF      GENERAL SORT ROUTINE FOR REAL ARRAY
C   PRGMMR: D. A. KEYSER     ORG: NP22       DATE: 1995-05-30
C
C ABSTRACT: USES EFFICIENT SORT ALGORITHM TO PRODUCE INDEX SORT LIST
C   FOR AN REAL ARRAY.  DOES NOT REARRANGE THE FILE.
C
C PROGRAM HISTORY LOG:
C 1993-06-05  R  KISTLER -- FORTRAN VERSION OF C-PROGRAM
C 1995-05-30  D. A. KEYSER - TESTS FOR < 2 ELEMENTS IN SORT LIST,
C             IF SO RETURNS WITHOUT SORTING (BUT FILLS INDX ARRAY)
C
C USAGE:    CALL INDEXF(N,ARRIN,INDX)
C   INPUT ARGUMENT LIST:
C     N        - SIZE OF ARRAY TO BE SORTED
C     ARRIN    - REAL ARRAY TO BE SORTED
C
C   OUTPUT ARGUMENT LIST:
C     INDX     - ARRAY OF POINTERS GIVING SORT ORDER OF ARRIN IN
C              - ASCENDING ORDER {E.G., ARRIN(INDX(I)) IS SORTED IN
C              - ASCENDING ORDER FOR ORIGINAL I = 1, ... ,N}
C
C REMARKS: NONE.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE:  IBM-SP
C
C$$$
      SUBROUTINE INDEXF(N,ARRIN,INDX)
      INTEGER  INDX(N)
      REAL  ARRIN(N)
      DO J = 1,N
         INDX(J) = J
      ENDDO
C MUST BE > 1 ELEMENT IN SORT LIST, ELSE RETURN
      IF(N.LE.1)  RETURN
      L = N/2 + 1
      IR = N
   33 CONTINUE
      IF(L.GT.1)  THEN
         L = L - 1
         INDXT = INDX(L)
         AA = ARRIN(INDXT)
      ELSE
         INDXT = INDX(IR)
         AA = ARRIN(INDXT)
         INDX(IR) = INDX(1)
         IR = IR - 1
         IF(IR.EQ.1)  THEN
            INDX(1) = INDXT
            RETURN
         END IF
      END IF
      I = L
      J = L * 2
   30 CONTINUE
      IF(J.LE.IR)  THEN
         IF(J.LT.IR)  THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))  J = J + 1
         END IF
         IF(AA.LT.ARRIN(INDX(J)))  THEN
            INDX(I) = INDX(J)
            I = J
            J = J + I
         ELSE
            J = IR + 1
         END IF
      END IF
      IF(J.LE.IR)  GO TO 30
      INDX(I) = INDXT
      GO TO 33
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      BLOCK DATA

      PARAMETER ( NHR8PM = 14 )    ! Actual number of BUFR parameters
                                   !  in header
      PARAMETER ( MXR8PM =  8 )    ! Maximum number of BUFR parameters
                                   !  in level data (non-radiance) or
                                   !  in channel data (radiance)
      PARAMETER ( MXR8LV = 255 )   ! Maximum number of BUFR levels/
                                   !                        channels
      PARAMETER ( MXR8VN = 10 )    ! Max. number of BUFR event sequences
                                   !  (non-radiance reports)
      PARAMETER ( MXR8VT = 17)     ! Max. number of BUFR variable types
                                   !  (non-radiance reports)

      REAL*8    hdr ( NHR8PM ), qkswnd_hdr(3), aircar_hdr(6),
     +          aircft_hdr(4),  adpupa_hdr(1), goesnd1_hdr(6),
     +          goesnd2_hdr(MXR8LV), adpsfc_hdr(5), sfcshp_hdr(2),
     +          satwnd_hdr(1), evns ( MXR8PM, MXR8LV, MXR8VN, MXR8VT ),
     +          drft ( 3, MXR8LV ), brts ( MXR8PM, MXR8LV ),
     +          toth ( 2, MXR8LV )

      INTEGER   indxvr1 ( 100:299 ), indxvr2 ( 100:299 )

      COMMON  / PREPBC /      hdr, evns, drft, toth, brts, qkswnd_hdr,
     +  aircar_hdr, aircft_hdr, adpupa_hdr, goesnd1_hdr, goesnd2_hdr,
     +  adpsfc_hdr, sfcshp_hdr, satwnd_hdr, indxvr1, indxvr2, nlev, nchn

      DATA      nlev /0/, nchn /0/

      DATA      indxvr1
     + /   20*1,    30*1,   10,  15,   9,    3*1,    4*11,    40*1 ,
C         100-119  120-149  150  151  152  153-155  156-159  160-199
     +     20*1,    40*1,    25*1,    1,   1,   13*1    /
C         200-219  220-259  260-284  285  286  287-299

      DATA      indxvr2
     + /   20*4,    30*8,   10,  17,   9,    3*6,    4*14,    40*6 ,
C         100-119  120-149  150  151  152  153-155  156-159  160-299
     +     20*6,    40*8,    25*6,    8,   8,   13*6    /
C         200-219  220-259  260-284  285  286  287-299

      END



