! $Id: bufr2ropp_mod.f90 6490 2020-09-10 17:36:27Z idculv $

MODULE bufr2ropp

!****m* bufr2ropp/bufr2ropp *
!
! NAME
!   bufr2ropp    (bufr2ropp_mod.f90)
!
! SYNOPSIS
!   Module defining fixed values & common subroutines/functions for the
!   bufr2ropp main program (ECMWF or MetDB versions)
!
!   USE bufr2ropp
!
! USED BY
!   bufr2ropp  (bufr2ropp_ec or bufr2ropp_mo)
!
! AUTHOR
!   Met Office, Exeter, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

! Modules

  USE messages
  USE typesizes, dp => EightByteReal

! Fixed values

  INTEGER, PARAMETER :: RODescr = 310026    ! Table D descriptor for RO

  INTEGER, PARAMETER :: nb      =   500000  ! Max. length of BUFR message (bytes)
  INTEGER, PARAMETER :: nd      =   200000  ! Max. no. of descriptors
  INTEGER, PARAMETER :: no      =        1  ! Max. no. of observations
  INTEGER, PARAMETER :: nv      =    nd*no  ! Max. no. of data values

! Default (missing) data values

  INTEGER,  PARAMETER :: NMDFV = -9999999     ! Integer missing data flag value (MetDB)
  INTEGER,  PARAMETER :: NVIND = 2147483647   ! Integer missing data flag value (ECMWF)
  REAL,     PARAMETER :: RMDFV = -9999999.0   ! Real    missing data flag value (MetDB)
  REAL(dp), PARAMETER :: RVIND = 1.7E38_dp    ! Real    missing data flag value (ECMWF)

! MetDB I/O interface

  INTEGER, PARAMETER  :: Input  = 1        ! I/O input mode (r)

CONTAINS
!-------------------------------------------------------------------------

SUBROUTINE Usage()

!****s* bufr2ropp/Usage *
!
! NAME
!   Usage
!
! SYNOPIS
!   USE bufr2ropp
!   CALL Usage()
!
! INPUTS
!   None
!
! OUTPUTS
!   Summary usage text to stdout
!
! CALLED BY
!   GetOptions
!
! DESCRIPTION
!   Prints a summary of program usage (help) to stdout.
!
! AUTHOR
!   Met Office, Exeter, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

  PRINT *, 'Purpose:'
  PRINT *, '  Decode RO BUFR messages to ROPP (netCDF) format.'
  PRINT *, 'Usage:'
  PRINT *, ' > bufr2ropp bufr_file [bufr_file...] [-o ropp_file] [-m] [-a]'
  PRINT *, '                       [-f first] [-n number] [-d] [-h|?] [-v]'
  PRINT *, 'Input:'
  PRINT *, '  One or more files containing BUFR messages or WMO bulletins.'
  PRINT *, '  Any non-RO messages will be skipped.'
  PRINT *, 'Output:'
  PRINT *, '  One or more ROPP files.'
  PRINT *, 'Options:'
  PRINT *, '  -o ROPP netCDF output file name. Mandatory argument if used'
  PRINT *, '  -m write ROPP output as a multifile, i.e. if there'
  PRINT *, '     is more than one RO BUFR message in the input, they'
  PRINT *, '     are decoded into a single netCDF output file'
  PRINT *, '  -a append new profiles to file specified by -o. (-a implies -m)'
  PRINT *, '  -f specify the first RO message to decode'
  PRINT *, '  -n specify the number of RO messages to decode'
  PRINT *, '  -d outputs additional diagnostics to stdout'
  PRINT *, '  -h this help'
  PRINT *, '  -v version information'
  PRINT *, 'Defaults:'
  PRINT *, '  Input  file name : ropp.bufr'
  PRINT *, '  Output file name : from (first) occultation ID <occid>.nc'
  PRINT *, '  Output mode      : one netCDF output file per input RO BUFR message'
  PRINT *, '  Append mode      : exising file is overwritten'
  PRINT *, '  first   : 1'
  PRINT *, '  number  : 999999'
  PRINT *, 'See bufr2ropp(1) for details.'
  PRINT *, ''
END SUBROUTINE Usage
!-------------------------------------------------------------------------

SUBROUTINE GetOptions ( BUFRdsn,      & ! (out)
                        nfiles,       & ! (out)
                        ROPPdsn,      & ! (out)
                        multi,        & ! (out)
                        newfile,      & ! (out)
                        fMsgToDecode, & ! (out)
                        nMsgToDecode )  ! (out)

!****s* bufr2ropp/GetOptions *
!
! NAME
!   GetOptions
!
! SYNOPSIS
!   Get command line information & options or set defaults
!
!    USE bufr2ropp
!    CHARACTER (LEN=100), DIMENSION(:), ALLOCATABLE :: BUFRdsn
!    CHARACTER (LEN=100) :: roppdsn
!    INTEGER :: nfiles, fMsgToDecode, nMsgToDecode
!    LOGICAL :: multi, newfile
!    nfiles=IARGC()
!    ALLOCATE ( bufrdsn(nfiles) )
!    CALL getoptions ( bufrdsn, nfiles, roppdsn, &
!                      multi, newfile, fMsgToDecode, nMsgToDecode )
!   On the command line:
!    >  bufr2ropp bufr_file [...bufr_file...] [-o ropp_file]
!                           [-m] [-a]   [-f first] [-n number]
!                           [-d] [-h|?] [-v]
!
! INPUTS
!   None
!
! OUTPUTS
!   BUFRdsn       chr  BUFR input file name(s) (required)
!   nfiles        int  No. of BUFR input files
!   ROPPdsn       chr  ROPP output file name (optional)
!   multi         log  Multifile output (netCDF only) (default: single files)
!   newfile       log  .T. to start a new file, else append (default: new)
!   fMsgToDecode  int  First  RO message  to decode (default: 1)
!   nMsgToDecode  int  No. of RO messages to decode (default: all)
!
! CALLS
!   Usage
!   message
!   GETARG
!   IARGC
!
! CALLED BY
!   bufr2ropp
!
! DESCRIPTION
!   Provides a command line interface for the BUFR-to-ROPP
!   decoder application. See comments for main program bufr2ropp
!   for the command line details.
!
! SEE ALSO
!   bufr2ropp(1)
!
! AUTHOR
!   Met Office, Exeter, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

  IMPLICIT NONE

