MODULE AMV_EN_FEATURE_TRACKING_UTILS_M
!-------------------------------------------------------------------------------
!> @file
! Name:
!   Feature_Tracking_Utils
!
! Function:
!   A collection of pattern-matching utilities used to derive atmospheric motion
!   vectors from an image sequence (3 or more images).
!
! Description:
!   The routines in this module are called by the main target selection program
!   to perform correlation matching of satellite features.
!
! References:
!   Velden et al (2005), "Recent Innovations in Deriving Tropospheric Winds From
!   Meteorological Satellites." Bulletin of the American Meteorological Society,
!   February 2005.
!
!   Nieman et al (1997), "Fully Automated Cloud-Drift Winds in NESDIS
!   Operations." Bulletin of the American Meteorological Society, June 1997.
!
! Calling sequence:
!   See description below
!
! Inputs:
!   See description below
!
! Outputs:
!   See description below
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!     07/2007 - Wayne Bresky - Created
!  04/27/2009 - Wayne Bresky - Nested tracking introduced
!     11/2009 - Wayne Bresky - DBSCAN cluster analysis implemented within
!                              nested tracking
!     04/2010 - Wayne Bresky - Code Unit Test Review (CUTR) changes
!> @ingroup WINDS_AMV_EN
!-------------------------------------------------------------------------------

 USE TYPE_KINDS_AIT
 USE GEOCAT_CONSTANTS
 USE Framework_Global_Variables_Module
 USE NF_PARM
 USE COMMON_VAR_VALUE

 USE AMV_EN_TargetBox_m
 USE AMV_EN_DBSCAN_m
 USE AMV_EN_Output_m
 USE AMV_EN_NAV_MSG_M
 USE AMV_EN_NAV_ABI_M
 USE AMV_EN_NAV_POLAR_M
 USE AMV_EN_Nav_M
 USE AMV_EN_WINDS_INC




IMPLICIT NONE

! sub-vector number (1=reverse, 2=forward)
INTEGER(LONG), PARAMETER :: REVERSE = 1
INTEGER(LONG), PARAMETER :: FORWARD = 2

INTEGER(SHORT), PARAMETER :: MATCH_ON_BOUNDARY_FAILURE  = 15
INTEGER(SHORT), PARAMETER :: MISSING_DATA_SEARCH        = 20
INTEGER(SHORT), PARAMETER :: ZERO_INITIAL_SAMPLE        = 21
INTEGER(SHORT), PARAMETER :: ZERO_FINAL_SAMPLE          = 22

! Target box dimensions
INTEGER(LONG) :: Target_Line_Size
INTEGER(LONG) :: Target_Element_Size

! Search box dimensions
INTEGER(LONG) :: Search_Line_Size
INTEGER(LONG) :: Search_Element_Size

! Correlation lag matrix dimensions
INTEGER(LONG) :: Lag_Line_Size
INTEGER(LONG) :: Lag_Element_Size
REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Lag_Matrix
REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: XCorr_Values

! 1-D data
REAL(SINGLE), DIMENSION(:), ALLOCATABLE :: Target_1D
REAL(SINGLE), DIMENSION(:), ALLOCATABLE :: Search_1D

! DBSCAN parameters
INTEGER(BYTE),PARAMETER :: APPEND_PT = 0
INTEGER(BYTE),PARAMETER :: DELETE_PT = 1
REAL(SINGLE), PARAMETER :: UNCLASSIFIED = 0.0
REAL(SINGLE), PARAMETER :: NOISE = 1.0
REAL(SINGLE), PARAMETER :: CLASSIFIED = 2.0


CONTAINS

!-------------------------------------------------------------------------------
!
! Name:
!   Compute_Correlation
!
! Function:
!   Finds a match in the brightness temperature field of a subsequent image.
!
! Description:
!   A pattern matching routine which utilizes one of two methods to find
!   a match in the brightness temperature field of a larger search region.
!   The match is found by either minimizing the sum-of-squared differences
!   (SSD) or maximizing the linear correlation (also called the cross
!   correlion). The routine returns a displacement of the match "relative"
!   to the center of the search region. The program that calls this routine
!   must convert the relative displacement to an actual displacement.
!
! References:
!   Nieman et al (1997), "Fully Automated Cloud-Drift Winds in NESDIS
!   Operations." Bulletin of the American Meteorological Society, June 1997.
!
!   Merrill (1989), "Advances in the automated production of wind estimates
!   from geostationary satellite imaging." Proc. Fourth Conf. Satellite
!   Meteorology, San Diego, CA, Amer. Meteor. Soc., 246-249.
!
! Calling sequence:
!   Call Compute_Correlation
!
! Input:
!   Box_Data                  Target box data structure
!   Search_BT_Values          Brightness temperature values for search region
!   Lag_Table                 List of positions in spiral
!
! Output:
!   Line_Disp                 Line displacement of match relative to center
!                             of search
!   Element_Disp              Element displacement of match  relative to
!                             center of search
!   Least_Square_Error        Minimum sum of squared differences value
!   Target_Scene_Mean         Mean brightness temperature value of target
!   Match_Scene_Mean          Mean brightness temperature value of match
!   Target_Scene_StdDev       Standard deviation of target scene brightness
!                             temperature values
!   Match_Scene_StdDev        Standard deviation of match scene brightness
!                             temperature values
!   Corr_Coeff                Linear correlation value of matching scene
!   Status_Flag               Success value:
!                                 0 - Match found
!                               200 - Match found on boundary of lag matrix
!                               300 - Missing data
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!  08/2007 - Wayne Bresky - Created
!  04/2010 - Wayne Bresky - CUTR changes
!
!-------------------------------------------------------------------------------
SUBROUTINE Compute_Correlation(Box_Data, Search_BT_Values, Lag_Table, &
                                          Delta_X_Fcst, Delta_Y_Fcst, &
                         Line_Disp, Element_Disp, Least_Square_Error, &
            Target_Scene_Mean, Match_Scene_Mean, Target_Scene_StdDev, &
                         Match_Scene_StdDev, Corr_Coeff, Status_Flag, &
                                                          SSD_Matrix)

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

  ! Lag matrix center point
  INTEGER(LONG) :: Lag_Center_Element
  INTEGER(LONG) :: Lag_Center_Line
  INTEGER(LONG) :: Total_Points_Lag

  ! Flag to indicate success or failure of correlation matching
  INTEGER(LONG), INTENT(OUT) :: Status_Flag

  ! Least square error and linear correlation coefficient associated match
  INTEGER(LONG), INTENT(OUT) :: Least_Square_Error
  REAL(SINGLE), INTENT(OUT) :: Corr_Coeff

  INTEGER(LONG), DIMENSION(:, :), INTENT(IN) :: Lag_Table

  REAL(SINGLE), INTENT(INOUT), DIMENSION(:) :: SSD_Matrix

  ! Minimum SSD value and location
  REAL(DOUBLE) :: Min_SSD_Value
  INTEGER(LONG) :: Min_SSD_Location
  INTEGER(LONG) :: Min_SSD_Element
  INTEGER(LONG) :: Min_SSD_Line

  INTEGER(LONG) :: Match_Lag_Element
  INTEGER(LONG) :: Match_Lag_Line

  ! Offset for location of match area within search grid
  INTEGER(LONG) :: Match_Offset
  INTEGER(LONG) Diff_In_Size

  INTEGER(LONG) :: Element_Counter_Search
  INTEGER(LONG) :: Element
  INTEGER(LONG) :: Line
  INTEGER(LONG) :: Search_Index
  INTEGER(LONG) :: Total_Points_Search

  ! Brightness Temperature data and statistics for target and match scenes
  REAL(SINGLE), INTENT(IN), DIMENSION(:,:) :: Search_BT_Values
  REAL(SINGLE), INTENT(IN) :: Delta_X_Fcst
  REAL(SINGLE), INTENT(IN) :: Delta_Y_Fcst

  REAL(DOUBLE), INTENT(OUT) :: Target_Scene_Mean
  REAL(DOUBLE), INTENT(OUT) :: Target_Scene_StdDev
  REAL(DOUBLE), INTENT(OUT) :: Match_Scene_Mean
  REAL(DOUBLE), INTENT(OUT) :: Match_Scene_StdDev

  REAL(DOUBLE) :: Sum_Target
  REAL(DOUBLE) :: Sum_Squared_Target
  REAL(DOUBLE) :: Sum_Squared_Match
  REAL(DOUBLE) :: Target_Scene_Variance
  REAL(DOUBLE) :: Match_Scene_Variance

  ! Vector displcement in satellite coordinates
  REAL(SINGLE), INTENT(OUT) :: Line_Disp
  REAL(SINGLE), INTENT(OUT) :: Element_Disp
  REAL(SINGLE) :: Fractional_Element
  REAL(SINGLE) :: Fractional_Line
  INTEGER(LONG) :: Total_Points_Target
  REAL(DOUBLE) :: Corr_Coeff_DB

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

  Target_Element_Size = SIZE(Box_Data%BrtTemp, 1)
  Target_Line_Size = SIZE(Box_Data%BrtTemp, 2)
  Total_Points_Target = Target_Element_Size * Target_Line_Size
  Search_Element_Size = SIZE(Search_BT_Values, 1)
  Search_Line_Size = SIZE(Search_BT_Values, 2)
  Total_Points_Search = Search_Element_Size * Search_Line_Size
  Line_Disp = 0.0
  Element_Disp = 0.0
  Least_Square_Error = 0
  Status_Flag = 0
  Corr_Coeff = 0.0
  Target_Scene_Mean = 0.0
  Target_Scene_StdDev = 0.0
  Match_Scene_Mean = 0.0
  Match_Scene_StdDev = 0.0
  Sum_Target = 0.0
  Sum_Squared_Target = 0.0
  Sum_Squared_Match = 0.0
  Target_Scene_Variance = 0.0
  Match_Scene_Variance = 0.0
  Diff_In_Size = 0



  ALLOCATE(Target_1D(Total_Points_Target), Search_1D(Total_Points_Search))

  !-----------------------------------------------------------------------------
  ! Reshape 2-D box data into vectors for pattern matching code.
  !-----------------------------------------------------------------------------

  Target_1D = RESHAPE(Box_Data%BrtTemp, (/Total_Points_Target/))
  Search_1D = RESHAPE(Search_BT_Values, (/Total_Points_Search/))

  ! Compute summations for target scene

  summation_loop: DO Element = 1, Total_Points_Target

    Sum_Target = Sum_Target + Target_1D(Element)
    Sum_Squared_Target = Sum_Squared_Target + Target_1D(Element)**2

  END DO summation_loop

  ! ----------------------------------------------------------------------------
  ! Specify "lag" coefficient matrix size. This array contains the SSD
  ! values for all possible lag positions of the template within the
  ! search region. This array is filled regardless of the tracking
  ! method because it is used by both the LP norm and CC methods.
  !-----------------------------------------------------------------------------

  Lag_Element_Size = Search_Element_Size - Target_Element_Size + 1
  Lag_Line_Size = Search_Line_Size - Target_Line_Size + 1
  Lag_Center_Element = (Lag_Element_Size + 1) / 2
  Lag_Center_Line = (Lag_Line_Size + 1) / 2
  Total_Points_Lag = Lag_Element_Size * Lag_Line_Size

  ALLOCATE(Lag_Matrix(Total_Points_Lag))

  ! Call pattern match code

  CALL Fill_Lag_Matrix_1D(Lag_Element_Size, Diff_In_Size, Target_1D, Search_1D, &
                                                         Lag_Table, Lag_Matrix)

  SSD_Matrix = Lag_Matrix

  ! ----------------------------------------------------------------------------
  ! Determine location of match within lag matrix.
  ! ----------------------------------------------------------------------------

  Min_SSD_Value = MINVAL(Lag_Matrix)
  Least_Square_Error = Min_SSD_Value
  Min_SSD_Location = MINLOC(Lag_Matrix, 1)
  Min_SSD_Element = MOD(Min_SSD_Location, Lag_Element_Size)
  IF (MOD(Min_SSD_Location, Lag_Element_Size) .EQ. 0) THEN
    Min_SSD_Element = Lag_Element_Size
    Min_SSD_Line = (Min_SSD_Location / Lag_Element_Size)
  ELSE
    Min_SSD_Line = (Min_SSD_Location / Lag_Element_Size) + 1
  ENDIF

  ! ----------------------------------------------------------------------------
  ! Flag match on boundary. This condition may indicate the true match
  ! is located beyond the domain of the search.
  ! ----------------------------------------------------------------------------

  IF (Min_SSD_Element .EQ. 1 .OR. Min_SSD_Element .EQ. Lag_Element_Size .OR. &
              Min_SSD_Line .EQ. 1 .OR. Min_SSD_Line .EQ. Lag_Line_Size) THEN

    !Status_Flag = 200
    Status_Flag = MATCH_ON_BOUNDARY_FAILURE
    DEALLOCATE(Target_1D, Search_1D, Lag_Matrix)
    RETURN

  ENDIF

  ! ----------------------------------------------------------------------------
  ! Compute displacement relative to center of lag array. The center of
  ! the lag array also corresponds to the center of the search box.
  ! ----------------------------------------------------------------------------

  Element_Disp = REAL(Min_SSD_Element - Lag_Center_Element)
  Line_Disp = REAL(Min_SSD_Line - Lag_Center_Line)

  !-----------------------------------------------------------------------------
  ! Offset for location of match area within search grid.
  !-----------------------------------------------------------------------------

  Match_Offset = (Min_SSD_Line - 1) * Search_Element_Size + Min_SSD_Element - 1
  Element_Counter_Search = - Search_Element_Size

  line_loop: DO Line = 1, Target_Line_Size

    Element_Counter_Search = Element_Counter_Search + Search_Element_Size

    element_loop: DO Element = 1, Target_Element_Size

      Search_Index = Element_Counter_Search + Element + Match_Offset
      Match_Scene_Mean = Match_Scene_Mean + Search_1D(Search_Index)
      Sum_Squared_Match = Sum_Squared_Match + Search_1D(Search_Index)**2

    END DO element_loop

  END DO line_loop

  ! ----------------------------------------------------------------------------
  ! Compute scene statistics.
  ! ----------------------------------------------------------------------------

  Target_Scene_Mean = Sum_Target / REAL(Total_Points_Target)
  Match_Scene_Mean = Match_Scene_Mean / REAL(Total_Points_Target)

  Target_Scene_Variance = Sum_Squared_Target / REAL(Total_Points_Target) - &
                                                      Target_Scene_Mean**2

  Match_Scene_Variance = Sum_Squared_Match / REAL(Total_Points_Target) - &
                                                     Match_Scene_Mean**2

  Target_Scene_StdDev = SQRT(Target_Scene_Variance)
  Match_Scene_StdDev = SQRT(Match_Scene_Variance)

  ! ----------------------------------------------------------------------------
  ! Compute linear correlation coefficient.
  ! ----------------------------------------------------------------------------
  Corr_Coeff_DB = (Sum_Squared_Match + Sum_Squared_Target - Min_SSD_Value) / &
       (2.0 * Total_Points_Target) - Target_Scene_Mean * Match_Scene_Mean

  Corr_Coeff = Corr_Coeff_DB / (Target_Scene_StdDev * Match_Scene_StdDev)


  ! ----------------------------------------------------------------------------
  ! Compute sub-pixel displacement. Interpolate lag matrix values to
  ! get fractional line and element displacements.
  ! ----------------------------------------------------------------------------


  Match_Lag_Element = Min_SSD_Element
  Match_Lag_Line = Min_SSD_Line

  CALL Interpolate_Disp(Match_Lag_Element, Match_Lag_Line, Fractional_Element, &
                                                              Fractional_Line)

  ! Add fractional values to integer displacement.

  Element_Disp = Element_Disp + Fractional_Element
  Line_Disp = Line_Disp + Fractional_Line

  ! convert to absolute displacement

  Element_Disp = Element_Disp + Delta_X_Fcst
  Line_Disp = Line_Disp + Delta_Y_Fcst

  DEALLOCATE(Target_1D, Search_1D, Lag_Matrix)

  IF (Corr_Coeff .LT. 0.) Corr_Coeff = 0.0

