! ------------------------------------------------------------------------------ ! ! 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) ! ! 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 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) !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, & 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, & 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 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) 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 ALLOCATE(SetOfPoints((Box_Size-4)**2, 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, & 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 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. ! -------------------------------------------------------------------- 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, & Least_Square_Error, Factor, Target_Scene_Mean, & Match_Scene_Mean, Target_Scene_StdDev, Match_Scene_StdDev, & CC_Test2, SF_Test2, DBSCAN_Out, SetOfPoints) 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 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 ! Find the median pressure associated with the combined sample ! (sample 1 is from the reverse vector while sample 2 is from ! the forward vector). CALL Find_Combined_Median(DBSCAN_Out, Fcst_Prof, & Output, Combined_Median, Combined_Median_Hgt, QC_Flag) !outputs 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