! Fixed parameters

  CHARACTER (LEN=*), PARAMETER :: Defdsn = "ropp.bufr" ! default input file name

! Argument list parameters

  CHARACTER (LEN=*), INTENT(OUT) :: BUFRdsn(:)    ! BUFR (input)  dataset name(s)
  CHARACTER (LEN=*), INTENT(OUT) :: ROPPdsn       ! ROPP (output) dataset name
  LOGICAL,           INTENT(OUT) :: multi         ! Single/Multifile output flag
  LOGICAL,           INTENT(OUT) :: newfile       ! New/Append output flag
  INTEGER,           INTENT(OUT) :: nfiles        ! No. of BUFR input files
  INTEGER,           INTENT(OUT) :: fMsgToDecode  ! First  ROPP message  to decode
  INTEGER,           INTENT(OUT) :: nMsgToDecode  ! No. of ROPP messages to decode

! Local parameters

  CHARACTER (LEN=256) :: carg    ! command line argument
  CHARACTER (LEN=256) :: Value   ! value extracted from argument
  INTEGER             :: ia, n   ! counters
  INTEGER             :: Narg    ! number of command line arguments
  INTEGER             :: ierr    ! error status code

  ! Some compilers may need the following declaration to be commented out
  INTEGER :: IARGC

!-------------------------------------------------------------
! 1. Initialise
!-------------------------------------------------------------

  BUFRdsn(:)   = Defdsn
  nfiles       = 0
  ROPPdsn      = " "
  multi        = .FALSE.
  newfile      = .TRUE.
  fMsgToDecode = 1
  nMsgToDecode = 999999

!-------------------------------------------------------------
! 2. Loop over all command line arguments.
!    If a switch has a trailing blank, then we need to get
!    the next string as its argument.
!-------------------------------------------------------------

  ia   = 1
  narg = IARGC()

  DO WHILE ( ia <= Narg )

    CALL GETARG ( ia, carg )
    IF ( carg(1:1) == "?" .OR. &
         carg(1:6) == "--help"    ) carg = "-h"
    IF ( carg(1:9) == "--version" ) carg = "-v"

    IF ( carg(1:1) == "-" ) THEN   ! is this an option introducer?
                                   ! If so, which one?
      SELECT CASE (carg(2:2))

        CASE ("a","A")             ! Append (multifile) output requested
          newfile = .FALSE.
          multi   = .TRUE.

        CASE ("d","D")             ! debug/diagnostics wanted
          msg_MODE = VerboseMode

        CASE ("f","F" )            ! Start Message (first)
          Value = carg(3:)
          IF ( Value == " " ) THEN
            ia = ia + 1
            CALL GETARG ( ia, Value )
          END IF
          READ ( Value, *, IOSTAT=ierr ) n
          IF ( ierr == 0 .AND. &
                n    >  1 ) fMsgToDecode = n

        CASE ("h","H")             ! Help wanted
          CALL Usage()
          CALL EXIT(msg_exit_ok)

        CASE ("m","M")             ! Multifile output requested
          multi = .TRUE.

        CASE ("n","N" )            ! Number of Messages
          Value = carg(3:)
          IF ( Value == " " ) THEN
            ia = ia + 1
            CALL GETARG ( ia, Value )
          END IF
          READ ( Value, *, IOSTAT=ierr ) n
          IF ( ierr == 0 .AND. &
                n   >  0 ) nMsgToDecode = n

        CASE ("o","O" )            ! Output file
          ROPPdsn = carg(3:)
          IF ( ROPPdsn == " " ) THEN
            ia = ia + 1
            CALL GETARG ( ia, ROPPdsn )
          END IF

        CASE ("v","V" )            ! Only program version ID wanted
          CALL version_info()
          CALL EXIT(msg_exit_ok)

! Ignore anything else

        CASE DEFAULT               ! unknown option
      END SELECT

    ELSE
      nfiles = nfiles + 1
      BUFRdsn(nfiles) = carg       ! not an option - must be an input name
    END IF

    ia = ia + 1
  END DO         ! argument loop

  IF ( nfiles == 0 ) nfiles = 1    ! No input files - try a default name

END SUBROUTINE GetOptions
!----------------------------------------------------------------------------

SUBROUTINE ConvertBUFRtoROPP ( Values, & ! (in)
                               NFreq,  & ! (in)
                               ROdata, & ! (inout)
                               ierr )    ! (out)
!
!****s* bufr2ropp/ConvertBUFRtoROPP *
!
! NAME
!   ConvertBUFRtoROPP
!
! SYNOPSIS
!   Convert BUFR data to ROPP specification
!
!    USE ropp_io_types
!    USE bufr2ropp
!    TYPE (ROprof) ROdata
!    REAL*8  :: values(nv)
!    INTEGER :: nfreq, ierr
!    CALL ConvertBUFRtoROPP ( values, nfreq, ROdata, ierr )
!   where
!    ne is the number of elements (data items from BUFR)
!
! INPUTS
!   Values  dflt  Array of values from BUFR decoder
!   NFreq    int  No. of frequency sets in BUFR (1 or 3)
!   ROdata  dtyp  ROPP data - derived type
!
! OUTPUTS
!   ROdata  dtyp  ROPP data - derived type, updated
!   ierr     int  Return code: 0=OK
!                              1=sync error (mis-matched replic. factor)
!                              2=wrong frequency value
!                              3=wrong vertical significance code
!
! USES
!   ropp_io_types  - ROPP file I/O support
!   ropp_io        - i/o routines
!   geodesy        - geodesy routines
!   messages       - messages interface
!
! CALLS
!   ConvertCodes
!   ropp_io_occid
!   geometric2geopotential
!   messages
!   message_get_routine
!   message_set_routine
!
! CALLED BY
!   bufr2ropp
!
! DESCRIPTION
!   Converts decoded BUFR plain array data to ROPP derived type units etc.
!   This procedure is mostly scaling and/or range changing (e.g. Pa to hPa).
!   This routine also performs gross error checking, so that if data is not
!   valid (normally indicated by the BUFR 'missing data' flag value) that data
!   value is left as default "missing" in the ROPP structure.
!
! REFERENCES
!   1) ROPP User Guide - Part I
!      SAF/ROM/METO/UG/ROPP/002
!   2) WMO FM94 (BUFR) Specification for ROM SAF Processed Radio
!      Occultation Data. SAF/ROM/METO/FMT/BUFR/001
!
! AUTHOR
!   Met Office, Exeter, UK.
!   Any comments on this software should be given via the ROM SAF
!   Helpdesk at http://www.romsaf.org
!
! COPYRIGHT
!   (c) EUMETSAT. All rights reserved.
!   For further details please refer to the file COPYRIGHT
!   which you should have received as part of this distribution.
!
!****

