! $Id: eum2bufr_ec.f90 2738 2011-03-22 15:34:18Z frdo $

PROGRAM eum2bufr_eccodes

!****x* Programs/eum2bufr_eccodes *
!
! NAME
!   eum2bufr    (eum2bufr_eccodes.f90)
!
! SYNOPSIS
!   Encode a EUMETSAT file to WMO FM-94 (BUFR) [ECMWF ecCodes library]
!
!   > export BUFR_TABLES=bufr_directory_path
!   > eum2bufr_eccodes eum_file [eum_file ...] [-o bufr_file]
!                         [-g[i]] [-s csn_file]
!                         [-p thin_file] [-t time] [-r resn]
!                         [-e] [-u] [-l] [-m] [-h|?] [-v] [-d]
! INPUTS
!   eum_file is the input file(s) which must be in EUM netCDF4 format
!
! OUTPUTS
!   bufr_file is the output file, which will contain one encoded
!             BUFR message per input profile (See Ref.2)
!             The output file name is optional, and if not specified,
!             is generated from the occulation ID.
!
! OPTIONS
!   Option switches can be in any order and are case-insensitive;
!   any space(s) between a switch and its (mandatory) argument is
!   optional.
!     -o  specifies the BUFR output file name
!     -s  specifies a channel sequence number file
!     -p  specifies a thinning control file or max. no. of levels
!     -t  specifies a time (age) rejection threshold
!     -u  leave profiles unordered - disables the default re-ordering
!         of output profiles to ascending.
!         NB: using -u, profiles thinned using one of the interpolation methods
!         will retain the order of the fixed levels in the control file; other
!         methods will retain the ordering of the input profiles.
!     -l  do not encode L1+L2 data (bending angle, Level 1b), if present
!     -m  do not encode met.  data (geophysical, Level 2b,c), if present
!     -e  take L2 extrapolated bending angle from netCDF4 EUM file
!     -r  resolution group of netCDF4 EUM file
!     -d  to output additional diagnostics
!     -h or ? to output summary help
!     -v  to output program version ID
!   Defaults
!     Input file name         : none - at least one required
!     Output file name        : <occid>.bufr
!     GTS routing headers     : not generated
!     Channel sequence        : initialised at 001
!     Time threshold          : 00:00 (no cut-off) unless one of
!                               -g options present, when 23:50
!     Encode                  : all available Level 1b, 2b & 2c data
!     Thinning                : sample to no more than 375 levels
!     Resolution              : 'thinned'
!     L2 bending angle        : not extrapolated
!
! CALLS
!   IARGC
!   ConvertROPPtoBUFR
!   EncodeBUFR
!   GetOptions
!   ropp_io_ascend
!   ropp_io_occid
!   ropp_io_read
!   ropp_io_thin
!   ropp_io_free
!   PBOPEN
!   PBWRITE
!   PBCLOSE
!   GTShdrCSN
!   GTShdrIPH
!   CalToJul
!   Date_and_Time_UTC
!   DateTimeOffset
!   To_Lower
!   File_Delete
!   message
!   message_set_routine
!
! MODULES
!   ropp2bufr           - Fixed parameter definitions & common routines
!   ropp_io             - ROPP I/O file support
!   ropp_io_types       - ROPP data type definitions
!   DateTimeProgs       - Date & Time conversion routines
!   DateTimeTypes       - Date & Time conversion definitions
!   messages            - ROPP message library
!   system              - System utility routines
!   GTShdrs             - Routines to add WMO/GTS routing header/trailer
!
! DEPENDENCIES
!   ECMWF BUFR package  - BUFR kernel routines
!   ROPP I/O library    - ROPP file I/O support
!   ROPP Utils library  - ROPP utility routines
!   netCDF library      - netCDF file support
!
! ENVIRONMENT VARIABLES
!   BUFR_TABLES         - Path for run-time files
!
! ERRORS
!   Program (shell) return codes:
!     0 = OK
!     1 = At least one Warning occurred
!     2 = At least one Error occurred
!     3 = A Fatal error occurred
!
! DESCRIPTION
!   A BUFR encoder for Radio Occultation data.
!   Reads from one or more EUM netCDF 4 files and encodes data therein
!   to one BUFR message per profile output to a single BUR file.
!   Various options are provided to control the generation of GTS
!   routing headers and rejection based on the age of the data
!   and to skip encoding certain profile subsets and thinning (see Refs.3,4).
!   BUFR tables and other run-time files are found via the environment
!   variable 'BUFR_TABLES'.
!
! 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
!   3) Monodimensional data thinning for GPS radio occultations
!      SAF/GRAS/METO/ALG/ROPP/001
!   4) ROPP thinner algorithm
!      SAF/GRAS/METO/REP/GSR/8
!
! SEE ALSO
!   ropp2bufr(1), 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.
!
!****

! Modules
  USE eccodes
  USE messages
  USE ropp2bufr
  USE ropp_io_types, ONLY: ROprof
  USE ropp_io,       ONLY: ropp_io_nrec,  &
                           ropp_io_occid, &
                           ropp_io_read,  &
                           ropp_io_thin,  &
                           ropp_io_free
  USE ropp_utils,    ONLY: File_Delete
  USE DateTimeProgs, ONLY: Date_and_Time_UTC, &
                           DateTimeOffset,    &
                           MonthOfYear,       &
                           TimeSince
  USE DateTimeTypes, ONLY: IdxYear, IdxMonth, IdxDay, &
                           IdxHour, IdxMinute
  USE GTShdrs,       ONLY: GTShdrCSN, &
                           GTShdrIPH, &
                           GTShdrGAD
  USE system,        ONLY: setenv

  IMPLICIT NONE

! Fixed values

  CHARACTER (LEN=*), PARAMETER :: DTfmt1 = & ! hh:mm dd-mm-yyyy
                     "(I2.2,':',I2.2,'UT ',I2.2,'-',A3,'-',I4.4)"
  CHARACTER (LEN=3), PARAMETER :: Month(0:12) = &
                     (/ "???", "Jan", "Feb", "Mar", "Apr", "May", "Jun", &
                               "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" /)

