!********************************************************************************
!> @file
!> NAME:         AMV_EN_DBSCAN_M
!>
!> FUNCTION:     DBSCAN routines
!>
!> DESCRIPTION:
!>
!>     Density based spatial clustering variables/routines.
!>
!> REFERENCE:
!
!> CALLING SEQUENCE:
!>
!>     USE AMV_EN_DBSCAN_M
!>
!> DEPENDENCIES:
!>
!> RESTRICTIONS:
!>
!> HISTORY:
!> - 01/2017 - Brian Helgans - Created.
!>
!> @ingroup WINDS_AMV_EN
!********************************************************************************
MODULE AMV_EN_DBSCAN_M

   USE TYPE_KINDS_AIT
   USE GEOCAT_CONSTANTS
   USE fw_log_mod

   IMPLICIT NONE

   PRIVATE

   PUBLIC :: DBSCAN_Output, DBSCAN_Output_Initialize

   !> @brief define DBSCAN_Output structure
   TYPE :: DBSCAN_Output

      INTEGER(LONG) :: Number_CTP_Sample1
      INTEGER(LONG) :: Number_CTP_Sample2
      !INTEGER(LONG) :: Number_CTP_Sample
      INTEGER(LONG) :: Sample1
      INTEGER(LONG) :: Sample2
      INTEGER(LONG) :: Initial_Sample
      INTEGER(LONG) :: Number_BT_Sample1
      INTEGER(LONG) :: Number_BT_Sample2
      !INTEGER(LONG) :: Number_BT_Sample
      INTEGER(LONG) :: NumberOfClusters1
      INTEGER(LONG) :: NumberOfClusters2
      INTEGER(LONG) :: MaxCluster1
      INTEGER(LONG) :: MaxCluster2
      INTEGER(LONG) :: MaxClusterSize1
      INTEGER(LONG) :: MaxClusterSize2
      REAL(SINGLE) :: Max1_PctOfSample
      REAL(SINGLE) :: Max2_PctOfSample
      REAL(SINGLE) :: Avg_Delta_X_Sample1
      REAL(SINGLE) :: Avg_Delta_Y_Sample1
      REAL(SINGLE) :: Avg_Delta_X_Sample2
      REAL(SINGLE) :: Avg_Delta_Y_Sample2
      ! add 2nd and 3rd cluster averages - wcb 9/7/12
      REAL(SINGLE) :: Avg_Delta_X_2ndCluster_Rev
      REAL(SINGLE) :: Avg_Delta_Y_2ndCluster_Rev
      REAL(SINGLE) :: Avg_Delta_X_2ndCluster_For
      REAL(SINGLE) :: Avg_Delta_Y_2ndCluster_For
      REAL(SINGLE) :: Avg_Delta_X_3rdCluster_Rev
      REAL(SINGLE) :: Avg_Delta_Y_3rdCluster_Rev
      REAL(SINGLE) :: Avg_Delta_X_3rdCluster_For
      REAL(SINGLE) :: Avg_Delta_Y_3rdCluster_For
      REAL(SINGLE) :: Median_Press_2ndCluster
      REAL(SINGLE) :: Median_Press_3rdCluster
      INTEGER(LONG) :: PointsIn2ndCluster
      INTEGER(LONG) :: PointsIn3rdCluster
      REAL(SINGLE) :: StdDev_MVD1
      REAL(SINGLE) :: StdDev_MVD2
      REAL(SINGLE) :: PctOfAvg1
      REAL(SINGLE) :: PctOfAvg2
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CTP_5x5_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CTP_5x5_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CTP_5x5_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CORR_5x5_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CORR_5x5_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CORR_5x5_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CldTemp_5x5_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CldTemp_5x5_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: CldTemp_5x5_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Cld_Hgt_5x5_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Cld_Hgt_5x5_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Cld_Hgt_5x5_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Press_Error_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Press_Error_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Press_Error_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Temp_Error_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Temp_Error_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Temp_Error_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: BT_5x5_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: BT_5x5_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: BT_5x5_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Cld_OD_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Cld_OD_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Cld_OD_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: ACHA_Cost_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: ACHA_Cost_Sample2
      !REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: ACHA_Cost_Sample
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Phase_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Phase_Sample2
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Type_Sample1
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Type_Sample2


      ! save displacements for the reverse and forward vectors
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Delta_X_Reverse
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Delta_Y_Reverse
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Delta_X_Forward
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: Delta_Y_Forward

      ! save 2D fields for plotting
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:,:) :: Delta_X_Init
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:,:) :: Delta_Y_Init
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:,:) :: Delta_X_Final
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:,:) :: Delta_Y_Final
      ! add array to hold sub-vector number (1=reverse time step, 2=forward time step)
      REAL(SINGLE), ALLOCATABLE, DIMENSION(:) :: vector
      CONTAINS
      PROCEDURE :: Create
      FINAL :: Destroy
   END TYPE DBSCAN_Output

   CONTAINS



   !> @brief allocates the structure
   FUNCTION Create(DBSCAN_Out, Box_Size) RESULT(res)
      CHARACTER(LEN=*), PARAMETER :: FUNC = "Create_DBSCAN"
      CLASS(DBSCAN_Output) :: DBSCAN_Out !< structure to be created
      INTEGER(LONG), INTENT(IN) :: Box_Size !< size of the box.
      LOGICAL :: res
      INTEGER(LONG) astatus

      ! aab Tandem (dbscan)   vvvvvvvvvvvvvvvvv
      integer comb_size

      !ALLOCATE(DBSCAN_Out%CTP_5x5_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%CTP_5x5_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%CORR_5x5_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%CORR_5x5_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%CldTemp_5x5_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%CldTemp_5x5_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%Cld_Hgt_5x5_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%Cld_Hgt_5x5_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%BT_5x5_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%BT_5x5_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%Press_Error_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%Press_Error_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%Temp_Error_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%Temp_Error_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%Cld_OD_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%Cld_OD_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%ACHA_Cost_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%ACHA_Cost_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%Phase_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%Phase_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%Type_Sample1((Box_Size-4)**2), &
      !   DBSCAN_Out%Type_Sample2((Box_Size-4)**2), &
      !   DBSCAN_Out%Delta_X_Reverse((Box_Size-4)**2), &
      !   DBSCAN_Out%Delta_Y_Reverse((Box_Size-4)**2), &
      !   DBSCAN_Out%Delta_X_Forward((Box_Size-4)**2), &
      !   DBSCAN_Out%Delta_Y_Forward((Box_Size-4)**2), &
      !   DBSCAN_Out%Delta_X_Init(Box_Size, Box_Size), &
      !   DBSCAN_Out%Delta_Y_Init(Box_Size, Box_Size), &
      !   DBSCAN_Out%Delta_X_Final(Box_Size, Box_Size), &
      !   DBSCAN_Out%Delta_Y_Final(Box_Size, Box_Size), stat=astatus)

      ! increase by a factor of 2 to hold samples from both time steps
      comb_size = 2*((Box_Size-4)**2) 
      ALLOCATE(DBSCAN_Out%CTP_5x5_Sample1(comb_size), &
         DBSCAN_Out%CTP_5x5_Sample2(comb_size), &
         DBSCAN_Out%CORR_5x5_Sample1(comb_size), &
         DBSCAN_Out%CORR_5x5_Sample2(comb_size), &
         DBSCAN_Out%CldTemp_5x5_Sample1(comb_size), &
         DBSCAN_Out%CldTemp_5x5_Sample2(comb_size), &
         DBSCAN_Out%Cld_Hgt_5x5_Sample1(comb_size), &
         DBSCAN_Out%Cld_Hgt_5x5_Sample2(comb_size), &
         DBSCAN_Out%BT_5x5_Sample1(comb_size), &
         DBSCAN_Out%BT_5x5_Sample2(comb_size), &
         DBSCAN_Out%Press_Error_Sample1(comb_size), &
         DBSCAN_Out%Press_Error_Sample2(comb_size), &
         DBSCAN_Out%Temp_Error_Sample1(comb_size), &
         DBSCAN_Out%Temp_Error_Sample2(comb_size), &
         DBSCAN_Out%Cld_OD_Sample1(comb_size), &
         DBSCAN_Out%Cld_OD_Sample2(comb_size), &
         DBSCAN_Out%ACHA_Cost_Sample1(comb_size), &
         DBSCAN_Out%ACHA_Cost_Sample2(comb_size), &
         DBSCAN_Out%Phase_Sample1(comb_size), &
         DBSCAN_Out%Phase_Sample2(comb_size), &
         DBSCAN_Out%Type_Sample1(comb_size), &
         DBSCAN_Out%Type_Sample2(comb_size), &
         DBSCAN_Out%Delta_X_Reverse(comb_size), &
         DBSCAN_Out%Delta_Y_Reverse(comb_size), &
         DBSCAN_Out%Delta_X_Forward(comb_size), &
         DBSCAN_Out%Delta_Y_Forward(comb_size), &
         DBSCAN_Out%vector(comb_size), &
         DBSCAN_Out%Delta_X_Init(2*Box_Size, 2*Box_Size), &
         DBSCAN_Out%Delta_Y_Init(2*Box_Size, 2*Box_Size), &
         DBSCAN_Out%Delta_X_Final(2*Box_Size, 2*Box_Size), &
         DBSCAN_Out%Delta_Y_Final(2*Box_Size, 2*Box_Size), stat=astatus)
      ! aab Tandem (dbscan)   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

      IF (astatus /= 0) THEN
         CALL fw_log_error(FUNC, 'failure to allocate')
         res = .FALSE.
      ELSE
         res = .TRUE.
      ENDIF

   END FUNCTION Create


   !> @brief Initialize variable DBSCAN_Out to default values
   SUBROUTINE DBSCAN_Output_Initialize(DBSCAN_Out)
      TYPE(DBSCAN_Output), INTENT(INOUT) :: DBSCAN_Out !< to be initialized

      DBSCAN_Out%Number_CTP_Sample1 = MISSING_VALUE_INT4
      DBSCAN_Out%Number_CTP_Sample2 = MISSING_VALUE_INT4
      DBSCAN_Out%Sample1 = MISSING_VALUE_INT4
      DBSCAN_Out%Sample2 = MISSING_VALUE_INT4
      DBSCAN_Out%Initial_Sample = 0
      DBSCAN_Out%Number_BT_Sample1 = MISSING_VALUE_INT4
      DBSCAN_Out%Number_BT_Sample2 = MISSING_VALUE_INT4
      DBSCAN_Out%NumberOfClusters1 = MISSING_VALUE_INT4
      DBSCAN_Out%NumberOfClusters2 = MISSING_VALUE_INT4
      DBSCAN_Out%MaxCluster1 = MISSING_VALUE_INT4
      DBSCAN_Out%MaxCluster2 = MISSING_VALUE_INT4
      DBSCAN_Out%MaxClusterSize1 = MISSING_VALUE_INT4
      DBSCAN_Out%MaxClusterSize2  = MISSING_VALUE_INT4
      DBSCAN_Out%Max1_PctOfSample  = MISSING_VALUE_REAL4
      DBSCAN_Out%Max2_PctOfSample  = MISSING_VALUE_REAL4
      DBSCAN_Out%StdDev_MVD1 = MISSING_VALUE_REAL4
      DBSCAN_Out%StdDev_MVD2 = MISSING_VALUE_REAL4
      DBSCAN_Out%PctOfAvg1 =MISSING_VALUE_REAL4
      DBSCAN_Out%PctOfAvg2 =MISSING_VALUE_REAL4
      DBSCAN_Out%CTP_5x5_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%CTP_5x5_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%CORR_5x5_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%CORR_5x5_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%CldTemp_5x5_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%CldTemp_5x5_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Cld_Hgt_5x5_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Cld_Hgt_5x5_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%BT_5x5_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%BT_5x5_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Press_Error_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Press_Error_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Temp_Error_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Temp_Error_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Cld_OD_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Cld_OD_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%ACHA_Cost_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%ACHA_Cost_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Phase_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Phase_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Type_Sample1(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Type_Sample2(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Delta_X_Reverse(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Delta_Y_Reverse(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Delta_X_Forward(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Delta_Y_Forward(:) = MISSING_VALUE_REAL4
      ! aab Tandem (dbscan)
      DBSCAN_Out%vector(:) = MISSING_VALUE_REAL4
      DBSCAN_Out%Delta_X_Init = 0.0
      DBSCAN_Out%Delta_Y_Init = 0.0
      DBSCAN_Out%Delta_X_Final = 0.0
      DBSCAN_Out%Delta_Y_Final = 0.0
      DBSCAN_Out%Avg_Delta_X_Sample1 = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_Y_Sample1 = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_X_Sample2 = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_Y_Sample2 = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_X_2ndCluster_Rev = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_Y_2ndCluster_Rev = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_X_2ndCluster_For = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_Y_2ndCluster_For = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_X_3rdCluster_Rev = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_Y_3rdCluster_Rev = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_X_3rdCluster_For = MISSING_VALUE_REAL4
      DBSCAN_Out%Avg_Delta_Y_3rdCluster_For = MISSING_VALUE_REAL4
      DBSCAN_Out%Median_Press_2ndCluster = MISSING_VALUE_REAL4
      DBSCAN_Out%Median_Press_3rdCluster = MISSING_VALUE_REAL4
      DBSCAN_Out%PointsIn2ndCluster = MISSING_VALUE_INT4
      DBSCAN_Out%PointsIn3rdCluster = MISSING_VALUE_INT4



   END SUBROUTINE DBSCAN_Output_Initialize

   !> @brief deallocate space for variable DBSCAN_Out
   !>      Should be called automatically as the structures FINAL method.
   SUBROUTINE Destroy(DBSCAN_Out)
      TYPE(DBSCAN_Output) :: DBSCAN_Out !< to be destroyed.

      IF(ALLOCATED(DBSCAN_Out%CTP_5x5_Sample1)) THEN
         DEALLOCATE(DBSCAN_Out%CTP_5x5_Sample1, &
         DBSCAN_Out%CTP_5x5_Sample2, &
         DBSCAN_Out%CORR_5x5_Sample1, &
         DBSCAN_Out%CORR_5x5_Sample2, &
         DBSCAN_Out%CldTemp_5x5_Sample1, &
         DBSCAN_Out%CldTemp_5x5_Sample2, &
         DBSCAN_Out%Cld_Hgt_5x5_Sample1, &
         DBSCAN_Out%Cld_Hgt_5x5_Sample2, &
         DBSCAN_Out%BT_5x5_Sample1, &
         DBSCAN_Out%BT_5x5_Sample2, &
         DBSCAN_Out%Press_Error_Sample1, &
         DBSCAN_Out%Press_Error_Sample2, &
         DBSCAN_Out%Temp_Error_Sample1, &
         DBSCAN_Out%Temp_Error_Sample2, &
         DBSCAN_Out%Cld_OD_Sample1, &
         DBSCAN_Out%Cld_OD_Sample2, &
         DBSCAN_Out%ACHA_Cost_Sample1, &
         DBSCAN_Out%ACHA_Cost_Sample2, &
         DBSCAN_Out%Phase_Sample1, &
         DBSCAN_Out%Phase_Sample2, &
         DBSCAN_Out%Type_Sample1, &
         DBSCAN_Out%Type_Sample2, &
         DBSCAN_Out%Delta_X_Reverse, &
         DBSCAN_Out%Delta_Y_Reverse, &
         DBSCAN_Out%Delta_X_Forward, &
         DBSCAN_Out%Delta_Y_Forward, &
         ! aab Tandem (dbscan)
         DBSCAN_Out%vector, &
         DBSCAN_Out%Delta_X_Init, &
         DBSCAN_Out%Delta_Y_Init, &
         DBSCAN_Out%Delta_X_Final, &
         DBSCAN_Out%Delta_Y_Final)
      ENDIF

   END SUBROUTINE Destroy

END MODULE AMV_EN_DBSCAN_M
