! $Id: grib2bgrasc.f90 4130 2014-05-01 08:05:56Z idculv $

PROGRAM grib2bgrasc

!****x* Programs/grib2bgrasc *
!
! NAME
!   grib2bgrasc    (grib2bgrasc.f90)
!
! SYNOPSIS
!   Extract RO background data from GRIB (ed 1 or 2) file(s)
!   and output in fortran namelist format.
!
!   > grib2bgrasc grib_file -lat <lat> -lon <lon>
!                 [-date <date>] [-time <time>]
!                 [-g grib_file2]
!                 [-z grib_filez]
!                 [-y obs_file]
!                 [-o <output namelist file>]
!                 [-x <output netCDF file>]
!                 [-r <radius_of_curvature>]
!                 [-u <undulation>]
!                 [-a <azimuth>]
!                 [-c <centre of curvature>]
!                 [-d] [-h] [-v]
!
! ARGUMENTS
!   grib_file  - GRIB file containing background fields.
!   lat        - Latitude (deg)
!   lon        - Longitude (deg)
!
! OPTIONS
!   Option switches can be in any order and are case-insensitive;
!   any space(s) between a switch and its (madatory) argument is
!   optional.
!     -date <date> (yyyymmdd). Date, if time interpolation needed.
!     -time <time> (HHMM or HHMMSS). Time, if time interpolation needed.
!     -g grib_file2. Second GRIB file containing background fields, if time interpolation needed.
!     -z grib_filez.  GRIB file containing surface geopotential, aka g_wmo*orography, if needed.
!     -y obs_file.  Observation file containing date, time, lat, lon, roc and und for which background is desired.
!     -o <namelist_file> Output file name.
!     -x to output ROPP netCDF file.
!     -r <roc>. Radius of curvature of tangent plane (m).
!     -u <und>. Undulation at tangent point (m).
!     -a <azi>. Azimuth at tangent point (m).
!     -c <coc>. Centre of curvature (3 components) (m).
!     -d to output additional diagnostics.
!     -h help.
!     -v version information.
!   Defaults:
!     Input file name : required
!     Output file name : $PWD/input_file.nml
! NOTE
!   The environment variables GEOPOT_COEF and GEOPOT_CORR must be
!   set appropriately if the user wants the undulation to be calculated.
!
! INPUTS
!   grib_files containing background fields in GRIB (ed 1 or 2) format.
!     Each input file must contain {p*, phi*, T, q, Ak, Bk}.
!   obs file must be in ROPP netCDF format - see ROPP UG I: I/O for details.
!
! OUTPUTS
!   Ascii file holding background profile in Fortan namelist format.
!     The output file name is optional; if unspecified, it is derived from the GRIB file name.
!
! CALLS
!   IARGC
!   message
!   extract_prof_from_GRIB
!   generate_GPH
!
! MODULES
!   grib_api
!   typesizes
!   messages
!   coordinates
!   DateTimeProgs
!   ropp_utils
!   ropp_io
!   ropp_io_ecmwf
!
! ERRORS
!   Program (shell) return codes:
!     0 = OK
!     1 = At least one Warning occurred
!     2 = At least one Error occurred
!     3 = A Fatal error occurred
!
! DEPENDENCIES
!   ECMWF GRIB package  - GRIB kernel routines
!   ROPP Utils library  - ROPP utility routines
!   The environment variables GEOPOT_COEF and GEOPOT_CORR must be 
!     set appropriately if the user wants the undulation to be calculated.
!
! DESCRIPTION
!   A GRIB background data translator for Radio Occultation data.
!   Generates one profile, in Fortran namelist format, from one GRIB file.
!
! SEE ALSO
!   bgrasc2ropp
!
! REFERENCES
!   1) GRIB API
!      https://software.ecmwf.int/wiki/display/GRIB/Home
!   2) ROPP User Guide - Part 1: IO module
!      SAF/ROM/METO/UG/ROPP/002
!
! AUTHOR
!   Met Office, Exeter, UK and DMI, Copenhagen, Denmark.
!   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,     wp => EightByteReal
  USE messages
  USE DateTimeProgs, ONLY: CalToJul
  USE ropp_utils,    ONLY: ropp_MDFV, ropp_MDTV, ropp_MIFV, Get_IO_Unit, isinrange
  USE ropp_io,       ONLY: ropp_io_read, ropp_io_write, ropp_io_free, &
                           ropp_io_nrec, ropp_io_rangecheck, ropp_io_version
  USE ropp_io_ecmwf
  USE ropp_io_types, ONLY: ROprof, PCD_occultation

  IMPLICIT NONE

! Input parameters
  CHARACTER(LEN=256)                :: file_in            ! Input file
  CHARACTER(LEN=256)                :: file_in2=' '       ! Optional input file (if interpolating in time)
  CHARACTER(LEN=256)                :: file_inz=' '       ! Optional input file (if external orography source needed)
  CHARACTER(LEN=256)                :: file_obs=' '       ! Optional input file (holding desired date, time, lat and lon)
  CHARACTER(LEN=256)                :: file_out_nml       ! Output namelist file
  CHARACTER(LEN=256)                :: file_out_nc        ! Output netCDF file
  CHARACTER(LEN=8)                  :: datestring=' '     ! Input date
  CHARACTER(LEN=8)                  :: timestring=' '     ! Input time
  INTEGER                           :: year_in=ropp_MIFV  ! Input year
  INTEGER                           :: mon_in=ropp_MIFV   ! Input month
  INTEGER                           :: day_in=ropp_MIFV   ! Input day
  INTEGER                           :: hour_in=ropp_MIFV  ! Input hour
  INTEGER                           :: min_in=ropp_MIFV   ! Input minute
  INTEGER                           :: sec_in=ropp_MIFV   ! Input second
  INTEGER                           :: msec_in=0          ! Input millisecond
  REAL(wp)                          :: lon_in=ropp_MDFV   ! Input lon (deg)
  REAL(wp)                          :: lat_in=ropp_MDFV   ! Input lat (deg)
  REAL(wp)                          :: roc_in=ropp_MDFV   ! Input radius of curvature (m)
  REAL(wp)                          :: und_in=ropp_MDFV   ! Input undulation (m)
  REAL(wp)                          :: azi_in=ropp_MDFV   ! Input azimuth (deg)
  REAL(wp), DIMENSION(3)            :: coc_in=(/ ropp_MDFV, &
                                                 ropp_MDFV, &
                                                 ropp_MDFV /) ! Input centre of curvature (m)

! Local variables
  TYPE(ROprof)                      :: data_obs           ! Holds file_obs data
  TYPE(ROprof)                      :: grib1              ! Holds file_in data
  TYPE(ROprof)                      :: grib2              ! Holds file_in2 data
  TYPE(ROprof)                      :: gribz              ! Holds file_inz data
  CHARACTER(LEN=4)                  :: syear_in           ! Input year
  CHARACTER(LEN=2)                  :: smon_in            ! Input month
  CHARACTER(LEN=2)                  :: sday_in            ! Input day
  CHARACTER(LEN=2)                  :: shour_in           ! Input hour
  CHARACTER(LEN=2)                  :: smin_in            ! Input minute
  CHARACTER(LEN=2)                  :: ssec_in            ! Input second
  CHARACTER(LEN=12)                 :: slon_in            ! Input lon (deg)
  CHARACTER(LEN=12)                 :: slat_in            ! Input lat (deg)
  CHARACTER(LEN=20)                 :: sroc_in            ! Input radius of curvature (m)
  CHARACTER(LEN=12)                 :: sund_in            ! Input undulation (m)
  CHARACTER(LEN=12)                 :: sazi_in            ! Input azimuth (deg)
  CHARACTER(LEN=12)                 :: scoc_in            ! Input centre of curvature (m)
  CHARACTER(LEN=4)                  :: sleo_id_in         ! Input leo_id
  CHARACTER(LEN=4)                  :: sstn_id_in         ! Input stn_id
  CHARACTER(LEN=4)                  :: sgns_id_in         ! Input gns_id
  CHARACTER(LEN=32)                 :: socc_id_in         ! Input occ_id

  CHARACTER(LEN=256)                :: arg                ! Command line argument
  CHARACTER(LEN=256)                :: GEOPOT_COEF_FILE   ! EGM96 coeff file name
  CHARACTER(LEN=256)                :: GEOPOT_CORR_FILE   ! EGM96 corrn file name
  LOGICAL                           :: exists             ! File existence flag
  INTEGER                           :: namelist_lun       ! Namelist logical unit
  INTEGER                           :: i, ii              ! Loop counters / indices
  INTEGER                           :: iarg               ! Command line argument index
  INTEGER                           :: narg               ! No. of command line arguments
  INTEGER                           :: iostatus           ! I/O Status
  INTEGER                           :: n_profs            ! No. of RO profiles
  INTEGER, DIMENSION(8)             :: CDT                ! Date/time combo for cal2jul
  REAL(wp), DIMENSION(2)            :: time_coeff         ! Time interpolation coefficients
  REAL(KIND(1.D0))                  :: jul_day1           ! Julian day of file_in
  REAL(KIND(1.D0))                  :: jul_day2           ! Julian day of file_in2
  REAL(KIND(1.D0))                  :: jul_day_in         ! Julian day of requested date/time
  REAL(wp)                          :: swver              ! Software version
  REAL(wp)                          :: dtor               ! pi/180
  LOGICAL                           :: ranchk = .TRUE.    ! Range check when reading/writing ROPP files

