MODULE AMV_EN_TARGET_SELECTION_UTILS_M

USE fw_log_mod

USE NF_PARM
USE GEOCAT_CONSTANTS

USE TYPE_KINDS_AIT
USE COMMON_VAR_VALUE
USE Framework_Global_Variables_Module
USE NUMERICAL_ROUTINES

USE AMV_EN_Forecast_M
USE AMV_EN_FEATURE_TRACKING_UTILS_M
USE AMV_EN_WINDS_INC


IMPLICIT NONE

! Variables accessible to all subroutines in this module
INTEGER(LONG) :: Good_Winds
INTEGER(LONG) :: Box_Size
INTEGER(LONG) :: Total_Points

INTEGER(LONG) :: Channel


! Satellite ID
INTEGER(LONG) :: SENSOR_SERIES

! Clear-sky water vapor flag
INTEGER(LONG) :: CSWV_Flag



REAL(SINGLE)   :: MIN_PERCENT_PRESS_VALUES
REAL(SINGLE)   :: IR_GRADIENT_THRESHOLD





CONTAINS



!-------------------------------------------------------------------------------
!
! Name:
!   Convolve_X
!
! Function:
!   Compute gradient in element direction using a 4-pt centered difference
!
! Description:
!   - Use 5-pt kernel to estimate derivative
!
! Reference:
!   Barron et al. (1994), "Performance of Optical Flow Techniques."
!   International Journal of Computer Vision, 12:1, 43-77.
!
! Inputs:
!   Kernel               5-pt kernel
!   Element              Image element location
!   Line                 Image line location
!   PIXEL_OFFSET         Pixel offset used to compute centered difference
!
! Output:
!   Diff_X               Derivative in element direction
!
! Dependencies:
!   None
!
! Restrictions:
!   - May convert this function to a generic interface
!
! History:
!   4/2007 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
FUNCTION Convolve_X(Kernel,Element,Line,BrtTemp,PIXEL_OFFSET) RESULT(Diff_X)

  INTEGER(LONG), INTENT(IN) :: Line
  INTEGER(LONG), INTENT(IN) :: Element
  INTEGER(LONG), INTENT(IN) :: PIXEL_OFFSET
  INTEGER(LONG) :: Pixel

  REAL(SINGLE), DIMENSION(:), INTENT(IN) :: Kernel
  REAL(SINGLE), DIMENSION(:,:), INTENT(IN) :: BrtTemp
  REAL(SINGLE) :: Diff_X

  Diff_X = 0.0

  sum_loop: DO Pixel=-PIXEL_OFFSET, PIXEL_OFFSET
    !what about missing values?
    Diff_X = Diff_X + Kernel(Pixel+3)*BrtTemp(Element+Pixel,Line)

  END DO sum_loop

END FUNCTION Convolve_X

!-------------------------------------------------------------------------------
!
! Name:
!   Convolve_Y
!
! Function:
!   Compute gradient in line direction using a 4-pt centered difference
!
! Description:
!   - Use 5-pt kernel to estimate derivative
!
! Reference:
!   Barron et al. (1994), "Performance of Optical Flow Techniques."
!   International Journal of Computer Vision, 12:1, 43-77.
!
! Inputs:
!   Kernel               5-pt kernel
!   Element              Image element location
!   Line                 Image line location
!   PIXEL_OFFSET         Pixel offset used to compute centered difference
!
! Output:
!   Diff_Y               Derivative in line direction
!
! Dependencies:
!   None
!
! Restrictions:
!   - May convert this function to a generic interface
!
! History:
!   4/2007 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
FUNCTION Convolve_Y(Kernel,Element,Line,BrtTemp,PIXEL_OFFSET) RESULT(Diff_Y)

  INTEGER(LONG), INTENT(IN) :: Line
  INTEGER(LONG), INTENT(IN) :: Element
  INTEGER(LONG), INTENT(IN) :: PIXEL_OFFSET
  INTEGER(LONG) :: Pixel

  REAL(SINGLE), DIMENSION(:), INTENT(IN) :: Kernel
  REAL(SINGLE), DIMENSION(:,:), INTENT(IN) :: BrtTemp
  REAL(SINGLE) :: Diff_Y

  Diff_Y = 0.0

  sum_loop: DO Pixel=-PIXEL_OFFSET, PIXEL_OFFSET
    !what about missing values?
    Diff_Y = Diff_Y + Kernel(Pixel+3)*BrtTemp(Element,Line+Pixel)

  END DO sum_loop

END FUNCTION Convolve_Y

! ------------------------------------------------------------------------------
!
! Name:
!   Target_QC
!
! Function:
!   Perform quality control of target scene.
!
! Description:
!   Apply various tests to determine if scene is a suitable target. Tests
!   include:
!
!               contrast (max BrtTemp - min BrtTemp)
!               earth edge (space)
!               cloud amount
!                 - cloud tracers; at least 10 percent cloud coverage
!                 - clear sky tracer; 100 percent clear coverage
!               bad data (reasonable data values)
!               spatial coherence
!               multi-layer check
!               day/night terminator check (visible and 3.9 channel only)
!
! Reference:
!   Nieman et al (1997), "Fully Automated Cloud-Drift Winds in NESDIS
!   Operations." Bulletin of the American Meteorological Society, June 1997.
!
!   Nieman et al (1993), "A Comparison of Several Techniques to Assign Heights
!   to Cloud Tracers." Journal of Applied Meteorology, September 1993.
!
!   Coakley and Bretherton (1982), "Cloud Cover From High-Resolution Scanner
!   Data: Detecting and Allowing for Partially Filled Fields of View." Journal
!   of Geophysical Research, June 1982.
!
! Calling sequence:
!   CALL Target_QC
!
! Inputs:
!   Box_Data                   Data structure holding target box data arrays
!
! Output:
!   QC_Flag              Quality control flag, possible values include:
!                          0 - Good target scene
!                          1 - Target failed contrast test
!                          2 - Target located on earth edge
!
!                          3 - Cloud amount below 10 percent
!                                         or
!                              Target scene not clear enough
!
!                          4 - Bad data (max BrtTemp <= 180 or min BrtTemp = 0)
!                          5 - Bad data (max BrtTemp > 340)
!                          6 - multi-layer scene
!                          7 - target too coherent
!                         13 - target too close to day/night terminator
!                              (0.6 and 3.9 channels only)
!
! Dependencies:
!   Relies on satellite data structure for brightness temperature data.
!   Requires spatial coherence arrays output from "compute_spatial_uniformity"
!   subroutine.
!
! Restrictions:
!   None
!
! History:
!   5/2007 - Wayne Bresky - Created
!  10/2008 - Wayne Bresky - Added check for day/night terminator
!  04/2010 - Wayne Bresky - Introduced data structure to hold target box data
!-------------------------------------------------------------------------------
SUBROUTINE Target_QC(NavType, Box_Data, QC_Flag)

  CHARACTER(LEN=*), INTENT(IN) :: NavType
  TYPE(Target_Box_Data), INTENT(INOUT) :: Box_Data
  INTEGER(LONG), INTENT(OUT) :: QC_Flag

  INTEGER(LONG) :: Number_Of_Lines
  INTEGER(LONG) :: Number_Of_Elements
  INTEGER(LONG) :: Total_Sample
  INTEGER(LONG) :: Line
  INTEGER(LONG) :: Element
  INTEGER(LONG) :: Peak_Sample_Size
  INTEGER(LONG) :: Cold_Sample_Size
  INTEGER(LONG) :: Filtered_Sample_Size
  INTEGER(LONG) :: Cloudy_Pixels
  INTEGER(LONG) :: Cloudy_Pixels_WV
  INTEGER(LONG) :: Cloudy_Threshold
  INTEGER(LONG) :: Sample
  INTEGER(LONG) :: Sum_Of_Samples
  INTEGER(LONG) :: Peak
  INTEGER(LONG) :: Peak_Cold
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: Peak_Slots
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: Cold_Slots

  REAL(SINGLE), DIMENSION(Total_Points) :: Mean_Values_Filtered_Sample
  REAL(SINGLE) :: BrtTemp_Gradient_Threshold
  REAL(SINGLE) :: Sample_Threshold
  REAL(SINGLE) :: Coherence_Threshold
  REAL(SINGLE) :: Min_Box_Value
  REAL(SINGLE) :: Max_Box_Value
  REAL(SINGLE) :: Range_Of_Values
  REAL(SINGLE) :: Solar_Zen_Center


  ! ----------------------------------------------------------------------------
  ! Initialize variables
  ! ----------------------------------------------------------------------------

  Number_Of_Elements = SIZE(Box_Data%Mean,1)
  Number_Of_Lines = SIZE(Box_Data%Mean,2)
  Peak_Sample_Size = 0
  Filtered_Sample_Size = 0
  Cold_Sample_Size = 0
  Total_Sample = 0
  QC_Flag = sym%SUCCESS

  ! The BT extremes are needed for the contrast test

  Max_Box_Value = MAXVAL(Box_Data%BrtTemp)
  Min_Box_Value = MINVAL(Box_Data%BrtTemp)

  Sample_Threshold = 0.0
  Cloudy_Pixels = 0
  Cloudy_Pixels_WV = 0
  Sum_Of_Samples = 0

  ! pick up solar zenith angle at center of target box

  Solar_Zen_Center = Box_Data%Sol_Zen(Number_Of_Elements/2 + 1, &
                                         Number_Of_Lines/2 + 1)

  ! Avoid the day/night terminator when processing visible and 3.9 channels

  IF (Channel .EQ. NATIVE_ABI_CHN2) THEN

    IF (Solar_Zen_Center .GT. MAX_SOLAR_ZEN_VISIBLE) THEN

      QC_Flag = TERMINATOR_FAILURE
      RETURN

    ENDIF

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN7) THEN

    IF (Solar_Zen_Center .LT. MIN_SOLAR_ZEN_SWIR .OR. &
        Solar_Zen_Center .GT. MAX_SOLAR_ZEN_SWIR) THEN

      QC_Flag = TERMINATOR_FAILURE
      RETURN

    ENDIF

  ENDIF

  ! Use a 10 percent cloud amount threshold for cloud-tracked winds
  ! this is the default

  Cloudy_Threshold = NINT(REAL(Number_Of_Elements * Number_Of_Lines) * 0.1)

  ! For the clear-sky WV winds require a completely clear scene
  ! (i.e., no cloudy pixels allowed)

  IF (CSWV_Flag .EQ. sym%YES) Cloudy_Threshold = 0

  line_loop: DO Line=1, Number_Of_Lines

    element_loop: DO Element=1, Number_Of_Elements

      ! ------------------------------------------------------------------------
      ! Discard target if any portion of the box is located in space
      ! ------------------------------------------------------------------------

      IF (Box_Data%Space(Element,Line) == SYM%SPACE) THEN
        QC_Flag = EARTH_EDGE_FAILURE
        RETURN
      ENDIF

      IF (Box_Data%StdDev(Element,Line) .NE. MISSING_VALUE_REAL4) &
                                  Total_Sample = Total_Sample + 1

      ! ------------------------------------------------------------------------
      ! Keep track of number of cloudy pixels in box. Accept all cloudy and
      ! probably cloudy pixels.
      ! ------------------------------------------------------------------------

      ! for polar winds only count cloudy pixels not probably cloudy pixels
      ! this increases the coverage of clear-sky WV winds - 1/25/16

      IF (TRIM(NavType) .EQ. 'PS') THEN

        IF (Box_Data%Cloud_Mask(Element,Line) == SYM%CLOUDY) THEN

          Cloudy_Pixels = Cloudy_Pixels + 1

          ! For the mid- and upper-level WV only count cloudy pixels above 600 mb. This
          ! allows winds to be generated in clear regions above low clouds.

          IF (Channel .EQ. NATIVE_ABI_CHN8 .AND. Box_Data%Cld_Top_Press(Element,Line) .LT. &
              CLR_CLDY_PRESS_CUTOFF) Cloudy_Pixels_WV = Cloudy_Pixels_WV + 1

          IF (Channel .EQ. NATIVE_ABI_CHN9 .AND. Box_Data%Cld_Top_Press(Element,Line) .LT. &
              CLR_CLDY_PRESS_CUTOFF) Cloudy_Pixels_WV = Cloudy_Pixels_WV + 1

        ENDIF

      ELSE

        IF (Box_Data%Cloud_Mask(Element,Line) == SYM%CLOUDY .OR. &
            Box_Data%Cloud_Mask(Element,Line) == SYM%PROB_CLOUDY) THEN

          Cloudy_Pixels = Cloudy_Pixels + 1

          ! For the mid- and upper-level WV only count cloudy pixels above 600 mb. This
          ! allows winds to be generated in clear regions above low clouds.

          IF (Channel .EQ. NATIVE_ABI_CHN8 .AND. Box_Data%Cld_Top_Press(Element,Line) .LT. &
              CLR_CLDY_PRESS_CUTOFF) Cloudy_Pixels_WV = Cloudy_Pixels_WV + 1

          IF (Channel .EQ. NATIVE_ABI_CHN9 .AND. Box_Data%Cld_Top_Press(Element,Line) .LT. &
              CLR_CLDY_PRESS_CUTOFF) Cloudy_Pixels_WV = Cloudy_Pixels_WV + 1

        ENDIF

      ENDIF !1/25/16

      ! Keep track of pixels that are spatially coherent.

      IF (Box_Data%StdDev(Element,Line) .GE. 0.0 .AND. Box_Data%StdDev &
              (Element,Line) .LE. STDDEV_THRESHOLD .AND. Box_Data%Mean &
                                        (Element, Line) .GT. 0.0) THEN

        Filtered_Sample_Size = Filtered_Sample_Size + 1
        Mean_Values_Filtered_Sample(Filtered_Sample_Size) = &
                             Box_Data%Mean(Element, Line)

      ENDIF

    END DO element_loop

  END DO line_loop

  ! ----------------------------------------------------------------------------
  ! Apply cloud amount threshold. For the cloud-top winds, exclude scenes which
  ! do not have at least 10 percent cloud cover. For the clear-sky WV winds,
  ! require that the entire scene be clear. Ignore low clouds when processing
  ! the 6.2 channel.
  ! ----------------------------------------------------------------------------

  ! save cloudy pixel count - wcb 08/31/12

  Box_Data%CloudyPixels = Cloudy_Pixels

  SELECT CASE(Channel)

  CASE (NATIVE_ABI_CHN2, NATIVE_ABI_CHN7, NATIVE_ABI_CHN14)  ! visible, SWIR and LWIR

    IF (Cloudy_Pixels .LT. Cloudy_Threshold) THEN

      QC_Flag = CLOUD_AMOUNT_FAILURE
      RETURN

    ENDIF

  CASE (NATIVE_ABI_CHN8 : NATIVE_ABI_CHN9)  ! mid- and upper-level WV

    ! These bands are unique because we use them to track cloud and clear sky
    ! features. This means we have to test for two conditions.

    IF ((CSWV_Flag .EQ. sym%NO .AND. Cloudy_Pixels_WV .LT. Cloudy_Threshold) .OR. &
        (CSWV_Flag .EQ. sym%YES .AND. Cloudy_Pixels_WV .GT. Cloudy_Threshold)) THEN

      QC_Flag = CLOUD_AMOUNT_FAILURE
      RETURN

    ENDIF

  CASE (NATIVE_ABI_CHN10)  ! low-level WV

    ! For this channel we can't ignore low-level clouds, so check total number
    ! of cloudy pixels and fail if box isn't completely clear. We are only
    ! tracking clear sky features with this band.

    IF (CSWV_Flag .EQ. sym%YES .AND. Cloudy_Pixels .GT. Cloudy_Threshold) THEN

      QC_Flag = CLOUD_AMOUNT_FAILURE
      RETURN

    ENDIF

  END SELECT

  ! ----------------------------------------------------------------------------
  ! Apply contrast check
  ! ----------------------------------------------------------------------------

  Range_Of_Values = Max_Box_Value - Min_Box_Value
  IF (Channel .EQ. NATIVE_ABI_CHN2) THEN

    !aab SWIR - Mods for VIIRS...use different threshold for VIIRS SWIR
    !           VIIRS SWIR is a reflectance value, so it's treated as a VIS channel
    !BrtTemp_Gradient_Threshold = VISIBLE_GRADIENT_THRESHOLD * Box_Size / &
    !                                                   NOMINAL_BOX_SIZE
    IF ( Sensor_Series == SERIES_JPSS_VIIRS ) THEN
      BrtTemp_Gradient_Threshold = VISIBLE_GRADIENT_THRESHOLD_VIIRS * Box_Size / &
                                                          NOMINAL_BOX_SIZE
    ELSE
      BrtTemp_Gradient_Threshold = VISIBLE_GRADIENT_THRESHOLD * Box_Size / &
                                                          NOMINAL_BOX_SIZE
    ENDIF

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN7) THEN

    ! aab Tandem (dbscan) 
    !BrtTemp_Gradient_Threshold = SWIR_GRADIENT_THRESHOLD * Box_Size / &
    !                                            NOMINAL_SWIR_BOX_SIZE
    ! now using lower threshold for SWIR winds 
    BrtTemp_Gradient_Threshold = SWIR_GRADIENT_THRESHOLD

  ! clear-sky WV

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN8  .AND. CSWV_Flag .EQ. sym%YES .OR. &
           Channel .EQ. NATIVE_ABI_CHN9  .AND. CSWV_Flag .EQ. sym%YES .OR. &
           Channel .EQ. NATIVE_ABI_CHN10 .AND. CSWV_Flag .EQ. sym%YES) THEN

    BrtTemp_Gradient_Threshold = WVCS_GRADIENT_THRESHOLD * Box_Size / &
                                                     NOMINAL_BOX_SIZE

  ! cloud-top WV (note: we do not use band 10 for cloud tracking)

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN8 .AND. CSWV_Flag .EQ. sym%NO .OR. &
           Channel .EQ. NATIVE_ABI_CHN9 .AND. CSWV_Flag .EQ. sym%NO) THEN

    BrtTemp_Gradient_Threshold = WVCT_GRADIENT_THRESHOLD * Box_Size / &
                                                     NOMINAL_BOX_SIZE

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN14) THEN

    ! aab Tandem (dbscan) 
    !BrtTemp_Gradient_Threshold = IR_GRADIENT_THRESHOLD * Box_Size / &
    !                                               NOMINAL_BOX_SIZE
    ! now using lower threshold for LWIR winds
    BrtTemp_Gradient_Threshold = IR_GRADIENT_THRESHOLD

  ENDIF


  IF (Range_Of_Values .LT. BrtTemp_Gradient_Threshold) THEN

    QC_Flag = GRADIENT_THRESHOLD_FAILURE
    RETURN

  ENDIF

  ! ----------------------------------------------------------------------------
  ! Flag target scenes with max/min values outside of the acceptable range
  ! ----------------------------------------------------------------------------

  IF (Channel .EQ. NATIVE_ABI_CHN2) THEN

    ! check reflectance values instead of BT but use same error code

    IF (Min_Box_Value .LT. MIN_VALID_REFLECTANCE .OR. &
        Max_Box_Value .GT. MAX_VALID_REFLECTANCE) THEN
      QC_Flag = BAD_BT_FAILURE
      RETURN
    ENDIF

  ELSE

    IF (Min_Box_Value .LT. MIN_VALID_BT .OR. &
        Max_Box_Value .GT. MAX_VALID_BT) THEN

      QC_Flag = BAD_BT_FAILURE
      RETURN

    ENDIF

  ENDIF

  ! ----------------------------------------------------------------------------
  ! For WV channels this is the extent of the target QC.
  ! ----------------------------------------------------------------------------

  ! polar winds test
  ! do not perform additional QC

  !021915
  !IF (Channel .EQ. NATIVE_ABI_CHN8  .OR. Channel .EQ. NATIVE_ABI_CHN9 .OR. &
  !    Channel .EQ. NATIVE_ABI_CHN10) RETURN

  IF (Channel .EQ. NATIVE_ABI_CHN2  .OR. Channel .EQ. NATIVE_ABI_CHN7 .OR. &
      Channel .EQ. NATIVE_ABI_CHN8  .OR. Channel .EQ. NATIVE_ABI_CHN9 .OR. &
      Channel .EQ. NATIVE_ABI_CHN10 ) RETURN


  IF (TRIM(NavType) .EQ. 'PS') RETURN

  ! ----------------------------------------------------------------------------
  ! Apply spatial coherence check
  ! ----------------------------------------------------------------------------

  Coherence_Threshold = COHERENCE_THRESHOLD_PCT * REAL(Total_Sample)

  ! aab Tandem (dbscan) 
  ! disable spatial coherency check
  !IF (REAL(Filtered_Sample_Size) .GT. Coherence_Threshold) THEN

  !  QC_Flag = COHERENCY_FAILURE
  !  RETURN

  !ENDIF

  IF (Filtered_Sample_Size > 0) THEN

    ! --------------------------------------------------------------------------
    ! Perform cluster analysis on filtered sample
    ! --------------------------------------------------------------------------

    CALL Cluster_Analysis(Mean_Values_Filtered_Sample(1:Filtered_Sample_Size), &
                                    Cold_Slots, Peak_Slots, Peak, Peak_Cold)

    ! --------------------------------------------------------------------------
    ! Count the number of warm and cold samples output from cluster analysis
    ! --------------------------------------------------------------------------

    sample_loop: DO Sample= 1, Filtered_Sample_Size

      IF (Cold_Slots(Sample) .EQ. 1) Cold_Sample_Size = Cold_Sample_Size + 1
      IF (Peak_Slots(Sample) .EQ. 1) Peak_Sample_Size = Peak_Sample_Size + 1

    END DO sample_loop

    ! --------------------------------------------------------------------------
    ! Filter out multi-layer cloud scenes (the warmest layer is assumed
    ! to be the surface).
    ! --------------------------------------------------------------------------

    Sum_Of_Samples = Cold_Sample_Size + Peak_Sample_Size
    Sample_Threshold = SAMPLE_THRESHOLD_PCT * Filtered_Sample_Size

    IF (REAL(Sum_Of_Samples) .LT. Sample_Threshold) THEN

      QC_Flag = MULTIPLE_LAYERS_FAILURE
      RETURN

    ENDIF

    DEALLOCATE(Cold_Slots, Peak_Slots)

  ENDIF

