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

SUBROUTINE ropp_apps_pblh_tdry(ro_data, diag)

!****s* PlanetaryBoundaryLayerHeight/ropp_apps_pblh_tdry *
!
! NAME
!   ropp_apps_pblh_tdry
!
! SYNOPSIS
!   Planetary boundary layer height diagnostic based on dry temperature
!
!   CALL ropp_apps_pblh_tdry(ro_data, diag)
!
! DESCRIPTION
!   Diagnose planetary boundary layer height from the kinks in the dry temperature profile,
!   using the gradient, as recommended by Xie, 2014 (SAF/ROM/DMI/REP/VS21/001).
!
! INPUTS
!   TYPE(ROprof), INTENT(INOUT)    :: ro_data    ! input RO profile containing lev2a data
!   LOGICAL, OPTIONAL, INTENT(IN)  :: diag       ! extra diagnostics required
!
! OUTPUTS
!   TYPE(ROprof), INTENT(INOUT)    :: ro_data    ! output RO profile containing lev2a data
!
! 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_constants, ONLY: dtor
  USE ropp_apps_utils, ONLY: ropp_apps_calc_tdry, ropp_apps_pblh_locate, &
                             ropp_apps_pblh_check_min_height, &
                             ropp_apps_pblh_check_max_height
  USE ropp_apps_types
  USE ropp_apps, ONLY: ropp_apps_pblh_region

  IMPLICIT NONE

  TYPE(ROprof), INTENT(INOUT)          :: ro_data    ! input RO profile containing lev2a data
  LOGICAL, OPTIONAL, INTENT(IN)        :: diag       ! extra diagnostics required

! Parameters in the dry temperature-based PBLH diagnosis algorithms

  REAL(wp), PARAMETER                  :: pblh_min=300.0_wp ! m

  REAL(wp), PARAMETER                  :: pblh_max=5000.0_wp ! m

  INTEGER,  PARAMETER                  :: min_no_points=4

! Local variables

  INTEGER, PARAMETER                   :: max_pblh_indexes=100
  INTEGER                              :: pblh_indexes(max_pblh_indexes)
  INTEGER                              :: pblh_index1,pblh_index2
  CHARACTER(LEN=3)                     :: smax_pblh_indexes

  REAL(wp)                             :: pblx1,pblh1
  REAL(wp)                             :: pblx2,pblh2
  REAL(wp)                             :: lat, lon, z_star
  REAL(wp)                             :: var_1,var_n

  REAL(wp), DIMENSION(:), ALLOCATABLE  :: var_x,var_y,var_grad,var_y1
  REAL(wp), DIMENSION(:), ALLOCATABLE  :: geom
  REAL(wp), DIMENSION(:), ALLOCATABLE  :: tdry

  LOGICAL                              :: pblh_possible
  LOGICAL                              :: l_diag

  INTEGER                              :: pblh_qc_flag,n_pblh
  INTEGER                              :: n_points,n_points1
  INTEGER                              :: i,i1
  INTEGER                              :: i_lon, i_lat

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

  LOGICAL, DIMENSION(:), ALLOCATABLE   :: lregion

! PBLH region variables
  INTEGER, PARAMETER                       :: PBLH_nlon=90
  REAL(wp), DIMENSION(PBLH_nlon)           :: PBLH_lon
  INTEGER, PARAMETER                       :: PBLH_nlat=46
  REAL(wp), DIMENSION(PBLH_nlat)           :: PBLH_lat
  INTEGER, DIMENSION(PBLH_nlon, PBLH_nlat) :: PBLH_region


!-------------------------------------------------------------------------------
! 2. Calculate PBLH based on dry temperature gradient
!-------------------------------------------------------------------------------

! 2.0 Initial messages
! --------------------

  CALL message_get_routine(routine)

  CALL message_set_routine('ropp_apps_pblh_tdry')

  CALL message(msg_info, "Calculating dry-temperature-based planetary boundary layer height \n")

! 2.1 Initialise variables
! ------------------------

  n_points = ro_data%lev2a%npoints

  pblh_qc_flag = 0

  l_diag = .FALSE.
  IF ( PRESENT(diag) ) l_diag = diag

  IF ((n_points == 0) .OR. ro_data%lev2a%missing) THEN

    CALL message(msg_info, "No (valid) Level 2a data in profile ... " //    &
                           "cannot calculate dry-temperature-based PBLH \n")

    ro_data%lev2c%pblh_tdry_flag = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)

    IF ( l_diag .AND. (SIZE(ro_data%lev2a%alt_refrac) > 0) ) &
      CALL pblh_tdry_dummy_diag(ro_data, SIZE(ro_data%lev2a%alt_refrac))

    CALL message_set_routine(routine)

    RETURN

  END IF

  pblh1 = ropp_MDFV  ;  pblx1 = ropp_MDFV
  pblh2 = ropp_MDFV  ;  pblx2 = ropp_MDFV

  pblh_possible = .TRUE.

  ALLOCATE ( tdry(n_points) )  ;  tdry = ro_data%lev2a%dry_temp