END SUBROUTINE Compute_Correlation

!-------------------------------------------------------------------------------
!
! Name:
!   Compute_UV
!
! Function:
!   Computes the u- and v- components of the wind from the beginning and ending
!   earth locations of the cloud feature. Computation is done in spherical
!   coordinates.
!
! Description:
!
! References:
!
! Calling sequence:
!   CALL Compute_UV
!
! Inputs:
!   Lat1                 Initial latitude of cloud feature
!   Lon1                 Initial longitude of cloud feature
!   Lat2                 Final latitude of cloud feature
!   Lon2                 Final longitude of cloud feature
!   Delta_Time           Time interval in fractional hours
!
! Outputs:
!   U                    U-component of wind (positive is west, negative is
!                        east) in meters per second
!   V                    V-component of wind (positive is south, negative is
!                        north) in meters per second
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!   7/2007 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE  Compute_UV(NavType, Lat1, Lon1, Lat2, Lon2, Delta_Time, U, V)

  REAL(SINGLE), INTENT(IN) :: Lat1
  REAL(SINGLE), INTENT(IN) :: Lon1
  REAL(SINGLE), INTENT(IN) :: Lat2
  REAL(SINGLE), INTENT(IN) :: Lon2
  REAL(SINGLE), INTENT(IN) :: Delta_Time
  Real(SINGLE) :: Delta_Lat
  Real(SINGLE) :: Delta_Lon
  Real(SINGLE) :: Cosine_Of_Lat
  Real(SINGLE) :: Avg_Lat
  REAL(SINGLE) :: Absolute_Delta_Time
  REAL(SINGLE) :: Delta_Time_Sec
  REAL(SINGLE), INTENT(OUT) :: U
  REAL(SINGLE), INTENT(OUT) :: V

  !  Earth equatorial radius (km)
  REAL(DOUBLE), PARAMETER :: R_EQ = 6378.137D0

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

  !
  ! Initialize wind components
  !

  U = 0.0
  V = 0.0

  !
  ! Convert hours to seconds
  !

  Delta_Time_Sec = Delta_Time * SEC_IN_A_HOUR
  Absolute_Delta_Time = ABS(Delta_Time) * SEC_IN_A_HOUR

  !
  ! Compute change in latitude and longitude
  !

  Delta_Lat = Lat2 - Lat1
  Delta_Lon = Lon2 - Lon1

  !
  ! Compute average latitude and cosine of average latitude
  !

  Avg_Lat = (Lat1 + Lat2) / 2.0
  Cosine_Of_Lat = COS(DTOR * Avg_Lat)

  !
  ! Convert longitude to range -180/+180 if necessary
  !

  IF (ABS(Delta_Lon) .GT. DEGREE180 .AND. Delta_Lon .GT. 0.0) Delta_Lon = &
                                                    Delta_Lon - DEGREE360
  IF (ABS(Delta_Lon) .GT. DEGREE180 .AND. Delta_Lon .LT. 0.0) Delta_Lon = &
                                                    Delta_Lon + DEGREE360
  !
  ! Convert change in degrees to radians
  !

  Delta_Lon = Delta_Lon * DTOR
  Delta_Lat = Delta_Lat * DTOR

  !
  ! Compute wind components
  !

  !Multiplying by 1000 converts the units from kilometers to meters


  U = R_EQ * KM2M * Cosine_Of_Lat * Delta_Lon / Delta_Time_Sec
  V = R_EQ * KM2M * Delta_Lat / Delta_Time_Sec

  IF (TRIM(NavType) .EQ. 'PS' ) THEN
      U = R_EQ * KM2M * Cosine_Of_Lat * Delta_Lon / Absolute_Delta_Time
      V = R_EQ * KM2M * Delta_Lat / Absolute_Delta_Time
  ENDIF

END SUBROUTINE Compute_UV

!-------------------------------------------------------------------------------
!
! Name:
!   Get_Fcst_Displacement
!
! Function:
!   Compute a (line and element) displacement from a starting position using the
!   input speed and time interval.
!
! Description:
!   This routine computes the expected displacement of a feature given a
!   starting location, a time interval and the forecast u- and v- wind
!   components. Note the time can be negative or positive. The displacement
!   returned by this routine is used to compute the ceter location of the
!   search region.


! Calling sequence:
!   CALL Get_Fcst_Displacement
!
! Inputs:
!   Starting_Lat         Initial latitude of cloud feature
!   Starting_Lon         Initial longitude of cloud feature
!   Delta_Time           Time interval in fractional hours (+ or -)
!   U                    Forecast u-component
!   V                    Forecast v-component
!   Sat_ID               Satellite ID
!   Target_Element       Center element of target box
!   Target_Line          Center line of target box
!
! Outputs:
!   Search_Element       Center element of search box
!   Search_Line          Center line of search box
!   Delta_X_Fcst         Forecast element displacement
!   Delta_Y_Fcst         Forecast line displacement
!
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!  04/2010 - Wayne Bresky - Created from existing "New Position" subroutine
!
!-------------------------------------------------------------------------------
SUBROUTINE Get_Fcst_Displacement(Nav, NavType, Starting_Lat, Starting_Lon, Delta_Time, U, V, &
                                 Sensor_Series, Target_Element, Target_Line, Delta_X_Fcst, Delta_Y_Fcst)
  TYPE(Nav_t), INTENT(IN):: Nav
  CHARACTER(LEN=*), INTENT(IN) :: NavType
  INTEGER(LONG) :: Col
  INTEGER(LONG) :: Row
  INTEGER(LONG), INTENT(IN) :: Sensor_Series

  REAL(SINGLE), INTENT(IN) :: Starting_Lat
  REAL(SINGLE), INTENT(IN) :: Starting_Lon
  REAL(SINGLE), INTENT(IN) :: Delta_Time
  REAL(SINGLE), INTENT(IN) :: U
  REAL(SINGLE), INTENT(IN) :: V
  REAL(SINGLE), INTENT(IN) :: Target_Element
  REAL(SINGLE), INTENT(IN) :: Target_Line
  REAL(SINGLE) :: Ending_Lat
  REAL(SINGLE) :: Ending_Lon
  REAL(SINGLE), INTENT(OUT) :: Delta_X_Fcst
  REAL(SINGLE), INTENT(OUT) :: Delta_Y_Fcst
  REAL(SINGLE) :: Cosine_Of_Lat
  REAL(SINGLE) :: Avg_Lat
  REAL(SINGLE) :: Col_Real
  REAL(SINGLE) :: Row_Real

  !  Earth equatorial radius (km)
  REAL(DOUBLE), PARAMETER :: R_EQ = 6378.137D0



  ! Multiplying by 60**2 converts time from hours to seconds
  ! Multiplying by 1000 converts the units from kilometers to meters
  Ending_Lat = Starting_Lat + (V * Delta_Time * SEC_IN_A_HOUR) &
                        / (R_EQ * KM2M * DTOR)

  Avg_Lat = (Starting_Lat + Ending_Lat) / 2.0
  Cosine_Of_Lat = COS(DTOR * Avg_Lat)

  SELECT CASE (Sensor_Series)

  CASE (SERIES_MSG_SEVIRI)

    ! positive u corresponds with positive change in longitude

    ! Multiplying by 60**2 converts time from hours to seconds
    ! Multiplying by 1000 converts the units from kilometers to meters
    Ending_Lon = Starting_Lon + (U * Delta_Time * SEC_IN_A_HOUR) / &
                (R_EQ * KM2M * Cosine_Of_Lat * DTOR)

    ! Convert earth location to image pixel

    CALL geocoord2pixcoord_dp(Nav, DBLE(Ending_Lat), DBLE(Ending_Lon), Col_Real, Row_Real)

    Col = NINT(Col_Real)
    Row = NINT(Row_Real)

!  CASE (SERIES_GOES_N_IMG)
!
!    ! positive u corresponds with positive change in longitude
!
!    ! Multiplying by 60**2 converts time from hours to seconds
!    ! Multiplying by 1000 converts the units from kilometers to meters
!    Ending_Lon = Starting_Lon + (U * Delta_Time * SEC_IN_A_HOUR) / &
!                (R_EQ * KM2M * Cosine_Of_Lat * DTOR)
!
!    ! Convert earth location to image pixel
!
!    CALL geocoord2pixcoord_gvar(Nav, Ending_Lat, Ending_Lon, 1, Col_Real, Row_Real)
!
!    Col = NINT(Col_Real)
!    Row = NINT(Row_Real)

  CASE (SERIES_GOESR_ABI,SERIES_HIMAWARI_AHI)

    ! Multiplying by 60**2 converts time from hours to seconds
    ! Multiplying by 1000 converts the units from kilometers to meters
    Ending_Lon = Starting_Lon + (U * Delta_Time * SEC_IN_A_HOUR) / &
                (R_EQ * KM2M * Cosine_Of_Lat * DTOR)

    ! Convert earth location to image pixel

    CALL geocoord2pixcoord_abi(Nav,  DBLE(Ending_Lat), DBLE(Ending_Lon), Col, Row)


  END SELECT

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

    ! positive u corresponds with positive change in longitude

    ! Multiplying by 60**2 converts time from hours to seconds
    ! Multiplying by 1000 converts the units from kilometers to meters
    Ending_Lon = Starting_Lon - (U * Delta_Time * SEC_IN_A_HOUR) / &
                (R_EQ * KM2M * Cosine_Of_Lat * DTOR)

    CALL geocoord2pixcoord_ps(Nav, Ending_Lat, Ending_Lon, Col_Real, Row_Real)

    Col = NINT(Col_Real)
    Row = NINT(Row_Real)

  ENDIF

  Delta_X_Fcst = REAL(Col) - Target_Element
  Delta_Y_Fcst = REAL(Row) - Target_Line

  ! EUMETSAT coordinate system is reversed from GEOCAT coordinate system,
  ! so we need to reverse the sign of the forecast displacement

  IF (Sensor_Series == SERIES_MSG_SEVIRI) THEN
    Delta_X_Fcst = -Delta_X_Fcst
    Delta_Y_Fcst = -Delta_Y_Fcst
  ENDIF

END SUBROUTINE Get_Fcst_Displacement

!-------------------------------------------------------------------------------
!
! Name:
!   Interpolate_Disp
!
! Function:
!   Compute sub-pixel displacement through interpolation.
!
! Description:
!   This routine computes the fractional displacement through interpolation
!   of the 4 points surrounding the match location.