! Local variables

  CHARACTER (LEN=256),  &
          DIMENSION(:), &
          ALLOCATABLE :: EUMdsn         ! Name(s) of input EUM file(s)
  CHARACTER (LEN=256) :: BUFRdsn        ! Output (BUFR) file name
  CHARACTER (LEN=256) :: CSNdsn         ! Channel sequence file name
  CHARACTER (LEN=256) :: Thindsn        ! Thinning control file name
  CHARACTER (LEN=4)   :: OrigICAO       ! Originating centre ICAO code
  INTEGER             :: BUFRunit       ! BUFR I/O unit (C FILE pointer)
  INTEGER             :: nFiles         ! No. of file names on command line
  INTEGER             :: OrigCentre     ! Originating centre BUFR code
  INTEGER             :: SubCentre      ! Originating sub-centre BUFR code
  INTEGER             :: GTShdrType     ! Code for GTS header generation
  INTEGER             :: RejTimeDiff    ! reject obs older than this
  LOGICAL             :: CorrOnly       ! Flag for Corrected only
  LOGICAL             :: nomet          ! Flag for no met. data
  LOGICAL             :: unordered      ! Enable ordering of profiles to ascending
  CHARACTER (LEN=20)  :: centre = 'EUM' ! Centre (EUM here)
  LOGICAL             :: getextrap      ! Take L2 extrapolated bending angle from input files
  CHARACTER (LEN=20)  :: resolution     ! Resolution group of EUM netCDF4 files

  INTEGER,  ALLOCATABLE :: Descr(:)     ! Descriptor sequence
  INTEGER,  ALLOCATABLE :: RepFac(:)    ! Replication Factors
  REAL(dp), ALLOCATABLE :: Values(:)    ! Data values for BUFR (d.p.)

  CHARACTER (LEN=10) :: cIPH            ! Dummy IPH end-of-file sequence
  INTEGER :: LenBUF = 0                 ! Length of BUFR message/bulletin
  INTEGER :: LenEOF                     ! Length of end-of-file message
  INTEGER :: LenIPH                     ! Length of IPH sequence
  INTEGER :: LenTot = 0                 ! Total length of all bulletins (BUFR
                                        ! messages + headers) written out to
                                        ! BUFR file
  INTEGER :: uBYTES(10), pBYTES(3)      ! unpacked & packed header byte containers

  INTEGER :: nElem                      ! No. of elements
  INTEGER :: nValues                    ! No. of elements after thinning
  INTEGER :: nRepFac                    ! No. of Replication Factors
  INTEGER :: CSN = 0                    ! Channel sequence number
  INTEGER :: ierr, status               ! File error & return status codes
  INTEGER :: nmsg = 0                   ! Count of BUFR messages
  INTEGER :: nfreq  = 0                 ! No. of Level 1b frequencies
  INTEGER :: nvalid = 0                 ! No. of valid L1 b/angles
  INTEGER :: i, in, iprof, ifile, n,j,k ! Loop counters
  LOGICAL :: exists                     ! File present flag
  LOGICAL :: first = .TRUE.             ! First profile flag
  LOGICAL :: lgetlevel1b=.TRUE.         ! Always get level 1b data
  CHARACTER(LEN=20)::lgetlevel1a='none' ! For bufr there is no need to get full level 1a data
  LOGICAL :: lgetdirect=.FALSE.         ! Without full level 1a data, this one doesn't matter

  TYPE(ROprof)          :: ROdata       ! ROPP data structure
  INTEGER               :: nprofs       ! No. of profiles in i/p file
  INTEGER               :: tprofs = 0   ! Total profiles in all files

  INTEGER, DIMENSION(8) :: DT8, Offset  ! Date/time arrays
  CHARACTER (LEN=80)    :: outmsg       ! Output text string
  CHARACTER (LEN=10)    :: MonthName    ! Month name
  CHARACTER (LEN=10)    :: number       ! Numeric value as string
  REAL(dp)              :: MinRej       ! Rejection threshold (minutes since 00:00UT 1-Jan-2000)
  REAL(dp)              :: MinObs       ! Ob time (minutes since 00:00UT 1-Jan-2000)

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

  ! New variables for ecCodes
  CHARACTER*80                       :: descriptor ! Holds descriptor string for codes_set e.g. #1#airTemperature
  INTEGER                            :: iret, imessage, ioutfile
  INTEGER, DIMENSION(:), ALLOCATABLE :: ivalues
  REAL, DIMENSION(:), ALLOCATABLE    :: rvalues
  CHARACTER (LEN=256)                :: outfile
  INTEGER                            :: istat ! For setting environment variable
  CHARACTER (LEN=18)                 :: GTShdrTxt ! The input to gtsheaders_bufr.pl
  CHARACTER (LEN=5)                  :: GTShdrOpt ! Either '-g' or '-gi'