! Namelist variables
  INTEGER, PARAMETER                :: nlevs_max=200      ! Maximum number of model levels
  INTEGER                           :: nlevs              ! Number of full model levels
  INTEGER                           :: year, mon, day     ! Background validity date
  INTEGER                           :: hour, min, sec     ! Background validity time
  REAL(wp)                          :: tlead=ropp_MDFV    ! Background forecast range (hr)
  REAL(wp)                          :: lon=ropp_MDFV      ! Longitude (deg)
  REAL(wp)                          :: lat=ropp_MDFV      ! Latitude (deg)
  REAL(wp)                          :: und=ropp_MDFV      ! Undulation (m)
  REAL(wp)                          :: roc=ropp_MDFV      ! Radius of curvature (m)
  REAL(wp), DIMENSION(3)            :: coc=(/ ropp_MDFV, &
                                              ropp_MDFV, &
                                              ropp_MDFV /)! Centre of curvature (m)
  REAL(wp)                          :: azi=ropp_MDFV      ! GNSS->LEO line of sight angle (degT)
  REAL(wp)                          :: Z0=ropp_MDFV       ! Surface geopotential height (m)
  REAL(wp)                          :: P0=ropp_MDFV       ! Surface pressure (Pa)
  REAL(wp), DIMENSION(nlevs_max)    :: P=ropp_MDFV        ! Pressure (Pa)
  REAL(wp), DIMENSION(nlevs_max)    :: T=ropp_MDFV        ! Temperature (K)
  REAL(wp), DIMENSION(nlevs_max)    :: Q=ropp_MDFV        ! Specific humidity (g/kg)
  REAL(wp), DIMENSION(nlevs_max)    :: Z=ropp_MDFV        ! Geopotential height (m)
  REAL(wp), DIMENSION(nlevs_max+1)  :: Ak=ropp_MDFV       ! Hybrid/Eta level A-coefficient (Pa)
  REAL(wp), DIMENSION(nlevs_max+1)  :: Bk=ropp_MDFV       ! Hybrid/Eta level B-coefficient
  NAMELIST / bgr_profile / year, mon, day, hour, min, sec, tlead, &
                           lon, lat, und, roc, coc, azi, &
                           Z0, P0, &
                           nlevs, &
                           P, T, Q, Z, Ak, Bk

! Functions

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


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

  CALL message_set_routine ( 'grib2bgrasc' )

  CALL message(msg_noin, '')
  CALL message(msg_noin, &
       '---------------------------------------------------------------------')
  CALL message(msg_noin, &
       '                     GRIB to background ascii converter'              )
  CALL message(msg_noin, &
       '---------------------------------------------------------------------')
  CALL message(msg_noin, '')

! Constants

  dtor = ATAN(1.0_wp) / 45.0_wp

  namelist_lun = Get_IO_Unit()

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

  narg = IARGC()

  file_in      = ' '       ! no default for i/p file name
  file_out_nml = ' '       ! assume a default generated from i/p file name
  file_out_nc  = ' '       ! assume a default generated from i/p file name

  lon_in = ropp_MDFV
  lat_in = ropp_MDFV

  iarg = 1

  DO WHILE ( iarg <= narg )

    CALL GETARG ( iarg, arg )

    SELECT CASE (arg)

      CASE ('-d', '-D', '--debug')
        msg_MODE = VerboseMode

      CASE ('-h', '-H', '--help' ,'?')
        narg = 0
        file_in = 'dummy'

      CASE ('-o', '-O', '--output')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        file_out_nml = arg

      CASE ('-lon', '-Lon', '-LON')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        READ (arg, *) lon_in

      CASE ('-lat', '-Lat', '-LAT')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        READ (arg, *) lat_in

      CASE ('-date')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        datestring = arg

      CASE ('-time')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        timestring = arg

      CASE ('-v', '-V', '--version')
        CALL version_info()
        CALL EXIT(msg_exit_ok)

      CASE ('-g', '-G')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        file_in2 = arg

      CASE ('-z', '-Z')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        file_inz = arg

      CASE ('-y', '-Y')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        file_obs = arg

      CASE ('-x', '-X')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        file_out_nc = arg

      CASE ('-r', '-R', '-roc', '-RoC', '-ROC')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        READ (arg, *) roc_in

      CASE ('-u', '-U', '-und', '-Und', '-UND')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        READ (arg, *) und_in

      CASE ('-a', '-A', '-azi', '-Azi', '-AZI')
        iarg = iarg + 1
        CALL GETARG ( iarg, arg )
        READ (arg, *) azi_in

      CASE ('-c', '-C', '-coc', '-CoC', '-COC')
        DO i=1,3
          iarg = iarg + 1
          CALL GETARG ( iarg, arg )
          READ (arg, *) coc_in(i)
        ENDDO

      CASE ('--no-ranchk')
        ranchk = .FALSE.

      CASE DEFAULT
         IF ( arg(1:1) /= '-' ) THEN
           file_in = arg
         END IF

    END SELECT

    iarg = iarg + 1

  END DO

  IF ( file_in == ' ' ) THEN
    CALL message ( msg_error, 'No input file(s) specified' )
    narg = 0
  END IF

  IF ( narg == 0 ) THEN
    CALL Usage
    CALL EXIT(msg_exit_status)
  ENDIF


!-------------------------------------------------------------
! 3. Input data preprocessing
!-------------------------------------------------------------

! 3.1  Check existence of input file(s) and environment variables
!----------------------------------------------------------------

  IF ( file_in == ' ' ) THEN
    CALL message ( msg_fatal, 'No input file specified' )
    narg = 0
  ENDIF

  INQUIRE ( FILE=file_in, EXIST=exists )
  IF ( .NOT. exists ) &
    CALL message ( msg_fatal, 'GRIB input file ' // TRIM(file_in) // &
                              ' not found' )

  IF (file_in2 /= ' ') THEN
    INQUIRE ( FILE=file_in2, EXIST=exists )
    IF ( .NOT. exists ) &
      CALL message ( msg_fatal, 'GRIB input file ' // TRIM(file_in2) // &
                                ' not found' )
  ENDIF

  IF (file_inz /= ' ') THEN
    INQUIRE ( FILE=file_inz, EXIST=exists )
    IF ( .NOT. exists ) &
      CALL message ( msg_fatal, 'GRIB input file ' // TRIM(file_inz) // &
                                ' not found' )
  ENDIF

  IF (file_obs /= ' ') THEN
    INQUIRE ( FILE=file_obs, EXIST=exists )
    IF ( .NOT. exists ) &
      CALL message ( msg_fatal, 'obs input file ' // TRIM(file_obs) // &
                                ' not found' )
  ENDIF

  IF (file_out_nml == ' ') THEN
    ii = INDEX(file_in, '/', BACK=.TRUE.)
    file_out_nml = file_in(ii+1:)
    i = INDEX(file_out_nml, '.', BACK=.TRUE.)
    IF (i > 0) THEN
      file_out_nml = TRIM(ADJUSTL(file_out_nml(1:i))) // 'nml'
    ELSE
      file_out_nml = TRIM(ADJUSTL(file_out_nml)) // '.nml'
    ENDIF
  ENDIF

  IF (file_out_nc == ' ') THEN
    ii = INDEX(file_in, '/', BACK=.TRUE.)
    file_out_nc = file_in(ii+1:)
    i = INDEX(file_out_nc, '.', BACK=.TRUE.)
    IF (i > 0) THEN
      file_out_nc = TRIM(ADJUSTL(file_out_nc(1:i))) // 'nc'
    ELSE
      file_out_nc = TRIM(ADJUSTL(file_out_nc)) // '.nc'
    ENDIF
  ENDIF

  GEOPOT_COEF_FILE = ' '
  CALL getenv ( 'GEOPOT_COEF', GEOPOT_COEF_FILE )
  IF (GEOPOT_COEF_FILE == ' ') &
    CALL message(msg_warn, '$GEOPOT_COEF environment file not set; ' // &
                           'undulation will not be calculable \n' // &
                           '(though may still be provided on command line or ' // &
                           'observation file)')

  GEOPOT_CORR_FILE = ' '
  CALL getenv ( 'GEOPOT_CORR', GEOPOT_CORR_FILE )
  IF (GEOPOT_CORR_FILE == ' ') &
    CALL message(msg_warn, '$GEOPOT_CORR environment file not set; ' // &
                           'undulation will not be calculable. \n' // &
                           '(though may still be provided on command line or ' // &
                           'observation file)')