END SUBROUTINE Target_QC

! ------------------------------------------------------------------------------
!
! Name:
!   Cluster_Analysis
!
! Function:
!   Identify clusters from spatially coherent (StdDev <= threshold)
!   sample of 3x3 mean values. The sample is filtered before calling this
!   routine.
!
! Description:
!   - A histogram is created from the input 3x3 mean data
!   - The peak frequency is identified
!   - A Gaussian histogram is generated centered on the peak location
!   - Points are selected by applying standard deviation thresholds
!   - The process above is repeated for the peak frequency of the cold end
!
! Reference:
!   Nieman et al (1993), "A Comparison of Several Techniques to Assign Heights
!   to Cloud Tracers." Journal of Applied Meteorology, September 1993.
!   Page 1563.
!
! Calling sequence:
!   CALL Cluster_Analysis
!
! Inputs:
!   Input_Data                Input 3x3 mean data
!   Number_Of_Points          Number of values in Input_Data array
!
! Outputs:
!   Cold_Slots                Selected points from the histogram analysis of the
!                             cold peak
!   Peak_Slots                Selected points from the histogram analysis of the
!                             overall peak
!   Max_Frequency_Slot        Radiance value of histogram peak
!   Max_Frequency_Slot_Cold   Radiance value of cold peak
!
! Dependencies:
!   Requires a spatially coherent filtered sample of 3x3 mean values
!
! Restrictions:
!
! History:
!   5/2007 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE Cluster_Analysis (Input_Data, Cold_Slots, &
                                       Peak_Slots, Max_Frequency_Slot, &
                                              Max_Frequency_Slot_Cold)
  REAL(SINGLE), INTENT(IN), DIMENSION(:) :: Input_Data
  INTEGER(LONG), INTENT(OUT), ALLOCATABLE, DIMENSION(:) :: Cold_Slots
  INTEGER(LONG), INTENT(OUT), ALLOCATABLE, DIMENSION(:) :: Peak_Slots
  INTEGER(LONG), INTENT(OUT) :: Max_Frequency_Slot
  INTEGER(LONG), INTENT(OUT) :: Max_Frequency_Slot_Cold

  INTEGER(LONG), PARAMETER :: Number_Of_StdDev_Cutoffs = 4
  INTEGER(LONG), PARAMETER :: Maximum_Number_Of_Cold_Clusters = 5
  REAL(SINGLE), PARAMETER :: Maximum_Variance = 25.0
  INTEGER(LONG), DIMENSION(Number_Of_StdDev_Cutoffs) :: StdDev_Cutoffs
  INTEGER(LONG):: Number_Of_Points
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: Slot_Values
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: Histogram
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: Gaussian_Histogram
  INTEGER(LONG) :: Slot_Number
  INTEGER(LONG) :: Max_Frequency
  INTEGER(LONG) :: Slot
  INTEGER(LONG) :: Point
  INTEGER(LONG) :: Frequency
  INTEGER(LONG) :: Cold_Counter
  INTEGER(LONG) :: Max_Frequency_Cold

  REAL(SINGLE) :: Variance_LHS
  REAL(SINGLE) :: Variance_RHS
  REAL(SINGLE) :: StdDev_LHS
  REAL(SINGLE) :: StdDev_RHS
  Number_Of_Points = SIZE(Input_Data)
  ALLOCATE(Cold_Slots(Number_Of_Points), Peak_Slots(Number_Of_Points))
  ALLOCATE(Slot_Values(Total_Points), Histogram(NUMBER_OF_SLOTS))

  ! ----------------------------------------------------------------------------
  ! Initialize variables
  ! ----------------------------------------------------------------------------

  StdDev_Cutoffs(1:Number_Of_StdDev_Cutoffs) = 0
  Cold_Slots(1:Number_Of_Points) = 0
  Peak_Slots(1:Number_Of_Points) = 0

  ! ----------------------------------------------------------------------------
  ! Create a 1-D histogram
  ! ----------------------------------------------------------------------------

  CALL Make_Histogram(Input_Data, Slot_Values, &
                   Max_Frequency_Slot, Max_Frequency, Histogram)

  ! ----------------------------------------------------------------------------
  ! Fit half-Gaussian curves to both sides of peak frequency using the three
  ! closest points on either side of the peak to estimate the average variance
  ! on each side.
  ! ----------------------------------------------------------------------------

  CALL Gaussian_Fit (Max_Frequency_Slot, Histogram, Variance_LHS, Variance_RHS)

  ! ----------------------------------------------------------------------------
  ! Constrain variances to an upper limit
  ! ----------------------------------------------------------------------------

  IF (Variance_LHS .GT. Maximum_Variance) Variance_LHS = Maximum_Variance
  IF (Variance_RHS .GT. Maximum_Variance) Variance_RHS = Maximum_Variance

  ! ----------------------------------------------------------------------------
  ! Call routine to compute Gaussian curve for all slots. This is the full
  ! Gaussian curve created from the two average variance values.
  ! ----------------------------------------------------------------------------

  ALLOCATE(Gaussian_Histogram(NUMBER_OF_SLOTS))

  CALL Gaussian_Curve (Max_Frequency_Slot, Max_Frequency, Variance_LHS, &
                                      Variance_RHS, Gaussian_Histogram)

  ! ----------------------------------------------------------------------------
  ! Compute indices corresponding to point selection thresholds. These
  ! thresholds will be used to filter the sample according to the following
  ! rules: Select all points within 1 standard deviation of peak. Between 1 - 3
  ! standard deviations only select as many points as suggested by the Gaussian
  ! fit. Exclude all points more than 3 standard deviations from peak.
  ! ----------------------------------------------------------------------------

  StdDev_LHS = SQRT(Variance_LHS) / INTERVAL_LENGTH
  StdDev_RHS = SQRT(Variance_RHS) / INTERVAL_LENGTH

  StdDev_Cutoffs(1) = Max_Frequency_Slot - SIGMA_3 * StdDev_LHS
  StdDev_Cutoffs(2) = Max_Frequency_Slot - SIGMA_1 * StdDev_LHS
  StdDev_Cutoffs(3) = Max_Frequency_Slot + SIGMA_1 * StdDev_RHS
  StdDev_Cutoffs(4) = Max_Frequency_Slot + SIGMA_3 * StdDev_RHS

  peak_loop: DO Point = 1, Number_Of_Points

    Slot = Slot_Values(Point)

    ! --------------------------------------------------------------------------
    ! Select all points within 1 S.D. of Peak.
    ! --------------------------------------------------------------------------

    IF (Slot .GE. StdDev_Cutoffs(2) .AND. Slot .LE. StdDev_Cutoffs(3)) THEN
      Peak_Slots(Point) = 1

    ! --------------------------------------------------------------------------
    ! On the interval 1 - 3 S.D. select only as many points as suggested by
    ! Gaussian fit.
    ! --------------------------------------------------------------------------

    ELSE IF ((Slot .GE. StdDev_Cutoffs(1) .AND. Slot .LT. StdDev_Cutoffs(2)) &
       .OR. (Slot .GT. StdDev_Cutoffs(3) .AND. Slot .LE. StdDev_Cutoffs(4))) &
                                                                        THEN

      IF (Gaussian_Histogram(Slot) .GT. 0) THEN
        Peak_Slots(Point) = 1
        Gaussian_Histogram(Slot) = Gaussian_Histogram(Slot) - 1
      ENDIF

    ENDIF

  END DO peak_loop

  ! ----------------------------------------------------------------------------
  ! Find peak frequency of 5 coldest clusters.
  ! ----------------------------------------------------------------------------

  Max_Frequency_Cold = 0
  Max_Frequency_Slot_Cold = 0
  Cold_Counter = 0

  ! search from cold end
  slot_loop: DO Slot_Number = 1, NUMBER_OF_SLOTS

    Frequency = Histogram(Slot_Number)

    IF (Frequency .GE. 1) Cold_Counter = Cold_Counter + 1

    IF (Frequency .GT. Max_Frequency_Cold) THEN
      Max_Frequency_Cold = Histogram(Slot_Number)
      Max_Frequency_Slot_Cold = Slot_Number
    ENDIF

    IF (Cold_Counter .GE. Maximum_Number_Of_Cold_Clusters) EXIT slot_loop

  END DO slot_loop

  ! ----------------------------------------------------------------------------
  ! Repeat process used above, this time for the cold sample.
  ! ----------------------------------------------------------------------------

  CALL Gaussian_Fit(Max_Frequency_Slot_Cold, Histogram, Variance_LHS, &
                                                        Variance_RHS)

  IF (Variance_LHS .GT. Maximum_Variance) Variance_LHS = Maximum_Variance
  IF (Variance_RHS .GT. Maximum_Variance) Variance_RHS = Maximum_Variance

  CALL Gaussian_Curve(Max_Frequency_Slot_Cold, Max_Frequency_Cold, &
                   Variance_LHS, Variance_RHS, Gaussian_Histogram)

  StdDev_LHS = SQRT(Variance_LHS) / INTERVAL_LENGTH
  StdDev_RHS = SQRT(Variance_RHS) / INTERVAL_LENGTH

  StdDev_Cutoffs(1) = Max_Frequency_Slot_Cold - SIGMA_3 * StdDev_LHS
  StdDev_Cutoffs(2) = Max_Frequency_Slot_Cold - SIGMA_1 * StdDev_LHS
  StdDev_Cutoffs(3) = Max_Frequency_Slot_Cold + SIGMA_1 * StdDev_RHS
  StdDev_Cutoffs(4) = Max_Frequency_Slot_Cold + SIGMA_3 * StdDev_RHS

  warm_loop: DO Point = 1, Number_Of_Points

    Slot = Slot_Values(Point)

    ! --------------------------------------------------------------------------
    ! Select all points within 1 S.D. of Peak.
    ! --------------------------------------------------------------------------

    IF (Slot .GE. StdDev_Cutoffs(2) .AND. Slot .LE. StdDev_Cutoffs(3)) THEN

      Cold_Slots(Point) = 1

    ! --------------------------------------------------------------------------
    ! On the interval 1 - 3 S.D. select only as many points as suggested by
    ! Gaussian fit.
    ! --------------------------------------------------------------------------

    ELSE IF ((Slot .GE. StdDev_Cutoffs(1) .AND. Slot .LT. StdDev_Cutoffs(2)) &
       .OR. (Slot .GT. StdDev_Cutoffs(3) .AND. Slot .LE. StdDev_Cutoffs(4))) &
                                                                        THEN

      IF (Gaussian_Histogram(Slot) .GT. 0) THEN
        Cold_Slots(Point) = 1
        Gaussian_Histogram(Slot) = Gaussian_Histogram(Slot) - 1
      ENDIF

    ENDIF

  END DO warm_loop

  DEALLOCATE(Histogram, Slot_Values, Gaussian_Histogram)

