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

PROGRAM ropp_fm_bg2ro_1d

!****p* Programs/ropp_fm_bg2ro_1d *
!
! NAME
!   ropp_fm_bg2ro_1d
!
! SYNOPSIS
!   Calculate radio occultation pseudo observation from background model
!   data using 1d forward models
!
!   > ropp_fm_bg2ro_1d <infile(s)> -o <outfile>
!                     [-l] <levels_file> [-c] <config_file>
!                     [-f] [-use_logp] [-use_logq]
!                     [-comp] [-check_qsat] [-nocheck_qmin] [-new_op] [-direct_ion] [-best]
!                     [-refrac_only] [-zmin <geop_min>] [-zmax <geop_max>] [-nz <n_geop>]
!                     [-bangle_only] [-247L] [-ih_from_geop] [-ihmin <ih_min>] [-ihmax <ih_max>] [-nih <n_ih>]
!                     [-d] [-h] [-v]
!
! ARGUMENTS
!   <infile(s)>         one (or more) input file names
!
! OPTIONS
!   -o <outfile>        name of output file (default: bg2ro.nc)
!   -l <levels_file>    name of file containing a set of altitudes and
!                       impact heights. This will typically be an
!                       observation file. Only works for single file.
!   -c <config_file>    name of file containing the configuration namelist
!   -f                  forward model only, no gradient calculation.
!   -use_logp           use log(pressure) for forward model
!   -use_logq           use log(spec humidity) for forward model
!   -comp               include non ideal gas compressibility
!   -check_qsat         check against supersaturation
!   -nocheck_qmin       do not check against superdryness
!   -new_op             use alternative refrac and bangle interpolation
!   -direct_ion         forward model L1 and L2 directly, using model ionosphere
!   -best               currently advised 'best' options
!   -refrac_only        do not generate bending angles
!   -zmin <geop_min>    minimum refractivity geopotential (gpm)
!   -zmax <geop_max>    maximum refractivity geopotential (gpm)
!   -nz <n_geop>        number of uniformly spaced refractivity geopotentials
!   -bangle_only        do not generate refractivities
!   -247L               output bangle on 'standard' 247 levels
!   -ih_from_geop       infer impact heights from refractivity geopotentials
!   -ihmin <ih_min>     minimum impact height (m)
!   -ihmax <ih_max>     maximum impact height (m)
!   -nih <n_ih>         number of uniformly spaced impact_heights
!   -d                  output additional diagnostics
!   -h                  help
!   -v                  version information
!
! DESCRIPTION
!   This program reads model data on model levels from the input data files
!   and calculates vertical profiles of bending angle and refractivity using
!   the 1d forward operators. The result is written to an ROPP formatted
!   output file.
!
! NOTES
!   If the input file is a multifile, or more than one input files are
!   specified, the output file is a multifile.
!
!   Already existing output files will be overwritten.
!
! EXAMPLES
! 1) To calculate bending angle and refractivity from one of the example
!    (single-) files in the data directory:
!
!      > ropp_fm_bg2ro_1d ../data/bgr20090401_000329_M02_2030337800_N0007_XXXX.nc
!
! 2) To calculate bending angle and refractivity profiles from all singlefiles
!    in the data directory:
!
!      > ropp_fm_bg2ro_1d ../data/bgr20090401*_N0007_XXXX.nc -o eg_02.nc
!
!    Note that the resulting eg_02.nc file contains forward modelled data from
!    all example profiles.
!
! 3) To calculate forward modelled bending angle and refractivity profiles from
!    all profiles contained in the multifile bgr20090401_multi.nc:
!
!      > ropp_fm_bg2ro_1d ../data/bgr20090401_multi.nc -o eg_03.nc
!
!    Since the ecmwf_multi_* file was generated by concatenating the other
!    files in the data directory, eg_02.nc and eg_03.nc should be identical
!    apart from the file names.
!
! 4) To calculate forward modelled L1 and L2 bending angles, assuming
!    a model ionosphere, from one of the example (single-) files
!    in the data directory:
!
!      > ropp_fm_bg2ro_1d --direct_ion ../data/bgr20090401_000329_M02_2030337800_N0007_YYYY.nc
!
!     The model ionosphere comprises a set of Chapman layers, each of
!     whose {Ne_max, H_peak, H_width} parameters are contained in the input file.
!     (The L2ctype and state1dFM structures have been suitably extended.)
!
! SEE ALSO
!   ropp_fm_bangle_1d
!   ropp_fm_refrac_1d
!   ropp_fm_refrac_1d_new
!
! 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.
!
!****