! Modules

  USE typesizes,     dp => EightByteReal
  USE ropp_io_types, ONLY: ROprof,          &
                           PCD_occultation
  USE ropp_io,       ONLY: ropp_io_occid
  USE geodesy,       ONLY: geometric2geopotential

  IMPLICIT NONE

! Fixed parameters

  CHARACTER (LEN=*), PARAMETER :: FmtNum = "(I20)"
  REAL(dp), PARAMETER :: MISSING = RVIND - 0.7E38_dp  ! Missing data check value (10^38)

! Argument list parameters

  REAL(dp),      INTENT(IN)    :: Values(:) ! Decoded values
  INTEGER,       INTENT(IN)    :: NFreq     ! No. of frequencies
  TYPE (ROprof), INTENT(INOUT) :: ROdata    ! RO data structure
  INTEGER,       INTENT(OUT)   :: ierr      ! Return code

! Local parameters

  CHARACTER (LEN=50) :: routine  ! saved routine name
  CHARACTER (LEN=80) :: outmsg   ! output message string
  CHARACTER (LEN=20) :: number   ! temporary string for numeric values
  CHARACTER (LEN=16) :: QC_str   ! temporary string to hold QC values
  CHARACTER (LEN=4)  :: Ccode    ! ICAO code associated with Ocode
  INTEGER            :: Gclass   ! GNSS satellite class value
  INTEGER            :: Gcode    ! GNSS PRN
  INTEGER            :: Lcode    ! LEO  code value
  INTEGER            :: Icode    ! Instrument code value
  INTEGER            :: Ocode    ! Origin. centre code value
  INTEGER            :: Scode    ! Sub-centre code value
  INTEGER            :: Bcode    ! B/G generator code value
  INTEGER            :: ProdType ! Product type code
  INTEGER            :: PCD      ! PCD bit flags (16-bit)
  INTEGER            :: in       ! loop counter for profile arrays
  INTEGER            :: IE       ! index to Values element
  INTEGER            :: RepFac   ! Replication Factor
  INTEGER            :: FOSsig   ! First-order statistics significance code
  INTEGER            :: TimeSig  ! Time signficance code
  INTEGER            :: ioerr    ! I/O error code
  REAL               :: SWver    ! Software version number
  REAL(dp)           :: lat      ! Nominal latitude of occultation
  REAL(dp)           :: ht       ! Sample geometric height
  REAL(dp)           :: MeanFreq ! Mean frequency (Hz)

!-------------------------------------------------
! 0. Initialise
!-------------------------------------------------

  CALL message_get_routine ( routine )
  CALL message_set_routine ( 'ConvertBUFRtoROPP' )

  ierr = 0