! 3.2  Read information from obs file
!------------------------------------

  IF (file_obs /= ' ') THEN

    n_profs = ropp_io_nrec(file_obs)

    IF (n_profs > 1) &
      CALL message(msg_warn, 'More than one profile detected in obs file; ' // &
                             'taking parameters from first.')

    CALL ropp_io_read(data_obs, file_obs, rec=1, ranchk=ranchk)

! Use the values in the file if they're not specified at the command line

    IF (lon_in < ropp_MDTV) lon_in = data_obs%georef%lon
    IF (lat_in < ropp_MDTV) lat_in = data_obs%georef%lat
    IF (roc_in < ropp_MDTV) roc_in = data_obs%georef%roc
    IF (und_in < ropp_MDTV) und_in = data_obs%georef%undulation
    IF (azi_in < ropp_MDTV) azi_in = data_obs%georef%azimuth
    DO i=1,3
      IF (ABS(coc_in(i) - ropp_MDFV) < 1.0e3_wp) coc_in(i) = data_obs%georef%r_coc(i)
    ENDDO

    sleo_id_in = data_obs%leo_id
    sstn_id_in = data_obs%stn_id
    sgns_id_in = data_obs%gns_id
    socc_id_in = data_obs%occ_id

    IF (datestring == ' ') THEN
      WRITE(syear_in, '(i4.4)') data_obs%dtocc%year
      WRITE(smon_in,  '(i2.2)') data_obs%dtocc%month
      WRITE(sday_in,  '(i2.2)') data_obs%dtocc%day
      datestring = syear_in // smon_in // sday_in
    ENDIF

    IF (timestring == ' ') THEN
      WRITE(shour_in, '(i2.2)') data_obs%dtocc%hour
      WRITE(smin_in,  '(i2.2)') data_obs%dtocc%minute
      WRITE(ssec_in,  '(i2.2)') data_obs%dtocc%second
      timestring = shour_in // smin_in // ssec_in
      IF ( isinrange(data_obs%dtocc%msec, data_obs%dtocc%range%msec) ) &
        msec_in = data_obs%dtocc%msec
    ENDIF

    CALL ropp_io_free(data_obs)

  ENDIF