! Calling sequence:
!   CALL Interpolate_Disp
!
! Inputs:
!   Match_Lag_Element    Element coordinate of match within lag matrix
!   Match_Lag_Line       Line coordinate of match within lag matrix
!
! Outputs:
!   Fractional_Element   Fractional element displacement
!   Fractional_Line      Fractional line displacement
!
! Dependencies:
!   Requires the lag matrix to be filled with values.
!
! Restrictions:
!   None
!
! History:
!   7/2007 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE Interpolate_Disp(Match_Lag_Element, Match_Lag_Line, &
                          Fractional_Element, Fractional_Line)

  INTEGER(LONG), INTENT(IN) :: Match_Lag_Element
  INTEGER(LONG), INTENT(IN) :: Match_Lag_Line
  REAL(SINGLE), INTENT(OUT) :: Fractional_Element
  REAL(SINGLE), INTENT(OUT) :: Fractional_Line
  INTEGER(LONG) :: First_Index
  INTEGER(LONG) :: Middle_Index
  INTEGER(LONG) :: Third_Index
  REAL(SINGLE) :: Denominator
  REAL(SINGLE) :: Middle_Point
  REAL(SINGLE) :: First_Point
  REAL(SINGLE) :: Third_Point

  Fractional_Element = 0.0
  Fractional_Line = 0.0

  ! Check if match is on boundary. If so, do not interpolate.

  IF (Match_Lag_Element .EQ. 1 .OR. Match_Lag_Element .EQ. Lag_Element_Size &
  .OR. Match_Lag_Line .EQ. 1 .OR. Match_Lag_Line .EQ. Lag_Line_Size) RETURN

  ! Line interpolation

  Middle_Index = Lag_Element_Size * (Match_Lag_Line - 1) + Match_Lag_Element
  First_Index = Middle_Index - Lag_Element_Size
  Third_Index = Middle_Index + Lag_Element_Size

  First_Point = Lag_Matrix(First_Index)
  Middle_Point = Lag_Matrix(Middle_Index)
  Third_Point = Lag_Matrix(Third_Index)

  ! Linear interpolation equation in 1D
  Denominator = 2.0 * (First_Point + Third_Point - 2.0 * Middle_Point)

  IF (Denominator .NE. 0.0) Fractional_Line = (First_Point - Third_Point) / &
                                                                Denominator

  ! Element interpolation

  First_Index = Middle_Index - 1
  Third_Index = Middle_Index +  1

  First_Point = Lag_Matrix(First_Index)
  Third_Point = Lag_Matrix(Third_Index)

  ! Linear interpolation equation in 1D
  Denominator = 2.0 * (First_Point + Third_Point - 2.0 * Middle_Point)

  IF (Denominator .NE. 0.0) Fractional_Element = (First_Point - Third_Point) / &
                                                                   Denominator

END SUBROUTINE Interpolate_Disp

!-------------------------------------------------------------------------------
!
! Name:
!   Generate_Spiral_Table
!
! Function:
!   Produces a table of positions that spirals out from from the search
!   center point.
!
! Description:
!   A simple routine for generating a prioritized list of positions to search.
!   The list of positions describes a clockwise spiral starting at the center
!   point of the search region. By using the spiral table, the search for
!   the match begins at the location predicted by the forecast instead of
!   beginning at the top left corner as is the usual case. This list,
!   when used with the strategy of terminating the SSD summation early
!   once a minimum has been exceeded, speeds up the computational performance
!   of the search process by setting the SSD minimum to a low initial value.
!
! Calling sequence:
!   Generate_Spiral_Table
!
! Inputs:
!   Search_Size          Size of Search Box
!   Lag_Size             Size of lag matrix
!   Size_Diff            Size difference in pixels between the nested box and
!                        the outer target box.
!
! Outputs:
!   Lag_Table2           List of positions in spiral
!
! Dependencies:
!   Assumes SSD similarity measure is being used for tracking.
!
! Restrictions:
!   None
!
! History:
!   06/01/2009 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE Generate_Spiral_Table(Search_Size, Lag_Size, Size_Diff, Lag_Table2)

  INTEGER(BYTE) Sign_Pixel
  INTEGER(LONG) I
  INTEGER(LONG) J
  INTEGER(LONG) Point
  INTEGER(LONG) Element
  INTEGER(LONG) Line
  INTEGER(LONG) Lag_Index
  INTEGER(LONG) Search_Index
  INTEGER(LONG) Small_Search_Element_Size
  INTEGER(LONG) Lag_Center_Element
  INTEGER(LONG) Lag_Center_Line
  INTEGER(LONG) Total_Points_Lag
  INTEGER(LONG), INTENT(IN) :: Search_Size
  INTEGER(LONG), INTENT(IN) :: Lag_Size
  INTEGER(LONG), INTENT(IN) :: Size_Diff
  INTEGER(LONG), INTENT(INOUT), DIMENSION(:,:) :: Lag_Table2

  Small_Search_Element_Size = Search_Size - Size_Diff
  Lag_Center_Element = (Lag_Size + 1) / 2
  Lag_Center_Line = Lag_Center_Element
  Total_Points_Lag = Lag_Size**2

  Point = 1
  Sign_Pixel = -1
  Element = Lag_Center_Element
  Line = Lag_Center_Line
  Lag_Table2(1,1) = Lag_Center_Element
  Lag_Table2(1,2) = Lag_Center_Line
  Lag_Table2(1,3) = (Total_Points_Lag) / 2 + 1
  Lag_Table2(1,4) = (Lag_Center_Line - 1) * Small_Search_Element_Size &
                                                + Lag_Center_Element

  DO I = 1, Total_Points_Lag

    Sign_Pixel = Sign_Pixel * (-1)

    DO J = 1,I

      Point = Point + 1

      IF (Point .GT. Total_Points_Lag) RETURN

      Element = Element + Sign_Pixel

      IF (Element .GT. 0 .AND. Line .GT. 0 .AND. Element .LE. &
         Lag_Size .AND. Line .LE. Lag_Size) THEN

        Lag_Index = (Line - 1) * Lag_Size + Element
        Search_Index = (Line - 1) * Small_Search_Element_Size + Element
        Lag_Table2(Point,1) = Element
        Lag_Table2(Point,2) = Line
        Lag_Table2(Point,3) = Lag_Index
        Lag_Table2(Point,4) = Search_Index

      ELSE

        Point = Point - 1

      ENDIF

    ENDDO

    DO J = 1,I

      Point = Point + 1

      IF (Point .GT. Total_Points_Lag) RETURN

      Line = Line + Sign_Pixel

      IF (Element .GT. 0 .AND. Line .GT. 0 .AND. Element .LE. &
         Lag_Size .AND. Line .LE. Lag_Size) THEN

        Lag_Index = (Line - 1) * Lag_Size + Element
        Search_Index = (Line - 1) * Small_Search_Element_Size + Element
        Lag_Table2(Point,1) = Element
        Lag_Table2(Point,2) = Line
        Lag_Table2(Point,3) = Lag_Index
        Lag_Table2(Point,4) = Search_Index

      ELSE

        Point = Point - 1

      ENDIF

    ENDDO

  ENDDO

END SUBROUTINE Generate_Spiral_Table