END SUBROUTINE Cluster_Analysis

! ------------------------------------------------------------------------------
!
! Name:
!   Make_Histogram
!
! Function:
!   Construct a 1-D histogram using image (radiance) data.
!
! Description:
!
! Reference:
!   None
!
! Calling sequence:
!   CALL Make_Histogram
!
! Inputs:
!   Input_Data           Input image data
!   Number_Of_Points     Number of values in Input_Data array
!
! Outputs:
!   Slot_Values          Output array marking each data point's resultant
!                        histogram slot
!   Histogram            Output array of # of points in each histogram
!                        slot (i.e. the histogram)
!
! Dependencies:
!
! Restrictions:
!
! History:
!   5/2007 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE Make_Histogram(Input_Data, Slot_Values, Peak_Slot, Peak_Value, Histogram)

  REAL(SINGLE), INTENT(IN), DIMENSION(:) :: Input_Data
  INTEGER(LONG), DIMENSION(:) :: Slot_Values
  INTEGER(LONG), INTENT(OUT) :: Peak_Slot
  INTEGER(LONG), INTENT(OUT) :: Peak_Value
  INTEGER(LONG), DIMENSION(:), INTENT(OUT)  :: Histogram
  INTEGER(LONG) :: Number_Of_Points
  INTEGER(LONG) :: Point
  INTEGER(LONG) :: Data_Value


  Number_Of_Points = SIZE(Input_Data)
  ! ----------------------------------------------------------------------------
  ! Initialize variables
  ! ----------------------------------------------------------------------------

  Histogram(1:NUMBER_OF_SLOTS) = 0
  Slot_Values(1:Total_Points) = 0
  Peak_Slot = 0
  Peak_Value = 0
  Data_Value = 0

  ! ----------------------------------------------------------------------------
  ! Loop through data and fill histogram
  ! ----------------------------------------------------------------------------

  data_loop: DO Point=1, Number_Of_Points

    ! --------------------------------------------------------------------------
    ! Calculate which histogram slot data point falls into
    ! --------------------------------------------------------------------------

    Slot_Values(Point) = NINT((Input_Data(Point) - MIN_HISTOGRAM_VALUE) &
                                                     / INTERVAL_LENGTH)
    Data_Value = Slot_Values(Point)

    ! --------------------------------------------------------------------------
    ! Add point to correct histogram slot and keep track of max frequency
    ! --------------------------------------------------------------------------
    !021915
    IF ( Data_Value .GT. 0 .AND. Data_Value .LE. NUMBER_OF_SLOTS) THEN
       Histogram(Data_Value) = Histogram(Data_Value) + 1

       IF (Histogram(Data_Value) .GT. Peak_Value) THEN
          Peak_Value = Histogram(Data_Value)
          Peak_Slot =  Data_Value
       ENDIF
    ENDIF

  END DO data_loop

END SUBROUTINE Make_Histogram

! ------------------------------------------------------------------------------
!
! Name:
!   Make_Histogram_Temp
!
! Function:
!   Construct a 1-D histogram using temperature data (either BT or CTT) then
!   determine temperature threshold using coldest 20 (WV) or 25 (IR) percent
!   of histogram.
!
! Description:
!
! Reference:
!   None
!
! Calling sequence:
!   CALL Make_Histogram_Temp
!
! Inputs:
!   Box_Data                      Data structure holding target box data
!   Number_Of_Points              Number of values in Temp_Data array
!
! Outputs:
!   Cold_Sample_Counter           Number of points in cold sample
!   Variance_Press                Variance of Cld_Top_Press values (cold sample)
!   Inversion_Sample_Flag         Flag indicating whether or not a low-level
!                                 inversion is present
!
! Dependencies:
!
! Restrictions:
!
! History:
!  08/2007 - Wayne Bresky - Created
!  03/2009 - Wayne Bresky - Added low-level inversion flag
!  07/2009 - Wayne Bresky - Added cloud height values
!  04/2010 - Wayne Bresky - Introduced target box data structure
!
!-------------------------------------------------------------------------------
SUBROUTINE Make_Histogram_Temp(NavType, Box_Data, Number_Of_Points, Cold_Sample_Counter, &
                                         Variance_Press, Inversion_Sample_Flag)
  CHARACTER(len=*), PARAMETER :: FUNC = "Make_Histogram_Temp"
  CHARACTER(LEN=*), INTENT(IN) :: NavType
  TYPE(Target_Box_Data), INTENT(INOUT) :: Box_Data
  INTEGER(LONG), INTENT(IN) :: Number_Of_Points
  INTEGER(LONG), INTENT(OUT) :: Cold_Sample_Counter
  REAL(SINGLE), INTENT(OUT) :: Variance_Press
  INTEGER(BYTE), INTENT(OUT) :: Inversion_Sample_Flag

  INTEGER(LONG) :: Point
  INTEGER(LONG) :: Data_Value
  INTEGER(LONG) :: BrtTemp
  INTEGER(LONG) :: Cold_Sample
  INTEGER(LONG) :: Number_Of_BrtTemp_Slots = BRTTEMP_MAX - BRTTEMP_MIN + 1
  INTEGER(LONG) :: Histogram_Points
  INTEGER(LONG) :: Point_Cutoff
  INTEGER(LONG) :: Number_Of_Bins
  INTEGER(LONG) :: Cold_Slot_Threshold
  INTEGER(LONG) :: Lower_Bound
  INTEGER(LONG) :: Upper_Bound
  INTEGER(SHORT), ALLOCATABLE, DIMENSION(:) :: Slot_Values
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: Histogram
  INTEGER(LONG) :: Count_Sample1
  INTEGER(LONG) :: Count_Sample2
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: Histogram_Sample1
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: Histogram_Sample2
  INTEGER(LONG), SAVE :: Num_Call = 0
  CHARACTER(len=100) :: message

  ! 1-D data
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Input_Data
  REAL(DOUBLE), ALLOCATABLE, DIMENSION(:) :: Cloud_Top_Pressure
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Cld_Hgt
  INTEGER(BYTE), ALLOCATABLE, DIMENSION(:) :: Cld_Top_Press_Quality
  INTEGER(BYTE), ALLOCATABLE, DIMENSION(:) :: Cloud_Mask
  INTEGER(BYTE), ALLOCATABLE, DIMENSION(:) :: LL_Inver_Flag
  !WCB - add cloud height quality flag
  INTEGER(BYTE), ALLOCATABLE, DIMENSION(:) :: CldHgtQF

  ! cold sample arrays
  REAL(DOUBLE), DIMENSION(:), ALLOCATABLE :: Pressure_Cold_Sample
  REAL(SINGLE), DIMENSION(:), ALLOCATABLE :: BrtTemp_Cold_Sample
  REAL(SINGLE), DIMENSION(:), ALLOCATABLE :: Cld_Hgt_Cold_Sample
  REAL(SINGLE) :: Cold_Threshold
  REAL(DOUBLE) :: Sum_Press
  REAL(DOUBLE) :: Sum_Press2



  ALLOCATE(Histogram(Number_Of_BrtTemp_Slots), source = 0)
  ALLOCATE(Slot_Values(Total_Points), source = MISSING_VALUE_INT2)
  ALLOCATE(Histogram_Sample1(Number_Of_BrtTemp_Slots), source = 0)
  ALLOCATE(Histogram_Sample2(Number_Of_BrtTemp_Slots), source = 0)
  ALLOCATE(Pressure_Cold_Sample(Total_Points), source = MISSING_VALUE_REAL8)
  ALLOCATE(Cld_Hgt_Cold_Sample(Total_Points), source = MISSING_VALUE_REAL4)
  ALLOCATE(BrtTemp_Cold_Sample(Total_Points), source = MISSING_VALUE_REAL4)

  Num_Call = Num_Call + 1

  ! ----------------------------------------------------------------------------
  ! Initialize variables
  ! ----------------------------------------------------------------------------

  IF (Channel .EQ. NATIVE_ABI_CHN8  .AND. CSWV_Flag .EQ. sym%YES .OR. &
      Channel .eq. NATIVE_ABI_CHN9  .AND. CSWV_Flag .EQ. sym%YES .OR. &
      Channel .EQ. NATIVE_ABI_CHN10 .AND. CSWV_Flag .EQ. sym%YES) THEN

    ! Clear-sky WV threshold

    Cold_Threshold = WVCS_COLD_THRESHOLD

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN8 .AND. CSWV_Flag .EQ. sym%NO) THEN

    ! Cloud-top WV threshold

    Cold_Threshold = WVCT_COLD_THRESHOLD

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN9 .AND. CSWV_Flag .EQ. sym%NO) THEN

    ! Cloud-top WV threshold

    Cold_Threshold = WVCT_COLD_THRESHOLD

  ELSE

    Cold_Threshold = IR_COLD_THRESHOLD

  ENDIF

  IF(Num_Call .EQ. 1) THEN
      WRITE(message, *) 'Cold_Threshold = ', Cold_Threshold
      CALL fw_log_debug(FUNC, message)
  ENDIF

  Inversion_Sample_Flag = sym%NO
  Data_Value = 0
  Histogram_Points = 0
  Count_Sample1 = 0
  Count_Sample2 = 0
  Point_Cutoff = 0
  Cold_Sample = 0
  Cold_Sample_Counter = 0
  Number_Of_Bins = 0
  Cold_Slot_Threshold = 0
  Sum_Press = 0.0
  Sum_Press2 = 0.0
  Variance_Press = 0.0

  ! Reshape 2-D arrays to 1-D.
  ! For the CSWV winds construct the histogram from brightness temperatures,
  ! otherwise use cloud-top temperatures.
  ALLOCATE(Input_Data(Total_Points))
  IF (CSWV_Flag .EQ. sym%YES) THEN
    Input_Data = RESHAPE(Box_Data%BrtTemp, (/Total_Points/))
  ELSE
    Input_Data = RESHAPE(Box_Data%Cld_Top_Temp, (/Total_Points/))
  ENDIF

  ALLOCATE(Cloud_Mask(Total_Points))
  Cloud_Mask = RESHAPE(Box_Data%Cloud_Mask, (/Total_Points/))

  ALLOCATE(Cloud_Top_Pressure(Total_Points))
  Cloud_Top_Pressure = RESHAPE(Box_Data%Cld_Top_Press, (/Total_Points/))

  ALLOCATE(Cld_Top_Press_Quality(Total_Points))
  Cld_Top_Press_Quality = RESHAPE(Box_Data%Cld_Top_Press_Quality, &
                                                (/Total_Points/))

  ALLOCATE(LL_Inver_Flag(Total_Points))
  LL_Inver_Flag = RESHAPE(Box_Data%LL_Inver, (/Total_Points/))

  ALLOCATE(Cld_Hgt(Total_Points))
  Cld_Hgt = RESHAPE(Box_Data%Cld_Hgt, (/Total_Points/))

  !WCB - add cloud height quality flag
  ALLOCATE(CldHgtQF(Total_Points))
  CldHgtQF = RESHAPE(Box_Data%CldHgtQF, (/Total_Points/))

  ! Set min/max bounds of histogram

  Lower_Bound = MAX(MINVAL(Input_Data), MIN_VALID_BT) *  SCALE_FACTOR
  Upper_Bound = MIN(MAXVAL(Input_Data), MAX_VALID_BT) *  SCALE_FACTOR

  ! ----------------------------------------------------------------------------
  ! Loop through data and fill histogram.
  ! ----------------------------------------------------------------------------

  data_loop: DO Point=1, Number_Of_Points

    ! filter out cloudy pixels for the clear-sky WV, otherwise
    ! check for clear pixels and missing CTP values

    IF (TRIM(NavType) .EQ. 'PS' .AND. CSWV_Flag .EQ. sym%YES) THEN

       ! for polar clear sky winds only check for cloudy pixels - 1/25/16
       IF (Cloud_Mask(Point) .EQ. SYM%CLOUDY) CYCLE

    ElSE IF (TRIM(NavType) .EQ. 'PS' .AND. CSWV_Flag .EQ. sym%NO) THEN

       ! for polar cloudy winds check for clear and prob clear pixels - 1/25/16
       IF (Cloud_Mask(Point) .EQ. SYM%CLEAR .OR. &
           Cloud_Mask(Point) .EQ. SYM%PROB_CLEAR .OR. &
           !WCB - add cloud height quality flag
           CldHgtQF(Point) .EQ. CTH_DQF_MARGINAL_RETREVIAL .OR.  &
           Cloud_Top_Pressure(Point) .EQ. MISSING_VALUE_REAL4) CYCLE

    ElSE IF (TRIM(NavType) .NE. 'PS' .AND. CSWV_Flag .EQ. sym%YES) THEN

       ! for geostationary clear sky winds check for both prob cloudy
       ! and cloudy pixels - 1/25/16

      IF (Cloud_Mask(Point) .EQ. SYM%CLOUDY .OR. &
          Cloud_Mask(Point) .EQ. SYM%PROB_CLOUDY) CYCLE

    ELSE IF (TRIM(NavType) .NE. 'PS' .AND. CSWV_Flag .EQ. sym%NO) THEN

      ! for geostationary cloudy winds check for both clear and prob clear
      ! pixels and missing CTP - 1/25/16

      IF (Cloud_Mask(Point) .EQ. SYM%CLEAR .OR. &
          Cloud_Mask(Point) .EQ. SYM%PROB_CLEAR  .OR. &
          !WCB - add cloud height quality flag
          CldHgtQF(Point) .EQ. CTH_DQF_MARGINAL_RETREVIAL .OR.  &
          Cloud_Top_Pressure(Point) .EQ. MISSING_VALUE_REAL4) CYCLE

    ENDIF

    ! Check for reasonable temperature values

    IF (Input_Data(Point) .EQ. MISSING_VALUE_REAL4 .OR. &
        Input_Data(Point) .LT. MIN_VALID_BT .OR. &
        Input_Data(Point) .GT. MAX_VALID_BT) CYCLE

    ! --------------------------------------------------------------------------
    ! Calculate which histogram slot data point falls into
    ! --------------------------------------------------------------------------

    Slot_Values(Point) = NINT(Input_Data(Point) * REAL( SCALE_FACTOR)) - &
                                                        BRTTEMP_MIN + 1

    Data_Value = Slot_Values(Point)

    ! --------------------------------------------------------------------------
    ! Check low-level inversion flag before adding point to correct histogram
    ! slot.
    ! --------------------------------------------------------------------------

    IF (LL_Inver_Flag(Point) .EQ. sym%YES) THEN

      Count_Sample2 = Count_Sample2 + 1
      Histogram_Sample2(Data_Value) = Histogram_Sample2(Data_Value) + 1

    ELSE

      Count_Sample1 = Count_Sample1 + 1
      Histogram_Sample1(Data_Value) = Histogram_Sample1(Data_Value) + 1

    ENDIF

  END DO data_loop

  ! Use sample with larger number of points

  IF (Count_Sample1 .GT. Count_Sample2) THEN

    Histogram(:) = Histogram_Sample1(:)
    Histogram_Points = Count_Sample1

  ELSE IF (Count_Sample1 .LE. Count_Sample2) THEN

    Histogram(:) = Histogram_Sample2(:)
    Histogram_Points = Count_Sample2
    Inversion_Sample_Flag = sym%YES

  ENDIF

  ! Set threshold number of points based on number of points in histogram
  ! and percentage.

  Point_Cutoff = NINT(REAL(Histogram_Points) * Cold_Threshold)

  ! Scan histogram from the cold end and determine cutoff slot

  threshold_loop: DO BrtTemp = Lower_Bound, Upper_Bound

    Cold_Sample = Cold_Sample + Histogram(BrtTemp)

    IF (Histogram(BrtTemp) .GT. 0) Number_Of_Bins = Number_Of_Bins + 1

    IF (Cold_Sample .GT. Point_Cutoff .AND. Number_Of_Bins .GT. 1) THEN

      Cold_Sample = Cold_Sample - Histogram(BrtTemp)
      Cold_Slot_Threshold = BrtTemp - 1
      EXIT

    ! Keep at least one histogram bin
    ELSE IF (Cold_Sample .GT. Point_Cutoff .AND. Number_Of_Bins .EQ. 1) THEN

      Cold_Sample = Cold_Sample
      Cold_Slot_Threshold = BrtTemp
      EXIT

    ELSE IF (Cold_Sample .LE. Point_Cutoff .AND. BrtTemp .EQ. Upper_Bound) THEN

      Cold_Sample = Cold_Sample
      Cold_Slot_Threshold = BrtTemp
      EXIT

    ENDIF

  END DO threshold_loop

  ! Identify cold sample. Keep track of temperature, pressure and cloud height.

  IF (Cold_Sample .GT. 0) THEN

    DO Point=1, Number_Of_Points

      IF (Slot_Values(Point) .NE. MISSING_VALUE_INT2 .AND. Slot_Values(Point) &
                     .LE. Cold_Slot_Threshold .AND. LL_Inver_Flag(Point) .EQ. &
                                                  Inversion_Sample_Flag) THEN

        ! Add point to cold sample array

        Cold_Sample_Counter = Cold_Sample_Counter + 1
        Pressure_Cold_Sample(Cold_Sample_Counter) = Cloud_Top_Pressure(Point)
        BrtTemp_Cold_Sample(Cold_Sample_Counter) = Input_Data(Point)
        Sum_Press2 = Sum_Press2 + Cloud_Top_Pressure(Point)**2
        Cld_Hgt_Cold_Sample(Cold_Sample_Counter) = Cld_Hgt(Point)

      ENDIF

    ENDDO

    ! Compute the variance of the cold sample

    IF (Cold_Sample_Counter .EQ. 1) THEN

      Sum_Press = Pressure_Cold_Sample(1)
      Variance_Press = 0.0

    ELSE IF (Cold_Sample_Counter .GT. 1) THEN

      Sum_Press = SUM(Pressure_Cold_Sample(1:Cold_Sample_Counter))
      Variance_Press = Sum_Press2 / Cold_Sample_Counter - &
                     (Sum_Press / Cold_Sample_Counter)**2

    ENDIF

    ! Copy array data to target box data structure

    Box_Data%Pressure_Cold_Sample(:) = Pressure_Cold_Sample(:)
    Box_Data%BrtTemp_Cold_Sample(:) = BrtTemp_Cold_Sample(:)
    Box_Data%Cld_Hgt_Cold_Sample(:) = Cld_Hgt_Cold_Sample(:)

  ENDIF

  DEALLOCATE(Histogram, Histogram_Sample1, Histogram_Sample2, Slot_Values, &
        Input_Data, Cloud_Mask, Cloud_Top_Pressure, Cld_Top_Press_Quality, &
                 LL_Inver_Flag, Pressure_Cold_Sample, BrtTemp_Cold_Sample, &
                                             Cld_Hgt_Cold_Sample, Cld_Hgt, &
                                             !WCB - add cloud height quality
                                             !flag
                                             CldHgtQF)


