! $Id: ropp_fm_bg2ro_1d.f90 4432 2015-01-19 14:26:46Z 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> [-f] [-use_logp] [-use_logq] 
!                     [-comp] [-check_qsat] [-new_op] [-direct_ion] [-247L]
!                     -zmin <geop_min> -zmax <geop_max> -nz <n_geop>
!                     [-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.
!   -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 saturation
!   -new_op             use alternative refrac and bangle interpolation
!   -direct_ion         forward model L1 and L2 directly, using model ionosphere
!   -247L               output refrac and bangle on 'standard' 247L
!   -zmin <geop_min>    minimum refractivity geopotential (gpm)'
!   -zmax <geop_max>    maximum refractivity geopotential (gpm)'
!   -nz <n_geop>        number of uniformly spaced refractivity geopotentials (gpm)'
!   -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), DIMENSION(:,:), ALLOCATABLE            :: gradient_bangle
  REAL(wp), DIMENSION(:,:), ALLOCATABLE            :: gradient_refrac
  REAL(wp)                                         :: geop_min
  REAL(wp)                                         :: geop_max
  INTEGER                                          :: n_geop

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

  LOGICAL                                          :: give_help   = .FALSE.
  LOGICAL                                          :: do_adj      = .TRUE.
  LOGICAL                                          :: lfile_exist = .FALSE.
  LOGICAL                                          :: ranchk      = .TRUE.
  LOGICAL                                          :: compress    = .FALSE.
  LOGICAL                                          :: checksat    = .FALSE.
  LOGICAL                                          :: new_op      = .FALSE.
  LOGICAL                                          :: direct_ion  = .FALSE.
  LOGICAL                                          :: use_247L    = .FALSE.

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

!-------------------------------------------------------------------------------
! 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. Command line arguments
!-------------------------------------------------------------------------------

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

  geop_min =   200.0_wp
  geop_max = 60000.0_wp
  n_geop   =   300

  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
           lfile_exist = .TRUE.
        CASE ('-f')                         ! Perform forward model only
           do_adj = .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 saturation
           checksat = .TRUE.
        CASE ('-new_op')                    ! New interpolation
           new_op = .TRUE.
        CASE ('-direct_ion')                ! Model L1 and L2 bangles directly
           direct_ion = .TRUE.
        CASE ('--no-ranchk')                ! Disable rangecheck; for experts only
           CALL message ( msg_warn, "Range checking is disabled" )
           ranchk = .FALSE.                 ! Do not document as a user option.
        CASE ('-247L', '-247l', '-247')     ! Output refrac and bangle on standard 247L
           use_247L = .TRUE.
        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('-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

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

  CALL file_delete(ofile, idummy)

!-------------------------------------------------------------------------------
! 4. Set default units
!-------------------------------------------------------------------------------

  CALL ropp_fm_set_units(ro_data)

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

  DO k = 1, n_files

!-------------------------------------------------------------------------------
! 6. 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 )