! 2.2 Check for numerically invalid data
! --------------------------------------

  IF (ALL((                    tdry < ropp_MDTV) .OR. &
          (ro_data%lev2a%alt_refrac < ropp_MDTV))) THEN

    CALL message(msg_diag, "No common non-missing dry temperatures and " // &
                           "refractivity altitudes ... will try to calculate Tdry from refractivity \n")

    CALL ropp_apps_calc_tdry(ro_data, tdry)

    IF (ALL((                    tdry < ropp_MDTV) .OR. &
            (ro_data%lev2a%alt_refrac < ropp_MDTV))) THEN

      CALL message(msg_info, "No common non-missing dry temperatures and " // &
                             "refractivity altitudes ... cannot calculate PBLH \n")
      pblh_possible = .FALSE.
      pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)

    END IF

  END IF

  IF (ro_data%GEOref%lat < ropp_MDTV) THEN
    CALL message(msg_warn, "Missing tangent point latitude ... " // &
                           "setting to zero to calculate latitude-dependent parameters \n")
    lat = 0.0_wp
    pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_missing_lat)
  ELSE
    lat = ro_data%GEOref%lat
  END IF

  IF (ro_data%GEOref%lon < ropp_MDTV) THEN
    CALL message(msg_warn, "Missing tangent point longitude ... " // &
                           "cannot calculate longitude-dependent parameters \n")
    pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_missing_lon)
  ELSE
    lon = ro_data%GEOref%lon
  END IF

  IF (ro_data%lev2c%geop_sfc < ropp_MDTV) THEN
    CALL message(msg_warn, "Missing surface geopotential ... assuming it to be zero \n")
    z_star = 0.0_wp
  ELSE
    z_star = ro_data%lev2c%geop_sfc
  END IF

  ALLOCATE(lregion(n_points))  ;  lregion = .FALSE.

  WHERE ((                    tdry > ropp_MDTV) .AND. &
         (ro_data%lev2a%alt_refrac > ropp_MDTV)) lregion = .TRUE.
  n_points1 = COUNT(lregion)

  IF (n_points1 < min_no_points) THEN
    CALL message(msg_info, "Not enough non-missing points for algorithm " // &
                           "to proceed ... cannot calculate PBLH\n")
    pblh_possible = .FALSE.
    pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
  END IF

  IF (BTEST(pblh_qc_flag, PBLH_QC_data_invalid)) THEN
    pblh_possible = .FALSE.
    CALL message(msg_info, "Numerically invalid data")
  END IF

! 2.3 Calculate geometric height above surface corresponding to the given impact parameter
! ----------------------------------------------------------------------------------------

  ALLOCATE (geom(n_points))

  geom = ro_data%lev2a%alt_refrac - geopotential2geometric(lat, z_star)

! 2.4 Check for scientifically invalid data
! -----------------------------------------

  IF (MINVAL(geom, mask=lregion) > pblh_min) THEN
    CALL message(msg_info, "Refractivity altitudes do not start deep enough")
    pblh_possible = .FALSE.
    pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_prof_depth)
  END IF

  IF (MAXVAL(geom, mask=lregion) < pblh_max) THEN
    CALL message(msg_info, "Refractivity altitudes do not reach high enough")
    pblh_possible = .FALSE.
    pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_prof_height)
  END IF

  IF (BTEST(pblh_qc_flag,  PBLH_QC_prof_depth) .OR. &
      BTEST(pblh_qc_flag, PBLH_QC_prof_height)) THEN
    pblh_possible = .FALSE.
    CALL message(msg_info, "Scientifically invalid data")
    pblh_qc_flag  = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
  END IF

! 2.5 Prune out missing data
! --------------------------

  IF (pblh_possible) THEN

    ALLOCATE(var_y(n_points1), var_x(n_points1), var_grad(n_points1)) ! Fill these with the non-missing data

    i1 = 1
    DO i=1,n_points
      IF (lregion(i)) THEN
        var_y(i1) = geom(i)
        var_x(i1) = tdry(i)
        i1        = i1 + 1
      END IF
    END DO

! Below this point, the length of all arrays should be n_points1

! 2.6 Smooth Tdry
! ---------------

    var_1 = var_x(1)   ;  var_n = var_x(n_points1)
    var_x = (CSHIFT(var_x, -1) + 2.0_wp*var_x + CSHIFT(var_x, 1)) / 4.0_wp
    var_x(1) = var_1   ;  var_x(n_points1) = var_n

