! $Id: ropp_pp_invert_tool.f90 2021 2009-01-16 10:49:04Z frhl $

PROGRAM ropp_pp_invert_tool

!****p* Programs/ropp_pp_invert_tool *
!
! NAME
!    ropp_pp_invert_tool
!
! SYNOPSIS
!   Calculate refractivity and dry temperature profile from L1 and L2
!   channel radio occultation bending angle data using ionospheric
!   correction, statistical optimization, Abel transform, and
!   hydrostatic equation.
!
!   > ropp_pp_invert_tool infile(s)
!                         [-c <config_file>]  [-o <output_file>[
!                         [-m <method>]  [-mfile <file>]  [-bfile <file>]
!                         [-d] [-h] [-v]
!
! ARGUMENTS
!   infile(s)   One or more input file names.
!
! OPTIONS
!   -c <config_file> name of configuration file
!   -o <output_file> name of ROPP netCDF output file
!   -m <method>      ionospheric correction method
!                    [NONE,MSIS,GMSIS,BG,BARO,GBARO], (default GMSIS)
!   -mfile <file>    model coefficients file
!                    [MSIS_coeff.nc,BAROCLIM_coeff.nc], (default MSIS_coeff.nc)
!   -bfile <file>    background model atmospheric profile file path
!                    (if using BG ionospheric correction method)
!   -d               output additional diagnostics
!   -h               help
!   -v               version information
!
!    -o <outfile>  Name of output file (default: ropp_pp_inv.nc).
!
! DESCRIPTION
!    This program reads RO L1 and L2 bending angle data on impact parameter
!    levels from the input data files and calculates vertical profiles of
!    ionospheric corrected bending angle, refractivity, and dry temperature.
!    The result is written to a 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.
!
! ERRORS
!   Program (shell) return codes:
!     0 = OK
!     1 = At least one Warning occurred
!     2 = At least one Error occurred
!     3 = A Fatal error occurred
!
! EXAMPLE
!   To calculate bending angle, refractivity, and dry temperature
!   from one of the example (single-) files in the data directory:
!     > ropp_pp_invert_tool ../data/input.nc -o example_01.nc
!
!   To calculate bending angle, refractivity and dry temperature profiles
!   from all singlefiles in the data directory:
!
!     > ropp_pp_occ_tool ../data/*.nc -o example_02.nc
!
!   Note that the resulting example_02.nc file contains processed data from
!   all example profiles.
!
!   To calculate bending angle, refractivity and dry temperature profiles from
!   all profiles contained in the multifile multi.nc:
!     > ropp_pp_invert_tool ../data/multi.nc -o example_03.nc
!
!   Since the multi_* file was generated by concatenating the other files
!   in the data directory, example_02.nc and example_03.nc should be identical
!   apart from the file names.
!
! 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, L1btype, L2atype
  USE ropp_pp
  USE ropp_pp_types, ONLY: PPConfig, PPDiag

  IMPLICIT NONE

  TYPE(ROprof)   :: ro_data      ! Input RO data
  TYPE(L1btype)  :: Lev1b
  TYPE(L1btype)  :: bangle       ! Input bending angle profiles
  TYPE(L1btype)  :: out_ba       ! Corrected bending angles
  TYPE(L1btype)  :: smt_ba       ! Smoothed bending angle
  TYPE(L1btype)  :: mod_ba       ! Model bending angle

  TYPE(L2atype)  :: refrac       ! Input refractivity profile
  TYPE(L2atype)  :: out_refrac   ! Retrieved refractivity profile
  TYPE(L2btype)  :: dum_meteo    ! Level2b structure for dummy meteo variables

  INTEGER        :: idummy, imin, imax
  INTEGER        :: i, iargc, argc, k
  INTEGER        :: n_files, n_profiles, nbi

  LOGICAL        :: give_help
  LOGICAL        :: ranchk = .TRUE.

  CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
  CHARACTER(len = 4096)                            :: ofile
  CHARACTER(len = 4096)                            :: cfg_file = " "
  CHARACTER(len =  256)                            :: buffer
  CHARACTER(len =  256)                            :: mfile = " "
  CHARACTER(len =  256)                            :: bfile = " "
  CHARACTER(len =   10)                            :: method = " "
  CHARACTER(len =   10)                            :: pstr1, pstr2
  CHARACTER(len =    4)                            :: istr
  CHARACTER(len =    6)                            :: nstr

  REAL(wp) :: Pmax, Pmin
  REAL(wp), DIMENSION(:),   ALLOCATABLE :: mod_refrac   ! climate model refractivity

  TYPE(ppConfig) :: config
  TYPE(ppDiag)   :: diag

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

  give_help = .FALSE.
  ofile     = "ropp_pp_inv.nc"

  CALL message(msg_noin, '')
  CALL message(msg_noin, &
       '----------------------------------------------------------------------')
  CALL message(msg_noin, &
       '                ROPP Pre-processor Invert Tool')
  CALL message(msg_noin, &
       '----------------------------------------------------------------------')
  CALL message(msg_noin, '')

!-------------------------------------------------------------------------------
! 3. Command line arguments
!-------------------------------------------------------------------------------

  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('-c')                          ! Configuration file name
           CALL getarg(i+1, buffer)
           cfg_file = buffer
           i = i + 1
        CASE('-m')                          ! Ionospheric correction method
           CALL getarg(i+1, buffer)
           method = buffer
           i = i + 1
        CASE('-mfile')                      ! Model coefficients file
           CALL getarg(i+1, buffer)
           mfile = buffer
           i = i + 1
        CASE('-bfile')                      ! Background atmosphere profile file
           CALL getarg(i+1, buffer)
           bfile = buffer
           i = i + 1
        CASE ('-no-ranchk')                 ! Use no rangecheck on output
           ranchk = .FALSE.
        CASE('-d')                          ! Additional diagnostic mode
           msg_MODE = VerboseMode
        CASE('-h', '--help', '?')           ! Give help
           give_help = .TRUE.
        CASE('-v', '-V', '--version')       ! Give 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 ) THEN
    CALL message ( msg_error, "No input file(s) specified" )
  END IF

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

  ! 3.1 Read configuration file (if exists), preserving command-line options

  IF (cfg_file /= " ") THEN
    CALL message(msg_info, &
          "Reading configuration file " // TRIM(cfg_file) // ".\n")
    CALL ropp_pp_read_config(cfg_file, config)
  ENDIF

  IF (method /= " ") config%method = method
  IF (mfile  /= " ") config%mfile  = mfile
  IF (bfile  /= " ") config%bfile  = bfile

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

  CALL file_delete(ofile, idummy)

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

  DO k = 1, n_files

    n_profiles = ropp_io_nrec(ifiles(k))

    DO i = 1, n_profiles

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

!-------------------------------------------------------------------------------
! 6. Read data
!-------------------------------------------------------------------------------

        ! 6.1 Read input file

        CALL ropp_io_read(ro_data, ifiles(k), rec=i, ranchk=.TRUE.)
        CALL message(msg_info, "(" // TRIM(ro_data%occ_id) // ") \n")

        IF (ro_data%lev1b%npoints == 0) THEN
          CALL message(msg_fatal, "No Level1b data in file " //    &
             TRIM(ifiles(k)) // ". No data to invert. \n")
        ENDIF

        ! 6.2 Shrink ro_data to correct size (not just for multiple profiles)

        WHERE ( ro_data%lev1b%impact_L1 < ropp_MDTV )
          ro_data%lev1b%impact_L2  = ropp_MDFV
          ro_data%lev1b%impact     = ropp_MDFV
          ro_data%lev1b%impact_opt = ropp_MDFV
        END WHERE

        WHERE ( ro_data%lev1b%bangle_L1 < ropp_MDTV )
          ro_data%lev1b%impact_L1  = ropp_MDFV
          ro_data%lev1b%impact_L2  = ropp_MDFV
          ro_data%lev1b%impact     = ropp_MDFV
          ro_data%lev1b%impact_opt = ropp_MDFV
        END WHERE

        WHERE ( ro_data%lev1b%impact_L2 < ropp_MDTV )
          ro_data%lev1b%impact_L1  = ropp_MDFV
          ro_data%lev1b%impact     = ropp_MDFV
          ro_data%lev1b%impact_opt = ropp_MDFV
        END WHERE

        WHERE ( ro_data%lev1b%bangle_L2 < ropp_MDTV )
          ro_data%lev1b%impact_L1  = ropp_MDFV
          ro_data%lev1b%impact_L2  = ropp_MDFV
          ro_data%lev1b%impact     = ropp_MDFV
          ro_data%lev1b%impact_opt = ropp_MDFV
        END WHERE

        CALL ropp_io_rangecheck ( ro_data )

        CALL ropp_io_init(Lev1b, ro_data%lev1b%npoints)
        Lev1b = ro_data%lev1b
        ro_data%lev1b = Lev1b
        CALL ropp_io_free(Lev1b)

       ! 6.3 Make sure there are no remaining missing values in L1 and L2

        IF (ANY(ro_data%lev1b%impact_L1 <= ropp_MDTV) .OR.  &
            ANY(ro_data%lev1b%impact_L2 <= ropp_MDTV) .OR.  &
            ANY(ro_data%lev1b%bangle_L1 <= ropp_MDTV) .OR.  &
            ANY(ro_data%lev1b%bangle_L2 <= ropp_MDTV)) THEN
          CALL message(msg_error, "Missing L1/L2 impact parameter or " // &
                                  "bending angle values detected.")
          CALL message(msg_info,  "No output generated for " // &
                                  TRIM(ro_data%occ_id) // "\n")
          CALL ropp_io_free(ro_data)
          CYCLE
        ENDIF

!-------------------------------------------------------------------------------
! 7. Copy data in RO structure to observation vectors
!-------------------------------------------------------------------------------

        ! 7.1 Copy valid data by assignment

        bangle = ro_data%lev1b
        refrac = ro_data%lev2a  !never used

        config%r_curve = ro_data%georef%roc
        config%npoints = ro_data%lev1b%npoints

        config%Pmax = MAXVAL(bangle%impact_L1(:))
        config%Pmin = MINVAL(bangle%impact_L1(:))

!-------------------------------------------------------------------------------
! 8. Sort impact parameters to monotonous profile
!-------------------------------------------------------------------------------

        CALL ropp_pp_monotonous(bangle%impact_L1, -1)
        CALL ropp_pp_monotonous(bangle%impact_L2, -1)

!-------------------------------------------------------------------------------
! 9. Interpolate input data to standard grid
!-------------------------------------------------------------------------------

        ! 9.1 Calculate size of standard grid

        IF ( INDEX(config%method, "NONE" ) == 1 ) THEN
           Pmax = MAXVAL(bangle%impact_L1(:)) - 5000.0_wp
        ELSE
           Pmax = config%ztop_invert + config%r_curve
        ENDIF
        Pmin = MINVAL(bangle%impact_L1(:))

        nbi = 1 + CEILING((Pmax - Pmin)/config%dpi)

        WRITE(pstr2, '(f10.3)') Pmax-config%r_curve
        WRITE(pstr1, '(f10.3)') Pmin-config%r_curve
        CALL message(msg_diag, "Pmin = " // pstr1 // " Pmax = " // pstr2)
        WRITE(nstr, '(I6)') nbi
        CALL message(msg_diag, "Standard grid size nbi = " // nstr)

        ! 9.2 Initialise standard grid data structures

        CALL ropp_io_init(out_ba, nbi)
        CALL ropp_io_init(out_refrac, nbi)

        ! 9.3 Interpolate input data onto standard grid

        CALL ropp_pp_merge_profile(bangle%impact_L1, bangle%bangle_L1, &
                                   bangle%impact_L2, bangle%bangle_L2, &
                                   out_ba%impact_L1, out_ba%bangle_L1, &
                                   out_ba%impact_L2, out_ba%bangle_L2, &
                                   Pmin, Pmax)

!-------------------------------------------------------------------------------
! 10. Ionospheric correction of bending angle profile by linear combination
!-------------------------------------------------------------------------------

        IF ( INDEX(config%method, "NONE" ) == 1 ) THEN

          CALL message(msg_info,   &
                    "Retrieving bending angle profile by LINEAR COMBINATION \n")

          CALL ropp_pp_linear_combination(out_ba%impact_L1, out_ba%bangle_L1, &
                                          out_ba%impact_L2, out_ba%bangle_L2, &
                                          out_ba%impact   , out_ba%bangle )
          out_ba%impact_opt = out_ba%impact
          out_ba%bangle_opt = out_ba%bangle

!-------------------------------------------------------------------------------
! 11. Ionospheric correction of bending angles by statistical optimization
!-------------------------------------------------------------------------------

        ELSE

          CALL message(msg_info,   &
              "Retrieving bending angle profile by STATISTICAL OPTIMISATION \n")

! 11.1 Smooth bending angle

          smt_ba = bangle
          CALL ropp_pp_linear_combination(bangle%impact_L1, bangle%bangle_L1, &
                                          bangle%impact_L2, bangle%bangle_L2, &
                                          bangle%impact   , bangle%bangle )
          CALL ropp_pp_smooth_profile(bangle%impact, bangle%bangle,           &
                                      smt_ba%bangle, config)
          smt_ba%impact = bangle%impact

! 11.2 Retrieve model bending angles on L1 impact parameter levels

          CALL ropp_io_init(mod_ba, nbi)

          IF (INDEX(config%method, "GMSIS" ) == 1 .OR. &
              INDEX(config%method, "GBARO" ) == 1) THEN
            IF (INDEX(config%sf_method, "serial" ) == 1) THEN
              CALL ropp_pp_search_model_refraction_new( config%mfile,   &
                                                 ro_data%dtocc%month,   &
                                                 ro_data%georef%lat,    &
                                                 ro_data%georef%lon,    &
                                                 smt_ba%impact,         &
                                                 smt_ba%bangle,         &
                                                 out_ba%impact_L1,      &
                                                 mod_ba%bangle,         &
                                                 config )
            ELSE IF (INDEX(config%sf_method, "parallel" ) == 1) THEN
              CALL ropp_pp_search_model_refraction( config%mfile,   &
                                             ro_data%dtocc%month,   &
                                             ro_data%georef%lat,    &
                                             ro_data%georef%lon,    &
                                             smt_ba%impact,         &
                                             smt_ba%bangle,         &
                                             out_ba%impact_L1,      &
                                             mod_ba%bangle,         &
                                             config )
            ELSE
              CALL message(msg_fatal, "Search and fit method " // &
                                  config%sf_method // " not supported")
              EXIT
            ENDIF
          ELSE IF (INDEX(config%method, "MSIS" ) == 1 .OR. &
                   INDEX(config%method, "BARO" ) == 1) THEN
            CALL ropp_pp_model_refraction(config%mfile, ro_data%dtocc%month, &
                                      ro_data%georef%lat,ro_data%georef%lon, &
                                      out_ba%impact_L1, mod_ba%bangle, config)
          ELSE IF (INDEX(config%method, "BG" ) == 1 ) THEN
            CALL ropp_pp_bg_refraction(config%bfile, ro_data%dtocc%month, &
                                      ro_data%georef%lat,ro_data%georef%lon, &
                                      out_ba%impact_L1, mod_ba%bangle, config)
          ELSE
            CALL message(msg_fatal, "Statistical optimisation method " //     &
                                  config%method // " not supported")
            EXIT
          ENDIF

! 11.3 Fit model bending angle with (smoothed) observed bending angles

          IF (INDEX(config%sf_method, "serial" ) == 1) THEN
            CALL ropp_pp_fit_model_refraction_new(bangle%impact, smt_ba%bangle,    &
                                                  out_ba%impact_L1, mod_ba%bangle, &
                                                  config)
          ELSE IF (INDEX(config%sf_method, "parallel" ) == 1) THEN
            CALL ropp_pp_fit_model_refraction(bangle%impact, smt_ba%bangle,    &
                                              out_ba%impact_L1, mod_ba%bangle, &
                                              config)
          ELSE
            CALL message(msg_fatal, "Search and fit method " // &
                                config%sf_method // " not supported")
            EXIT
          ENDIF

! 11.4 Perform ionospheric correction with statistical optimization
          Imax = SUM(MINLOC(ABS(out_ba%impact_L1(:)-(config%Pmax-5000.0_wp))))
          Imin = SUM(MINLOC(ABS(out_ba%impact_L1(:)-config%Pmin)))

          !! Merge observed bending angles with model above observation top
          WHERE ( out_ba%impact_L1 > out_ba%impact_L1(imax) )
             out_ba%bangle_L1(:) = mod_ba%bangle(:)
             out_ba%bangle_L2(:) = mod_ba%bangle(:)
          END WHERE

          out_ba%impact_opt(:) = out_ba%impact_L1(:)
          out_ba%bangle_opt(:) = out_ba%bangle_L1(:)
          CALL ropp_pp_ionospheric_correction(out_ba%impact_L1(imin:imax),  &
                                              out_ba%bangle_L1(imin:imax),  &
                                              out_ba%impact_L2(imin:imax),  &
                                              out_ba%bangle_L2(imin:imax),  &
                                              out_ba%impact_L1(imin:imax),  &
                                              mod_ba%bangle(imin:imax),     &
                                              config,                       &
                                              out_ba%impact_opt(imin:imax), &
                                              out_ba%bangle_opt(imin:imax), &
                                              diag)

        ENDIF  ! INDEX(config%method, "NONE" ) == 1

! 11.5 Interpolate LC bending angles to output grid

        out_ba%impact = out_ba%impact_opt
        CALL ropp_pp_interpol(bangle%impact, out_ba%impact, &
                              bangle%bangle, out_ba%bangle)

!-------------------------------------------------------------------------------
! 12. Perform inverse Abel transform to compute refractivity
!-------------------------------------------------------------------------------

    ! 12.1 Abel inversion of corrected bending angle profile
        IF ( INDEX(config%method, "NONE" ) == 1 ) THEN

          CALL message(msg_info,   &
                "Retrieving refractivity profile by " // TRIM(config%abel) //  &
                " ABEL TRANSFORM \n")

          IF ( INDEX(config%abel, "LIN" ) == 1 ) THEN
             CALL ropp_pp_invert_LIN(out_ba%impact_opt, out_ba%bangle_opt,    &
                                     out_ba%impact_opt, out_refrac%refrac)
          ELSE IF ( INDEX(config%abel, "EXP" ) == 1 ) THEN
             CALL ropp_pp_invert_EXP(out_ba%impact_opt, out_ba%bangle_opt,    &
                                     out_ba%impact_opt, out_refrac%refrac)
          ELSE
             CALL message(msg_fatal, "Abel integral method " // config%abel // &
                                     " not supported")
          ENDIF

        out_refrac%alt_refrac = (out_ba%impact_opt /                           &
                                (1.0_wp + out_refrac%refrac * 1.e-6_wp))   &
                                - (ro_data%GEOref%roc+ro_data%GEOref%undulation)
        out_refrac%geop_refrac =    &
              geometric2geopotential(ro_data%georef%lat, out_refrac%alt_refrac)

    ! 12.2 Abel inversion of corrected bending angle profile with
    !      statistical optimization

        ELSE
  
          ALLOCATE(mod_refrac(nbi) )

          CALL message(msg_info,   &
              "Retrieving refractivity profile by " // TRIM(config%abel) //  &
              " ABEL TRANSFORM with STAT OPT \n")

          CALL ropp_pp_invert_refraction(config%mfile,             &
                                         ro_data%dtocc%month,      &
                                         ro_data%georef%lat,       &
                                         ro_data%georef%lon,       &
                                         out_ba%impact_opt,        &
                                         out_ba%bangle_opt,        &
                                         out_refrac%geop_refrac,   &
                                         out_refrac%refrac,        &
                                         mod_refrac,              &
                                         config)
          DEALLOCATE(mod_refrac)

          IF (ro_data%georef%undulation > ropp_MDTV) THEN
            out_refrac%alt_refrac =                                           &
               ((out_ba%impact_opt / (1.0_wp + out_refrac%refrac * 1.e-6_wp)) &
                 - (ro_data%GEOref%roc + ro_data%GEOref%undulation))
          ELSE
            CALL message(msg_warn, "Invalid undulation calculated. " //      &
               "Check for valid EGM geoid coefficient and correction file." // &
               "\n    Writing output altitude scales with respect to ELLIPSOID.")
            out_refrac%alt_refrac =                                         &
               ((out_ba%impact_opt / (1.0_wp + out_refrac%refrac * 1.e-6_wp))   &
               - ro_data%GEOref%roc)
          ENDIF

          out_refrac%geop_refrac =     &
               geometric2geopotential(ro_data%georef%lat, out_refrac%alt_refrac)

        ENDIF

!-------------------------------------------------------------------------------
! 13. Compute Tdry and Pdry - though only Tdry is written to ro_data
!-------------------------------------------------------------------------------

        CALL message(msg_info, "Computing dry temperature \n")
        CALL ropp_io_init(dum_meteo, nbi) ! easy way to get dummy storage
        dum_meteo%shum = 0.0_wp
        CALL ropp_pp_tdry(ro_data%georef%lat, out_refrac%alt_refrac,   &
             out_refrac%refrac, dum_meteo%shum, out_refrac%dry_temp,   &
             dum_meteo%press)

!-------------------------------------------------------------------------------
! 14. Copy retrieved profiles to RO structure 
!     (only output data within observed range and according to output flags)
!-------------------------------------------------------------------------------

        IF (config%output_lev1b) THEN
          out_ba%npoints = imax - imin + 1
          CALL ropp_io_roprof2roprof(out_ba, ro_data%lev1b)
        ELSE
          CALL ropp_io_free(ro_data%lev1b)
        END IF

        IF (config%output_lev2a) THEN
          out_refrac%npoints = imax - imin + 1
          CALL ropp_io_roprof2roprof(out_refrac, ro_data%lev2a)
        ELSE
          CALL ropp_io_free(ro_data%lev2a)
        END IF

        IF (.NOT.config%output_lev1a) CALL ropp_io_free(ro_data%lev1a)

        CALL ropp_io_free(ro_data%vlist)   !! interim, avoid writing lcf data
        IF (config%output_diag) THEN
          CALL message(msg_info, "Writing additional diagnostic output \n")
          CALL ropp_pp_diag2roprof(diag, ro_data)
        ENDIF

!-------------------------------------------------------------------------------
! 15. Write data
!-------------------------------------------------------------------------------

        CALL message(msg_info, "Writing to output file " // TRIM(ofile) // "\n")
        CALL ropp_io_write(ro_data, ofile, append=.TRUE., ranchk=ranchk )

!-------------------------------------------------------------------------------
! 16. Clean up
!-------------------------------------------------------------------------------

        CALL ropp_io_free(ro_data)
        CALL ropp_io_free(bangle)
        CALL ropp_io_free(refrac)
        CALL ropp_io_free(out_ba)
        CALL ropp_io_free(out_refrac)
        CALL ropp_io_free(dum_meteo)

        IF (ASSOCIATED(diag%ba_ion)) DEALLOCATE(diag%ba_ion)
        IF (ASSOCIATED(diag%err_ion)) DEALLOCATE(diag%err_ion)
        IF (ASSOCIATED(diag%err_neut)) DEALLOCATE(diag%err_neut)
        IF (ASSOCIATED(diag%wt_data)) DEALLOCATE(diag%wt_data)
        IF (ASSOCIATED(diag%ba_model)) DEALLOCATE(diag%ba_model)

    END DO  ! profiles

  END DO  ! files

  CALL EXIT(msg_exit_status)

  CONTAINS

!-------------------------------------------------------------------------------
! 17. Usage information
!-------------------------------------------------------------------------------

  SUBROUTINE usage()
    PRINT *, 'Purpose:'
    PRINT *, '  Calculate corrected bending angle, refractivity'
    PRINT *, '  and dry temperature from L1, L2 bending angle data'
    PRINT *, 'Usage:'
    PRINT *, '  > ropp_pp_invert_tool [<options>] <input_file(s)>'
    PRINT *, 'Options:'
    PRINT *, '  -c <config_file> name of configuration file'
    PRINT *, '  -o <output_file> name of ROPP netCDF output file'
    PRINT *, '  -m <method>      ionospheric correction method'
    PRINT *, '                    [NONE,MSIS,GMSIS,BG,BARO,GBARO], (default GMSIS)'
    PRINT *, '  -mfile <file>    model coefficients file'
    PRINT *, '                   [MSIS_coeff.nc,BAROCLIM_coeff.nc], (default MSIS_coeff.nc)'
    PRINT *, '  -bfile <file>    background model atmospheric profile file path'
    PRINT *, '                    (if using BG ionospheric correction method)'
    PRINT *, '  -d               output additional diagnostics'
    PRINT *, '  -h               this help'
    PRINT *, '  -v               version information'
    PRINT *, ''
  END SUBROUTINE usage

!-------------------------------------------------------------------------------
! 18. Version information
!-------------------------------------------------------------------------------

  SUBROUTINE version_info()
    CHARACTER (LEN=40) :: version
    version = ropp_pp_version()
    PRINT *, 'ropp_pp_invert_tool - Pre-processor tool:'
    PRINT *, 'Calculate corrected bending angles, refractivity, '
    PRINT *, 'and dry temperature.'
    PRINT *, ''
    PRINT *, 'This program is part of ROPP (PP) Release ' // TRIM(version)
    PRINT *, ''
  END SUBROUTINE version_info

END PROGRAM ropp_pp_invert_tool