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

PROGRAM ropp_apps_pblh_tool

!****p* Programs/ropp_apps_pblh_tool *
!
! NAME
!   ropp_apps_pblh_tool
!
! SYNOPSIS
!   Planetary boundary layer height (PBLH) diagnostic
!
!   > ropp_apps_pblh_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 PBLH
!   -n                calculate refractivity-based PBLH
!   -y                calculate dry-temperature-based PBLH
!   -t                calculate temperature-based PBLH
!   -q                calculate specific-humidity-based PBLH
!   -r                calculate relative-humidity-based PBLH
!   -h                help
!   -d                output additional diagnostics
!   -v                version information
!
! DESCRIPTION
!   Diagnose planetary boundary layer 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 PBLHs.
!
! 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 ncdf

  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, ro_data1

! Local variables

  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_shum, l_rhum, &
                                                      l_all_pblh
  LOGICAL                                          :: l_diag

  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)                            :: pblh_str, pblx_str
  CHARACTER(len = 4096)                            :: flag_meaning

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

  CALL message_set_routine('ropp_apps_pblh_tool')

  give_help  = .FALSE.

  l_bangle   = .FALSE.
  l_refrac   = .FALSE.
  l_tdry     = .FALSE.
  l_temp     = .FALSE.
  l_shum     = .FALSE.
  l_rhum     = .FALSE.
  l_all_pblh = .FALSE.

  l_diag     = .FALSE.

  ofile      = "ropp_apps_pblh.nc"

  CALL message(msg_noin, '')
  CALL message(msg_noin, &
       '----------------------------------------------------------------------')
  CALL message(msg_noin, &
       '        ROPP Applications Planetary Boundary Layer 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 PBLHs
         l_bangle  = .TRUE.
      CASE('-n')                          ! Calculate refractivity-based PBLHs
         l_refrac  = .TRUE.
      CASE('-y')                          ! Calculate dry-temperature-based PBLHs
         l_tdry    = .TRUE.
      CASE('-t')                          ! Calculate temperature-based PBLHs
         l_temp    = .TRUE.
      CASE('-q')                          ! Calculate specific-humidity-based PBLHs
         l_shum    = .TRUE.
      CASE('-r')                          ! Calculate relative-humidity-based PBLHs
         l_rhum    = .TRUE.
      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   .OR. l_shum   .OR. l_rhum)) l_all_pblh = .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
!-------------------------------------------------------------------------------

      CALL ropp_io_read(ro_data1,  ifiles(k), rec=i, ranchk=.TRUE.)

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

