SUBROUTINE ComputeBdryTangentNormal( Bdry, BotTop )
! Does some pre-processing on the boundary points to pre-compute segment
! lengths (%Len),
! tangents (%t, %nodet),
! normals (%n, %noden), and
! curvatures (%kappa)
!
! The boundary is also extended with a constant depth to infinity to cover cases where the ray
! exits the domain defined by the user
INTEGER, SAVE :: NPts = 0
REAL (KIND=8), ALLOCATABLE :: phi( : )
REAL (KIND=8) :: sss
TYPE(BdryPt), INTENT( INOUT ) :: Bdry( : )
CHARACTER (LEN=3), INTENT( IN ) :: BotTop ! Flag indicating bottom or top reflection
CHARACTER (LEN=2), SAVE :: CurvilinearFlag = '-'
SELECT CASE ( BotTop )
CASE ( 'Bot' )
NPts = NbtyPts
CurvilinearFlag = btyType
CASE ( 'Top' )
NPts = NatiPts
CurvilinearFlag = atiType
END SELECT
! extend the bathymetry to +/- infinity in a piecewise constant fashion
Bdry( 1 )%x( 1 ) = -sqrt( huge( Bdry( 1 )%x( 1 ) ) ) / 1.0d5
Bdry( 1 )%x( 2 ) = Bdry( 2 )%x( 2 )
Bdry( 1 )%HS = Bdry( 2 )%HS
Bdry( NPts )%x( 1 ) = +sqrt( huge( Bdry( 1 )%x( 1 ) ) ) / 1.0d5
Bdry( NPts )%x( 2 ) = Bdry( NPts - 1 )%x( 2 )
Bdry( NPts )%HS = Bdry( NPts - 1 )%HS
! compute tangent and outward-pointing normal to each bottom segment
! tBdry( 1, : ) = xBdry( 1, 2:NPts ) - xBdry( 1, 1:NPts - 1 )
! tBdry( 2, : ) = xBdry( 2, 2:NPts ) - xBdry( 2, 1:NPts - 1 )
! above caused compiler problems
BoundaryPt: DO ii = 1, NPts - 1
Bdry( ii )%t = Bdry( ii + 1 )%x - Bdry( ii )%x
Bdry( ii )%Dx = Bdry( ii )%t( 2 ) / Bdry( ii )%t( 1 ) ! first derivative
! write( *, * ) 'Dx, t', Bdry( ii )%Dx, Bdry( ii )%x, 1 / ( Bdry( ii )%x( 2 ) / 500 )
! normalize the tangent vector
Bdry( ii )%Len = NORM2( Bdry( ii )%t )
Bdry( ii )%t = Bdry( ii )%t / Bdry( ii )%Len
SELECT CASE ( BotTop )
CASE ( 'Bot' )
Bdry( ii )%n( 1 ) = -Bdry( ii )%t( 2 )
Bdry( ii )%n( 2 ) = +Bdry( ii )%t( 1 )
CASE ( 'Top' )
Bdry( ii )%n( 1 ) = +Bdry( ii )%t( 2 )
Bdry( ii )%n( 2 ) = -Bdry( ii )%t( 1 )
END SELECT
END DO BoundaryPt
IF ( CurvilinearFlag( 1 : 1 ) == 'C' ) THEN ! curvilinear option: compute tangent and normal at node by averaging normals on adjacent segments
! averaging two centered differences is equivalent to forming a single centered difference of two steps ...
DO ii = 2, NPts - 1
sss = Bdry( ii - 1 )%Len / ( Bdry( ii - 1 )%Len + Bdry( ii )%Len )
sss = 0.5 ! LP: BUG? Line above is overwritten.
Bdry( ii )%Nodet = ( 1.0 - sss ) * Bdry( ii - 1 )%t + sss * Bdry( ii )%t
END DO
Bdry( 1 )%Nodet = [ 1.0, 0.0 ] ! tangent left-end node
Bdry( NPts )%Nodet = [ 1.0, 0.0 ] ! tangent right-end node
SELECT CASE ( BotTop )
CASE ( 'Bot' )
Bdry( : )%Noden( 1 ) = -Bdry( : )%Nodet( 2 )
Bdry( : )%Noden( 2 ) = +Bdry( : )%Nodet( 1 )
CASE ( 'Top' )
Bdry( : )%Noden( 1 ) = +Bdry( : )%Nodet( 2 )
Bdry( : )%Noden( 2 ) = -Bdry( : )%Nodet( 1 )
END SELECT
! compute curvature in each segment
! LP: This allocation is not necessary, could just have two variables for
! current and next phi. Operating on the whole array can trigger compiler
! SIMD parallelism (AVX-512 etc.), but this is unlikely to happen for
! atan2, and this is in one-time setup code anyway.
ALLOCATE( phi( NPts ), Stat = IAllocStat )
phi = atan2( Bdry( : )%Nodet( 2 ), Bdry( : )%Nodet( 1 ) ) ! this is the angle at each node
DO ii = 1, NPts - 1
Bdry( ii )%kappa = ( phi( ii + 1 ) - phi( ii ) ) / Bdry( ii )%Len ! this is curvature = dphi/ds
Bdry( ii )%Dxx = ( Bdry( ii + 1 )%Dx - Bdry( ii )%Dx ) / & ! second derivative
( Bdry( ii + 1 )%x( 1 ) - Bdry( ii )%x( 1 ) )
Bdry( ii )%Dss = Bdry( ii )%Dxx * Bdry( ii )%t( 1 ) ** 3 ! derivative in direction of tangent
!write( *, * ) 'kappa, Dss, Dxx', Bdry( ii )%kappa, Bdry( ii )%Dss, Bdry( ii )%Dxx, &
! 1 / ( ( 8 / 1000 ** 2 ) * ABS( Bdry( ii )%x( 2 ) ) ** 3 ), Bdry( ii )%x( 2 )
! -1 / ( 4 * ( Bdry( ii )%x( 2 ) ) ** 3 / 1000000 ), Bdry( ii )%x( 2 )
Bdry( ii )%kappa = Bdry( ii )%Dss !over-ride kappa !!!!!
END DO
DEALLOCATE( phi ) ! LP: was missing deallocation
ELSE
Bdry%kappa = 0
END IF
END SUBROUTINE ComputeBdryTangentNormal