!--------------------------------------------------------------
!  1. Initalise
!--------------------------------------------------------------

  CALL message_set_routine ( "eum2bufr_eccodes" )

  CALL message(msg_noin, '')
  CALL message(msg_noin, &
       '---------------------------------------------------------------------')
  CALL message(msg_noin, &
       '                  EUMETSAT to BUFR (ecCodes) Encoder'                 )
  CALL message(msg_noin, &
       '---------------------------------------------------------------------')
  CALL message(msg_noin, '')

  BUFRunit = 0

  ! Set environment variable needed for out-of-bounds BUFR encoding
  istat = setenv('ECCODES_BUFR_SET_TO_MISSING_IF_OUT_OF_RANGE'//ACHAR(0), &
                 '1'//ACHAR(0), 1)

!--------------------------------------------------------------
! 2. Parse command line options
!--------------------------------------------------------------

  nFiles = MAX ( IARGC(), 1 )
  ALLOCATE ( EUMdsn(nFiles), STAT=status )
  IF ( status /= 0 ) &
    CALL message ( msg_fatal, "Failed to allocate memory for ipfile array" )
  CALL GetOptions ( centre,      &
                    EUMdsn,      &
                    nFiles,      &
                    BUFRdsn,     &
                    CSNdsn,      &
                    Thindsn,     &
                    GTShdrType,  &
                    RejTimeDiff, &
                    CorrOnly,    &
                    nomet,       &
                    unordered,   &
                    getextrap,   &
                    resolution   )

!--------------------------------------------------------------
! 3. If time rejection on, set time rejection threshold in
!    minutes since 00:00UT 1-Jan-2000 for specified period
!    back from 'now'.
!--------------------------------------------------------------

  IF ( RejTimeDiff > 0 ) THEN
    Offset = (/0,0,0,0,0,RejTimeDiff,0,0/)
    CALL DateTimeOffset ( DT8, "-", Offset )
    CALL TimeSince      ( DT8, MinRej, 1, Base="JM2000" )
    CALL MonthOfYear    ( DT8(IdxMonth), MonthName, 1 )
    WRITE ( outmsg, FMT=DTfmt1 ) DT8(IdxHour),   &
                                 DT8(IdxMinute), &
                                 DT8(IdxDay),    &
                                 MonthName(1:3), &
                                 DT8(IdxYear)
    CALL message ( msg_diag, " Rejecting occultations older than "// &
                             TRIM(outmsg) )
  ELSE
    MinRej = 0.0_dp
  END IF

!--------------------------------------------------------------
! 4. If GTS headers to be generated, issue a warning and advice
!    on how to do it
!--------------------------------------------------------------

  IF ( GTShdrType /= NOhdrs ) THEN

    GTShdrOpt = "'-g'"
    IF ( GTShdrType == IPhdrs ) GTShdrOpt = "'-gi'"

    CALL message(msg_warn, &
         TRIM(GTShdrOpt) // ' requested, but GTS header encoding not ' // &
         'currently available in ecCodes.')
    CALL message(msg_cont, &
         'For this functionality either use the perl script at ' // &
         '$ROPP_SRC/ropp_io/tools/gtsheaders_bufr.pl with the instructions below,')
    CALL message(msg_cont, &
         'or use the MetDB or ECMWF BUFR packages.')
    CALL message(msg_noin, '')

  ENDIF

!--------------------------------------------------------------
! 5. Loop over input files
!--------------------------------------------------------------

  DO ifile = 1, nFiles

    INQUIRE ( FILE=EUMdsn(ifile), EXIST=exists )
    IF ( .NOT. exists ) THEN
      CALL message ( msg_error, "EUM input file '"// &
                                TRIM(EUMdsn(ifile))//"' not found" )
      CYCLE
    ENDIF

    CALL message ( msg_info, "Reading  EUM data from "// &
                             TRIM(EUMdsn(ifile)) )

    nprofs = ropp_io_nrec ( EUMdsn(ifile) )
    IF ( nprofs < 0 ) nprofs = 1            ! assume 1 profile for text type
    tprofs = tprofs + nprofs

!--------------------------------------------------------------
! 6. Loop over occultations from current file
!    (If a read error, skip to next file)
!--------------------------------------------------------------

    DO iprof = 1, nprofs

      CALL ropp_io_read ( ROdata,                 &
                          file=EUMdsn(ifile),     &
                          rec=iprof,              &
                          ierr=status,            &
                          centre=centre,          &
                          getlevel1b=lgetlevel1b, &
                          getextrap=getextrap,    &
                          resolution=resolution,  &
                          getlevel1a=lgetlevel1a, &
                          getdirect=lgetdirect)

      IF ( status /= 0 ) &
        CALL message ( msg_fatal, "Failed to read file" )

!--------------------------------------------------------------
! 7. On first profile, open output file for BUFR (default name
!    from first occultation ID)
!--------------------------------------------------------------

      CALL ropp_io_occid ( ROdata )
      IF ( first ) THEN
        IF ( BUFRdsn == " " ) THEN
          BUFRdsn = TRIM(ROdata%Occ_id) // ".bufr"
          CALL To_Lower ( BUFRdsn )
        END IF

        CALL codes_open_file( BUFRunit,TRIM(BUFRdsn) ,'w' , ierr)

        IF ( ierr /= 0 ) THEN
          CALL message ( msg_fatal, "Failed to open BUFR output file "// &
                                    TRIM(BUFRdsn) )
        END IF
        first = .FALSE.
      END IF

      WRITE( number, FMT="(I4)") iprof
      CALL message ( msg_info, "Encoding profile "//TRIM(number)// &
                               " : "//TRIM(ROdata%Occ_ID) )

!--------------------------------------------------------------
! 7.1 If GTS time rejection on, skip if occultation time
!     is too old
!--------------------------------------------------------------

      IF ( MinRej > 0.5_dp ) THEN
        DT8 = (/ROdata%DTocc%Year,   ROdata%DTocc%Month,  &
                ROdata%DTocc%Day,    0,                   &
                ROdata%DTocc%Hour,   ROdata%DTocc%Minute, &
                ROdata%DTocc%Second, 0/)
        CALL TimeSince ( DT8, MinObs, 1, Base="JM2000" )
        IF ( MinObs < MinRej ) THEN
          CALL message ( msg_warn, "Occultation is too old for GTS "// &
                                   "- not encoded." )
          CYCLE
        END IF
      END IF

!--------------------------------------------------------------
! 7.2 Use (at most) one Level 1a data for nominal POD
!     Only encode Level 1b L1+L2 data if:
!     a) 'Corrected only' option not taken and
!     b) there is at least one valid L1 bending angle value present
!--------------------------------------------------------------

      ROdata%Lev1a%Npoints = MIN ( 1, ROdata%Lev1a%Npoints )

      IF ( .NOT. CorrOnly ) THEN
        nvalid = 0
        DO in = 1, ROdata%Lev1b%Npoints
          IF ( ROdata%Lev1b%BAngle_L1(in) > 0.0 .AND. &
               ROdata%Lev1b%BAngle_L1(in) < 0.082 ) THEN
            nvalid = nvalid + 1
          END IF
        END DO
        IF ( nvalid < 1 ) CorrOnly = .TRUE.
      END IF

      IF ( ROdata%Lev1b%Npoints > 0 ) THEN
        IF ( CorrOnly ) THEN
          nfreq = 1
        ELSE
          nfreq = 3
        END IF
      END IF

!--------------------------------------------------------------
! 7.3 Only encode 'Met' data if:
!     a) 'No Met' option not taken and
!     b) there is at least one valid temperature value present.
!     Level 2c (Surface) data is always encoded, but show as '0'
!     if not valid.
!     Ignore any level 2d data
!--------------------------------------------------------------

      nvalid = 0
      DO in = 1, ROdata%Lev2b%Npoints
        IF ( ROdata%Lev2b%Temp(in) > 150.0 .AND. &
             ROdata%Lev2b%Temp(in) < 350.0 ) THEN
          nvalid = nvalid + 1
        END IF
      END DO
      IF ( nvalid < 1 ) nomet = .TRUE.

      IF ( nomet ) THEN
        IF ( ROdata%Lev2b%Npoints > 0 ) THEN
          ROdata%Lev2b%Npoints = 0
          ROdata%Lev2c%Npoints = 0
        END IF
      END IF

      IF ( ROdata%Lev2c%Geop_Sfc < -1000.0 .OR. &
           ROdata%Lev2c%Geop_Sfc > 10000.0 )    &
           ROdata%Lev2c%Npoints = 0

      ROdata%Lev2d%Npoints = 0

