SUBROUTINE ReduceStep3D( x0, urayt, iSegx0, iSegy0, iSegz0, h )
INTEGER, INTENT( IN ) :: iSegx0, iSegy0, iSegz0
REAL (KIND=8), INTENT( IN ) :: x0( 3 ), urayt( 3 ) ! ray coordinate and tangent
REAL (KIND=8), INTENT( INOUT ) :: h ! reduced step size
REAL (KIND=8) :: hInt, hBoxx, hBoxy, hBoxz, hTop, hBot, hxSeg, hySeg, hTopDiag, hBotDiag ! step sizes
REAL (KIND=8) :: d( 3 ), d0( 3 ), tri_n( 3 )
REAL (KIND=8) :: x( 3 ) ! ray coordinate after full trial step
REAL (KIND=8) :: xSeg( 2 ), ySeg( 2 )
REAL (KIND=8) :: over_diag_amount
LOGICAL :: newSide, newOnEdge
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'ReduceStep3D'
! Detect SSP interface or boundary crossing and reduce step, if necessary, to land on that crossing.
! Keep in mind possibility that user put source right on an interface
! and that multiple events can occur (crossing interface, top, and bottom in a single step).
x = x0 + h * urayt ! make a trial step
! interface crossing in depth
! Step reduction is not done for the top or bottom layer
! Instead the SSP is extrapolated
! This prevents problems when the boundaries are outside the domain of the SSP
hInt = huge( hInt )
IF ( ABS( urayt( 3 ) ) > EPSILON( hInt ) ) THEN
IF ( SSP%z( iSegz0 ) > x( 3 ) .AND. iSegz0 > 1 ) THEN
hInt = ( SSP%z( iSegz0 ) - x0( 3 ) ) / urayt( 3 )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Shallower bound SSP Z > z; hInt', SSP%z( iSegz0 ), x( 3 ), hInt
ELSE IF ( SSP%z( iSegz0 + 1 ) < x( 3 ) .AND. iSegz0 + 1 < SSP%Nz ) THEN
hInt = ( SSP%z( iSegz0 + 1 ) - x0( 3 ) ) / urayt( 3 )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Deeper bound SSP Z < z; hInt', SSP%z( iSegz0 + 1 ), x( 3 ), hInt
END IF
END IF
! ray mask using a box centered at ( xs_3D( 1 ), xs_3D( 2 ), 0 )
hBoxx = huge( hBoxx )
IF ( ABS( x( 1 ) - xs_3D( 1 ) ) > Beam%Box%x ) THEN
hBoxx = ( Beam%Box%x - ABS( ( x0( 1 ) - xs_3D( 1 ) ) ) ) / ABS( urayt( 1 ) )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Beam box crossing X, hBoxx', Beam%Box%x, hBoxx
END IF
hBoxy = huge( hBoxy )
IF ( ABS( x( 2 ) - xs_3D( 2 ) ) > Beam%Box%y ) THEN
hBoxy = ( Beam%Box%y - ABS( ( x0( 2 ) - xs_3D( 2 ) ) ) ) / ABS( urayt( 2 ) )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Beam box crossing Y, hBoxy', Beam%Box%y, hBoxy
END IF
hBoxz = huge( hBoxz )
IF ( ABS( x( 3 ) ) > Beam%Box%z ) THEN
hBoxz = ( Beam%Box%z - ABS( x0( 3 ) ) ) / ABS( urayt( 3 ) )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Beam box crossing Z, hBoxz', Beam%Box%z, hBoxz
END IF
! top crossing
hTop = huge( hTop )
d = x - Topx ! vector from top to ray
! LP: Changed from EPSILON( hTop ) to 0 in conjunction with StepToBdry3D change here.
IF ( DOT_PRODUCT( Topn, d ) >= 0.0D0 ) THEN
d0 = x0 - Topx ! vector from top node to ray origin
hTop = -DOT_PRODUCT( d0, Topn ) / DOT_PRODUCT( urayt, Topn )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Top crossing hTop', hTop
END IF
! bottom crossing
hBot = huge( hBot )
d = x - Botx ! vector from bottom to ray
! LP: Changed from EPSILON( h1 ) (should have been h3!) to 0 in conjunction
! with StepToBdry3D change here.
IF ( DOT_PRODUCT( Botn, d ) >= 0.0D0 ) THEN
d0 = x0 - Botx ! vector from bottom node to ray origin
hBot = -DOT_PRODUCT( d0, Botn ) / DOT_PRODUCT( urayt, Botn )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Bottom crossing hBot', hBot
END IF
! top/bottom/ocean segment crossing in x
hxSeg = huge( hxSeg )
xSeg( 1 ) = MAX( xTopSeg( 1 ), xBotSeg( 1 ) )
xSeg( 2 ) = MIN( xTopSeg( 2 ), xBotSeg( 2 ) )
IF ( SSP%Type == 'H' ) THEN ! ocean segment
xSeg( 1 ) = MAX( xSeg( 1 ), SSP%Seg%x( iSegx0 ) )
xSeg( 2 ) = MIN( xSeg( 2 ), SSP%Seg%x( iSegx0 + 1 ) )
END IF
IF ( ABS( urayt( 1 ) ) > EPSILON( hxSeg ) ) THEN
IF ( x( 1 ) < xSeg( 1 ) ) THEN
hxSeg = -( x0( 1 ) - xSeg( 1 ) ) / urayt( 1 )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Min bound SSP X > x; hxSeg', xSeg( 1 ), x( 1 ), hxSeg
ELSE IF ( x( 1 ) > xSeg( 2 ) ) THEN
hxSeg = -( x0( 1 ) - xSeg( 2 ) ) / urayt( 1 )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Max bound SSP X < x; hxSeg', xSeg( 2 ), x( 1 ), hxSeg
END IF
END IF
! top/bottom/ocean segment crossing in y
hySeg = huge( hySeg )
ySeg( 1 ) = MAX( yTopSeg( 1 ), yBotSeg( 1 ) )
ySeg( 2 ) = MIN( yTopSeg( 2 ), yBotSeg( 2 ) )
IF ( SSP%Type == 'H' ) THEN ! ocean segment
ySeg( 1 ) = MAX( ySeg( 1 ), SSP%Seg%y( iSegy0 ) )
ySeg( 2 ) = MIN( ySeg( 2 ), SSP%Seg%y( iSegy0 + 1 ) )
END IF
! LP: removed 1000 * epsilon which mbp had comment questioning also
IF ( ABS( urayt( 2 ) ) > EPSILON( hySeg ) ) THEN
IF ( x( 2 ) < ySeg( 1 ) ) THEN
hySeg = -( x0( 2 ) - ySeg( 1 ) ) / urayt( 2 )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Min bound SSP Y > y; hySeg', ySeg( 1 ), x( 2 ), hySeg
ELSE IF ( x( 2 ) > ySeg( 2 ) ) THEN
hySeg = -( x0( 2 ) - ySeg( 2 ) ) / urayt( 2 )
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Max bound SSP Y < y; hySeg', ySeg( 2 ), x( 2 ), hySeg
END IF
END IF
! triangle crossing within a top segment
hTopDiag = huge( hTopDiag )
IF ( .NOT. Top_td_onEdge ) THEN
d = x - Topxmid ! vector from top center to ray end
d0 = x0 - Topxmid ! vector from top center to ray origin
tri_n = [ -( yTopSeg( 2 ) - yTopSeg( 1 ) ), xTopSeg( 2 ) - xTopSeg( 1 ), 0.0d0 ]
tri_n = tri_n / NORM2( tri_n )
over_diag_amount = DOT_PRODUCT( d, tri_n )
newSide = ( over_diag_amount >= 0.0D0 )
IF ( newSide .NEQV. Top_td_side ) THEN
newOnEdge = ( ABS( over_diag_amount ) < TRIDIAG_THRESH )
IF ( newOnEdge ) THEN
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'ReduceStep3D Top naturally stepped to tri diag edge, h over_diag_amount', &
over_diag_amount
ELSE
hTopDiag = -DOT_PRODUCT( d0, tri_n ) / DOT_PRODUCT( urayt, tri_n )
IF ( hTopDiag < 0.0D0 ) THEN
WRITE( PRTFile, * ) 'ReduceStep3D Top BHC_WARN_TRIDIAG_H_NEGATIVE'
hTopDiag = 0.0D0
END IF
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Top tri diag crossing hTopDiag dot(n, d0) dot(n, d)', &
hTopDiag, DOT_PRODUCT( tri_n, d0 ), over_diag_amount
END IF
END IF
END IF
! triangle crossing within a bottom segment
hBotDiag = huge( hBotDiag )
IF ( .NOT. Bot_td_onEdge ) THEN
d = x - Botxmid ! vector from bottom center to ray end
d0 = x0 - Botxmid ! vector from bottom center to ray origin
tri_n = [ -( yBotSeg( 2 ) - yBotSeg( 1 ) ), xBotSeg( 2 ) - xBotSeg( 1 ), 0.0d0 ]
tri_n = tri_n / NORM2( tri_n )
over_diag_amount = DOT_PRODUCT( d, tri_n )
newSide = ( over_diag_amount >= 0.0D0 )
IF ( newSide .NEQV. Bot_td_side ) THEN
newOnEdge = ( ABS( over_diag_amount ) < TRIDIAG_THRESH )
IF ( newOnEdge ) THEN
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'ReduceStep3D Bot naturally stepped to tri diag edge, h over_diag_amount', &
over_diag_amount
ELSE
hBotDiag = -DOT_PRODUCT( d0, tri_n ) / DOT_PRODUCT( urayt, tri_n )
IF ( hBotDiag < 0.0D0 ) THEN
WRITE( PRTFile, * ) 'ReduceStep3D Bot BHC_WARN_TRIDIAG_H_NEGATIVE'
hBotDiag = 0.0D0
END IF
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'Bot tri diag crossing hBotDiag dot(n, d0) dot(n, d)', &
hBotDiag, DOT_PRODUCT( tri_n, d0 ), over_diag_amount
END IF
END IF
END IF
h = MIN( h, hInt, hBoxx, hBoxy, hBoxz, hTop, hBot, hxSeg, hySeg, hTopDiag, hBotDiag ) ! take limit set by shortest distance to a crossing
IF ( h < -1d-4 ) THEN
WRITE( PRTFile, * ) 'ReduceStep3D: negative h'
! CALL ERROUT( 'ReduceStep3D', 'negative h' )
END IF
IF ( h < INFINITESIMAL_STEP_SIZE * Beam%deltas ) THEN ! is it taking an infinitesimal step?
h = INFINITESIMAL_STEP_SIZE * Beam%deltas ! make sure we make some motion
iSmallStepCtr = iSmallStepCtr + 1 ! keep a count of the number of sequential small steps
ELSE
iSmallStepCtr = 0 ! didn't do a small step so reset the counter
END IF
END SUBROUTINE ReduceStep3D