!-------------------------------------------------------------------------------
!
! Name:
!   Nested_Tracking
!
! Function:
!   Version of Compute_Correlation designed to compute multiple displacement
!   estimates for small sub-regions (5x5 pixel box) of larger target box.
!   Essentially, the smaller box is "nested" within the larger box and shifted
!   one pixel at a time in order to compute the local motion for all interior
!   pixel locations.
!
! Description:
!   This is a modified version of Compute_Correlation that computes multiple
!   displacement estimates for small sub-regions of the larger target box.
!   A cluster analysis algorithm called DBSCAN is then used to find the
!   largest cluster in line/element space. The displacements associated with
!   the points in the largest cluster are then averaged to yield a final motion
!   estimate for the larger box. Likewise, the cloud top pressure values (or BT
!   when processing the WV) associated with the points in the largest cluster
!   are retained so that they can be used by target_selection to assign a median
!   pressure for the wind.
!
! References:
!   Daniels and Bresky (2010), "A New Nested Tracking Approach for Reducing the
!   Slow Speed Bias Associated with Atmospheric Motion Vectors (AMVs)."
!   Proceedings of the Tenth International Winds Workshop (IWW10).
!
! Calling sequence:
!   Call Nested_Tracking
!
! Input:
!   Vector               Vector number (1=reverse, 2=forward)
!   CSWV_Flag            Clear-sky WV flag (0=no, 1=yes)
!   Box_Data             Target box data structure
!   Search_BT_Values     Brightness temperature values for search region
!   Delta_X_Sat          Distance between target center and search center
!   Delta_Y_Sat          Distance between target center and search center
!   Band                 Satellite band
!   Lag_Table            List of positions in spiral search
!   Factor               pixel factor (ratio of E-W distance to N-S distance)
!   Output               output variables data structure of target
!
! Input/Output:
!   DBSCAN_Out           DBSCAN output structure
!   Status_Flag          Success value:
!                            0   - Match found
!                            300 - Missing data
!                            400 - Zero good sample
!
! Output:
!   Line_Disp            Line displacement of match relative to center of search
!   Element_Disp         Element displacement of match relative to center
!                        of search
!   Least_Square_Error   Minimum sum of squared differences value
!   Target_Scene_Mean    Mean brightness temperature value of target
!   Match_Scene_Mean     Mean brightness temperature value of match
!   Target_Scene_StdDev  Standard deviation of target scene brightness
!                        temperature values
!   Match_Scene_StdDev   Standard deviation of match scene brightness
!                        temperature values
!   Corr_Coeff           Linear correlation value of matching scene
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!   04/27/2009 - Wayne Bresky - Created
!   11/01/2009 - Wayne Bresky - DBSCAN cluster analysis added
!   08/10/2010 - Wayne Bresky - DBSCAN output structure created
!   07/27/2011 - Wayne Bresky - output variables data structure included in call
!   09/30/2025 - Wayne Bresky - modified to combine displacements from both time steps 
!                               and only call DBSCAN once 
!
!-------------------------------------------------------------------------------
SUBROUTINE Nested_Tracking(NavType, Vector, CSWV_Flag, Box_Data, Search_BT_Values, &
                                           Delta_X_Sat, Delta_Y_Sat, Band, Lag_Table2, &
                                          Line_Disp, Element_Disp, TandemTriplet, Least_Square_Error, & !<-- aab Tandem mod
                                          Factor, Target_Scene_Mean, Match_Scene_Mean, &
                                  Target_Scene_StdDev, Match_Scene_StdDev, Corr_Coeff, &
                                    Status_Flag, DBSCAN_Out, SetOfPoints_Init)
  CHARACTER(LEN=*), PARAMETER :: FUNC ="Nested_Tracking"
  CHARACTER(LEN=*), INTENT(IN) :: NavType
  ! Vector number (1=reverse, 2=forward)
  INTEGER(LONG), INTENT(IN) :: Vector
  ! clear sky flag
  INTEGER(LONG), INTENT(IN) :: CSWV_Flag
  TYPE(Target_Box_Data), INTENT(IN) :: Box_Data
  REAL(SINGLE), INTENT(IN), DIMENSION(:,:) :: Search_BT_Values !< Brightness Temperature for target and match scenes
  REAL(SINGLE), INTENT(IN) :: Delta_X_Sat !< Distance in pixels between target center and search center
  REAL(SINGLE), INTENT(IN) :: Delta_Y_Sat !< Distance in pixels between target center and search center
  INTEGER(LONG), INTENT(IN) :: Band !< Satellite band
  INTEGER(LONG), DIMENSION(:, :), INTENT(IN) :: Lag_Table2
  REAL(SINGLE), INTENT(OUT) :: Line_Disp !< Vector displcement in satellite coordinates
  REAL(SINGLE), INTENT(OUT) :: Element_Disp !< Vector displcement in satellite coordinates
  CHARACTER(LEN=*), INTENT(IN) :: TandemTriplet   !<-- aab Tandem mod
  INTEGER(LONG), INTENT(OUT) :: Least_Square_Error!< Least square error associated match
  REAL(SINGLE), INTENT(IN) :: Factor !< Normalization factor
  REAL(DOUBLE), INTENT(OUT) :: Target_Scene_Mean
  REAL(DOUBLE), INTENT(OUT) :: Match_Scene_Mean
  REAL(DOUBLE), INTENT(OUT) :: Target_Scene_StdDev
  REAL(DOUBLE), INTENT(OUT) :: Match_Scene_StdDev
  REAL(SINGLE), INTENT(OUT) :: Corr_Coeff !< linear correlation coefficient associated match
  INTEGER(LONG), INTENT(INOUT) :: Status_Flag !< Flag to indicate success or failure of correlation matching
  TYPE(DBSCAN_Output), INTENT(INOUT) :: DBSCAN_Out
  REAL(SINGLE), DIMENSION(:,:), INTENT(INOUT) :: SetOfPoints_Init


  INTEGER(LONG) :: Lag_Center_Element !< Lag matrix center point
  INTEGER(LONG) :: Lag_Center_Line!< Lag matrix center point
  INTEGER(LONG) :: Total_Points_Lag



  ! DBSCAN variables
  INTEGER(LONG) MinPts
  INTEGER(LONG) :: NumberOfClusters
  INTEGER(LONG) :: MaxCluster
  INTEGER(LONG) :: MaxClusterSize
  REAL(SINGLE) :: Eps
  REAL(SINGLE), ALLOCATABLE, DIMENSION(:,:) :: SetOfPoints

  ! Minimum SSD value and location
  REAL(DOUBLE) :: Min_SSD_Value
  INTEGER(LONG) :: Min_SSD_Location
  INTEGER(LONG) :: Min_SSD_Element
  INTEGER(LONG) :: Min_SSD_Line

  INTEGER(LONG) :: Match_Lag_Element
  INTEGER(LONG) :: Match_Lag_Line

  ! Offset for location of match area within search grid
  INTEGER(LONG) :: Match_Offset

  INTEGER(LONG) :: Element_Counter_Target
  INTEGER(LONG) :: Target_Index
  INTEGER(LONG) :: Element_Counter_Search
  INTEGER(LONG) :: Element
  INTEGER(LONG) :: Line
  INTEGER(LONG) :: Search_Index


  REAL(DOUBLE) :: Sum_Target
  REAL(DOUBLE) :: Sum_Squared_Target
  REAL(DOUBLE) :: Sum_Squared_Match
  REAL(DOUBLE) :: Target_Scene_Variance
  REAL(DOUBLE) :: Match_Scene_Variance
  INTEGER(LONG) :: Total_Points_Search

  ! statistics for 5x5 target and match scenes
  REAL(DOUBLE) :: Sum_Product
  REAL(DOUBLE) :: Sum_Diff_Squared_Target
  REAL(DOUBLE) :: Sum_Diff_Squared_Match
  REAL(DOUBLE) :: Small_Target_Mean
  REAL(DOUBLE) :: Small_Match_Mean
  REAL(DOUBLE) :: Small_Target_StdDev
  REAL(DOUBLE) :: Small_Match_StdDev
  REAL(DOUBLE) :: Corr_Coeff_5x5

  REAL(DOUBLE) :: Corr_Coeff_5x5_Thresh

  REAL(SINGLE) Relative_Ele_Disp
  REAL(SINGLE) Relative_Line_Disp
  REAL(SINGLE) :: Fractional_Element
  REAL(SINGLE) :: Fractional_Line


  INTEGER(LONG) :: Total_Points_Target
  INTEGER(LONG) Small_Target_Element_Size
  INTEGER(LONG) Small_Target_Line_Size
  INTEGER(LONG) Number_Of_Loops
  INTEGER(LONG) Small_Search_Element_Size
  INTEGER(LONG) Small_Search_Line_Size
  INTEGER(LONG) Points_Small_Target
  INTEGER(LONG) Points_Small_Search
  INTEGER(LONG) Counter
  INTEGER(LONG) Good_Sample, good_sample_rev, good_sample_x2, idx1, idx2, mxsize_rev, mxsize_for  ! aab Tandem (dbscan)
  INTEGER(LONG) Initial_Sample
  INTEGER(LONG) Good_CTP_Sample
  INTEGER(LONG) sample_size_rev, sample_size_for  ! aab Tandem (dbscan) 
  INTEGER(LONG) Good_BT_Sample
  INTEGER(LONG) Center_Element
  INTEGER(LONG) Center_Line
  INTEGER(LONG) Center_Element_FS
  INTEGER(LONG) Center_Line_FS
  INTEGER(LONG) Center_Element_SS
  INTEGER(LONG) Center_Line_SS
  INTEGER(LONG) Point
  INTEGER(LONG) Diff_In_Size
  INTEGER(LONG) Pixel_Offset
  INTEGER(LONG) Pixel_Offset_Search
  ! Offset between small search center and large search center
  INTEGER(LONG) Ele_Shift_Search
  INTEGER(LONG) Line_Shift_Search


  INTEGER(LONG) I_Upper, I1, I2, I3, I4, L1

  ! Array holding the phase type counts of largest cluster
  ! (0=clear, 1=liquid, 2=supercooled, 3=mixed, 4=ice)

  INTEGER phase_counts(4), dominant_phase

  ! Add variables for initial sample
  REAL(DOUBLE) Sum_Delta_X_Initial
  REAL(DOUBLE) Sum_Delta_Y_Initial
  REAL(DOUBLE) Avg_Delta_X_Initial
  REAL(DOUBLE) Avg_Delta_Y_Initial

  REAL(DOUBLE), DIMENSION(:), ALLOCATABLE :: Delta_X
  REAL(DOUBLE), DIMENSION(:), ALLOCATABLE :: Delta_Y
  REAL(DOUBLE), DIMENSION(:), ALLOCATABLE :: Delta_X_Final
  REAL(DOUBLE), DIMENSION(:), ALLOCATABLE :: Delta_Y_Final
  REAL(DOUBLE), DIMENSION(:), ALLOCATABLE :: Dist
  REAL(SINGLE), DIMENSION(:), ALLOCATABLE :: Small_Target
  REAL(SINGLE), DIMENSION(:), ALLOCATABLE :: Small_Search
  REAL(DOUBLE), DIMENSION(:), ALLOCATABLE :: Small_Match

  REAL(DOUBLE) Sum_Delta_X
  REAL(DOUBLE) Sum_Delta_Y
  REAL(DOUBLE) Sum_Dist
  REAL(DOUBLE) Sum_Dist_Squared
  REAL(DOUBLE) Avg_Delta_X
  REAL(DOUBLE) Avg_Delta_Y
  REAL(DOUBLE) Avg_Dist
  REAL(DOUBLE) StdDev_Delta_X
  REAL(DOUBLE) StdDev_Delta_Y
  REAL(DOUBLE) StdDev
  REAL(DOUBLE) PctOfAvg
  REAL(DOUBLE) Mean_Abs_Dev
  REAL(DOUBLE) Cutoff_1SD
  REAL(DOUBLE) Median_Delta_X
  REAL(DOUBLE) Median_Delta_Y

  REAL(DOUBLE) :: Corr_Coeff_DB

  integer i,j



  ! Make sure search array is larger than target array. Return control if
  ! it isn't.

  Target_Element_Size = SIZE(Box_Data%BrtTemp, 1)
  Target_Line_Size = SIZE(Box_Data%BrtTemp, 2)
  Total_Points_Target = Target_Element_Size * Target_Line_Size
  Search_Element_Size = SIZE(Search_BT_Values, 1)
  Search_Line_Size = SIZE(Search_BT_Values, 2)
  Total_Points_Search = Search_Element_Size * Search_Line_Size


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

  Line_Disp = 0.0
  Element_Disp = 0.0
  Least_Square_Error = 0
  Corr_Coeff = 0.0
  Target_Scene_Mean = 0.0
  Target_Scene_StdDev = 0.0
  Match_Scene_Mean = 0.0
  Match_Scene_StdDev = 0.0
  Sum_Target = 0.0
  Sum_Squared_Target = 0.0
  Sum_Squared_Match = 0.0
  Sum_Delta_X_Initial = 0.0_DOUBLE
  Sum_Delta_Y_Initial = 0.0_DOUBLE
  Sum_Delta_X = 0.0
  Sum_Delta_Y = 0.0
  Avg_Delta_X = 0.0
  Avg_Delta_Y = 0.0
  Median_Delta_X = 0.0
  Median_Delta_Y = 0.0
  StdDev_Delta_X = 0.
  StdDev_Delta_Y = 0.
  Target_Scene_Variance = 0.0
  Match_Scene_Variance = 0.0
  Small_Target_Element_Size = 5
  Small_Target_Line_Size = 5
  Points_Small_Target = Small_Target_Element_Size * Small_Target_Line_Size
  Counter = 0
  Good_Sample = 0
  good_sample_x2 = 0
  
  ! aab Tandem (dbscan): keep track of combined sample if processing forward time step
  if (vector .eq. 2) then
          good_sample_rev = DBSCAN_Out%Initial_Sample  ! number of initial points from reverse time step
          good_sample_x2 = DBSCAN_Out%Initial_Sample   ! cumulative total from both time steps
  endif
  if (vector .eq. 1) Status_Flag = ZERO_INITIAL_SAMPLE

  StdDev = 0.
  PctOfAvg = 0.
  Sum_Dist = 0.
  Sum_Dist_Squared = 0.
  Avg_Dist = 0.
  Mean_Abs_Dev = 0.
  Cutoff_1SD = 0.
  Center_Element = (Target_Element_Size + 1) / 2
  Center_Line = (Target_Line_Size + 1) / 2
  Pixel_Offset = Small_Target_Element_Size / 2
  Center_Element_FS = (Search_Element_Size + 1) / 2
  Center_Line_FS = (Search_Line_Size + 1) / 2

  ALLOCATE(Target_1D(Total_Points_Target), Search_1D(Total_Points_Search))

  !-----------------------------------------------------------------------------
  ! Reshape 2-D box data into vectors for pattern matching code.
  !-----------------------------------------------------------------------------

  Target_1D = RESHAPE(Box_Data%BrtTemp, (/Total_Points_Target/))
  Search_1D = RESHAPE(Search_BT_Values, (/Total_Points_Search/))

  ! ----------------------------------------------------------------------------
  ! Specify "lag" coefficient matrix size. This array contains the SSD
  ! values for all possible lag positions of the template within the
  ! search region. This array is filled regardless of the tracking
  ! method because it is used by both the LP norm and CC methods.
  !
  ! Note: The lag size is the same whether or not nested tracking is used.
  !-----------------------------------------------------------------------------

  Lag_Element_Size = Search_Element_Size - Target_Element_Size + 1
  Lag_Line_Size = Search_Line_Size - Target_Line_Size + 1
  Lag_Center_Element = (Lag_Element_Size + 1) / 2
  Lag_Center_Line = (Lag_Line_Size + 1) / 2
  Total_Points_Lag = Lag_Element_Size * Lag_Line_Size

  Small_Search_Element_Size = Small_Target_Element_Size + Lag_Element_Size - 1
  Small_Search_Line_Size = Small_Target_Line_Size + Lag_Line_Size - 1
  Points_Small_Search = Small_Search_Element_Size * Small_Search_Line_Size
  Number_Of_Loops = (Target_Element_Size - 2*Pixel_Offset) * (Target_Line_Size - 2*Pixel_Offset)
  Pixel_Offset_Search = Small_Search_Element_Size / 2

  ! ----------------------------------------------------------------------------
  ! Need to compute difference in size between the large target box and the
  ! small target box. This is also the difference in size between the large
  ! search region and the small search region.
  ! ----------------------------------------------------------------------------

  Diff_In_Size = Target_Element_Size - Small_Target_Element_Size

  ALLOCATE(Lag_Matrix(Total_Points_Lag))
  ALLOCATE(Small_Target(Points_Small_Target))
  ALLOCATE(Small_Search(Points_Small_Search))
  ALLOCATE(Small_Match(Points_Small_Target), source=MISSING_VALUE_DOUBLE)
  ALLOCATE(Delta_X(Number_Of_Loops), Delta_Y(Number_Of_Loops))
  !ALLOCATE(SetOfPoints(Number_Of_Loops, DBSCAN_Output_Variables))
  ALLOCATE(SetOfPoints(size(SetOfPoints_Init, dim=1), size(SetOfPoints_Init, dim=2)), source=MISSING_VALUE_SINGLE)

  ! For each interior pixel of the larger target box, find the local motion
  ! associated with a smaller 5x5 scene. In effect, we are "nesting" a 5x5
  ! box within the larger box containing the entire target scene. This
  ! produces a sample of local motion vectors.

  DO  Line = Pixel_Offset + 1, Target_Line_Size - Pixel_Offset

    DO Element = Pixel_Offset + 1, Target_Element_Size - Pixel_Offset

      Counter = Counter + 1
      Delta_X(Counter) = MISSING_VALUE_REAL4
      Delta_Y(Counter) = MISSING_VALUE_REAL4
      Sum_Product = 0.0
      Small_Target_Mean = 0.0
      Small_Match_Mean = 0.0
      Small_Target_StdDev = 0.0
      Small_Match_StdDev = 0.0
      Corr_Coeff_5x5 = 0.0

      ! compute center element of small search box
      Center_Element_SS = Element + ((Small_Search_Element_Size + 1)/2) - 3
      ! center line of small search box
      Center_Line_SS = Line + ((Small_Search_Line_Size + 1)/2) - 3

      Ele_Shift_Search = Center_Element_SS - Center_Element_FS
      Line_Shift_Search = Center_Line_SS - Center_Line_FS

      ! wcb 11/19/13
      ! filter out missing values
      ! wcb 11/19/15
      ! modified check to filter out values exceeding 1100
      ! associated with negative heights
      ! wcb - 1/25/16
      ! only apply check if processing cloudy targets

      IF (CSWV_Flag .EQ. sym%NO) THEN

        IF (Box_Data%Cld_Top_Press(Element, Line) .EQ. MISSING_VALUE_REAL4 .OR. &
            Box_Data%Cld_Top_Press(Element, Line) .GT. CTP_MAX) CYCLE

        ! WCB - also check cloud height quality flag
        ! AAB - Now skipping this check. Allow AMVs from the below QF values through.

        !IF (Box_Data%CldHgtQF(Element, Line) .EQ. MISSING_VALUE_INT1 .OR. &

        !    Box_Data%CldHgtQF(Element, Line) .EQ. CTH_DQF_MARGINAL_RETREVIAL .OR.  &

        !    Box_Data%CldHgtQF(Element, Line) .EQ. CTH_DQF_RETREVIAL_ATTEMPTED .OR.  &

        !    Box_Data%CldHgtQF(Element, Line) .EQ. CTH_DQF_BAD_RETREVIAL) CYCLE

      ENDIF

      ! fill small target

      I_Upper = Small_Target_Element_Size - 1
      DO I = 0, I_Upper

        I1 = I*Small_Target_Element_Size+1
        I2 = I1 + Small_Target_Element_Size - 1
        I3 = Element-Pixel_Offset
        I4 = Element+Pixel_Offset
        L1 = Line-Pixel_Offset+I
        Small_Target(I1:I2) = Box_Data%BrtTemp(I3:I4, L1)

      ENDDO

      ! Subsample the large search region to fill small search region.

      I_Upper = Small_Search_Element_Size - 1
      DO I = 0, I_Upper

        I1 = I*Small_Search_Element_Size+1
        I2 = I1 + Small_Search_Element_Size - 1
        I3 = Center_Element_SS-Pixel_Offset_Search
        I4 = Center_Element_SS+Pixel_Offset_Search
        L1 = Center_Line_SS-Pixel_Offset_Search+I
        Small_Search(I1:I2) = Search_BT_Values(I3:I4, L1)

      ENDDO

      CALL Fill_Lag_Matrix_1D(Lag_Element_Size, Diff_In_Size, Small_Target, &
                                      Small_Search, Lag_Table2, Lag_Matrix)

      ! ----------------------------------------------------------------------------
      ! Determine location of match within lag matrix.
      ! ----------------------------------------------------------------------------

      Min_SSD_Value = MINVAL(Lag_Matrix)
      Least_Square_Error = Min_SSD_Value
      Min_SSD_Location = MINLOC(Lag_Matrix, 1)
      Min_SSD_Element = MOD(Min_SSD_Location, Lag_Element_Size)

      IF (MOD(Min_SSD_Location, Lag_Element_Size) .EQ. 0) THEN
        Min_SSD_Element = Lag_Element_Size
        Min_SSD_Line = (Min_SSD_Location / Lag_Element_Size)
      ELSE
        Min_SSD_Line = (Min_SSD_Location / Lag_Element_Size) + 1
      ENDIF

      ! ----------------------------------------------------------------------------
      ! Do not accept match if located on boundary. This condition may indicate the
      ! true match is located beyond the domain of the search.
      ! ----------------------------------------------------------------------------

      IF (Min_SSD_Element .EQ. 1 .OR. Min_SSD_Element .EQ. Lag_Element_Size .OR. &
                 Min_SSD_Line .EQ. 1 .OR. Min_SSD_Line .EQ. Lag_Line_Size) CYCLE

      ! Fill small match array and compute correlation

      Match_Offset = (Min_SSD_Line - 1) * Small_Search_Element_Size &
                                                  + Min_SSD_Element

      I_Upper = Small_Target_Element_Size - 1
      DO I = 0, I_Upper

        I1 = I*Small_Target_Element_Size + 1 !row1,column1
        I2 = I1 + Small_Target_Element_Size - 1 !row2,column2
        I3 = Match_Offset
        I4 = I3 + Small_Target_Element_Size - 1
        Small_Match(I1:I2) = Small_Search(I3:I4)

        Match_Offset = Match_Offset + Small_Search_Element_Size

      ENDDO

      ! discard match if missing values are present

      IF (ANY(Small_Match == MISSING_VALUE_REAL4)) CYCLE

      Small_Target_Mean = SUM(Small_Target) / REAL(Points_Small_Target)
      Small_Match_Mean = SUM(Small_Match) / REAL(Points_Small_Target)

      Sum_Product = 0.0
      Sum_Diff_Squared_Target = 0.0
      Sum_Diff_Squared_Match = 0.0

      DO Point=1,Points_Small_Target

        Sum_Product = Sum_Product + (Small_Target(Point) - Small_Target_Mean) &
                                    * (Small_Match(Point) - Small_Match_Mean)

        Sum_Diff_Squared_Target = Sum_Diff_Squared_Target &
           + (Small_Target(Point) - Small_Target_Mean)**2

        Sum_Diff_Squared_Match = Sum_Diff_Squared_Match &
           + (Small_Match(Point) - Small_Match_Mean)**2

      ENDDO

      Small_Target_StdDev = SQRT(Sum_Diff_Squared_Target)
      Small_Match_StdDev = SQRT(Sum_Diff_Squared_Match)

      Corr_Coeff_5x5= Sum_Product / (Small_Target_StdDev * Small_Match_StdDev)

      ! use appropriate correlation threshold - 1/25/16

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

        Corr_Coeff_5x5_Thresh = CORR_COEFF_5X5_POLAR_CSWV

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

        ! aab Tandem mod
        !Corr_Coeff_5x5_Thresh = CORR_COEFF_5X5_POLAR_CLD

        ! Tandem VIIRS Polar winds Corr Coeff is set based on triplet configuration.
        SELECT CASE(TandemTriplet)
           CASE ('S') ! Single-Sat (default)
              Corr_Coeff_5x5_Thresh = CORR_COEFF_5X5_POLAR_CLD_S
           CASE ('D') ! Dual-Sat
              Corr_Coeff_5x5_Thresh = CORR_COEFF_5X5_POLAR_CLD_D
           CASE ('T') ! Three-Sat
              Corr_Coeff_5x5_Thresh = CORR_COEFF_5X5_POLAR_CLD_T
           CASE DEFAULT
              Corr_Coeff_5x5_Thresh = CORR_COEFF_5X5_POLAR_CLD_S
        END SELECT

      ELSE

        Corr_Coeff_5x5_Thresh = CORR_COEFF_5X5_GEO

      ENDIF

      ! filter out matches having low correlation scores
      ! introduced Corr_Coeff_5x5_Thresh - 1/25/16

      IF (Corr_Coeff_5x5 .LT. Corr_Coeff_5x5_Thresh)   CYCLE

      Good_Sample = Good_Sample + 1
      ! aab Tandem (dbscan): increment combined sample at the same time
      good_sample_x2 = good_sample_x2 + 1

      ! Compute displacement of match relative to lag center point

      Element_Disp = REAL(Min_SSD_Element - Lag_Center_Element)
      Line_Disp = REAL(Min_SSD_Line - Lag_Center_Line)
      Match_Lag_Element = Min_SSD_Element
      Match_Lag_Line = Min_SSD_Line

      CALL Interpolate_Disp(Match_Lag_Element, Match_Lag_Line, &
                          Fractional_Element, Fractional_Line)

      ! Add fractional values to integer displacement.

      Element_Disp = Element_Disp + Fractional_Element
      Line_Disp = Line_Disp + Fractional_Line

      ! Compute absolute displacement taking into account the displacement of
      ! the small target box relative to the center of the larger box. The
      ! absolute displacement is needed for the DBSCAN analysis.

      Delta_X(Counter) = REAL(Center_Element) + Delta_X_Sat &
                           + Element_Disp - Element

      Delta_Y(Counter) = REAL(Center_Line) + Delta_Y_Sat &
                              + Line_Disp - Line

      ! account for offset between small search and large search

      Delta_X(Counter) = Delta_X(Counter) + Ele_Shift_Search
      Delta_Y(Counter) = Delta_Y(Counter) + Line_Shift_Search

      ! ------------------------------------------------------------------------
      ! Array of points for DBSCAN
      ! This is the initial set of points that is analyzed for motion clusters
      ! ------------------------------------------------------------------------

      !SetOfPoints(Good_Sample, 1) = Delta_X(Counter)

      ! Multiply element displacement by normalization factor

      SetOfPoints(Good_Sample, 1) = Delta_X(Counter) * Factor
      SetOfPoints(Good_Sample, 2) = Delta_Y(Counter)
      SetOfPoints(Good_Sample, 5) = REAL(Element)
      SetOfPoints(Good_Sample, 6) = REAL(Line)
      SetOfPoints(Good_Sample, 7) = Box_Data%Cld_Top_Press(Element, Line)
      SetOfPoints(Good_Sample, 8) = Box_Data%Cld_Hgt(Element, Line)
      SetOfPoints(Good_Sample, 9) = Box_Data%Press_Error(Element, Line)
      SetOfPoints(Good_Sample, 10) = Box_Data%Temp_Error(Element, Line)
      SetOfPoints(Good_Sample, 11) = Box_Data%Cld_Top_Temp(Element, Line)
      SetOfPoints(Good_Sample, 12) = Corr_Coeff_5x5
      SetOfPoints(Good_Sample, 13) = Box_Data%Cloud_Phase(Element, Line)
      SetOfPoints(Good_Sample, 14) = Box_Data%Cloud_Type(Element, Line)
      SetOfPoints(Good_Sample, 15) = Box_Data%Cld_OD(Element, Line)
      SetOfPoints(Good_Sample, 16) = Box_Data%ACHA_Cost(Element, Line)
      SetOfPoints(Good_Sample, 17) = vector  ! aab Tandem (dbscan)

      ! fill 2-D arrays for plotting

      DBSCAN_Out%Delta_X_Init(Element, Line) = Delta_X(Counter)
      DBSCAN_Out%Delta_Y_Init(Element, Line) = Delta_Y(Counter)

    END DO        ! element loop

  END DO          ! line loop

  ! Call DBSCAN algorithm

  MinPts = 4       ! minimum number of points in a cluster

  ! aab Tandem (dbscan)
  ! Specify neighborhood radius parameter (epsilon), use larger value for ABI visible winds
  IF (Band .EQ. NATIVE_ABI_CHN2) THEN
    Eps = EPS_GEO_VZ
  ELSE 
    Eps = EPS_GEO    ! neighborhood radius
  ENDIF

  ! use larger radius with polar winds
  ! EPS is the neighborhood radius parameter in pixels
  ! A larger value (3.3) is used in polar winds
  ! processing because of the larger time interval involved

  IF (TRIM(NavType) .EQ. 'PS') Eps = EPS_POLAR

  if (vector .eq. 1 .and. good_sample .gt. 1) then  ! aab Tandem (dbscan)
          status_flag = sym%SUCCESS 
          setofpoints_init = setofpoints
          dbscan_out%initial_sample = good_sample
  endif

  ! aab Tandem (dbscan): Only call DBSCAN after forward time step is processed (i.e., when vector=2)

  if (vector .eq. 2) then  ! aab Tandem (dbscan)
    idx1 = good_sample + 1
    idx2 = good_sample_x2

    ! aab Tandem (dbscan): reverse line/ele displacements to account for reverse time step
    setofpoints_init(1:good_sample_rev,1) = -setofpoints_init(1:good_sample_rev,1)
    setofpoints_init(1:good_sample_rev,2) = -setofpoints_init(1:good_sample_rev,2)

    ! aab Tandem (dbscan): combine reverse time step sample with forward time step sample before calling DBSCAN
    setofpoints(idx1:idx2,:) = setofpoints_init(1:good_sample_rev,:)
    good_sample = good_sample_x2

    CALL DBSCAN(Good_Sample, SetOfPoints, MinPts, Eps, &
           NumberOfClusters, MaxCluster, MaxClusterSize)

    ! save initial sample of points

    IF ( SIZE(SetOfPoints_Init, 1) .NE. SIZE(SetOfPoints, 1) .OR. &
         SIZE(SetOfPoints_Init, 2) .NE. SIZE(SetOfPoints, 2) ) THEN
       CALL fw_log_error(FUNC, 'Shape of SetOfPoints_Init does not match SetOfPoints')
    ENDIF

    SetOfPoints_Init = SetOfPoints

    IF (Good_Sample .GT. 1) THEN

      Avg_Delta_X = 0.0
      Avg_Delta_Y = 0.0
      Sum_Delta_X = 0.0
      Sum_Delta_Y = 0.0
      Initial_Sample = Good_Sample
      DBSCAN_Out%Initial_Sample = Good_Sample

      Good_Sample = 0
      Good_CTP_Sample = 0
      Good_BT_Sample = 0
      sample_size_rev = 0
      sample_size_for = 0

      ! Allocate memory for arrays to hold final sample
  
      !ALLOCATE(Delta_X_Final(Number_Of_Loops), Delta_Y_Final(Number_Of_Loops))  ! aab Tandem (dbscan)
      !ALLOCATE(Dist(Number_Of_Loops))  ! aab Tandem (dbscan)
      ! aab Tandem (dbscan): double size for combined sample processing
      ALLOCATE(Delta_X_Final(Number_Of_Loops*2), Delta_Y_Final(Number_Of_Loops*2))
      ALLOCATE(Dist(Number_Of_Loops*2))

      ! determine dominant phase of largest cluster
      phase_counts(1) = COUNT(SetOfPoints(:,13) == 1. .AND. SetOfPoints(:,4) == REAL(MaxCluster))
      phase_counts(2) = COUNT(SetOfPoints(:,13) == 2. .AND. SetOfPoints(:,4) == REAL(MaxCluster))
      phase_counts(3) = COUNT(SetOfPoints(:,13) == 3. .AND. SetOfPoints(:,4) == REAL(MaxCluster))
      phase_counts(4) = COUNT(SetOfPoints(:,13) == 4. .AND. SetOfPoints(:,4) == REAL(MaxCluster))
      dominant_phase = MAXLOC(phase_counts, DIM=1)

      ! aab Tandem (dbscan): count the number of samples for each vector
      mxsize_rev = COUNT(SetOfPoints(:,17) == 1. .AND. SetOfPoints(:,4) == REAL(MaxCluster))
      mxsize_for = COUNT(SetOfPoints(:,17) == 2. .AND. SetOfPoints(:,4) == REAL(MaxCluster))

      ! use output of DBSCAN analysis to determine final sample

      DO Point = 1, Initial_Sample

        ! Add initial sample variables
        Sum_Delta_X_Initial = Sum_Delta_X_Initial + SetOfPoints(Point,1) / Factor
        Sum_Delta_Y_Initial = Sum_Delta_Y_Initial + SetOfPoints(Point,2)

        ! if point is not part of largest cluster or is missing BT or CTP, skip it

        IF (CSWV_Flag .EQ. sym%YES) THEN

          ! do not check CTP value for CSWV winds only max cluster - 1/25/16

          IF (SetOfPoints(Point,4) .NE. REAL(MaxCluster)) CYCLE

        ELSE

          ! check CTP value for all cloudy winds and max cluster - 1/25/16

          IF (SetOfPoints(Point,4) .NE. REAL(MaxCluster) .OR. &
                SetOfPoints(Point,7) .EQ. MISSING_VALUE_REAL4) CYCLE

        ENDIF

        ! if channel is long-wave IR and cloud top height is the surface, skip it

        IF (Band .EQ. 14 .AND. SetOfPoints(Point,8) .EQ. 0.0) CYCLE

        ! check phase of SWIR, WV and LWIR bands

        IF (Band .GE. 7) THEN

          ! if dominant phase is ice, exclude liquid pixels and vice versa

          IF (dominant_phase .EQ. sym%ICE_PHASE .AND. INT(SetOfPoints(Point,13)) .EQ. sym%WATER_PHASE) CYCLE
          IF (dominant_phase .EQ. sym%WATER_PHASE .AND. INT(SetOfPoints(Point,13)) .EQ. sym%ICE_PHASE) CYCLE

        ENDIF

        Good_Sample = Good_Sample + 1

        ! Remove pixel factor before computing final solution

        Sum_Delta_X = Sum_Delta_X + SetOfPoints(Point,1) / Factor
        Sum_Delta_Y = Sum_Delta_Y + SetOfPoints(Point,2)
        Delta_X_Final(Good_Sample) = SetOfPoints(Point,1) / Factor
        Delta_Y_Final(Good_Sample) = SetOfPoints(Point,2)

        ! Fill output structure

        Good_CTP_Sample = Good_CTP_Sample + 1   ! now represents total from both time steps 

        ! Save final sample of points in output structure

        !IF (Vector .EQ. REVERSE) THEN  ! aab Tandem (dbscan)
        if(setofpoints(point,17) == 1.0) then

          sample_size_rev = sample_size_rev + 1
          DBSCAN_Out%Number_CTP_Sample1 = sample_size_rev 
          DBSCAN_Out%CTP_5x5_Sample1(sample_size_rev) = SetOfPoints(Point,7)
          DBSCAN_Out%CldTemp_5x5_Sample1(sample_size_rev) = SetOfPoints(Point,11)
          DBSCAN_Out%Cld_Hgt_5x5_Sample1(sample_size_rev) = SetOfPoints(Point,8)
          DBSCAN_Out%Press_Error_Sample1(sample_size_rev) = SetOfPoints(Point,9)
          DBSCAN_Out%Temp_Error_Sample1(sample_size_rev) = SetOfPoints(Point,10)
          DBSCAN_Out%CORR_5x5_Sample1(sample_size_rev) = SetOfPoints(Point,12)
          DBSCAN_Out%Cld_OD_Sample1(sample_size_rev) = SetOfPoints(Point,15)
          DBSCAN_Out%ACHA_Cost_Sample1(sample_size_rev) = SetOfPoints(Point,16)
          DBSCAN_Out%Delta_X_Reverse(sample_size_rev) = SetOfPoints(Point,1) / Factor
          DBSCAN_Out%Delta_Y_Reverse(sample_size_rev) = SetOfPoints(Point,2)

          DBSCAN_Out%Phase_Sample1(sample_size_rev) = SetOfPoints(Point, 13)
          DBSCAN_Out%Type_Sample1(sample_size_rev) = SetOfPoints(Point, 14)

        !ELSE IF (Vector .EQ. FORWARD) THEN  ! aab Tandem (dbscan)
        else if(setofpoints(point,17) == 2.0) then

          sample_size_for = sample_size_for + 1
          DBSCAN_Out%Number_CTP_Sample2 = sample_size_for 
          DBSCAN_Out%CTP_5x5_Sample2(sample_size_for) = SetOfPoints(Point,7)
          DBSCAN_Out%CldTemp_5x5_Sample2(sample_size_for) = SetOfPoints(Point,11)
          DBSCAN_Out%Cld_Hgt_5x5_Sample2(sample_size_for) = SetOfPoints(Point,8)
          DBSCAN_Out%Press_Error_Sample2(sample_size_for) = SetOfPoints(Point,9)
          DBSCAN_Out%Temp_Error_Sample2(sample_size_for) = SetOfPoints(Point,10)
          DBSCAN_Out%CORR_5x5_Sample2(sample_size_for) = SetOfPoints(Point,12)
          DBSCAN_Out%Cld_OD_Sample2(sample_size_for) = SetOfPoints(Point,15)
          DBSCAN_Out%ACHA_Cost_Sample2(sample_size_for) = SetOfPoints(Point,16)
          DBSCAN_Out%Delta_X_Forward(sample_size_for) = SetOfPoints(Point,1) / Factor
          DBSCAN_Out%Delta_Y_Forward(sample_size_for) = SetOfPoints(Point,2)

          DBSCAN_Out%Phase_Sample2(sample_size_for) = SetOfPoints(Point, 13)
          DBSCAN_Out%Type_Sample2(sample_size_for) = SetOfPoints(Point, 14)

        ENDIF

        ! keep track of sub-vector number
        DBSCAN_Out%vector(good_ctp_sample) = SetOfPoints(Point, 17)

        ! fill 2D arrays for plotting

        I = INT(SetOfPoints(Point,5))
        J = INT(SetOfPoints(Point,6))

        DBSCAN_Out%Delta_X_Final(I,J) = SetOfPoints(Point,1) / Factor
        DBSCAN_Out%Delta_Y_Final(I,J) = SetOfPoints(Point,2)

        ! sum up CTP sample for statistics