!-------------------------------------------------------------------------------
! 7. Read data
!-------------------------------------------------------------------------------

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

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

      IF (lfile_exist) THEN

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

        CALL ropp_io_read(l_data, lfile, ranchk=ranchk)

        ro_data%dtocc  = l_data%dtocc
        ro_data%georef = l_data%georef
        ro_data%lev1b  = l_data%lev1b
        ro_data%lev2a  = l_data%lev2a

        WRITE (slev1b, '(i4)') ro_data%lev1b%Npoints
        CALL message(msg_info, "There are " // slev1b // &
                               " level 1b points in " // lfile // " \n")

        WRITE (slev2a, '(i4)') ro_data%lev2a%Npoints
        CALL message(msg_info, "There are " // slev2a // &
                               " level 2a points in " // lfile // " \n")

        CALL ropp_io_free(l_data)

      END IF

      IF (direct_ion) THEN

        ro_data%lev2c%direct_ion = .TRUE.

        CALL ropp_fm_iono_set_default(ro_data)

      END IF

!-------------------------------------------------------------------------------
! 8. Copy data in RO structure to state and refrac obs vectors
!-------------------------------------------------------------------------------

      IF (checksat) state%check_qsat = .TRUE. ! Switch on check against saturation

      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

      IF (ro_data%lev2a%Npoints > 0) THEN
        CALL ropp_fm_roprof2obs(ro_data, obs_refrac)    ! Set to obs levels
        obs_refrac%refrac(:) = ropp_MDFV
      ELSE
        IF (use_247L) THEN
          CALL set_obs_levels_refrac_247(ro_data, obs_refrac, state) ! Pre-defined levels
        ELSE
          CALL set_obs_levels_refrac(ro_data, obs_refrac, &
               geop_min=geop_min, geop_max=geop_max, n_geop=n_geop)  ! Uniformly spaced levels
        END IF
      END IF

!-------------------------------------------------------------------------------
! 9. Calculate refractivity and its gradient
!-------------------------------------------------------------------------------

      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

        IF (do_adj) THEN
          ALLOCATE(gradient_refrac(SIZE(obs_refrac%refrac),SIZE(state%state)))
          CALL ropp_fm_refrac_1d_grad(state, obs_refrac, gradient_refrac)
        END IF

      END IF

!-------------------------------------------------------------------------------
! 10. Copy data in RO and refrac structure to bending angle obs vector
!-------------------------------------------------------------------------------

      IF (ro_data%lev1b%Npoints > 0) THEN
        CALL ropp_fm_roprof2obs(ro_data, obs_bangle)   ! Set to obs levels
        obs_bangle%bangle(:) = ropp_MDFV
      ELSE
        IF (use_247L) THEN
          CALL set_obs_levels_bangle_247(ro_data, obs_bangle) ! Pre-defined levels
        ELSE
          CALL set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle) !default
        END IF
      END IF

!-------------------------------------------------------------------------------
! 11. Calculate bending angle and its gradient
!-------------------------------------------------------------------------------

      IF (state%state_ok) THEN

        CALL ropp_fm_bangle_1d(state, obs_bangle)

        IF (do_adj .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

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

      CALL ropp_fm_obs2roprof(obs_refrac, ro_data)

      CALL ropp_fm_obs2roprof(obs_bangle, ro_data)

      IF (state%state_ok) THEN

        IF (do_adj) THEN
          CALL ropp_io_addvar(ro_data,                                       &
                              name      = "gradient_refrac",                 &
                              long_name =                                    &
                              "Gradient of the refractivity forward model",  &
                              units     = "1",                               &
                              range     = (/MINVAL(gradient_refrac),         &
                                            MAXVAL(gradient_refrac)/),       &
                              DATA      = gradient_refrac)
        END IF

        IF (do_adj .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

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

        CALL ropp_fm_state2roprof(state, ro_data)

      END IF

!-------------------------------------------------------------------------------
! 14. Write data
!-------------------------------------------------------------------------------

      IF (direct_ion) CALL ropp_fm_iono_unpack(ro_data)

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

!-------------------------------------------------------------------------------
! 15. Clean up
!-------------------------------------------------------------------------------

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

      IF (state%state_ok .AND. do_adj .AND. (.NOT. direct_ion)) 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

!-------------------------------------------------------------------------------
! 16a. Calculate uniformly spaced observation levels for refractivity
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_refrac(ro_data, obs_refrac, &
             geop_min, geop_max, n_geop)
!
! n_geop uniformly spaced geopotential heights between geop_min and geop_max (gpm).
!

!   16a.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
    REAL(wp), OPTIONAL               :: geop_min
    REAL(wp), OPTIONAL               :: geop_max
    INTEGER, OPTIONAL                :: n_geop

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

    CALL message_get_routine(routine)
    CALL message_set_routine('set_obs_levels_refrac')

!   16a. Vertical geopotential height levels between geop_min and geop_max
!   ----------------------------------------------------------------------

    zmin = 200.0_wp
    IF ( PRESENT(geop_min) ) zmin = geop_min

    zmax = 60.0e3_wp
    IF ( PRESENT(geop_max) ) zmax = geop_max

    nz = 300
    IF ( PRESENT(n_geop) ) nz = n_geop

    IF ( (zmin >= zmax) .OR. (n_geop <= 1) ) THEN
      CALL message(msg_warn, 'Invalid geop_min, geop_max or n_geop' // &
                             ' ... resetting all to default values')
      zmin =   200.0_wp
      zmax = 60000.0_wp
      nz   =   300
    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+(i-1)*(zmax-zmin)/REAL(nz-1, wp), i=1,nz) /)
    obs_refrac%weights(:) = 1.0_wp
    obs_refrac%lon        = ro_data%GEOref%lon
    obs_refrac%lat        = ro_data%GEOref%lat

    CALL message_set_routine(routine)

  END SUBROUTINE set_obs_levels_refrac


!-------------------------------------------------------------------------------
! 16b. Calculate standard 247 observation levels for refractivity
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_refrac_247(ro_data, obs_refrac, state)
!
! 247 geopotential height levels, going up to about 60 km, chosen in such a
! way that the resulting impact heights match the 'standard' ones in 
! impact_height_eum_247 in ropp_fm/common/ropp_fm_levels.f90
!
!   16b.1 Declarations
!   -----------------

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

    IMPLICIT NONE

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

    INTEGER, PARAMETER                   :: n=247
    REAL(wp), DIMENSION(n)               :: impact_height
    REAL(wp), DIMENSION(n)               :: tmp
    REAL(wp), PARAMETER                  :: refrac_surface_default=300.0_wp
    REAL(wp)                             :: refrac_surface
    REAL(wp), PARAMETER                  :: tol=1.0e-3_wp ! Demand < 1mm diffs
    REAL(wp)                             :: max_refrac
    REAL(wp)                             :: min_geop
    REAL(wp)                             :: roc, und, lat
    INTEGER, PARAMETER                   :: max_iter=1000
    INTEGER                              :: iter, imiss
    INTEGER                              :: lowest_refrac_index

    CHARACTER(LEN=3)                     :: smax_iter
    CHARACTER(LEN=256)                   :: routine

!   16b.2 Initialisations
!   ---------------------

    CALL message_get_routine(routine)
    CALL message_set_routine('set_obs_levels_refrac_247')

    impact_height = impact_height_eum_247  ! The particular 247L set

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

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

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

!   16b.3 Vertical geopotential height levels
!   -----------------------------------------

! We seek solutions h_i of 
! (nr = ) (1 + 1e-6 N(Z(h_i)))*(h_i + RoC + und) = a_i + RoC
! where a_i are the given impact heights, RoC is the radius of curvature 
! and und is the undulation.  h_i is the geometric height; 
! Z_i is the geopotential height.  Non-linear because N depends on Z. 
! When geop_refrac is below the lowest state refrac, the refracs are
! (now) missing data.  In such cases, replace the missing refracs with ones
! calculated with the old refrac scheme with this nullification switched off.
! Because the same principle is applied in set_obs_levels_bangle, the
! resulting impact heights are still good.

!   16b.4 Generate radius of curvature
!   ----------------------------------

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

!   16b.5 Generate undulation
!   -------------------------

    IF ( ro_data%GEOref%undulation > ropp_MDTV ) THEN
      und = ro_data%GEOref%undulation
    ELSE
      CALL message(msg_warn, "Undulation missing ... will assume to " // &
                             "be zero when calculating geopotential heights.")
      und = 0.0_wp
    END IF

!   16b.6 Generate latitude
!   -----------------------

    IF ( ro_data%GEOref%lat > ropp_MDTV ) THEN
      lat = ro_data%GEOref%lat
    ELSE
      CALL message(msg_warn, "Latitude missing ... will assume to " // &
                             "be zero when calculating geopotential heights.")
      lat = 0.0_wp
    END IF

!   16b.7 Initialise
!   ----------------

    IF ( (state%temp(1) > ropp_MDTV) .AND. (state%pres(1) > ropp_MDTV) ) THEN
      refrac_surface = kappa1 * state%pres(1) / state%temp(1) ! roughly 
    ELSE
      refrac_surface = refrac_surface_default
    END IF

    obs_refrac%geop = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*refrac_surface)) - & 
                       (roc + und) ! roughly

    CALL ropp_fm_refrac_1d(state, obs_refrac, nonull=.TRUE.) ! Initialise with the old interp scheme sans missing data

