GetTopSeg3D Subroutine

  SUBROUTINE GetTopSeg3D( x, t, isInit )
!! Gets the Top segment info

    ! Get the Top segment info (index and range interval) for XY position, x
    ! sets Topx and Topn
    ! LP: See discussion of changes in README.md.

    INTEGER, PARAMETER :: PRTFile = 6
    REAL (KIND=8), INTENT( IN ) :: x( 3 ), t( 3 )
    LOGICAL,       INTENT( IN ) :: isInit
    INTEGER :: nx, ny
    REAL (KIND=8) :: Top_tri_n( 2 )   ! triangle normals
    REAL (KIND=8) :: over_diag_amount

    nx = NatiPts( 1 )
    ny = NatiPts( 2 )

    IsegTopx = MIN( MAX( IsegTopx, 1 ), nx - 1 )
    IsegTopy = MIN( MAX( IsegTopy, 1 ), ny - 1 )

    IF ( t( 1 ) >= 0.0 ) THEN
       DO WHILE ( IsegTopx >= 1 .AND. Top( IsegTopx, 1 )%x( 1 ) > x( 1 ) )
          IsegTopx = IsegTopx - 1
       END DO
       DO WHILE ( IsegTopx >= 1 .AND. IsegTopx < nx .AND. Top( IsegTopx + 1, 1 )%x( 1 ) <= x( 1 ) )
          IsegTopx = IsegTopx + 1
       END DO
    ELSE
      DO WHILE ( IsegTopx < nx .AND. Top( IsegTopx + 1, 1 )%x( 1 ) < x( 1 ) )
         IsegTopx = IsegTopx + 1
      END DO
      DO WHILE ( IsegTopx >= 1 .AND. IsegTopx < nx .AND. Top( IsegTopx, 1 )%x( 1 ) >= x( 1 ) )
         IsegTopx = IsegTopx - 1
      END DO
    END IF
    IF ( t( 2 ) >= 0.0 ) THEN
      DO WHILE ( IsegTopy >= 1 .AND. Top( 1, IsegTopy )%x( 2 ) > x( 2 ) )
         IsegTopy = IsegTopy - 1
      END DO
      DO WHILE ( IsegTopy >= 1 .AND. IsegTopy < ny .AND. Top( 1, IsegTopy + 1 )%x( 2 ) <= x( 2 ) )
         IsegTopy = IsegTopy + 1
      END DO
    ELSE
      DO WHILE ( IsegTopy < ny .AND. Top( 1, IsegTopy + 1 )%x( 2 ) < x( 2 ) )
         IsegTopy = IsegTopy + 1
      END DO
      DO WHILE ( IsegTopy >= 1 .AND. IsegTopy < ny .AND. Top( 1, IsegTopy )%x( 2 ) >= x( 2 ) )
         IsegTopy = IsegTopy - 1
      END DO
    END IF

    IF ( IsegTopx ==  0 .AND. Top(  1,  1 )%x( 1 ) == x( 1 ) ) IsegTopx = 1
    IF ( IsegTopx == nx .AND. Top( nx,  1 )%x( 1 ) == x( 1 ) ) IsegTopx = nx - 1
    IF ( IsegTopy ==  0 .AND. Top(  1,  1 )%x( 2 ) == x( 2 ) ) IsegTopy = 1
    IF ( IsegTopy == ny .AND. Top(  1, ny )%x( 2 ) == x( 2 ) ) IsegTopy = ny - 1

    IF ( IsegTopx <= 0 .OR. IsegTopx >= nx .OR. IsegTopy <= 0 .OR. IsegTopy >= ny ) THEN
       WRITE( PRTFile, * ) 'Warning: GetTopSeg3D: Top altimetry undefined above the ray, x', x
       IsegTopx = MIN( MAX( IsegTopx, 1 ), nx - 1 )
       IsegTopy = MIN( MAX( IsegTopy, 1 ), ny - 1 )
    END IF


    xTopSeg  = [ Top( IsegTopx, 1 )%x( 1 ), Top( IsegTopx + 1, 1 )%x( 1 ) ]   ! segment limits in range
    yTopSeg  = [ Top( 1, IsegTopy )%x( 2 ), Top( 1, IsegTopy + 1 )%x( 2 ) ]   ! segment limits in range

    Topx = Top( IsegTopx, IsegTopy )%x
    Topxmid = ( Topx + Top( IsegTopx + 1, IsegTopy + 1 )%x ) * 0.5D0

    ! WRITE( PRTFile, * ) 'IsegTop', IsegTopx, IsegTopy
    ! WRITE( PRTFile, * ) 'Topx x', Topx, x

    ! identify the normal based on the active triangle of a pair
    ! normal of triangle side pointing up and to the left
    Top_tri_n = [ -( yTopSeg( 2 ) - yTopSeg( 1 ) ), xTopSeg( 2 ) - xTopSeg( 1 ) ]
    Top_tri_n = Top_tri_n / NORM2( Top_tri_n )
    over_diag_amount = DOT_PRODUCT( x( 1 : 2 ) - Topxmid( 1 : 2 ), Top_tri_n )
    Top_td_onEdge = ( ABS( over_diag_amount ) < TRIDIAG_THRESH )
    ! WRITE( PRTFile, * ) 'Top_tri_n over_diag_amount', Top_tri_n, over_diag_amount
    IF ( .NOT. isInit .AND. Top_td_justSteppedTo ) THEN
       Top_td_side = Top_td_outgoingSide
    ELSE IF ( .NOT. isInit .AND. Top_td_onEdge ) THEN
       Top_td_side = ( DOT_PRODUCT( t( 1 : 2 ), Top_tri_n ) >= 0.0D0 )
    ELSE
       Top_td_side = ( over_diag_amount >= 0.0D0 )
    END IF
    Top_td_justSteppedTo = .FALSE.
    IF ( .NOT. Top_td_side ) THEN
       Topn = Top( IsegTopx, IsegTopy )%n1
    ELSE
       Topn = Top( IsegTopx, IsegTopy )%n2
    END IF

    ! if the Bot depth is bad (a NaN) then error out
    ! LP: Originally set segment to invalid and relied on later code to catch this.
    IF ( ISNAN( Topx( 3 ) ) .OR. ANY( ISNAN( Topn ) ) ) THEN
       WRITE( PRTFile, * ) 'Error: Boundary segment contains NaN!'
       CALL ERROUT( 'BELLHOP: GetTopSeg3D', 'Boundary segment contains NaN' )
    END IF
  END SUBROUTINE GetTopSeg3D