! ------------------------------------------------------------------------------
!
! Name:
!   Target_Selection
!
! Function:
!   Select targets to be used for deriving atmospheric motion vectors (AMVs).
!
! Description:
!   The purpose of this module is to select image features that will
!   subsequently be tracked to derive atmospheric motion. The image data is
!   provided by the main GEOCAT program in the form of line segments. The
!   line segment data is subdivided into smaller regions (typically 15 lines
!   by 15 elements) called target scenes and a variety of tests are performed
!   on each scene. 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
!  Target QC    bad data (reasonable data values)
!               spatial coherence
!               multi-layer check
!               day/night terminator check (visible and 3.9 channel only)
!               valid height assignment (at least 2 percent of pixels)
!               pressure thresholds (band specific)
!   ----------------------------------------------------------------------------
!                max correlation on boundary of search
!                minimum correlation check (current threshold is 0.6)
!   Tracking QC  acceleration check + gross error checks
!                EUMETSAT Quality Indicator (QI)
!                Expected Error (EE)
!
!   Additonal checks have been introduced with the advent of nested tracking,
!   these include:
!
!               zero sample of local motion vectors
!               change in median pressure between sample 1 (reverse vector)
!               and sample 2 (forward vector) exceeds 100 mb
!
!   Notes:
!
!   1. This routine uses a buffer (nominally 105 lines) to hold image data
!   (brightness temperatures, reflectances and radiances) as well as output
!   from the cloud mask and cloud height algorithms, which must run prior to
!   this routine being called. No AMVs are generated until the data buffer
!   is full.
!
!   2. Only the middle portion of the buffer is processed. This is to allow
!   cloud features to move outside the strip being processed for targets.
!
!
!   3. The target box is positioned over the maximum BrtTemp gradient within
!   the original target scene. This requires that the box data be re-sampled
!   because the initial target box may or may not be centered on the maximum
!   gradient. The above tests are applied to the scene centered on the max
!   gradient.
!
!   4. The spatial coherence analysis is used to filter out scenes which are too
!   coherent (and thus difficult to track by correlation methods) or that
!   contain more than two layers (height assignment of multi-layered target
!   scenes is not attempetd).
!
! 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).
!
!   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:
!   Use GEOCAT command line with cloud mask, cloud type, cloud height and the
!   AMV algorithms all specified.
!
!   (e.g., geocat -akey baseline_cmask_seviri baseline_ctype_seviri
!                       baseline_cld_hght_seviri winds ....)
!
! Inputs:
!   All input data is passed through GEOCAT structures (satellite data, cloud
!   mask, space mask, cloud top pressure). The file "geocat.default" should be
!   used to specify default configuration values. Most importantly, the number
!   of scan lines processed at one time should be set (usually 15 lines).
!
!   The include file "winds.inc" is used to specify certain parameters specific
!   to the AMV algorithm.
!
! Outputs:
!   1. ASCII file - to be replaced with netCDF file.
!
! Dependencies:
!   GEOCAT satellite data structure must be populated for this segment as well
!   as the cloud height error from the out2 data structure. The configutration
!   file "winds.config" controls several processing configuration settings
!   including the size of the target box, the satellite band to process, the
!   loop time interval and whether or not to attempt nested tracking. The
!   include file "winds.inc" is also required.
!
!   The cloud mask, cloud type and cloud height algorithms must run prior
!   to calling this routine.
!
! Restrictions:
!
! History:
!  04/2007 - Wayne Bresky - Created.
!  11/2008 - Wayne Bresky - New bands added. Removed constraint that line
!                           segment had to be multiple of target box size.
!  02/2009 - Wayne Bresky - Gross error checks introduced with input from
!                           Chris Velden and Steve Wanzong at CIMSS.
!  03/2009 - Wayne Bresky - Low level inversion flag, output from Cloud
!                           height algorithm, added to buffer.
!  03/2010 - Wayne Bresky - Nested tracking introduced.
!  04/2010 - Wayne Bresky - Code Unit Test Review (CUTR) issues addressed.
!  08/2010 - Wayne Bresky / Steve Wanzong - Weighting function-derived heights
!                           added for clear-sky WV winds, additional diagnostic
!                           output variables added.
!
!-------------------------------------------------------------------------------
MODULE AMV_EN_TARGET_SELECTION_M
      USE ds_data_mod
      USE cf_data_mod
      USE fw_log_mod
      USE fw_scene_mod
      USE fw_str_util_mod
      USE URIHelper_mod

      USE AMV_EN_input_m
      USE AMV_EN_Quality_m
      USE AMV_EN_BestFitPress_m
      USE AMV_EN_Wind_Data_m
      USE AMV_EN_Error_Counts_m
      USE AMV_EN_XPATH_M
      USE AMV_EN_TARGET_SELECTION_UTILS_M
      USE AMV_EN_FEATURE_TRACKING_UTILS_M
      USE AMV_EN_WINDS_INC
   IMPLICIT NONE

   PRIVATE
   PUBLIC :: Target_Selection_Main