!   16b.8 Iterate on GPH
!   --------------------

    iter = 0

    DO

      tmp = ((impact_height + roc)/(1.0_wp + 1.0e-6_wp*obs_refrac%refrac)) - &
             (roc + und)  ! geometric height

      obs_refrac%geop = geometric2geopotential(lat, tmp)

      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( (/ (i, i=1,n) /), mask=obs_refrac%refrac>ropp_MDTV ) )

        IF ( (lowest_refrac_index > 1) .AND. (lowest_refrac_index < n) ) 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

          tmp(1:lowest_refrac_index-1) = ((impact_height(1:lowest_refrac_index-1) + roc) / &
                                          (1.0_wp + 1.0e-6_wp*obs_refrac%refrac(1:lowest_refrac_index-1))) - &
                                          (roc + und)  ! geometric height exactly right in this region
          obs_refrac%geop(1:lowest_refrac_index-1) = geometric2geopotential(lat, tmp(1:lowest_refrac_index-1))

         END IF

      END IF

      tmp = (1.0_wp + 1.0e-6_wp*obs_refrac%refrac)*(tmp + roc + und) - &
            (impact_height + roc) ! residual

      IF ( MAXVAL(ABS(tmp)) < tol ) EXIT

      iter = iter + 1

      IF ( iter > max_iter ) THEN
        WRITE (smax_iter, '(i3)') max_iter
        CALL message(msg_fatal, 'GPH iteration failed to converge in ' // &
                                 smax_iter // ' iterations')
      END IF

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

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

    END DO

    WRITE (smax_iter, '(i3)') iter
    CALL message(msg_diag, 'GPH iteration converged in ' // &
                                 smax_iter // ' iterations')

!   16b.9 Set other elements of obs_refrac structure
!   ------------------------------------------------

    obs_refrac%refrac     = 0.0_wp ! Refrac will be calculated in the main program
    obs_refrac%weights    = 1.0_wp
    obs_refrac%lon        = ro_data%GEOref%lon
    obs_refrac%lat        = ro_data%GEOref%lat

!   16b.10 Clean up
!   ---------------

    CALL message_set_routine(routine)

  END SUBROUTINE set_obs_levels_refrac_247


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

  SUBROUTINE set_obs_levels_bangle(state, ro_data, obs_refrac, obs_bangle)

!   17.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                             :: n, imiss
    INTEGER                             :: lowest_refrac_index

    CHARACTER(LEN=256)                  :: routine

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

!   17.2 Allocate arrays
!   --------------------

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

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

    n = SIZE(obs_refrac%geop)

    obs_bangle%n_L1 = n

    ALLOCATE(tmp(n))

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

!   17.3 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( (/ (i, i=1,n) /), mask=obs_refrac%refrac>ropp_MDTV ) )

      IF ( (lowest_refrac_index > 1) .AND. (lowest_refrac_index < n) ) 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

!   17.4 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

!   17.5 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:n)     = tmp
      obs_bangle%impact(n+1:2*n) = tmp
    END IF