!-------------------------------------------------------------------------------
! 1. Declarations
!-------------------------------------------------------------------------------

  USE typesizes, ONLY: wp => EightByteReal
  USE ropp_utils
  USE ropp_io
  USE ropp_io_types, ONLY: ROprof
  USE ropp_fm
  USE ropp_fm_types
  USE ropp_fm_copy
  USE ropp_fm_iono
  USE ropp_fm_levels

  IMPLICIT NONE

  TYPE(ROprof)                                     :: ro_data, l_data
  TYPE(State1dFM)                                  :: state
  TYPE(Obs1dRefrac)                                :: obs_refrac
  TYPE(Obs1dBangle)                                :: obs_bangle

  REAL(wp), PARAMETER                              :: zmin_default=200.0_wp
  REAL(wp), PARAMETER                              :: zmax_default=60000.0_wp
  INTEGER, PARAMETER                               :: nz_default=300

  REAL(wp), PARAMETER                              :: amin_default=2000.0_wp
  REAL(wp), PARAMETER                              :: amax_default=60000.0_wp
  INTEGER, PARAMETER                               :: na_default=291

  REAL(wp), DIMENSION(:,:), ALLOCATABLE            :: gradient_bangle
  REAL(wp), DIMENSION(:,:), ALLOCATABLE            :: gradient_refrac
  REAL(wp), DIMENSION(:), ALLOCATABLE              :: tdry
  REAL(wp)                                         :: geop_min=zmin_default
  REAL(wp)                                         :: geop_max=zmax_default
  INTEGER                                          :: n_geop=nz_default
  REAL(wp)                                         :: ih_min=amin_default
  REAL(wp)                                         :: ih_max=amax_default
  INTEGER                                          :: n_ih=na_default

  INTEGER                                          :: idummy
  INTEGER                                          :: i, iargc, argc, k
  INTEGER                                          :: n_files, n_profiles
  INTEGER                                          :: namelist_lun
  INTEGER                                          :: iostatus

  LOGICAL                                          :: refrac_wanted = .TRUE.
  LOGICAL                                          :: bangle_wanted = .TRUE.
  LOGICAL                                          :: ih_from_geop  = .TRUE.
  LOGICAL                                          :: use_logp      = .FALSE.
  LOGICAL                                          :: use_logq      = .FALSE.
  LOGICAL                                          :: calc_grad     = .TRUE.
  LOGICAL                                          :: lfile_exists  = .FALSE.
  LOGICAL                                          :: cfile_exists  = .FALSE.
  LOGICAL                                          :: ranchk        = .TRUE.
  LOGICAL                                          :: compress      = .FALSE.
  LOGICAL                                          :: check_qsat    = .FALSE.
  LOGICAL                                          :: check_qmin    = .TRUE.
  LOGICAL                                          :: new_op        = .FALSE.
  LOGICAL                                          :: direct_ion    = .FALSE.
  LOGICAL                                          :: use_247L      = .FALSE.
  LOGICAL                                          :: give_help     = .FALSE.

  CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
  CHARACTER(len = 4096)                            :: ofile='bg2ro.nc'
  CHARACTER(len = 4096)                            :: lfile='' ! levels file
  CHARACTER(len = 4096)                            :: cfile='' ! config file
  CHARACTER(len =  256)                            :: buffer
  CHARACTER(len =    4)                            :: istr
  CHARACTER(len =    6)                            :: nstr
  CHARACTER(len =    4)                            :: slev1b, slev2a

! Configuration namelist
  NAMELIST / config_fm / refrac_wanted, geop_min, geop_max, n_geop, &
                         bangle_wanted, use_247L, ih_from_geop, ih_min, ih_max, n_ih, &
                         lfile, &
                         calc_grad, use_logp, use_logq, compress, &
                         check_qsat, check_qmin, new_op, direct_ion

!-------------------------------------------------------------------------------
! 2. Default settings
!-------------------------------------------------------------------------------

  CALL message_set_routine ( 'ropp_fm_bg2ro_1d' )

  CALL message(msg_noin, '')
  CALL message(msg_noin, &
      '-----------------------------------------------------------------------')
  CALL message(msg_noin, &
      '                     ROPP Forward Model'                                )
  CALL message(msg_noin, &
      '-----------------------------------------------------------------------')
  CALL message(msg_noin, '')

!-------------------------------------------------------------------------------
! 3. Read configuration file, if given and possible
!-------------------------------------------------------------------------------

! 3.1 Identify the configuration file
! -----------------------------------

  i = 1

  DO WHILE ( i <= iargc() )

    CALL getarg(i, buffer)

    IF ( buffer == '-c' ) THEN
      CALL getarg(i+1, buffer)
      cfile = buffer
    END IF

    IF ( buffer == '-d' ) THEN
      msg_MODE = VerboseMode    ! Needed in the next section
    END IF

    i = i + 1

  END DO