! 2.7 Compute gradient at half points (i+1/2), assuming T varies linearly with z
! ------------------------------------------------------------------------------

!  NB: var_grad(i) = dTdry/dz(i+1/2)

    var_grad = (CSHIFT(var_x, 1) - var_x) / &
               (CSHIFT(var_y, 1) - var_y + TINY(1.0_wp))

    var_grad(n_points1) = var_grad(n_points1-1)

! 2.8 Append gradient and height to the ROprof structure and thence the output file
! ---------------------------------------------------------------------------------

    IF (l_diag) THEN

      CALL ropp_io_addvar_rodataD1d( ro_data, &
           name      = "pblh_diag_tdry", &
           long_name = "PBLH diagnostic: tdry", &
           units     = "K", &
           range     = (/ -1000000.0_wp, 1000000.0_wp /), &
           DATA      = UNPACK(var_x, MASK=lregion, FIELD=(/ (ropp_MDFV, i=1,n_points) /)) )

      CALL ropp_io_addvar_rodataD1d( ro_data, &
           name      = "pblh_diag_tdry_gradient", &
           long_name = "PBLH diagnostic: tdry gradient", &
           units     = "K/m", &
           range     = (/ -1000000.0_wp, 1000000.0_wp /), &
           DATA      = UNPACK(var_grad, MASK=lregion, FIELD=(/ (ropp_MDFV, i=1,n_points) /)) )

      ALLOCATE(var_y1(n_points1))  ! This is where gradient is defined (approximately)
      var_y1 = (var_y + CSHIFT(var_y, 1)) / 2.0_wp
      var_y1(n_points1) = var_y(n_points1)
      CALL ropp_io_addvar_rodataD1d( ro_data, &
           name      = "pblh_diag_tdry_height", &
           long_name = "PBLH diagnostic: tdry height", &
           units     = "m", &
           range     = (/ -1000000.0_wp, 1000000.0_wp /), &
           DATA      = UNPACK(var_y1, MASK=lregion, FIELD=(/ (ropp_MDFV, i=1,n_points) /)) )

      DEALLOCATE(var_y1)

    END IF

    DEALLOCATE(lregion)

! 2.9 Find indices and number of local maximum vertical gradients
! ----------------------------------------------------------------

    n_pblh = 0

    pblh_indexes(:) = -1

    DO i=n_points1-1,2,-1

      IF ( var_y(i) > pblh_max ) CYCLE

      IF ( var_y(i) < pblh_min ) EXIT

      IF ( (var_grad(i+1) < var_grad(i)) .AND. &
           (var_grad(i-1) < var_grad(i)) ) THEN

        n_pblh = n_pblh + 1

        IF ( n_pblh > max_pblh_indexes ) THEN
          WRITE ( smax_pblh_indexes, FMT='(I3)' ) max_pblh_indexes
          CALL message( msg_info, 'More than ' // smax_pblh_indexes // &
                                  ' possible PBLHs found ... ' // &
                                  ' will examine the highest ones.' )
          EXIT
        END IF

        pblh_indexes(n_pblh) = i

      END IF

    END DO

! 2.10 Estimate location of PBLH(s) by interpolating dlapse/dz
! ------------------------------------------------------------

    SELECT CASE (n_pblh)

      CASE (0)

        CALL message(msg_info, "Could not detect PBLH ... " // &
                               "leaving as missing data \n")

        pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)

      CASE (1)

        CALL message(msg_info, "  1 possible PBLH found \n")

        pblh_index1 = pblh_indexes(n_pblh)

        CALL ropp_apps_pblh_locate(pblh_index1, var_x, var_y, var_grad, pblh1, pblx1, pblh_qc_flag)

      CASE (2:)

        WRITE ( str_n_pblh, FMT='(I3)' ) n_pblh
        CALL message(msg_diag, str_n_pblh // &
                               " possible PBLHs found \n")

        IF ( n_pblh == 2 ) THEN
          pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_double_pblh)
        ELSE
          pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_multiple_pblh)
        END IF

        pblh_index1 = pblh_indexes(SUM(MAXLOC(var_grad(pblh_indexes(1:n_pblh)))))
        pblh_index2 = pblh_indexes(SUM(MAXLOC(var_grad(pblh_indexes(1:n_pblh)), &
                                              MASK=(pblh_indexes(1:n_pblh) /= pblh_index1))))

        CALL ropp_apps_pblh_locate(pblh_index1, var_x, var_y, var_grad, pblh1, pblx1, pblh_qc_flag)
        CALL ropp_apps_pblh_locate(pblh_index2, var_x, var_y, var_grad, pblh2, pblx2, pblh_qc_flag)

    END SELECT