! aab Tandem (dbscan)
!        Mean_CTP = Mean_CTP + SetOfPoints(Point,7)
!        Mean_CTP2 = Mean_CTP2 + SetOfPoints(Point,7)**2

      ENDDO ! Loop over initial sample

      ! Check final sample size. Fail if zero, otherwise compute
      ! average displacement.

      IF (Good_Sample .GT. 0) THEN

        Avg_Delta_X = Sum_Delta_X / REAL(Good_Sample)
        Avg_Delta_Y = Sum_Delta_Y / REAL(Good_Sample)

        ! Add initial sample variables
        Avg_Delta_X_Initial = Sum_Delta_X_Initial / REAL(Initial_Sample)
        Avg_Delta_Y_Initial = Sum_Delta_Y_Initial / REAL(Initial_Sample)

! aab Tandem (dbscan) 
!        Mean_CTP =  Mean_CTP / Good_CTP_Sample
!        Mean_CTP2 =  Mean_CTP2 / Good_CTP_Sample
!        SD_CTP = SQRT(Mean_CTP2 - Mean_CTP**2)

        Status_Flag = sym%SUCCESS

        ! Compute spread about mean solution

        DO Point = 1, Good_Sample

          Dist(Point) = SQRT((Delta_X_Final(Point) - Avg_Delta_X)**2 &
                          + (Delta_Y_Final(Point) - Avg_Delta_Y)**2)

          Sum_Dist = Sum_Dist + Dist(Point)

        ENDDO

        Avg_Dist = Sum_Dist / REAL(Good_Sample)

        DO Point = 1, Good_Sample

          Sum_Dist_Squared = Sum_Dist_Squared + (Dist(Point) - Avg_Dist)**2

        ENDDO

        StdDev = SQRT(Sum_Dist_Squared / REAL(Good_Sample))
        PctOfAvg = StdDev / SQRT(Avg_Delta_X**2 + Avg_Delta_Y**2)
        ! keep original code -- wcb 04/19/13
        ! save ratio of final sample to initial sample instead
        !PctOfAvg = REAL(Good_Sample) / REAL(Box_Data%CloudyPixels)

        ! Fill scalar values of output structure

        !IF (Vector .EQ. REVERSE) THEN  ! aab Tandem (dbscan)

          !DBSCAN_Out%StdDev_MVD1 = StdDev  ! aab Tandem (dbscan)
          !DBSCAN_Out%PctOfAvg1 = PctOfAvg  ! aab Tandem (dbscan)
          DBSCAN_Out%NumberOfClusters1 = NumberOfClusters
          DBSCAN_Out%MaxCluster1 = MaxCluster
          DBSCAN_Out%MaxClusterSize1 = MaxClusterSize
          !DBSCAN_Out%Max1_PctOfSample = REAL(MaxClusterSize) / REAL(Initial_Sample)  ! aab Tandem (dbscan)


        !ELSE IF (Vector .EQ. FORWARD) THEN  ! aab Tandem (dbscan)

          !DBSCAN_Out%StdDev_MVD2 = StdDev  ! aab Tandem (dbscan)
          !DBSCAN_Out%PctOfAvg2 = PctOfAvg  ! aab Tandem (dbscan)
          DBSCAN_Out%NumberOfClusters2 = NumberOfClusters
          DBSCAN_Out%MaxCluster2 = MaxCluster
          DBSCAN_Out%MaxClusterSize2 = MaxClusterSize
          !DBSCAN_Out%Max2_PctOfSample = REAL(MaxClusterSize) / REAL(Initial_Sample)  ! aab Tandem (dbscan)

        !ENDIF  ! aab Tandem (dbscan)

        ! deallocate arrays holding final sample of points

        DEALLOCATE(Delta_X_Final, Delta_Y_Final, Dist)

      ELSE  ! If final sample is zero flag target as failure

        Status_Flag = ZERO_FINAL_SAMPLE
        DEALLOCATE(Delta_X_Final, Delta_Y_Final, Dist)

      ENDIF

    ELSE  ! If initial sample is zero flag target as failure

      Status_Flag = ZERO_INITIAL_SAMPLE

    ENDIF
  endif ! combined sample testing

  ! Convert absolute displacement back to relative displacement from center.

  Element_Disp = Avg_Delta_X
  Line_Disp = Avg_Delta_Y
  Relative_Ele_Disp = Avg_Delta_X - Delta_X_Sat
  Relative_Line_Disp = Avg_Delta_Y - Delta_Y_Sat

  !-----------------------------------------------------------------------------
  ! Offset for location of match area within search grid.
  !-----------------------------------------------------------------------------

  Min_SSD_Line = Lag_Center_Line + NINT(Relative_Line_Disp)
  Min_SSD_Element = Lag_Center_Element + NINT(Relative_Ele_Disp)

  ! guard against exceeding bounds of search region
  IF (Min_SSD_Line .LT. 1) Min_SSD_Line = 1
  IF (Min_SSD_Element .LT. 1) Min_SSD_Element = 1
  IF (Min_SSD_Line .GT. Lag_Line_Size) Min_SSD_Line = Lag_Line_Size
  IF (Min_SSD_Element .GT. Lag_Element_Size) Min_SSD_Element = Lag_Element_Size

  Match_Offset = (Min_SSD_Line - 1) * Search_Element_Size + Min_SSD_Element - 1
  Min_SSD_Value = SSD_1D2(0, Target_1D, Search_1D, 1, Match_Offset, INT(1,BYTE))

  Element_Counter_Search = - Search_Element_Size
  Element_Counter_Target = - Target_Element_Size

  line_loop: DO Line = 1, Target_Line_Size

    Element_Counter_Search = Element_Counter_Search + Search_Element_Size
    Element_Counter_Target = Element_Counter_Target + Target_Element_Size

    element_loop: DO Element = 1, Target_Element_Size

      Search_Index = Element_Counter_Search + Element + Match_Offset
      Target_Index = Element_Counter_Target + Element
      Sum_Target = Sum_Target + Target_1D(Target_Index)
      Sum_Squared_Target = Sum_Squared_Target + Target_1D(Target_Index)**2
      Match_Scene_Mean = Match_Scene_Mean + Search_1D(Search_Index)
      Sum_Squared_Match = Sum_Squared_Match + Search_1D(Search_Index)**2

    END DO element_loop

  END DO line_loop

  ! ----------------------------------------------------------------------------
  ! Compute scene statistics.
  ! ----------------------------------------------------------------------------

  Target_Scene_Mean = Sum_Target / REAL(Total_Points_Target)
  Match_Scene_Mean = Match_Scene_Mean / REAL(Total_Points_Target)

  Target_Scene_Variance = Sum_Squared_Target / REAL(Total_Points_Target) - &
                                                      Target_Scene_Mean**2

  Match_Scene_Variance = Sum_Squared_Match / REAL(Total_Points_Target) - &
                                                     Match_Scene_Mean**2

  Target_Scene_StdDev = SQRT(Target_Scene_Variance)
  Match_Scene_StdDev = SQRT(Match_Scene_Variance)

  ! ----------------------------------------------------------------------------
  ! Compute linear correlation coefficient.
  ! ----------------------------------------------------------------------------
  Corr_Coeff_DB = (Sum_Squared_Match + Sum_Squared_Target - Min_SSD_Value) / &
       (2.0 * Total_Points_Target) - Target_Scene_Mean * Match_Scene_Mean

  Corr_Coeff = Corr_Coeff_DB / (Target_Scene_StdDev * Match_Scene_StdDev)


  IF (Corr_Coeff .LT. 0.) Corr_Coeff = 0.0

  ! Deallocate arrays used by nested tracking and DBSCAN

  DEALLOCATE(Target_1D, Search_1D, Lag_Matrix)
  DEALLOCATE(Small_Target, Small_Search, Small_Match)
  DEALLOCATE(Delta_X, Delta_Y)
  DEALLOCATE(SetOfPoints)

