!!$/=============================================================================/
!!$Copyright (c) 2007, The University of Massachusetts Dartmouth 
!!$Produced at the School of Marine Science & Technology 
!!$Marine Ecosystem Dynamics Modeling group
!!$All rights reserved.
!!$
!!$The FVCOM Offline Lagrangian Model has been developed by the joint UMASSD-WHOI
!!$research team.   For details of authorship and attribution of credit please see
!!$the FVCOM technical manual or contact the MEDM group.
!!$
!!$ 
!!$This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu/ The
!!$full copyright notice is contained in the file COPYRIGHT located in the root
!!$directory of the FVCOM code. This original header must be maintained in all
!!$distributed versions.
!!$
!!$THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
!!$ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
!!$IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
!!$ARE  DISCLAIMED.  
!!$
!!$/-----------------------------------------------------------------------------/
!!$CVS VERSION INFORMATION
!!$$Id: $
!!$$Name: $
!!$$Revision: $
!!$/=============================================================================/
SUBROUTINE TRIANGLE_GRID_EDGE

  !==============================================================================!
  !  This program is used to define the non-overlapped, unstructured             !
  !  triangular meshes used for flux computations. The mesh could be             !
  !  created using the commerical software called "sms8.0" or other              !
  !  mesh generation programs. The mesh file generated by sms8.0 can             !
  !  be directly used for this subroutine, while the mesh file                   !
  !  generated using other programs must be coverted the data format             !
  !  to meet the required format used here.                                      !
  !==============================================================================!
  !     variable list:							         !
  !  vx(m)    :: vx(i) = x-coordinate of node i (input from mesh)	         !
  !  vy(m)    :: vy(i) = y-coordinate of node i (input from mesh)                !
  !  nv(n,3)  :: nv(i:1-3) = 3 node numbers of element i                         !
  !  xc(n)    :: xc(i) = x-coordinate of element i (calculated from vx)          !
  !  yc(n)    :: yc(i) = y-coordinate of element i (calculated from vy)          !
  !                                                                              !
  !  nbe(n,3) :: nbe(i,1->3) = element index of 1->3 neighbors of element i      !
  !  isbce(n) :: flag if element is on the boundary, see below for values        !
  !  isonb(m) :: flag is node is on the boundary, see below for values           !
  !                                                                              !
  !  ntve(m)  :: the number of neighboring elements of node m                    !
  !  nbve(m,ntve(m)) :: nbve(i,1->ntve(i)) = ntve elements containing node i     !
  !  nbvt(m,ntve(m)) :: nbvt(i,j) = the node number of node i in element         ! 
  !                     nbve(i,j) (has a value of 1,2,or 3)                      ! 
  !                                                                              !
  !==============================================================================!
  !     classification of the triangles nodes, and edges                         !
  !                                                                              !
  !     isonb(i)=0:  node in the interior computational domain                   !
  !     isonb(i)=1:  node on the solid boundary                                  !
  !     isonb(i)=2:  node on the open boundary                                   !
  !                                                                              !
  !     isbce(i)=0:  element in the interior computational domain                !
  !     isbce(i)=1:  element on the solid boundary                               !
  !     isbce(i)=2:  element on the open boundary                                !
  !     isbce(i)=3:  element with 2 solid boundary edges                         !
  !                                                                              !
  !==============================================================================!

  USE ALL_VARS
  IMPLICIT NONE

  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: TEMP,NB_TMP,CELLS,NBET
  INTEGER, ALLOCATABLE, DIMENSION(:)   :: CELLCNT
  INTEGER                              :: I,J,II,JJ,NTMP,NCNT,NFLAG,JJB
  INTEGER                              :: N1,N2,N3,J1,J2,J3


  !==============================================================================|
  !   SET UP MESH (HORIZONTAL COORDINATES)                                       |
  !==============================================================================|

  !--------------CALCULATE GLOBAL MINIMUMS AND MAXIMUMS--------------------------!

  VXMIN = MINVAL(VX(1:M)) ; VXMAX = MAXVAL(VX(1:M))
  VYMIN = MINVAL(VY(1:M)) ; VYMAX = MAXVAL(VY(1:M))
  
  VXMIN = 0.0; VXMAX = 0.0
  VYMIN = 0.0; VYMAX = 0.0
  !--------------SHIFT GRID TO UPPER RIGHT CARTESIAN-----------------------------!

  VX = VX - VXMIN
  VY = VY - VYMIN

  !--------------CALCULATE GLOBAL ELEMENT CENTER GRID COORDINATES----------------!

  DO I=1,N   
     XC(I)  = (VX(NV(I,1)) + VX(NV(I,2)) + VX(NV(I,3)))/3.0_SP
     YC(I)  = (VY(NV(I,1)) + VY(NV(I,2)) + VY(NV(I,3)))/3.0_SP
  END DO

  XC(0) = 0.0_SP ; YC(0) = 0.0_SP


  !----------------------------INITIALIZE----------------------------------------!

  ISBCE = 0
  ISONB = 0
  NBE   = 0

  !
  !----DETERMINE NBE(i=1:n,j=1:3): INDEX OF 1 to 3 NEIGHBORING ELEMENTS----------!
  !
  ALLOCATE(NBET(N,3)) ; NBET = 0
  ALLOCATE(CELLS(M,50)) ; CELLS = 0
  ALLOCATE(CELLCNT(M))  ; CELLCNT = 0
  DO I=1,N
     N1 = NV(I,1) ; CELLCNT(N1) = CELLCNT(N1)+1
     N2 = NV(I,2) ; CELLCNT(N2) = CELLCNT(N2)+1
     N3 = NV(I,3) ; CELLCNT(N3) = CELLCNT(N3)+1
     CELLS(NV(I,1),CELLCNT(N1)) = I
     CELLS(NV(I,2),CELLCNT(N2)) = I
     CELLS(NV(I,3),CELLCNT(N3)) = I
  END DO
  if(maxval(cellcnt) > 50)write(*,*)'bad',maxval(cellcnt)
  DO I=1,N
     N1 = NV(I,1)
     N2 = NV(I,2)
     N3 = NV(I,3)
     DO J1 = 1,CELLCNT(N1) 
        DO J2 = 1,CELLCNT(N2) 
           IF((CELLS(N1,J1) == CELLS(N2,J2)).AND. CELLS(N1,J1) /= I)NBE(I,3) = CELLS(N1,J1)
        END DO
     END DO
     DO J2 = 1,CELLCNT(N2) 
        DO J3 = 1,CELLCNT(N3) 
           IF((CELLS(N2,J2) == CELLS(N3,J3)).AND. CELLS(N2,J2) /= I)NBE(I,1) = CELLS(N2,J2)
        END DO
     END DO
     DO J1 = 1,CELLCNT(N1) 
        DO J3 = 1,CELLCNT(N3) 
           IF((CELLS(N1,J1) == CELLS(N3,J3)).AND. CELLS(N1,J1) /= I)NBE(I,2) = CELLS(N3,J3)
        END DO
     END DO
  END DO
  DEALLOCATE(CELLS,CELLCNT)
  !   IF(MSR)WRITE(IPT,*)  '!  NEIGHBOR FINDING      :    COMPLETE'
  !
  !--ENSURE ALL ELEMENTS HAVE AT LEAST ONE NEIGHBOR------------------------------!
  !
  NFLAG = 0
  DO I=1,N
     IF(SUM(NBE(I,1:3))==0)THEN 
        NFLAG = 1
        WRITE(*,*)'ELEMENT ',I,' AT ',XC(I),YC(I),' HAS NO NEIGHBORS'
        STOP
     END IF
  END DO
  IF(NFLAG == 1) STOP

  !
  !----IF ELEMENT ON BOUNDARY SET ISBCE(I)=1 AND ISONB(J)=1 FOR BOUNDARY NODES J-!
  !

  DO I=1,N 
     IF(MIN(NBE(I,1),NBE(I,2),NBE(I,3))==0)THEN    !!ELEMENT ON BOUNDARY
        ISBCE(I) = 1
        IF(NBE(I,1) == 0)THEN 
           ISONB(NV(I,2)) = 1 ; ISONB(NV(I,3)) = 1
        END IF
        IF(NBE(I,2) ==0) THEN
           ISONB(NV(I,1)) = 1 ; ISONB(NV(I,3)) = 1
        END IF
        IF(NBE(I,3) ==0) THEN
           ISONB(NV(I,1)) = 1 ; ISONB(NV(I,2)) = 1
        END IF
     END IF
  END DO

  !==============================================================================|
  !             DEFINE NTVE, NBVE, NBVT                                          !
  !                                                                              !
  ! ntve(1:m):           total number of the surrounding triangles               !
  !                      connected to the given node                             !
  ! nbve(1:m, 1:ntve+1): the identification number of surrounding                !
  !                      triangles with a common node (counted clockwise)        !
  ! nbvt(1:m,ntve(1:m)): the idenfication number of a given node over            !
  !                      each individual surrounding triangle(counted            !
  !                      clockwise)                                              !
  !==============================================================================|

  !
  !----DETERMINE MAX NUMBER OF SURROUNDING ELEMENTS------------------------------!
  !
  MX_NBR_ELEM = 0
  DO I=1,M
     NCNT = 0
     DO J=1,N
        IF( FLOAT(NV(J,1)-I)*FLOAT(NV(J,2)-I)*FLOAT(NV(J,3)-I) == 0.0_SP) &
             NCNT = NCNT + 1
     END DO
     MX_NBR_ELEM = MAX(MX_NBR_ELEM,NCNT)
  END DO

  !
  !----ALLOCATE ARRAYS BASED ON MX_NBR_ELEM--------------------------------------!
  ! 
  ALLOCATE(NBVE(M,MX_NBR_ELEM+1))
  ALLOCATE(NBVT(M,MX_NBR_ELEM+1))
  !
  !--DETERMINE NUMBER OF SURROUNDING ELEMENTS FOR NODE I = NTVE(I)---------------!
  !--DETERMINE NBVE - INDICES OF NEIGHBORING ELEMENTS OF NODE I------------------!
  !--DETERMINE NBVT - INDEX (1,2, or 3) OF NODE I IN NEIGHBORING ELEMENT---------!
  !
  DO I=1,M
     NCNT=0
     DO J=1,N
        IF (FLOAT(NV(J,1)-I)*FLOAT(NV(J,2)-I)*FLOAT(NV(J,3)-I) == 0.0_SP)THEN
           NCNT = NCNT+1
           NBVE(I,NCNT)=J
           IF((NV(J,1)-I) == 0) NBVT(I,NCNT)=1
           IF((NV(J,2)-I) == 0) NBVT(I,NCNT)=2
           IF((NV(J,3)-I) == 0) NBVT(I,NCNT)=3
        END IF
     ENDDO
     NTVE(I)=NCNT
  ENDDO

  ALLOCATE(NB_TMP(M,MX_NBR_ELEM+1))
  DO I=1,M
     IF(ISONB(I) == 0) THEN
        NB_TMP(1,1)=NBVE(I,1)
        NB_TMP(1,2)=NBVT(I,1)
        DO J=2,NTVE(I)+1
           II=NB_TMP(J-1,1)
           JJ=NB_TMP(J-1,2)
           NB_TMP(J,1)=NBE(II,JJ+1-INT((JJ+1)/4)*3)
           JJ=NB_TMP(J,1)
           IF((NV(JJ,1)-I) == 0) NB_TMP(J,2)=1
           IF((NV(JJ,2)-I) == 0) NB_TMP(J,2)=2
           IF((NV(JJ,3)-I) == 0) NB_TMP(J,2)=3
        ENDDO

        DO J=2,NTVE(I)+1
           NBVE(I,J)=NB_TMP(J,1)
        ENDDO

        DO J=2,NTVE(I)+1
           NBVT(I,J)=NB_TMP(J,2)
        ENDDO

        NTMP=NTVE(I)+1
        IF(NBVE(I,1) /= NBVE(I,NTMP)) THEN
           PRINT*, I,'NBVE(I) NOT CORRECT!!'
           STOP
        ENDIF
        IF(NBVT(I,1) /= NBVT(I,NTMP)) THEN
           PRINT*, I,'NBVT(I) NOT CORRECT!!'
           STOP
        END IF

     ELSE 
        JJB=0

        DO J=1,NTVE(I)
           JJ=NBVT(I,J)
           IF(NBE(NBVE(I,J),JJ+2-INT((JJ+2)/4)*3) == 0) THEN
              JJB=JJB+1
              NB_TMP(JJB,1)=NBVE(I,J)
              NB_TMP(JJB,2)=NBVT(I,J)
           END IF
        ENDDO

        IF(JJB /= 1) THEN
           PRINT*, 'ERROR IN ISONB !,I,J', I,J
           PAUSE
        END IF

        DO J=2,NTVE(I)
           II=NB_TMP(J-1,1)
           JJ=NB_TMP(J-1,2)
           NB_TMP(J,1)=NBE(II,JJ+1-INT((JJ+1)/4)*3)
           JJ=NB_TMP(J,1)
           IF((NV(JJ,1)-I) == 0) NB_TMP(J,2)=1
           IF((NV(JJ,2)-I) == 0) NB_TMP(J,2)=2
           IF((NV(JJ,3)-I) == 0) NB_TMP(J,2)=3
        ENDDO

        DO J=1,NTVE(I)
           NBVE(I,J)=NB_TMP(J,1)
           NBVT(I,J)=NB_TMP(J,2)
        ENDDO
        NBVE(I,NTVE(I)+1)=0

     END IF
  END DO
  DEALLOCATE(NB_TMP)

  RETURN
END SUBROUTINE TRIANGLE_GRID_EDGE

   !==============================================================================|