CONTAINS

   FUNCTION get_instrument_parameters(Nav, SatName, IntrumentName, resolution, ChannelName, &
                                      ImageStartElem_native, ImageStartLine_native, SatLat, SatLon) RESULT(ret)
      CHARACTER(len=*), PARAMETER :: FUNC = "get_instrument_parameters"
      TYPE(Nav_t), INTENT(OUT) :: Nav !< contains everything needed for pixcoord2geocoord and geocoord2pixcoord.
      CHARACTER(LEN=:), ALLOCATABLE, INTENT(OUT) :: SatName !< GOES16, GOES17, NPP, NOAA20, etc.
      CHARACTER(LEN=:), ALLOCATABLE, INTENT(OUT) :: IntrumentName !< ABI, VIIRS, etc.
      CHARACTER(LEN=:), ALLOCATABLE, INTENT(OUT) :: resolution !< ABI_2KM, ABI_500M, VIIRS_750, etc.
      CHARACTER(LEN=:), ALLOCATABLE, INTENT(OUT) :: ChannelName !< AMV Channel
      type(fw_data_context) :: ctx_line  !< context. Describes location within the larger picture
      type(fw_data_context):: ctx_elem !< context. Describes location within the larger picture
      INTEGER(DLONG), INTENT(OUT) :: ImageStartElem_native !< segment start x (element)
                                                           !< in the native image coordinates (0 or 1 at origin)
      INTEGER(DLONG), INTENT(OUT) :: ImageStartLine_native !< segment start y (line)
                                                           !< in the native image coordinates (0 or 1 at origin)
      REAL(SINGLE), INTENT(OUT)  :: SatLat !< satellite latitude subpoint used in computing wind azimuth/zenith angle
      REAL(SINGLE), INTENT(OUT)  :: SatLon !< satellite longitude subpoint used in computing wind azimuth/zenith angle
      INTEGER(LONG):: Sat_ID

      LOGICAL :: ret

      !TYPE(Nav_t) :: Nav
      INTEGER(SHORT), DIMENSION(:), POINTER :: X !< indicates the navigation x/y. used to determine first x/y (0 or 1 origin)
      INTEGER(SHORT), DIMENSION(:), POINTER :: Y !< indicates the navigation x/y. used to determine first x/y (0 or 1 origin)
      LOGICAL :: has_INPUT_X !< certain sensors might not store x and y
      type(fw_data_id) :: radiometricID !< describes the input BT/Radiance
      ret = .FALSE.



      NULLIFY(X,Y)


      !
      ! Get a descriptor for the incoming BT or Radiance.
      !
      IF(.NOT. cf_get(INPUT_MIDDLE_RADIOMETRIC_DATA_LBL, radiometricID, FUNC)) RETURN
      SatName = radiometricID%satellite_id
      IntrumentName = radiometricID%instrument_id
      resolution = radiometricID%resolution_id
      ChannelName = radiometricID%channel_id
      IF(.NOT. ds_get(INPUT_SENSOR_SERIES_LBL, Sensor_Series, FUNC)) RETURN


      IF(.NOT. cf_get(PAR_NAV_TYPE_LBL, Nav%NavType, FUNC)) RETURN
      IF(.NOT. ds_get(INPUT_PROJECTED_LAT_LBL, SatLat, FUNC)) RETURN ! inaccurate to use projection origin
                                                                     ! we should injest both projection and nadir.
      IF(.NOT. ds_get(INPUT_PROJECTED_LON_LBL, SatLon, FUNC)) RETURN ! inaccurate to use projection origin
                                                                     ! we should injest both projection and nadir.
      IF(.NOT. ds_get(INPUT_MCIDAS_ID_LBL, Sat_ID, FUNC)) RETURN
      Nav%Sat_ID = Sat_ID
      Nav%Sensor_Series = Sensor_Series ! this is a global right now.




      IF(IntrumentName == "SEVIRI") THEN
         IF(.NOT.ds_get(INPUT_ORBITAL_RADIUS_LBL, Nav%SEVIRI%SAT_HEIGHT_MSG, FUNC)) RETURN
         IF(.NOT.ds_get(INPUT_COLUMN_FACTOR_LBL, Nav%SEVIRI%CFAC, FUNC)) RETURN
         IF(.NOT.ds_get(INPUT_LINE_FACTOR_LBL, Nav%SEVIRI%LFAC, FUNC)) RETURN
         IF(.NOT.ds_get(INPUT_COLUMN_OFFSET_LBL, Nav%SEVIRI%COFF, FUNC)) RETURN
         IF(.NOT.ds_get(INPUT_LINE_OFFSET_LBL, Nav%SEVIRI%LOFF, FUNC)) RETURN
         Nav%SEVIRI%SubSatLonRad = REAL(SatLon/DTOR, DOUBLE)
      ELSEIF(IntrumentName == "IMAGER") THEN
         !BH CALL Build_GVAR_Nav_Struct(Ctxt)
         CALL fw_log_error(FUNC, "Not yet supported.")
         RETURN
      ELSEIF (TRIM(Nav%NavType) .EQ. 'PS' ) THEN
         IF(.NOT.FillNav_Polar(Nav)) RETURN
      ELSE
         Nav%ABI%SubSatLon = SatLon
         Nav%ABI%SubSatLat = SatLat
         IF(.NOT.ds_get(INPUT_X_SCALE_LBL, Nav%ABI%XScale, FUNC)) RETURN
         IF(.NOT.ds_get(INPUT_Y_SCALE_LBL, Nav%ABI%YScale, FUNC)) RETURN
         IF(.NOT.ds_get(INPUT_X_OFFSET_LBL, Nav%ABI%XOffset, FUNC)) RETURN
         IF(.NOT.ds_get(INPUT_Y_OFFSET_LBL, Nav%ABI%YOffset, FUNC)) RETURN
      ENDIF


      IF(.NOT. fw_uri_segment_ctx(INPUT_MIDDLE_RADIOMETRIC_DATA_LBL, FW_DIM_ROW, ctx_line, FUNC)) RETURN
      IF(.NOT. fw_uri_segment_ctx(INPUT_MIDDLE_RADIOMETRIC_DATA_LBL, FW_DIM_COLUMN, ctx_elem, FUNC)) RETURN
      IF (.NOT.cf_has(INPUT_X_LBL, has_INPUT_X, FUNC)) RETURN
      IF( has_INPUT_X ) THEN
         IF(.NOT.ds_get(INPUT_X_LBL, X, FUNC)) RETURN
         IF(.NOT.ds_get(INPUT_Y_LBL, Y, FUNC)) RETURN
         ImageStartElem_native = INT(X(1), KIND=DLONG)
         ImageStartLine_native = INT(Y(1), KIND=DLONG)
      ELSE
         ImageStartElem_native = ctx_elem%position%start_pos
         ImageStartLine_native = ctx_line%position%start_pos
      ENDIF


      ret = .TRUE.
   END FUNCTION get_instrument_parameters


   !> this is the main subroutine to perform wind retrieval
   SUBROUTINE Target_Selection_Main(Stat)


      IMPLICIT NONE
      CHARACTER(LEN=*), PARAMETER :: FUNC = "Target_Selection_NOP_Main"

      TYPE(Nav_t) :: Nav !< stores all the needed navigation parameters for any given sensor/projection
      INTEGER(DLONG) :: ImageStartElem_native !< segment start x (element) in the native image coordinates (0 or 1 at origin)
      INTEGER(DLONG) :: ImageStartLine_native !< segment start y (line) in the native image coordinates (0 or 1 at origin)
      CHARACTER(LEN=:), ALLOCATABLE :: SatName !< GOES16, GOES17, METOPC, etc.
      CHARACTER(LEN=:), ALLOCATABLE :: IntrumentName !< ABI, SEVIRI, AVHRR, etc.
      CHARACTER(LEN=:), ALLOCATABLE :: resolution !< ABI_2KM, ABI_500M, etc.
      CHARACTER(LEN=:), ALLOCATABLE :: ChannelName !< The AMV Channel, before potential degrading
      INTEGER(LONG) :: Stat

      !input descriptions
      CHARACTER(len=128) :: string
      TYPE(fw_time_stamp) :: middle_time_stamp !< nominal timestamp of the middle image
      INTEGER(LONG), DIMENSION(NUMBER_TIMESTEPS) :: EpochImageTimes !< granule times of 3 images [seconds since 1/1/1970]
      REAL(SINGLE)  :: Time_Interval  !< average time between images [min]
      REAL(SINGLE)  :: Time_Interval_minus !< time between middle and first image [min]
      REAL(SINGLE)  :: Time_Interval_plus !< time between middle and last image [min]


      !Navigation inputs
      REAL(SINGLE) :: SatLat  !< satellite latitude subpoint used in computing wind azimuth/zenith angle
      REAL(SINGLE) :: SatLon  !< satellite longitude subpoint used in computing wind azimuth/zenith angle
      type(fw_data_context):: ctx_line !< context. Describes location within the larger picture
      type(fw_data_context):: ctx_elem !< context. Describes location within the larger picture
      type(fw_data_range):: lineROI !< the region of interest within the padded segment
      type(fw_data_range):: elemROI !< the region of interest within the padded segment


      !
       !inputs
      !
      TYPE(input_t) :: input !< stores all the inputs
      INTEGER(BYTE), DIMENSION(:,:), POINTER :: Space_Mask !< space mask of the middle granule
      TYPE(ForecastProfilesGrid_t) :: forecastProfiles !< gfs grid/profiles for the middle granule

      INTEGER(SHORT) :: BUFFER_SEGMENTS !< the original winds implementation had this.
                                        !!  Tracking starts at the center line of a buffer.
                                        !!  The buffer is roughly BUFFER_SEGMENTS times the nominal box size.
                                        !!  The buffer is then incremented towards the bottom of the image
                                        !!  until there is no more padding remaining in the segment.
                                        !!  All it's doing now is preventing tracking from starting
                                        !!  at the beginning and ending of an image
                                        !!  and that isn't especially desirable anyway.

      CHARACTER(len=:), ALLOCATABLE :: targetType  !< Cloud Top ('CT') or Clear Sky ('CS')
      INTEGER(LONG) :: Lag_Size !< Lag Size used in the search
      INTEGER(LONG) :: Nested_Tracking_Flag !< flag indicating nested tracking is enabled (sym%YES or sym%NO)

      ! aab Tandem mod
      CHARACTER(len=:), ALLOCATABLE :: TandemTriplet  !< Tandem-VPW Triplet Configuration 

      !
      ! Outputs
      !
      INTEGER(LONG), DIMENSION(41) :: Error_Counts(-10:30)
      REAL(SINGLE), DIMENSION(:, :), ALLOCATABLE :: Wind_Data !< per target data storage
                                                              !! where the second index is used to store different variables.
      TYPE(TargetData_t) :: TargetData !< per target data storage. can replace Wind_data.
      INTEGER(LONG) :: Target_Number !< number of targets
      INTEGER(LONG) :: Good_Targets !< number of targets with no errors

      INTEGER(LONG) :: MAX_WIND_RECORDS !< must be above the maximum number of targets otherwise wind_data runs out of room to store results.
      REAL(SINGLE)   :: CORR_MIN  !< correlation must be greater than this value.


      LOGICAL :: isTrackingSkipped !<  tracking can be skipped when inputs are insufficient quality

      Stat = RETURN_FAIL

      Target_Number = 0
      Good_Targets = 0
      Good_Winds = 0
      Error_Counts(-10:30) = 0


      IF(.NOT. cf_get(PAR_AMV_IDX_LBL, Channel, FUNC)) RETURN

      IF(.NOT. get_instrument_parameters(Nav, SatName, IntrumentName, resolution, ChannelName, &
                                         ImageStartElem_native, ImageStartLine_native, SatLat, SatLon)) RETURN

      !
      ! Obtain the position and region of interest within the BT/Radiance image.
      ! The image has padding and the region of interest is the part we wish to keep.
      IF(.NOT. fw_uri_segment_ctx(INPUT_MIDDLE_RADIOMETRIC_DATA_LBL, FW_DIM_ROW, ctx_line, FUNC)) RETURN
      IF(.NOT. fw_uri_segment_ctx(INPUT_MIDDLE_RADIOMETRIC_DATA_LBL, FW_DIM_COLUMN, ctx_elem, FUNC)) RETURN
      lineROI = fw_roi_in_ctx(ctx_line)
      elemROI = fw_roi_in_ctx(ctx_elem)


      IF(.NOT. ds_get(INPUT_SPACE_MSK_LBL, Space_Mask, FUNC)) RETURN


      IF(.NOT. ds_get(INPUT_MIDDLE_TIME_STAMP_LBL, middle_time_stamp, FUNC)) RETURN



      IF(.NOT. cf_get(PAR_BOX_SIZE_LBL, Box_Size, FUNC)) RETURN
      IF(.NOT. cf_get(PAR_LAG_SIZE_LBL, Lag_Size, FUNC)) RETURN !< this is a function of time delta
      IF(.NOT. cf_get(PAR_NESTED_TRACKING_FLAG_LBL, Nested_Tracking_Flag, FUNC)) RETURN
      IF(.NOT. cf_get(PAR_TARGET_TYPE_LBL, targetType, FUNC)) RETURN

      ! aab Tandem mod
      IF(.NOT. cf_get(PAR_TANDEM_TRIPLET_LBL, TandemTriplet, FUNC)) RETURN

      IF(.NOT. cf_get(PAR_MAX_WINDS_PER_SEGMENT_LBL, MAX_WIND_RECORDS, FUNC)) RETURN


      IF(.NOT. ds_get(FIRST_GMT_EPOCH1970_LBL, EpochImageTimes(FIRST), FUNC)) RETURN
      IF(.NOT. ds_get(MIDDLE_GMT_EPOCH1970_LBL,EpochImageTimes(MIDDLE), FUNC)) RETURN
      IF(.NOT. ds_get(LAST_GMT_EPOCH1970_LBL, EpochImageTimes(LAST), FUNC)) RETURN

      IF(.NOT. ds_get(TIME_INTERVAL_MINUS_LBL, Time_Interval_minus, FUNC)) RETURN
      IF(.NOT. ds_get(TIME_INTERVAL_PLUS_LBL, Time_Interval_plus, FUNC)) RETURN
      IF(.NOT. ds_get(TIME_INTERVAL_AVERAGE_LBL, Time_Interval,FUNC)) RETURN
      !---------------------------------------
      ! set parameters according to instrument
      !---------------------------------------
      IF (TRIM(Nav%NavType) .EQ. 'PS' ) THEN
         BUFFER_SEGMENTS = BUFFER_SEGMENTS_POLAR
         MIN_PERCENT_PRESS_VALUES = MIN_PERCENT_PRESS_VALUES_POLAR
         IR_GRADIENT_THRESHOLD = IR_GRADIENT_THRESHOLD_POLAR
         CORR_MIN = CORR_MIN_POLAR
      ELSE
         BUFFER_SEGMENTS = BUFFER_SEGMENTS_GEO
         MIN_PERCENT_PRESS_VALUES = MIN_PERCENT_PRESS_VALUES_GEO
         IR_GRADIENT_THRESHOLD = IR_GRADIENT_THRESHOLD_GEO
         CORR_MIN = CORR_MIN_GEO
      ENDIF





      IF (targetType == TARGET_CLEAR) CALL fw_log_info(FUNC, 'Clear sky processing active')
      ! This check was added to prevent nested tracking from being used with the clear-sky WV winds
      ! but, allow nested tracking to be used for polar CSWV processing
      IF ((targetType == TARGET_CLEAR) .AND. Nested_Tracking_Flag .EQ. sym%YES  .AND. TRIM(Nav%NavType) .NE. 'PS') THEN

         CALL fw_log_error(FUNC, 'Nested tracking should be off when processing clear-sky WV winds')
         RETURN

      ENDIF


      CALL fw_log_debug(FUNC, 'satellite: '//SatName)
      CALL fw_log_debug(FUNC, 'year: '//fw_to_str(middle_time_stamp%Year))
      CALL fw_log_debug(FUNC, 'day: '//fw_to_str(middle_time_stamp%julian_day))
      CALL fw_log_debug(FUNC, 'hour: '//fw_to_str(middle_time_stamp%Hour))
      CALL fw_log_debug(FUNC, 'minute: '//fw_to_str(middle_time_stamp%Minute))
      CALL fw_log_debug(FUNC, 'AMV Channel Index: '//fw_to_str(Channel))
      CALL fw_log_info(FUNC, 'Instrument Channel: '//ChannelName)
      WRITE(string,'(f10.2)') Time_Interval_minus
      CALL fw_log_info(FUNC, 'time interval backward (min): '//string)
      WRITE(string,'(f10.2)') Time_Interval_plus
      CALL fw_log_info(FUNC, 'time interval forward (min): '//string)
      WRITE(string,'(f10.2)') Time_Interval
      CALL fw_log_info(FUNC, 'Time_Interval: '//string)
      CALL fw_log_info(FUNC, 'target box size: '//fw_to_str(Box_Size))
      CALL fw_log_debug(FUNC, 'lag size: '//fw_to_str(Lag_Size))
      CALL fw_log_info(FUNC, 'Nested tracking flag: '//fw_to_str(Nested_Tracking_Flag))
      CALL fw_log_info(FUNC, "Product Resolution: "//resolution)
      ! aab Tandem mod
      IF (TRIM(Nav%NavType) .EQ. 'PS' ) THEN
         CALL fw_log_info(FUNC, "Tandem Triplet Configuration: "//TandemTriplet)
      ENDIF
      !CALL fw_log_debug(FUNC, 'Flag to read cloud file: '//fw_to_str(Temporal_Cloud_Flag))


      IF(.NOT.ds_get(INPUT_NWP_PRESSURE_LEVELS_LBL,forecastProfiles%PressureProfiles,FUNC)) RETURN
      forecastProfiles%Pressure => forecastProfiles%PressureProfiles(1,1,:)
      IF(.NOT.ds_get(INPUT_NWP_WND_U_PROF_LBL,forecastProfiles%U_Wind,FUNC)) RETURN
      IF(.NOT.ds_get(INPUT_NWP_WND_V_PROF_LBL,forecastProfiles%V_Wind,FUNC)) RETURN
      IF(.NOT.ds_get(INPUT_NWP_TEMPERATURE_PROF_LBL,forecastProfiles%Temperature,FUNC)) RETURN

      IF(.NOT. input%Load(isTrackingSkipped)) RETURN



      IF(.NOT. isTrackingSkipped) THEN
         ALLOCATE(Wind_Data(MAX_WIND_RECORDS, MAX_PARAMETERS))
         Wind_Data = MISSING_VALUE_REAL4
         IF(.NOT. TargetData%Create(MAX_WIND_RECORDS)) RETURN

         IF(.NOT. Target_Selection_Loop(&
                                    CORR_MIN, BUFFER_SEGMENTS, MAX_WIND_RECORDS, &
                                    Nested_Tracking_Flag,targetType, Lag_Size, TandemTriplet, &  !<--- aab Tandem mod
                                    ctx_line, &
                                    elemROI, lineROI, & !inputs
                                    ImageStartElem_native, ImageStartLine_native, & !inputs
                                    EpochImageTimes,Time_Interval_minus,Time_Interval_plus,&
                                    Nav, input, forecastProfiles, & !inputs
                                    SatLat, SatLon, Space_Mask, & !inputs
                                    Target_Number, Good_Targets,& !outputs ! Good_Winds is a global
                                    TargetData, Error_Counts, Wind_Data)& ! outputs
                                    &) RETURN
      ELSE
         CALL fw_log_warning(FUNC, 'Expected error encountered: skipping processing of channel .')
      ENDIF




      
      ! Put the values into cache
      
      IF(.NOT.ds_put(makeSegmentURI(ADHOC_URI_SCHEME, ERROR_COUNTS_ADHOC), Error_Counts,FUNC)) RETURN

      IF(Target_Number > 0) THEN
         IF(.NOT.ds_put(makeSegmentURI(ADHOC_URI_SCHEME, WIND_DATA_ADHOC), Wind_Data(1:Target_Number,:),FUNC)) RETURN
         IF(.NOT.ds_put(makeSegmentURI(ADHOC_URI_SCHEME, FIRST_TARGET_EPOCH_LBL),&
                                       TargetData%EpochTimes(1:Target_Number,FIRST),FUNC)) RETURN
         IF(.NOT.ds_put(makeSegmentURI(ADHOC_URI_SCHEME, MIDDLE_TARGET_EPOCH_LBL),&
                                       TargetData%EpochTimes(1:Target_Number,MIDDLE),FUNC)) RETURN
         IF(.NOT.ds_put(makeSegmentURI(ADHOC_URI_SCHEME, LAST_TARGET_EPOCH_LBL),&
                                       TargetData%EpochTimes(1:Target_Number,LAST),FUNC)) RETURN
      ENDIF

      IF(.NOT.ds_put(makeSegmentURI(ADHOC_URI_SCHEME, STAT_NUM_TARGETS_LBL), Target_Number,FUNC)) RETURN
      IF(.NOT.ds_put(makeSegmentURI(ADHOC_URI_SCHEME, STAT_NUM_GOOD_WINDS_LBL), Good_Winds,FUNC)) RETURN
      IF(.NOT.ds_put(makeSegmentURI(ADHOC_URI_SCHEME, STAT_NUM_GOOD_TARGETS_LBL), Good_Targets, FUNC)) RETURN

      CALL ErrorCountsLog(Error_Counts,Target_Number,Good_Targets,Good_Winds)


      IF(ALLOCATED(Wind_Data)) DEALLOCATE(Wind_Data)


      Stat=RETURN_SUCCESS

   END SUBROUTINE Target_Selection_Main

   !> targeting loop that will step through the segment. targeting will not be recorded if the start is outside the region of interest
   FUNCTION Target_Selection_Loop(&
   CORR_MIN, BUFFER_SEGMENTS, MAX_WIND_RECORDS, &
   Nested_Tracking_Flag,targetType, Lag_Size, TandemTriplet, &  !<-- aab Tandem mod
   ctx_line, &
   elemROI, lineROI, & !inputs
   ImageStartElem_native, ImageStartLine_native, & !inputs
   EpochImageTimes,Time_Interval_minus,Time_Interval_plus,&
   Nav, input, forecastProfiles, & !inputs
   SatLat, SatLon, Space_Mask, & !inputs
   Target_Number, Good_Targets,& !outputs ! Good_Winds is a global
   TargetData, Error_Counts, Wind_Data) RESULT(ret) ! outputs
      CHARACTER(len=*), PARAMETER :: FUNC = "Target_Selection_Loop"
      !Navigation inputs
      REAL(SINGLE), INTENT(IN)   :: CORR_MIN !< correlation must be greater than this value.
      INTEGER(SHORT), INTENT(IN) :: BUFFER_SEGMENTS !< the original winds implementation had this.
                                                    !  Tracking starts at the center line of a buffer.
                                                    !  The buffer is roughly BUFFER_SEGMENTS times the nominal box size.
                                                    !  The buffer is then incremented towards the bottom of the image
                                                    !  until there is no more padding remaining in the segment.
                                                    !  All it's doing now is preventing tracking from starting
                                                    ! at the beginning and ending of an image
                                                    !  and that isn't especially desirable anyway.
      INTEGER(LONG), INTENT(IN) :: MAX_WIND_RECORDS !< The winds record has a buffer size equal to this
      INTEGER(LONG), INTENT(IN) :: Nested_Tracking_Flag !< flag indicating nested tracking is enabled
      CHARACTER(len=*), INTENT(IN) :: targetType !< Cloud Top ('CT') or Clear Sky ('CS')
      INTEGER(LONG), INTENT(IN) :: Lag_Size !< Lag Size used in the search
      !aab Tandem mod
      CHARACTER(len=*), INTENT(IN) :: TandemTriplet !< Tandem-VPW Triplet Configuration
      type(fw_data_range), INTENT(IN):: elemROI !< the region of interest within the padded segment
      type(fw_data_range), INTENT(IN):: lineROI !< the region of interest within the padded segment
      type(fw_data_context), INTENT(IN):: ctx_line !< context. Describes location within the larger picture
      INTEGER(DLONG), INTENT(IN) :: ImageStartElem_native !< segment start x (element)
                                                          !< in the native image coordinates (0 or 1 at origin)
      INTEGER(DLONG), INTENT(IN) :: ImageStartLine_native !< segment start y (line)
                                                          !< in the native image coordinates (0 or 1 at origin)
      TYPE(Nav_t), INTENT(IN) :: Nav !< saves all the navigation parameters needed to map between x/y and lat/lon
      TYPE(input_t), INTENT(IN) :: input !< stores all the inputs
      REAL(SINGLE), INTENT(IN) :: SatLat !< satellite latitude subpoint
      REAL(SINGLE), INTENT(IN) :: SatLon !< satellite longitude subpoint
      INTEGER(BYTE), DIMENSION(:,:), INTENT(IN) :: Space_Mask !< space mask of the middle granule
      TYPE(ForecastProfilesGrid_t), INTENT(IN) :: forecastProfiles
      INTEGER(LONG), DIMENSION(NUMBER_TIMESTEPS), INTENT(IN) :: EpochImageTimes !< granule image times.
      REAL(SINGLE), INTENT(IN)  :: Time_Interval_minus !< granule image time interval minus (units: seconds)
      REAL(SINGLE), INTENT(IN)  :: Time_Interval_plus !< granule image time interval plus (units: seconds)

      ! Outputs Arrays
      INTEGER(LONG), INTENT(OUT) :: Target_Number
      INTEGER(LONG), INTENT(OUT) :: Good_Targets
      !INTEGER(LONG), INTENT(OUT) :: Good_winds ! global currently.
      INTEGER(LONG), INTENT(INOUT) :: Error_Counts(-10:30)
      REAL(SINGLE), INTENT(INOUT) :: Wind_Data(:, :)
      TYPE(TargetData_t), INTENT(INOUT) :: TargetData
      LOGICAL :: ret

      INTEGER(LONG) :: Elems
      INTEGER(LONG) :: Lines


      INTEGER(LONG):: Press_Levels
      INTEGER(LONG) :: X_NWP_Target
      INTEGER(LONG) :: Y_NWP_Target
      INTEGER(BYTE) :: LL_Inver_Flag !< low-level inversion flag


      TYPE(DBSCAN_Output) :: DBSCAN_Out


      ! position variables
      INTEGER(LONG) :: gradientLine !< line within the gradient loop
      INTEGER(LONG) :: gradientElem !< element within the gradient loop
      INTEGER(LONG) :: targetingLines !< lines available in the scene for targeting
      INTEGER(LONG) :: cutoff !< the lines to cut off before targetSceneStartLine
      INTEGER(LONG) :: targetSceneStartLine !< Line at the start of the "buffer" scene
      INTEGER(LONG) :: targetSceneEndLine !< Line at the end of the "buffer" scene
      INTEGER(LONG) :: gradientTrackingCenterLine !< line center of the "buffer" scene
      INTEGER(LONG) :: gradientTrackingCenterElem !< element center of the "buffer" scene
      INTEGER(LONG) :: gradientTrackingStartLine !< start of the gradient tracking
      INTEGER(LONG) :: gradientTrackingEndLine !< end of the the gradient tracking
      INTEGER(LONG) :: gradientTrackingStartElem !< start of the gradient tracking
      INTEGER(LONG) :: gradientTrackingEndElem !< end of the the gradient tracking
      INTEGER(LONG) :: Target_Center_Elem !< target center element relative to the segment
      INTEGER(LONG) :: Target_Center_Line !< target center line relative to the segment
      INTEGER(LONG), DIMENSION(Location_Dimensions) :: Max_Location !< location of maximum gradient magnitude in target box
      REAL(SINGLE) :: Target_Elem_native !< target x in the native image coordinates (0 or 1 based)
      REAL(SINGLE) :: Target_Line_native !< target y in the native image coordinates (0 or 1 based)


      TYPE(Output_Variables) :: Output


      REAL(SINGLE) :: Delta_Time_minus
      REAL(SINGLE) :: Delta_Time_plus
      INTEGER(LONG), DIMENSION(NUMBER_TIMESTEPS) :: EpochScanTimes


      ! image coordinates
      INTEGER(LONG) :: Target_Half_Width
      INTEGER(LONG) :: Search_Half_Width_Element
      INTEGER(LONG) :: Search_Half_Width_Line

      INTEGER(LONG) :: Diff_In_Size


      INTEGER(LONG) :: QC_Flag !< target quality control flag


      ! BrtTemp gradient information
      INTEGER(LONG) :: Gradient_Offset
      REAL(SINGLE) :: Gradient_X
      REAL(SINGLE) :: Gradient_Y
      REAL(SINGLE) :: Gradient_Max
      REAL(SINGLE) :: Gradient_Magnitude


      ! displacement of search box center from target box center
      REAL(SINGLE) :: Delta_X_Fcst
      REAL(SINGLE) :: Delta_Y_Fcst


      ! pixel factor
      REAL(SINGLE) :: Factor

      ! earth coordinates of target, search and match regions
      REAL(DOUBLE) :: Latitude  !< Lat_Target in double precision
      REAL(DOUBLE) :: Longitude !< Lon_Target in double precision
      REAL(SINGLE) :: Lat_Target
      REAL(SINGLE) :: Lon_Target
      REAL(SINGLE) :: Lat_Match
      REAL(SINGLE) :: Lon_Match
      REAL(SINGLE) :: Lat_Match2
      REAL(SINGLE) :: Lon_Match2

      ! end points for nested tracking estimate
      REAL(SINGLE) :: Lat_Match3
      REAL(SINGLE) :: Lon_Match3
      REAL(SINGLE) :: Lat_Match4
      REAL(SINGLE) :: Lon_Match4


      ! tracking variables
      INTEGER(LONG) :: Status_Flag1
      INTEGER(LONG) :: Status_Flag2
      INTEGER(LONG) :: Least_Square_Error
      REAL(SINGLE) :: Line_Disp
      REAL(SINGLE) :: Element_Disp
      REAL(SINGLE) :: Line_Disp2
      REAL(SINGLE) :: Element_Disp2

      REAL(DOUBLE) :: Target_Scene_Mean
      REAL(DOUBLE) :: Target_Scene_StdDev
      REAL(DOUBLE) :: Match_Scene_Mean
      REAL(DOUBLE) :: Match_Scene_StdDev


      REAL(SINGLE) :: U_Avg
      REAL(SINGLE) :: V_Avg
      REAL(SINGLE) :: U_Component !< u wind of reverse vector (negative time interval)
      REAL(SINGLE) :: V_Component !< v wind of reverse vector (negative time interval)
      REAL(SINGLE) :: U_Component2 !< u wind of forward vector (positive time interval)
      REAL(SINGLE) :: V_Component2 !< v wind of forward vector (positive time interval)

      ! Data from largest clusters found through DBSCAN analysis
      ! First sample is for reverse vector, second sample is for forward vector
      REAL(SINGLE) :: Median_CTP_Sample1
      REAL(SINGLE) :: Median_CTP_Sample2
      REAL(SINGLE) :: Combined_Median
      REAL(SINGLE) :: Combined_Median_Hgt


      ! height assignment variables
      REAL(SINGLE) :: Median_Press2
      REAL(SINGLE) :: Median_BrtTemp


      INTEGER(LONG), DIMENSION(:,:), ALLOCATABLE :: Lag_Table
      INTEGER(LONG), DIMENSION(:,:), ALLOCATABLE :: Lag_Table2
      REAL(SINGLE), DIMENSION(:,:), ALLOCATABLE :: Search_BrtTemp_Values
      REAL(SINGLE), DIMENSION(:,:), ALLOCATABLE :: Search_BrtTemp_Values2
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:,:) :: SetOfPoints
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: SSD_Matrix

      TYPE(Forecast_Profiles) :: Fcst_Prof
      TYPE(Target_Box_Data) :: Box_Data
      INTEGER(LONG) :: Box_Element_Start
      INTEGER(LONG) :: Box_Element_End
      INTEGER(LONG) :: Box_Line_Start
      INTEGER(LONG) :: Box_Line_End

      INTEGER(LONG) :: Search_Size
      INTEGER(LONG) :: Lag_Array_Size
      INTEGER(LONG) :: Start_Element
      INTEGER(LONG) :: End_Element
      INTEGER(LONG) :: Level
      TYPE(Uniformity_t) :: Unif

      REAL(SINGLE) LD_Test  !< avg line displacement of pts in largest cluster (reverse vector)
      REAL(SINGLE) ED_Test  !< avg ele displacement of pts in largest cluster (reverse vector)
      REAL(SINGLE) CC_Test  !< linear correlation of matching scene (reverse vector)
      INTEGER(LONG) SF_Test  !< quality flag (reverse vector)

      REAL(SINGLE) LD_Test2  !< avg line displacement of pts in largest cluster (forward vector)
      REAL(SINGLE) ED_Test2  !< avg ele displacement of pts in largest cluster (forward vector)
      REAL(SINGLE) CC_Test2  !< linear correlation of matching scene (forward vector)
      INTEGER(LONG) SF_Test2  !< quality flag (forward vector)

      ! aab Tandem (dbscan)
      INTEGER(LONG) comb_size
      REAL(SINGLE) Sum_DX1, Sum_DY1, Avg_DX1, Avg_DY1, Sum_DX2, Sum_DY2, Avg_DX2, Avg_DY2

      ret = .FALSE.

      Target_Number = 0
      Good_Targets = 0



      SELECT CASE(targetType)
         CASE (TARGET_CLEAR)
            CSWV_Flag = CSWV_CLEAR
         CASE (TARGET_CLOUD)
            CSWV_Flag = CSWV_CLOUD
         CASE DEFAULT
            CALL fw_log_error(FUNC, "Target Type is not properly assigned")
            RETURN
      END SELECT


      !For most satellites the scan time of target is assumed constant
      Delta_Time_minus = Time_Interval_minus / MIN_IN_A_HOUR
      Delta_Time_plus = Time_Interval_plus / MIN_IN_A_HOUR

      Lines = SIZE(Space_Mask,2)
      Elems = SIZE(Space_Mask,1)

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

      Press_Levels = size(forecastProfiles%Pressure)

      Total_Points = Box_Size**2
      Target_Half_Width = Box_Size / 2

      ! the lag array has one additional element to account for center position (lag position zero)
      Lag_Array_Size = Lag_Size + MOD(Lag_Size+1,2)
      Search_Size = Box_Size + (Lag_Size / 2) * 2
      Search_Half_Width_Element = (Search_Size - 1) / 2
      Search_Half_Width_Line = Search_Half_Width_Element



      CALL fw_log_info(FUNC, 'Lag Array Size: '//fw_to_str(Lag_Array_Size))
      CALL fw_log_info(FUNC, 'Search Size: '//fw_to_str(Search_Size))
      CALL fw_log_info(FUNC, 'Search Half Width Element: '//fw_to_str(Search_Half_Width_Element))
      CALL fw_log_info(FUNC, 'Dimensions of search (elements x lines): '//&
                             fw_to_str(2 * Search_Half_Width_Element + 1)//&
                             " x "//fw_to_str(2 * Search_Half_Width_Line + 1))




      ! Allocate memory for the target box data structure
      IF(.NOT. Box_Data%Create(Box_Size)) RETURN



      ! Allocate memory for the forecast profiles data structure
      IF(.NOT. Fcst_Prof%Create(Press_Levels)) RETURN

      ! Allocate memory for the arrays holding the search image data
      ! (previous time and later time)
      ALLOCATE(Search_BrtTemp_Values(Search_Size, Search_Size))
      ALLOCATE(Search_BrtTemp_Values2(Search_Size, Search_Size))
      IF (SIZE(Search_BrtTemp_Values) .LE. Total_Points) THEN
            CALL fw_log_error(FUNC, 'Search box must be larger than target scene')
         RETURN
      ENDIF

      IF(.NOT. DBSCAN_Out%Create(Box_Size)) RETURN

      ! aab Tandem (dbscan)
      !ALLOCATE(SetOfPoints((Box_Size-4)**2, DBSCAN_Output_Variables))
      ! Double size of first dimension of SetOfPoints array for combined sample testing
      comb_size = ((Box_Size-4)**2) * 2
      ALLOCATE(SetOfPoints(comb_size, DBSCAN_Output_Variables))
      
      SetOfPoints = MISSING_VALUE_REAL4

      IF (.NOT. ALLOCATED(SSD_Matrix)) ALLOCATE(SSD_Matrix(Lag_Array_Size**2))
      SSD_Matrix = MISSING_VALUE_REAL4

      ! --------------------------------------------------------------------------
      ! Generate tables of positions to conduct a "spiral" search beginning at the
      ! center of the lag matrix (which corresponds to the forecast position) and
      ! spiraling outward. This is a more efficient search for the SSD minimum.
      ! --------------------------------------------------------------------------
      ALLOCATE(Lag_Table(Lag_Array_Size**2, 4))
      ALLOCATE(Lag_Table2(Lag_Array_Size**2, 4))

      Diff_In_Size = 0
      CALL Generate_Spiral_Table(Search_Size, Lag_Array_Size, Diff_In_Size, Lag_Table)

      !the difference between the target box size and the size of the nested tracking box
      !(5 pixels)
      Diff_In_Size = Box_Size - 5
      CALL Generate_Spiral_Table(Search_Size, Lag_Array_Size, Diff_In_Size, Lag_Table2)

      !---------------------------------------------------------------------------
      ! Call GEOCAT routine to compute 3x3 mean values for entire buffer.
      ! Although only computing the gradient for the middle segment, we will
      ! need the extra data when the target box is re-centered on the max
      ! gradient.
      !---------------------------------------------------------------------------
      CALL Compute_Spatial_Uniformity(1, 1, input%Space%d, input%Radiance%d, &
         Unif%mean, Unif%max, Unif%min, Unif%uni)



      ! Begin processing segment data
      ! the segment loop
      !More targets can be found if it is just BUFFER_SEGMENTS*NOMINAL_BOX_SIZE
      targetingLines = (BUFFER_SEGMENTS - 1)*NOMINAL_BOX_SIZE + Box_Size
      !we want to cut off the beginning if its not a multiple of NOMINAL_BOX_SIZE
      cutoff = NOMINAL_BOX_SIZE - MOD(ctx_line%position%start_pos - 1, NOMINAL_BOX_SIZE)
      IF(cutoff == NOMINAL_BOX_SIZE) cutoff = 0
      targetSceneStartLine = 1 + cutoff
      targetSceneEndLine = targetSceneStartLine + targetingLines - 1
      IF(targetSceneEndLine > lines) THEN
         CALL fw_log_error(FUNC, "Framework Segmentation is insufficient for targeting")
      ENDIF




      !Initialize the scan line times with the nominal image time.
      !Some satellites just use the image time as the scan line time.
      EpochScanTimes(FIRST) = EpochImageTimes(FIRST)
      EpochScanTimes(MIDDLE) = EpochImageTimes(MIDDLE)
      EpochScanTimes(LAST) = EpochImageTimes(LAST)



      ! pixel offset for gradient calculation loop
      IF (Box_Size .GE. NOMINAL_BOX_SIZE) THEN
         Gradient_Offset = (Box_Size - NOMINAL_BOX_SIZE) / 2
      ELSE
         Gradient_Offset = 0
      ENDIF
      Gradient_X = 0.0
      Gradient_Y = 0.0

      ! The segment loop (line loop)
      DO WHILE(targetSceneEndLine <= lines)

         ! Initialize target box element bounds
         Start_Element = 1 + PIXEL_OFFSET
         !---------------------------------------------------------------------------
         ! Beginning of main target selection loop. The middle portion of
         ! the buffer is essentially divided into small sub-regions, the size
         ! of which is controlled by the Box_Size variable. At each pixel
         ! of each box the brightness temperature gradient is computed using a
         ! 4-pt centered difference convolution kernel. The box is then
         ! recentered on the maximum gradient location. Once this is done a
         ! variety of tests are performed to determine if the target is a
         ! suitable tracer. The loop is terminated when the edge of the data
         ! is reached.
         !---------------------------------------------------------------------------

         ! The main loop (element loop)
         DO

            !-------------------------------------------------------------------------
            ! Initialize data structures and target box variables
            !-------------------------------------------------------------------------
            CALL Output_Initialize(Output)
            CALL DBSCAN_Output_Initialize(DBSCAN_Out)
            CALL Box_Data_Initialize(Box_Data)

            Output%Channel = Channel
            Gradient_Max = 0.0
            Max_Location(1) = 0
            Max_Location(2) = 0

            QC_Flag = sym%SUCCESS
            Lat_Target = MISSING_VALUE_REAL4
            Lon_Target = MISSING_VALUE_REAL4
            Lat_Match = MISSING_VALUE_REAL4
            Lon_Match = MISSING_VALUE_REAL4
            Lat_Match2 = MISSING_VALUE_REAL4
            Lon_Match2 = MISSING_VALUE_REAL4
            Lat_Match3 = MISSING_VALUE_REAL4
            Lon_Match3 = MISSING_VALUE_REAL4
            Lat_Match4 = MISSING_VALUE_REAL4
            Lon_Match4 = MISSING_VALUE_REAL4
            U_Avg = MISSING_VALUE_REAL4
            V_Avg = MISSING_VALUE_REAL4
            U_Component = MISSING_VALUE_REAL4
            V_Component = MISSING_VALUE_REAL4
            U_Component2 = MISSING_VALUE_REAL4
            V_Component2 = MISSING_VALUE_REAL4
            Median_Press2 = MISSING_VALUE_REAL4
            Median_BrtTemp = MISSING_VALUE_REAL4
            X_NWP_Target = MISSING_VALUE_INT4
            Y_NWP_Target = MISSING_VALUE_INT4
            Delta_X_Fcst = MISSING_VALUE_REAL4
            Delta_Y_Fcst = MISSING_VALUE_REAL4
            Median_CTP_Sample1 = MISSING_VALUE_REAL4
            Median_CTP_Sample2 = MISSING_VALUE_REAL4
            Combined_Median = MISSING_VALUE_REAL4
            Combined_Median_Hgt = MISSING_VALUE_REAL4
            Output%Flag_5x5 = Nested_Tracking_Flag
            SF_Test = sym%FAILURE
            SF_Test2 = sym%FAILURE


            !-------------------------------------------------------------------------
            ! Determine the end element for the current target box. Exit loop
            ! if beyond the edge of the image.
            !-------------------------------------------------------------------------
            End_Element = Start_Element + Box_Size - 1

            !-------------------------------------------------------------------------
            ! Loop over every pixel in target scene and compute gradient.
            !-------------------------------------------------------------------------

            gradientTrackingCenterLine = (targetSceneEndLine - targetSceneStartLine) / 2 + targetSceneStartLine
            gradientTrackingStartLine = gradientTrackingCenterLine - Target_Half_Width + Gradient_Offset
            gradientTrackingEndLine = gradientTrackingStartLine + (Box_Size - Gradient_Offset*2) - 1
            gradientTrackingStartElem = Start_Element + Gradient_Offset
            gradientTrackingEndElem = End_Element - Gradient_Offset
            gradientTrackingCenterElem = &
		(gradientTrackingEndElem - gradientTrackingStartElem) / 2 + gradientTrackingStartElem


            IF (gradientTrackingEndElem .GT. (Elems - PIXEL_OFFSET)) EXIT !This one is needed to prevent a possible segfault.
            IF (End_Element .GT. (Elems - PIXEL_OFFSET)) EXIT  !FIXME:  put check before convolve like the line check.

            !if the tracking starts after the ROI then a later segment will pick it up.
            IF( gradientTrackingStartElem < elemROI%start_pos .OR. gradientTrackingStartElem > elemROI%end_pos .OR. &
                gradientTrackingStartLine < lineROI%start_pos .OR. gradientTrackingStartLine > lineROI%end_pos ) THEN
                !increment as if successful. a bit faster to just exit the element loop in the event of a line boundary issue.
                CALL IncrementElement(sym%SUCCESS, Start_Element)
                CYCLE
            ENDIF


            DO gradientLine = gradientTrackingStartLine, gradientTrackingEndLine
               DO gradientElem = gradientTrackingStartElem, gradientTrackingEndElem

                  ! Initialize gradient magnitude to zero at each pixel.
                  Gradient_Magnitude = 0.0

                  ! Only proceed if pixel is non-space look.
                  IF (Space_Mask(gradientElem,gradientLine) == SYM%NO_SPACE) THEN

                     !-----------------------------------------------------------------
                     ! Calculate the IR brightness temperature gradient.
                     ! Using buffer data will allow us to compute a 4-pt centered
                     ! difference for every pixel in the target box.
                     !-----------------------------------------------------------------

                     IF ( (gradientLine + PIXEL_OFFSET) > lines ) CYCLE !alternatively alter the convolve
                     Gradient_X = Convolve_X(Kernel,gradientElem,gradientLine, input%BrtTemp%d, INT(PIXEL_OFFSET,LONG))

                     Gradient_Y = Convolve_Y(Kernel,gradientElem,gradientLine, input%BrtTemp%d, INT(PIXEL_OFFSET,LONG))

                     Gradient_Magnitude = SQRT(Gradient_X**2 + Gradient_Y**2)

                     IF (Gradient_Max < Gradient_Magnitude) THEN
                        Gradient_Max = Gradient_Magnitude
                        Max_Location(1) = gradientElem
                        Max_Location(2) = gradientLine
                     ENDIF

                  ENDIF

               END DO
            END DO
            ! If gradient max is zero, go to next target.
            IF (Gradient_Max .EQ. 0.0) THEN

               !We already checked if the start of the tracking is within the ROI so we can store this as a failure.
               QC_Flag = ZERO_GRADIENT_FAILURE
               !Update_Wind_Data is here in case we want to start storing ZERO_GRADIENT_FAILURE failures.
               !right now it does not even enter the missing value lat/lon if block inside Update_Wind_data.
               CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,& 
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
               CALL IncrementElement(QC_Flag, Start_Element)
               CYCLE
            ENDIF

            ! ------------------------------------------------------------------------
            ! Center the target on the maximum gradiant
            ! ------------------------------------------------------------------------
            Target_Center_Elem = Max_Location(1)
            Target_Center_Line = Max_Location(2)


            ! Convert to the native coordinates needed for output and pixcoord2geocoord (e.g. ABI is 0 indexed.)
            Target_Elem_native = REAL(ImageStartElem_native + Target_Center_Elem - 1)
            Target_Line_native = REAL(ImageStartLine_native + Target_Center_Line - 1)

            ! Recompute boundaries. Do not exceed image limits.
            Box_Element_Start = MAX(1, Target_Center_Elem - Target_Half_Width)
            Box_Element_End = MIN(Elems, Target_Center_Elem + Target_Half_Width)
            Box_Line_Start = MAX(1, Target_Center_Line - Target_Half_Width)
            Box_Line_End = MIN(Lines, Target_Center_Line + Target_Half_Width)



            ! Call navigation to get the earth location of target center pixel.
            IF (TRIM(Nav%NavType) .EQ. 'PS' ) THEN

               CALL pixcoord2geocoord_ps(Nav, Target_Elem_native, Target_Line_native, Lat_Target, Lon_Target)

               Lon_Target = -Lon_Target
               Output%Lon_Target = Lon_Target

            ELSE

               SELECT CASE (Sensor_Series)

                  CASE (SERIES_MSG_SEVIRI)

                     ! Convert to EUMETSAT coord
                     ! 3713 is the number of lines  (or elements) in a SEVIRI full disk image plus 1.
                     !Target_Elem_native = REAL(3713 - (ImageStartElem_native + Target_Center_Elem))
                     !Target_Line_native = REAL(3713 - (Center_Of_Buffer + Lines_From_Center))

                     ! compute earth coordinates of target center
                     CALL pixcoord2geocoord_dp(Nav, Target_Elem_native, Target_Line_native, Latitude, Longitude)

                     Lat_Target = REAL(Latitude)
                     Lon_Target = REAL(Longitude)

                     Output%Lon_Target = Lon_Target

                  !      CASE (SERIES_GOES_N_IMG)
                  !
                  !        CALL pixcoord2geocoord_gvar(Target_Elem_native, Target_Line_native, &
                  !                            sat_xstride, Lat_Target, Lon_Target)
                  !
                  !        Output%Lon_Target = Lon_Target

                  !        EpochScanTimes(FIRST) = input%Line_Time_Image0(Channel,Max_Location(2))
                  !        EpochScanTimes(MIDDLE) = input%Line_Time_Image1(Channel,Max_Location(2))
                  !        EpochScanTimes(LAST) = input%Line_Time_Image2(Channel,Max_Location(2))
                  !        Delta_Time_minus = seconds_to_hours(EpochScanTimes(MIDDLE) - EpochScanTimes(FIRST))
                  !        Delta_Time_plus  = seconds_to_hours(EpochScanTimes(LAST) - EpochScanTimes(MIDDLE))

                  CASE (SERIES_GOESR_ABI, SERIES_HIMAWARI_AHI)

                     CALL pixcoord2geocoord_abi_real(Nav, Target_Elem_native, Target_Line_native, &
                        Latitude, Longitude)


                     Lat_Target = REAL(Latitude)
                     Lon_Target = REAL(Longitude)

                     Output%Lon_Target = Lon_Target

                     ! 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_Target > 180.0) Lon_Target = Lon_Target - 360.0


               END SELECT
            ENDIF

            ! Compute ratio of 1-pixel element displacement over 1-pixel line
            ! displacement. This ratio will be passed to DBSCAN and used when
            ! computing distance in pixel space.

            CALL Pixel_Factor(Nav, Target_Elem_native, Target_Line_native, Factor)

            Output%Lat_Target = Lat_Target

            ! add image coordinates of target center - wcb 12/2/15
            Output%Target_Element = Target_Elem_native
            Output%Target_Line = Target_Line_native

            ! Check the Land Mask and determine if the target is over land or water.
            ! If it's over land, set output flag to 1.
            ! If it's over water, set output flag to 0.

            IF (input%Land%d(Max_Location(1),Max_Location(2)) == sym%Land .OR. &
               input%Land%d(Max_Location(1),Max_Location(2)) == sym%COASTLINE) THEN

               Output%Land_Flag = LAND_BINARY_MASK

            ELSE

               Output%Land_Flag = WATER_BINARY_MASK

            ENDIF

            IF (TRIM(Nav%NavType) .NE. 'PS' ) THEN

               ! Compute satellite zenith angle (geostationary processing only)
               CALL Compute_Zenith_Angle(SatLon, SatLat, Lon_Target, &
                  Lat_Target, Output%Sat_Zenith_Angle)

               ! WCB
               ! add zenith angle cut off
               IF (Output%Sat_Zenith_Angle .GT. ZENITH_ANGLE_CUTOFF) THEN
                  QC_Flag = EARTH_EDGE_FAILURE
                  CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                    EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,& 
                                    Box_Data) ! aab added Box_Data for CloudyPixels output
                  CALL IncrementElement(QC_Flag, Start_Element)
                  CYCLE
               ENDIF
            ENDIF

            ! Fill target box with data from buffer.
            CALL Fill_Box_Data(Box_Data, Box_Size, Box_Element_Start, Box_Element_End, &
               Box_Line_Start, Box_Line_End, input, unif%mean, unif%uni, QC_Flag)

            IF (QC_Flag .NE. sym%SUCCESS) THEN

               CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
               CALL IncrementElement(QC_Flag, Start_Element)
               CYCLE

            ENDIF

            !-------------------------------------------------------------------------
            ! Apply QC to target box. Keep track of failures.
            !-------------------------------------------------------------------------

            CALL Target_QC(Nav%NavType, Box_Data, QC_Flag)

            IF (QC_Flag .NE. sym%SUCCESS) THEN

               CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
               CALL IncrementElement(QC_Flag, Start_Element)
               CYCLE

            ENDIF

            !-------------------------------------------------------------------------
            ! Beginning of height assignment section. First step is to construct
            ! a 1-D histogram of cloud top temperature values and find the temperature
            ! threshold of the coldest 25 percent. The next step is to loop through
            ! the pixels and save the Cld_Top_Press values of those pixels colder than
            ! the threshold value. The final step is to compute the median pressure
            ! of the cold sample.
            !-------------------------------------------------------------------------

            IF (QC_Flag .EQ. sym%SUCCESS .AND. Gradient_Max > 0.0) THEN

               ! Make histogram of temperature data and return cold sample

               CALL Make_Histogram_Temp(Nav%NavType, Box_Data, Total_Points, Output%Point_Index, &
                  Output%Variance_Press, LL_Inver_Flag)

               ! Copy forecast wind and temperature profiles

               X_NWP_Target = input%X_NWP%d(Max_Location(1), Max_Location(2))
               Y_NWP_Target = input%Y_NWP%d(Max_Location(1), Max_Location(2))
               Output%X_NWP_Target = X_NWP_Target
               Output%Y_NWP_Target = Y_NWP_Target

               ! Set Inversion Flag for output file.

               IF (LL_Inver_Flag .EQ. sym%YES) THEN

                  Output%Inversion_Flag = sym%YES

               ELSE

                  Output%Inversion_Flag = sym%NO

               ENDIF

               ! copy forecast profiles at max gradient location
               CALL Forecast_Profiles_Set(Fcst_Prof, forecastProfiles,Press_Levels, X_NWP_Target, Y_NWP_Target)

               ! Find the dominant cloud phase and cloud type of the target scene.  This
               ! will simply be stored as a diagnostic output variable.

               CALL Dominant_Cloud_Phase(Box_Data, Output%Cloud_Phase)

               CALL Dominant_Cloud_Type(Box_Data, Output%Cloud_Type)

               ! ----------------------------------------------------------------------
               ! Sort pressures from target scene into ascending order and find median.
               ! We need the forecast wind at this pressure level to center the search
               ! box further below. Fail target if height assignment fails. This can
               ! happen if the 1-DVAR cloud height algorithm fails to produce a height.
               ! ----------------------------------------------------------------------

               IF (Output%Point_Index .GT. 0) THEN

                  CALL Assign_Height(Press_Levels, Output, Box_Data, Fcst_Prof, QC_Flag)

                  IF (QC_Flag .NE. sym%SUCCESS) THEN

                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ENDIF

               ENDIF ! Number of points check

               ! Check the AMV height assignment. Skip to next target if a failure
               ! is encountered. The primary role of this subroutine is to apply
               ! channel-specific pressure thresholds.

               CALL Check_Height_Assignment(Output%Point_Index, Output%Median_Press, &
                  Fcst_Prof%Pressure, QC_Flag)

               IF (QC_Flag .NE. sym%SUCCESS) THEN

                  CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                  CALL IncrementElement(QC_Flag, Start_Element)
                  CYCLE

               ENDIF

               ! ----------------------------------------------------------------------
               ! At this point the target has passed all quality control and is
               ! suitable for tracking.
               ! ----------------------------------------------------------------------

               Good_Targets = Good_Targets + 1

               Box_Data%Target_ID = Target_Number

               ! Interpolate forecast to AMV height.

               CALL Interpolate_To_AMV_Press(Output, Press_Levels, Fcst_Prof)

               ! Compute speed and direction of forecast at original pressure
               ! assignment. The pressure assignment will change if nested
               ! tracking is being used, so save the original value.

               CALL UV2SPD(Output%Spd_Fcst_Original_Press, &
                  Output%Dir_Fcst_Original_Press, &
                  Output%U_Fcst, Output%V_Fcst)

               !-----------------------------------------------------------------------
               ! Begin tracking portion of AMV algorithm. The forecast wind is used as
               ! a "guide" to position the search region in the subsequent image. To
               ! do this we need to compute the displacement predicted by the model.
               !-----------------------------------------------------------------------

               ! Retrieve forecast position of feature in previous image (negative
               ! time step)

               CALL Get_Fcst_Displacement(Nav,Nav%NavType, Lat_Target, Lon_Target, -Delta_Time_minus, &
                  Output%U_Fcst, Output%V_Fcst, Sensor_Series, Target_Elem_native, &
                  Target_Line_native, Delta_X_Fcst, Delta_Y_Fcst)

               !-----------------------------------------------------------------------
               ! Use predicted displacement to center search box and fill with data
               ! from buffer.
               !
               ! This routine starts with the target position (max gradient location)
               ! in buffer coordinates and adds the predicted line and element
               ! displacement to determine where to center the search box in buffer
               ! coordinates. It then fills the search box with data from the buffer.
               !-----------------------------------------------------------------------

               CALL Fill_Search_Box(Elems, Max_Location(1), Max_Location(2), Delta_X_Fcst, &
                  Delta_Y_Fcst, Search_Half_Width_Line, Search_Half_Width_Element, &
                  Lines, input%Search%d, Search_BrtTemp_Values, QC_Flag)
               ! superficially restrict the search boundaries to the "buffer only"
               !CALL Fill_Search_Box(Elems, Max_Location(1), Max_Location(2)-targetSceneStartLine + 1, Delta_X_Fcst, &
               !   Delta_Y_Fcst, Search_Half_Width_Line, Search_Half_Width_Element, &
               !   targetingLines, input%Search%d(:,targetSceneStartLine:targetSceneEndLine), Search_BrtTemp_Values, QC_Flag)


               IF (QC_Flag .NE. sym%SUCCESS) THEN

                  CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                  CALL IncrementElement(QC_Flag, Start_Element)
                  CYCLE

               ENDIF

               ! check to see if we are tracking with the traditional method

               IF (Nested_Tracking_Flag .EQ. sym%NO) THEN

                  ! Find matching scene.

                  CALL Compute_Correlation(Box_Data, Search_BrtTemp_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, Output%Corr_Coeff, &
                     Status_Flag1, SSD_Matrix)

                  SSD_Matrix = MISSING_VALUE_REAL4

                  ! Check success of tracking

                  IF (Status_Flag1 .NE. sym%SUCCESS) THEN

                     QC_Flag = Status_Flag1
                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ENDIF

                  ! ----------------------------------------------------------------------
                  ! Check for valid correlation value. There is no point in continuing
                  ! if the value for the 1st vector is bad.
                  ! ----------------------------------------------------------------------

                  IF (Output%Corr_Coeff .LT. CORR_MIN .OR. &
                     Output%Corr_Coeff .GT. 1.0) THEN

                     QC_Flag = CORRELATION_FAILURE
                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ENDIF

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

                  ! ----------------------------------------------------------------------
                  ! Attempt nested 5x5 tracking to improve estimate
                  ! ----------------------------------------------------------------------

                  ! pass in clear sky flag - 1/25/16

                  CALL Nested_Tracking(Nav%NavType, REVERSE, CSWV_Flag, Box_Data, Search_BrtTemp_Values,Delta_X_Fcst, &
                     Delta_Y_Fcst, Channel, Lag_Table2, LD_Test, ED_Test, TandemTriplet, & !<-- aab Tandem mod
                     Least_Square_Error, Factor, Target_Scene_Mean, &
                     Match_Scene_Mean, Target_Scene_StdDev, Match_Scene_StdDev, &
                     CC_Test, SF_Test, DBSCAN_Out, SetOfPoints)

                  Output%Corr_Coeff = CC_Test
                  ! aab Tandem (dbscan)
                  ! combined sample test - keep data from reverse
                  !SetOfPoints = MISSING_VALUE_REAL4

                  ! Check success of tracking
                  IF (SF_Test .NE. sym%SUCCESS) THEN

                     QC_Flag = SF_Test
                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ENDIF

                  ! do not use CTP for CSWV - 1/25/16
                  IF (CSWV_Flag .EQ. sym%NO) THEN

                     ! --------------------------------------------------------------------
                     ! Need to find median pressure associated with the cluster of pixels
                     ! returned from the nested tracking process. We need this value
                     ! because we compare it to the value derived from the second sample
                     ! further below. A limit of 100 mb is placed on the difference.
                     ! --------------------------------------------------------------------

                     ! aab Tandem (dbscan)
                     ! DBSCAN is now not called until after the forward time step is processed, so remove this code
                     !CALL Sort_Shell_NP(DBSCAN_Out%Number_CTP_Sample1, &
                     !   DBSCAN_Out%CTP_5x5_Sample1, &
                     !   Median_CTP_Sample1)

                  ENDIF  ! CSWV check - 1/25/16


               ENDIF  ! Nested tracking check (first vector)

               ! ----------------------------------------------------------------------
               ! Beginning of logic to track feature forward in time. This logic
               ! mirrors what is done above for the reverse vector.
               ! ----------------------------------------------------------------------

               ! Retrieve forecast position of feature in next image (+Delta_Time)

               CALL Get_Fcst_Displacement(Nav,Nav%NavType, Lat_Target, Lon_Target, Delta_Time_plus, &
                  Output%U_Fcst, Output%V_Fcst, Sensor_Series, Target_Elem_native, &
                  Target_Line_native, Delta_X_Fcst, Delta_Y_Fcst)

               !-----------------------------------------------------------------------
               ! Fill search box with buffer data
               !-----------------------------------------------------------------------
               CALL Fill_Search_Box(Elems, Max_Location(1), Max_Location(2), Delta_X_Fcst, &
                  Delta_Y_Fcst, Search_Half_Width_Line, Search_Half_Width_Element, &
                  Lines, input%Search2%d, Search_BrtTemp_Values2, QC_Flag)
               ! superficially restrict the search boundaries to the "buffer only"
               !CALL Fill_Search_Box(Elems, Max_Location(1), Max_Location(2)-targetSceneStartLine+1, Delta_X_Fcst, &
               !   Delta_Y_Fcst, Search_Half_Width_Line, Search_Half_Width_Element, &
               !   targetingLines, input%Search2%d(:,targetSceneStartLine:targetSceneEndLine), Search_BrtTemp_Values2, QC_Flag)

               ! Check if search box was filled properly
               IF (QC_Flag .NE. sym%SUCCESS) THEN

                  CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                  CALL IncrementElement(QC_Flag, Start_Element)
                  CYCLE

               ENDIF

               IF (Nested_Tracking_Flag .EQ. sym%NO) THEN

                  ! Compute correlation for second vector (the forward vector)
                  CALL Compute_Correlation(Box_Data, Search_BrtTemp_Values2, &
                     Lag_Table, Delta_X_Fcst, Delta_Y_Fcst, &
                     Line_Disp2, Element_Disp2, &
                     Least_Square_Error, Target_Scene_Mean, &
                     Match_Scene_Mean, Target_Scene_StdDev, &
                     Match_Scene_StdDev, Output%Corr_Coeff2, &
                     Status_Flag2, SSD_Matrix)

                  ! Check success of tracking
                  IF (Status_Flag2 .NE. sym%SUCCESS) THEN

                     QC_Flag = Status_Flag2
                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ENDIF

                  ! Check for valid correlation value
                  IF (Output%Corr_Coeff2 .LT. CORR_MIN .OR. Output%Corr_Coeff2 .GT. 1.0) THEN

                     QC_Flag = CORRELATION_FAILURE
                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ENDIF

                  ! ----------------------------------------------------------------------
                  ! Compute vector end points for backward and forward vectors
                  ! ----------------------------------------------------------------------

                  CALL Compute_End_Points(Nav, Target_Elem_native, Element_Disp, &
                     Target_Line_native, Line_Disp,  &
                     Lat_Match, Lon_Match)

                  CALL Compute_End_Points(Nav, Target_Elem_native, Element_Disp2, &
                     Target_Line_native, Line_Disp2, &
                     Lat_Match2, Lon_Match2)

                  ! ----------------------------------------------------------------------
                  ! Compute u- and v- wind components. The backward vector is computed
                  ! first and then the forward vector is computed.
                  ! ----------------------------------------------------------------------

                  ! For GOES-NOP save scan time of match locations and update time
                  ! interval

                  !IF (Sensor_Series == SERIES_GOES_N_IMG) THEN

                    !Store in epoch 1970, not hour in day.
                    !EpochScanTimes(FIRST) = input%Line_Time_Image1(Channel, (Max_Location(2)+INT(Line_Disp)))
                    !EpochScanTimes(LAST) = input%Line_Time_Image2(Channel, (Max_Location(2)+INT(Line_Disp2)))

                    !Science:CALL Check_Scan_Times(GOESScanTimeFirst, GOESScanTimeMiddle, GOESScanTimeLast, &
                    !Science:                                      Delta_Time_minus, Delta_Time_plus)
                    !Delta_Time_minus = seconds_to_hours(EpochScanTimes(MIDDLE) - EpochScanTimes(FIRST))
                    !Delta_Time_plus  = seconds_to_hours(EpochScanTimes(LAST) - EpochScanTimes(MIDDLE))
                  !ENDIF

                  CALL Compute_UV(Nav%NavType, Lat_Target, Lon_Target, Lat_Match, Lon_Match, &
                     -Delta_Time_minus, U_Component, V_Component)

                  ! For the reverse vector change sign to account for negative time
                  ! interval. Note, however, that the sign is changed only for the
                  ! v-component. This is because decreasing longitude values are
                  ! associated with positve u-component values (McIDAS convention).

                  IF (TRIM(Nav%NavType) .EQ. 'PS' ) V_Component = -V_Component

                  CALL Compute_UV(Nav%NavType, Lat_Target, Lon_Target, Lat_Match2, Lon_Match2, &
                     Delta_Time_plus, U_Component2, V_Component2)
                  IF (TRIM(Nav%NavType) .EQ. 'PS' ) U_Component2 = -U_Component2

                  ! save data to output structure
                  Output%Lat_Match = Lat_Match
                  Output%Lon_Match = Lon_Match
                  Output%U_Component = U_Component
                  Output%V_Component = V_Component
                  Output%Lat_Match2 = Lat_Match2
                  Output%Lon_Match2 = Lon_Match2
                  Output%U_Component2 = U_Component2
                  Output%V_Component2 = V_Component2


               ! Attempt nested 5x5 tracking to improve estimate
               ELSE IF (Nested_Tracking_Flag .EQ. sym%YES) THEN

                  ! pass in clear sky flag - 1/25/16

                  CALL Nested_Tracking(Nav%NavType, FORWARD, CSWV_Flag, Box_Data,Search_BrtTemp_Values2,Delta_X_Fcst, &
                     Delta_Y_Fcst, Channel, Lag_Table2, LD_Test2, ED_Test2, TandemTriplet, & !<--- aab Tandem mod
                     Least_Square_Error, Factor, Target_Scene_Mean, &
                     Match_Scene_Mean, Target_Scene_StdDev, Match_Scene_StdDev, &
                     CC_Test2, SF_Test2, DBSCAN_Out, SetOfPoints)

                  ! aab Tandem (dbscan)
                  ! new code to compute average displacements from DBSCAN output structure
                  Sum_DX1 = SUM(DBSCAN_Out%Delta_X_Reverse(1:DBSCAN_Out%Number_CTP_Sample1))
                  Sum_DY1 = SUM(DBSCAN_Out%Delta_Y_Reverse(1:DBSCAN_Out%Number_CTP_Sample1))
                  Avg_DX1 = Sum_DX1/DBSCAN_Out%Number_CTP_Sample1
                  Avg_DY1 = Sum_DY1/DBSCAN_Out%Number_CTP_Sample1
                  ED_Test = -Avg_DX1
                  LD_Test = -Avg_DY1
                  Sum_DX2 = SUM(DBSCAN_Out%Delta_X_Forward(1:DBSCAN_Out%Number_CTP_Sample2))
                  Sum_DY2 = SUM(DBSCAN_Out%Delta_Y_Forward(1:DBSCAN_Out%Number_CTP_Sample2))
                  Avg_DX2 = Sum_DX2/DBSCAN_Out%Number_CTP_Sample2
                  Avg_DY2 = Sum_DY2/DBSCAN_Out%Number_CTP_Sample2
                  ED_Test2 = Avg_DX2
                  LD_Test2 = Avg_DY2
                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                  Output%Corr_Coeff2 = CC_Test2

                  ! Check success of tracking
                  IF (SF_Test2 .NE. sym%SUCCESS) THEN

                     QC_Flag = SF_Test2
                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ENDIF

                  ! do not use CTP for CSWV - 1/25/16
                  IF (CSWV_Flag .EQ. sym%NO) THEN

                     ! aab Tandem (dbscan) - comment out vvvvvvvvvvvvvvvv
                     !CALL Sort_Shell_NP(DBSCAN_Out%Number_CTP_Sample2, &
                     !   DBSCAN_Out%CTP_5x5_Sample2, &
                     !   Median_CTP_Sample2)

                     ! Check consistency of median CTP samples, fail if difference is 100 mb or greater.
                     !IF (ABS(Median_CTP_Sample2 - Median_CTP_Sample1) .GE. &
                     !   PRESSURE_CONSISTENCY_THRESHOLD) THEN

                     !   QC_Flag = CTP_CONSISTENCY_FAILURE
                     !   CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                     !                EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                     !                Box_Data) ! aab added Box_Data for CloudyPixels output
                     !   CALL IncrementElement(QC_Flag, Start_Element)
                     !   CYCLE

                     !ENDIF
                     ! aab Tandem (dbscan) - comment out ^^^^^^^^^^^^^^^^^^


                     ! Find the median pressure associated with the combined sample
                     ! (sample 1 is from the reverse vector while sample 2 is from
                     ! the forward vector).
                     ! aab Tandem (dbscan) vvvvvvvvvvvvvvvvvvvvv
                     !CALL Find_Combined_Median(DBSCAN_Out, Fcst_Prof, &
                     !                          Output, Combined_Median, Combined_Median_Hgt, QC_Flag) !outputs
                     if (DBSCAN_Out%Number_CTP_Sample1 > 0 .AND. DBSCAN_Out%Number_CTP_Sample2 > 0) then
                       CALL Find_Combined_Median(DBSCAN_Out, Fcst_Prof, &
                                                 Output, Combined_Median, Combined_Median_Hgt, QC_Flag) !outputs
                     else
                       QC_Flag = 17
                     endif
                     ! aab Tandem (dbscan) ^^^^^^^^^^^^^^^^^^^^^^

                     IF (QC_Flag .NE. sym%SUCCESS) THEN

                        CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                        CALL IncrementElement(QC_Flag, Start_Element)
                        CYCLE

                     ENDIF

                  ENDIF  ! CSWV check - 1/25/16

                  ! Set Inversion Flag for output file.
                  IF (LL_Inver_Flag .EQ. sym%YES) THEN

                     Output%Inversion_Flag = sym%YES

                  ELSE

                     Output%Inversion_Flag = sym%NO

                  ENDIF

                  ! Add absolute displacement to location of target box center
                  ! to get endpoint of motion vector, then return earth location
                  ! of the endpoint.

                  CALL Compute_End_Points(Nav, Target_Elem_native, ED_Test, Target_Line_native, &
                     LD_Test, Lat_Match3, Lon_Match3)

                  CALL Compute_End_Points(Nav, Target_Elem_native, ED_Test2, Target_Line_native, &
                     LD_Test2,  Lat_Match4, Lon_Match4)

                  ! For GOES-NOP save scan time of match locations and update time
                  ! interval
                  !IF (Sensor_Series == SERIES_GOES_N_IMG) THEN
                  !
                  !  EpochScanTimes(FIRST) = input%Line_Time_Image1(Channel, (Max_Location(2)+INT(LD_Test)))
                  !  EpochScanTimes(LAST) = input%Line_Time_Image2(Channel, (Max_Location(2)+INT(LD_Test2)))
                  !
                  !  !Science:CALL Check_Scan_Times(GOESScanTimeFirst, GOESScanTimeMiddle, GOESScanTimeLast, &
                  !  !Science:                                      Delta_Time_minus, Delta_Time_plus)
                  !  Delta_Time_minus = seconds_to_hours(EpochScanTimes(MIDDLE) - EpochScanTimes(FIRST))
                  !  Delta_Time_plus  = seconds_to_hours(EpochScanTimes(LAST) - EpochScanTimes(MIDDLE))
                  !ENDIF

                  ! Compute the U- and V- wind components for the nested tracking
                  ! estimate.

                  CALL Compute_UV(Nav%NavType, Lat_Target, Lon_Target, Lat_Match3, Lon_Match3, &
                     -Delta_Time_minus, U_Component, V_Component)

                  IF (TRIM(Nav%NavType) .EQ. 'PS')  V_Component = -V_Component

                  CALL Compute_UV(Nav%NavType, Lat_Target, Lon_Target, Lat_Match4, Lon_Match4, &
                     Delta_Time_plus, U_Component2, V_Component2)

                  IF (TRIM(Nav%NavType) .EQ. 'PS') U_Component2 = -U_Component2

                  Output%Lat_Match = Lat_Match3
                  Output%Lon_Match = Lon_Match3
                  Output%Lat_Match2 = Lat_Match4
                  Output%Lon_Match2 = Lon_Match4
                  Output%U_Component = U_Component
                  Output%V_Component = V_Component
                  Output%U_Component2 = U_Component2
                  Output%V_Component2 = V_Component2

               ENDIF ! Check for 5x5 tracking

               ! ----------------------------------------------------------------------
               ! Average the u- and v- components to obtain the average vector.
               ! ----------------------------------------------------------------------

               IF (Nested_Tracking_Flag .EQ. sym%NO) THEN

                  U_Avg = (U_Component + U_Component2) / 2.0
                  V_Avg = (V_Component + V_Component2) / 2.0
                  Output%U_Avg = U_Avg
                  Output%V_Avg = V_Avg

                  ! Call routine to check accelerations and apply gross error checks.
                  !
                  CALL Gross_Error_Checks(Nav%NavType, U_Component, U_Component2, V_Component, &
                  ! WCB
                  !   V_Component2, Output%U_Fcst, Output%V_Fcst, &
                     V_Component2, Output%U_Fcst, Output%V_Fcst, Output%Median_Press, Output%Cloud_Type, &
                     QC_Flag)

                  IF (QC_Flag .NE. sym%SUCCESS) THEN

                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ENDIF

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

                  ! If 5x5 tracking is successful we need to update forecast wind
                  ! before checking for gross errors because the pressure has changed.
                  ! Follow same logic as before but with combined median pressure.

                  CALL Locate(Fcst_Prof%Pressure, Press_Levels, Combined_Median, Level)

                  ! Interpolate forecast wind to new median pressure level.
                  IF ( (Level+1) .LE. Press_Levels) THEN

                     Output%U_Fcst = ((Fcst_Prof%Pressure(Level+1) - Combined_Median) * &
                        Fcst_Prof%U_Wind(Level) +  (Combined_Median - &
                        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) - Combined_Median) * &
                        Fcst_Prof%V_Wind(Level) +  (Combined_Median - &
                        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

                  ! Check 5x5 estimate for gross errors
                  !
                  CALL Gross_Error_Checks(Nav%NavType, U_Component, U_Component2, V_Component, &
                  ! WCB
                  !   V_Component2, Output%U_Fcst, Output%V_Fcst, &
                     V_Component2, Output%U_Fcst, Output%V_Fcst, Output%Median_Press, Output%Cloud_Type, &
                     QC_Flag)

                  Output%QC_Flag = QC_Flag

                  IF (QC_Flag .NE. sym%SUCCESS) THEN

                     CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
                     CALL IncrementElement(QC_Flag, Start_Element)
                     CYCLE

                  ELSE

                     ! If the estimate passes the gross checks compute average value
                     U_Avg = (U_Component + U_Component2) / 2.0
                     V_Avg = (V_Component + V_Component2) / 2.0
                     Output%U_Avg = U_Avg
                     Output%V_Avg = V_Avg

                  ENDIF

               ENDIF  ! done checking 5x5 estimate for gross errors


              ! ----------------------------------------------------------------------
              ! Count number of winds without errors.
              ! ----------------------------------------------------------------------

            ENDIF   ! max gradient check

            ! ------------------------------------------------------------------------
            ! Add target to array of wind data, then process next box
            ! ------------------------------------------------------------------------

            CALL Update_Wind_Data(Wind_Data, TargetData, Error_Counts, Target_Number, Good_Winds,& !updates
                                     EpochScanTimes, QC_Flag, Output, DBSCAN_Out, CSWV_Flag, MAX_WIND_RECORDS,&
                                     Box_Data) ! aab added Box_Data for CloudyPixels output
            CALL IncrementElement(QC_Flag, Start_Element)

         END DO



         targetSceneStartLine = targetSceneStartLine + NOMINAL_BOX_SIZE
         targetSceneEndLine = targetSceneEndLine + NOMINAL_BOX_SIZE
      END DO



      ! --------------------------------------------------------------------------
      ! Deallocate arrays.
      ! --------------------------------------------------------------------------

      !CALL Destroy_Spatial_Uniformity(Mean_Values, Max_Values, Min_Values, StdDev_Values)

      DEALLOCATE(Search_BrtTemp_Values, Search_BrtTemp_Values2)
      DEALLOCATE(Lag_Table, Lag_Table2)
      DEALLOCATE(SetOfPoints)
      DEALLOCATE(SSD_Matrix)


      ret = .TRUE.
   END FUNCTION Target_Selection_Loop





END MODULE AMV_EN_TARGET_SELECTION_M



