! $Id: ropp_apps_tph_tool.f90 3491 2013-02-06 12:43:43Z idculv $

PROGRAM ropp_apps_tph_tool

!****p* Programs/ropp_apps_tph_tool *
!
! NAME
!   ropp_apps_tph_tool
!
! SYNOPSIS
!   Tropopause height (TPH) diagnostic
!
!   > ropp_apps_tph_tool [<options>] <infile(s)>
!
! ARGUMENTS
!   <infile(s)>   One (or more) input file names.
!
! OPTIONS
!   -o <output_file>  name of ROPP netCDF output file
!   -b                calculate bending-angle-based TPH
!   -n                calculate refractivity-based TPH
!   -y                calculate dry-temperature-based TPH
!   -t                calculate temperature-based TPH
!   --no-ranchk       no range check performed on input
!                     (output is never range checked)
!   -h                help
!   -d                output additional diagnostics
!   -v                version information
!
! DESCRIPTION
!   Diagnose tropopause height from the kinks in one or more RO profiles,
!   using the covariance transform method of Lewis (GRL 2009).
!
! NOTES
!   If the input file is a multifile, or more than one input files are
!   specified, the output file is a multifile.
!
!   Existing output files will be overwritten.
!
!   If none of the {-b, -n, -y, -t} options is specified, the tool will
!   attempt to calculate all four TPHs.
!
! ERRORS
!   Program (shell) return codes:
!     0 = OK
!     1 = At least one Warning occurred
!     2 = At least one Error occurred
!     3 = A Fatal error occurred
!
! 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_apps_types
  USE ropp_apps

  IMPLICIT NONE

  TYPE(ROprof)                                     :: ro_data