END SUBROUTINE Nested_Tracking

!-------------------------------------------------------------------------------
!
! Name:
!   Fill_Lag_Matrix_1D
!
! Function:
!   Coputes the sum of squared difference values for all target-sized
!   scenes in search region.
!
! Description:
!   This subroutine fills the lag matrix, the array of SSD values of all
!   potential matching scenes in the search region.
!
! Calling sequence:
!   CALL Fill_Lag_Matrix_1D
!
! Inputs:
!   Size_Diff            Size difference in pixels between the nested box and
!                        the outer target box.
!   Target_1D            One-dimensional array of target box values
!   Search_1D            One-dimensional array of search box values
!   Lag_Table2           prioritized list of positions to search (a spiral
!                        centered on the middle of the search region)
!
! Outputs:
!   Lag_Matrix           SSD values of potential matching scenes
!
! Dependencies:
!   Assumes SSD similarity measure is being used for tracking.
!
! Restrictions:
!   None
!
! History:
!   06/01/2009 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE Fill_Lag_Matrix_1D(Lag_Size, Size_Diff, Target_1D, Search_1D, &
                                                 Lag_Table2, Lag_Matrix)

  INTEGER(BYTE) Flag

  INTEGER(LONG) :: Point
  INTEGER(LONG) :: Lag_Index
  INTEGER(LONG) :: Search_Index
  INTEGER(LONG) :: Min_SSD_Location
  INTEGER(LONG) :: Min_SSD_Element
  INTEGER(LONG) :: Min_SSD_Line
  INTEGER(LONG) :: Small_Search_Element_Size
  INTEGER(LONG) :: Total_Points_Lag
  INTEGER(LONG), INTENT(IN) :: Size_Diff
  INTEGER(LONG), INTENT(IN) :: Lag_Size

  ! List of positions to use in spiral search for minimum SSD
  INTEGER(LONG), DIMENSION(:,:), INTENT(IN) :: Lag_Table2

  REAL(SINGLE), DIMENSION(:), INTENT(INOUT) :: Lag_Matrix
  REAL(SINGLE), DIMENSION(:), INTENT(IN) :: Target_1D
  REAL(SINGLE), DIMENSION(:), INTENT(IN) :: Search_1D

  ! ----------------------------------------------------------------------------
  ! Initialize all values to -1
  ! ----------------------------------------------------------------------------

  Lag_Matrix(:) = -1.0
  Small_Search_Element_Size = Search_Element_Size - Size_Diff
  Total_Points_Lag = Lag_Size**2

  ! Initialize flag to terminate SSD early if current minimum is exceeded

  Flag = 0

  lag_point_loop: DO Point = 1, Total_Points_Lag

    Lag_Index = Lag_Table2(Point, 3)
    Search_Index = Lag_Table2(Point, 4)
    Lag_Matrix(Lag_Index) = SSD_1D2(Size_Diff, Target_1D, Search_1D, Point, &
                                                        Search_Index, Flag)

  END DO lag_point_loop

  ! ----------------------------------------------------------------------------
  ! Compute SSD value for 4 pts surrounding minimum. These values are necessary
  ! for the sub-pixel interpolation.
  ! ----------------------------------------------------------------------------

  Flag = 1
  Point = 1
  Min_SSD_Location = MINLOC(Lag_Matrix, 1)

  Min_SSD_Element = MOD(Min_SSD_Location, Lag_Element_Size)
  IF (MOD(Min_SSD_Location, Lag_Element_Size) .EQ. 0) THEN
    Min_SSD_Element = Lag_Element_Size
    Min_SSD_Line = (Min_SSD_Location / Lag_Element_Size)
  ELSE
    Min_SSD_Line = (Min_SSD_Location / Lag_Element_Size) + 1
  ENDIF

  IF (Min_SSD_Element .EQ. 1 .OR. Min_SSD_Element .EQ. Lag_Element_Size .OR. &
            Min_SSD_Line .EQ. 1 .OR. Min_SSD_Line .EQ. Lag_Line_Size) RETURN

  ! Top point

  Lag_Index = Min_SSD_Location - Lag_Element_Size
  Search_Index = (Min_SSD_Line - 2) * Small_Search_Element_Size + Min_SSD_Element
  Lag_Matrix(Lag_Index) = SSD_1D2(Size_Diff, Target_1D, Search_1D, Point, &
                                                      Search_Index, Flag)

  ! Left point

  Lag_Index = Min_SSD_Location - 1
  Search_Index = (Min_SSD_Line - 1) * Small_Search_Element_Size + Min_SSD_Element - 1
  Lag_Matrix(Lag_Index) = SSD_1D2(Size_Diff, Target_1D, Search_1D, Point, &
                                                      Search_Index, Flag)

  ! Right point

  Lag_Index = Min_SSD_Location + 1
  Search_Index = (Min_SSD_Line - 1) * Small_Search_Element_Size + Min_SSD_Element + 1
  Lag_Matrix(Lag_Index) = SSD_1D2(Size_Diff, Target_1D, Search_1D, Point, &
                                                      Search_Index, Flag)

  ! Bottom point

  Lag_Index = Min_SSD_Location + Lag_Element_Size
  Search_Index = Min_SSD_Line * Small_Search_Element_Size + Min_SSD_Element
  Lag_Matrix(Lag_Index) = SSD_1D2(Size_Diff, Target_1D, Search_1D, Point, &
                                                      Search_Index, Flag)