! Skip this profile if no valid bending angles, refractivity,
! met. or surface met. present

      IF ( ROdata%Lev1b%Npoints <= 0 .AND. &
           ROdata%Lev2a%Npoints <= 0 .AND. &
           ROdata%Lev2b%Npoints <= 0 .AND. &
           ROdata%Lev2c%Npoints <= 0 ) THEN
        CALL message ( msg_warn, "No. of L1b,2a,2b,2c samples"// &
                                 " all zero - skipping this profile" )
        CYCLE
      END IF

!--------------------------------------------------------------
! 7.4 Thin BA, N & T,q,p profiles as required; ensure all
!     profiles to be encoded are in ascending height order
!--------------------------------------------------------------

      CALL ropp_io_thin ( ROdata, Thindsn )

      IF ( .NOT. unordered ) THEN
        CALL message ( msg_diag, "Ensuring all profiles are in "// &
                                 "ascending height order..." )
        CALL ropp_io_ascend ( ROdata )
      END IF

!--------------------------------------------------------------
! 7.5 Calculate total number of BUFR elements for this profile
!     and allocate working arrays for BUFR-interface data values
!--------------------------------------------------------------

! No. of BUFR elements expected

      nElem = 37                                            &  ! Header
            +  1 + ROdata%Lev1b%Npoints * ( 5 + nfreq * 6 ) &  ! Level 1b
            +  1 + ROdata%Lev2a%Npoints * 6                 &  ! Level 2a
            +  1 + ROdata%Lev2b%Npoints * 10                &  ! Level 2b
            +  7                                               ! Level 2c

      ALLOCATE ( Values(1:nElem), STAT=status )
      IF ( status /= 0 ) THEN
        CALL message ( msg_fatal, "Failed to allocate memory for "// &
                                  "Values array" )
      END IF
      Values(:)  = 0.0_dp

! No. of BUFR Replication Factors expected

      nRepFac = ROdata%Lev1b%Npoints*nfreq + 3
      ALLOCATE ( RepFac(1:nRepFac), STAT=status )
      IF ( status /= 0 ) THEN
        CALL message ( msg_fatal, "Failed to allocate memory for "// &
                                  "RepFac array" )
      END IF

!--------------------------------------------------------------
! 7.6 Convert RO data to BUFR array
!--------------------------------------------------------------

      CALL ConvertROPPtoBUFR ( ROdata,                          &
                               CorrOnly,                        &
                               OrigICAO, OrigCentre, SubCentre, &
                               Values,   nValues,               &
                               RepFac,   nRepFac )

!--------------------------------------------------------------
! 7.7 Allocate working array for descriptors, based on thinned
!     number of data values, and allowing headroom for expansion.
!--------------------------------------------------------------

!      nExtra = MAX(minExtra,(nValues*5/10))
!      ALLOCATE ( Descr(1:nValues+nExtra), STAT=status )
!      IF ( status /= 0 ) THEN
!        CALL message ( msg_fatal, "Failed to allocate memory for "// &
!                                  "Descr array" )
!      ELSE
!        WRITE ( number, FMT="(I6)" ) nValues+nExtra
!        CALL message ( msg_diag, "  Allocated Descriptor space   :"// &
!                                 TRIM(number) )
!      END IF
!      Descr(:) = 0