END SUBROUTINE Make_Histogram_Temp

! ------------------------------------------------------------------------------
!
! Name:
!   Gaussian_Fit
!
! Function:
!   Fit two half-Gaussian curves to a histogram peak
!
! Description:
!
! Reference:
!   None
!
! Calling sequence:
!   CALL Gaussian_Fit
!
! Inputs:
!   Center_Point         Center histogram point (i.e., peak) to fit curve to
!   Histogram            Array of histogram values
!
! Outputs:
!   Avg_Variance_LHS     Variance describing left-hand-side Gaussian curve
!   Avg_Variance_RHS     Variance describing right-hand-side Gaussian curve
!
! Dependencies:
!
! Restrictions:
!
! History:
!   5/2007 - Wayne Bresky - Created
!   9/2016 - Updated for boundary check
!
!-------------------------------------------------------------------------------
SUBROUTINE Gaussian_Fit(Center_Point, Histogram, Avg_Variance_LHS, &
                                                 Avg_Variance_RHS)

  INTEGER(LONG), INTENT(IN) :: Center_Point
  INTEGER(LONG), INTENT(IN), DIMENSION(:) :: Histogram
  INTEGER(LONG)  :: Second_Point
  INTEGER(LONG)  :: Non_Zero_Count_LHS
  INTEGER(LONG)  :: Non_Zero_Count_RHS
  INTEGER(LONG)  :: Lower_Bound
  INTEGER(LONG)  :: Upper_Bound

  REAL(SINGLE), INTENT(OUT) :: Avg_Variance_LHS
  REAL(SINGLE), INTENT(OUT) :: Avg_Variance_RHS
  REAL(SINGLE) :: Value1_Center
  REAL(SINGLE) :: Value2_Center
  REAL(SINGLE) :: Average_Value_Center
  REAL(SINGLE) :: Value1_Second
  REAL(SINGLE) :: Value2_Second
  REAL(SINGLE) :: Average_Value_Second
  REAL(SINGLE) :: Denominator
  REAL(SINGLE) :: Numerator
  REAL(SINGLE) :: Central_Frequency
  REAL(SINGLE) :: Second_Frequency
  REAL(SINGLE) :: Variance
  REAL(SINGLE) :: Sum_Variance_LHS
  REAL(SINGLE) :: Sum_Variance_RHS

  ! ----------------------------------------------------------------------------
  ! Initialize variables
  ! ----------------------------------------------------------------------------

  Non_Zero_Count_LHS = 0
  Non_Zero_Count_RHS = 0
  Value1_Center = 0.0
  Value2_Center = 0.0
  Average_Value_Center = 0.0
  Central_Frequency = 0.0
  Value1_Second = 0.0
  Value2_Second = 0.0
  Average_Value_Second = 0.0
  Second_Frequency = 0.0
  Numerator = 0.0
  Denominator = 0.0
  Avg_Variance_LHS = 0.0
  Avg_Variance_RHS = 0.0
  Sum_Variance_LHS = 0.0
  Sum_Variance_RHS = 0.0
  Lower_Bound = 0
  Upper_Bound = NUMBER_OF_SLOTS

  ! ----------------------------------------------------------------------------
  ! Compute average value of center point.
  ! Note: Averaging may be unneccessary. Might use actual value of peak.
  ! ----------------------------------------------------------------------------

  Value1_Center = MIN_HISTOGRAM_VALUE + REAL(Center_Point - 1) * INTERVAL_LENGTH
  Value2_Center = Value1_Center + INTERVAL_LENGTH
  Average_Value_Center = (Value1_Center + Value2_Center) / 2.0
  Central_Frequency = REAL(Histogram(Center_Point))

  IF (Central_Frequency .LT. 0.0) Central_Frequency = 0.0

  point_loop: DO Second_Point = Center_Point - 3, Center_Point + 3

    ! --------------------------------------------------------------------------
    ! Avoid center point and bounds of histogram
    ! --------------------------------------------------------------------------

    !IF (Second_Point .EQ. Center_Point) CYCLE
    IF (Second_Point .EQ. Center_Point .OR. &
        Second_Point .LE. Lower_Bound .OR. &
        Second_Point .GT. Upper_Bound) CYCLE

    Value1_Second = MIN_HISTOGRAM_VALUE + REAL(Second_Point - 1) &
                                               * INTERVAL_LENGTH
    Value2_Second = Value1_Second + INTERVAL_LENGTH
    Average_Value_Second = (Value1_Second + Value2_Second) / 2.0
    Second_Frequency = REAL(Histogram(Second_Point))

    IF (Second_Frequency .LT. 0.0) Second_Frequency = 0.0

    IF (Second_Frequency .NE. 0.0 .AND. Central_Frequency .NE. 0.0 .AND. &
                           Second_Frequency .NE. Central_Frequency) THEN
      Numerator = -1.0 * (Average_Value_Second - Average_Value_Center)**2
      Denominator = 2.0 * ALOG(Second_Frequency / Central_Frequency)
      Variance = Numerator/Denominator
    ELSE
      Variance = 0.0
    ENDIF

    ! --------------------------------------------------------------------------
    ! Count number of non-zero variance estimates.
    ! --------------------------------------------------------------------------

    IF (Variance .GT. 0.0 .AND. Second_Point .LT. Center_Point) THEN
      Non_Zero_Count_LHS = Non_Zero_Count_LHS + 1
      Sum_Variance_LHS = Sum_Variance_LHS + Variance
    ELSE IF (Variance .GT. 0.0 .AND. Second_Point .GT. Center_Point) THEN
      Non_Zero_Count_RHS = Non_Zero_Count_RHS + 1
      Sum_Variance_RHS = Sum_Variance_RHS + Variance
    ENDIF

  END DO point_loop

  ! ----------------------------------------------------------------------------
  ! Determine average variance for each half curve.
  ! ----------------------------------------------------------------------------

  IF (Non_Zero_Count_LHS .GT. 0) Avg_Variance_LHS = Sum_Variance_LHS / &
                                              REAL(Non_Zero_Count_LHS)
  IF (Non_Zero_Count_RHS .GT. 0) Avg_Variance_RHS = Sum_Variance_RHS / &
                                              REAL(Non_Zero_Count_RHS)

END SUBROUTINE Gaussian_Fit
! ------------------------------------------------------------------------------
!
! Name:
!   Gaussian_Curve
!
! Function:
!  Compute the full Gaussian curve from two half-Gaussian estimates.
!
! Description:
!   - Determine interval over which the Gaussian is to be computed (-/+ 5 S.D.)
!   - Evaluate Gaussian function
!
! Reference:
!   None
!
! Calling sequence:
!   CALL Gaussian_Curve
!
! Inputs:
!   Max_Frequency_Slot   Center point of Gaussian curve.
!   Max_Frequency        Frequency at center of histogram.
!   Variance_LHS         Variance estimate for left side of curve.
!   Variance_RHS         Variance estimate for right side of curve.
!
! Outputs:
!   Gaussian_Histogram   Full Gaussian histogram.
!
! Dependencies:
!
! Restrictions:
!
! History:
!   5/2007 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE Gaussian_Curve(Max_Frequency_Slot, Max_Frequency, Variance_LHS, &
                                         Variance_RHS, Gaussian_Histogram)

  INTEGER(LONG), INTENT(IN) :: Max_Frequency, Max_Frequency_Slot
  INTEGER(LONG), INTENT(OUT), DIMENSION(:) :: Gaussian_Histogram
  INTEGER(LONG) :: Slot_High, Slot_Low, Slot_Number

  REAL(SINGLE), INTENT(IN) :: Variance_LHS, Variance_RHS
  REAL(SINGLE) :: Radiance_Center_Avg, Radiance_Center_High, Radiance_Center_Low
  REAL(SINGLE) :: Radiance_Point_Avg, Radiance_Point_High, Radiance_Point_Low
  REAL(SINGLE) :: Gaussian_Exponent, Variance

  ! ----------------------------------------------------------------------------
  ! Initialize variables
  ! ----------------------------------------------------------------------------

  Radiance_Center_Low = 0.0
  Radiance_Center_High = 0.0
  Radiance_Center_Avg = 0.0
  Radiance_Point_Low = 0.0
  Radiance_Point_High = 0.0
  Radiance_Point_Avg = 0.0
  Gaussian_Histogram(1:NUMBER_OF_SLOTS) = 0

  Radiance_Center_Low = MIN_HISTOGRAM_VALUE + REAL(Max_Frequency_Slot - 1) &
                                                         * INTERVAL_LENGTH
  Radiance_Center_High = Radiance_Center_Low + INTERVAL_LENGTH
  Radiance_Center_Avg = (Radiance_Center_Low + Radiance_Center_High) / 2.0

  ! ----------------------------------------------------------------------------
  ! Set interval as -/+ 5 standard deviation units
  ! ----------------------------------------------------------------------------

  Slot_Low = Max_Frequency_Slot - NINT(5.0 * SQRT(Variance_LHS))
  Slot_High = Max_Frequency_Slot + NINT(5.0 * SQRT(Variance_RHS))
  IF (Slot_Low .LT. 1) Slot_Low = 1
  IF (Slot_High .GT. NUMBER_OF_SLOTS) Slot_High = NUMBER_OF_SLOTS

  slot_loop: DO Slot_Number = Slot_Low, Slot_High

    Gaussian_Exponent = MISSING_VALUE_REAL4

    IF (Slot_Number .LE. Max_Frequency_Slot) Variance = Variance_LHS
    IF (Slot_Number .GT. Max_Frequency_Slot) Variance = Variance_RHS

    Radiance_Point_Low = MIN_HISTOGRAM_VALUE + REAL(Slot_Number - 1) * &
                                                       INTERVAL_LENGTH
    Radiance_Point_High = Radiance_Point_Low + INTERVAL_LENGTH
    Radiance_Point_Avg = (Radiance_Point_Low + Radiance_Point_High) / 2.0
    IF (Variance .GT. 0.0) Gaussian_Exponent = -1.0 * (Radiance_Point_Avg &
                             - Radiance_Center_Avg)**2 / (2.0 * Variance)

    ! --------------------------------------------------------------------------
    ! Evaluate Gaussian function
    ! --------------------------------------------------------------------------

    IF (Gaussian_Exponent .LT. LOWER_GAUSSIAN_EXPONENT) THEN
      Gaussian_Histogram(Slot_Number) = 0
    ELSE
      Gaussian_Histogram(Slot_Number) = INT(Max_Frequency * &
                                    EXP(Gaussian_Exponent))
    ENDIF

  END DO slot_loop

END SUBROUTINE Gaussian_Curve