END SUBROUTINE Fill_Lag_Matrix_1D

!-------------------------------------------------------------------------------
!
! Name:
!   SSD_1D2
!
! Function:
!   Coputes the sum of squared difference value at a single location of
!   the lag matrix.
!
! Description:
!   This function computes the SSD value at a single location of the lag
!   matrix, the array of SSD values of all potential matching scenes in
!   the search region.
!
! Calling sequence:
!   SSD_1D2
!
! Inputs:
!   Size_Diff            Size difference in pixels between the nested box and
!                        the outer target box.
!   Target_1D            One-dimensional array of target box values
!   Search_1D            One-dimensional array of search box values
!   Point                Location in lag matrix
!   Search_Index         Starting pixel in search region
!   Flag                 Flag for terminating the SSD calculation if the
!                        current minimum has been exceeded
!
! Outputs/Result:
!   Sum_Of_Squared_Differences
!
! Dependencies:
!   Assumes SSD similarity measure is being used for tracking.
!
! Restrictions:
!   None
!
! History:
!   06/01/2009 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
FUNCTION SSD_1D2(Size_Diff, Target_1D, Search_1D, Point, Search_Index, Flag) &
                                         RESULT (Sum_Of_Squared_Differences)

  INTEGER(BYTE), INTENT(IN) :: Flag
  INTEGER(LONG), INTENT(IN) :: Point
  INTEGER(LONG), INTENT(IN) :: Search_Index
  INTEGER(LONG) :: Search_Index2
  INTEGER(LONG) :: Target_Index
  INTEGER(LONG) :: Element
  INTEGER(LONG) :: Line
  INTEGER(LONG) :: Element_Counter_Target
  INTEGER(LONG) :: Element_Counter_Search
  INTEGER(LONG) :: Small_Target_Element_Size
  INTEGER(LONG) :: Small_Target_Line_Size
  INTEGER(LONG) :: Small_Search_Element_Size
  INTEGER(LONG), INTENT(IN) :: Size_Diff
  REAL(SINGLE), DIMENSION(:), INTENT(IN) :: Target_1D
  REAL(SINGLE), DIMENSION(:), INTENT(IN) :: Search_1D
  REAL(SINGLE) :: Sum_Of_Squared_Differences
  REAL(SINGLE), SAVE :: Current_Minimum

  Sum_Of_Squared_Differences = 0.0
  Small_Target_Element_Size = 5
  Small_Target_Line_Size = 5

  ! If size difference is zero this indicates that nested tracking is not being
  ! used so we need to set the small target values to the larger size

  IF (Size_Diff .EQ. 0) THEN
    Small_Target_Element_Size = Target_Element_Size
    Small_Target_Line_Size = Target_Line_Size
  ENDIF

  Small_Search_Element_Size = Search_Element_Size - Size_Diff
  Element_Counter_Target = - Small_Target_Element_Size
  Element_Counter_Search = Search_Index - Small_Search_Element_Size - 1

  IF (Point .EQ. 1) Current_Minimum = HUGE(REAL(SHORT))

  line_loop: DO Line = 1, Small_Target_Line_Size

    Element_Counter_Target = Element_Counter_Target + Small_Target_Element_Size
    Element_Counter_Search = Element_Counter_Search + Small_Search_Element_Size

    element_loop: DO Element = 1, Small_Target_Element_Size

      Target_Index = Element_Counter_Target + Element
      Search_Index2 = Element_Counter_Search + Element

      ! perform this check to avoid a zero index
      ! wcb 2/28/12

      IF (Search_Index2 .EQ. 0) Search_Index2 = 1

      Sum_Of_Squared_Differences = Sum_Of_Squared_Differences + &
         (Target_1D(Target_Index) - Search_1D(Search_Index2))**2

      ! ------------------------------------------------------------------------
      ! Terminate summation if current value exceeds minimum
      ! ------------------------------------------------------------------------

      IF (Sum_Of_Squared_Differences .GT. Current_Minimum &
                                .AND. Flag .EQ. 0) RETURN

    END DO element_loop

  END DO line_loop

  ! ----------------------------------------------------------------------------
  ! Reset current minimum for new minimum
  ! ----------------------------------------------------------------------------

  IF (Sum_Of_Squared_Differences .LT. Current_Minimum) Current_Minimum = &
                                              Sum_Of_Squared_Differences

END FUNCTION SSD_1D2

! ------------------------------------------------------------------------------
!
! Name:
!   Sort_Shell_NP
!
! Function:
!  Sort array of values into ascending order and find median. Differs from
!  original version in that it accepts as one of the arguments the number
!  of points to sort. This is useful if the user does not wish to sort the
!  entire array which may have missing elements.
!
! Description:
!   - Uses Shell's method
!
! Reference:
!   Numerical Recipes
!
! Calling sequence:
!   CALL Sort_Shell
!
! Inputs:
!   Number_Of_Points     Number of points to be sorted.
!   Array                1-D array of values to be sorted.
!
! Outputs:
!   Median               Median of array values.
!
! Dependencies:
!
! Restrictions:
!
! History:
!   8/2008 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE Sort_Shell_NP(Number_Of_Points, Array, Median)

  INTEGER(LONG) :: I
  INTEGER(LONG) :: J
  INTEGER(LONG) :: Increment
  INTEGER(LONG), INTENT(IN) :: Number_Of_Points
  INTEGER(LONG) :: Mid_Point
  REAL(SINGLE) :: Value
  REAL(SINGLE), DIMENSION(:), INTENT(INOUT) :: Array
  REAL(SINGLE), INTENT(OUT) :: Median
  REAL(SINGLE) Mean

  Increment = 1

  IF (Number_Of_Points .EQ. 1) THEN

    Median = Array(1)
    RETURN

  ENDIF

  DO
    Increment = 3 * Increment + 1
    IF (Increment > Number_Of_Points) EXIT
  END DO

  DO
    Increment = Increment / 3
    DO I = Increment + 1, Number_Of_Points
      Value = Array(I)
      J = I
      DO
        IF (Array(J-Increment) <= Value) EXIT
        Array(J) = Array(J-Increment)
        J = J - Increment
        IF (J <= Increment) EXIT
      END DO

      Array(J) = Value

    END DO

    IF (Increment <= 1) EXIT

  END DO


  Mid_Point = Number_Of_Points / 2

  IF (MOD(Number_Of_Points, 2) .EQ. 0) THEN
    Median = 0.5*(Array(Mid_Point) + Array(Mid_Point + 1))
  ELSE
    Median = Array(Mid_Point + 1)
  ENDIF

! Test computing average instead of median

  Mean = SUM(Array, MASK = Array > 0.0) / REAL(COUNT(Array > 0.0))

END SUBROUTINE Sort_Shell_NP