! 3.3  Check validity of input parameters
!----------------------------------------

  WRITE(slon_in, '(f12.5)') lon_in
  IF ( ABS(lon_in) > 360.0_wp ) THEN  ! Also catches lon_in = ropp_MDFV
    CALL message ( msg_fatal, 'Unset or invalid input longitude ' // TRIM(slon_in) )
  ENDIF

  WRITE(slat_in, '(f12.5)') lat_in
  IF ( ABS(lat_in) > 90.0_wp ) THEN  ! Also catches lat_in = ropp_MDFV
    CALL message ( msg_fatal, 'Unset or invalid input latitude ' // TRIM(slat_in) )
  ENDIF

  WRITE(sroc_in, '(f20.10)') roc_in
  IF ( (ABS(roc_in - 6.4e6_wp) > 0.2e6_wp) .AND. (roc_in > ropp_MDTV) ) THEN ! Usual valid range for roc
    CALL message ( msg_fatal, 'Invalid input radius of curvature ' // TRIM(sroc_in) )
  ENDIF

  WRITE(sund_in, '(f12.5)') und_in
  IF ( (ABS(und_in) > 150.0_wp) .AND. (und_in > ropp_MDTV) ) THEN ! Usual valid range for und
    CALL message ( msg_fatal, 'Invalid input undulation ' // TRIM(sund_in) )
  ENDIF

  WRITE(sazi_in, '(f12.5)') azi_in
  IF ( (ABS(azi_in) > 360.0_wp) .AND. (azi_in > ropp_MDTV) ) THEN
    CALL message ( msg_fatal, 'Invalid input azimuth ' // TRIM(sazi_in) )
  ENDIF

  DO i=1,3
    WRITE(scoc_in, '(f12.5)') coc_in(i)
    IF ( (ABS(coc_in(i)) > 5.0e4_wp) .AND. (ABS(coc_in(i) - ropp_MDFV) > 1.0e3_wp) ) THEN ! Usual valid range for coc
      CALL message ( msg_fatal, 'Invalid input centre of curvature component ' // TRIM(scoc_in) )
    ENDIF
  ENDDO

  IF (LEN(TRIM(datestring)) == 8) THEN
    READ(datestring, *) ii
    year_in = ii/10000
    mon_in  = (ii - year_in*10000)/100
    day_in  = ii - year_in*10000 - mon_in*100
  ELSE
    IF (datestring /= ' ') &
      CALL message ( msg_fatal, 'Input date ' // TRIM(datestring) // ' not valid.\n' )
  ENDIF

  IF (LEN(TRIM(timestring)) == 6) THEN
    READ(timestring, *) ii
    hour_in = ii/10000
    min_in  = (ii - hour_in*10000)/100
    sec_in  = ii - hour_in*10000 - min_in*100
  ELSE IF (LEN(TRIM(timestring)) == 4) THEN
    READ(timestring,*) ii
    hour_in   = ii/100
    min_in  = ii - hour_in*100
    sec_in  = 0
  ELSE
    IF (timestring /= ' ') &
      CALL message ( msg_fatal, 'Input time ' // TRIM(timestring) // &
                                ' not valid.\n' )
  ENDIF

  IF ( (file_in2 /= ' ') .AND. (datestring == ' ') .AND. (timestring == ' ') ) THEN
    CALL message ( msg_warn, 'Time interpolation implied by use of second GRIB file, ' // &
                             'but date and time not specified. \n' // &
                             'Will use fields from first file only, ' // &
                             'and no time interpolation will be carried out.')
    file_in2 = ' '
  ENDIF


! 3.4  Write out input parameters
!--------------------------------

  IF (msg_MODE  == VerboseMode) THEN
    CALL message ( msg_diag, 'file_in = ' // file_in )
    CALL message ( msg_diag, 'file_out_nml = ' // file_out_nml )
    CALL message ( msg_diag, 'file_out_nc = ' // file_out_nc )
    CALL message ( msg_diag, 'lon_in = ' // slon_in )
    CALL message ( msg_diag, 'lat_in = ' // slat_in )
    IF (datestring /= ' ')  CALL message ( msg_diag, 'datestring = ' // datestring )
    IF (timestring /= ' ')  CALL message ( msg_diag, 'timestring = ' // timestring )
    IF (file_in2 /= ' ')    CALL message ( msg_diag, 'file_in2 = ' // file_in2 )
    IF (file_inz /= ' ')    CALL message ( msg_diag, 'file_inz = ' // file_inz )
    IF (file_obs /= ' ')    CALL message ( msg_diag, 'file_obs = ' // file_obs )
    IF (roc_in > ropp_MDTV) CALL message ( msg_diag, 'roc_in = ' // sroc_in )
    IF (und_in > ropp_MDTV) CALL message ( msg_diag, 'und_in = ' // sund_in )
    IF (azi_in > ropp_MDTV) CALL message ( msg_diag, 'azi_in = ' // sazi_in )
    DO i=1,3
      WRITE(scoc_in, '(f12.5)') coc_in(i)
      IF (ABS(coc_in(i) - ropp_MDFV) > 1.0e3_wp) &
        CALL message ( msg_diag, 'coc_in = ' // scoc_in )
    ENDDO
  ENDIF


!-------------------------------------------------------------
! 4. Read GRIB file(s)
!-------------------------------------------------------------

! 4.1  Read variables from file_in
! --------------------------------

  CALL message ( msg_info, 'Reading file ' // TRIM(file_in) )

  CALL extract_prof_from_GRIB(file_in, lon_in, lat_in, grib1)

! 4.1.1 Check consistency with input data
! ---------------------------------------

  IF (grib1%lev2b%npoints > nlevs_max) &
    CALL message ( msg_fatal, 'Too many levels. ' // &
                              'Increase nlevs_max and recompile. \n' )

  IF ( (ABS(COS(dtor*lon_in)-COS(dtor*grib1%georef%lon)) >= 1.e-6_wp) .OR. &
       (ABS(SIN(dtor*lon_in)-SIN(dtor*grib1%georef%lon)) >= 1.e-6_wp) ) &
    CALL message ( msg_fatal, 'Output and input longitudes do not match. \n' )

  IF (ABS(lat_in-grib1%georef%lat) >= 1.0e-6_wp) &
    CALL message ( msg_fatal, 'Output and input latitudes do not match. \n' )

  CDT = (/ grib1%bg%year, grib1%bg%month, grib1%bg%day, 0, grib1%bg%hour, grib1%bg%minute, 0, 0/)
  CALL CalToJul( CDT,  jul_day1, 1 )

  IF ((datestring /= ' ') .AND. (timestring /= ' ')) THEN ! Compare input (if specified) and bgr times
    CDT = (/ year_in,  mon_in,  day_in,  0,  hour_in,  min_in,  sec_in,  msec_in/)
    CALL CalToJul( CDT,  jul_day_in, 1 )
    IF (ABS(jul_day1-jul_day_in) > 0.25d0) &
      CALL message ( msg_warn, 'Validity time of background data in ' // &
                                TRIM(file_in) // ' differs by ' // &
                               'more than 6 hours from input time. \n' )
  ENDIF

! 4.1.2 Populate the namelist
! ---------------------------

  year          = grib1%bg%year  ;  mon = grib1%bg%month   ;  day = grib1%bg%day
  hour          = grib1%bg%hour  ;  min = grib1%bg%minute  ;  sec = 0
  tlead         = grib1%bg%fcperiod

  lon           = grib1%georef%lon  ! = the input lon
  lat           = grib1%georef%lat  ! = the input lat

  Z0            = grib1%lev2c%geop_sfc
  P0            = grib1%lev2c%press_sfc*1.0e2_wp

  nlevs         = grib1%lev2b%npoints
  P(1:nlevs)    = grib1%lev2b%press*1.0e2_wp
  T(1:nlevs)    = grib1%lev2b%temp
  Q(1:nlevs)    = grib1%lev2b%shum

  Ak(1:nlevs+1) = grib1%lev2d%level_coeff_a*1.0e2_wp
  Bk(1:nlevs+1) = grib1%lev2d%level_coeff_b

  IF (roc_in < ropp_MDTV) THEN ! It's not specified at the command line or in the obs file
    roc = grib1%georef%roc
  ELSE
    roc = roc_in               ! Use the value specified by the user
  ENDIF

  IF (und_in < ropp_MDTV) THEN ! It's not specified at the command line or in the obs file
    und = grib1%georef%undulation
  ELSE
    und = und_in               ! Use the value specified by the user
  ENDIF

  IF (azi_in < ropp_MDTV) THEN ! It's not specified at the command line or in the obs file
    azi = grib1%georef%azimuth
  ELSE
    azi = azi_in               ! Use the value specified by the user
  ENDIF

  DO i=1,3
    IF (ABS(coc_in(i) - ropp_MDFV) < 1.0e3_wp) THEN ! It's not specified at the command line or in the obs file
      coc(i) = grib1%georef%r_coc(i)
    ELSE
      coc(i) = coc_in(i)       ! Use the value specified by the user
    ENDIF
  ENDDO

! 4.2  Read variables from file_in2
! ---------------------------------

  IF (file_in2 /= ' ') THEN

    CALL message ( msg_info, 'Reading file ' // TRIM(file_in2) )

    CALL extract_prof_from_GRIB(file_in2, lon_in, lat_in, grib2)

! 4.2.1 Check consistency with input data and that in file1
! ---------------------------------------------------------

    IF (grib1%lev2b%npoints > nlevs_max) &
      CALL message ( msg_fatal, 'Too many levels. ' // &
                                'Increase nlevs_max and recompile. \n' )

    IF ( (ABS(COS(dtor*lon_in)-COS(dtor*grib2%georef%lon)) >= 1.e-6_wp) .OR. &
         (ABS(SIN(dtor*lon_in)-SIN(dtor*grib2%georef%lon)) >= 1.e-6_wp) ) &
      CALL message ( msg_fatal, 'Output and input longitudes do not match. \n' )

    IF (ABS(lat_in-grib2%georef%lat) >= 1.0e-6_wp) &
      CALL message ( msg_fatal, 'Output and input latitudes do not match. \n' )

    IF (grib2%lev2b%npoints /= grib1%lev2b%npoints) &
      CALL message ( msg_fatal, 'Different number of levels in the two grib files. \n' )

    IF (ANY(ABS(grib2%lev2d%level_coeff_a-grib1%lev2d%level_coeff_a) >= 1.0e-6_wp)) &
      CALL message ( msg_fatal, 'Different Aks in the two grib files. \n' )

    IF (ANY(ABS(grib2%lev2d%level_coeff_b-grib1%lev2d%level_coeff_b) >= 1.0e-6_wp)) &
      CALL message ( msg_fatal, 'Different Bks in the two grib files. \n' )

    IF ( (ABS(COS(dtor*grib1%georef%lon)-COS(dtor*grib2%georef%lon)) >= 1.e-6_wp) .OR. &
         (ABS(SIN(dtor*grib1%georef%lon)-SIN(dtor*grib2%georef%lon)) >= 1.e-6_wp) ) &
      CALL message ( msg_fatal, 'Different longitudes in the two grib files. \n' )

    IF (ABS(grib2%georef%lat-grib1%georef%lat) >= 1.0e-6_wp) &
      CALL message ( msg_fatal, 'Different latitudes in the two grib files. \n' )

    IF (ABS(grib2%georef%undulation-grib1%georef%undulation) >= 1.0e-6_wp) &
      CALL message ( msg_fatal, 'Different undulations in the two grib files. \n' )

    IF (ABS(grib2%georef%roc-grib1%georef%roc) >= 1.0e-6_wp) &
      CALL message ( msg_fatal, 'Different radii of curvature in the two grib files. \n' )

    IF (ANY(ABS(grib2%georef%r_coc-grib1%georef%r_coc) >= 1.0e-6_wp)) &
      CALL message ( msg_fatal, 'Different centres of curvature in the two grib files. \n' )

    IF (ABS(grib2%georef%azimuth-grib1%georef%azimuth) >= 1.0e-6_wp) &
      CALL message ( msg_fatal, 'Different azimuths in the two grib files. \n' )

    IF ((grib1%lev2c%geop_sfc > ropp_MDTV) .AND. (grib2%lev2c%geop_sfc > ropp_MDTV) .AND. &
        (ABS(grib2%lev2c%geop_sfc-grib1%lev2c%geop_sfc) >= 1.0e-6_wp)) &
      CALL message ( msg_fatal, 'Different surface geopotentials in the two grib files. \n' )

! 4.2.2 Interpolate in time
! -------------------------

    CDT = (/ grib2%bg%year, grib2%bg%month, grib2%bg%day, 0, grib2%bg%hour, grib2%bg%minute, 0, 0/)
    CALL CalToJul( CDT,  jul_day2, 1 )

    IF (ABS(jul_day2 - jul_day1) <= 1.0d-6) THEN
      CALL message ( msg_warn, 'Same background times in the two grib files. \n'// &
                               'Using 50-50 mixture of fields from each. \n' )
      time_coeff(2) = 0.5_wp
    ELSE
      time_coeff(2) = REAL((jul_day_in - jul_day1) / (jul_day2 - jul_day1), wp)  ! time interpolation coefficient
    ENDIF
    time_coeff(1) = 1.0_wp - time_coeff(2)

    P0          = (time_coeff(1)*grib1%lev2c%press_sfc + time_coeff(2)*grib2%lev2c%press_sfc)*1.0e2_wp
    P(1:nlevs)  = (time_coeff(1)*grib1%lev2b%press     + time_coeff(2)*grib2%lev2b%press)*1.0e2_wp
    T(1:nlevs)  = time_coeff(1)*grib1%lev2b%temp       + time_coeff(2)*grib2%lev2b%temp
    Q(1:nlevs)  = time_coeff(1)*grib1%lev2b%shum       + time_coeff(2)*grib2%lev2b%shum

    IF ((grib1%lev2c%geop_sfc > ropp_MDTV) .AND. (grib2%lev2c%geop_sfc > ropp_MDTV)) THEN
      Z0        = time_coeff(1)*grib1%lev2c%geop_sfc + time_coeff(2)*grib2%lev2c%geop_sfc  ! (Should return Z0 anyway)
    ELSE
      IF (grib2%lev2c%geop_sfc > ropp_MDTV) Z0 = grib2%lev2c%geop_sfc  ! Take Z0 from the second file
    ENDIF

! In the case of time interpolation, set hour etc to be the _requested_ date/time.
! (This also avoids having to interpolate integers.)

    year = year_in  ;  mon = mon_in  ;  day = day_in
    hour = hour_in  ;  min = min_in  ;  sec = sec_in
    tlead = ropp_MDFV  ! Since the lead time is apparently not well defined for interpolated fields.

    CALL ropp_io_free ( grib2 )

  ENDIF  ! Interpolating in time


! 4.3  Read surface geopotential Z0 (only) from file_inz, if necessary
! --------------------------------------------------------------------

  IF (Z0 < ropp_MDTV) THEN

    IF (file_inz == ' ') &
      CALL message ( msg_fatal, 'Surface geopotential not present in input file(s), ' // &
                                'and an auxiliary file containing it is not specified.\n' )

    CALL message ( msg_info, 'Reading file ' // TRIM(file_inz) )

    CALL extract_prof_from_GRIB(file_inz, lon_in, lat_in, gribz, lgeop_sfc=.TRUE.)

    IF (gribz%lev2c%geop_sfc < ropp_MDTV) &
      CALL message ( msg_fatal, 'Surface geopotential not present in the auxiliary file ' // &
                                'which is required to contain it.\n' )

    Z0 = gribz%lev2c%geop_sfc ! Ignore the other gribz variables.

    CALL ropp_io_free ( gribz )

  ENDIF

  IF (Z0 < ropp_MDTV) &
    CALL message ( msg_fatal, 'Surface geopotential not present in any input file. \n' )


! 4.4 Generate Z at the profile location
! --------------------------------------

  CALL generate_GPH(Z0, P0, T(1:nlevs), Q(1:nlevs), Ak(1:nlevs+1), Bk(1:nlevs+1), Z(1:nlevs))


!-------------------------------------------------------------------------------
! 5. Define and write background profile namelist and ROPP netCDF file
!-------------------------------------------------------------------------------

! 5.1 Write out namelist for diagnostic purposes
! ----------------------------------------------

  IF (msg_MODE  == VerboseMode) THEN
    CALL message( msg_diag, 'Contents of bgr_profile namelist:' )
    WRITE (*, NML=bgr_profile)  ! For diagnostic purposes
  ENDIF

! 5.2 Write out namelist to ascii file
! ------------------------------------

  CALL message ( msg_info, 'Writing file ' // TRIM(ADJUSTL(file_out_nml)) )

  OPEN (UNIT=namelist_lun, FILE=TRIM(ADJUSTL(file_out_nml)), STATUS='replace', ACTION='write', IOSTAT=iostatus)

  IF (iostatus > 0) &
    CALL message ( msg_fatal, 'I/O error when opening namelist file' )

  WRITE (namelist_lun, NML=bgr_profile, IOSTAT=iostatus) 

  IF (iostatus > 0) &
    CALL message ( msg_fatal, 'I/O error while writing namelist file' )

  CLOSE (UNIT=namelist_lun)

! 5.3 Write out bg structure to ROPP file
! ---------------------------------------

! 5.3.1 Transfer data from obs file to ROprof structure grib1
! -----------------------------------------------------------

    IF (file_obs /= ' ') THEN
      grib1%leo_id = TRIM(sleo_id_in)
      grib1%gns_id = TRIM(sgns_id_in)
      grib1%stn_id = TRIM(sstn_id_in)
      grib1%occ_id = TRIM(socc_id_in)
      grib1%occ_id(1:2) = 'BG'
    ENDIF

    grib1%dtocc%year   = year_in
    grib1%dtocc%month  = mon_in
    grib1%dtocc%day    = day_in
    grib1%dtocc%hour   = hour_in
    grib1%dtocc%minute = min_in
    grib1%dtocc%second = sec_in
    grib1%dtocc%msec   = msec_in

    grib1%georef%lon = ATAN2(SIN(dtor*grib1%georef%lon), COS(dtor*grib1%georef%lon)) / dtor  ! Between -180 and 180
    IF ( ABS(grib1%georef%lon) < 1.0E-6_wp ) grib1%georef%lon = 0.0_wp

! 5.3.2 Copy namelist variables that might have changed
!       by time interpolation, to ROprof structure grib1
! ------------------------------------------------------

    grib1%lev2c%geop_sfc  = Z0
    grib1%lev2c%press_sfc = P0*1.0e-2_wp

    grib1%lev2b%geop      = Z(1:nlevs)
    grib1%lev2b%press     = P(1:nlevs)*1.0e-2_wp
    grib1%lev2b%temp      = T(1:nlevs)
    grib1%lev2b%shum      = Q(1:nlevs)

    grib1%bg%year         = year
    grib1%bg%month        = mon
    grib1%bg%day          = day
    grib1%bg%hour         = hour
    grib1%bg%minute       = min
    grib1%bg%fcperiod     = tlead

    grib1%PCD = 0
    grib1%PCD = IBSET(grib1%PCD, PCD_occultation)  ! set PCD bit 15 to 1 (background data)

    grib1%software_version = ropp_io_version()
    IF (INDEX(grib1%software_version, ' ') > 1) THEN
      READ(grib1%software_version(2:INDEX(grib1%software_version, ' ')-1), *) swver
      IF (INT(swver) < 10) THEN
        grib1%software_version = 'V0' // grib1%software_version(2:INDEX(grib1%software_version, ' ')-1)
      ELSE
        grib1%software_version = 'V'  // grib1%software_version(2:INDEX(grib1%software_version, ' ')-1)
      ENDIF
      IF (LEN_TRIM(grib1%software_version) < 8) THEN
        DO i=LEN_TRIM(grib1%software_version),8
          grib1%software_version(i:i) = '0'
        ENDDO
      ENDIF
    ENDIF

    CALL ropp_io_write(grib1, file_out_nc, IERR=iostatus, ranchk=ranchk)

    IF ( iostatus > 0 ) THEN
      CALL message ( msg_fatal, "Error while writing output ROPP file" )
    END IF

! 5.4 Tidy up - deallocate structures & free memory
!--------------------------------------------------

    CALL ropp_io_free ( grib1 )
    CALL message ( msg_noin, ' ' )
    CALL EXIT(msg_exit_status)

!=============================================================================

CONTAINS


!-------------------------------------------------------------------------------
! 6. Extract profile from GRIB file
!-------------------------------------------------------------------------------

  SUBROUTINE extract_prof_from_GRIB(gribfile, lon_in, lat_in, grib, lgeop_sfc)

    USE grib_api
    USE EarthMod, ONLY: Datum_HMSL
    USE ropp_io,  ONLY: ropp_io_init

    IMPLICIT NONE

! Input variables
    CHARACTER(LEN=256), INTENT(IN)          :: gribfile  ! Input grib file
    REAL(wp), INTENT(IN)                    :: lon_in    ! Input longitude (deg)
    REAL(wp), INTENT(IN)                    :: lat_in    ! Input latitude (deg)
    LOGICAL, OPTIONAL                       :: lgeop_sfc ! Stop looking after geop_sfc has been found

! Output variables
    TYPE(ROprof), INTENT(OUT)               :: grib      ! ROprof structure holding bg profile

! Local variables
! GRIB-specific
    INTEGER                                 :: nlon, nlat, nlev, numberOfPoints
    REAL(wp)                                :: lon0, dlon, lon1, lon2
    REAL(wp)                                :: lat0, dlat, lat1, lat2
    INTEGER                                 :: ilon, ineg
    INTEGER                                 :: ilat, jpos, jcons, jalt
    INTEGER                                 :: date, time, tlead
    INTEGER                                 :: gribfile_lun, iret, igrib, nmess, ilev
    INTEGER                                 :: centre, gridsize, pvsize
    INTEGER                                 :: pvpresent, grib_ed, grib_procid
    CHARACTER(LEN=4)                        :: slev
    CHARACTER(LEN=8)                        :: shortname
    CHARACTER(LEN=256)                      :: gridtype, sprocid
    REAL(wp), DIMENSION(:),  ALLOCATABLE    :: pv       ! Ak and Bk.
    REAL(wp), DIMENSION(4)                  :: lats, lons, distances, values
    INTEGER, DIMENSION(4)                   :: indices
! Other
    CHARACTER(LEN=256)                      :: sret, routine
    REAL(wp), DIMENSION(4)                  :: Z0_grib, P0_grib
    REAL(wp), DIMENSION(:, :), ALLOCATABLE  :: T_grib, Q_grib
    REAL(wp), DIMENSION(2)                  :: lon_coeff    ! Interpolation coeff
    REAL(wp), DIMENSION(2)                  :: lat_coeff    ! Interpolation coeff
    REAL(wp), DIMENSION(4)                  :: lonlat_coeff ! Interpolation coeff
    LOGICAL                                 :: l_geop_sfc
! Constants
    REAL(wp), PARAMETER                     :: g_wmo=9.80665_wp   ! Standard gravity (m/s2)

! 6.1 Initialise
! --------------

  gribfile_lun = Get_IO_Unit()

  l_geop_sfc = .FALSE.
  IF ( PRESENT(lgeop_sfc) ) l_geop_sfc = lgeop_sfc

! 6.1.1 Set routine for error messages
! ------------------------------------

    CALL message_get_routine(routine)
    CALL message_set_routine('extract_prof_from_GRIB')

! 6.1.2 Open file for reading metadata from first message;
!       check that it is readable as a GRIB file.
! --------------------------------------------------------

    CALL grib_open_file(gribfile_lun, TRIM(gribfile), 'R', iret)

    IF (iret /= GRIB_SUCCESS) THEN
      IF (msg_MODE  == VerboseMode) THEN
        CALL grib_get_error_string (iret, sret)
        CALL message( msg_diag, 'GRIB error message: ' // TRIM(sret) )
      ENDIF
      CALL message ( msg_fatal, 'Error opening ' // TRIM(gribfile) // &
                                ' - check file existence and format' )
    ENDIF

    CALL grib_new_from_file(gribfile_lun, igrib, iret)

    IF (iret /= GRIB_SUCCESS) THEN
      IF (msg_MODE  == VerboseMode) THEN
        CALL grib_get_error_string (iret, sret)
        CALL message( msg_diag, 'GRIB error message: ' // TRIM(sret) )
      ENDIF
      CALL message ( msg_fatal, 'Error reading ' // TRIM(gribfile) // &
                                ' - check file format' )
    ENDIF

! 6.1.3  Read header info from first message in GRIB file
! -------------------------------------------------------

    CALL grib_get(igrib, 'editionNumber', grib_ed)
    IF (grib_ed == 1) THEN
      CALL message ( msg_diag, 'GRIB Edition 1' )
    ELSE IF (grib_ed == 2) THEN
      CALL message ( msg_diag, 'GRIB Edition 2' )
    ELSE
      CALL message ( msg_fatal, 'Invalid GRIB edition (not 1 or 2)' )
    ENDIF

    CALL grib_get(igrib, 'generatingProcessIdentifier', grib_procid)
    WRITE (sprocid, '(I4)') grib_procid
    CALL message ( msg_diag, 'ECMWF generatingProcessIdentifier = ' // sprocid )

    CALL grib_get(igrib, 'centre', centre)
    IF (centre /= 98) &
      CALL message ( msg_fatal, 'Invalid processing centre (not 98 = ECMWF)' )

    CALL grib_get(igrib, 'gridType', gridtype)
    IF ( (TRIM(ADJUSTL(gridtype)) /= 'regular_ll') .AND. &
         (TRIM(ADJUSTL(gridtype)) /= 'regular_gg') ) &
      CALL message ( msg_fatal, 'Not a regular lat-lon or Gaussian grid' )

    CALL grib_get(igrib, 'numberOfPoints', gridsize) ! Alias 'numberOfDataPoints'

    CALL grib_get(igrib, 'Ni', nlon)
    CALL grib_get(igrib, 'Nj', nlat)

    IF ( (nlat < 1) .OR. (nlon < 1) .OR. (gridsize /= nlat*nlon) ) &
      CALL message ( msg_fatal, &
        'Error: invalid or inconsistent nlat, nlon or gridsize' )

    CALL grib_get(igrib, 'PVPresent', pvpresent)
    IF ( pvpresent == 1 ) THEN
      CALL grib_get_size(igrib, 'pv', pvsize)
    ELSE
      CALL message ( msg_fatal, 'pv not in file ... cannot calculate number of levels' )
    ENDIF

    nlev = (pvsize - 2) / 2       ! Number of full model levels

    IF ( (nlev /= 60) .AND. (nlev /= 91) .AND. (nlev /= 137) ) THEN
      WRITE ( slev, '(I4)' ) nlev
      CALL message ( msg_warn, 'A non-standard number of levels' // &
                               ' (' // slev // ') in the GRIB file' )
    ENDIF

! 6.1.4  Initialise ROprof structure that will hold the GRIB data
! ---------------------------------------------------------------

    CALL ropp_io_init (grib, 0, 0, 0, nlev, 1, nlev+1)

! 6.1.5  Extract verification (or validity) date/time info from GRIB file.
!        (This is how bg%hour is defined in the ROPP I/O UG.)
! -----------------------------------------------------------

    CALL grib_get(igrib, 'validityDate', date)
    CALL grib_get(igrib, 'validityTime', time)
    CALL grib_get(igrib, 'stepRange', tlead)

    grib%bg%year     = date / 10000
    grib%bg%month    = (date - grib%bg%year*10000) / 100
    grib%bg%day      = date - grib%bg%year*10000 - grib%bg%month*100
    grib%bg%hour     = time / 100
    grib%bg%minute   = time - grib%bg%hour*100
    grib%bg%fcperiod = tlead

    grib%bg%source        = 'ECMWF'
    grib%lev2d%level_type = 'ECMWF'

! 6.1.6 Read pv and extract Ak and Bk from it, if possible.
!       (The ECMWF pv values changed with the introduction
!       of GRIB2 output at Cy37r2 on 18 May 2011.  In GRIB1 files
!       they are held less accurately.  The same inaccuracy is inherited
!       by nominally GRIB2 files which are converted from GRIB1 format.)
! ----------------------------------------------------------------

    ALLOCATE (pv(pvsize))

    CALL grib_get(igrib, 'pv', pv)

    grib%Lev2d%level_coeff_a = pv(nlev+1:1:-1)*1.0e-2_wp  ! Pa to hPa for ROPP
    grib%Lev2d%level_coeff_b = pv(2*(nlev+1):nlev+2:-1)

    DEALLOCATE (pv)

    IF ((grib_ed == 1) .OR. ((grib_ed == 2) .AND. (grib_procid < 140))) THEN

      CALL message( msg_info, 'Ak and Bk in GRIB Ed 1 file are inaccurate ... ' // &
                              'extracting from known standard values (if possible)' )

      SELECT CASE (nlev)

        CASE (60)
          grib%Lev2d%level_coeff_a = Ak60*1.0e-2_wp
          grib%Lev2d%level_coeff_b = Bk60

        CASE (91)
          grib%Lev2d%level_coeff_a = Ak91*1.0e-2_wp
          grib%Lev2d%level_coeff_b = Bk91

        CASE (137)
          grib%Lev2d%level_coeff_a = Ak137*1.0e-2_wp
          grib%Lev2d%level_coeff_b = Bk137

        CASE DEFAULT

          CALL message( msg_warn, 'Non-standard number of levels ... ' // &
                                  'using valudes in GRIB file anyway' )

      END SELECT

    ENDIF

    grib%Lev2d%level_coeff_a(nlev+1) = 1.0e-34_wp  ! always

! 6.1.7 Find indices of the nearest points to (lon_in, lat_in)
! ------------------------------------------------------------

    CALL grib_find_nearest(igrib, .FALSE., lat_in, lon_in, &
                           lats, lons, values, distances, indices) ! .FALSE. to find nearest point whether or not it's land.

    IF ( ANY(indices < 0)                                            .OR. &
         ANY(indices > (gridsize-1))                                 .OR. &
         ANY((/ indices(2), indices(3), indices(4) /) == indices(1)) .OR. &
         ANY((/ indices(3), indices(4), indices(1) /) == indices(2)) .OR. &
         ANY((/ indices(4), indices(1), indices(2) /) == indices(3)) .OR. &
         ANY((/ indices(1), indices(2), indices(3) /) == indices(4)) ) THEN
      CALL message( msg_warn, 'grib_find_nearest has returned invalid indices ... ' // &
                              'will try to generate them' )

      IF ( TRIM(ADJUSTL(gridtype)) /= 'regular_ll' ) &
        CALL message( msg_fatal, 'Cannot generate indices for gridtype ' // gridtype )

      CALL grib_get(igrib, 'longitudeOfFirstGridPointInDegrees', lon0)
      CALL grib_get(igrib, 'iDirectionIncrementInDegrees', dlon)
      CALL grib_get(igrib, 'iScansNegatively', ineg)

      CALL grib_get(igrib, 'latitudeOfFirstGridPointInDegrees', lat0)
      CALL grib_get(igrib, 'jDirectionIncrementInDegrees', dlat)
      CALL grib_get(igrib, 'jScansPositively', jpos)

      CALL grib_get(igrib, 'jPointsAreConsecutive', jcons)
      IF (jcons == 1) &
        CALL message (msg_fatal, 'Cannot (yet) handle grids with consecutive j points' )
! Fixing this should be basically a matter of swapping over the two loops below.

      CALL grib_get(igrib, 'alternativeRowScanning', jalt)
      IF (jalt == 1) &
        CALL message (msg_fatal, 'Cannot (yet) handle grids with alternative row scanning' )
! This would be harder to fix.

      indices = (/ 0, 1, nlon, nlon+1 /)  ! Initialise to the corners of the first gridblock.  NB: 0-based indexing.

      lat1 = lat0  ;  lat2 = lat1 + dlat*(2*jpos - 1)
      DO ilat=1,nlat
        IF ( (lat2-lat_in)*(lat1-lat_in) <= 0 ) EXIT
        indices = indices + (/ nlon, nlon, nlon, nlon /)
        lat1 = lat2  ;  lat2 = lat1 + dlat*(2*jpos - 1)
      ENDDO
      IF (ilat == nlat) indices(3:4) = indices(3:4) - 2*nlon  ! Bounce backwards

      lon1 = lon0  ;  lon2 = lon1 + dlon*(1 - 2*ineg)
      DO ilon=1,nlon
        IF (((SIN(dtor*lon2)-SIN(dtor*lon_in))*(SIN(dtor*lon1)-SIN(dtor*lon_in)) <= 0) .AND. &
            ((COS(dtor*lon2)-COS(dtor*lon_in))*(COS(dtor*lon1)-COS(dtor*lon_in)) <= 0)) EXIT
        indices = indices + (/ 1, 1, 1, 1 /)
        lon1 = lon2  ;  lon2 = lon1 + dlon*(1 - 2*ineg)
      ENDDO
      IF (ilon == nlon) indices(2:4:2) = indices(2:4:2) - nlon  ! Cycle back

      CALL grib_get_element(igrib, 'longitudes', indices, lons)
      CALL grib_get_element(igrib, 'latitudes', indices, lats)

    ENDIF

    CALL grib_close_file(gribfile_lun)

! 6.1.8 Calculate interpolation coefficients
! ------------------------------------------
! Currently set to bilinear, but we could use the distances above in 
! an 'inverse distance' method (eg as in SCRIP).  This would allow the code
! to be used with non-uniform grids, such as reduced gaussian grids, 
! and might be better near the pole, where bilinear could have problems.

! Prototype code for inverse-distance interpolation method:
!   lonlat_coeff = (1.0_wp/(distances + TINY(1.0_wp))) / &
!              SUM((1.0_wp/(distances + TINY(1.0_wp))))

! Prototype code for nearest neighbour interpolation method:
! (Doesn't cope with degenerate cases, but could be made to do so, which
!  CALL grib_find_nearest with scalar values, distance, indices etc, couldn't.)
!   lonlat_coeff(:) = 0.0_wp
!   lonlat_coeff(SUM(MINLOC(distances))) = 1.0_wp

! Bilinear interpolation in (lon, lat).
! This formula should be safe against lons crossing the cut, and if lons(2) < lons(1)
    lon_coeff(2) = ATAN2(SIN(dtor*(lon_in -lons(1))), COS(dtor*(lon_in -lons(1)))) / &
                   ATAN2(SIN(dtor*(lons(2)-lons(1))), COS(dtor*(lons(2)-lons(1))))
    lon_coeff(1) = 1.0_wp - lon_coeff(2)

! This formula should be safe against lats(3) < lats(1)
    lat_coeff(2) = ATAN2(SIN(dtor*(lat_in -lats(1))), COS(dtor*(lat_in -lats(1)))) / &
                   ATAN2(SIN(dtor*(lats(3)-lats(1))), COS(dtor*(lats(3)-lats(1))))
    lat_coeff(1) = 1.0_wp - lat_coeff(2)

! Calculate bilinear interpolation coefficients
    lonlat_coeff(1) = lon_coeff(1)*lat_coeff(1)
    lonlat_coeff(2) = lon_coeff(2)*lat_coeff(1)
    lonlat_coeff(3) = lon_coeff(1)*lat_coeff(2)
    lonlat_coeff(4) = lon_coeff(2)*lat_coeff(2)

! 6.2 Reopen file and read all required grib messages
! ---------------------------------------------------

    Z0_grib = ropp_MDFV
    P0_grib = ropp_MDFV
    ALLOCATE(T_grib(4, nlev), Q_grib(4, nlev))
    T_grib = ropp_MDFV
    Q_grib = ropp_MDFV

    nmess = 0

    CALL grib_open_file(gribfile_lun, TRIM(gribfile), 'R', iret)

    CALL grib_new_from_file(gribfile_lun, igrib, iret)

    DO WHILE (iret /= GRIB_END_OF_FILE)  ! Loop over all fields in grib file

      CALL grib_get(igrib, 'level', ilev)

      ilev = nlev - ilev + 1  ! Reset level index

      CALL grib_get(igrib, 'shortName', shortname)

      IF (TRIM(shortName) == 't') THEN  ! Temperature in K

        CALL grib_get_element(igrib, 'values', indices, values)
        T_grib(:, ilev) = values
        nmess = nmess + 1

      ELSEIF (TRIM(shortName) == 'q') THEN  ! Specific humidity in kg/kg

        CALL grib_get_element(igrib, 'values', indices, values)
        Q_grib(:, ilev) = values*1000.0_wp  ! Convert to specific humidity in g/kg
        nmess = nmess + 1

      ELSEIF ((TRIM(shortName) == 'z') .AND. (ilev == nlev)) THEN  ! Surface geopotential in m2/s2

        CALL grib_get_element(igrib, 'values', indices, values)
        Z0_grib(:) = values / g_wmo  ! Convert to height in m
        nmess = nmess + 1

      ELSEIF ((TRIM(shortName) == 'lnsp') .AND. (ilev == nlev)) THEN  ! Log(surface pressure in Pa)

        CALL grib_get_element(igrib, 'values', indices, values)
        P0_grib(:) = values  ! Interpolate log(surface pressure in Pa), then exponentiate.
        nmess = nmess + 1

      ENDIF

      CALL grib_release(igrib)

      CALL grib_new_from_file(gribfile_lun, igrib, iret)

      IF (l_geop_sfc .AND. ALL(Z0_grib > ropp_MDTV)) iret = GRIB_END_OF_FILE ! Don't look any further

    ENDDO  ! Loop over fields in grib file

    CALL grib_close_file(gribfile_lun)

    IF ((.NOT. l_geop_sfc) .AND. (nmess < (2*nlev+1))) &
      CALL message( msg_fatal, &
        'Not enough T, q, lnsp or Z0 messages in GRIB file' )

    IF ((.NOT. l_geop_sfc) .AND. (nmess > (2*nlev+2))) &
      CALL message( msg_warn, &
        'More than enough T, q, lnsp or Z0 messages in GRIB file ' // &
        '... beware of inconsistencies' )

! 6.3 Interpolate to tangent point location
! -----------------------------------------

    grib%georef%lon = ATAN2(SIN(dtor*lon_in), COS(dtor*lon_in)) / dtor  ! Between -180 and 180
    IF ( ABS(grib%georef%lon) < 1.0E-6_wp ) grib%georef%lon = 0.0_wp

    grib%georef%lat = lat_in

    grib%lev2c%geop_sfc = SUM(lonlat_coeff*Z0_grib)

    grib%lev2c%press_sfc = EXP(SUM(lonlat_coeff*P0_grib))*1.0e-2_wp ! Pa to hPa for ROPP

    DO ilev=1,nlev
      grib%lev2b%temp(ilev)  = SUM(lonlat_coeff(:)*T_grib(:, ilev))
      grib%lev2b%shum(ilev)  = SUM(lonlat_coeff(:)*Q_grib(:, ilev))
      grib%lev2b%press(ilev) = 0.5_wp*( &
        SUM(grib%Lev2d%level_coeff_a(ilev:ilev+1)) + &
        SUM(grib%Lev2d%level_coeff_b(ilev:ilev+1))*grib%lev2c%press_sfc )
    ENDDO

    DEALLOCATE(T_grib, Q_grib)

    grib%lev2c%press_sfc_qual = 100.0_wp
    grib%lev2b%meteo_qual     = 100.0_wp

! 6.4 Generate undulation at the tangent point
! --------------------------------------------
!     This routine requires the user to have exported the environment variables:
! (1) GEOPOT_COEF, which is the full pathname of the file containing the egm96 geoid coefficients
!     (eg ...ropp_pp/data/egm96.dat); and
! (2) GEOPOT_CORR, which is the full pathname of the file containing the correction coefficients
!     (eg ...ropp_pp/data/corrcoef.dat)
! -----------------------------------------------------------------------------------------------

    CALL datum_HMSL('WGS84', (/lat_in, lon_in, 0.0_wp/), grib%georef%undulation)

    IF (grib%georef%undulation > ropp_MDTV) &
      grib%georef%undulation = -grib%georef%undulation ! Because datum_HMSL actually returns height of ellipsoid (0 here) wrt MSL.

! 6.5 Clean up
! ------------

    CALL message_set_routine(routine)

  END SUBROUTINE extract_prof_from_GRIB


!-------------------------------------------------------------------------------
! 7. Generate geopotential height at the interpolated point by integrating density
!-------------------------------------------------------------------------------

  SUBROUTINE generate_GPH(Z0, P0, T, Q, Ak, Bk, Z)

    IMPLICIT NONE

! Input variables
    REAL(wp),               INTENT(IN)    :: Z0, P0
    REAL(wp), DIMENSION(:), INTENT(IN)    :: T, Q
    REAL(wp), DIMENSION(:), INTENT(IN)    :: Ak, Bk

! Output variables
    REAL(wp), DIMENSION(:), INTENT(OUT)   :: Z

! Local
    REAL(wp), DIMENSION(:), ALLOCATABLE   :: Tvirt, Phlv, deltaP
    REAL(wp), DIMENSION(:), ALLOCATABLE   :: lnP, Zhlv, dZ, alpha
    INTEGER                               :: ihlev, Nz
    CHARACTER(LEN=256)                    :: routine

! Constants
    REAL(wp), PARAMETER                   :: g_wmo=9.80665_wp   ! Standard gravity (m/s2)
    REAL(wp), PARAMETER                   :: R_dry=287.0597_wp  ! Dry gas constant (K/kg/K)

! 7.1 Set subroutine name for messages
! ------------------------------------

    CALL message_get_routine(routine)
    CALL message_set_routine('generate_GPH')

! 7.2 Initialise
! --------------

    Nz = SIZE(ak) - 1
    ALLOCATE(Tvirt(Nz), Phlv(Nz+1), deltaP(Nz), lnP(Nz), Zhlv(Nz+1), dZ(Nz), alpha(Nz))

! 7.3 Integrate density to get GPH
! --------------------------------

    !--- Virtual temperatures
    Tvirt = (1.0_wp + 0.61_wp * Q/1000.0_wp) * T

    !--- Pressure differences
    Phlv = ak + bk*P0
    deltaP = Phlv(1:Nz) - Phlv(2:Nz+1)

    !--- Log of pressure ratio
    lnP = LOG(Phlv(1:Nz)/Phlv(2:Nz+1))

    !--- Interpolation coefficients
    alpha = 1.0_wp - Phlv(2:Nz+1)/deltaP * lnP
    alpha(Nz) = REAL(LOG(2.0d0), wp)

    !--- Function to be integrated
    dZ = R_dry * Tvirt * lnP / g_wmo

    !--- Calculate geopotential height integral
    Zhlv(1) = 0.0_wp
    DO ihlev = 2,Nz+1
      Zhlv(ihlev) = SUM(dZ(1:ihlev-1))
    ENDDO
    IF (Z0 > ropp_MDTV) Zhlv(:) = Zhlv(:) + Z0

    !--- Interpolate onto full levels
    Z = Zhlv(1:Nz) + alpha * R_dry * Tvirt / g_wmo

    DEALLOCATE(Tvirt, Phlv, deltaP, lnP, Zhlv, dZ, alpha)

! 7.4 Clean up
! ------------

    CALL message_set_routine(routine)

  END SUBROUTINE generate_GPH


!-------------------------------------------------------------------------------
! 8. Usage (help) information
!-------------------------------------------------------------------------------

  SUBROUTINE Usage()
    PRINT*, 'Purpose:'
    PRINT*, '  Extract a background fields profile, in fortran namelist format,'
    PRINT*, '  from a GRIB (edition 1 or 2) file containing background model data.'
    PRINT*, 'Usage:'
    PRINT*, '  > grib2bgrasc grib_file -lat <lat> -lon <lon>'
    PRINT*, '                [-date <date> ] [-time <time>]'
    PRINT*, '                [-g <grib_file2>]'
    PRINT*, '                [-z <orog_file>]'
    PRINT*, '                [-y <obs_file>]'
    PRINT*, '                [-o <namelist_file>]'
    PRINT*, '                [-d] [-h] [-v]'
    PRINT*, '  where:'
    PRINT*, '    grib_file is a GRIB file holding {p*, phi*, T, q};'
    PRINT*, '    -lat <lat> is the latitude of the required profile (deg);'
    PRINT*, '    -lon <lon> is the longitude of the required profile (deg).'
    PRINT*, 'Options:'
    PRINT*, '  -date <date> is the date required (yyyymmdd)'
    PRINT*, '  -time <time> is the time required (HHMM or HHMMSS)'
    PRINT*, '  -g <grib_file2> is a second GRIB file, if time interpolation is desired'
    PRINT*, '  -z <grib_filez> is an additional GRIB file containing surface geopotential,'
    PRINT*, '                  if it is not in either grib_file or grib_file2'
    PRINT*, '  -o <namelist_file> is the output namelist file name'
    PRINT*, '  -y <obs_file>   is a ROPP_format observation file holding the'
    PRINT*, '                  desired date, time, lat, lon, roc, undulation and'
    PRINT*, '                  satellite and occultation IDs for the background'
    PRINT*, '                  (Command line values take precendence.)'
    PRINT*, '  -x <netCDF_file> is the name of the output ROPP netCDF file'
    PRINT*, '  -r <radius_of_curvature> is the radius of curvature (in m)'
    PRINT*, '  -u <undulation> is the undulation (in m)'
    PRINT*, '  -a <azimuth> is the azimuth (in deg)'
    PRINT*, '  -c <centre of curvature> is the x,y,z components of the '
    PRINT*, '     centre of curvature in the ECF frame (in m)'
    PRINT*, '  -d prints out some additional diagnostics to stdout'
    PRINT*, '  -h this help'
    PRINT*, '  -v version information'
    PRINT*, 'Defaults:'
    PRINT*, '  Input  file name : required'
    PRINT*, '  Output file names : $PWD/<input_file>.nml and $PWD/<input_file>.nc'
    PRINT*, 'Note that the environment variables GEOPOT_COEF and GEOPOT_CORR must be'
    PRINT*, '  set appropriately if the user wants the undulation to be calculated.'
    PRINT*, 'See grib2bgrasc(1) for details.'
    PRINT*, ''
  END SUBROUTINE Usage


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

  SUBROUTINE version_info()
    USE ropp_io, ONLY : ropp_io_version
    CHARACTER (LEN=40) :: version
    version = ropp_io_version()
    PRINT*, 'grib2bgrasc -  GRIB to background ASCII decoder [ECMWF library]'
    PRINT*, ''
    PRINT*, 'This program is part of ROPP (IO) Release ' // TRIM(version)
    PRINT*, ''
  END SUBROUTINE version_info


END PROGRAM grib2bgrasc