!   17.6 Fill other arrays
!   ----------------------

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

!   17.7 Clean up
!   -------------

    IF ( state%new_ref_op ) THEN ! Regenerate avec missing data
      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


!-------------------------------------------------------------------------------
! 17b. Calculate standard 247 observation levels for bending angle
!-------------------------------------------------------------------------------

  SUBROUTINE set_obs_levels_bangle_247(ro_data, obs_bangle)

!   17b.1 Declarations
!   ------------------

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

    IMPLICIT NONE

    TYPE(ROprof)                        :: ro_data
    TYPE(Obs1dbangle)                   :: obs_bangle

    INTEGER, PARAMETER                  :: n=247
    REAL(wp), DIMENSION(n)              :: impact_height

    CHARACTER(LEN=256)                  :: routine

    CALL message_get_routine(routine)
    CALL message_set_routine('set_obs_levels_bangle_247')

!   17b.2 Allocate arrays
!   ---------------------

    obs_bangle%n_L1 = n

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

!   17b.3 Define standard impact heights
!   ------------------------------------

    impact_height = impact_height_eum_247  ! The particular 247L set

!   17b.4 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

!   17b.5 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:n)     = impact_height + obs_bangle%r_curve
      obs_bangle%impact(n+1:2*n) = impact_height + obs_bangle%r_curve
    END IF

!   17b.6 Fill other arrays
!   -----------------------

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

!   17b.7 Clean up
!   --------------

    CALL message_set_routine(routine)

  END SUBROUTINE set_obs_levels_bangle_247


!-------------------------------------------------------------------------------
! 18. 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 *, '  -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 *, '  -new_op          use alternative refrac and bangle interpolation'
    PRINT *, '  -direct_ion      model L1 and L2 directly'
    PRINT *, '  -247L            output refrac and bangle on standard 247L'
    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 (gpm)'
    PRINT *, '  -d               output additional diagnostics'
    PRINT *, '  -h               this help'
    PRINT *, '  -v               version information'
    PRINT *, ''
  END SUBROUTINE usage


!-------------------------------------------------------------------------------
! 19. 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