!-------------------------------------------------
! 1. ROPP Header block
!-------------------------------------------------

  CALL message ( msg_diag, '==========================' )
  CALL message ( msg_diag, 'Contents of BUFR Section 4' )
  CALL message ( msg_diag, '==========================' )

  CALL message ( msg_diag, '------------------------' )
  CALL message ( msg_diag, 'Radio Occultation header' )
  CALL message ( msg_diag, '------------------------' )

  Lcode  = NINT(Values(1))                      ! [001007] LEO ID
  WRITE (number, *) Lcode
  CALL message ( msg_diag, '[001007] LEO ID = ' // TRIM(ADJUSTL(number)) )

  Icode  = NINT(Values(2))                      ! [002019] RO Instrument
  WRITE (number, *) Icode
  CALL message ( msg_diag, '[002019] RO instrument = ' // TRIM(ADJUSTL(number)) )

  Ocode  = NINT(Values(3))                      ! [001033] Proc.centre
  WRITE (number, *) Ocode
  CALL message ( msg_diag, '[001033] Orig. centre = ' // TRIM(ADJUSTL(number)) )

  Scode  = 0
  IF ( Ocode == 78 ) Scode = 173
  WRITE (number, *) Scode
  CALL message ( msg_diag, '[------] Orig. sub-centre = ' // TRIM(ADJUSTL(number)) )

  Bcode  = Ocode
  WRITE (number, *) Bcode
  CALL message ( msg_diag, '[------] Back. centre = ' // TRIM(ADJUSTL(number)) )

  Gclass = NINT(Values(21))                     ! [002020] GNSS class
  WRITE (number, *) Gclass
  CALL message ( msg_diag, '[002020] GNSS class = ' // TRIM(ADJUSTL(number)) )

  Gcode  = NINT(Values(22))                     ! [001050] GNSS PRN
  WRITE (number, *) Gcode
  CALL message ( msg_diag, '[001050] GNSS PRN = ' // TRIM(ADJUSTL(number)) )

  CALL ConvertCodes ( ROdata,        &
                      Gclass, Gcode, &
                      Lcode,  Icode, &
                      Ocode,  Scode, &
                      Ccode,  Bcode, &
                      -1 )

  ProdType = NINT(Values(4))                    ! [002172] Product type
  WRITE (number, *) ProdType
  CALL message ( msg_diag, '[002172] Product type = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(5) < MISSING ) THEN               ! [025060] Software version
    WRITE (number, '(F10.3)') Values(5)
    CALL message ( msg_diag, '[025060] Software version = ' // TRIM(ADJUSTL(number)) )
    SWver = Values(5) * 1E-3
    WRITE ( number,        &
            FMT="(F10.3)", &
            IOSTAT=ioerr ) SWver
    IF ( ioerr == 0 ) THEN
      IF ( SWver < 10.0 ) THEN
        ROdata%Software_Version = "V0" // ADJUSTL ( number )
      ELSE
        ROdata%Software_Version = "V"  // ADJUSTL ( number )
      END IF
    END IF
  ELSE
    WRITE (number, *) '- - - - -'
    CALL message ( msg_diag, '[025060] Software version = ' // TRIM(ADJUSTL(number)) )
  END IF

! Date/time of start of occultation

  CALL message ( msg_diag, '-------------------------' )
  CALL message ( msg_diag, 'Time of occultation start' )
  CALL message ( msg_diag, '-------------------------' )

  TimeSig             = NINT(Values(6))         ! [008021] Time.sig (17=start)
  WRITE (number, *) TimeSig
  CALL message ( msg_diag, '[008021] Time significance (17==>start) = ' // TRIM(ADJUSTL(number)) )

  ROdata%DTocc%Year   = NINT(Values(7))         ! [004001] Year
  WRITE (number, *) ROdata%DTocc%Year
  CALL message ( msg_diag, '[004001] Year = ' // TRIM(ADJUSTL(number)) )

  ROdata%DTocc%Month  = NINT(Values(8))         ! [004002] Month
  WRITE (number, *) ROdata%DTocc%Month
  CALL message ( msg_diag, '[004002] Month = ' // TRIM(ADJUSTL(number)) )

  ROdata%DTocc%Day    = NINT(Values(9))         ! [004003] Day
  WRITE (number, *) ROdata%DTocc%Day
  CALL message ( msg_diag, '[004003] Day = ' // TRIM(ADJUSTL(number)) )

  ROdata%DTocc%Hour   = NINT(Values(10))        ! [004004] Hour
  WRITE (number, *) ROdata%DTocc%Hour
  CALL message ( msg_diag, '[004004] Hour = ' // TRIM(ADJUSTL(number)) )

  ROdata%DTocc%Minute = NINT(Values(11))        ! [004005] Minute
  WRITE (number, *) ROdata%DTocc%Minute
  CALL message ( msg_diag, '[004005] Minute = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(12) < MISSING ) THEN              ! [004006] Seconds & MSecs
    ROdata%DTocc%Second = INT(Values(12))
    ROdata%DTocc%MSec   = NINT(MOD(Values(12),1.0_dp)*1E3)
    WRITE (number, '(F10.3)') Values(12)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[004006] Second = ' // TRIM(ADJUSTL(number)) )

! Summary quality information. Only use 1st 16 bits in unswapped bit order

  CALL message ( msg_diag, '------------------------------' )
  CALL message ( msg_diag, 'RO summary quality information' )
  CALL message ( msg_diag, '------------------------------' )

  IF ( Values(13) < MISSING ) THEN              ! [033039] Quality flags for RO data
    PCD = NINT(Values(13))
    WRITE (number, *) PCD
    ROdata%PCD = 0
    QC_str = ''
    DO in = 0, 15
      IF ( BTEST(PCD, in) ) THEN
        ROdata%PCD = IBSET(ROdata%PCD, 15-in)
        QC_str = '1' // QC_str
      ELSE
        QC_str = '0' // QC_str
      END IF
    END DO
  ELSE
    WRITE (number, *) '- - - - -'
    WRITE (QC_str, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[033039] Quality flags for RO data = ' // TRIM(ADJUSTL(number)) // &
                           ' (base 10) = ' // QC_str // ' (base 2)' )

  IF ( Values(14) < MISSING ) THEN              ! [033007] Percent confidence
    ROdata%Overall_Qual = Values(14)
    WRITE (number, '(F10.3)') Values(14)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[033007] Percent confidence = ' // TRIM(ADJUSTL(number)) )

! We now have enough header info to generate the occultation ID

  CALL ropp_io_occid ( ROdata )

! Background information

  IF ( BTEST(ROdata%PCD,PCD_occultation) ) THEN
    ROdata%bg%Year   = ROdata%DTocc%Year
    ROdata%bg%Month  = ROdata%DTocc%Month
    ROdata%bg%Day    = ROdata%DTocc%Day
    ROdata%bg%Hour   = ROdata%DTocc%Hour
    ROdata%bg%Minute = ROdata%DTocc%Minute
  ELSE
    ROdata%bg%Source = "NONE"
  END IF

! Local Earth parameters and reference POD

  CALL message ( msg_diag, '----------------' )
  CALL message ( msg_diag, 'LEO and GNSS POD' )
  CALL message ( msg_diag, '----------------' )

  IE = 14
  DO in=1,3
    IF ( Values(IE+in) < MISSING ) THEN         ! [027031, 028031, 010031] = LEO (X, Y, Z) pos (m)
      ROdata%GeoRef%LEO_Pod%Pos(in) = Values(IE+in)
      WRITE (number, '(F20.3)') Values(IE+in)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    SELECT CASE (in)
      CASE (1)
        CALL message ( msg_diag, '[027031] LEO X pos (m) = ' // TRIM(ADJUSTL(number)) )
      CASE (2)
        CALL message ( msg_diag, '[028031] LEO Y pos (m) = ' // TRIM(ADJUSTL(number)) )
      CASE (3)
        CALL message ( msg_diag, '[010031] LEO Z pos (m) = ' // TRIM(ADJUSTL(number)) )
    END SELECT
  END DO

  IE = 17
  DO in=1,3
    IF ( Values(IE+in) < MISSING ) THEN         ! [001041, 001042, 001043] = LEO (X, Y, Z) vel (m/s)
      ROdata%GeoRef%LEO_Pod%Vel(in) = Values(IE+in)
      WRITE (number, '(F20.5)') Values(IE+in)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    SELECT CASE (in)
      CASE (1)
        CALL message ( msg_diag, '[001041] LEO X vel (m/s) = ' // TRIM(ADJUSTL(number)) )
      CASE (2)
        CALL message ( msg_diag, '[001042] LEO Y vel (m/s) = ' // TRIM(ADJUSTL(number)) )
      CASE (3)
        CALL message ( msg_diag, '[001043] LEO Z vel (m/s) = ' // TRIM(ADJUSTL(number)) )
    END SELECT
  END DO

  IE = 22
  DO in=1,3
    IF ( Values(IE+in) < MISSING ) THEN         ! [027031, 028031, 010031] = GNSS (X, Y, Z) pos (m)
      ROdata%GeoRef%GNS_Pod%Pos(in) = Values(IE+in)
      WRITE (number, '(F20.3)') Values(IE+in)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    SELECT CASE (in)
      CASE (1)
        CALL message ( msg_diag, '[027031] GNSS X pos (m) = ' // TRIM(ADJUSTL(number)) )
      CASE (2)
        CALL message ( msg_diag, '[028031] GNSS Y pos (m) = ' // TRIM(ADJUSTL(number)) )
      CASE (3)
        CALL message ( msg_diag, '[010031] GNSS Z pos (m) = ' // TRIM(ADJUSTL(number)) )
    END SELECT
  END DO

  IE = 25
  DO in=1,3
    IF ( Values(IE+in) < MISSING ) THEN         ! [001041, 001042, 001043] = GNSS (X, Y, Z) vel (m/s)
      ROdata%GeoRef%GNS_Pod%Vel(in) = Values(IE+in)
      WRITE (number, '(F20.5)') Values(IE+in)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    SELECT CASE (in)
      CASE (1)
        CALL message ( msg_diag, '[001041] GNSS X vel (m/s) = ' // TRIM(ADJUSTL(number)) )
      CASE (2)
        CALL message ( msg_diag, '[001042] GNSS Y vel (m/s) = ' // TRIM(ADJUSTL(number)) )
      CASE (3)
        CALL message ( msg_diag, '[001043] GNSS Z vel (m/s) = ' // TRIM(ADJUSTL(number)) )
    END SELECT
  END DO

  CALL message ( msg_diag, '----------------------' )
  CALL message ( msg_diag, 'Local Earth parameters' )
  CALL message ( msg_diag, '----------------------' )

  IF ( Values(29) < MISSING ) THEN              ! [004016] Time/start (s)
    ROdata%GeoRef%Time_Offset = Values(29)
    WRITE (number, '(F10.3)') Values(29)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[004016] Time increment since start (s) = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(30) < MISSING ) THEN              ! [005001] Latitude (deg)
    ROdata%GeoRef%Lat = Values(30)
    WRITE (number, '(F15.5)') Values(30)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[005001] Latitude (deg) = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(31) < MISSING ) THEN              ! [006001] Longitude (deg)
    ROdata%GeoRef%Lon = Values(31)
    WRITE (number, '(F15.5)') Values(31)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[006001] Longitude (deg) = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(32) < MISSING ) THEN              ! [027031] CofC X (m)
    ROdata%GeoRef%r_CoC(1) = Values(32)
    WRITE (number, '(F10.3)') Values(32)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[027031] Centre of Curvature X (m) = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(33) < MISSING ) THEN              ! [028031] CofC Y (m)
    ROdata%GeoRef%r_CoC(2) = Values(33)
    WRITE (number, '(F10.3)') Values(33)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[028031] Centre of Curvature Y (m) = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(34) < MISSING ) THEN              ! [010031] CofC Z (m)
    ROdata%GeoRef%r_CoC(3) = Values(34)
    WRITE (number, '(F10.3)') Values(34)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[010031] Centre of Curvature Z (m) = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(35) < MISSING ) THEN              ! [010035] Radius value (m)
    ROdata%GeoRef%RoC = Values(35)
    WRITE (number, '(F15.3)') Values(35)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[010035] Radius of Curvature (m) = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(36) < MISSING ) THEN              ! [005021] Line of sight bearing (degT)
    ROdata%GeoRef%Azimuth = Values(36)
    WRITE (number, '(F10.3)') Values(36)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[005021] GNSS->LEO azimuth (degT) = ' // TRIM(ADJUSTL(number)) )

  IF ( Values(37) < MISSING ) THEN              ! [010036] Geoid undulation (m)
    ROdata%GeoRef%Undulation = Values(37)
    WRITE (number, '(F10.3)') Values(37)
  ELSE
    WRITE (number, *) '- - - - -'
  END IF
  CALL message ( msg_diag, '[010036] Geoid undulation (m) = ' // TRIM(ADJUSTL(number)) )

  IE = 37

!-------------------------------------------------
! 2. Level 1b data (bending angle profiles)
!-------------------------------------------------

  CALL message ( msg_diag, "-----------------" )
  CALL message ( msg_diag, "RO 'Step 1b' data" )
  CALL message ( msg_diag, "-----------------" )

  RepFac = NINT(Values(IE+1))                   ! [031002] Sample Replication factor
  WRITE (number, *) RepFac
  CALL message ( msg_diag, '[031002] Replication factor (Lev 1b) = ' // TRIM(ADJUSTL(number)) )

  IF ( RepFac /= ROdata%Lev1b%Npoints ) THEN
    WRITE ( outmsg, FMT="(A,2I10)" ) "Sync error: L1b RepFac /= NPoints", &
                                     repfac, ROdata%Lev1b%Npoints
    CALL message ( msg_error, TRIM(outmsg) )
    ierr = 1
  END IF

  DO in = 1, ROdata%Lev1b%Npoints

! Coordinates

    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) THEN
      WRITE (number, *) in
      CALL message ( msg_diag, 'For point i = ' // TRIM(ADJUSTL(number)) // ':' )
    END IF

    IF ( Values(IE+2) < MISSING ) THEN          ! [005001] Latitude (deg)
      ROdata%Lev1b%Lat_tp(in) = Values(IE+2)
      WRITE (number, '(F15.5)') Values(IE+2)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '  [005001] Latitude (deg) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+3) < MISSING ) THEN          ! [006001] Longitude (deg)
      ROdata%Lev1b%Lon_tp(in) = Values(IE+3)
      WRITE (number, '(F15.5)') Values(IE+3)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '  [006001] Longitude (deg) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+4) < MISSING ) THEN          ! [005021] Line of sight bearing (degT)
      ROdata%Lev1b%Azimuth_tp(in) = Values(IE+4)
      WRITE (number, '(F10.3)') Values(IE+4)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '  [005021] GNSS->LEO azimuth (degT) = ' // TRIM(ADJUSTL(number)) )

! Extract L1+L2 if they were encoded

    RepFac = NINT(Values(IE+5))                 ! [031001] Frequency Replication factor
    WRITE (number, *) RepFac
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '  [031001] Frequency replication factor = ' // TRIM(ADJUSTL(number)) )

    IF ( RepFac /= nFreq ) THEN
      WRITE ( outmsg, FMT="(A,I12,3I6)" ) "Sync error: L1b RepFac /= nFreq ", &
                                           repfac, nFreq, in, ie+5
      CALL message ( msg_error, TRIM(outmsg) )
      ierr = 1
    END IF

    IF ( NFreq == 3 ) THEN

! L1 data

      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) THEN
        CALL message ( msg_diag, '  L1 data:' )
      END IF

      MeanFreq = Values(IE+6)                   ! [002121] Mean frequency (L1=1.5GHz)
      WRITE (number, '(ES15.5)') MeanFreq
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [002121] Mean frequency = ' // TRIM(ADJUSTL(number)) )

      IF ( ABS(MeanFreq - 1.5E9) > 2.0E8 ) THEN
        WRITE ( outmsg, FMT="(A,F5.1,A)") "Wrong L1 frequency (expected around 1.57 GHz): ", &
                                           MeanFreq/1E9," GHz"
        CALL message ( msg_error, TRIM(outmsg) )
        ierr = 2
      END IF

      IF ( Values(IE+7) < MISSING ) THEN        ! [007040] Impact parameter (m)
        ROdata%Lev1b%Impact_L1(in) = Values(IE+7)
        WRITE (number, '(F15.3)') Values(IE+7)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [007040] L1 Impact parameter (m) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+8) < MISSING ) THEN        ! [015037] B/angle (rad)
        ROdata%Lev1b%BAngle_L1(in) = Values(IE+8)
        WRITE (number, '(F20.8)') Values(IE+8)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [015037] L1 Bending angle (rad) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+9) < MISSING ) THEN        ! [008023] First-order stats signif. (missing==>off)
        FOSsig = NINT(Values(IE+9))
        WRITE (number, *) FOSsig
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [008023] L1 First-order stats signif. (13==>RMS) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+10) < MISSING ) THEN       ! [015037] B/angle error (rad)
        ROdata%Lev1b%BAngle_L1_Sigma(in) = Values(IE+10)
        WRITE (number, '(F20.8)') Values(IE+10)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [015037] L1 Bending angle error (rad) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+11) < MISSING ) THEN       ! [008023] First-order stats signif. (missing==>off)
        FOSsig = NINT(Values(IE+11))
        WRITE (number, *) FOSsig
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [008023] L1 First-order stats signif. (missing==>off) = ' // TRIM(ADJUSTL(number)) )

! L2 data

      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) THEN
        CALL message ( msg_diag, '  L2 data:' )
      END IF

      MeanFreq = Values(IE+12)                  ! [002121] Mean frequency (L2=1.2GHz)
      WRITE (number, '(ES15.5)') MeanFreq
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [002121] Mean frequency = ' // TRIM(ADJUSTL(number)) )

      IF ( ABS(MeanFreq - 1.2E9) > 2.0E8 ) THEN
        WRITE ( outmsg, FMT="(A,F5.1,A)") "Wrong L2/L5 frequency (expected around 1.2 GHz): ", &
                                           MeanFreq/1E9," GHz"
        CALL message ( msg_error, TRIM(outmsg) )
        ierr = 2
      END IF

      IF ( Values(IE+13) < MISSING ) THEN       ! [007040] Impact parameter (m)
        ROdata%Lev1b%Impact_L2(in) = Values(IE+13)
        WRITE (number, '(F15.3)') Values(IE+13)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [007040] L2 Impact parameter (m) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+14) < MISSING ) THEN       ! [015037] B/angle (rad)
        ROdata%Lev1b%BAngle_L2(in) = Values(IE+14)
        WRITE (number, '(F20.8)') Values(IE+14)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [015037] L2 Bending angle (rad) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+15) < MISSING ) THEN       ! [008023] First-order stats signif. (missing==>off)
        FOSsig = NINT(Values(IE+15))
        WRITE (number, *) FOSsig
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [008023] L2 First-order stats signif. (13==>RMS) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+16) < MISSING ) THEN       ! [015037] B/angle error (rad)
        ROdata%Lev1b%BAngle_L2_Sigma(in) = Values(IE+16)
        WRITE (number, '(F20.8)') Values(IE+16)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [015037] L2 Bending angle error (rad) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+17) < MISSING ) THEN       ! [008023] First-order stats signif. (missing==>off)
        FOSsig = NINT(Values(IE+17))
        WRITE (number, *) FOSsig
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
        CALL message ( msg_diag, '    [008023] L2 First-order stats signif. (missing==>off) = ' // TRIM(ADJUSTL(number)) )

    ELSE

      IE = IE - 12

    END IF

! Corrected bending angle (always encoded)

    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) THEN
      CALL message ( msg_diag, '  LC data:' )
    END IF

    MeanFreq = Values(IE+18)                    ! [002121] Mean frequency (Corrected=0)
    WRITE (number, '(ES15.5)') MeanFreq
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '    [002121] Mean frequency = ' // TRIM(ADJUSTL(number)) )

    IF ( NINT(MeanFreq) /= 0 ) THEN
      WRITE ( outmsg, FMT="(A,F5.1,A)" ) "Wrong Corr (0.0) frequency: ", &
                                         MeanFreq/1E9," GHz"
      CALL message ( msg_error, TRIM(outmsg) )
      ierr = 2
    END IF

    IF ( Values(IE+19) < MISSING ) THEN         ! [007040] Impact parameter (m)
      ROdata%Lev1b%Impact(in) = Values(IE+19)
      WRITE (number, '(F15.3)') Values(IE+19)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '    [007040] Impact parameter (m) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+20) < MISSING ) THEN         ! [015037] B/angle (rad)
      ROdata%Lev1b%BAngle(in) = Values(IE+20)
      WRITE (number, '(F20.8)') Values(IE+20)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '    [015037] Bending angle (rad) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+21) < MISSING ) THEN         ! [008023] First-order stats signif. (missing==>off)
      FOSsig = NINT(Values(IE+21))
      WRITE (number, *) FOSsig
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '    [008023] First-order stats signif. (13==>RMS) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+22) < MISSING ) THEN         ! [015037] B/angle error (rad)
      ROdata%Lev1b%BAngle_Sigma(in) = Values(IE+22)
      WRITE (number, '(F20.8)') Values(IE+22)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '    [015037] Bending angle error (rad) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+23) < MISSING ) THEN         ! [008023] First-order stats signif. (missing==>off)
      FOSsig = NINT(Values(IE+23))
      WRITE (number, *) FOSsig
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '    [008023] First-order stats signif. (missing==>off) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+24) < MISSING ) THEN         ! [033007] Percent confidence
      ROdata%Lev1b%Bangle_Qual(in) = Values(IE+24)
      WRITE (number, '(F10.3)') Values(IE+24)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev1b%Npoints) ) &
      CALL message ( msg_diag, '    [033007] Percent confidence = ' // TRIM(ADJUSTL(number)) )

    IE = IE + 23  ! L1b values per sample

  END DO

  IE = IE + 1     ! L1b RepFac

!-------------------------------------------------
! 3. Level 2a data (derived refractivity profile)
!-------------------------------------------------

  CALL message ( msg_diag, "-----------------" )
  CALL message ( msg_diag, "RO 'Step 2a' data" )
  CALL message ( msg_diag, "-----------------" )

  lat = 0.0
  IF ( ROdata%GEOref%lat >= -90.0_dp ) &
    lat = ROdata%GEOref%lat

  RepFac = NINT(Values(IE+1))                   ! [031002] Sample Replication factor
  WRITE (number, *) RepFac
  CALL message ( msg_diag, '[031002] Replication factor (Lev 2a) = ' // TRIM(ADJUSTL(number)) )

  IF ( RepFac /= ROdata%Lev2a%Npoints ) THEN
    WRITE ( outmsg, FMT="(A,2I10)" ) "Sync error: L2a RepFac /= NPoints", &
                                      repfac, ROdata%Lev2a%Npoints
    CALL message ( msg_error, TRIM(outmsg) )
    ierr = 1
  END IF

  DO in = 1, ROdata%Lev2a%Npoints

    IF ( (in == 1) .OR. (in == ROdata%Lev2a%Npoints) ) THEN
      WRITE (number, *) in
      CALL message ( msg_diag, 'For point i = ' // TRIM(ADJUSTL(number)) // ':' )
    END IF

    IF ( Values(IE+2) < MISSING ) THEN          ! [007007] Height amsl (m)
      ROdata%Lev2a%Alt_Refrac(in) = Values(IE+2)
! Geopot ht is not in BUFR: re-create a value by applying
! latitude- & height-dependent gravity model
      ht = ROdata%Lev2a%Alt_Refrac(in)
      ROdata%Lev2a%Geop_Refrac(in) = geometric2geopotential(lat,ht) ! Geopot ht (gpm)
      WRITE (number, '(F10.3)') Values(IE+2)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2a%Npoints) ) &
      CALL message ( msg_diag, '  [007007] Height AMSL (m) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+3) < MISSING ) THEN          ! [015036] Refractivity (N-units)
      ROdata%Lev2a%Refrac(in) = Values(IE+3)
      WRITE (number, '(F10.3)') Values(IE+3)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2a%Npoints) ) &
      CALL message ( msg_diag, '  [015036] Refractivity (N-units) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+4) < MISSING ) THEN          ! [008023] First-order stats signif. (missing==>off)
      FOSsig = NINT(Values(IE+4))
      WRITE (number, *) FOSsig
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2a%Npoints) ) &
      CALL message ( msg_diag, '  [008023] First-order stats signif. (13==>RMS) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+5) < MISSING ) THEN          ! [015036] Refractivity error (N-units)
      ROdata%Lev2a%Refrac_Sigma(in) = Values(IE+5)
      WRITE (number, '(F10.3)') Values(IE+5)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2a%Npoints) ) &
      CALL message ( msg_diag, '  [015036] Refractivity error (N-units) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+6) < MISSING ) THEN          ! [008023] First-order stats signif. (missing==>off)
      FOSsig = NINT(Values(IE+6))
      WRITE (number, *) FOSsig
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2a%Npoints) ) &
      CALL message ( msg_diag, '  [008023] First-order stats signif. (missing==>off) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+7) < MISSING ) THEN          ! [033007] Percent confidence
      ROdata%Lev2a%Refrac_Qual(in) = Values(IE+7)
      WRITE (number, '(F10.3)') Values(IE+7)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2a%Npoints) ) &
      CALL message ( msg_diag, '  [033007] Percent confidence = ' // TRIM(ADJUSTL(number)) )

    IE = IE + 6   ! L2a values per sample

  END DO

  IE = IE + 1     ! L2a RepFac

!-------------------------------------------------
! 4. Level 2b data (retrieved P,T,q profile)
!-------------------------------------------------

  CALL message ( msg_diag, "-----------------" )
  CALL message ( msg_diag, "RO 'Level2b' data" )
  CALL message ( msg_diag, "-----------------" )

  RepFac = NINT(Values(IE+1))                   ! [031002] Sample Replication factor
  WRITE (number, *) RepFac
  CALL message ( msg_diag, '[031002] Replication factor (Lev 2b) = ' // TRIM(ADJUSTL(number)) )

  IF ( RepFac /= ROdata%Lev2b%Npoints ) THEN
    WRITE ( outmsg, FMT="(A,2I10)" ) "Sync error: L2b RepFac /= NPoints", &
                                     repfac, ROdata%Lev2b%Npoints
    CALL message ( msg_error, TRIM(outmsg) )
    ierr = 1
  END IF

  DO in = 1, ROdata%Lev2b%Npoints

    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) THEN
      WRITE (number, *) in
      CALL message ( msg_diag, 'For point i = ' // TRIM(ADJUSTL(number)) // ':' )
    END IF

    IF ( Values(IE+2) < MISSING ) THEN             ! [007009] Geopot ht (gpm)
      ROdata%Lev2b%Geop(in) = Values(IE+2)
      WRITE (number, '(F10.3)') Values(IE+2)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [007009] Geopotential height (gpm) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+3) < MISSING ) THEN          ! [010004] Pressure (Pa-->hPa)
      ROdata%Lev2b%Press(in) = Values(IE+3) * 1E-2
      WRITE (number, '(F10.3)') Values(IE+3)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [010004] Pressure (Pa) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+4) < MISSING ) THEN          ! [012001] Temperature (K)
      ROdata%Lev2b%Temp(in) = Values(IE+4)
      WRITE (number, '(F10.3)') Values(IE+4)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [012001] Temperature (K) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+5) < MISSING ) THEN          ! [013001] Spec/humidity (Kg/Kg)
      ROdata%Lev2b%SHum(in) = Values(IE+5) * 1E3
      WRITE (number, '(F10.5)') Values(IE+5)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [013001] Specific humidity (Kg/Kg) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+6) < MISSING ) THEN          ! [008023] First-order stats signif. (missing==>off)
      FOSsig = NINT(Values(IE+6))
      WRITE (number, *) FOSsig
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [008023] First-order stats signif. (13==>RMS) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+7) < MISSING ) THEN          ! [010004] Pressure error (Pa-->hPa)
      ROdata%Lev2b%Press_Sigma(in) = Values(IE+7) * 1E-2
      WRITE (number, '(F10.3)') Values(IE+7)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [010004] Pressure error (Pa) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+8) < MISSING ) THEN          ! [012001] Temperature error (K)
      ROdata%Lev2b%Temp_Sigma(in) = Values(IE+8)
      WRITE (number, '(F10.3)') Values(IE+8)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [012001] Temperature error (K) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+9) < MISSING ) THEN          ! [013001] S/Hum error (Kg/Kg)
      ROdata%Lev2b%SHum_Sigma(in) = Values(IE+9) * 1E3
      WRITE (number, '(F10.5)') Values(IE+9)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
       CALL message ( msg_diag, '  [013001] Spec/humidity error (Kg/Kg) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+10) < MISSING ) THEN         ! [008023] First-order stats signif. (missing==>off)
      FOSsig = NINT(Values(IE+10))
      WRITE (number, *) FOSsig
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [008023] First-order stats signif. (missing==>off) = ' // TRIM(ADJUSTL(number)) )

    IF ( Values(IE+11) < MISSING ) THEN        ! [033007] Percent confidence
      ROdata%Lev2b%Meteo_Qual(in) = Values(IE+11)
      WRITE (number, '(F10.3)') Values(IE+11)
    ELSE
      WRITE (number, *) '- - - - -'
    END IF
    IF ( (in == 1) .OR. (in == ROdata%Lev2b%Npoints) ) &
      CALL message ( msg_diag, '  [033007] Percent confidence = ' // TRIM(ADJUSTL(number)) )

    IE = IE + 10   ! L2b values per sample

  END DO

  IE = IE + 1     ! L2b RepFac

!-------------------------------------------------
! 5. Level 2c data (retrieved surface params)
!-------------------------------------------------

  CALL message ( msg_diag, "-----------------" )
  CALL message ( msg_diag, "RO 'Level2c' data" )
  CALL message ( msg_diag, "-----------------" )

  WRITE (number, *) ROdata%Lev2c%Npoints
  CALL message ( msg_diag, '[------] Replication factor (Lev 2c) = ' // TRIM(ADJUSTL(number)) )

  IF ( NINT(Values(IE+1)) == 0 ) THEN           ! [008003] Vertical signif. (0=surface)

    IF ( ROdata%Lev2c%Npoints == 1 ) THEN

      WRITE (number, *) NINT(Values(IE+1)) ! Must be zero
      CALL message ( msg_diag, '[008003] Vertical significance (0==>surface) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+2) < MISSING ) THEN        ! [007009] Geopot.Ht. (of surf) (gpm)
        ROdata%Lev2c%Geop_Sfc = Values(IE+2)
        WRITE (number, '(F10.3)') Values(IE+2)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      CALL message ( msg_diag, '[007009] Geopotential height (of surface) (gpm) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+3) < MISSING ) THEN        ! [010004] Surface pressure (Pa-->hPa)
        ROdata%Lev2c%Press_Sfc = Values(IE+3) * 1E-2
        WRITE (number, '(F10.3)') Values(IE+3)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      CALL message ( msg_diag, '[010004] Surface pressure (Pa) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+4) < MISSING ) THEN        ! [008023] First-order stats signif. (missing==>off)
        FOSsig = NINT(Values(IE+4))
        WRITE (number, *) FOSsig
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      CALL message ( msg_diag, '[008023] First-order stats signif. (13==>RMS) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+5) < MISSING ) THEN        ! [010004] S/press error (Pa-->hPa)
        ROdata%Lev2c%Press_Sfc_Sigma = Values(IE+5) * 1E-2
        WRITE (number, '(F10.3)') Values(IE+5)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      CALL message ( msg_diag, '[010004] Surface pressure error (Pa) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+6) < MISSING ) THEN        ! [008023] First-order stats signif. (missing==>off)
        FOSsig = NINT(Values(IE+6))
        WRITE (number, *) FOSsig
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      CALL message ( msg_diag, '[008023] First-order stats signif. (missing==>off) = ' // TRIM(ADJUSTL(number)) )

      IF ( Values(IE+7) < MISSING ) THEN        ! [033007] Percent confidence
        ROdata%Lev2c%Press_Sfc_Qual = Values(IE+7)
        WRITE (number, '(F10.3)') Values(IE+7)
      ELSE
        WRITE (number, *) '- - - - -'
      END IF
      CALL message ( msg_diag, '[033007] Percent confidence = ' // TRIM(ADJUSTL(number)) )

    END IF

  ELSE IF ( Values(IE+1) < MISSING ) THEN

    WRITE ( number, FMT=FmtNum ) NINT(Values(IE+1))
    CALL message ( msg_error, "Wrong vertical signficance, not surface (=0): "// &
                               TRIM(ADJUSTL(number)) )
    ierr = 3

  END IF

  CALL message_set_routine ( routine )

END SUBROUTINE ConvertBUFRtoROPP
!------------------------------------------------------------------------

END MODULE bufr2ropp