! 2.11 Check that PBLH is not too low
! -----------------------------------

    IF (.NOT. BTEST(pblh_qc_flag, PBLH_QC_data_invalid)) &
      CALL ropp_apps_pblh_check_min_height(n_pblh, pblh_min, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)

! 2.12 Check that PBLH is not too high
! ------------------------------------

    IF (.NOT. BTEST(pblh_qc_flag, PBLH_QC_data_invalid)) &
      CALL ropp_apps_pblh_check_max_height(n_pblh, pblh_max, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)

! 2.13 Calculate region for QC purposes
! -------------------------------------

    IF ( BTEST(pblh_qc_flag, PBLH_QC_missing_lon) .OR. &
         BTEST(pblh_qc_flag, PBLH_QC_missing_lat) ) THEN

      CALL message(msg_diag, "Cannot calculate PBLH geographical region")

    ELSE

      CALL message(msg_diag, "Calculating PBLH geographical region")

      CALL ropp_apps_pblh_region(pblh_lon, pblh_lat, pblh_region)

! Find index of nearest longitude
      i_lon = SUM ( MINLOC ( ABS ( ATAN2 ( SIN(dtor*(lon-pblh_lon)), &
                                           COS(dtor*(lon-pblh_lon)) ) ) ) )

! Find index of nearest latitude
      i_lat = SUM ( MINLOC ( ABS ( ATAN2 ( SIN(dtor*(lat-pblh_lat)), &
                                           COS(dtor*(lat-pblh_lat)) ) ) ) )

! Update QC flag accordingly
      pblh_qc_flag = IBSET(pblh_qc_flag, pblh_region(i_lon, i_lat) + PBLH_QC_land) ! Land points are zero in the region mask

    END IF

! 2.14 Clean up
! -------------

    DEALLOCATE(var_y, var_x, var_grad)

  ELSE ! pblh_possible = .FALSE.

    IF ( l_diag .AND. (SIZE(ro_data%lev2a%alt_refrac) > 0) ) &
      CALL pblh_tdry_dummy_diag(ro_data, SIZE(ro_data%lev2a%alt_refrac))

  END IF ! pblh_possible = .TRUE.

  DEALLOCATE (tdry, geom)


!-------------------------------------------------------------------------------
! 3. Copy PBLH variables to ROprof structure
!-------------------------------------------------------------------------------

  ro_data%lev2c%pblh_tdry      = pblh1     ! Return the PBLH refractivity altitude above surface
  ro_data%lev2c%pblt_tdry      = pblx1     ! Return the PBLH dry temperature
  ro_data%lev2c%pblh_tdry2     = pblh2     ! Return the 2nd PBLH refractivity altitude above surface
  ro_data%lev2c%pblt_tdry2     = pblx2     ! Return the 2nd PBLH dry temperature
  ro_data%lev2c%pblh_tdry_flag = pblh_qc_flag


!-------------------------------------------------------------------------------
! 4. Reset routine name
!-------------------------------------------------------------------------------

  CALL message_set_routine(routine)


CONTAINS


!-------------------------------------------------------------------------------
! 5. Dummy diagnostics
!-------------------------------------------------------------------------------

  SUBROUTINE pblh_tdry_dummy_diag (rodata, n)

! Add 1D-arrays, of length n, of missing data to the ROprof structure rodata.

    USE ropp_io_types, ONLY: ROprof

    IMPLICIT NONE

    TYPE(ROprof), INTENT(INOUT)   :: rodata
    INTEGER                       :: i, n

    CALL ropp_io_addvar_rodataD1d( rodata, &
         name      = "pblh_diag_tdry", &
         long_name = "PBLH diagnostic: tdry", &
         units     = "K", &
         range     = (/ -1000000.0_wp, 1000000.0_wp /), &
         DATA      = (/ (ropp_MDFV, i=1,n) /))

    CALL ropp_io_addvar_rodataD1d( rodata, &
         name      = "pblh_diag_tdry_gradient", &
         long_name = "PBLH diagnostic: tdry gradient", &
         units     = "K/m", &
         range     = (/ -1000000.0_wp, 1000000.0_wp /), &
         DATA      = (/ (ropp_MDFV, i=1,n) /))

    CALL ropp_io_addvar_rodataD1d( rodata, &
         name      = "pblh_diag_tdry_height", &
         long_name = "PBLH diagnostic: tdry height", &
         units     = "m", &
         range     = (/ -1000000.0_wp, 1000000.0_wp /), &
         DATA      = (/ (ropp_MDFV, i=1,n) /))

  END SUBROUTINE pblh_tdry_dummy_diag


END SUBROUTINE ropp_apps_pblh_tdry