!--------------------------------------------------------------
! 7.8 Diagnostics of what we're about to encode
!--------------------------------------------------------------

      CALL message ( msg_diag, "Encoding the following data: " )

      CALL MonthOfYear ( ROdata%DTocc%Month, MonthName, 1 )
      WRITE ( outmsg, FMT=DTfmt1 ) ROdata%DTocc%Hour,   &
                                   ROdata%DTocc%Minute, &
                                   ROdata%DTocc%Day,    &
                                   MonthName(1:3),      &
                                   ROdata%DTocc%Year
      CALL message ( msg_diag, "  Nominal time of occultation  : "//TRIM(outmsg) )

      IF ( ROdata%GeoRef%Lon > 180.0 ) &
             ROdata%GeoRef%Lon = ROdata%GeoRef%Lon - 360.0
      WRITE ( outmsg, FMT="(F6.2,',',F7.2)") ROdata%GeoRef%Lat, ROdata%GeoRef%Lon
      CALL message ( msg_diag, "  Nominal occ lat/lon location : "//TRIM(outmsg) )
      WRITE ( number, FMT="(I6)") ROdata%Lev1a%Npoints
      CALL message ( msg_diag, "  No. of orbit state vectors   : "//TRIM(number) )
      WRITE ( number, FMT="(I6)") ROdata%Lev1b%Npoints
      CALL message ( msg_diag, "  No. of bending angle samples : "//TRIM(number) )
      IF ( ROdata%Lev1b%Npoints > 0 ) THEN
        IF ( CorrOnly) THEN
          CALL message ( msg_diag, "  Bending angles present       : " // &
                                   "Corrected only" )
        ELSE
          CALL message ( msg_diag, "  Bending angles present       : " // &
                                   "L1+L2+Corrected" )
        ENDIF
      ENDIF
      WRITE ( number, FMT="(I6)") ROdata%Lev2a%Npoints
      CALL message ( msg_diag, "  No. of refractivity samples  : "//TRIM(number) )
      WRITE ( number, FMT="(I6)") ROdata%Lev2b%Npoints
      CALL message ( msg_diag, "  No. of geophysical samples   : "//TRIM(number) )
      WRITE ( number, FMT="(I6)") ROdata%Lev2c%Npoints
      CALL message ( msg_diag, "  No. of surface geo. samples  : "//TRIM(number) )
      WRITE ( number, FMT="(I6)") ROdata%Lev2d%Npoints
      CALL message ( msg_diag, "  No. of model coeff. levels   : "//TRIM(number) )

      IF ( nValues == nElem ) THEN
        WRITE ( number, FMT="(I6)" ) nElem
        CALL message ( msg_diag, "  Total no. of BUFR elements   :"//  &
                                 TRIM(number) )
      ELSE
        WRITE ( number, FMT="(I6)" ) nValues
        CALL message ( msg_diag, "  Thinned no. of BUFR elements :"// &
                                 TRIM(number) )
      END IF
!      WRITE ( number, FMT="(I6)" ) nValues+nExtra
!      CALL message ( msg_diag, "  Allocated Descriptor space   :"// &
!                               TRIM(number) )

!--------------------------------------------------------------
! 7.9 Encode this occultation & write it to output BUFR file
!--------------------------------------------------------------

      ! Have to include this if writing a BUFR file
      CALL codes_bufr_new_from_samples(BUFRunit,'BUFR4',iret)
      IF (iret/=CODES_SUCCESS) THEN
        STOP 1
      ENDIF

      ! Set replication factors
      IF(ALLOCATED(ivalues)) DEALLOCATE(ivalues)
      ALLOCATE(ivalues(ROdata%Lev1b%Npoints))
      ivalues(:) = nFreq
      CALL codes_set(BUFRunit,'inputDelayedDescriptorReplicationFactor',ivalues)

      IF(ALLOCATED(ivalues)) DEALLOCATE(ivalues)
      ALLOCATE(ivalues(3))

      ivalues=(/ ROdata%Lev1b%Npoints, ROdata%Lev2a%Npoints, ROdata%Lev2b%Npoints /)
      CALL codes_set(BUFRunit,'inputExtendedDelayedDescriptorReplicationFactor',ivalues)

      ! Set sections 0 to 3
      CALL codes_set(BUFRunit,'edition',Edition)               ! BUFR edition (3 or 4)
      CALL codes_set(BUFRunit,'masterTableNumber',MasterTable) ! BUFR Master Table used
      CALL codes_set(BUFRunit,'bufrHeaderCentre',OrigCentre)   ! Originating Centre code
      CALL codes_set(BUFRunit,'bufrHeaderSubCentre',SubCentre) ! Originating sub-centre
      CALL codes_set(BUFRunit,'updateSequenceNumber',0)        ! Update sequence number (original)
      CALL codes_set(BUFRunit,'dataCategory',DataType)         ! BUFR message type (Table A code)
      CALL codes_set(BUFRunit,'internationalDataSubCategory',IntlSubType) ! International sub-category

      CALL codes_set(BUFRunit,'dataSubCategory',LoclSubType) ! BUFR message subtype (local use)

      CALL codes_set(BUFRunit,'masterTablesVersionNumber',VerMasTable) ! Version no. of Master Table
      CALL codes_set(BUFRunit,'localTablesVersionNumber',VerLocTable)  ! Version no. of local table used

      CALL codes_set(BUFRunit,'typicalYear',ROdata%DTocc%Year)
      CALL codes_set(BUFRunit,'typicalMonth',ROdata%DTocc%Month)
      CALL codes_set(BUFRunit,'typicalDay',ROdata%DTocc%Day)
      CALL codes_set(BUFRunit,'typicalHour',ROdata%DTocc%Hour)
      CALL codes_set(BUFRunit,'typicalMinute',ROdata%DTocc%Minute)
      CALL codes_set(BUFRunit,'typicalSecond',NINT(ROdata%DTocc%Second + (ROdata%DTocc%msec / 1000.0)))

      CALL codes_set(BUFRunit,'numberOfSubsets',1) ! No. of subsets (observations)

      CALL codes_set(BUFRunit,'observedData',1) ! Data flags
      CALL codes_set(BUFRunit,'compressedData',0)

      ! Create the structure of the data section
      CALL codes_set(BUFRunit,'unexpandedDescriptors',310026)

      !!!!!!!!!!!!! Set Section 4 !!!!!!!!!!!!!!!!!!
      !-------------------------------------------------------------
      ! Satellite data introducer
      !-------------------------------------------------------------

      CALL codes_set(BUFRunit,'satelliteIdentifier',Values(1)) ! [001007] LEO ID

      CALL codes_set(BUFRunit,'satelliteInstruments',Values(2)) ! [002019] RO Instrument

      CALL codes_set(BUFRunit,'centre',Values(3)) ! [001033] B/g or Proc. centre 

      CALL codes_set(BUFRunit,'productTypeForRetrievedAtmosphericGases',Values(4))! [002172] Product type (limb sounding)

      CALL codes_set(BUFRunit,'softwareIdentification',Values(5)) ! [025060] Software version

      ! Date/time of start of occultation (or background profile)

      CALL codes_set(BUFRunit,'timeSignificance',Values(6)) ! [008021] Time.sig (start)
      CALL codes_set(BUFRunit,'year',Values(7))    ! [004001] Year
      CALL codes_set(BUFRunit,'month',Values(8))   ! [004002] Month
      CALL codes_set(BUFRunit,'day',Values(9))     ! [004003] Day
      CALL codes_set(BUFRunit,'hour',Values(10))   ! [004004] Hour
      CALL codes_set(BUFRunit,'minute',Values(11)) ! [004005] Minute
      CALL codes_set(BUFRunit,'second',Values(12)) ! [004006] Seconds & MSecs
      ! Summary quality information

      CALL codes_set(BUFRunit,'radioOccultationDataQualityFlags',Values(13))
      CALL codes_set(BUFRunit,'#1#percentConfidence',Values(14)) ! [033007] Percent confidence

      ! location and velocity of LEO and GNSS
      ! [027031, 028031, 010031] = LEO (X, Y, Z) posn (m)
      CALL codes_set(BUFRunit,'#1#DistanceFromEarthCentreInDirectionOf0DegreesLongitude',Values(15))
      CALL codes_set(BUFRunit,'#1#DistanceFromEarthCentreInDirection90DegreesEast',Values(16))
      CALL codes_set(BUFRunit,'#1#DistanceFromEarthCentreInDirectionOfNorthPole',Values(17))
      ! [001041, 001042, 001043] = LEO (X, Y, Z) vely (m/s)
      CALL codes_set(BUFRunit,'#1#absolutePlatformVelocityFirstComponent',Values(18))
      CALL codes_set(BUFRunit,'#1#absolutePlatformVelocitySecondComponent',Values(19))
      CALL codes_set(BUFRunit,'#1#absolutePlatformVelocityThirdComponent',Values(20))

      CALL codes_set(BUFRunit,'satelliteClassification',Values(21))  ! [002020] GNSS class
      CALL codes_set(BUFRunit,'platformTransmitterIdNumber',Values(22)) ! [001050] GNSS PRN

      ! [027031, 028031, 010031] = GNSS (X, Y, Z) posn (m)
      CALL codes_set(BUFRunit,'#2#DistanceFromEarthCentreInDirectionOf0DegreesLongitude',Values(23))
      CALL codes_set(BUFRunit,'#2#DistanceFromEarthCentreInDirection90DegreesEast',Values(24))
      CALL codes_set(BUFRunit,'#2#DistanceFromEarthCentreInDirectionOfNorthPole',Values(25))
      ! [001041, 001042, 001043] = GNSS (X, Y, Z) vely (m/s)
      CALL codes_set(BUFRunit,'#2#absolutePlatformVelocityFirstComponent',Values(26))
      CALL codes_set(BUFRunit,'#2#absolutePlatformVelocitySecondComponent',Values(27))
      CALL codes_set(BUFRunit,'#2#absolutePlatformVelocityThirdComponent',Values(28))

      ! Local Earth parameters
      CALL codes_set(BUFRunit,'timeIncrement',Values(29))  ! [004016] Time/start (s)
      CALL codes_set(BUFRunit,'#1#latitude',Values(30)) ! [005001] Latitude (deg)
      CALL codes_set(BUFRunit,'#1#longitude',Values(31)) ! [006001] Longitude (deg)
      ! [027031] CofC X, [028031] CofC Y, [010031] CofC Z (m)
      CALL codes_set(BUFRunit,'#3#DistanceFromEarthCentreInDirectionOf0DegreesLongitude',Values(32))
      CALL codes_set(BUFRunit,'#3#DistanceFromEarthCentreInDirection90DegreesEast',Values(33))
      CALL codes_set(BUFRunit,'#3#DistanceFromEarthCentreInDirectionOfNorthPole',Values(34))
      CALL codes_set(BUFRunit,'earthLocalRadiusOfCurvature',Values(35))  ! [010035] Radius value (m)
      CALL codes_set(BUFRunit,'#1#bearingOrAzimuth',Values(36)) ! [005021] Line of sight bearing (degT)
      CALL codes_set(BUFRunit,'geoidUndulation',Values(37)) ! [010036] Geoid undulation (m)

      !-------------------------------------------------------------
      ! Level 1b data (bending angle profile)
      !-------------------------------------------------------------

      DO i=0,ROdata%Lev1b%Npoints-1 !,43+(6*nfreq*ROdata%Lev1b%Npoints),6

        n = 39 + i * (5 + (6 * nFreq))

        ! [005001] Latitude (deg)
        WRITE(descriptor,'(A,I0,A)') '#',i+2,'#latitude' ! I0 adjusts width for int
        CALL codes_set(BUFRunit,descriptor,Values(n))

        ! [006001] Longitude (deg)
        WRITE(descriptor,'(A,I0,A)') '#',i+2,'#longitude'
        CALL codes_set(BUFRunit,descriptor,Values(n+1))

        ! [005021]  GNSS-to-LEO line of sight bearing
        WRITE(descriptor,'(A,I0,A)') '#',i+2,'#bearingOrAzimuth'
        CALL codes_set(BUFRunit,descriptor,Values(n+2))

        DO j=0,nFreq-1 ! nFreq can be 3 (L1, L2, Corrected) or 1 (Corrected only)

          ! meanFrequency
          WRITE(descriptor,'(A,I0,A)') '#',nFreq*i+j+1,'#meanFrequency'
          CALL codes_set(BUFRunit,descriptor,Values(n+4+6*j)) ! [002121] L1=1.5Ghz

          ! impactParameter
          WRITE(descriptor,'(A,I0,A)') '#',nFreq*i+j+1,'#impactParameter'
          CALL codes_set(BUFRunit,descriptor,Values(n+5+6*j)) ! [007040] Impact parameter (m)

          ! bendingAngle
          WRITE(descriptor,'(A,I0,A)') '#',2*nFreq*i+2*j+1,'#bendingAngle'
          CALL codes_set(BUFRunit,descriptor,Values(n+6+6*j)) ! [015037] B/angle (rad)

          ! firstOrderStatistics
          WRITE(descriptor,'(A,I0,A)') '#',2*nFreq*i+2*j+1,'#firstOrderStatistics'
          CALL codes_set(BUFRunit,descriptor,Values(n+7+6*j)) ! [008023] 1st order stats (rms)

          ! error in bendingAngle
          WRITE(descriptor,'(A,I0,A)') '#',2*nFreq*i+2*j+2,'#bendingAngle'
          CALL codes_set(BUFRunit,descriptor,Values(n+8+6*j)) ! [015037] B/angle error (rad)

          ! firstOrderStatistics
          WRITE(descriptor,'(A,I0,A)') '#',2*nFreq*i+2*j+2,'#firstOrderStatistics'
          CALL codes_set(BUFRunit,descriptor,Values(n+9+6*j)) ! [008023] 1st order stats (off)

        END DO

        ! percentConfidence
        WRITE(descriptor,'(A,I0,A)') '#',i+2,'#percentConfidence'  
        CALL codes_set(BUFRunit,descriptor,Values(n + 4 + 6*nFreq)) ! [033007] Percent confidence

      END DO

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

      DO i=0,ROdata%Lev2a%Npoints-1

        n = ROdata%Lev1b%Npoints * (6*nFreq + 5) + 40 + 6*i

        WRITE(descriptor,'(A,I0,A)') '#',i+1,'#height'
        CALL codes_set(BUFRunit,descriptor,Values(n)) ! [007007] Height amsl (m)

        WRITE(descriptor,'(A,I0,A)') '#',2*i+1,'#atmosphericRefractivity'
        CALL codes_set(BUFRunit,descriptor,Values(n+1)) ! [015036] Refrac (N-units)

        WRITE(descriptor,'(A,I0,A)') '#',(ROdata%Lev1b%Npoints * nFreq * 2) + 2*i + 1 , &
                                     '#firstOrderStatistics'
        CALL codes_set(BUFRunit,descriptor,Values(n+2)) ! [008023] 1st order stats (rms)

        WRITE(descriptor,'(A,I0,A)') '#',2*i+2,'#atmosphericRefractivity'
        CALL codes_set(BUFRunit,descriptor,Values(n+3)) ! [015036] Refrac error (N-units)

        WRITE(descriptor,'(A,I0,A)') '#',(ROdata%Lev1b%Npoints * nFreq * 2) + 2*i + 2, &
                                     '#firstOrderStatistics'
        CALL codes_set(BUFRunit,descriptor,Values(n+4)) ! [008023] 1st order stats (off)

        WRITE(descriptor,'(A,I0,A)') '#',ROdata%Lev1b%Npoints+2+i,'#percentConfidence'
        CALL codes_set(BUFRunit,descriptor,Values(n+5)) ! [033007] Percent confidence

      END DO

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

      DO i=0,ROdata%Lev2b%Npoints-1

        n =  ROdata%Lev1b%Npoints * (6*nFreq + 5) + 41 + ROdata%Lev2a%Npoints*6 + 10*i

        WRITE(descriptor,'(A,I0,A)') '#',i+1,'#geopotentialHeight'
        CALL codes_set(BUFRunit,descriptor,Values(n)) ! [007009] Geopot ht (gpm)

        WRITE(descriptor,'(A,I0,A)') '#',2*i+1,'#nonCoordinatePressure'
        CALL codes_set(BUFRunit,descriptor,Values(n+1)) ! [010004] Pressure (Pa)

        WRITE(descriptor,'(A,I0,A)') '#',2*i+1,'#airTemperature'
        CALL codes_set(BUFRunit,descriptor,Values(n+2)) ! [012001] Temperature (K)

        WRITE(descriptor,'(A,I0,A)') '#',2*i+1,'#specificHumidity'
        CALL codes_set(BUFRunit,descriptor,Values(n+3)) ! [013001] Spec/humidity (Kg/Kg)

        WRITE(descriptor,'(A,I0,A)') '#',2*nFreq*ROdata%Lev1b%Npoints+2*ROdata%Lev2a%Npoints+2*i+1, &
                                     '#firstOrderStatistics'
        CALL codes_set(BUFRunit,descriptor,Values(n+4)) ! [008023] 1st order stats (rms)

        WRITE(descriptor,'(A,I0,A)') '#',2*i+2,'#nonCoordinatePressure'
        CALL codes_set(BUFRunit,descriptor,Values(n+5)) ! [010004] Pressure error (Pa)

        WRITE(descriptor,'(A,I0,A)') '#',2*i+2,'#airTemperature'
        CALL codes_set(BUFRunit,descriptor,Values(n+6)) ! [012001] Temperature error (K)

        WRITE(descriptor,'(A,I0,A)') '#',2*i+2,'#specificHumidity'
        CALL codes_set(BUFRunit,descriptor,Values(n+7)) ! [013001] S/Hum error (Kg/Kg)

        WRITE(descriptor,'(A,I0,A)') '#',2*nFreq*ROdata%Lev1b%Npoints+2*ROdata%Lev2a%Npoints+2*i+2, &
                                     '#firstOrderStatistics'
        CALL codes_set(BUFRunit,descriptor,Values(n+8)) ! [008023] 1st order stats (off)

        WRITE(descriptor,'(A,I0,A)') '#',ROdata%Lev1b%Npoints+ROdata%Lev2a%Npoints+2+i, &
                                     '#percentConfidence'
        CALL codes_set(BUFRunit,descriptor,Values(n+9)) ! [033007] Percent confidence

      END DO

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

      n = 39 + &                                   ! header
          ROdata%Lev1b%Npoints * (6*nFreq + 5) + & ! num lev 1b values
          ROdata%Lev2a%Npoints *6 + 1 + &          ! num lev 2a values
          ROdata%Lev2b%Npoints *10 + 1             ! num lev 2b values

      ! [008003] Vertical sig. (surf)
      CALL codes_set(BUFRunit,'verticalSignificanceSatelliteObservations',Values(n))

      WRITE(descriptor,'(A,I0,A)') '#',ROdata%Lev2b%Npoints+1,'#geopotentialHeight'
      CALL codes_set(BUFRunit,descriptor,Values(n+1)) ! [007009] Geoptot.Ht. (of surf)

      WRITE(descriptor,'(A,I0,A)') '#',2*ROdata%Lev2b%Npoints+1,'#nonCoordinatePressure'
      CALL codes_set(BUFRunit,descriptor,Values(n+2)) ! [010004] Surface pressure (Pa)

      WRITE(descriptor,'(A,I0,A)') '#',2*nFreq*ROdata%Lev1b%Npoints+2*ROdata%Lev2a%Npoints+2*ROdata%Lev2b%Npoints+1, &
                                   '#firstOrderStatistics'
      CALL codes_set(BUFRunit,descriptor,Values(n+3)) ! [008023] 1st order stats (rms)

      WRITE(descriptor,'(A,I0,A)') '#',2*ROdata%Lev2b%Npoints+2,'#nonCoordinatePressure'
      CALL codes_set(BUFRunit,descriptor,Values(n+4)) ! [010004] S/press error (Pa)

      WRITE(descriptor,'(A,I0,A)') '#',2*nFreq*ROdata%Lev1b%Npoints+2*ROdata%Lev2a%Npoints+2*ROdata%Lev2b%Npoints+2, &
                                   '#firstOrderStatistics'
      CALL codes_set(BUFRunit,descriptor,Values(n+5)) ! [008023] 1st order stats (off)

      WRITE(descriptor,'(A,I0,A)') '#',ROdata%Lev1b%Npoints+ROdata%Lev2a%Npoints+ROdata%Lev2b%Npoints+2, &
                                   '#percentConfidence'
      CALL codes_set(BUFRunit,descriptor,Values(n+6)) ! [033007] Percent confidence

      !-------------------------------------------------------------
      ! Write to file and deallocate arrays
      !-------------------------------------------------------------

      outfile=TRIM(BUFRdsn)
      CALL codes_set(BUFRunit,'pack',1)
      IF (imessage .eq. 1) then
        CALL codes_open_file(ioutfile,outfile,'w')
      ELSE
        CALL codes_open_file(ioutfile,outfile,'a')
      ENDIF
      CALL codes_write(BUFRunit,ioutfile)
      CALL codes_close_file(ioutfile)
      CALL codes_release(BUFRunit)

      LenTot = LenTot + LenBUF
      IF ( LenBUF > 0 ) nmsg = nmsg + 1

      IF ( ALLOCATED ( Values ) ) DEALLOCATE ( Values )
      IF ( ALLOCATED ( Descr  ) ) DEALLOCATE ( Descr )
      IF ( ALLOCATED ( RepFac ) ) DEALLOCATE ( RepFac )

      IF ( GTShdrType /= NOhdrs ) THEN

        WRITE ( GTShdrTxt, '(A3,A1,A2,1x,A4,1x,I0.2,I0.2,I0.2)' ) &
               'IUT', TRIM(GTShdrGAD(INT(ROdata%GeoRef%Lat), INT(ROdata%GeoRef%Lon))), &
               '14', OrigICAO, ROdata%DTocc%Day, ROdata%DTocc%Hour, ROdata%DTocc%Minute

        IF ( GTShdrType == ARhdrs ) THEN

          CALL message(msg_noin, '')
          CALL message(msg_info, &
          'To wrap BUFR file in GTS ARH headers, run:')
          CALL message(msg_cont, &
            'perl $ROPP_SRC/ropp_io/tools/gtsheaders_bufr.pl ' &
            // TRIM(ADJUSTL(BUFRdsn)) // ' ' &
            // TRIM(ADJUSTL(BUFRdsn)) // '_gts (for example) ' // GTShdrTxt )
          CALL message(msg_noin, '')

        ELSE IF ( GTShdrType == IPhdrs ) THEN

          CALL message(msg_noin, '')
          CALL message(msg_info, &
          'To wrap BUFR file in GTS ARH + IPH headers, run:')
          CALL message(msg_cont, &
            'perl $ROPP_SRC/ropp_io/tools/gtsheaders_bufr.pl --iph ' &
            // TRIM(ADJUSTL(BUFRdsn)) // ' ' &
            // TRIM(ADJUSTL(BUFRdsn)) // '_gts (for example) ' // GTShdrTxt )
          CALL message(msg_noin, '')

        ENDIF

      ENDIF

    END DO                       ! end of profiles loop

!--------------------------------------------------------------
! 7.10 Free memory ready for next file
!--------------------------------------------------------------

    CALL ropp_io_free ( ROdata )

  END DO                       ! end of file loop

!--------------------------------------------------------------
! 8. Tidy up & finish
!--------------------------------------------------------------

  IF ( ALLOCATED(EUMdsn) ) DEALLOCATE (EUMdsn)

  CALL message ( msg_noin, " " )
  CALL EXIT(msg_exit_status)

END PROGRAM eum2bufr_eccodes

!-------------------------------------------------------------------------------
! 9. Version information
!-------------------------------------------------------------------------------

SUBROUTINE version_info()
  USE ropp_io, ONLY : ropp_io_version
  CHARACTER (LEN=40) :: version
  version = ropp_io_version()
  PRINT *, 'eum2bufr - EUM netCDF4 to BUFR encoder [ECMWF ecCodes library]'
  PRINT *, ''
  PRINT *, 'This program is part of ROPP (IO) Release ' // TRIM(version)
  PRINT *, ''
END SUBROUTINE version_info