!-------------------------------------------------------------------------------
!
! Name:
!   Gross_Error_Checks
!
! Function:
!   Apply gross error checks to AMV and flag suspect winds
!
! Description:
!   Check AMV wind components for large accelerations and apply channel-specific
!   tests using the forecast wind estimate.
!
! Reference:
!
! Inputs:
!   U1               First vector u-component
!   U2               Second vector u-component
!   V1               First vector v-component
!   V2               Second vector v-component
!   U_Fcst           Forecast u-component
!   V_Fcst           Forecast v-component
!   Median_P         Median pressure assignment of AMV
!
! Output:
!   QC_Flag          Quality flag. A value greater than 0 indicates the AMV
!                    has failed some test.
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!   02/09/2009 - Wayne Bresky - Created with input from Chris Velden
!                               and Steve Wanzong of CIMSS
!   11/22/2021 - Wayne Bresky - Implement pressure cutoff thresholds for various 
!                               cloud types
!
!-------------------------------------------------------------------------------


SUBROUTINE Gross_Error_Checks(NavType, U1, U2, V1, V2, U_Fcst, V_Fcst, Median_PW, DomCldType, QC_Flag)

  ! AMV quality control flag
  INTEGER(LONG), INTENT(INOUT) :: QC_Flag
  ! dominant cloud type of target scene
  INTEGER(LONG), INTENT(IN) :: DomCldType

  ! AMV vector components
  REAL(SINGLE), INTENT(IN) :: U1
  REAL(SINGLE), INTENT(IN) :: U2
  REAL(SINGLE), INTENT(IN) :: V1
  REAL(SINGLE), INTENT(IN) :: V2
  REAL(SINGLE), INTENT(IN) :: Median_PW
  REAL(SINGLE) ::  U_Avg
  REAL(SINGLE) ::  V_Avg
  REAL(SINGLE) ::  Speed_AMV
  REAL(SINGLE) ::  Dir_AMV

  REAL(SINGLE)  Vec_Diff
  REAL(SINGLE)  Vec_Diff_Cutoff

  ! Forecast vector components
  REAL(SINGLE), INTENT(IN) :: U_Fcst
  REAL(SINGLE), INTENT(IN) :: V_Fcst


  ! Max acceleration threshold
  REAL(SINGLE) :: Max_Acceleration

  ! Component accelerations
  REAL(SINGLE) :: U_Acceleration
  REAL(SINGLE) :: V_Acceleration


  CHARACTER(LEN=*), INTENT(IN) :: NavType

  ! ----------------------------------------------------------------------------
  ! Check wind component accelerations and flag those exceeding max limit
  ! ----------------------------------------------------------------------------

  Max_Acceleration = ACCEL_THRESHOLD
  IF (Channel .EQ. NATIVE_ABI_CHN2) Max_Acceleration = ACCEL_THRESHOLD_VZ

  U_Acceleration = ABS(U2 - U1)
  V_Acceleration = ABS(V2 - V1)

  IF (U_Acceleration .GT. Max_Acceleration .AND. &
      V_Acceleration .LE. Max_Acceleration) THEN

    QC_Flag = U_ACCELERATION_FAILURE
    RETURN

  ELSE IF (U_Acceleration .LE. Max_Acceleration .AND. &
           V_Acceleration .GT. Max_Acceleration) THEN

    QC_Flag = V_ACCELERATION_FAILURE
    RETURN

  ELSE IF (U_Acceleration .GT. Max_Acceleration .AND. &
           V_Acceleration .GT. Max_Acceleration) THEN

    QC_Flag = UV_ACCELERATION_FAILURE
    RETURN

  ENDIF

  ! ----------------------------------------------------------------------------
  ! Compute average wind.
  ! ----------------------------------------------------------------------------

  U_Avg = (U1 + U2) / 2.0
  V_Avg = (V1 + V2) / 2.0

  CALL UV2Spd(Speed_AMV, Dir_AMV, U_Avg, V_Avg)

  ! ----------------------------------------------------------------------------
  ! Flag slow winds.
  ! ----------------------------------------------------------------------------

  IF (Speed_AMV .LT. 3.0) THEN

    QC_Flag = SLOW_WIND_FAILURE
    RETURN

  ENDIF

  ! WCB
  ! Implement pressure limit on thick ice, cirrus, supercooled and liquid cloud scenes and
  ! discard all mixed cloud types. Applied only to cloud-track winds.
  IF (TRIM(NavType) .NE. 'PS' ) THEN

    IF (CSWV_Flag .EQ. sym%NO) THEN

      IF (DomCldType .EQ. sym%MIXED_TYPE) THEN

        QC_Flag = GROSS_ERROR_FAILURE
        RETURN

      ENDIF

      IF (Median_PW .GT. 500. .AND. (DomCldType .EQ. sym%TICE_TYPE .OR. DomCldType .EQ. sym%CIRRUS_TYPE)) THEN
        QC_Flag = GROSS_ERROR_FAILURE
        RETURN

      ENDIF

      IF (Median_PW .LT. 500. .AND. (DomCldType .EQ. sym%WATER_TYPE .OR. DomCldType .EQ. sym%SUPERCOOLED_TYPE)) THEN
        QC_Flag = GROSS_ERROR_FAILURE
        RETURN

      ENDIF

    ENDIF

  ENDIF

  ! ----------------------------------------------------------------------------
  ! Convert forecast u-/v-components to speed and direction
  ! ----------------------------------------------------------------------------

  !For additional possible tests:
  !CALL UV2Spd(Speed_Fcst, Dir_Fcst, U_Fcst, V_Fcst)
  !Speed_Diff = Speed_AMV - Speed_Fcst
  !Abs_Dir_Diff1 = ABS( Dir_AMV - Dir_Fcst )
  !Abs_Dir_Diff2 = ABS( ABS( (Dir_Fcst - Dir_AMV) ) - DEGREE360 )

  Vec_Diff = SQRT( (U_Avg - U_Fcst)**2 + (V_Avg - V_Fcst)**2 )

  ! ----------------------------------------------------------------------------
  ! Apply channel-specific gross error checks using the forecast wind.
  ! ----------------------------------------------------------------------------

  IF (Channel .EQ. NATIVE_ABI_CHN2) THEN

    Vec_Diff_Cutoff = VEC_DIFF_VZ

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN7) THEN

    Vec_Diff_Cutoff = VEC_DIFF_SWIR

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN8 .AND. CSWV_Flag .EQ. sym%NO) THEN

    Vec_Diff_Cutoff = VEC_DIFF_WVCT8

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN8 .AND. CSWV_Flag .EQ. sym%YES) THEN

    Vec_Diff_Cutoff = VEC_DIFF_WVCS8

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN9) THEN

    Vec_Diff_Cutoff = VEC_DIFF_WVCS9

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN10) THEN

    Vec_Diff_Cutoff = VEC_DIFF_WVCS10

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN14 .AND. TRIM(NavType) .NE. 'PS') THEN

    Vec_Diff_Cutoff = VEC_DIFF_LWIR

  ELSE IF (Channel .EQ. NATIVE_ABI_CHN14 .AND. TRIM(NavType) .EQ. 'PS') THEN

    Vec_Diff_Cutoff = VEC_DIFF_LWIR_POLAR

  ENDIF

  IF (Vec_Diff .GT. Vec_Diff_Cutoff) THEN

    QC_Flag = GROSS_ERROR_FAILURE
    RETURN

  ENDIF

  QC_Flag = GOOD_TARGET

END SUBROUTINE Gross_Error_Checks




!-------------------------------------------------------------------------------
!
! Name:
!   UV2Spd
!
! Function:
!   Converts the u- and v-components of a wind vector to speed and direction.
!
! Description:
!   This routine takes the u- and v-components of a wind vector and converts
!   them to speed and direction.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL UV2Spd
!
! Inputs:
!   U                        U-component wind (a west wind is positive)
!   V                        V-component wind (a south wind is positive)
!
! Outputs:
!   Spd                      Wind speed (m/s)
!   Dir                      Wind direction (degrees)
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 08/2009 - Wayne Bresky - Created
! ------------------------------------------------------------------------------
SUBROUTINE UV2Spd(Spd,Dir,U,V)

  REAL(SINGLE), PARAMETER :: RADDEG = DTOR**(-1)
  REAL(SINGLE) :: Spd
  REAL(SINGLE) :: Dir
  REAL(SINGLE), INTENT(IN) :: U
  REAL(SINGLE), INTENT(IN) :: V
  REAL(SINGLE) :: Ang

  IF (U .EQ. 0. .AND. V .EQ. 0.) THEN
    Spd = 0.
    Dir = 0.
    RETURN
  ENDIF

  Spd = SQRT(U*U + V*V)
  IF (U .EQ. 0. .AND. V .LT. 0.) THEN
    Dir = DEGREE0
  ELSE IF (V .EQ. 0. .AND. U .LT. 0.) THEN
    Dir = DEGREE90
  ELSE IF (U .EQ. 0. .AND. V .GT. 0.) THEN
    Dir = DEGREE180
  ELSE IF (V .EQ. 0. .AND. U .GT. 0.) THEN
    Dir = DEGREE270
  ELSE IF (U .LT. 0. .AND. V .LT. 0.) THEN
    Ang = ASIN(ABS(U)/Spd)
    Dir = Ang*RADDEG
  ELSE IF (U .LT. 0. .AND. V .GT. 0.) THEN
    Ang = ASIN(ABS(V)/Spd)
    Dir = Ang*RADDEG + DEGREE90
  ELSE IF (U .GT. 0. .AND. V .GT. 0.) THEN
    Ang = ASIN(ABS(U)/Spd)
    Dir = Ang*RADDEG + DEGREE180
  ELSE IF (U .GT. 0. .AND. V .LT. 0.) THEN
    Ang = ASIN(ABS(V)/Spd)
    Dir = Ang*RADDEG + DEGREE270
  ENDIF

  RETURN

END SUBROUTINE UV2Spd

!-------------------------------------------------------------------------------
!
! Name:
!   Spd2UV
!
! Function:
!   Converts speed and direction to u- and v-components.
!
! Description:
!   This routine decomposes a wind vector into u- and v-components.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Spd2UV
!
! Inputs:
!   Spd                      Wind speed (m/s)
!   Dir                      Wind direction (degrees)
!
! Outputs:
!   U                        U-component wind (a west wind is positive)
!   V                        V-component wind (a south wind is positive)
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 08/2009 - Wayne Bresky - Created
! ------------------------------------------------------------------------------
SUBROUTINE Spd2UV(Spd,Dir,U,V)

  REAL(SINGLE) :: Spd
  REAL(SINGLE) :: Dir
  REAL(SINGLE) :: U
  REAL(SINGLE) :: V
  REAL(SINGLE) :: Ang

  IF (Spd .EQ. 0.) THEN
    U = 0.
    V = 0.
    RETURN
  ENDIF

  IF (Dir .EQ. DEGREE0 .OR. Dir .EQ. DEGREE360) THEN
    U = 0.
    V = - Spd
  ELSE IF (Dir .EQ. DEGREE90) THEN
    U = -Spd
    V = 0.
  ELSE IF (Dir .EQ. DEGREE180) THEN
    U = 0.
    V = Spd
  ELSE IF (Dir .EQ. DEGREE270) THEN
    U = Spd
    Dir = 0.
  ELSE IF (Dir .GT. DEGREE0 .AND. Dir .LT. DEGREE90) THEN
    Ang = Dir*DTOR
    U = - Spd*SIN(Ang)
    V = - Spd*COS(Ang)
  ELSE IF (Dir .GT. DEGREE90 .AND. Dir .LT. DEGREE180) THEN
    Ang = (Dir- DEGREE90)*DTOR
    U = - Spd*COS(Ang)
    V =   Spd*SIN(Ang)
  ELSE IF (Dir .GT. DEGREE180 .AND. Dir .LT. DEGREE270) THEN
    Ang = (Dir - DEGREE180)*DTOR
    U =   Spd*SIN(Ang)
    V =   Spd*COS(Ang)
  ELSE IF (Dir .GT. DEGREE270 .AND. Dir .LT. DEGREE360) THEN
    Ang = (Dir - DEGREE270)*DTOR
    U =   Spd*COS(Ang)
    V = - Spd*SIN(Ang)
  ENDIF

  RETURN

END SUBROUTINE Spd2UV

!-------------------------------------------------------------------------------
!
! Name:
!   Compute_Zenith_Angle
!
! Function:
!   Computes the satellite zenith angle for a given earth location
!
! Description:
!   Computes the satellite zenith angle for a given latitude and longitude.
!   The satellite location is also required.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Compute_Zenith_Angle
!
! Inputs:
!   glon                     Satellite longitude
!   glat                     Satellite latitude
!   xlon                     Longitude of pixel
!   xlat                     Latitude of pixel
!
! Outputs:
!   zenith                   Satellite zenith angle
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 08/2009 - Wayne Bresky - Created
! ------------------------------------------------------------------------------
SUBROUTINE Compute_Zenith_Angle(glon, glat, xlon, xlat, zenith)

  REAL(SINGLE), INTENT(IN)  :: glon
  REAL(SINGLE), INTENT(IN)  :: glat
  REAL(SINGLE), INTENT(IN)  :: xlon
  REAL(SINGLE), INTENT(IN)  :: xlat
  REAL(SINGLE), INTENT(OUT) :: zenith

  ! Local variables
  REAL(DOUBLE) :: satlon
  REAL(DOUBLE) :: satlat
  REAL(DOUBLE) :: lat
  REAL(DOUBLE) :: lon
  REAL(DOUBLE) :: beta
  REAL(DOUBLE) :: sin_beta

  zenith = MISSING_VALUE_REAL4

  satlon = glon
  satlat = glat


  lon = (xlon - satlon) * DTOR   ! in radians
  lat = (xlat - satlat) * DTOR   ! in radians

  beta = ACOS( COS(lat) * COS(lon) )
  sin_beta = SIN(beta)

  ! zenith angle
  ! R**2 + (R+h)**2 = 1.808e09
  ! 2R(R+h) = 5.37e08
  ! R is the radius of the earth and h is the height of the satellite

  zenith = ASIN(MAX(-1.0_real8, MIN(1.0_real8, 42164.0* sin_beta / &
                                 SQRT(1.808e09 - 5.3725e08*COS(beta)))))
  zenith = zenith / DTOR

END SUBROUTINE Compute_Zenith_Angle

