GetBotSeg3D Subroutine

  SUBROUTINE GetBotSeg3D( x, t, isInit )
!! Gets the Bottom segment info

    ! Get the Bottom segment info (index and range interval) for XY position, x
    ! sets Botx and Botn
    ! LP: See comment in GetTopSeg3D.

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

    nx = NbtyPts( 1 )
    ny = NbtyPts( 2 )

    IsegBotx = MIN( MAX( IsegBotx, 1 ), nx - 1 )
    IsegBoty = MIN( MAX( IsegBoty, 1 ), ny - 1 )

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

    IF ( IsegBotx ==  0 .AND. Bot(  1,  1 )%x( 1 ) == x( 1 ) ) IsegBotx = 1
    IF ( IsegBotx == nx .AND. Bot( nx,  1 )%x( 1 ) == x( 1 ) ) IsegBotx = nx - 1
    IF ( IsegBoty ==  0 .AND. Bot(  1,  1 )%x( 2 ) == x( 2 ) ) IsegBoty = 1
    IF ( IsegBoty == ny .AND. Bot(  1, ny )%x( 2 ) == x( 2 ) ) IsegBoty = ny - 1

    IF ( IsegBotx <= 0 .OR. IsegBotx >= nx .OR. IsegBoty <= 0 .OR. IsegBoty >= ny ) THEN
       WRITE( PRTFile, * ) 'Warning: GetBotSeg3D: Bottom bathymetry undefined below the ray, x', x
       IsegBotx = MIN( MAX( IsegBotx, 1 ), nx - 1 )
       IsegBoty = MIN( MAX( IsegBoty, 1 ), ny - 1 )
    END IF

    ! WRITE( PRTFile, * ) 'IsegBot', IsegBotx, IsegBoty

    xBotSeg  = [ Bot( IsegBotx, 1 )%x( 1 ), Bot( IsegBotx + 1, 1 )%x( 1 ) ]   ! segment limits in range
    yBotSeg  = [ Bot( 1, IsegBoty )%x( 2 ), Bot( 1, IsegBoty + 1 )%x( 2 ) ]   ! segment limits in range

    Botx = Bot( IsegBotx, IsegBoty )%x
    Botxmid = ( Botx + Bot( IsegBotx + 1, IsegBoty + 1 )%x ) * 0.5D0

    ! identify the normal based on the active triangle of a pair
    ! normal of triangle side pointing up and to the left
    Bot_tri_n = [ -( yBotSeg( 2 ) - yBotSeg( 1 ) ), xBotSeg( 2 ) - xBotSeg( 1 ) ]
    Bot_tri_n = Bot_tri_n / NORM2( Bot_tri_n )
    over_diag_amount = DOT_PRODUCT( x( 1 : 2 ) - Botxmid( 1 : 2 ), Bot_tri_n )
    Bot_td_onEdge = ( ABS( over_diag_amount ) < TRIDIAG_THRESH )
    IF ( .NOT. isInit .AND. Bot_td_justSteppedTo ) THEN
       Bot_td_side = Bot_td_outgoingSide
    ELSE IF ( .NOT. isInit .AND. Bot_td_onEdge ) THEN
       Bot_td_side = ( DOT_PRODUCT( t( 1 : 2 ), Bot_tri_n ) >= 0.0D0 )
    ELSE
       Bot_td_side = ( over_diag_amount >= 0.0D0 )
    END IF
    Bot_td_justSteppedTo = .FALSE.
    IF ( .NOT. Bot_td_side ) THEN
       Botn = Bot( IsegBotx, IsegBoty )%n1
    ELSE
       Botn = Bot( IsegBotx, IsegBoty )%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( Botx( 3 ) ) .OR. ANY( ISNAN( Botn ) ) ) THEN
       WRITE( PRTFile, * ) 'Error: Boundary segment contains NaN!'
       CALL ERROUT( 'BELLHOP: GetBotSeg3D', 'Boundary segment contains NaN' )
    END IF
  END SUBROUTINE GetBotSeg3D