! 3.2 Read the configuration file
! -------------------------------

  IF ( TRIM(ADJUSTL(cfile)) /= '' ) THEN

    INQUIRE ( FILE=cfile, EXIST=cfile_exists )

    IF ( cfile_exists ) THEN

      CALL message ( msg_info, 'Reading config namelist from  ' // &
                               TRIM(ADJUSTL(cfile)) )

      namelist_lun = Get_IO_Unit()  ! From ropp_utils

      OPEN ( UNIT=namelist_lun, FILE=cfile, STATUS='OLD', ACTION='READ', IOSTAT=iostatus )

      IF ( iostatus > 0 ) &
        CALL message ( msg_fatal, 'I/O error while opening config file ' // &
                                  TRIM(ADJUSTL(cfile)) )

      READ ( UNIT=namelist_lun, NML=config_fm, IOSTAT=iostatus )

      IF ( iostatus > 0 ) &
        CALL message ( msg_fatal, 'I/O error while reading config namelist ' // &
                                  TRIM(ADJUSTL(cfile)) )

      IF ( msg_MODE == VerboseMode ) THEN

        CALL message( msg_diag, 'Contents of config_fm namelist:' )
        WRITE ( *, NML=config_fm )

      ENDIF

      CLOSE ( UNIT=namelist_lun )

    ELSE

      CALL message( msg_fatal, 'Requested configuration file ' // &
                               TRIM(ADJUSTL(cfile)) // ' does not exist' )

    END IF

  END IF

!-------------------------------------------------------------------------------
! 4. Command line arguments, which therefore override those in the namelist
!-------------------------------------------------------------------------------

  argc = iargc()
  i = 1
  n_files = 0
  ALLOCATE(ifiles(argc))

  DO WHILE(i <= argc)

    CALL getarg(i, buffer)

    SELECT CASE (buffer)

        CASE('-o')                          ! Output file name (netCDF output)
           CALL getarg(i+1, buffer)
           ofile = buffer
           i = i + 1

        CASE ('-l')                         ! Use level structure defined here
           CALL getarg(i+1, buffer)
           lfile = buffer
           i = i + 1

        CASE ('-f')                         ! Perform forward model only
           calc_grad = .FALSE.

        CASE ('-use_logp')                  ! Use log(pressure) for FM
           state%use_logp = .TRUE.

        CASE ('-use_logq')                  ! Use log(shum) for FM
           state%use_logq = .TRUE.

        CASE ('-comp')                      ! Non ideal gas switch
           compress = .TRUE.

        CASE ('-check_qsat')                ! Check against supersaturation
           check_qsat = .TRUE.

        CASE ('-nocheck_qmin')              ! Do not check against superdryness
           check_qmin = .FALSE.

        CASE ('-new_op')                    ! New interpolation
           new_op = .TRUE.

        CASE ('-direct_ion')                ! Model L1 and L2 bangles directly
           direct_ion = .TRUE.

        CASE ('-best', '-Best', '-BEST')    ! 'Best' options, namely ...
           calc_grad = .FALSE.
           compress = .TRUE.
           new_op = .TRUE.
           use_247L = .TRUE.
           ih_from_geop = .FALSE.

        CASE ('--no-ranchk')                ! Disable rangecheck; for experts only
                                            ! Do not document as a user option.
           CALL message ( msg_warn, 'Range checking is disabled' )
           ranchk = .FALSE.

        CASE ('-refrac_only' )              ! Don't generate bending angles
           bangle_wanted = .FALSE.

        CASE ('-zmin', '-Zmin', '-ZMIN')    ! Min refrac geopotential
           CALL getarg(i+1, buffer)
           READ(buffer, *) geop_min
           i = i + 1

        CASE ('-zmax', '-Zmax', '-ZMAX')    ! Max refrac geopotential
           CALL getarg(i+1, buffer)
           READ(buffer, *) geop_max
           i = i + 1

        CASE ('-nz', '-Nz', '-NZ')          ! Number of refractivity levels
           CALL getarg(i+1, buffer)
           READ(buffer, *) n_geop
           i = i + 1

        CASE ('-bangle_only' )              ! Don't generate refractivities
           refrac_wanted = .FALSE.

        CASE ('-247L', '-247l', '-247')     ! Output bending angles on standard 247L
           use_247L = .TRUE.
           ih_from_geop = .FALSE.

        CASE ('-ih_from_geop')              ! Make impact heights correspond to refractivity geopotentials
           ih_from_geop = .TRUE.

        CASE ('-ihmin', '-IHmin', '-IHMIN') ! Min impact height
           CALL getarg(i+1, buffer)
           READ(buffer, *) ih_min
           i = i + 1

        CASE ('-ihmax', '-IHmax', '-IHMAX') ! Max impact height
           CALL getarg(i+1, buffer)
           READ(buffer, *) ih_max
           i = i + 1

        CASE ('-nih', '-Nih', '-NIH')       ! Number of impact heights
           CALL getarg(i+1, buffer)
           READ(buffer, *) n_ih
           i = i + 1

        CASE('-d')                          ! Additional diagnostic mode
           msg_MODE = VerboseMode

        CASE('-h', '--help', '?')           ! Help
           give_help = .TRUE.

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

        CASE default                        ! Input file name
           IF ( buffer(1:1) /= '-' ) THEN
             n_files = n_files + 1
             ifiles(n_files) = buffer
           END IF

    END SELECT

    i = i + 1

  END DO

  IF ( (n_files == 0) .AND. (.NOT. give_help) ) &
    CALL message ( msg_error, 'No input file(s) specified' )

  IF ( argc == 0 .OR. n_files == 0 .OR. give_help ) THEN
    CALL usage()
    CALL EXIT(msg_exit_status)
  END IF

!-------------------------------------------------------------------------------
! 5. Remove pre-existing output file
!-------------------------------------------------------------------------------

  CALL file_delete(ofile, idummy)

!-------------------------------------------------------------------------------
! 6. Set default units
!-------------------------------------------------------------------------------

  CALL ropp_fm_set_units(ro_data)

!-------------------------------------------------------------------------------
! 7. Loop over all input files
!-------------------------------------------------------------------------------

  DO k = 1, n_files

!-------------------------------------------------------------------------------
! 8. Loop over all profiles
!-------------------------------------------------------------------------------

    n_profiles = ropp_io_nrec(ifiles(k))

    DO i = 1, n_profiles

      WRITE(istr, '(i4)') i
      WRITE(nstr, '(i6)') n_profiles
      CALL message(msg_noin, '')
      CALL message(msg_info, 'Processing profile ' // istr // ' of ' // nstr )

!-------------------------------------------------------------------------------
! 9. Read data
!-------------------------------------------------------------------------------

! 9.1 Background data
! -------------------

      CALL ropp_io_read(ro_data, ifiles(k), rec=i, ranchk=ranchk)

      CALL message(msg_info, '(' // TRIM(ro_data%occ_id) // ') \n')

! 9.2 Levels data, if given and possible
! --------------------------------------

      IF ( TRIM(ADJUSTL(lfile)) /= '' ) THEN

        INQUIRE ( FILE=lfile, EXIST=lfile_exists )

        IF ( lfile_exists ) THEN

          CALL message(msg_info, 'Reading levels information from ' // lfile // ' \n')

          CALL ropp_fm_set_units(l_data) ! l_data must use the same units as ro_data

          CALL ropp_io_read(l_data, lfile, ranchk=ranchk)

          ro_data%georef = l_data%georef

          IF ( bangle_wanted ) THEN  ! only use the heights from the levels file
            CALL ropp_io_init(ro_data%lev1b, l_data%lev1b%npoints)
            ro_data%lev1b%impact = l_data%lev1b%impact
          END IF

          IF ( refrac_wanted ) THEN  ! only use the heights from the levels file
            CALL ropp_io_init(ro_data%lev2a, l_data%lev2a%npoints)
            ro_data%lev2a%geop_refrac = l_data%lev2a%geop_refrac
          END IF

          CALL ropp_io_free(l_data)

          WRITE (slev1b, FMT='(I4)') ro_data%lev1b%Npoints
          CALL message(msg_info, 'There are ' // TRIM(ADJUSTL(slev1b)) // &
                                 ' level 1b points in ' // lfile // ' \n')

          WRITE (slev2a, FMT='(I4)') ro_data%lev2a%Npoints
          CALL message(msg_info, 'There are ' // TRIM(ADJUSTL(slev2a)) // &
                                 ' level 2a points in ' // lfile // ' \n')

        ELSE

          CALL message(msg_fatal, 'Requested levels file ' // &
                                   TRIM(ADJUSTL(lfile)) // ' does not exist')

        END IF

      END IF

!-------------------------------------------------------------------------------
! 10. Set default ionospheric state vector elements if modelling L1 and L2 directly
!-------------------------------------------------------------------------------

      IF ( direct_ion ) THEN

        ro_data%lev2c%direct_ion = .TRUE.

        CALL ropp_fm_iono_set_default(ro_data)

      END IF

!-------------------------------------------------------------------------------
! 11. Copy data in RO structure to state, and set relevant flags
!-------------------------------------------------------------------------------

!  state%check_qsat and state%check_qmin are both used by ropp_fm_roprof2state,
!  and must therefore be set before calling it.

      IF ( check_qsat ) state%check_qsat = .TRUE.

      IF ( .NOT. check_qmin ) state%check_qmin = .FALSE.

      CALL ropp_fm_roprof2state(ro_data, state)

      IF ( compress ) state%non_ideal = .TRUE.  ! Switch on non-ideal gas compressibility option

      IF ( new_op ) THEN
        state%new_ref_op    = .TRUE.            ! Use new refrac interpolation
        state%new_bangle_op = .TRUE.            ! Use new bangle interpolation
      END IF

!-------------------------------------------------------------------------------
! 12. Generate refractivities, if desired or required
!-------------------------------------------------------------------------------

      IF ( refrac_wanted .OR. ih_from_geop ) THEN

! 12.1 Levels
! -----------

        IF ( ro_data%lev2a%Npoints > 0 ) THEN   ! Set the refrac levels to those in the input file or levels file

          IF ( lfile_exists ) THEN ! Take from levels file

            CALL set_obs_levels_refrac_file(ro_data, obs_refrac)

          ELSE ! Take from input file

            CALL ropp_fm_roprof2obs(ro_data, obs_refrac)

          END IF

        ELSE                                    ! Generate uniformly spaced refrac levels

          CALL set_obs_levels_refrac_uni(ro_data, obs_refrac)

        END IF

! 12.2 Refractivities
! -------------------

        IF ( state%state_ok ) THEN

          IF ( state%new_ref_op ) THEN
            CALL ropp_fm_refrac_1d_new(state, obs_refrac)
          ELSE
            CALL ropp_fm_refrac_1d(state, obs_refrac)
          END IF

! 12.3 Refractivity gradients
! ---------------------------

          IF ( calc_grad .AND. refrac_wanted ) THEN

            ALLOCATE(gradient_refrac(SIZE(obs_refrac%refrac), SIZE(state%state)))

            CALL ropp_fm_refrac_1d_grad(state, obs_refrac, gradient_refrac)

          END IF

! 12.4 Dry temperatures
! ---------------------

          ALLOCATE( tdry( SIZE(obs_refrac%refrac) ) )

          tdry = ropp_MDFV

          CALL ropp_fm_tdry_1d(state, obs_refrac, tdry)


        END IF  ! state%state_ok

      END IF  ! refrac_wanted

!-------------------------------------------------------------------------------
! 13. Generate bending angles, if desired
!-------------------------------------------------------------------------------

      IF ( bangle_wanted ) THEN

! 13.1 Levels
! -----------

        IF ( ro_data%lev1b%Npoints > 0 ) THEN   ! Set the bangle levels to those in the input file or levels file

          IF ( lfile_exists ) THEN ! Take from levels file

            CALL set_obs_levels_bangle_file(ro_data, obs_bangle)

          ELSE ! Take from input file

            CALL ropp_fm_roprof2obs(ro_data, obs_bangle)

          END IF

        ELSE

          IF ( ih_from_geop .AND. (.NOT. use_247L) ) THEN

            CALL set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle) ! Calculate from refrac levels

          ELSE

            CALL set_obs_levels_bangle_uni(ro_data, obs_bangle) ! 247L or uniform

          END IF

        END IF

! 13.2 Bending angles
!--------------------

        IF ( state%state_ok ) THEN

          CALL ropp_fm_bangle_1d(state, obs_bangle)

! 13.3 Bending angle gradients
!-----------------------------

          IF ( calc_grad .AND. (.NOT. direct_ion) ) THEN

            ALLOCATE(gradient_bangle(SIZE(obs_bangle%bangle),SIZE(state%state)))

            CALL ropp_fm_bangle_1d_grad(state, obs_bangle, gradient_bangle)

          END IF

        END IF  ! state%state_ok = .TRUE.

      END IF  ! bangle_wanted = .TRUE.

!-------------------------------------------------------------------------------
! 14. Copy simulated observations to RO structure and write data
!-------------------------------------------------------------------------------

      IF ( refrac_wanted ) THEN

        CALL ropp_fm_obs2roprof(obs_refrac, ro_data)

        ro_data%Lev2a%dry_temp      = tdry
        ro_data%Lev2a%dry_temp_qual = 100.0_wp

      END IF

      IF ( bangle_wanted ) CALL ropp_fm_obs2roprof(obs_bangle, ro_data)

      IF (state%state_ok) THEN

        IF ( calc_grad .AND. refrac_wanted ) THEN

          CALL ropp_io_addvar(ro_data,                                       &
                              name      = 'gradient_refrac',                 &
                              long_name =                                    &
                              'Gradient of the refractivity forward model',  &
                              units     = 'N-unit',                          &
                              range     = (/MINVAL(gradient_refrac),         &
                                            MAXVAL(gradient_refrac)/),       &
                              DATA      = gradient_refrac)

        END IF


        IF ( calc_grad .AND. bangle_wanted .AND. (.NOT. direct_ion) ) THEN

          CALL ropp_io_addvar(ro_data,                                       &
                              name      = 'gradient_bangle',                 &
                              long_name =                                    &
                              'Gradient of the bending angle forward model', &
                              units     = 'rad',                             &
                              range     = (/MINVAL(gradient_bangle),         &
                                            MAXVAL(gradient_bangle)/),       &
                              DATA      = gradient_bangle)

        END IF

!-------------------------------------------------------------------------------
! 15. Update RO structure with computed state vector variables
!-------------------------------------------------------------------------------

        CALL ropp_fm_state2roprof(state, ro_data)

      END IF

!-------------------------------------------------------------------------------
! 16. If modelling L1 and L2 directly, unpack concatenated bending angles
!-------------------------------------------------------------------------------

      IF ( direct_ion ) CALL ropp_fm_iono_unpack(ro_data)

!-------------------------------------------------------------------------------
! 17. Write data
!-------------------------------------------------------------------------------

      CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk)

!-------------------------------------------------------------------------------
! 18. Clean up
!-------------------------------------------------------------------------------

      IF ( ALLOCATED(tdry) ) THEN
        DEALLOCATE(tdry)
      END IF

      IF ( state%state_ok .AND. calc_grad .AND. refrac_wanted ) THEN
        DEALLOCATE (gradient_refrac)
      END IF

      IF ( state%state_ok .AND. calc_grad .AND. (.NOT. direct_ion) .AND. bangle_wanted ) THEN
        DEALLOCATE (gradient_bangle)
      END IF

      CALL ropp_fm_free(state)

      CALL ropp_fm_free(obs_refrac)

      CALL ropp_fm_free(obs_bangle)

      CALL ropp_io_free(ro_data)

    END DO

  END DO

  CALL EXIT(msg_exit_status)


CONTAINS


!-------------------------------------------------------------------------------
! 19. Get observation levels for refractivity from file
!     (Basically does the same job as roprof2obs, but setting refrac to ropp_MDFV
!      and without copying covariances.)
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_refrac_file(ro_data, obs_refrac)

!   19.1 Declarations
!   -----------------

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_io_types
    USE ropp_fm

    IMPLICIT NONE

!   Input/Output variables
    TYPE(ROprof)                     :: ro_data
    TYPE(Obs1dRefrac)                :: obs_refrac

!   Local variables
    CHARACTER(len=256)               :: routine

    CALL message_get_routine(routine)
    CALL message_set_routine('set_obs_levels_refrac_file')

!   19.2 Copy obs levels directly from ro_data structure
!   ----------------------------------------------------

    ALLOCATE(obs_refrac%refrac(ro_data%Lev2a%Npoints))
    ALLOCATE(obs_refrac%geop(ro_data%Lev2a%Npoints))
    ALLOCATE(obs_refrac%weights(ro_data%Lev2a%Npoints))

    obs_refrac%refrac(:)  = ropp_MDFV
    obs_refrac%geop(:)    = ro_data%Lev2a%geop_refrac
    obs_refrac%weights(:) = 1.0_wp
    obs_refrac%lon        = ro_data%GEOref%lon
    obs_refrac%lat        = ro_data%GEOref%lat
    obs_refrac%obs_ok     = .TRUE.
    obs_refrac%cov_ok     = .FALSE.

    CALL message_set_routine(routine)

  END SUBROUTINE set_obs_levels_refrac_file


!-------------------------------------------------------------------------------
! 20. Calculate uniformly spaced observation levels for refractivity
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_refrac_uni(ro_data, obs_refrac)

!   Generate n_geop uniformly spaced geopotential heights between
!   geop_min and geop_max (gpm).  These parameters are inherited from the
!   main program.

!   20.1 Declarations
!   -----------------

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_io_types
    USE ropp_fm

    IMPLICIT NONE

!   Input/Output variables
    TYPE(ROprof)                     :: ro_data
    TYPE(Obs1dRefrac)                :: obs_refrac

!   Local variables
    REAL(wp)                         :: zmin, zmax
    INTEGER                          :: j, nz
    CHARACTER(len=256)               :: routine

!   20.2 Define routine name for messaging
!   --------------------------------------

    CALL message_get_routine(routine)
    CALL message_set_routine('set_obs_levels_refrac_uni')

!   20.3 Generate the geopotential heights
!   --------------------------------------

    zmin = geop_min
    zmax = geop_max
    nz   = n_geop

    IF ( (zmin >= zmax) .OR. (nz <= 1) ) THEN
      CALL message(msg_warn, 'Invalid geop_min, geop_max or n_geop' // &
                             ' ... resetting all to default values')
      zmin = zmin_default
      zmax = zmax_default
      nz   = nz_default
    END IF

    ALLOCATE(obs_refrac%refrac(nz))
    ALLOCATE(obs_refrac%geop(nz))
    ALLOCATE(obs_refrac%weights(nz))

    obs_refrac%refrac(:)  = 0.0_wp
    obs_refrac%geop(:)    = (/ (zmin+(j-1)*(zmax-zmin)/REAL(nz-1, wp), j=1,nz) /)
    obs_refrac%weights(:) = 1.0_wp
    obs_refrac%lon        = ro_data%GEOref%lon
    obs_refrac%lat        = ro_data%GEOref%lat
    obs_refrac%obs_ok     = .TRUE.

!   20.4 Clean up
!   -------------

    CALL message_set_routine(routine)

  END SUBROUTINE set_obs_levels_refrac_uni


!-------------------------------------------------------------------------------
! 21. Get obs levels for bending angle from file
!     (Basically does the same job as roprof2obs, but setting bangle to ropp_MDFV
!      and without copying covariances.)
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_bangle_file(ro_data, obs_bangle)

!   21.1 Declarations
!   -----------------

    USE typesizes, ONLY: wp => EightByteReal
    USE ropp_io_types
    USE ropp_fm

    IMPLICIT NONE

!   Input/Output variables
    TYPE(ROprof)                     :: ro_data
    TYPE(Obs1dBangle)                :: obs_bangle

!   Local variables
    CHARACTER(len=256)               :: routine

    CALL message_get_routine(routine)
    CALL message_set_routine('set_obs_levels_bangle_file')

!   21.2 Copy obs levels directly from ro_data (holding data from file)
!  --------------------------------------------------------------------

    ALLOCATE(obs_bangle%bangle(ro_data%Lev1b%Npoints))
    ALLOCATE(obs_bangle%impact(ro_data%Lev1b%Npoints))
    ALLOCATE(obs_bangle%weights(ro_data%Lev1b%Npoints))

    obs_bangle%bangle(:)  = ropp_MDFV
    obs_bangle%impact(:)  = ro_data%Lev1b%impact
    obs_bangle%weights(:) = 1.0_wp
    obs_bangle%lon        = ro_data%GEOref%lon
    obs_bangle%lat        = ro_data%GEOref%lat
    obs_bangle%r_curve    = ro_data%GEOref%roc
    obs_bangle%undulation = ro_data%GEOref%undulation
    obs_bangle%azimuth    = ro_data%GEOref%azimuth
    obs_bangle%obs_ok     = .TRUE.
    obs_bangle%cov_ok     = .FALSE.

    IF (isinrange(ro_data%GEOref%roc, ro_data%GEOref%range%roc)) THEN
      obs_bangle%r_curve = ro_data%GEOref%roc

    ELSE
      CALL message(msg_warn, "Radius of curvature out of range or missing.")
      CALL message(msg_warn, "Check input data and valid_range attributes")
      obs_bangle%obs_ok = .FALSE.
    ENDIF

    IF (isinrange(ro_data%GEOref%undulation, ro_data%GEOref%range%undulation)) THEN
      obs_bangle%undulation = ro_data%GEOref%undulation
    ELSE
      CALL message(msg_warn, "Undulation out of range or missing.")
      CALL message(msg_warn, "Check input data and valid_range attributes")
      obs_bangle%obs_ok = .FALSE.
    ENDIF

    CALL message_set_routine(routine)

  END SUBROUTINE set_obs_levels_bangle_file


!-------------------------------------------------------------------------------
! 22. Calculate obs levels for bending angle (consistent with obs_refrac)
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle)

!   Generate SIZE(obs_refrac%geop) impact parameters corresponding to
!   given refracs and their geopotential heights.

!   22.1 Declarations
!   -----------------

    USE typesizes, ONLY: wp => EightByteReal
    USE geodesy
    USE ropp_io_types
    USE ropp_fm

    IMPLICIT NONE

    TYPE(ROprof)                        :: ro_data
    TYPE(Obs1dRefrac)                   :: obs_refrac
    TYPE(Obs1dbangle)                   :: obs_bangle
    TYPE(State1dFM)                     :: state

    REAL(wp), DIMENSION(:), ALLOCATABLE :: tmp
    REAL(wp)                            :: max_refrac
    REAL(wp)                            :: min_geop

    INTEGER                             :: j, nlev, imiss
    INTEGER                             :: lowest_refrac_index

    CHARACTER(LEN=256)                  :: routine

!   22.2 Define routine name for messaging
!   ---------------------------------------

    CALL message_get_routine(routine)
    CALL message_set_routine('set_obs_levels_bangle')

!   22.3 Allocate arrays
!   --------------------

    max_refrac = ro_data%lev2a%range%refrac(2)

    min_geop   = ro_data%lev2b%range%geop(1)

    nlev = SIZE(obs_refrac%geop)

    obs_bangle%n_L1 = nlev

    ALLOCATE(tmp(nlev))

    IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
      ALLOCATE(obs_bangle%bangle(nlev))
      ALLOCATE(obs_bangle%impact(nlev))
      ALLOCATE(obs_bangle%weights(nlev))
    ELSE
      ALLOCATE(obs_bangle%bangle(2*nlev))
      ALLOCATE(obs_bangle%impact(2*nlev))
      ALLOCATE(obs_bangle%weights(2*nlev))
    END IF

!   22.4 Calculate refractivities
!   -----------------------------

    IF ( state%new_ref_op ) THEN
      CALL ropp_fm_refrac_1d_new(state, obs_refrac)
    ELSE
      CALL ropp_fm_refrac_1d(state, obs_refrac)
    END IF

    IF ( ANY(obs_refrac%refrac <= ropp_MDTV) ) THEN ! rederive approximately the zapped refracs

      lowest_refrac_index = SUM( MINLOC( (/ (j, j=1,nlev) /), &
                                         MASK=obs_refrac%refrac>ropp_MDTV ) )

      IF ( (lowest_refrac_index > 1) .AND. (lowest_refrac_index < nlev) ) THEN

        DO imiss=lowest_refrac_index-1,1,-1
          obs_refrac%refrac(imiss) = obs_refrac%refrac(imiss+1) + &
            (obs_refrac%refrac(lowest_refrac_index) - obs_refrac%refrac(lowest_refrac_index+1))
        END DO

      END IF

    END IF

    WHERE ( obs_refrac%refrac > max_refrac ) obs_refrac%refrac = max_refrac

    WHERE ( obs_refrac%geop   < min_geop   ) obs_refrac%geop   = min_geop

!   22.5 Set scalar arguments of the observation vector
!   ---------------------------------------------------

    obs_bangle%g_sfc        = gravity(ro_data%GEOref%lat)
    obs_bangle%r_earth      = R_eff(ro_data%GEOref%lat)
    obs_bangle%undulation   = ro_data%GEOref%undulation
    obs_bangle%lon          = ro_data%GEOref%lon
    obs_bangle%lat          = ro_data%GEOref%lat
    obs_bangle%azimuth      = ro_data%GEOref%azimuth

    IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
      obs_bangle%r_curve = ro_data%GEOref%roc
    ELSE
      CALL message (msg_warn, 'RoC missing from data structure ... ' // &
                              'setting equal to effective radius')
      obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
    END IF

    IF ( ro_data%Lev2c%direct_ion ) CALL calculate_r_leo(ro_data, obs_bangle)

!   22.6 Calculate levels to coincide with the geopotential levels
!   --------------------------------------------------------------

    tmp = geopotential2geometric(ro_data%GEOref%lat, obs_refrac%geop) + &
                                 obs_bangle%r_curve

    IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
      tmp = tmp + obs_bangle%undulation
    ELSE
      CALL message(msg_warn, 'Undulation missing ... will assume to ' // &
                             'be zero when calculating impact parameters.')
    END IF

    tmp = (1.0_wp + obs_refrac%refrac*1.e-6_wp) * tmp

    IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
      obs_bangle%impact                = tmp
    ELSE
      obs_bangle%impact(1:nlev)        = tmp
      obs_bangle%impact(nlev+1:2*nlev) = tmp
    END IF

!   22.7 Fill other arrays
!   ----------------------

    obs_bangle%bangle(:)  = 0.0_wp
    obs_bangle%weights(:) = 1.0_wp
    obs_bangle%obs_ok     = .TRUE.

!   22.8 Clean up
!   -------------

!   Regenerate refracs with missing data
    IF ( state%new_ref_op ) THEN
      CALL ropp_fm_refrac_1d_new(state, obs_refrac)
    ELSE
      CALL ropp_fm_refrac_1d(state, obs_refrac)
    END IF

    DEALLOCATE(tmp)

    CALL message_set_routine(routine)

  END SUBROUTINE set_obs_levels_bangle


!-------------------------------------------------------------------------------
! 23. Calculate obs levels for bending angle (independent of obs_refrac)
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_bangle_uni(ro_data, obs_bangle)

!   Define impact heights that are independent of the refractivity.
!   If use_247L is .TRUE., these are the 'standard' 247 impact heights.
!   If not, the n_ih impact heights are uniformly spaced between ih_min and ih_max.
!   These parameters are inherited from the main program.


!   23.1 Declarations
!   -----------------

    USE typesizes, ONLY: wp => EightByteReal
    USE geodesy
    USE ropp_io_types
    USE ropp_fm

    IMPLICIT NONE

!   Input/output
    TYPE(ROprof)                        :: ro_data
    TYPE(Obs1dbangle)                   :: obs_bangle

!   Local
    INTEGER                             :: j, nlev
    REAL(wp)                            :: amin, amax
    INTEGER                             :: na
    REAL(wp), DIMENSION(:), ALLOCATABLE :: impact_height

    CHARACTER(LEN=256)                  :: routine

!   23.2 Define routine name for messaging
!   --------------------------------------

    CALL message_get_routine(routine)
    CALL message_set_routine('set_obs_levels_bangle_uni')

!   23.3 Calculate number of levels, nlev
!   -------------------------------------

    IF ( use_247L ) THEN

      nlev = 247

    ELSE

      amin = ih_min
      amax = ih_max
      na   = n_ih

      IF ( (amin >= amax) .OR. (na <= 1) ) THEN
        CALL message(msg_warn, 'Invalid ih_min, ih_max or n_ih' // &
                               ' ... resetting all to default values')
        amin = amin_default
        amax = amax_default
        na   = na_default
      END IF

      nlev = na

    END IF

!   23.4 Allocate arrays
!   --------------------

    obs_bangle%n_L1 = nlev

    ALLOCATE(impact_height(nlev))

    IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
      ALLOCATE(obs_bangle%bangle(nlev))
      ALLOCATE(obs_bangle%impact(nlev))
      ALLOCATE(obs_bangle%weights(nlev))
    ELSE
      ALLOCATE(obs_bangle%bangle(2*nlev))
      ALLOCATE(obs_bangle%impact(2*nlev))
      ALLOCATE(obs_bangle%weights(2*nlev))
    END IF

!   23.5 Define impact heights
!   --------------------------

    IF ( use_247L ) THEN  ! The 'standard' 247 impact heights

      impact_height = impact_height_eum_247

    ELSE                  ! Define uniform set of impact heights

      impact_height = amin + ((amax - amin)/REAL(na-1, wp)) * (/ (j, j=0,na-1) /)

    END IF

!   23.6 Set scalar arguments of the observation vector
!   ---------------------------------------------------

    obs_bangle%g_sfc      = gravity(ro_data%GEOref%lat)
    obs_bangle%r_earth    = R_eff(ro_data%GEOref%lat)
    obs_bangle%undulation = ro_data%GEOref%undulation
    obs_bangle%lon        = ro_data%GEOref%lon
    obs_bangle%lat        = ro_data%GEOref%lat
    obs_bangle%azimuth    = ro_data%GEOref%azimuth

    IF ( ro_data%GEOref%roc > 0.0_wp ) THEN
      obs_bangle%r_curve = ro_data%GEOref%roc
    ELSE
      CALL message (msg_warn, 'RoC missing from data structure ... ' // &
                              'setting equal to effective radius')
      obs_bangle%r_curve = obs_bangle%r_earth ! Maybe +21km? See RSR 14.
    END IF

    IF ( ro_data%Lev2c%direct_ion ) CALL calculate_r_leo(ro_data, obs_bangle)

!   23.7 Calculate impact parameters
!   --------------------------------

    IF ( .NOT. ro_data%Lev2c%direct_ion ) THEN
      obs_bangle%impact                = impact_height + obs_bangle%r_curve
    ELSE
      obs_bangle%impact(1:nlev)        = impact_height + obs_bangle%r_curve
      obs_bangle%impact(nlev+1:2*nlev) = impact_height + obs_bangle%r_curve
    END IF

!   23.8 Fill other arrays
!   ----------------------

    obs_bangle%bangle(:)  = 0.0_wp
    obs_bangle%weights(:) = 1.0_wp

!   23.9 Clean up
!   -------------

    DEALLOCATE (impact_height)

    CALL message_set_routine(routine)

  END SUBROUTINE set_obs_levels_bangle_uni

!-------------------------------------------------------------------------------
! 24. Calculate r_LEO
!-------------------------------------------------------------------------------

  SUBROUTINE calculate_r_LEO(ro_data, y)

!   Derive r_LEO from LEO POD (ideally), or estimate from satellite ID.

!   24.1 Declarations
!   -----------------

    USE ropp_utils, ONLY: ropp_MDFV
    USE messages
    USE ropp_fm_types
    USE ropp_io_types

    IMPLICIT NONE

!   Input/output
    TYPE(ROprof), INTENT(INOUT)      :: ro_data
    TYPE(Obs1dbangle), INTENT(INOUT) :: y

!   Local
    CHARACTER(LEN=6)                 :: sr_leo
    CHARACTER(LEN=256)               :: routine

!   24.2 Define routine name for messaging
!   --------------------------------------

    CALL message_get_routine(routine)
    CALL message_set_routine('calculate_r_LEO')

!   24.3 Calculate r_LEO
!   --------------------

    IF ( ANY ( ABS(ro_data%georef%leo_pod%pos - ropp_MDFV) < 1.0_wp ) .OR. &
         ANY ( ABS(ro_data%georef%r_coc       - ropp_MDFV) < 1.0_wp ) ) THEN

      CALL message (msg_warn, 'Cannot find LEO position and/or centre of curvature ... ' // &
                              'will try to estimate r_LEO from satellite ID.')

      SELECT CASE (ro_data%leo_id(1:2))  ! Data from https://www.wmo-sat.info/oscar/

        CASE ('C0')
          CALL message (msg_info, 'Assuming LEO height of 800 km (COSMIC profile)')
          y%r_leo = y%r_curve + 800.E3_wp

        CASE ('CH')
          CALL message (msg_info, 'Assuming LEO height of 470 km (CHAMP profile)')
          y%r_leo = y%r_curve + 470.E3_wp

        CASE ('FY')
          CALL message (msg_info, 'Assuming LEO height of 840 km (FY-3C/D/... profile)')
          y%r_leo = y%r_curve + 840.E3_wp

        CASE ('GR')
          CALL message (msg_info, 'Assuming LEO height of 340 km (GRACE-A/B profile)')
          y%r_leo = y%r_curve + 340.E3_wp

        CASE ('ME')
          CALL message (msg_info, 'Assuming LEO height of 830 km (Metop-A/B/C profile)')
          y%r_leo = y%r_curve + 830.E3_wp

        CASE DEFAULT
          CALL message(msg_warn, 'Unrecognised LEO satellite ID (' // ro_data%leo_id // &
                                 ') ... r_LEO cannot be estimated.')
          CALL message(msg_warn, 'Check that LEO POD and CoC positions, or LEO ID, are valid.')
          CALL message(msg_warn, 'Setting status flag obs%ok to FALSE')
          y%r_leo =  ropp_MDFV

      END SELECT

    ELSE

      SELECT CASE ( TRIM(ADJUSTL(ro_data%georef%units%r_coc)) )

        CASE ('km', 'Km', 'KM', 'kilometres', 'kilometre', 'Kilometres', 'Kilometre') ! Likely
          y%r_leo = SQRT ( SUM ( (ro_data%georef%leo_pod%pos - &
                                  1.E3_wp*ro_data%georef%r_coc)**2 ) )

        CASE DEFAULT  ! ro_data%georef%r_coc assumed to be in m, like ro_data%georef%leo_pod%pos
          y%r_leo = SQRT ( SUM ( (ro_data%georef%leo_pod%pos - &
                                  ro_data%georef%r_coc)**2 ) )

      END SELECT

      WRITE (sr_leo, FMT='(F6.1)') (y%r_leo - y%r_curve) * 1.0E-3_wp
      CALL message (msg_info, 'Using a LEO altitude of ' // sr_leo // ' km.')

    END IF

!   24.4 Clean up
!   -------------

    CALL message_set_routine(routine)

  END SUBROUTINE calculate_r_LEO


!-------------------------------------------------------------------------------
! 25. Usage information
!-------------------------------------------------------------------------------

  SUBROUTINE usage()
    PRINT *, 'Purpose:'
    PRINT *, '  Bending angles and refractivity forward model'
    PRINT *, 'Usage:'
    PRINT *, '  > ropp_fm_bg2ro_1d [<options>] <input_file(s)>'
    PRINT *, 'Options:'
    PRINT *, '  -o <output_file>  name of ROPP netCDF output file'
    PRINT *, '  -l <levels_file>  optional name of (observation) file'
    PRINT *, '                    (non-multi) containing output level information'
    PRINT *, '  -c <config_file>  name of file containing the configuration namelist'
    PRINT *, '  -f                forward model only, no gradient calculation'
    PRINT *, '  -use_logp         use log(pressure) for forward model'
    PRINT *, '  -use_logq         use log(spec humidity) for forward model'
    PRINT *, '  -comp             include non ideal gas compressibility'
    PRINT *, '  -check_qsat       include check against saturation'
    PRINT *, '  -nocheck_qmin     do not include check against dryness'
    PRINT *, '  -new_op           use alternative refrac and bangle interpolation'
    PRINT *, '  -direct_ion       model L1 and L2 directly'
    PRINT *, '  -best             use currently advised "best" options (see FM User Guide)'
    PRINT *, '  -refrac_only      do not generate bending angles'
    PRINT *, '  -zmin <geop_min>  minimum refractivity geopotential (gpm)'
    PRINT *, '  -zmax <geop_max>  maximum refractivity geopotential (gpm)'
    PRINT *, '  -nz <n_geop>      number of uniformly spaced refractivity geopotentials'
    PRINT *, '  -bangle_only      do not generate refractivities'
    PRINT *, '  -247L             output bending angles on standard 247 levels'
    PRINT *, '  -ih_from_geop     infer impact heights from refractivity geopotentials'
    PRINT *, '  -ihmin <ih_min>   minimum impact height (m)'
    PRINT *, '  -ihmax <ih_max>   maximum impact height (m)'
    PRINT *, '  -nih <n_ih>       number of uniformly spaced impact heights'
    PRINT *, '  -d                output additional diagnostics'
    PRINT *, '  -h                this help'
    PRINT *, '  -v                version information'
    PRINT *, ''
  END SUBROUTINE usage


!-------------------------------------------------------------------------------
! 26. Version information
!-------------------------------------------------------------------------------

  SUBROUTINE version_info()
    CHARACTER (LEN=40) :: version
    version = ropp_fm_version()
    PRINT *, 'ropp_fm_bg2ro_1d - Bending angles and refractivity forward model.'
    PRINT *, ''
    PRINT *, 'This program is part of ROPP (FM) Release ' // TRIM(version)
    PRINT *, ''
  END SUBROUTINE version_info

END PROGRAM ropp_fm_bg2ro_1d