!-------------------------------------------------------------------------------
!
! Name:
!   Fill_Search_Box
!
! Function:
!   Loads data from the temporal data buffer into search box.
!
! Description:
!   This routine loads brightness temperature data from the buffer into the
!   search box. It requires the input of the target center point and the
!   forecast displacement in pixels.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Fill_Search_Box
!
! Inputs:
!   Target_Center_Element      Image element of target box center point
!   Target_Center_Line         Image line of target box center point
!   Delta_X                    Forecast element displacement
!   Delta_Y                    Forecast line displacement
!   Search_Half_Width_Line     Half the number of lines in the search region
!   Search_Half_Width_Element  Half the number of elements in the search region
!   Lines_In_Buffer            Total number of lines in the data buffer
!   Search_Buffer              Temporal brightness temperature data
!
! Outputs:
!   Search_BrtTemp_Values      Search box brightness temperature values
!   QC_Flag                    Status flag
!
! Dependencies:
!   Requires the temporal data buffers (first and third images) be full.
!
! Restrictions:
!   None.
!
! History:
! 08/2009 - Wayne Bresky - Created
! ------------------------------------------------------------------------------
SUBROUTINE  Fill_Search_Box(Elems, Target_Center_Element, Target_Center_Line, &
                             Delta_X, Delta_Y, Search_Half_Width_Line, &
                           Search_Half_Width_Element, Lines_In_Buffer, &
                        Search_Buffer, Search_BrtTemp_Values, QC_Flag)

  ! buffer coordinates of target box center
  INTEGER(LONG), INTENT(IN) :: Target_Center_Line
  INTEGER(LONG), INTENT(IN) :: Target_Center_Element
  ! buffer coordinates of search box center
  INTEGER(LONG) Search_Center_Line
  INTEGER(LONG) Search_Center_Element
  INTEGER(LONG), INTENT(IN) :: Search_Half_Width_Element
  INTEGER(LONG), INTENT(IN) :: Search_Half_Width_Line
  INTEGER(LONG), INTENT(IN) :: Lines_In_Buffer
  INTEGER(LONG), INTENT(INOUT) :: QC_Flag
  INTEGER(LONG) Search_Start_Line
  INTEGER(LONG) Search_End_Line
  INTEGER(LONG) Search_Start_Element
  INTEGER(LONG) Search_End_Element
  INTEGER(LONG) :: Search_Line_Size
  INTEGER(LONG) :: Search_Element_Size

  ! distance between target box center and search box center
  REAL(SINGLE), INTENT(IN) :: Delta_X
  REAL(SINGLE), INTENT(IN) :: Delta_Y
  REAL(SINGLE), DIMENSION(:,:), INTENT(IN) :: Search_Buffer
  REAL(SINGLE), DIMENSION(:,:), INTENT(INOUT) :: Search_BrtTemp_Values

  INTEGER(LONG), INTENT(IN) :: Elems


  ! Compute nominal bounds of search region in buffer

  Search_Center_Element = Target_Center_Element + INT(Delta_X)
  Search_Center_Line = Target_Center_Line + INT(Delta_Y)

  Search_Start_Line = Search_Center_Line - Search_Half_Width_Line
  Search_End_Line = Search_Center_Line + Search_Half_Width_Line
  Search_Start_Element = Search_Center_Element - Search_Half_Width_Element
  Search_End_Element = Search_Center_Element + Search_Half_Width_Element


  ! Center search on expected location. Do not exceed bounds of buffer.

  Search_Start_Line = MAX(1, Search_Start_Line)
  Search_End_Line = MIN(Lines_In_Buffer, Search_End_Line)
  Search_Line_Size = Search_End_Line - Search_Start_Line + 1
  Search_Start_Element = MAX(1, Search_Start_Element)
  Search_End_Element = MIN(Elems, Search_End_Element)
  Search_Element_Size = Search_End_Element - Search_Start_Element + 1

  ! Fail if we do not have a full-sized search box.

  IF (Search_Line_Size .LT. (Search_Half_Width_Line * 2 + 1) .OR. &
    Search_Element_Size .LT. (Search_Half_Width_Element * 2 + 1)) THEN

    QC_Flag = SEARCH_NOT_FULL_FAILURE
    RETURN

  ELSE

    Search_BrtTemp_Values(1:Search_Element_Size, 1:Search_Line_Size) = &
                Search_Buffer(Search_Start_Element:Search_End_Element, &
                                    Search_Start_Line:Search_End_Line)
  ENDIF

  !IF (ANY(Search_BrtTemp_Values == MISSING_VALUE_REAL4)) QC_Flag = MISSING_DATA_SEARCH
  IF (ALL(Search_BrtTemp_Values == MISSING_VALUE_REAL4)) QC_Flag = MISSING_DATA_SEARCH

END SUBROUTINE Fill_Search_Box

!-------------------------------------------------------------------------------
!
! Name:
!   Compute_End_Points
!
! Function:
!   Computes the end point of a vector in earth coordinates
!
! Description:
!   Given a starting point and displacement expressed in pixel coordinates,
!   this routine computes the end point in earth coordinates.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Compute_End_Points
!
! Inputs:
!   Search_Element           Image element of starting point
!   Element_Disp             Element displacement
!   Search_Line              Image line of starting point
!   Line_Disp                Line displacement
!   Sat_ID                   Satellite ID
!
! Outputs:
!   Lat_Match                Latitude of vector end point
!   Lon_Match                Longitude of vector end point
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 08/2009 - Wayne Bresky - Created
! 01/2012 - Wayne Bresky - Added PS map projection
! ------------------------------------------------------------------------------
SUBROUTINE Compute_End_Points(Nav, Search_Element, Element_Disp, Search_Line, &
                              Line_Disp, Lat_Match, Lon_Match)
  TYPE(Nav_t), INTENT(IN):: Nav

  REAL(SINGLE), INTENT(IN) :: Search_Element
  REAL(SINGLE), INTENT(IN) :: Search_Line
  REAL(SINGLE), INTENT(IN) :: Element_Disp
  REAL(SINGLE), INTENT(IN) :: Line_Disp
  REAL(SINGLE) Match_Element
  REAL(SINGLE) Match_Line
  REAL(SINGLE), INTENT(OUT) :: Lat_Match
  REAL(SINGLE), INTENT(OUT) :: Lon_Match

  REAL(DOUBLE) :: Latitude
  REAL(DOUBLE) :: Longitude




  ! Compute the endpoints of both vectors in image coordinates and convert
  ! to earth location.

  SELECT CASE (Nav%Sensor_Series)

  CASE (SERIES_MSG_SEVIRI)

    ! reverse sign of element displacement to account for EUMETSAT
    ! coordinate system

    Match_Element = Search_Element - Element_Disp
    Match_Line = Search_Line - Line_Disp
    CALL pixcoord2geocoord_dp(Nav, Match_Element, Match_Line, Latitude, Longitude)

    Lat_Match = REAL(Latitude)
    Lon_Match = REAL(Longitude)

!  CASE (SERIES_GOES_N_IMG)
!
!    Match_Element = Search_Element + Element_Disp
!    Match_Line = Search_Line + Line_Disp
!    CALL pixcoord2geocoord_gvar(Match_Element, Match_Line, &
!                          sat_xstride, Lat_Match, Lon_Match)

  CASE (SERIES_GOESR_ABI, SERIES_HIMAWARI_AHI)

    Match_Element = Search_Element + Element_Disp
    Match_Line = Search_Line + Line_Disp
    CALL pixcoord2geocoord_abi_real(Nav, Match_Element, Match_Line, &
                                         Latitude, Longitude)

    Lat_Match = REAL(Latitude)
    Lon_Match = REAL(Longitude)

    ! Because JMA sets their longitudes from 0 to 360, and
    ! we want 180 to -180, we need to switch here for longitudes
    ! that are greater than 180.

    IF (Lon_Match > 180.0) Lon_Match = Lon_Match - 360.0

  END SELECT

  ! polar winds test
  IF (TRIM(Nav%NavType) .EQ. 'PS' ) THEN

    Match_Element = Search_Element + Element_Disp
    Match_Line = Search_Line + Line_Disp
    CALL pixcoord2geocoord_ps(Nav, Match_Element, Match_Line, Lat_Match, Lon_Match)

    Lon_Match = -Lon_Match

  ENDIF

END SUBROUTINE Compute_End_Points

!-------------------------------------------------------------------------------
!
! Name:
!   Pixel_Factor
!
! Function:
!   Computes the ratio of E-W distance to N-S distance for a pixel
!
! Description:
!   Computes the ratio of x-distance to y-distance for a given pixel location.
!   This ratio is useful for normalizing a displacement expressed in elements
!   to a displacement expressed in lines without having to repeatedly transform
!   between pixel coordinates and earth coordinates.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Pixel_Factor
!
! Inputs:
!   Sat_ID                   Satellite ID
!   Point1_Element           Image element
!   Point1_Line              Image Line
!
! Outputs:
!   Factor                   Ratio of E-W distance to N-S distance
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 1/2010 - Wayne Bresky - Created
! ------------------------------------------------------------------------------
SUBROUTINE Pixel_Factor(Nav, Point1_Element, Point1_Line, Factor)
  CHARACTER(LEN=*), PARAMETER :: ROUTINE_NAME = "Pixel_Factor"
  TYPE(Nav_t), INTENT(IN):: Nav
  REAL(SINGLE), INTENT(IN) :: Point1_Element
  REAL(SINGLE), INTENT(IN) :: Point1_Line
  REAL(SINGLE), INTENT(OUT) :: Factor
  REAL(DOUBLE) Point1_Lat
  REAL(DOUBLE) Point1_Lon
  REAL(DOUBLE) Point2_Lat
  REAL(DOUBLE) Point2_Lon
  REAL(DOUBLE) Point3_Lat
  REAL(DOUBLE) Point3_Lon
  REAL(SINGLE) Pt1_Lat_Single
  REAL(SINGLE) Pt1_Lon_Single
  REAL(SINGLE) Pt2_Lat_Single
  REAL(SINGLE) Pt2_Lon_Single
  REAL(SINGLE) Pt3_Lat_Single
  REAL(SINGLE) Pt3_Lon_Single
  REAL(SINGLE) zlat1
  REAL(SINGLE) zlon1
  REAL(SINGLE) zlat2
  REAL(SINGLE) zlon2
  REAL(SINGLE) clat1
  REAL(SINGLE) clon1
  REAL(SINGLE) slat1
  REAL(SINGLE) slon1
  REAL(SINGLE) clat2
  REAL(SINGLE) clon2
  REAL(SINGLE) slat2
  REAL(SINGLE) slon2
  REAL(SINGLE) xx
  REAL(SINGLE) yy
  REAL(SINGLE) zz
  REAL(SINGLE) ARCL    !arclength
  REAL(SINGLE) DIST
  REAL(SINGLE) XDIST
  REAL(SINGLE) YDIST

  REAL(SINGLE), PARAMETER :: D2RAD = 0.017453292
  REAL(SINGLE), PARAMETER :: R = 6.378E3




  SELECT CASE(Nav%Sensor_Series)

  CASE (SERIES_MSG_SEVIRI)

    CALL pixcoord2geocoord_dp(Nav, Point1_Element, Point1_Line, Point1_Lat, Point1_Lon)

    CALL pixcoord2geocoord_dp(Nav, Point1_Element+1., Point1_Line, Point2_Lat, Point2_Lon)

    CALL pixcoord2geocoord_dp(Nav, Point1_Element, Point1_Line+1., Point3_Lat, Point3_Lon)

    Pt1_Lat_Single = REAL(Point1_Lat)
    Pt1_Lon_Single = REAL(Point1_Lon)
    Pt2_Lat_Single = REAL(Point2_Lat)
    Pt2_Lon_Single = REAL(Point2_Lon)
    Pt3_Lat_Single = REAL(Point3_Lat)
    Pt3_Lon_Single = REAL(Point3_Lon)

!  CASE (SERIES_GOES_N_IMG)
!
!    CALL pixcoord2geocoord_gvar(Point1_Element, Point1_Line, sat_xstride, &
!                                 Pt1_Lat_Single, Pt1_Lon_Single)
!    CALL pixcoord2geocoord_gvar(Point1_Element+1., Point1_Line, sat_xstride, &
!                                    Pt2_Lat_Single, Pt2_Lon_Single)
!    CALL pixcoord2geocoord_gvar(Point1_Element, Point1_Line+1., sat_xstride, &
!                                    Pt3_Lat_Single, Pt3_Lon_Single)

  CASE (SERIES_GOESR_ABI)

    CALL pixcoord2geocoord_abi_real(Nav, Point1_Element, Point1_Line, &
                                    Point1_Lat, Point1_Lon)
    CALL pixcoord2geocoord_abi_real(Nav, Point1_Element+1., Point1_Line, &
                                    Point2_Lat, Point2_Lon)
    CALL pixcoord2geocoord_abi_real(Nav, Point1_Element, Point1_Line+1., &
                                    Point3_Lat, Point3_Lon)

    Pt1_Lat_Single = REAL(Point1_Lat)
    Pt1_Lon_Single = REAL(Point1_Lon)
    Pt2_Lat_Single = REAL(Point2_Lat)
    Pt2_Lon_Single = REAL(Point2_Lon)
    Pt3_Lat_Single = REAL(Point3_Lat)
    Pt3_Lon_Single = REAL(Point3_Lon)

  CASE (SERIES_HIMAWARI_AHI)

    CALL pixcoord2geocoord_abi_real(Nav, Point1_Element, Point1_Line, &
                                    Point1_Lat, Point1_Lon)
    Pt1_Lat_Single = REAL(Point1_Lat)
    Pt1_Lon_Single = REAL(Point1_Lon)
    IF (Pt1_Lon_Single > 180.0) Pt1_Lon_Single = Pt1_Lon_Single - 360.0

    CALL pixcoord2geocoord_abi_real(Nav, Point1_Element+1., Point1_Line, &
                                    Point2_Lat, Point2_Lon)

    Pt2_Lat_Single = REAL(Point2_Lat)
    Pt2_Lon_Single = REAL(Point2_Lon)
    IF (Pt2_Lon_Single > 180.0) Pt2_Lon_Single = Pt2_Lon_Single - 360.0

    CALL pixcoord2geocoord_abi_real(Nav, Point1_Element, Point1_Line+1., &
                                    Point3_Lat, Point3_Lon)

    Pt3_Lat_Single = REAL(Point3_Lat)
    Pt3_Lon_Single = REAL(Point3_Lon)
    IF (Pt3_Lon_Single > 180.0) Pt3_Lon_Single = Pt3_Lon_Single - 360.0


  END SELECT

  ! polar winds test
  IF (TRIM(Nav%NavType) .EQ. 'PS') THEN

    CALL pixcoord2geocoord_ps(Nav, Point1_Element, Point1_Line, Pt1_Lat_Single, Pt1_Lon_Single)
    CALL pixcoord2geocoord_ps(Nav, Point1_Element+1., Point1_Line, Pt2_Lat_Single, Pt2_Lon_Single)
    CALL pixcoord2geocoord_ps(Nav, Point1_Element, Point1_Line+1., Pt3_Lat_Single, Pt3_Lon_Single)

  END IF

  ! Calculate earth distance between lines and elements
  ! This code comes from McIDAS subroutine m0distnc.for

  ! Convert degrees to radians
  zlat1 = Pt1_Lat_Single * D2RAD
  zlon1 = Pt1_Lon_Single * D2RAD
  zlat2 = Pt2_Lat_Single * D2RAD
  zlon2 = Pt2_Lon_Single * D2RAD

  ! Compute sines and cosines
  clat1=COS(zlat1)
  clon1=COS(zlon1)
  slat1=SIN(zlat1)
  slon1=SIN(zlon1)
  clat2=COS(zlat2)
  clon2=COS(zlon2)
  slat2=SIN(zlat2)
  slon2=SIN(zlon2)

  ! set up trig to get distance

  xx=clat2*clon2-clat1*clon1
  yy=clat2*slon2-clat1*slon1
  zz=slat2-slat1
  DIST=SQRT(xx*xx+yy*yy+zz*zz)
  ARCL=2.*ASIN(DIST/2.)*R
  XDIST=ARCL

  zlat2 = Pt3_Lat_Single * D2RAD
  zlon2 = Pt3_Lon_Single * D2RAD

  clat2=COS(zlat2)
  clon2=COS(zlon2)
  slat2=SIN(zlat2)
  slon2=SIN(zlon2)

  xx=clat2*clon2-clat1*clon1
  yy=clat2*slon2-clat1*slon1
  zz=slat2-slat1
  DIST=SQRT(xx*xx+yy*yy+zz*zz)
  ARCL=2.*ASIN(DIST/2.)*R
  YDIST=ARCL

  Factor = XDIST/YDIST

END SUBROUTINE Pixel_Factor