!-------------------------------------------------------------------------------
!
! Name:
!   DBSCAN
!
! Function:
!   Identifies clusters in the presence of noise.
!
! Description:
!   An efficient cluster analysis algorithm that separates areas of higher
!   density from areas of lower density. The algorithm is very effective at
!   discriminating between coherent clusters and noise.
!
! References:
!   Ester et al (1996), "A Density-Based Algorithm for Discovering Clusters in
!   Large Spatial Databases with Noise." Proceedings of 2nd International
!   Conference on Knowledge Discovery and Data Mining (KDD-96), Portland,
!   Oregon, USA, 226-231.
!
! Calling sequence:
!   CALL DBSCAN
!
! Inputs:
!   GoodSample           Number of points in initial sample
!   SetOfPoints          Set of x and y displacements
!   MinPts               Minimum number of points for cluster
!   Eps                  Radius around point to search for cluster
!   Median_Press         Median pressure of cold sample
!
! Outputs:
!   NumberOfClusters     Number of clusters identified
!   MaxCluster           Cluster ID of largest cluster
!   MaxClusterSize       Number of points in largest cluster
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!   11/2009 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE DBSCAN(GoodSample, SetOfPoints, MinPts, Eps, &
          NumberOfClusters, MaxCluster, MaxClusterSize)

  INTEGER(LONG), INTENT(IN) :: GoodSample
  INTEGER(LONG), INTENT(IN) :: MinPts
  INTEGER(LONG) ClusterID
  INTEGER(LONG) Point
  INTEGER(LONG) ClusterSize
  INTEGER(LONG), INTENT(OUT) :: NumberOfClusters
  INTEGER(LONG), INTENT(OUT) :: MaxClusterSize
  INTEGER(LONG), INTENT(OUT) :: MaxCluster
  INTEGER(LONG) Dim1
  INTEGER(LONG) Dim2
  INTEGER(LONG) I
  INTEGER(LONG) J
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: SeedList
  INTEGER(LONG), ALLOCATABLE, DIMENSION(:) :: ResultList

  REAL(SINGLE), INTENT(IN) :: Eps
  REAL(SINGLE), DIMENSION(:,:), INTENT(INOUT) :: SetOfPoints

! Initialize variables

  ClusterID = 0
  Dim1 = SIZE(SetOfPoints, 1)
  Dim2 = SIZE(SetOfPoints, 2)
  SetOfPoints(1:Dim1,3) = UNCLASSIFIED
  SetOfPoints(1:Dim1,4) = MISSING_VALUE_REAL4
  NumberOfClusters = 0
  MaxClusterSize = 0
  MaxCluster = 0

  ALLOCATE(SeedList(Dim1))
  ALLOCATE(ResultList(Dim1))


  DO I = 1,GoodSample

    Point = I

    ! Check if point has been classified already. If it has not already been
    ! classsified, then check to see if it is part of a cluster.

    IF (SetOfPoints(Point, 3) .EQ. UNCLASSIFIED) THEN

      IF (ExpandCluster(GoodSample, SetOfPoints, Point, &
          SeedList, ResultList, ClusterID, MinPts, Eps)) THEN

        ! If the point is part of a cluster, increment the ID counter

        ClusterID = ClusterID + 1

      ENDIF

    ENDIF

  ENDDO

  ! All points have been classified at this stage.

  MaxClusterSize = 0
  ClusterSize = 0

  ! If at least one cluster exists, go back and count the number of points
  ! in each cluster.

  IF(ClusterID .GT. 0) THEN

    DO J = 1, ClusterID

      DO I = 1, GoodSample

        IF (SetOfPoints(I,4) .EQ. REAL(J-1)) THEN

          ClusterSize = ClusterSize + 1

        ENDIF

      ENDDO

      IF (ClusterSize .GT. MaxClusterSize) THEN

        MaxClusterSize = ClusterSize
        MaxCluster = J - 1

      ENDIF

      ! reset counter for next cluster

      ClusterSize = 0

    ENDDO

  ENDIF

   ! Return the number of clusters found

   NumberOfClusters = ClusterID

  DEALLOCATE(SeedList, ResultList)

END SUBROUTINE DBSCAN

!-------------------------------------------------------------------------------
!
! Name:
!   ExpandCluster
!
! Function:
!   Function used to classify a given point as either noise or part of a cluster.
!
! Description:
!   Most important function in DBSCAN algorithm. It determines if the current
!   point is part of a cluster or noise.
!
! References:
!   Ester et al (1996), "A Density-Based Algorithm for Discovering Clusters in
!   Large Spatial Databases with Noise." Proceedings of 2nd International
!   Conference on Knowledge Discovery and Data Mining (KDD-96), Portland,
!   Oregon, USA, 226-231.
!
! Calling sequence:
!   ExpandCluster
!
! Inputs:
!   SampleSize           Number of points in initial sample
!   Point                Point being classified
!   Median_Press         Median pressure of cold sample
!   ClusterID            ID of current cluster
!   MinPts               Minimum number of points for cluster
!   Eps                  Radius around point to search for cluster
!
! Input/Output
!   SetOfPoints          Set of x and y displacements
!   SeedList             List of points within 1-Eps radius of current point
!   ResultList           List of points within 1-Eps radius of a seed point
!
! Outputs:
!   TRUE or FALSE        TRUE if point is located within cluster, FALSE if not
!
! Note:   A point can have one of 3 status values:
!
!         UNCLASSIFIED - point has not been classified yet
!         NOISE - point that doesn't belong to any cluster
!         CLASSIFIED - point has been classified previously
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!   11/2009 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
LOGICAL FUNCTION ExpandCluster(SampleSize, SetOfPoints, Point, SeedList, &
                                         ResultList, ClusterID, MinPts, Eps)

  INTEGER(LONG), INTENT(IN) :: SampleSize
  INTEGER(LONG), INTENT(IN) :: MinPts
  INTEGER(LONG), INTENT(IN) :: ClusterID
  INTEGER(LONG) NumberOfSeeds
  INTEGER(LONG), INTENT(IN) ::  Point
  INTEGER(LONG) CurrentP
  INTEGER(LONG) ResultP
  INTEGER(LONG) NumberOfResults
  INTEGER(LONG) I
  INTEGER(LONG) J
  INTEGER(LONG), INTENT(INOUT), DIMENSION(:) :: SeedList
  INTEGER(LONG), INTENT(INOUT), DIMENSION(:) :: ResultList

  REAL(SINGLE), INTENT(IN) :: Eps

  ! set of points being analyzed - these are actually displacements in
  ! line/element space
  REAL(SINGLE), DIMENSION(:,:), INTENT(INOUT) :: SetOfPoints

  ! set default value of function to false
  ExpandCluster = .FALSE.

  ! get list of points within 1 epsilon (Eps) of point under consideration

  CALL GetSeedList(SampleSize, SetOfPoints, Point, Eps, &
                               NumberOfSeeds, SeedList)

  ! if the number of points in the neighborhood is less than the minimum
  ! number for a cluster, then classify the point as noise

  IF (NumberOfSeeds .LT. MinPts) THEN

    SetOfPoints(Point, 3) = NOISE
    RETURN

  ELSE

    ! core point found

    SetOfPoints(Point, 3) = CLASSIFIED
    SetOfPoints(Point, 4) = ClusterID

    ! add the list of points to the current cluster and
    ! set the status of each to classified

    DO I = 1,NumberOfSeeds

      SetOfPoints(SeedList(I), 3) = CLASSIFIED
      SetOfPoints(SeedList(I), 4) = ClusterID

    ENDDO

    ! now check the neighborhood of each point in the seed list
    ! to see if it is "connected" to the current cluster

    DO WHILE (NumberOfSeeds .GT. 0)

      CurrentP = SeedList(1)

      CALL GetSeedList(SampleSize, SetOfPoints, CurrentP, Eps, &
                                  NumberOfResults, ResultList)

      IF (NumberOfResults .GE. MinPts) THEN

        DO J = 1,NumberOfResults

          ResultP = ResultList(J)

          IF (SetOfPoints(ResultP, 3) .EQ. UNCLASSIFIED .OR. &
                    SetOfPoints(ResultP, 3) .EQ. NOISE) THEN

            IF (SetOfPoints(ResultP, 3) .EQ. UNCLASSIFIED) THEN

!             Add point to list of seeds
              CALL UpdateSeedList(APPEND_PT, ResultP, NumberOfSeeds, SeedList)

            ENDIF

            SetOfPoints(ResultP, 3) = CLASSIFIED
            SetOfPoints(ResultP, 4) = ClusterID

          ENDIF

        ENDDO

      ENDIF

!     Delete point from list of seeds
      CALL UpdateSeedList(DELETE_PT, CurrentP, NumberOfSeeds, SeedList)

    ENDDO

    ExpandCluster = .TRUE.

  ENDIF

END FUNCTION ExpandCluster

!-------------------------------------------------------------------------------
!
! Name:
!   GetSeedList
!
! Function:
!   Returns a list of points within 1 Eps radius of current point.
!
! Description:
!   A routine for returning a list of points within 1 Eps of a given point.
!   The entire sample of points is searched.
!
!
! References:
!   Ester et al (1996), "A Density-Based Algorithm for Discovering Clusters in
!   Large Spatial Databases with Noise." Proceedings of 2nd International
!   Conference on Knowledge Discovery and Data Mining (KDD-96), Portland,
!   Oregon, USA, 226-231.
!
! Calling sequence:
!   CALL GetSeedList
!
! Inputs:
!   SampleSize           Number of points in initial sample
!   Points               Original sample of points
!   Point                Point being searched
!   Eps                  Radius around point to search for neighbors
!   Median_Press         Median pressure of cold sample
!
! Outputs:
!   NumberOfSeeds        Number of points within 1-Eps radius of current point
!   SeedList             List of points within 1-Eps radius of current point
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!   11/2009 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE GetSeedList(SampleSize, Points, Point, Eps, NumberOfSeeds, SeedList)

  INTEGER(LONG), INTENT(IN) :: SampleSize
  INTEGER(LONG), INTENT(IN) :: Point
  INTEGER(LONG), INTENT(OUT) :: NumberOfSeeds
  INTEGER(LONG) I
  INTEGER(LONG), INTENT(INOUT), DIMENSION(:) :: SeedList

  REAL(SINGLE), INTENT(IN) :: Eps
  REAL(SINGLE) Distance
  REAL(SINGLE), DIMENSION(:,:), INTENT(IN) :: Points
  REAL(SINGLE) X1
  REAL(SINGLE) Y1
  REAL(SINGLE) X2
  REAL(SINGLE) Y2
  ! test adding third dimension
  REAL(SINGLE) Z1
  REAL(SINGLE) Z2
  REAL(SINGLE) ZDIFF

  ! Find all points within 1 epsilon of point under consideration.

  NumberOfSeeds = 0
  X1 = Points(Point,1)
  Y1 = Points(Point,2)
  ! add vertical dimension (height in km)
  ! Division by 1000 converts units from meters to kilometers
  Z1 = Points(Point,8)/KM2M

  DO I = 1,SampleSize

    IF (I .NE. Point) THEN

      X2 = Points(I,1)
      Y2 = Points(I,2)
      ! add vertical dimension (height in km)
      ! Division by 1000 converts units from meters to kilometers
      Z2 = Points(I,8)/KM2M

      ZDIFF = Z2-Z1

      Distance = SQRT((X2-X1)**2 + (Y2-Y1)**2 )

      IF (Distance .LE. Eps) THEN

        NumberOfSeeds = NumberOfSeeds + 1
        SeedList(NumberOfSeeds) = I

      ENDIF

    ENDIF

  ENDDO

END SUBROUTINE GetSeedList

!-------------------------------------------------------------------------------
!
! Name:
!   UpdateSeedList
!
! Function:
!   Update the list of seed points.
!
! Description:
!   This routine adds or removes a point from the list of seeds.
!
! References:
!   Ester et al (1996), "A Density-Based Algorithm for Discovering Clusters in
!   Large Spatial Databases with Noise." Proceedings of 2nd International
!   Conference on Knowledge Discovery and Data Mining (KDD-96), Portland,
!   Oregon, USA, 226-231.
!
! Calling sequence:
!   CALL UpdateSeedList
!
! Inputs:
!   Operation            The operation to perform on the list
!   CurrentP             Point to operate on
!   NumberOfSeeds        Number of points in seed list
!
! Input/Output:
!   SeedList             List of seed points to be updated
!
! Dependencies:
!   None
!
! Restrictions:
!   None
!
! History:
!   11/2009 - Wayne Bresky - Created
!
!-------------------------------------------------------------------------------
SUBROUTINE UpdateSeedList(Operation, CurrentP, NumberOfSeeds, SeedList)

  INTEGER(LONG), INTENT (INOUT) :: NumberOfSeeds
  INTEGER(LONG), INTENT(IN) :: CurrentP
  INTEGER(LONG) I
  INTEGER(BYTE), INTENT(IN) :: Operation
  INTEGER(LONG), INTENT(INOUT), DIMENSION(:) :: SeedList

  IF (Operation .EQ. DELETE_PT) THEN

    NumberOfSeeds = NumberOfSeeds - 1

    IF (NumberOfSeeds .GT. 0) THEN

      DO I = 1,NumberOfSeeds
         SeedList(I) = SeedList(I+1)
      ENDDO

    ENDIF

  ELSE IF(Operation .EQ. APPEND_PT) THEN

    NumberOfSeeds = NumberOfSeeds + 1
    SeedList(NumberOfSeeds) = CurrentP

  ENDIF

END SUBROUTINE UpdateSeedList

END MODULE AMV_EN_FEATURE_TRACKING_UTILS_M