! Local variables

  REAL(wp)                                         :: und

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

  LOGICAL                                          :: give_help
  LOGICAL                                          :: l_bangle, l_refrac, &
                                                      l_tdry, l_temp, l_all_tph
  LOGICAL                                          :: l_diag
  LOGICAL                                          :: ranchk = .TRUE.

  CHARACTER(len = 4096), DIMENSION(:), ALLOCATABLE :: ifiles
  CHARACTER(len = 4096)                            :: ofile
  CHARACTER(len =  256)                            :: buffer
  CHARACTER(len =    4)                            :: istr
  CHARACTER(len =    6)                            :: nstr
  CHARACTER(len =   10)                            :: tph_str
  CHARACTER(len =   10)                            :: tpt_str
  CHARACTER(len = 4096)                            :: flag_meaning

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

  CALL message_set_routine('ropp_apps_tph_tool')

  give_help = .FALSE.

  l_bangle  = .FALSE.
  l_refrac  = .FALSE.
  l_tdry    = .FALSE.
  l_temp    = .FALSE.
  l_all_tph = .FALSE.

  l_diag    = .FALSE.

  ofile     = "ropp_apps_tph.nc"

  CALL message(msg_noin, '')
  CALL message(msg_noin, &
       '----------------------------------------------------------------------')
  CALL message(msg_noin, &
       '             ROPP Applications Tropopause Height 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('-b')                          ! Calculate bending-angle-based TPHs
         l_bangle  = .TRUE.
      CASE('-n')                          ! Calculate refractivity-based TPHs
         l_refrac  = .TRUE.
      CASE('-y')                          ! Calculate dry-temperature-based TPHs
         l_tdry    = .TRUE.
      CASE('-t')                          ! Calculate temperature-based TPHs
         l_temp    = .TRUE.
      CASE ('--no-ranchk')                ! Don't rangecheck the input
         ranchk = .FALSE.
      CASE('-h', '--help', '?')           ! Give some help
         give_help = .TRUE.
      CASE('-d')                          ! Output more diagnostic information
         msg_MODE = VerboseMode
         l_diag = .TRUE.
      CASE('-v', '-V', '--version')       ! Output version info
         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)
  END IF

  IF (.NOT. (l_bangle .OR. l_refrac .OR. l_tdry .OR. l_temp)) l_all_tph = .TRUE.

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

  CALL file_delete(ofile, idummy)

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

  DO k=1,n_files

    CALL message(msg_info, "Processing file " // TRIM(ADJUSTL(ifiles(k))))

!-------------------------------------------------------------------------------
! 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_info, "Processing profile " // istr // " of " // nstr )

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

! 7.1 Read data from file into an ROprof structure.

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

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

! 7.2 Make sure the ROprof structure includes a lev2c sub-structure.

      IF ( (ro_data%lev2c%geop_sfc  < ropp_MDTV) .AND. &
           (ro_data%lev2c%press_sfc < ropp_MDTV) ) THEN
        CALL ropp_io_init(ro_data%lev2c, 1)
      END IF

      ro_data%lev2c%npoints = 1
      ro_data%lev2c%missing = .FALSE.

! 7.3 Delete the extra data if the TPH routines would try to generate them

      IF ( l_diag ) CALL ropp_io_free_vlistD1d(ro_data%vlist%vlistD1d)

! 7.4 Ensure that profile's heights are increasing - 1st element nearest the surface

      CALL ropp_io_ascend(ro_data)

!-------------------------------------------------------------------------------
! 8. Calculate TPH based on bending angle, if possible
!-------------------------------------------------------------------------------

      IF (l_bangle .OR. l_all_tph) THEN

        ro_data%lev2c%tph_bangle      = ropp_MDFV
        ro_data%lev2c%tpa_bangle      = ropp_MDFV
        ro_data%lev2c%tph_bangle_flag = ropp_MIFV

        CALL ropp_apps_tph_bangle(ro_data, diag=l_diag)

        CALL tph_flag_decode(ro_data%lev2c%tph_bangle_flag, flag_meaning)

        CALL message(msg_diag, "Bending angle-based TPH QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%tph_bangle > ropp_MDTV) THEN

          IF (ro_data%GEOref%undulation < ropp_MDTV) THEN
            und = 0.0_wp
          ELSE
            und = ro_data%GEOref%undulation
          END IF

          WRITE (tph_str, FMT='(F10.5)') &
            (ro_data%lev2c%tph_bangle - ro_data%GEOref%roc - und)*1.0e-3_wp
          WRITE (tpt_str, FMT='(F10.5)') ro_data%lev2c%tpa_bangle*1.0e3_wp
          CALL message(msg_info, "Bending angle-based TPH: " // &
                                  tpt_str // " mrad at " // tph_str // " km \n")

        END IF

      END IF

!-------------------------------------------------------------------------------
! 9. Calculate TPH based on refractivity, if possible
!-------------------------------------------------------------------------------

      IF (l_refrac .OR. l_all_tph) THEN

        ro_data%lev2c%tph_refrac      = ropp_MDFV
        ro_data%lev2c%tpn_refrac      = ropp_MDFV
        ro_data%lev2c%tph_refrac_flag = ropp_MIFV

        CALL ropp_apps_tph_refrac(ro_data, diag=l_diag)

        CALL tph_flag_decode(ro_data%lev2c%tph_refrac_flag, flag_meaning)

        CALL message(msg_diag, "Refractivity-based TPH QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%tph_refrac > ropp_MDTV) THEN

          WRITE (tph_str, FMT='(F10.5)') ro_data%lev2c%tph_refrac*1.0e-3_wp
          WRITE (tpt_str, FMT='(F10.5)') ro_data%lev2c%tpn_refrac
          CALL message(msg_info, "Refractivity-based TPH: " // &
                                  tpt_str // " N-unit at " // tph_str // " km \n")

        END IF

      END IF

!-------------------------------------------------------------------------------
! 10. Calculate TPH based on dry temperature, if possible
!-------------------------------------------------------------------------------

      IF (l_tdry .OR. l_all_tph) THEN

        ro_data%lev2c%tph_tdry_lrt      = ropp_MDFV
        ro_data%lev2c%tpt_tdry_lrt      = ropp_MDFV
        ro_data%lev2c%tph_tdry_lrt_flag = ropp_MIFV

        ro_data%lev2c%tph_tdry_cpt      = ropp_MDFV
        ro_data%lev2c%tpt_tdry_cpt      = ropp_MDFV
        ro_data%lev2c%tph_tdry_cpt_flag = ropp_MIFV

        ro_data%lev2c%prh_tdry_cpt      = ropp_MDFV
        ro_data%lev2c%prt_tdry_cpt      = ropp_MDFV
        ro_data%lev2c%prh_tdry_cpt_flag = ropp_MIFV

        CALL ropp_apps_tph_tdry(ro_data, diag=l_diag)

! 10.1 Lapse rate

        CALL tph_flag_decode(ro_data%lev2c%tph_tdry_lrt_flag, flag_meaning)

        CALL message(msg_diag, "Dry temperature lapse rate-based TPH QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%tph_tdry_lrt > ropp_MDTV) THEN

          WRITE (tph_str, FMT='(F10.5)') ro_data%lev2c%tph_tdry_lrt*1.0e-3_wp
          WRITE (tpt_str, FMT='(F10.5)') ro_data%lev2c%tpt_tdry_lrt
          CALL message(msg_info, "Dry temperature lapse rate-based TPH: " // &
                                  tpt_str // " K at " // tph_str // " km \n")

        END IF

! 10.2 Cold point

        CALL tph_flag_decode(ro_data%lev2c%tph_tdry_cpt_flag, flag_meaning)

        CALL message(msg_diag, "Dry temperature cold point-based TPH QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%tph_tdry_cpt > ropp_MDTV) THEN

          WRITE (tph_str, FMT='(F10.5)') ro_data%lev2c%tph_tdry_cpt*1.0e-3_wp
          WRITE (tpt_str, FMT='(F10.5)') ro_data%lev2c%tpt_tdry_cpt
          CALL message(msg_info, "Dry temperature cold point-based TPH: " // &
                                  tpt_str // " K at " // tph_str // " km \n")

        END IF

! 10.3 Profile Tmin

        CALL tph_flag_decode(ro_data%lev2c%prh_tdry_cpt_flag, flag_meaning)

        CALL message(msg_diag, "Entire profile dry temperature cold point QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%prh_tdry_cpt > ropp_MDTV) THEN

          WRITE (tph_str, FMT='(F10.5)') ro_data%lev2c%prh_tdry_cpt*1.0e-3_wp
          WRITE (tpt_str, FMT='(F10.5)') ro_data%lev2c%prt_tdry_cpt
          CALL message(msg_diag, "Cold point of entire dry temperature profile: " // &
                                  tpt_str // " K at " // tph_str // " km \n")

        END IF

      END IF

!-------------------------------------------------------------------------------
! 11. Calculate TPH based on temperature, if possible
!-------------------------------------------------------------------------------

      IF (l_temp .OR. l_all_tph) THEN

        ro_data%lev2c%tph_temp_lrt      = ropp_MDFV
        ro_data%lev2c%tpt_temp_lrt      = ropp_MDFV
        ro_data%lev2c%tph_temp_lrt_flag = ropp_MIFV

        ro_data%lev2c%tph_temp_cpt      = ropp_MDFV
        ro_data%lev2c%tpt_temp_cpt      = ropp_MDFV
        ro_data%lev2c%tph_temp_cpt_flag = ropp_MIFV

        ro_data%lev2c%prh_temp_cpt      = ropp_MDFV
        ro_data%lev2c%prt_temp_cpt      = ropp_MDFV
        ro_data%lev2c%prh_temp_cpt_flag = ropp_MIFV

        CALL ropp_apps_tph_temp(ro_data, diag=l_diag)

! 11.1 Lapse rate

        CALL tph_flag_decode(ro_data%lev2c%tph_temp_lrt_flag, flag_meaning)

        CALL message(msg_diag, "Temperature lapse rate-based TPH QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%tph_temp_lrt > ropp_MDTV) THEN

          WRITE (tph_str, FMT='(F10.5)') ro_data%lev2c%tph_temp_lrt*1.0e-3_wp
          WRITE (tpt_str, FMT='(F10.5)') ro_data%lev2c%tpt_temp_lrt
          CALL message(msg_info, "Temperature lapse rate-based TPH: " // &
                                  tpt_str // " K at " // tph_str // " km \n")

        END IF

! 11.2 Cold point

        CALL tph_flag_decode(ro_data%lev2c%tph_temp_cpt_flag, flag_meaning)

        CALL message(msg_diag, "Temperature cold point-based TPH QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%tph_temp_cpt > ropp_MDTV) THEN

          WRITE (tph_str, FMT='(F10.5)') ro_data%lev2c%tph_temp_cpt*1.0e-3_wp
          WRITE (tpt_str, FMT='(F10.5)') ro_data%lev2c%tpt_temp_cpt
          CALL message(msg_info, "Temperature cold point-based TPH: " // &
                                  tpt_str // " K at " // tph_str // " km \n")

        END IF

! 11.3 Profile Tmin

        CALL tph_flag_decode(ro_data%lev2c%prh_temp_cpt_flag, flag_meaning)

        CALL message(msg_diag, "Entire profile temperature cold point QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%prh_temp_cpt > ropp_MDTV) THEN

          WRITE (tph_str, FMT='(F10.5)') ro_data%lev2c%prh_temp_cpt*1.0e-3_wp
          WRITE (tpt_str, FMT='(F10.5)') ro_data%lev2c%prt_temp_cpt
          CALL message(msg_diag, "Cold point of entire temperature profile: " // &
                                 tpt_str // " K at " // tph_str // " km \n")

        END IF

      END IF

!-------------------------------------------------------------------------------
! 12. Write data
!-------------------------------------------------------------------------------

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

!-------------------------------------------------------------------------------
! 13. Clean up
!-------------------------------------------------------------------------------

      CALL ropp_io_free(ro_data)

    END DO  ! loop over profiles

  END DO  ! loop over files

  CALL EXIT(msg_exit_status)

CONTAINS

!-------------------------------------------------------------------------------
! 14. TPH QC flag decoding
!-------------------------------------------------------------------------------

  SUBROUTINE tph_flag_decode(tph_flag, meaning)

    INTEGER                                        :: tph_flag
    CHARACTER (LEN=5)                              :: flag_str
    CHARACTER (LEN=4096)                           :: meaning

    WRITE (flag_str, '(i5)') tph_flag

    meaning = flag_str // " ==> "

    IF (tph_flag < ropp_MITV) THEN
      meaning = TRIM(meaning) // " Initial value or incalculable ... \n"
      RETURN
    END IF

    IF (tph_flag == 0) THEN
      meaning = TRIM(meaning) // " Diagnosed value OK ..."
      RETURN
    END IF

    IF (BTEST(tph_flag, TPH_QC_data_invalid)) &
      meaning = TRIM(meaning) // " Invalid profile data ... \n"

    IF (BTEST(tph_flag, TPH_QC_prof_depth)) &
      meaning = TRIM(meaning) // " Profile not deep enough ... \n"

    IF (BTEST(tph_flag, TPH_QC_prof_height)) &
      meaning = TRIM(meaning) // " Profile not high enough ... \n"

    IF (BTEST(tph_flag, TPH_QC_CT_smooth_above)) &
      meaning = TRIM(meaning) // " Cov transform too smooth above TPH ..."

    IF (BTEST(tph_flag, TPH_QC_CT_smooth_below)) &
      meaning = TRIM(meaning) // " Cov transform too smooth below TPH ..."

    IF (BTEST(tph_flag, TPH_QC_double_trop)) &
      meaning = TRIM(meaning) // " Double tropopause detected ..."

    IF (BTEST(tph_flag, TPH_QC_too_low)) &
      meaning = TRIM(meaning) // " Diagnosed TPH too low ..."

    IF (BTEST(tph_flag, TPH_QC_too_high)) &
      meaning = TRIM(meaning) // " Diagnosed TPH too high ..."

  END SUBROUTINE tph_flag_decode

!-------------------------------------------------------------------------------
! 15. Usage information
!-------------------------------------------------------------------------------

  SUBROUTINE usage()
    PRINT *, 'Purpose:'
    PRINT *, '  Tropopause height diagnostic'
    PRINT *, 'Usage:'
    PRINT *, '  > ropp_apps_tph_tool [<options>] <input_file(s)>'
    PRINT *, 'Options:'
    PRINT *, '  -o <output_file>  name of ROPP netCDF output file'
    PRINT *, '  -b                calculate bending-angle-based TPH'
    PRINT *, '  -n                calculate refractivity-based TPH'
    PRINT *, '  -y                calculate dry-temperature-based TPH'
    PRINT *, '  -t                calculate temperature-based TPH'
    PRINT *, '  -h                this help'
    PRINT *, '  -d                output additional diagnostics'
    PRINT *, '  -v                version information'
    PRINT *, ''
  END SUBROUTINE usage

!-------------------------------------------------------------------------------
! 16. Version information
!-------------------------------------------------------------------------------

  SUBROUTINE version_info()
    CHARACTER (LEN=40) :: version
    version = ropp_apps_version()
    PRINT *, 'ropp_apps_tph_tool - Tropopause height diagnostic.'
    PRINT *, ''
    PRINT *, 'This program is part of ROPP (PP) Release ' // TRIM(version)
    PRINT *, ''
  END SUBROUTINE version_info

END PROGRAM ropp_apps_tph_tool