!-------------------------------------------------------------------------------
!
! Name:
!   Check_Height_Assignment
!
! Function:
!   Checks the AMV height assignment for errors.
!
! Description:
!   Checks if the AMV height assignment is missing or based on a
!   small (less than 2%) sample of cloudy pixels. The subroutine
!   also applies channel-specific pressure thresholds.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Check_Height_Assignment
!
! Inputs:
!   Point_Index              Number of pixels used to assign height
!   Median_Press             Pressure assignment of AMV
!   Press_Profile            Forecast pressure levels
!
! Outputs:
!   QC_Flag                  Quality control flag
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 04/2010 - Wayne Bresky - Created
! ------------------------------------------------------------------------------
SUBROUTINE Check_Height_Assignment(Point_Index, Median_Press, &
                                      Press_Profile, QC_Flag)

  INTEGER(LONG), INTENT(INOUT) :: QC_Flag
  INTEGER(LONG), INTENT(IN) :: Point_Index
  INTEGER(LONG) Min_Sample_Size
  REAL(SINGLE), INTENT(IN) :: Median_Press
  REAL(SINGLE), DIMENSION(:), INTENT(IN) :: Press_Profile
  REAL(SINGLE) PMIN
  REAL(SINGLE) PMAX

  ! First check if median pressure is missing or less than the value of
  ! the first pressure level. This applies to all bands.

  IF (Median_Press .EQ. MISSING_VALUE_REAL4 .OR. &
      Median_Press .LT. Press_Profile(1) ) THEN

    QC_Flag = BAD_PRESS_FAILURE
    RETURN

  ENDIF

  ! Check to see if we are processing clear-sky WV winds and set
  ! pressure thresholds accordingly

  IF (CSWV_Flag .EQ. sym%NO) THEN

    ! For cloud-track winds also check size of CTP sample.

    Min_Sample_Size = CEILING(REAL(Total_Points) * MIN_PERCENT_PRESS_VALUES)

    IF (Point_Index .LT. Min_Sample_Size) THEN

      QC_Flag = BAD_PRESS_FAILURE
      RETURN

    ENDIF

    SELECT CASE(Channel)

    ! set upper and lower limits for all cloud-track channels

    CASE(NATIVE_ABI_CHN2)
      !aab SWIR - Mods for VIIRS...use different threshold for VIIRS SWIR
      !           VIIRS SWIR is a reflectance value, so it's treated as a VIS channel
      !PMIN = PMIN_VISIBLE
      IF ( Sensor_Series == SERIES_JPSS_VIIRS ) THEN
        PMIN = PMIN_VISIBLE_VIIRS
      ELSE
        PMIN = PMIN_VISIBLE
      ENDIF
      PMAX = PMAX_ALL_WINDS

    CASE(NATIVE_ABI_CHN7)
      PMIN = PMIN_SWIR
      PMAX = PMAX_ALL_WINDS
    CASE(NATIVE_ABI_CHN8 : NATIVE_ABI_CHN10)
      PMIN = PMIN_CTWV
      IF ( Sensor_Series == SERIES_EOS_MODIS ) THEN
         PMAX = PMAX_CTWV_MODIS
      ELSE
         PMAX = PMAX_CTWV
      ENDIF
    CASE(NATIVE_ABI_CHN14)
      PMIN = PMIN_LWIR
      PMAX = PMAX_ALL_WINDS
    END SELECT

  ELSE IF (CSWV_Flag .EQ. sym%YES) THEN

    SELECT CASE(Channel)

    ! set upper and lower limits for clear-sky WV winds

    CASE(NATIVE_ABI_CHN8 : NATIVE_ABI_CHN9)
      PMIN = PMIN_CTWV
      PMAX = PMAX_ALL_WINDS
    CASE(NATIVE_ABI_CHN10)
      PMIN = PMIN_CSWV10
      PMAX = PMAX_CSWV10
    END SELECT

  ENDIF

  ! Check min/max limits

  IF (Median_Press .LT. PMIN .OR. Median_Press .GT. PMAX) THEN

    QC_Flag = PRESS_THRESHOLD_FAILURE
    RETURN

  ENDIF

END SUBROUTINE Check_Height_Assignment

!-------------------------------------------------------------------------------
!
! Name:
!   Interpolate_To_AMV_Press
!
! Function:
!   Interpolates the forecast wind to the AMV pressure.
!
! Description:
!   Performs a linear interpolation of the forecast wind to the AMV
!   pressure assignment.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Interpolate_To_AMV_Press
!
! Inputs:
!   Press_Levels             Forecast pressure levels
!   Fcst_Prof                Forecast profile data structure
!
! Input/Output:
!   Output                   Output data structure
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 04/2010 - Wayne Bresky - Created
! ------------------------------------------------------------------------------
SUBROUTINE Interpolate_To_AMV_Press(Output, Press_Levels, Fcst_Prof)

  TYPE(Output_Variables), INTENT(INOUT) :: Output
  TYPE(Forecast_Profiles), INTENT(IN) :: Fcst_Prof
  INTEGER(LONG), INTENT(IN) :: Press_Levels
  INTEGER(LONG) :: Level

  ! Locate nearest forecast level corresponding to median pressure.

  CALL Locate(Fcst_Prof%Pressure, Press_Levels, Output%Median_Press, Level)

  ! Interpolate forecast wind to median pressure level.

  IF ( (Level+1) .LE. Press_Levels) THEN

    Output%U_Fcst = ((Fcst_Prof%Pressure(Level+1) - Output%Median_Press) * &
                         Fcst_Prof%U_Wind(Level) +  (Output%Median_Press - &
                 Fcst_Prof%Pressure(Level)) * Fcst_Prof%U_Wind(Level+1)) / &
                (Fcst_Prof%Pressure(Level+1) - Fcst_Prof%Pressure(Level))

    Output%V_Fcst = ((Fcst_Prof%Pressure(Level+1) - Output%Median_Press) * &
                         Fcst_Prof%V_Wind(Level) +  (Output%Median_Press - &
                 Fcst_Prof%Pressure(Level)) * Fcst_Prof%V_Wind(Level+1)) / &
                (Fcst_Prof%Pressure(Level+1) - Fcst_Prof%Pressure(Level))

  ELSE

    Output%U_Fcst = Fcst_Prof%U_Wind(Level)
    Output%V_Fcst = Fcst_Prof%V_Wind(Level)

  ENDIF

END SUBROUTINE Interpolate_To_AMV_Press

!-------------------------------------------------------------------------------
!
! Name:
!   Assign Height
!
! Function:
!   Assigns a representative height (pressure) to the AMV.
!
! Description:
!   Assigns a representative height to the AMV by either finding
!   the median cloud-top pressure of a sample of cloudy pixels
!   or by finding the forecast level corresponding to a median
!   brightness temperature.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Assign Height
!
! Inputs:
!   Box_Data                 Target box data structure
!   Fcst_Prof                Forecast profile data structure
!
! Input/Output:
!   Output                   Output data structure of current target scene
!
! Outputs:
!   QC_Flag                  Quality control flag
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 04/2010 - Wayne Bresky - Created
! ------------------------------------------------------------------------------
SUBROUTINE Assign_Height(Press_Levels, Output, Box_Data, Fcst_Prof, QC_Flag)

  INTEGER(LONG), INTENT(IN) :: Press_Levels
  TYPE(Output_Variables), INTENT(INOUT) :: Output
  TYPE(Target_Box_Data), INTENT(INOUT) :: Box_Data
  TYPE(Forecast_Profiles), INTENT(IN) :: Fcst_Prof

  INTEGER(LONG), INTENT(INOUT) :: QC_Flag
  INTEGER(LONG) Median_Level

  REAL(SINGLE) a

  ! find median values of pressure, temperature, brightness temperature
  ! and cloud height

  CALL Sort_Shell_NP(Output%Point_Index, Box_Data%Pressure_Cold_Sample, &
                                                   Output%Median_Press)

  ! need to add cloud top temperature here

  CALL Sort_Shell_NP(Output%Point_Index, Box_Data%BrtTemp_Cold_Sample, &
                                                Output%Median_BrtTemp)

  CALL Sort_Shell_NP(Output%Point_Index, Box_Data%Cld_Hgt_Cold_Sample, &
                                                    Output%Median_Hgt)

  ! For the clear-sky water vapor winds find level in forecast profile
  ! corresponding to median BrtTemp.

  IF (CSWV_Flag .EQ. sym%YES .AND. Output%Median_BrtTemp .NE. MISSING_VALUE_REAL4) THEN

    CALL Locate(Fcst_Prof%Temp, Press_levels, Output%Median_BrtTemp, &
                                                       Median_Level)

    ! Check for level index of zero which can happen if median BT is
    ! colder than coldest temperature in (26-level) profile. May need
    ! to access higher resolution profile to correct this. For now we
    ! will flag the target as bad. If we don't flag the target the
    ! result is a pressure value of NaN.

    IF (Median_Level .EQ. 0) THEN

      QC_Flag = BAD_PRESS_FAILURE
      RETURN

    ENDIF

    ! Interpolate to find pressure associated with median BrtTemp

    a = LOG(Output%Median_BrtTemp/Fcst_Prof%Temp(Median_Level)) / &
                          LOG(Fcst_Prof%Temp(Median_Level + 1)  / &
                                 Fcst_Prof%Temp(Median_Level))

    IF (a .NE. 0.0) Output%Median_Press =  Fcst_Prof%Pressure(Median_Level)* &
                         EXP(a * LOG (Fcst_Prof%Pressure(Median_Level + 1) / &
                                        Fcst_Prof%Pressure(Median_Level)))

  ENDIF ! clear-sky WV check

END SUBROUTINE Assign_Height


SUBROUTINE IncrementElement(QC_Flag, Start_Element)
  CHARACTER(LEN=*), PARAMETER :: FUNC = "IncrementElement"
  INTEGER(LONG), INTENT(IN) :: QC_Flag
  INTEGER(LONG), INTENT(INOUT) :: Start_Element

  INTEGER(LONG) :: Element_Shift

  ! ------------------------------------------------------------------------
  ! Reset starting element for next box
  ! ------------------------------------------------------------------------

  ! Range of flag values (1-7) indicating target failures
  IF (QC_Flag .GT. 0 .AND. QC_Flag .LE. 7) THEN

    Element_Shift = Box_Size / 2

    ! if fixed spacing of output is desired use this setting
    !Element_Shift = NOMINAL_BOX_SIZE

  ELSE

    ! use nominal box size instead of target box size
    ! this results in slight overlap for 11 micron channel which uses
    ! box size of 19 (nominal size is 15)
    !Element_Shift = Box_Size

    Element_Shift = NOMINAL_BOX_SIZE

  ENDIF

  Start_Element = Start_Element + Element_Shift

END SUBROUTINE IncrementElement