! 7.1 Generate an Roprof structure
! (a) includes a lev2c sub-structure (to hold the PBLHs) if necessary,
! (b) contains (eg) a lev2a substructure even if this is full of missing data originally, and
! (c) is absent of 'extra_data', which prevent variables of the same name being added later

      CALL ropp_io_init(ro_data, &
                        ro_data1%lev1a%npoints, &
                        ro_data1%lev1b%npoints, &
                        ro_data1%lev2a%npoints, &
                        ro_data1%lev2b%npoints, &
                                             1, &
                        ro_data1%lev2d%npoints)

      ro_data = ro_data1

      IF ( ro_data1%lev2c%npoints == 0 ) CALL ropp_io_init(ro_data%lev2c, 1)

      IF ( ro_data1%lev2c%geop_sfc > ropp_MDTV ) THEN 
        ro_data%lev2c%geop_sfc = ro_data1%lev2c%geop_sfc
      ELSE
        ro_data%lev2c%geop_sfc = 0.0_wp  ! To stop the range-checking zapping the lev2c component
      END IF

      IF ( ro_data1%lev2c%press_sfc > ropp_MDTV ) THEN 
        ro_data%lev2c%press_sfc = ro_data1%lev2c%press_sfc
      ELSE
        ro_data%lev2c%press_sfc = 1000.0_wp  ! To stop the range-checking zapping the lev2c component
      END IF

      ro_data%lev2d%level_type = ro_data1%lev2d%level_type  ! Not yet done in ropp_io_assign.f90

      CALL ropp_io_free(ro_data1)

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

      CALL ropp_io_ascend(ro_data)

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

      IF (l_bangle .OR. l_all_pblh) THEN

        ro_data%lev2c%pblh_bangle      = ropp_MDFV
        ro_data%lev2c%pblh_bangle2     = ropp_MDFV
        ro_data%lev2c%pbla_bangle      = ropp_MDFV
        ro_data%lev2c%pbla_bangle2     = ropp_MDFV
        ro_data%lev2c%pblh_bangle_flag = ropp_MIFV

        CALL ropp_apps_pblh_bangle(ro_data, diag=l_diag)

        CALL pblh_flag_decode(ro_data%lev2c%pblh_bangle_flag, flag_meaning)

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

        IF (ro_data%lev2c%pblh_bangle > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_bangle*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pbla_bangle*1.0e3_wp
          CALL message(msg_info, "Bending-angle-based PBLH: " // &
                                  pblx_str // " mrad at " // pblh_str // " km \n")

        END IF

        IF (ro_data%lev2c%pblh_bangle2 > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_bangle2*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pbla_bangle2*1.0e3_wp
          CALL message(msg_info, "2nd Bending-angle-based PBLH: " // &
                                  pblx_str // " mrad at " // pblh_str // " km \n")

        END IF

      END IF

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

      IF (l_refrac .OR. l_all_pblh) THEN

        ro_data%lev2c%pblh_refrac      = ropp_MDFV
        ro_data%lev2c%pblh_refrac2     = ropp_MDFV
        ro_data%lev2c%pbln_refrac      = ropp_MDFV
        ro_data%lev2c%pbln_refrac2     = ropp_MDFV
        ro_data%lev2c%pblh_refrac_flag = ropp_MIFV

        CALL ropp_apps_pblh_refrac(ro_data, diag=l_diag)

        CALL pblh_flag_decode(ro_data%lev2c%pblh_refrac_flag, flag_meaning)

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

        IF (ro_data%lev2c%pblh_refrac > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_refrac*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pbln_refrac
          CALL message(msg_info, "Refractivity-based PBLH: " // &
                                  pblx_str // " N-unit at " // pblh_str // " km \n")

        END IF

        IF (ro_data%lev2c%pblh_refrac2 > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_refrac2*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pbln_refrac2
          CALL message(msg_info, "2nd Refractivity-based PBLH: " // &
                                  pblx_str // " N-unit at " // pblh_str // " km \n")

        END IF

      END IF

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

      IF (l_tdry .OR. l_all_pblh) THEN

        ro_data%lev2c%pblh_tdry      = ropp_MDFV
        ro_data%lev2c%pblh_tdry2     = ropp_MDFV
        ro_data%lev2c%pblt_tdry      = ropp_MDFV
        ro_data%lev2c%pblt_tdry2     = ropp_MDFV
        ro_data%lev2c%pblh_tdry_flag = ropp_MIFV

        CALL ropp_apps_pblh_tdry(ro_data, diag=l_diag)

        CALL pblh_flag_decode(ro_data%lev2c%pblh_tdry_flag, flag_meaning)

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

        IF (ro_data%lev2c%pblh_tdry > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_tdry*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pblt_tdry
          CALL message(msg_info, "Dry-temperature-based PBLH: " // &
                                  pblx_str // " K at " // pblh_str // " km \n")

        END IF

        IF (ro_data%lev2c%pblh_tdry2 > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_tdry2*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pblt_tdry2
          CALL message(msg_info, "2nd Dry-temperature-based PBLH: " // &
                                  pblx_str // " K at " // pblh_str // " km \n")

        END IF

      END IF

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

      IF (l_temp .OR. l_all_pblh) THEN

        ro_data%lev2c%pblh_temp      = ropp_MDFV
        ro_data%lev2c%pblh_temp2     = ropp_MDFV
        ro_data%lev2c%pblt_temp      = ropp_MDFV
        ro_data%lev2c%pblt_temp2     = ropp_MDFV
        ro_data%lev2c%pblh_temp_flag = ropp_MIFV

        CALL ropp_apps_pblh_temp(ro_data, diag=l_diag)

        CALL pblh_flag_decode(ro_data%lev2c%pblh_temp_flag, flag_meaning)

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

        IF (ro_data%lev2c%pblh_temp > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_temp*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pblt_temp
          CALL message(msg_info, "Temperature-based PBLH: " // &
                                  pblx_str // " K at " // pblh_str // " km \n")

        END IF

        IF (ro_data%lev2c%pblh_temp2 > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_temp2*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pblt_temp2
          CALL message(msg_info, "2nd Temperature-based PBLH: " // &
                                  pblx_str // " K at " // pblh_str // " km \n")

        END IF

      END IF

!-------------------------------------------------------------------------------
! 12. Calculate PBLH based on specific humidity, if possible
!-------------------------------------------------------------------------------

      IF (l_shum .OR. l_all_pblh) THEN

        ro_data%lev2c%pblh_shum      = ropp_MDFV
        ro_data%lev2c%pblh_shum2     = ropp_MDFV
        ro_data%lev2c%pblq_shum      = ropp_MDFV
        ro_data%lev2c%pblq_shum2     = ropp_MDFV
        ro_data%lev2c%pblh_shum_flag = ropp_MIFV

        CALL ropp_apps_pblh_shum(ro_data, diag=l_diag)

        CALL pblh_flag_decode(ro_data%lev2c%pblh_shum_flag, flag_meaning)

        CALL message(msg_diag, "Specific-humidity-based PBLH QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%pblh_shum > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_shum*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pblq_shum
          CALL message(msg_info, "Specific-humidity-based PBLH: " // &
                                  pblx_str // " g/kg at " // pblh_str // " km \n")

        END IF

        IF (ro_data%lev2c%pblh_shum2 > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_shum2*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pblq_shum2
          CALL message(msg_info, "2nd Specific-humidity-based PBLH: " // &
                                  pblx_str // " g/kg at " // pblh_str // " km \n")

        END IF

      END IF

!-------------------------------------------------------------------------------
! 13. Calculate PBLH based on relative humidity, if possible
!-------------------------------------------------------------------------------

      IF (l_rhum .OR. l_all_pblh) THEN

        ro_data%lev2c%pblh_rhum      = ropp_MDFV
        ro_data%lev2c%pblh_rhum2     = ropp_MDFV
        ro_data%lev2c%pblr_rhum      = ropp_MDFV
        ro_data%lev2c%pblr_rhum2     = ropp_MDFV
        ro_data%lev2c%pblh_rhum_flag = ropp_MIFV

        CALL ropp_apps_pblh_rhum(ro_data, diag=l_diag)

        CALL pblh_flag_decode(ro_data%lev2c%pblh_rhum_flag, flag_meaning)

        CALL message(msg_diag, "Relative-humidity-based PBLH QC flag = " // &
                                TRIM(ADJUSTL(flag_meaning)))

        IF (ro_data%lev2c%pblh_rhum > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_rhum*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pblr_rhum
          CALL message(msg_info, "Relative-humidity-based PBLH: " // &
                                  pblx_str // " % at " // pblh_str // " km \n")

        END IF

        IF (ro_data%lev2c%pblh_rhum2 > ropp_MDTV) THEN

          WRITE (pblh_str, FMT='(F10.5)') ro_data%lev2c%pblh_rhum2*1.0e-3_wp
          WRITE (pblx_str, FMT='(F10.5)') ro_data%lev2c%pblr_rhum2
          CALL message(msg_info, "2nd Relative-humidity-based PBLH: " // &
                                  pblx_str // " % at " // pblh_str // " km \n")

        END IF

      END IF

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

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

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

      CALL ropp_io_free(ro_data)

    END DO  ! loop over profiles

  END DO  ! loop over files

  CALL EXIT(msg_exit_status)

CONTAINS

!-------------------------------------------------------------------------------
! 16. PBLH QC flag decoding
!-------------------------------------------------------------------------------

  SUBROUTINE pblh_flag_decode(pblh_flag, meaning)

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

    WRITE (flag_str, FMT='(I5)') pblh_flag

    meaning = flag_str // " ==> "

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

    IF (MODULO(pblh_flag, 2**PBLH_QC_missing_lon) == 0) THEN
      meaning = TRIM(meaning) // " Diagnosed value OK ..."
      IF (pblh_flag == 0) RETURN  ! Won't ever be true ... yet
    END IF

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

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

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

    IF (BTEST(pblh_flag, PBLH_QC_too_low)) &
      meaning = TRIM(meaning) // " Diagnosed PBLH too low ..."

    IF (BTEST(pblh_flag, PBLH_QC_too_high)) &
      meaning = TRIM(meaning) // " Diagnosed PBLH too high ..."

    IF (BTEST(pblh_flag, PBLH_QC_missing_lon)) &
      meaning = TRIM(meaning) // " Missing longitude ..."

    IF (BTEST(pblh_flag, PBLH_QC_missing_lat)) &
      meaning = TRIM(meaning) // " Missing latitude ..."

    IF (BTEST(pblh_flag, PBLH_QC_double_pblh)) &
      meaning = TRIM(meaning) // " Double PBLH ..."

    IF (BTEST(pblh_flag, PBLH_QC_multiple_pblh)) &
      meaning = TRIM(meaning) // " Multiple PBLH ..."

    IF (BTEST(pblh_flag, PBLH_QC_land)) &
      meaning = TRIM(meaning) // " Land point ..."

    IF (BTEST(pblh_flag, PBLH_QC_coast)) &
      meaning = TRIM(meaning) // " Coast point ..."

    IF (BTEST(pblh_flag, PBLH_QC_polar)) &
      meaning = TRIM(meaning) // " Polar point ..."

    IF (BTEST(pblh_flag, PBLH_QC_subtropics)) &
      meaning = TRIM(meaning) // " Subtropical point ..."

    IF (BTEST(pblh_flag, PBLH_QC_tropics)) &
      meaning = TRIM(meaning) // " Tropical point ..."

  END SUBROUTINE pblh_flag_decode

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

  SUBROUTINE usage()
    PRINT *, 'Purpose:'
    PRINT *, '  Planetary boundary layer height diagnostic'
    PRINT *, 'Usage:'
    PRINT *, '  > ropp_apps_pblh_tool [<options>] <input_file(s)>'
    PRINT *, 'Options:'
    PRINT *, '  -o <output_file>  name of ROPP netCDF output file'
    PRINT *, '  -b                calculate bending-angle-based PBLH'
    PRINT *, '  -n                calculate refractivity-based PBLH'
    PRINT *, '  -y                calculate dry-temperature-based PBLH'
    PRINT *, '  -t                calculate temperature-based PBLH'
    PRINT *, '  -q                calculate specific-humidity-based PBLH'
    PRINT *, '  -r                calculate relative-humidity-based PBLH'
    PRINT *, '  -h                this help'
    PRINT *, '  -d                output additional diagnostics'
    PRINT *, '  -v                version information'
    PRINT *, ''
  END SUBROUTINE usage

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

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

END PROGRAM ropp_apps_pblh_tool