!-------------------------------------------------------------------------------
!
! Name:
!   Find_Combined_Median
!
! Function:
!   Combines two samples and finds median value.
!
! Description:
!   This subroutine combines two samples of data (either cloud-top
!   pressure values or brightness temperature values) and finds the
!   median of the combined sample. The two samples come from the
!   nested tracking output of the reverse and forward vectors.
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Find_Combined_Median
!
! Inputs:
!   DBSCAN_Out               DBSCAN output data structure
!   LL_Inver_Flag            Low-level inversion flag
!   Fcst_Prof                Forecast profile data structure
!   Output                   Output data structure of current target scene
!
! Input/Output
!   QC_Flag                  Quality control flag
!
! Outputs
!   Combined_Median          Median pressure associated with combined sample
!   Combined_Median_Hgt      Median height associated with combined sample
!
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 04/2010 - Wayne Bresky - Created
! 08/2010 - Wayne Bresky - Added DBSCAN output structure
! ------------------------------------------------------------------------------
SUBROUTINE Find_Combined_Median(DBSCAN_Out, Fcst_Prof, Output, &
                               Combined_Median, Combined_Median_Hgt, QC_Flag)

  TYPE(DBSCAN_Output), INTENT(IN) :: DBSCAN_Out
  TYPE(Forecast_Profiles), INTENT(IN) :: Fcst_Prof
  TYPE(Output_Variables), INTENT(INOUT) :: Output
  REAL(SINGLE), INTENT(OUT) :: Combined_Median
  REAL(SINGLE), INTENT(OUT) :: Combined_Median_Hgt
  INTEGER(LONG), INTENT(INOUT) :: QC_Flag

  INTEGER(LONG) :: Combined_Sum
  INTEGER(LONG) :: Filtered_Sum
  INTEGER(LONG) :: Point


  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Combined_CTP_Sample
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Filtered_CTP_Sample
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Combined_Hgt_Sample
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Combined_CldTemp_Sample
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Filtered_CldTemp_Sample
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Combined_HgtErr_Sample
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Combined_TempErr_Sample
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: OD_Sample
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: ACHA_Cost_Sample
  REAL(SINGLE) Min_CTP
  REAL(SINGLE) Max_CTP
  REAL(SINGLE) Min_CTT
  REAL(SINGLE) Max_CTT
  REAL(SINGLE) :: Combined_Median_CldTemp
  REAL(SINGLE) :: Combined_Median_BrtTemp
  REAL(SINGLE) Combined_Median_HgtErr
  REAL(SINGLE) Combined_Median_TempErr
  REAL(SINGLE) SumOfValues
  REAL(SINGLE) SumOfWeights
  REAL(SINGLE) WAvg
  REAL(SINGLE) Mean_CTP
  REAL(SINGLE) Mean_CTT
  REAL(SINGLE) Variance_CTP
  REAL(SINGLE) SD_CTP
  REAL(SINGLE) Filtered_Median
  REAL(SINGLE) Filtered_Median_CldTemp
  ! optical depth
  REAL(SINGLE) Min_OD
  REAL(SINGLE) Max_OD
  REAL(SINGLE) Median_OD
  ! ACHA cost
  REAL(SINGLE) Min_Cost
  REAL(SINGLE) Max_Cost
  REAL(SINGLE) Median_Cost
  INTEGER(LONG) :: i
  INTEGER(LONG) :: counts
  INTEGER(LONG) :: maxcount
  INTEGER(LONG) :: dominant_phase
  INTEGER(LONG) :: dominant_type



  ! Initialize output value to missing

  Combined_Median = MISSING_VALUE_REAL4
  Filtered_Median = MISSING_VALUE_REAL4
  Combined_Median_CldTemp = MISSING_VALUE_REAL4
  Filtered_Median_CldTemp = MISSING_VALUE_REAL4
  Combined_Median_Hgt = MISSING_VALUE_REAL4
  Combined_Median_BrtTemp = MISSING_VALUE_REAL4
  Combined_Median_HgtErr = MISSING_VALUE_REAL4
  Combined_Median_TempErr = MISSING_VALUE_REAL4
  Median_OD = MISSING_VALUE_REAL4
  Median_Cost = MISSING_VALUE_REAL4

  ! added for Weighted CTP
  SumOfValues = 0.0
  SumOfWeights = 0.0
  WAvg = 0.0

  ! added for filtering of outliers
  Mean_CTP = 0.0
  Mean_CTT = 0.0
  Variance_CTP = 0.0
  SD_CTP = 0.0
  Filtered_Sum = 0

  ! combine the sample from the reverse vector with the sample
  ! from the forward vector

  Combined_Sum = DBSCAN_Out%Number_CTP_Sample1 + DBSCAN_Out%Number_CTP_Sample2

  counts = 0
  maxcount = 0
  dominant_phase = MISSING_VALUE_LONG
  DO I = sym%PHASE_MIN, sym%PHASE_MAX
     counts = COUNT(DBSCAN_Out%Phase_Sample1(1:DBSCAN_Out%Number_CTP_Sample1) == I)
     counts = counts + COUNT(DBSCAN_Out%Phase_Sample2(1:DBSCAN_Out%Number_CTP_Sample2) == I)
     IF(counts > maxcount) THEN
        dominant_phase = I
        maxcount = counts
     ENDIF
  ENDDO
  Output%Cloud_Phase = dominant_phase

  counts = 0
  maxcount = 0
  dominant_type = MISSING_VALUE_LONG
  DO I = sym%TYPE_MIN, sym%TYPE_MAX
     counts = COUNT(DBSCAN_Out%Type_Sample1(1:DBSCAN_Out%Number_CTP_Sample1) == I)
     counts = counts + COUNT(DBSCAN_Out%Type_Sample2(1:DBSCAN_Out%Number_CTP_Sample2) == I)
     IF(counts > maxcount) THEN
        dominant_type = I
        maxcount = counts
     ENDIF
  ENDDO
  Output%Cloud_Type = dominant_type


  ALLOCATE(Combined_CTP_Sample(Combined_Sum))
  ALLOCATE(Filtered_CTP_Sample(Combined_Sum))
  ALLOCATE(Combined_CldTemp_Sample(Combined_Sum))
  ALLOCATE(Filtered_CldTemp_Sample(Combined_Sum))
  ALLOCATE(Combined_Hgt_Sample(Combined_Sum))
  ALLOCATE(Combined_HgtErr_Sample(Combined_Sum))
  ALLOCATE(Combined_TempErr_Sample(Combined_Sum))
  ALLOCATE(OD_Sample(Combined_Sum))
  ALLOCATE(ACHA_Cost_Sample(Combined_Sum))

  Combined_CTP_Sample(1:DBSCAN_Out%Number_CTP_Sample1) = DBSCAN_Out% &
                    CTP_5x5_Sample1(1:DBSCAN_Out%Number_CTP_Sample1)

  Combined_CTP_Sample(DBSCAN_Out%Number_CTP_Sample1+1:) = DBSCAN_Out% &
                     CTP_5x5_Sample2(1:DBSCAN_Out%Number_CTP_Sample2)

  CALL Sort_Shell_NP(Combined_Sum, Combined_CTP_Sample, &
                                       Combined_Median)

  ! Find combined median cloud temperature

  Combined_CldTemp_Sample(1:DBSCAN_Out%Number_CTP_Sample1) = DBSCAN_Out% &
                    CldTemp_5x5_Sample1(1:DBSCAN_Out%Number_CTP_Sample1)

  Combined_CldTemp_Sample(DBSCAN_Out%Number_CTP_Sample1+1:) = DBSCAN_Out% &
                     CldTemp_5x5_Sample2(1:DBSCAN_Out%Number_CTP_Sample2)

  CALL Sort_Shell_NP(Combined_Sum, Combined_CldTemp_Sample, &
                                   Combined_Median_CldTemp)

  ! Find combined median cloud height

  Combined_Hgt_Sample(1:DBSCAN_Out%Number_CTP_Sample1) = DBSCAN_Out% &
                 Cld_Hgt_5x5_Sample1(1:DBSCAN_Out%Number_CTP_Sample1)

  Combined_Hgt_Sample(DBSCAN_Out%Number_CTP_Sample1+1:) = DBSCAN_Out% &
                 Cld_Hgt_5x5_Sample2(1:DBSCAN_Out%Number_CTP_Sample2)

  CALL Sort_Shell_NP(Combined_Sum, Combined_Hgt_Sample, &
                                   Combined_Median_Hgt)

  ! Find combined median height and temperature error

  Combined_HgtErr_Sample(1:DBSCAN_Out%Number_CTP_Sample1) = DBSCAN_Out% &
                   Press_Error_Sample1(1:DBSCAN_Out%Number_CTP_Sample1)

  Combined_HgtErr_Sample(DBSCAN_Out%Number_CTP_Sample1+1:) = DBSCAN_Out% &
                    Press_Error_Sample2(1:DBSCAN_Out%Number_CTP_Sample2)

  CALL Sort_Shell_NP(Combined_Sum, Combined_HgtErr_Sample, &
                                   Combined_Median_HgtErr)

  Combined_TempErr_Sample(1:DBSCAN_Out%Number_CTP_Sample1) = DBSCAN_Out% &
                     Temp_Error_Sample1(1:DBSCAN_Out%Number_CTP_Sample1)

  Combined_TempErr_Sample(DBSCAN_Out%Number_CTP_Sample1+1:) = DBSCAN_Out% &
                      Temp_Error_Sample2(1:DBSCAN_Out%Number_CTP_Sample2)

  CALL Sort_Shell_NP(Combined_Sum, Combined_TempErr_Sample, &
                                   Combined_Median_TempErr)

  ! Find median optical depth

  OD_Sample(1:DBSCAN_Out%Number_CTP_Sample1) = DBSCAN_Out% &
              Cld_OD_Sample1(1:DBSCAN_Out%Number_CTP_Sample1)

  OD_Sample(DBSCAN_Out%Number_CTP_Sample1+1:) = DBSCAN_Out% &
              Cld_OD_Sample2(1:DBSCAN_Out%Number_CTP_Sample2)

  CALL Sort_Shell_NP(Combined_Sum, OD_Sample, Median_OD)

  ! Find median cost

  ACHA_Cost_Sample(1:DBSCAN_Out%Number_CTP_Sample1) = DBSCAN_Out% &
              ACHA_Cost_Sample1(1:DBSCAN_Out%Number_CTP_Sample1)

  ACHA_Cost_Sample(DBSCAN_Out%Number_CTP_Sample1+1:) = DBSCAN_Out% &
              ACHA_Cost_Sample2(1:DBSCAN_Out%Number_CTP_Sample2)

  CALL Sort_Shell_NP(Combined_Sum, ACHA_Cost_Sample, Median_Cost)

  ! added for Weighted CTP (weighted by correlation)

  !DO Point=1,DBSCAN_Out%Number_CTP_Sample1

  !  ! use cold branch
  !  IF (DBSCAN_Out%CldTemp_5x5_Sample1(Point) .GT. Mean_CTT) CYCLE

  !  SumOfValues = SumOfValues + DBSCAN_Out%CTP_5x5_Sample1(Point) * DBSCAN_Out%CORR_5x5_Sample1(Point)
  !  SumOfWeights = SumOfWeights + DBSCAN_Out%CORR_5x5_Sample1(Point)

  !ENDDO

  !DO Point=1,DBSCAN_Out%Number_CTP_Sample2

  !  ! use cold branch
  !  IF (DBSCAN_Out%CldTemp_5x5_Sample2(Point) .GT. Mean_CTT) CYCLE

  !  SumOfValues = SumOfValues + DBSCAN_Out%CTP_5x5_Sample2(Point) * DBSCAN_Out%CORR_5x5_Sample2(Point)
  !  SumOfWeights = SumOfWeights + DBSCAN_Out%CORR_5x5_Sample2(Point)

  !ENDDO

  !WAvg = SumOfValues / SumOfWeights

  !IF (ISNAN(WAvg)) THEN

  !  WAvg = MISSING_VALUE_REAL4

  !ENDIF

  ! keep track of min and max values

  Min_CTP = MINVAL(Combined_CTP_Sample)
  Max_CTP = MAXVAL(Combined_CTP_Sample)
  Min_CTT = MINVAL(Combined_CldTemp_Sample)
  Max_CTT = MAXVAL(Combined_CldTemp_Sample)
  Min_OD = MINVAL(OD_Sample)
  Max_OD = MAXVAL(OD_Sample)
  Min_Cost = MINVAL(ACHA_Cost_Sample)
  Max_Cost = MAXVAL(ACHA_Cost_Sample)


  ! save values in output structure

  Output%Min_CTP = Min_CTP
  Output%Max_CTP = Max_CTP
  Output%Min_CTT = Min_CTT
  Output%Max_CTT = Max_CTT
  Output%Combined_Median = Combined_Median
  Output%Min_OD = Min_OD
  Output%Max_OD = Max_OD
  Output%Median_OD = Median_OD
  Output%Min_Cost = Min_Cost
  Output%Max_Cost = Max_Cost
  Output%Median_Cost = Median_Cost

  !use weighted pressure
  !Combined_Median = WAvg
  !Output%Combined_Median = WAvg
  !

  Output%Combined_Median_CldTemp = Combined_Median_CldTemp
  Output%Combined_Median_Hgt = Combined_Median_Hgt
  Output%Combined_Median_HgtErr = Combined_Median_HgtErr
  Output%Combined_Median_TempErr = Combined_Median_TempErr
  Output%Weighted_PW = WAvg

  ! update MedianPress
  Output%Median_Press = Combined_Median

  !deallocate combined sample arrays

  DEALLOCATE(Combined_CTP_Sample)
  DEALLOCATE(Filtered_CTP_Sample)
  DEALLOCATE(Combined_CldTemp_Sample)
  DEALLOCATE(Filtered_CldTemp_Sample)
  DEALLOCATE(Combined_Hgt_Sample)
  DEALLOCATE(Combined_HgtErr_Sample)
  DEALLOCATE(Combined_TempErr_Sample)
  DEALLOCATE(OD_Sample)
  DEALLOCATE(ACHA_Cost_Sample)

  CALL Check_Height_Assignment(Combined_Sum, Combined_Median, &
                                 Fcst_Prof%Pressure, QC_Flag)

END SUBROUTINE Find_Combined_Median

!-------------------------------------------------------------------------------
!
! Name:
!   Dominant_Cloud_Phase
!
! Function:
!   Find dominant cloud phase in target scene. Clear phase is ignored.
!
! Description:
!   Returns the dominant cloud phase in the target scene.  Will become
!   a diagnostic output variable.
!     CLEAR_PHASE       = 0     clear cloud mask
!     WATER_PHASE       = 1     liquid water clouds warmer than 273.16 K
!     SUPERCOOLED_PHASE = 2     liquid water clouds colder than 273.16 K
!     MIXED_PHASE       = 3     very cold liquid clouds with some ice
!     ICE_PHASE         = 4     ice clouds
!     UNKNOWN_PHASE     = 5     bad input data
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Dominant_Cloud_Phase
!
! Inputs:
!   Box_Data                 Target box data structure
!
! Input/Output:
!   Output                   Output data structure of current target scene
!
! Outputs:
!   Main_Cloud_Phase         Dominant cloud phase of target scene.
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 08/2010 - Steve Wanzong - Created
! ------------------------------------------------------------------------------
SUBROUTINE Dominant_Cloud_Phase(Box_Data, Main_Cloud_Phase)

  TYPE(Target_Box_Data), INTENT(IN) :: Box_Data

  INTEGER(LONG), INTENT(OUT) :: Main_Cloud_Phase

  INTEGER(LONG), PARAMETER :: PHASE_LOOP_NUM = 6
  INTEGER(LONG), DIMENSION(PHASE_LOOP_NUM) :: Phase_Type
  INTEGER(LONG) :: Loop_Index
  INTEGER(LONG) :: Phase_Count

  ! How many clear phase pixels in target scene?

  Phase_Type(1) = COUNT(Box_Data%Cloud_Phase == 0)

  ! How many water phase pixels in target scene?

  Phase_Type(2) = COUNT(Box_Data%Cloud_Phase == 1)

  ! How many supercooled phase pixels in target scene?

  Phase_Type(3) = COUNT(Box_Data%Cloud_Phase == 2)

  ! How many mixed phase pixels in target scene?

  Phase_Type(4) = COUNT(Box_Data%Cloud_Phase == 3)

  ! How many ice phase pixels in target scene?

  Phase_Type(5) = COUNT(Box_Data%Cloud_Phase == 4)

  ! How many unknown phase pixels in target scene?

  Phase_Type(6) = COUNT(Box_Data%Cloud_Phase == 5)

  ! Find main phase type

  Phase_Count = 0
  Main_Cloud_Phase = 0

  ! Skip clear phase, so begin loop at Index 2.

  phase_loop: DO Loop_Index = 2, PHASE_LOOP_NUM

    IF (Phase_Type(Loop_Index) > Phase_Count) THEN

      Phase_Count = Phase_Type(Loop_Index)
      Main_Cloud_Phase = Loop_Index - 1

    ENDIF

  END DO phase_loop

END SUBROUTINE Dominant_Cloud_Phase

!-------------------------------------------------------------------------------
!
! Name:
!   Dominant_Cloud_Type
!
! Function:
!   Find dominant cloud type in target scene. Clear type is ignored.
!
! Description:
!   Returns the dominant cloud type in the target scene.  Will become
!   a diagnostic output variable.  Technically, FOG_TYPE and OVERSHOOTING_TYPE
!   are not baseline cloud type's, but for ease of coding, they are included.
!     CLEAR_TYPE        = 0     clear cloud mask
!     FOG_TYPE          = 1
!     WATER_TYPE        = 2     liquid water clouds, warmer than 273.16 K
!     SUPERCOOLED_TYPE  = 3     liquid water clouds, colder than 273.16 K
!     MIXED_TYPE        = 4     very cold liquid water clouds with some ice
!     TICE_TYPE         = 5     thick, opaque ice clouds
!     CIRRUS_TYPE       = 6     thin, semi-transparent ice clouds
!     OVERLAP_TYPE      = 7     multi-layer clouds, ice over water
!     OVERSHOOTING_TYPE = 8
!     UNKNOWN_TYPE      = 9     bad input data
!
! Reference:
!   None.
!
! Calling Sequence:
!   CALL Dominant_Cloud_Type
!
! Inputs:
!   Box_Data                 Target box data structure
!
! Input/Output:
!   Output                   Output data structure of current target scene
!
! Outputs:
!   Main_Cloud_Phase         Dominant cloud type of target scene.
!
! Dependencies:
!   None.
!
! Restrictions:
!   None.
!
! History:
! 08/2010 - Steve Wanzong - Created
! ------------------------------------------------------------------------------
SUBROUTINE Dominant_Cloud_Type(Box_Data, Main_Cloud_Type)

  TYPE(Target_Box_Data), INTENT(IN) :: Box_Data

  INTEGER(LONG), INTENT(OUT) :: Main_Cloud_Type

  INTEGER(LONG), PARAMETER :: TYPE_LOOP_NUM = 10
  INTEGER(LONG), DIMENSION(TYPE_LOOP_NUM) :: Cloud_Type
  INTEGER(LONG) :: Loop_Index
  INTEGER(LONG) :: Type_Count


  ! How many clear type pixels in target scene?

  Cloud_Type(1) = COUNT(Box_Data%Cloud_Type == 0)

  ! How many fog type pixels in target scene?

  Cloud_Type(2) = COUNT(Box_Data%Cloud_Type == 1)

  ! How many water type pixels in target scene?

  Cloud_Type(3) = COUNT(Box_Data%Cloud_Type == 2)

  ! How many supercooled type pixels in target scene?

  Cloud_Type(4) = COUNT(Box_Data%Cloud_Type == 3)

  ! How many mixed type pixels in target scene?

  Cloud_Type(5) = COUNT(Box_Data%Cloud_Type == 4)

  ! How many tice type pixels in target scene?

  Cloud_Type(6) = COUNT(Box_Data%Cloud_Type == 5)

  ! How many cirrus type pixels in target scene?

  Cloud_Type(7) = COUNT(Box_Data%Cloud_Type == 6)

  ! How many overlap type pixels in target scene?

  Cloud_Type(8) = COUNT(Box_Data%Cloud_Type == 7)

  ! How many overshooting type pixels in target scene?

  Cloud_Type(9) = COUNT(Box_Data%Cloud_Type == 8)

  ! How many unknown type pixels in target scene?

  Cloud_Type(10) = COUNT(Box_Data%Cloud_Type == 9)

  ! Find main cloud type

  Type_Count = 0
  Main_Cloud_Type = 0

  ! Skip clear type, so begin loop at Index 2.

  type_loop: DO Loop_Index = 2, TYPE_LOOP_NUM

    IF (Cloud_Type(Loop_Index) > Type_Count) THEN

      Type_Count = Cloud_Type(Loop_Index)
      Main_Cloud_Type = Loop_Index - 1

    ENDIF

  END DO type_loop

END SUBROUTINE Dominant_Cloud_Type







END MODULE AMV_EN_TARGET_SELECTION_UTILS_M
