!! Provides BELLHOP3D program definition

PROGRAM BELLHOP3D
!! Gaussian beam tracing for ocean acoustics in three dimensions

  ! Gaussian beam tracing in three dimensions
  ! Michael B. Porter

  ! Copyright (C) 2009 Michael B. Porter

  ! This program is free software: you can redistribute it and/or modify
  ! it under the terms of the GNU General Public License as published by
  ! the Free Software Foundation, either version 3 of the License, or
  ! (at your option) any later version.

  ! This program is distributed in the hope that it will be useful,
  ! but WITHOUT ANY WARRANTY; without even the implied warranty of
  ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  ! GNU General Public License for more details.

  ! You should have received a copy of the GNU General Public License
  ! along with this program.  If not, see <http://www.gnu.org/licenses/>.

  ! Original version done at NRL in 1986
  ! Published at SACLANT Undersea Research Center 1989
  !
  ! Converted to Fortran 2003 and greatly enhanced for outside users in 2010
  ! with the support of the Agency for Defense Development, Republic of Korea
  ! Supported also by the U.S. Office of Naval Research
  !
  ! Changes included:
  !    Use of standard BELLHOP input files
  !    Inclusion of multiple sources and receivers
  !    Integrated option for Nx2D or full 3D
  !    Implementation of several standard beam options for both the Nx2D and 3D cases
  !    Reading in:
  !       oceanography from a SSP file
  !       bathymetry or altimetry  from a BTY or ATI file
  !       a reflection coefficient from a BRC or TRC file
  !       a source beam pattern from a SBP file

  ! Loose ends:
  !    Nx2D version does not handle jumps in cx, cy (routine step2d)
  !    cannot specify isingle( 2 ) for alpha and beta
  !    Trilinear interpolation (hex) ignores cxy values.
  !    If the water depth for the lower half space is much larger than that for the SSP
  !    in the water column, it will use a step that is too large and exit the ray box.

  !    Cerveny beams (rarely used):
  !       influenceC calls SSP; need to select EvaluateSSP2D or EvaluateSSP3D for that to work in BELLHOP3D
  !       efficiency changes for Cerveny beams as well (and in BELLHOP)
  !       space filling or minimum width options are applied to both alpha and beta--- often only want space filling in azimuth

  !    Influend3DGeoHat writes no eigenray info if number of receiver ranges NR=1
  !    r( 1 ) = 1 m in BELLHOP plus logic for reversing
  !    geogaussian should pre-calculate dtauds and dqds like geohat?
  !    fix automatic deltas selection
  !    detect Sz below bottom with bty file

  !    Should probably use a Structure of Arrays instead of an Array of Structures for speed

  !    Desired additional features:
  !       Terrain following option for receiver depth
  !       Variable bottom type vs.lat/long

  USE bellhopMod
  USE ReadEnvironmentBell
  USE RefCoef
  USE bdry3DMod
  USE BeamPattern
  USE sspMod
  USE influence
  USE Influence3D
  USE FatalError
  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT

  IMPLICIT NONE

  CHARACTER ( LEN=80 ) :: FileRoot

  ThreeD = .TRUE.

  ! get the file root for naming all input and output files
  ! should add some checks here ...

  CALL GET_COMMAND_ARGUMENT( 1, FileRoot )
  ! Open the print file
  OPEN( UNIT = PRTFile, FILE = TRIM( FileRoot ) // '.prt', STATUS = 'UNKNOWN', IOSTAT = iostat )

  ! Read in control data

  CALL ReadEnvironment(    FileRoot, ThreeD )
  CALL ReadATI3D( FileRoot, Bdry%Top%HS%Opt( 5 : 5 ), Bdry%Top%HS%Depth, PRTFile )    ! AlTImetry
  CALL ReadBTY3D( FileRoot, Bdry%Bot%HS%Opt( 2 : 2 ), Bdry%Bot%HS%Depth, PRTFile )    ! BaThYmetry

  CALL ReadReflectionCoefficient(  FileRoot, Bdry%Bot%HS%Opt( 1 : 1 ), Bdry%Top%HS%Opt( 2 : 2 ), PRTFile )    ! (top and bottom)
  SBPFlag = Beam%RunType( 3 : 3 )
  CALL ReadPAT( FileRoot,                                 PRTFile )    ! Source Beam Pattern
  CALL OpenOutputFiles( FileRoot, ThreeD )

  CALL BellhopCore

  CONTAINS

! **********************************************************************!

SUBROUTINE BellhopCore
  !! Core subroutine to run Bellhop algorithm

  USE angleMod
  USE SourceReceiverPositions
  USE ArrMod
  USE WriteRay

  ! @note "LP"
  ! Was 400M; this translates to 16 GB of arrivals, which is half my computer's
  ! memory, so I can't run this and the C++/CUDA version in parallel. This large
  ! a size should not be necessary for most runs anyway, and it's not great for
  ! a program to allocate this kind of memory just as a default (rather than in
  ! response to the user specifying a very large problem size). So, cut in half.
  ! @endnote
  INTEGER,   PARAMETER :: ArrivalsStorage = 200000000
  INTEGER              :: IBPvec( 1 ), ibp, iBeamWindow2, irz, itheta, isx, isy, isz, iRec, ir
  REAL      ( KIND=8 ) :: Tstart, Tstop
  REAL      ( KIND=8 ) :: Amp0, RadMax, S
  REAL      ( KIND=8 ) :: c0, cimag0, gradc( 3 ), cxx, cyy, czz, cxy, cxz, cyz, rho
  REAL      ( KIND=8 ) :: tdummy( 3 )
  REAL      ( KIND=8 ), ALLOCATABLE :: x_rcvrMat( :, :, : ), t_rcvr( :, : )
  COMPLEX   ( KIND=8 ) :: epsilon( 2 )
  COMPLEX, ALLOCATABLE :: P( :, :, : ), U( :, : )

  CALL CPU_TIME( Tstart )

  omega = 2.0 * pi * freq

  IF ( Beam%deltas == 0.0 ) Beam%deltas = ( Bdry%Bot%HS%Depth - Bdry%Top%HS%Depth ) / 10.0   ! Automatic step size selection

  Angles%alpha  = DegRad * Angles%alpha   ! convert to radians
  Angles%Dalpha = 0.0
  IF ( Angles%Nalpha /= 1 ) THEN
     Angles%Dalpha = ( Angles%alpha( Angles%Nalpha ) - Angles%alpha( 1 ) ) / ( Angles%Nalpha - 1 )
     ! angular spacing between beams
  END IF
  SELECT CASE ( Beam%RunType( 5 : 5 ) )
  CASE ( 'I' )
     NRz_per_range = 1         ! irregular grid
  CASE ( 'R' )
     NRz_per_range = Pos%NRz   ! rectilinear grid
  CASE DEFAULT
     WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown grid RunType: ', Beam%RunType( 5 : 5 )
     CALL ERROUT( 'BELLHOP3D', 'Unknown grid RunType' )
  END SELECT

  ! for a TL calculation, allocate space for the pressure matrix
  SELECT CASE ( Beam%RunType( 1 : 1 ) )
  CASE ( 'C', 'S', 'I' )        ! TL calculation
     ALLOCATE ( P( Pos%Ntheta, NRz_per_range, Pos%NRr ), Stat = IAllocStat )
     ALLOCATE ( U( NRz_per_range, Pos%NRr ), Stat = IAllocStat )    ! used for 2D option
     IF ( IAllocStat /= 0 ) THEN
          CALL ERROUT( 'BELLHOP', 'Insufficient memory for TL matrix: reduce Nr * NRz'  )
     END IF
  CASE ( 'A', 'a', 'R', 'E' )   ! Arrivals calculation
     ALLOCATE ( P( 1, 1, 1 ), Stat = IAllocStat )   ! open a dummy variable
     ALLOCATE ( U( 1, 1 ),    Stat = IAllocStat )   ! open a dummy variable
  CASE DEFAULT
     WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown RunType for allocation: ', Beam%RunType( 1 : 1 )
     CALL ERROUT( 'BELLHOP3D', 'Unknown RunType in allocation' )
  END SELECT

  ! for an arrivals run, allocate space for arrivals matrices
  SELECT CASE ( Beam%RunType( 1 : 1 ) )
  CASE ( 'A', 'a' )
     MaxNArr = MAX( ArrivalsStorage / ( Pos%Ntheta * NRz_per_range * Pos%NRr ), 10 )   ! allow space for at least 10 arrivals
     WRITE( PRTFile, * )
     WRITE( PRTFile, * ) '( Maximum # of arrivals = ', MaxNArr, ')'

     ALLOCATE ( Arr3D( Pos%Ntheta, NRz_per_range, Pos%NRr, MaxNArr ), &
               NArr3D( Pos%Ntheta, NRz_per_range, Pos%NRr ), Stat = IAllocStat )
     IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', &
          'Insufficient memory to allocate arrivals matrix; reduce parameter ArrivalsStorage' )

     ! For a 2D Arrivals run, also need to allocate a structure for that
     IF ( Beam%RunType( 6 : 6 ) == '2' ) THEN
         ALLOCATE ( Arr( NRz_per_range, Pos%NRr, MaxNArr ), &
                   NArr( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
         IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', &
            'Insufficient memory to allocate arrivals matrix; reduce parameter ArrivalsStorage' )
         NArr( 1 : NRz_per_range, 1 : Pos%NRr ) = 0
     END IF
  CASE DEFAULT
     MaxNArr = 1
     ALLOCATE ( Arr3D( Pos%Ntheta, NRz_per_range, Pos%NRr, 1 ), NArr3D( Pos%Ntheta, NRz_per_range, Pos%NRr ), Stat = IAllocStat )
  END SELECT

  ALLOCATE( x_rcvrMat( 2, Pos%Ntheta, Pos%NRr ), t_rcvr( 2, Pos%Ntheta ), Stat = IAllocStat )
  IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', &
       'Insufficient memory to x_rcvrMat; reduce Nr or Ntheta' )

  ! tangent along receiver bearing line
  t_rcvr( 1, : ) = COS( DegRad * Pos%theta( 1 : Pos%Ntheta ) )
  t_rcvr( 2, : ) = SIN( DegRad * Pos%theta( 1 : Pos%Ntheta ) )

  WRITE( PRTFile, * )

  Source_z: DO isz = 1, Pos%NSz         ! loop over source z-coordinate

     Source_x: DO isx = 1, Pos%Nsx      ! loop over source x-coordinate

        Source_y: DO isy = 1, Pos%Nsy   ! loop over source y-coordinate
           P = 0.0 !  Zero out field matrix
           U = 0.0
           NArr3D = 0 ! LP: Zero out 3D arrival matrix

           ! IF ( r( 1 ) == 0.0 ) r( 1 ) = 1.0
           xs_3D = [ Pos%sx( isx ), Pos%sy( isy ), Pos%sz( isz ) ]
           WRITE( PRTFile, * )
           WRITE( PRTFile, "( 'xs = ', G11.3, 2X, G11.3, 2X, G11.3 )" ) xs_3D

           ! positions of rcvrs in the x-y plane; this is pre-calculated for InfluenceGeoHatCart
           ! It is not clear that the pre-calculation saves time ...
           DO ir = 1, Pos%NRr
              DO itheta = 1, Pos%Ntheta
                 x_rcvrMat( 1 : 2, itheta, ir ) = xs_3D( 1 : 2 ) + Pos%Rr( ir ) * t_rcvr( :, itheta )
                     ! x-y coordinate of the receiver
              END DO
           END DO

           ! *** Compute 'optimal' beam constant ***

           tdummy = [ 0.0, 0.0, 1.0 ]
           CALL EvaluateSSP3D( xs_3D, tdummy, c0, cimag0, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, 'TAB' )
           ray2D( 1 )%c = c0
           ray3D( 1 )%c = c0
           CALL PickEpsilon( Beam%Type( 1 : 2 ), omega, c0, Angles%Dalpha, Angles%Dbeta, &
                             Beam%rLoop, Beam%epsMultiplier, epsilon ) ! beam constant

           ! *** Trace successive beams ***

           AzimuthalAngle: DO ibeta = 1, Angles%Nbeta ! this is also the receiver bearing angle for a 2D run
              SrcAzimAngle = RadDeg * Angles%beta( ibeta )           ! take-off azimuthal   angle in degrees
              !if ( ibeta /= 134 ) cycle AzimuthalAngle
              IF ( Angles%iSingle_beta == 0 .OR. ibeta == Angles%iSingle_beta ) THEN    ! Single beam run?
              !IF ( ibeta == 2 ) THEN    ! Single beam run?
              !IF ( mod( ibeta+1, 2 ) == 0 ) THEN    ! Single beam run?
                WRITE( PRTFile, FMT = "( 'Tracing azimuthal beam ', I4, F10.2 )" ) ibeta, SrcAzimAngle
                FLUSH( PRTFile )
                 ! WRITE( *,       FMT = "( 'Tracing beam ', I4, F10.2 )" ) ibeta, RadDeg * Angles%beta( ibeta )

                 DeclinationAngle: DO ialpha = 1, Angles%Nalpha
                    SrcDeclAngle = RadDeg *Angles%alpha( ialpha )          ! take-off declination angle in degrees

                    IF ( Angles%iSingle_alpha == 0 .OR. ialpha == Angles%iSingle_alpha ) THEN    ! Single beam run?
                    !IF ( ialpha  == 11 .or. ialpha == 12 ) THEN    ! Single beam run?

                       IF ( Angles%iSingle_alpha > 0 .OR. Angles%iSingle_beta > 0 .OR. &
                           Angles%Nalpha * Angles%Nbeta <= 20000 .OR. MOD( ialpha, 20 ) == 1 ) THEN
                           WRITE( PRTFile, FMT = "( '   Tracing declination beam ', I4, F10.2 )" ) ialpha, SrcDeclAngle
                       END IF
                       !flush( prtfile )

                       IBPvec = maxloc( SrcBmPat( :, 1 ), mask = SrcBmPat( :, 1 ) < SrcDeclAngle )
                           ! index of ray angle in beam pattern
                       IBP    = IBPvec( 1 )
                       IBP    = MAX( IBP, 1 )               ! don't go before beginning of table
                       IBP    = MIN( IBP, NSBPPts - 1 )     ! don't go past end of table

                       ! linear interpolation to get amplitude
                       s    = ( SrcDeclAngle  - SrcBmPat( IBP, 1 ) ) / ( SrcBmPat( IBP + 1, 1 ) - SrcBmPat( IBP, 1 ) )
                       Amp0 = ( 1             -                    s ) * SrcBmPat( IBP, 2 ) + s * SrcBmPat( IBP + 1, 2 )

                       ! Lloyd mirror pattern for semi-coherent option
                       IF ( Beam%RunType( 1 : 1 ) == 'S' ) THEN
                          Amp0 = Amp0 * SQRT( 2.0 ) * ABS( SIN( omega / c0 * xs_3D( 3 ) * SIN( Angles%alpha( ialpha ) ) ) )
                       END IF

                       SELECT CASE ( Beam%RunType( 6 : 6 ) )   ! flag for 2D or 3D calculation
                       CASE ( '2' )   ! Nx2D calculation, neglecting horizontal refraction
                          CALL TraceRay2D(  Angles%alpha( ialpha ), Angles%beta( ibeta ), Amp0 )
                          IF ( Beam%RunType( 1 : 1 ) /= 'R' ) THEN     ! If not a ray trace run, calculate the field
                             SELECT CASE ( Beam%Type( 1 : 1 ) )
                             CASE ( 'R' )
                                iBeamWindow2 = Beam%iBeamWindow ** 2
                                RadMax       = 50 * c0 / freq  ! 50 wavelength max radius
                                CALL InfluenceCervenyRayCen(   U, epsilon( 1 ), Angles%alpha( ialpha ), IBeamWindow2, RadMax )
                             CASE ( 'C' )
                                CALL ERROUT( 'BELLHOP3D', 'Run Type ''C'' not supported at this time' )
                                iBeamWindow2 = Beam%iBeamWindow ** 2
                                RadMax       = 50 * c0 / freq  ! 50 wavelength max radius
                                CALL InfluenceCervenyCart(     U, epsilon( 1 ), Angles%alpha( ialpha ), IBeamWindow2, RadMax )
                             CASE ( 'g' )
                                CALL InfluenceGeoHatRayCen(    U,       Angles%alpha( ialpha ), Angles%Dalpha )
                             CASE ( 'S' )
                                RadMax       = 50 * c0 / freq  ! LP: Was missing / uninitialized.
                                CALL InfluenceSGB(             U,       Angles%alpha( ialpha ), Angles%Dalpha, RadMax )
                             CASE ( 'B' )
                                CALL InfluenceGeoGaussianCart( U,       Angles%alpha( ialpha ), Angles%Dalpha )
                             CASE ( 'G', '^', ' ' )
                                CALL InfluenceGeoHatCart(      U,       Angles%alpha( ialpha ), Angles%Dalpha )
                             CASE DEFAULT
                                CALL ERROUT( 'BELLHOP3D', 'Invalid Run Type' )
                             END SELECT
                          END IF

                       CASE ( '3' )   ! full 3D calculation
                          CALL TraceRay3D( Angles%alpha( ialpha ), Angles%beta( ibeta ), epsilon, Amp0 )

                          IF ( Beam%RunType( 1 : 1 ) /= 'R' ) THEN     ! If not a ray trace run, calculate the field

                             SELECT CASE ( Beam%RunType( 2 : 2 ) )
                             CASE ( 'C' )   ! Cerveny style beams
                                CALL ERROUT( 'BELLHOP3D', 'Run Type ''C'' not supported at this time' )

                                ! option: assemble f, g, h from p-q
                                !ray3D( 1 : Beam%Nsteps )%f    =    ray3D( 1 : Beam%Nsteps )%p_tilde( 1 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_hat(   2 ) - &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%p_hat(   1 )
                                !ray3D( 1 : Beam%Nsteps )%g    = -( ray3D( 1 : Beam%Nsteps )%p_tilde( 2 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_hat(   1 ) - &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_tilde( 1 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%p_hat(   2 ) )
                                !ray3D( 1 : Beam%Nsteps )%h    =    ray3D( 1 : Beam%Nsteps )%p_tilde( 2 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_hat(   2 ) - &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%p_hat(   2 )
                                !ray3D( 1 : Beam%Nsteps )%DetP =    ray3D( 1 : Beam%Nsteps )%p_tilde( 1 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%p_hat(   2 ) - &
                                !                                   ray3D( 1 : Beam%Nsteps )%p_tilde( 2 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%p_hat(   1 )
                                !ray3D( 1 : Beam%Nsteps )%DetQ =    ray3D( 1 : Beam%Nsteps )%q_tilde( 1 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_hat(   2 ) - &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) * &
                                !                                   ray3D( 1 : Beam%Nsteps )%q_hat(   1 )

                                !CALL Influence3D_3D( xs_3D, Angles%alpha( ialpha ), iBeamWindow, P )
                             CASE ( 'g' )        ! Geometric-beams with hat-shape
                                CALL Influence3DGeoHatRayCen(       Angles%alpha( ialpha ), Angles%beta( ibeta ), &
                                                              Angles%Dalpha, Angles%Dbeta, P )
                             CASE ( 'G', '^', ' ' )   ! Geometric-beams with hat-shape in Cartesian coordinates
                                CALL Influence3DGeoHatCart(         Angles%alpha( ialpha ), Angles%beta( ibeta ), &
                                                              Angles%Dalpha, Angles%Dbeta, P, x_rcvrMat, t_rcvr )
                             CASE ( 'b' )        ! Geometric-beams with Gaussian-shape
                                CALL Influence3DGeoGaussianRayCen(  Angles%alpha( ialpha ), Angles%beta( ibeta ), &
                                                              Angles%Dalpha, Angles%Dbeta, P )
                             CASE ( 'B' )        ! Geometric-beams with Gaussian-shape in Cartesian coordiantes
                                CALL Influence3DGeoGaussianCart(    Angles%alpha( ialpha ), Angles%beta( ibeta ), &
                                                              Angles%Dalpha, Angles%Dbeta, P, x_rcvrMat, t_rcvr )
                            CASE DEFAULT
                                CALL ERROUT( 'BELLHOP3D', 'Invalid Run Type' )
                             END SELECT
                          END IF
                       CASE DEFAULT
                          WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown 2D/3D flag: ', Beam%RunType( 6 : 6 )
                          CALL ERROUT( 'BELLHOP3D', 'Unknown 2D/3D flag' )
                       END SELECT

                       ! Optionally dump rays to a disk file
                       IF ( Beam%RunType( 1 : 1 ) == 'R' ) THEN
                          CALL WriteRay3D( Angles%alpha( ialpha ), Angles%beta( ibeta ), Beam%Nsteps )
                       END IF
                    END IF   ! closes iSingle test
                 END DO DeclinationAngle

                 ! for a 2D TL run, scale the pressure and copy the 2D slice into the radial of the 3D field
                 IF ( Beam%RunType( 6 : 6 ) == '2' ) THEN  ! 2D calculation
                    SELECT CASE ( Beam%RunType( 1 : 1 ) )
                    CASE ( 'C', 'S', 'I' )   ! TL calculation
                       CALL ScalePressure( Angles%Dalpha, ray2D( 1 )%c, Pos%Rr, U, NRz_per_range, Pos%NRr, Beam%RunType, freq )
                       P( ibeta, :, : ) = U
                       U = 0   ! clear out the pressure field on the radial in prep for next radial
                    CASE ( 'A' )             ! arrivals calculation, ascii
                       NArr3D( ibeta, :, :    ) = NArr( :, : )
                       Arr3D(  ibeta, :, :, : ) = Arr(  :, :, : )
                       Arr3D(  ibeta, :, :, : )%SrcAzimAngle  = SNGL( SrcAzimAngle )   ! angle
                       Arr3D(  ibeta, :, :, : )%RcvrAzimAngle = SNGL( SrcAzimAngle )
                           ! angle (rcvr angle is same as source angle)
                       Narr = 0   ! this clears out the 2D arrival structure
                    CASE ( 'a' )             ! arrivals calculation, binary
                       NArr3D( ibeta, :, :    ) = NArr( :, : )
                       Arr3D(  ibeta, :, :, : ) = Arr(  :, :, : )
                       Arr3D(  ibeta, :, :, : )%SrcAzimAngle  = SNGL( SrcAzimAngle )   ! angle
                       Arr3D(  ibeta, :, :, : )%RcvrAzimAngle = SNGL( SrcAzimAngle )
                           ! angle (rcvr angle is same as source angle)
                       Narr = 0   ! this clears out the 2D arrival structure
                    CASE ( 'R', 'E' )         ! ray trace or eigenrays
                       CONTINUE
                    CASE DEFAULT
                       WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown RunType in 2D merge: ', Beam%RunType( 1 : 1 )
                       CALL ERROUT( 'BELLHOP3D', 'Unknown RunType in 2D merge' )
                    END SELECT
                 END IF
              END IF   ! closes iSingle test
           END DO AzimuthalAngle

           ! *** Scale the complex pressure field ***

           IF ( Beam%RunType( 6 : 6 ) == '3' ) THEN  ! 3D calculation
              SELECT CASE ( Beam%RunType( 1 : 1 ) )
              CASE ( 'C', 'S', 'I' )   ! TL calculation
                 CALL ScalePressure3D( Angles%Dalpha, Angles%Dbeta, ray2D( 1 )%c, epsilon, P, &
                                       Pos%Ntheta, NRz_per_range, Pos%NRr, Beam%RunType, freq )
              CASE ( 'A', 'a', 'R', 'E' )
                 CONTINUE
              CASE DEFAULT
                 WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown RunType in 3D scale: ', Beam%RunType( 1 : 1 )
                 CALL ERROUT( 'BELLHOP3D', 'Unknown RunType in 3D scale' )
              END SELECT
           END IF

           ! Write out the field

           SELECT CASE ( Beam%RunType( 1 : 1 ) )
           CASE ( 'C', 'S', 'I' )   ! TL calculation
              DO irz = 1, Pos%NRz
                 RcvrBearing: DO itheta = 1, Pos%Ntheta
                    iRec = 10 + ( isx    - 1 ) * Pos%Nsy * Pos%Ntheta * Pos%NSz * Pos%NRz + &
                                ( isy    - 1 ) *           Pos%Ntheta * Pos%NSz * Pos%NRz + &
                                ( itheta - 1 ) *                        Pos%NSz * Pos%NRz + &
                                ( isz    - 1 ) *                                  Pos%NRz + irz
                    WRITE( SHDFile, REC = IRec ) P( itheta, irz, 1 : Pos%NRr )

                 END DO RcvrBearing
              END DO
              CASE ( 'A' )             ! arrivals calculation, ascii
                 CALL WriteArrivalsASCII3D(  Pos%Rr, Pos%Ntheta, NRz_per_range, Pos%NRr )
              CASE ( 'a' )             ! arrivals calculation, binary
                 CALL WriteArrivalsBinary3D( Pos%Rr, Pos%Ntheta, NRz_per_range, Pos%NRr )
              CASE ( 'R', 'E' )         ! ray trace or eigenrays
                 CONTINUE
              CASE DEFAULT
                 WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown RunType when writing results: ', Beam%RunType( 1 : 1 )
                 CALL ERROUT( 'BELLHOP3D', 'Unknown RunType when writing results' )
           END SELECT
        END DO Source_y
     END DO Source_x
  END DO Source_z

  ! close all files
  SELECT CASE ( Beam%RunType( 1 : 1 ) )
  CASE ( 'C', 'S', 'I' )      ! TL calculation
     CLOSE( SHDFile )
  CASE ( 'A', 'a' )           ! arrivals calculation
     CLOSE( ARRFile )
  CASE ( 'R', 'E' )            ! ray trace or eigenrays
     CLOSE( RAYFile )
  CASE DEFAULT
     WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown RunType when closing files: ', Beam%RunType( 1 : 1 )
     CALL ERROUT( 'BELLHOP3D', 'Unknown RunType when closing files' )
  END SELECT

  ! Display run time
  CALL CPU_TIME( Tstop )
  WRITE( PRTFile, "( /, ' CPU Time = ', G15.3, 's' )" ) Tstop - Tstart

END SUBROUTINE BellhopCore

! **********************************************************************!

SUBROUTINE PickEpsilon( BeamType, omega, c, Dalpha, Dbeta, rLoop, EpsMultiplier, epsilon )

  ! Picks the optimum value for epsilon

  REAL      (KIND=8), INTENT( IN  ) :: omega, c             ! angular frequency, sound speed
  REAL      (KIND=8), INTENT( IN  ) :: Dalpha, Dbeta        ! angular spacing for ray fan
  REAL      (KIND=8), INTENT( IN  ) :: epsMultiplier, Rloop ! multiplier, loop range
  COMPLEX   (KIND=8), INTENT( OUT ) :: epsilon( 2 )         ! beam initial conditions
  CHARACTER (LEN= 2), INTENT( IN  ) :: BeamType
  LOGICAL, SAVE      :: INIFlag = .TRUE.
  REAL      (KIND=8), SAVE :: HalfWidth( 2 ) = [ 0.0, 0.0 ]
  COMPLEX   (KIND=8) :: epsilonOpt( 2 )
  CHARACTER (LEN=80) :: TAG

  ! LP: BUG: Multiple codepaths do not set epsilonOpt, leads to UB

  SELECT CASE ( BeamType( 1 : 1 ) )
  CASE ( 'C', 'R' )   ! Cerveny beams
     TAG    = 'Cerveny style beam'

     SELECT CASE ( BeamType( 2 : 2 ) )
     CASE ( 'F' )
        TAG            = 'Space filling beams'
        HalfWidth( 1 ) = 2.0 / ( ( omega / c ) * Dalpha )
        HalfWidth( 2 ) = 0.0
        IF ( Dbeta /= 0.0 ) HalfWidth( 2 ) = 2.0 / ( ( omega / c ) * Dbeta )
        epsilonOpt     = i * 0.5 * omega * HalfWidth ** 2
     CASE ( 'M' )
        TAG             = 'Minimum width beams'
        HalfWidth(  1 ) = SQRT( 2.0 * c * 1000.0 * rLoop / omega )
        HalfWidth(  2 ) = HalfWidth( 1 )
        epsilonOPT      = i * 0.5 * omega * HalfWidth **2
     CASE ( 'C' )
        TAG    = 'Cerveny style beam'
     CASE DEFAULT
        WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown beam width option: ', BeamType( 2 : 2 )
        CALL ERROUT( 'BELLHOP3D', 'Unknown beam width option' )
     END SELECT

  CASE ( 'g' )
     TAG             = 'Geometric beam, hat-shaped, Ray coord.'
     epsilonOPT      = 0.0
  CASE ( 'G', '^' )
     TAG             = 'Geometric beam, hat-shaped, Cart. coord.'
     epsilonOPT      = 0.0
  CASE ( 'b' )
     TAG             = 'Geometric beam, Gaussian-shaped, Ray coord.'
     epsilonOPT      = 0.0
  CASE ( 'B' )
     TAG             = 'Geometric beam, Gaussian-shaped, Cart. coord.'
     epsilonOPT      = 0.0
  CASE ( 'S' )
     TAG        = 'Simple Gaussian beams'
     halfwidth  = 2.0 / ( ( omega / c ) * Dalpha )
     epsilonOpt = i * 0.5 * omega * halfwidth ** 2
  CASE DEFAULT
     WRITE( ERROR_UNIT, * ) 'BELLHOP3D: Unknown beam type: ', BeamType( 1 : 1 )
     CALL ERROUT( 'BELLHOP3D', 'Unknown beam type' )
  END SELECT

  epsilon = epsMultiplier * epsilonOPT

  ! On first call write info to prt file
  IF ( INIFlag ) THEN
     WRITE( PRTFile, * )
     WRITE( PRTFile, * ) TAG
     WRITE( PRTFile, * ) 'HalfWidth1  = ', HalfWidth( 1 )
     WRITE( PRTFile, * ) 'HalfWidth2  = ', HalfWidth( 2 )
     WRITE( PRTFile, * ) 'epsilonOPT1 = ', epsilonOPT( 1 )
     WRITE( PRTFile, * ) 'epsilonOPT2 = ', epsilonOPT( 2 )
     WRITE( PRTFile, * ) 'EpsMult     = ', EpsMultiplier
     WRITE( PRTFile, * )
     INIFlag = .FALSE.
  END IF

END SUBROUTINE PickEpsilon


! **********************************************************************!

SUBROUTINE TraceRay2D( alpha, beta, Amp0 )

  ! Traces the beam corresponding to a particular take-off angle

  USE ReflectMod

  REAL     (KIND=8), INTENT( IN ) :: alpha, beta, Amp0 ! initial angles, amplitude
  INTEGER           :: is, is1                   ! index for a step along the ray
  REAL     (KIND=8) :: x( 3 )                    ! ray coordinate
  REAL     (KIND=8) :: t_o( 3 )                  ! tangent in ocean space
  !REAL     (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho
  REAL     (KIND=8) :: DistBegTop, DistEndTop, DistBegBot, DistEndBot ! Distances from ray beginning, end to top and bottom
  REAL     (KIND=8) :: tinit( 2 ), tradial( 2 ), BotnInt( 3 ), TopnInt( 3 ), s1, s2
  REAL     (KIND=8) :: z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy
  REAL     (KIND=8) :: term_minx, term_miny, term_maxx, term_maxy
  LOGICAL           :: topRefl, botRefl

  ! *** Initial conditions ***

  iSmallStepCtr = 0
  tinit = [ COS( alpha ), SIN( alpha ) ]
  tradial = [ COS( beta ), SIN( beta ) ]
  ray2D( 1 )%x = [ 0.0D0, xs_3D( 3 ) ]

  ! CALL EvaluateSSP2D( ray2D( 1 )%x, tinit, c, cimag, gradc, crr, crz, czz, rho, xs_3D, tradial, freq )
  ray2D( 1 )%t         = tinit / ray2D( 1 )%c ! LP: Set before PickEpsilon independent of ray angles; never overwritten
  ray2D( 1 )%p         = [ 1.0, 0.0 ]
  ray2D( 1 )%q         = [ 0.0, 1.0 ]
  ray2D( 1 )%tau       = 0.0
  ray2D( 1 )%Amp       = Amp0
  ray2D( 1 )%Phase     = 0.0
  ray2D( 1 )%NumTopBnc = 0
  ray2D( 1 )%NumBotBnc = 0

  IsegTopx = 1
  IsegTopy = 1
  IsegBotx = 1
  IsegBoty = 1
  t_o = RayToOceanT( ray2D( 1 )%t, tradial )
  CALL GetTopSeg3D( xs_3D, t_o, .TRUE. )   ! identify the top    segment above the source
  CALL GetBotSeg3D( xs_3D, t_o, .TRUE. )   ! identify the bottom segment below the source

  ! Trace the beam (note that Reflect alters the step index is)
  is = 0
  CALL Distances3D( xs_3D, Topx, Botx, Topn, Botn, DistBegTop, DistBegBot )

  IF ( DistBegTop <= 0 .OR. DistBegBot <= 0 ) THEN
     Beam%Nsteps = 1
     RETURN       ! source must be within the medium
  END IF

  Stepping: DO WHILE ( .TRUE. )
     is  = is + 1
     is1 = is + 1

     CALL Step2D( ray2D( is ), ray2D( is1 ), tradial, topRefl, botRefl )

     ! convert polar coordinate of ray to x-y coordinate
     x = RayToOceanX( ray2D( is1 )%x, xs_3D, tradial )
     t_o = RayToOceanT( ray2D( is1 )%t, tradial )

     CALL GetTopSeg3D( x, t_o, .FALSE. )    ! identify the top    segment above the source
     CALL GetBotSeg3D( x, t_o, .FALSE. )    ! identify the bottom segment below the source

     IF ( IsegTopx == 0 .OR. IsegTopy == 0 .OR. IsegBotx == 0 .OR. IsegBoty == 0 ) THEN ! we escaped the box
        Beam%Nsteps = is
        EXIT Stepping
     END IF

     ! Reflections?
     ! Tests that ray at step is is inside, and ray at step is+1 is outside
     ! to detect only a crossing from inside to outside
     ! DistBeg is the distance at step is,   which is saved
     ! DistEnd is the distance at step is+1, which needs to be calculated

     CALL Distances3D( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )

     IF ( topRefl ) THEN
        IF ( atiType == 'C' ) THEN

           ! LP: This is superfluous, it's the same x calculated above.
           x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial )

           s1     = ( x( 1 ) - Topx( 1 ) ) / ( xTopSeg( 2 ) - xTopSeg( 1 ) )   ! proportional distance along segment
           s2     = ( x( 2 ) - Topx( 2 ) ) / ( yTopSeg( 2 ) - yTopSeg( 1 ) )   ! proportional distance along segment

           TopnInt = Top( IsegTopx,     IsegTopy     )%Noden * ( 1 - s1 ) * ( 1 - s2 ) +  &
                     Top( IsegTopx + 1, IsegTopy     )%Noden * ( s1     ) * ( 1 - s2 ) +  &
                     Top( IsegTopx + 1, IsegTopy + 1 )%Noden * ( s1     ) * ( s2     ) +  &
                     Top( IsegTopx,     IsegTopy + 1 )%Noden * ( 1 - s1 ) * ( s2     )

           z_xx = Top( IsegTopx, IsegTopy )%z_xx
           z_xy = Top( IsegTopx, IsegTopy )%z_xy
           z_yy = Top( IsegTopx, IsegTopy )%z_yy

           kappa_xx = Top( IsegTopx, IsegTopy )%kappa_xx
           kappa_xy = Top( IsegTopx, IsegTopy )%kappa_xy
           kappa_yy = Top( IsegTopx, IsegTopy )%kappa_yy
        ELSE
           TopnInt = Topn   ! normal is constant in a segment
           z_xx = 0
           z_xy = 0
           z_yy = 0
           kappa_xx = 0
           kappa_xy = 0
           kappa_yy = 0
        END IF

        CALL Reflect2D( is, Bdry%Top%HS, 'TOP', TopnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RTop, NTopPTS, tradial)
        ray2D( is + 1 )%NumTopBnc = ray2D( is )%NumTopBnc + 1

        x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial )

        CALL Distances3D( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )

     ELSE IF ( botRefl ) THEN
        ! write( *, * ) 'Reflecting', x, Botx
        ! write( *, * ) 'Botn', Botn
        ! write( *, * ) 'Distances', DistEndTop, DistEndBot
        IF ( btyType == 'C' ) THEN

           x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial )

           s1     = ( x( 1 ) - Botx( 1 ) ) / ( xBotSeg( 2 ) - xBotSeg( 1 ) )   ! proportional distance along segment
           s2     = ( x( 2 ) - Botx( 2 ) ) / ( yBotSeg( 2 ) - yBotSeg( 1 ) )   ! proportional distance along segment

           BotnInt = Bot( IsegBotx,     IsegBoty     )%Noden * ( 1 - s1 ) * ( 1 - s2 ) +  &
                     Bot( IsegBotx + 1, IsegBoty     )%Noden * ( s1     ) * ( 1 - s2 ) +  &
                     Bot( IsegBotx + 1, IsegBoty + 1 )%Noden * ( s1     ) * ( s2     ) +  &
                     Bot( IsegBotx,     IsegBoty + 1 )%Noden * ( 1 - s1 ) * ( s2     )

           z_xx = Bot( IsegBotx, IsegBoty )%z_xx
           z_xy = Bot( IsegBotx, IsegBoty )%z_xy
           z_yy = Bot( IsegBotx, IsegBoty )%z_yy

           kappa_xx = Bot( IsegBotx, IsegBoty )%kappa_xx
           kappa_xy = Bot( IsegBotx, IsegBoty )%kappa_xy
           kappa_yy = Bot( IsegBotx, IsegBoty )%kappa_yy
        ELSE
           BotnInt = Botn   ! normal is constant in a segment
           z_xx = 0
           z_xy = 0
           z_yy = 0
           kappa_xx = 0
           kappa_xy = 0
           kappa_yy = 0
        END IF

        CALL Reflect2D( is, Bdry%Bot%HS, 'BOT', BotnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RBot, NBotPTS, tradial)
        ray2D( is + 1 )%NumBotBnc = ray2D( is )%NumBotBnc + 1

        x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial )

        CALL Distances3D( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )

     END IF

     ! Has the ray left the box, lost its energy, escaped the boundaries, or exceeded storage limit?
     ! LP: See explanation for changes in bdry3DMod: GetTopSeg3D.
     t_o = RayToOceanT( ray2D( is + 1 )%t, tradial )
     term_minx = MAX( Bot( 1,            1 )%x( 1 ), Top( 1,            1 )%x( 1 ) )
     term_miny = MAX( Bot( 1,            1 )%x( 2 ), Top( 1,            1 )%x( 2 ) )
     term_maxx = MIN( Bot( NBTYPts( 1 ), 1 )%x( 1 ), Top( NATIPts( 1 ), 1 )%x( 1 ) )
     term_maxy = MIN( Bot( 1, NBTYPts( 2 ) )%x( 2 ), Top( 1, NATIPts( 2 ) )%x( 2 ) )
     ! LP: The Beam%Box conditions were inexplicably commented out in 2022 revision of Nx2D, see README.
     IF ( ABS( x( 1 ) - xs_3D( 1 ) ) >= Beam%Box%x .OR. &
          ABS( x( 2 ) - xs_3D( 2 ) ) >= Beam%Box%y .OR. &
          ABS( x( 3 )              ) >= Beam%Box%z .OR. &  ! LP: Removed xs( 3 ) for consistency
          x( 1 ) < term_minx .OR. &
          x( 2 ) < term_miny .OR. &
          x( 1 ) > term_maxx .OR. &
          x( 2 ) > term_maxy .OR. &
          ( x( 1 ) == term_minx .AND. t_o( 1 ) < 0.0 ) .OR. &
          ( x( 2 ) == term_miny .AND. t_o( 2 ) < 0.0 ) .OR. &
          ( x( 1 ) == term_maxx .AND. t_o( 1 ) > 0.0 ) .OR. &
          ( x( 2 ) == term_maxy .AND. t_o( 2 ) > 0.0 ) .OR. &
          ray2D( is + 1 )%Amp < 0.005 .OR. &
          ! ray2D( is + 1 )%t( 1 ) < 0  .OR. & ! kills off a backward traveling ray
          iSmallStepCtr > 50 ) THEN
        Beam%Nsteps = is + 1
        EXIT Stepping
     ELSE IF ( is >= MaxN - 3 ) THEN
        WRITE( PRTFile, * ) 'Warning in TraceRay2D : Insufficient storage for ray trajectory'
        WRITE( PRTFile, * ) 'Angles are  alpha = ', alpha * RadDeg, '    beta = ', beta * RadDeg
        Beam%Nsteps = is
        EXIT Stepping
     END IF

     DistBegTop = DistEndTop
     DistBegBot = DistEndBot

  END DO Stepping

END SUBROUTINE TraceRay2D

! **********************************************************************!

SUBROUTINE Step2D( ray0, ray2, tradial, topRefl, botRefl )

  ! Does a single step along the ray
  ! x denotes the ray coordinate, (r,z)
  ! t denotes the scaled tangent to the ray (previously (rho, zeta))
  ! c * t would be the unit tangent

  USE Step3DMod

  REAL (KIND=8), INTENT( IN ) :: tradial( 2 )   ! coordinate of source and ray bearing angle
  LOGICAL, INTENT( OUT ) :: topRefl, botRefl
  TYPE( ray2DPt )    :: ray1
  TYPE( ray2DPt ), INTENT( IN )    :: ray0
  TYPE( ray2DPt ), INTENT( INOUT ) :: ray2
  INTEGER            :: iSegx0, iSegy0, iSegz0, snapDim
  REAL     (KIND=8 ) :: gradc0( 2 ), gradc1( 2 ), gradc2( 2 ), rho, &
                        c0, cimag0, crr0, crz0, czz0, csq0, cnn0_csq0, &
                        c1, cimag1, crr1, crz1, czz1, csq1, cnn1_csq1, &
                        c2, cimag2, crr2, crz2, czz2, urayt0( 2 ), urayt1( 2 ), urayt2( 2 ), &
                        h, halfh, hw0, hw1, ray2n( 2 ), RM, RN, gradcjump( 2 ), cnjump, csjump, w0, w1, &
                        rayx3D( 3 ), rayt3D( 3 ), x2_o( 3 ), x2_o_out( 3 )

  IF ( STEP_DEBUGGING ) THEN
     WRITE( PRTFile, * )
     WRITE( PRTFile, * ) 'ray0 x t', ray0%x, ray0%t
     WRITE( PRTFile, * ) 'iSegr iSegz', iSegr, iSegz
  END IF

  ! The numerical integrator used here is a version of the polygon (a.k.a. midpoint, leapfrog, or Box method), and similar
  ! to the Heun (second order Runge-Kutta method).
  ! However, it's modified to allow for a dynamic step change, while preserving the second-order accuracy).

  ! *** Phase 1 (an Euler step)

  CALL EvaluateSSP2D( ray0%x, ray0%t, c0, cimag0, gradc0, crr0, crz0, czz0, rho, xs_3D, tradial, freq )

  csq0      = c0 * c0
  cnn0_csq0 = crr0 * ray0%t( 2 )**2 - 2.0 * crz0 * ray0%t( 1 ) * ray0%t( 2 ) + czz0 * ray0%t( 1 )**2

  iSegx0 = iSegx     ! make note of current layer
  iSegy0 = iSegy
  iSegz0 = iSegz

  h = Beam%deltas            ! initially set the step h, to the basic one, deltas

  urayt0 = c0 * ray0%t  ! unit tangent
  rayx3D = RayToOceanX( ray0%x, xs_3D, tradial )
  rayt3D = RayToOceanT( urayt0, tradial )

  CALL ReduceStep3D( rayx3D, rayt3D, iSegx0, iSegy0, iSegz0, h ) ! reduce h to land on boundary

  halfh = 0.5 * h   ! first step of the modified polygon method is a half step

  ray1%x = ray0%x + halfh * urayt0
  ray1%t = ray0%t - halfh * gradc0 / csq0
  ray1%p = ray0%p - halfh * cnn0_csq0 * ray0%q
  ray1%q = ray0%q + halfh * c0        * ray0%p

  ! *** Phase 2

  CALL EvaluateSSP2D( ray1%x, ray1%t, c1, cimag1, gradc1, crr1, crz1, czz1, rho, xs_3D, tradial, freq )
  csq1      = c1 * c1
  cnn1_csq1 = crr1 * ray1%t( 2 )**2 - 2.0 * crz1 * ray1%t( 1 ) * ray1%t( 2 ) + czz1 * ray1%t( 1 )**2

  ! The Munk test case with a horizontally launched ray caused problems.
  ! The ray vertexes on an interface and can ping-pong around that interface.
  ! Have to be careful in that case about big changes to the stepsize (that invalidate the leap-frog scheme) in phase II.
  ! A modified Heun or Box method could also work.

  urayt1 = c1 * ray1%t   ! unit tangent
  rayt3D = RayToOceanT( urayt1, tradial )

  CALL ReduceStep3D( rayx3D, rayt3D, iSegx0, iSegy0, iSegz0, h ) ! reduce h to land on boundary

  ! use blend of f' based on proportion of a full step used.
  w1  = h / ( 2.0d0 * halfh )
  w0  = 1.0d0 - w1
  urayt2 = w0 * urayt0 + w1 * urayt1
  rayt3D = RayToOceanT( urayt2, tradial )
  ! Take the blended ray tangent ( urayt2 ) and find the minimum step size ( h )
  ! to put this on a boundary, and ensure that the resulting position
  ! ( ray2%x ) gets put precisely on the boundary.
  CALL StepToBdry3D( rayx3D, x2_o, rayt3D, iSegx0, iSegy0, iSegz0, h, &
     topRefl, botRefl, snapDim )
  ray2%x = OceanToRayX( x2_o, xs_3D, tradial, urayt2, snapDim )
  IF ( STEP_DEBUGGING ) THEN
     x2_o_out = RayToOceanX( ray2%x, xs_3D, tradial )
     WRITE( PRTFile, * ) 'OceanToRayX in ', x2_o
     WRITE( PRTFile, * ) '           ==> ', ray2%x
     WRITE( PRTFile, * ) '           ==> ', x2_o_out
  END IF
  !write( *, * ) 'final coord ', x2_o, ray2%x, w0, w1

  hw0 = h * w0
  hw1 = h * w1
  ray2%t   = ray0%t   - hw0 * gradc0 / csq0      - hw1 * gradc1 / csq1
  ray2%p   = ray0%p   - hw0 * cnn0_csq0 * ray0%q - hw1 * cnn1_csq1 * ray1%q
  ray2%q   = ray0%q   + hw0 * c0        * ray0%p + hw1 * c1        * ray1%p
  ray2%tau = ray0%tau + hw0 / CMPLX( c0, cimag0, KIND=8 ) + hw1 / CMPLX( c1, cimag1, KIND=8 )

  ray2%Amp       = ray0%Amp
  ray2%Phase     = ray0%Phase
  ray2%NumTopBnc = ray0%NumTopBnc
  ray2%NumBotBnc = ray0%NumBotBnc

  ! If we crossed an interface, apply jump condition

  CALL EvaluateSSP2D( ray2%x, ray2%t, c2, cimag2, gradc2, crr2, crz2, czz2, rho, xs_3D, tradial, freq )
  ray2%c = c2

  !!! this needs modifying like the full 3D version to handle jumps in the x-y direction
  IF ( iSegz /= iSegz0 ) THEN
     gradcjump =  gradc2 - gradc0  ! this is precise only for c-linear layers
     ray2n     = [ -ray2%t( 2 ), ray2%t( 1 ) ]

     cnjump    = DOT_PRODUCT( gradcjump, ray2n  )
     csjump    = DOT_PRODUCT( gradcjump, ray2%t )

     RM        = ray2%t( 1 ) / ray2%t( 2 )
     RN        = RM * ( 2 * cnjump - RM * csjump ) / c2
     RN        = -RN
     ray2%p    = ray2%p + ray2%q * RN

  END IF

END SUBROUTINE Step2D

FUNCTION RayToOceanX( x, xs, tradial )
!! Transform ray coordinates to ocean coordinates
  REAL ( KIND=8 ) :: RayToOceanX( 3 )
  REAL ( KIND=8 ), INTENT( IN ) :: x( 2 ), xs( 3 ), tradial( 2 )
  RayToOceanX = [ xs( 1 ) + x( 1 ) * tradial( 1 ), &
                  xs( 2 ) + x( 1 ) * tradial( 2 ), &
                  x( 2 ) ]
END FUNCTION RayToOceanX

FUNCTION RayToOceanT( t, tradial )
!! Transform ray tangent to ocean tangent coordinates
  REAL ( KIND=8 ) :: RayToOceanT( 3 )
  REAL ( KIND=8 ), INTENT( IN  ) :: t( 2 ), tradial( 2 )
  RayToOceanT = [ t( 1 ) * tradial( 1 ), &
                  t( 1 ) * tradial( 2 ), &
                  t( 2 ) ]
END FUNCTION RayToOceanT

FUNCTION OceanToRayX( x, xs, tradial, t, snapDim )
!! Transform ocean coordinates to ray coordinates

  ! LP: Going back and forth through the coordinate transform won't
  ! always keep the precise value, so we may have to finesse the floats.
  ! Valid values of snapDim:
  ! -2: Snap to X or Y unspecified
  ! -1: No snap
  !  0: Snap to X
  !  1: Snap to Y
  !  2: Snap to Z

  REAL ( KIND=8 ) :: OceanToRayX( 2 )
  REAL ( KIND=8 ), INTENT( IN ) :: x( 3 ), xs( 3 ), tradial( 2 ), t( 2 )
  INTEGER, INTENT( IN ) :: snapDim
  REAL ( KIND=8 ) :: ret( 2 ), x_back( 3 ), wantdir( 2 ), errdir( 2 )
  LOGICAL :: correctdir( 2 )
  INTEGER :: i

  ! Depth always transfers perfectly--not changed.
  ret( 2 ) = x( 3 )
  ! For range, use larger dimension--this avoids divide-by-zero or divide
  ! by a small number causing accuracy problems.
  IF ( ABS( tradial( 1 ) ) >= ABS( tradial( 2 ) ) ) THEN
     ret( 1 ) = ( x( 1 ) - xs( 1 ) ) / tradial( 1 )
  ELSE
     ret( 1 ) = ( x( 2 ) - xs( 2 ) ) / tradial( 2 )
  END IF
  IF ( snapDim < -2 .OR. snapDim == -1 .OR. snapDim >= 2 ) THEN
     ! Either:
     ! snapDim out of range (won't happen, but want to help compiler)
     ! No snap selected--this is the best estimate
     ! Snap to Z--Z already perfect, this is the best estimate
     OceanToRayX = ret
     RETURN
  END IF
  ! Only do this iteration a few times, then give up.
  DO i = 1, 4
     ! Go back from 2D to 3D, compare to original x.
     x_back = RayToOceanX( ret, xs, tradial )
     ! If we can't be on the boundary, we want to be slightly forward of
     ! the boundary, measured in terms of the ray tangent range. This also
     ! encompasses cases where one component exactly matches (errdir.x_or_y
     ! == RL(0.0)). For both of these values, only the sign matters.
     wantdir = tradial * t( 1 )
     errdir = x_back( 1 : 2 ) - x( 1 : 2 )
     correctdir = ( wantdir * errdir ) >= 0.0D0
     IF ( ( snapDim == 0 .AND. correctdir( 1 ) ) .OR. &
          ( snapDim == 1 .AND. correctdir( 2 ) ) .OR. &
          ( correctdir( 1 ) .AND. correctdir( 2 ) ) ) THEN
        OceanToRayX = ret
        RETURN
     END IF
     ! Move to the next floating point value for ret, in the direction
     ! of the ray tangent. How do we know this will actually produce a
     ! new value for x_back? If the scale of the floating point steps
     ! around x_back is larger than that of ret (several values of ret
     ! map to the same value of x_back after the transform), there
     ! should be no problem in finding a value that maps exactly, so
     ! we would not get here.
     ret( 1 ) = NEAREST( ret( 1 ), t( 1 ) )
  END DO
  WRITE( PRTFile, * ) 'Warning in OceanToRayX: Failed to transform 3D -> 2D -> 3D consistently'
  OceanToRayX = ret
END FUNCTION OceanToRayX

! **********************************************************************!

SUBROUTINE TraceRay3D( alpha, beta, epsilon, Amp0 )
!! Traces the beam corresponding to a particular take off angle

  USE Step3DMod
  USE Reflect3DMod

  REAL     ( KIND=8 ), INTENT( IN ) :: Amp0  ! initial amplitude
  REAL     ( KIND=8 ), INTENT( IN ) :: alpha, beta    ! take-off angles of the ray
  COMPLEX  ( KIND=8 ), INTENT( IN ) :: epsilon( 2 )   ! beam initial conditions
  INTEGER             :: is, is1
  REAL     ( KIND=8 ) :: DistBegTop, DistEndTop, DistBegBot, DistEndBot!, &
                         !c, cimag, gradc( 3 ), cxx, cyy, czz, cxy, cxz, cyz, rho   ! soundspeed derivatives
  REAL     (KIND=8) :: TopnInt( 3 ), BotnInt( 3 )
  REAL     (KIND=8) :: s1, s2
  REAL     (KIND=8) :: z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy
  REAL     (KIND=8) :: tinit( 3 )
  LOGICAL           :: topRefl, botRefl
  REAL     (KIND=8) :: term_minx, term_miny, term_maxx, term_maxy, term_x( 3 ), term_t( 3 )

  ! *** Initial conditions ***

  iSmallStepCtr = 0
  ray3D( 1 )%x    = xs_3D
  tinit = [ COS( alpha ) * COS( beta ), COS( alpha ) * SIN( beta ), SIN( alpha ) ]
  !CALL EvaluateSSP3D(  ray3D( 1 )%x, tinit, c, cimag, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, 'TAB' )

  ray3D( 1 )%t    = tinit / ray3D( 1 )%c ! LP: Set before PickEpsilon independent of ray angles; never overwritten
  !ray3D( 1 )%f    = epsilon( 2 )
  !ray3D( 1 )%g    = epsilon( 1 )
  !ray3D( 1 )%h    = 0.0
  !ray3D( 1 )%DetP = 1.0
  !ray3D( 1 )%DetQ = epsilon( 1 ) * epsilon( 2 )
  !ray3D( 1 )%c    = c   ! LP
  ray3D( 1 )%phi  = 0.0

  ray3D( 1 )%tau       = 0.0
  ray3D( 1 )%Amp       = Amp0
  ray3D( 1 )%Phase     = 0.0
  ray3D( 1 )%NumTopBnc = 0
  ray3D( 1 )%NumBotBnc = 0

  !ray3D( 1 )%p_tilde = [ 1.0, 0.0 ]   use these for complex Cerveny beams
  !ray3D( 1 )%q_tilde = [ epsilon( 1 ), CMPLX( 0.0, KIND=8 ) ]
  !ray3D( 1 )%p_hat   = [ 0.0, 1.0 ]
  !ray3D( 1 )%q_hat   = [ CMPLX( 0.0, KIND=8 ), epsilon( 2 ) ]

  ray3D( 1 )%p_tilde = [ 1.0, 0.0 ]
  ray3D( 1 )%q_tilde = [ 0.0, 0.0 ]
  ray3D( 1 )%p_hat   = [ 0.0, 1.0 ]
  ray3D( 1 )%q_hat   = [ 0.0, 0.0 ]

  ! dummy BotSeg info to force GetBotSeg to search for the active segment on first call
  xTopSeg = [ +big, -big ]
  yTopSeg = [ +big, -big ]
  xBotSeg = [ +big, -big ]
  yBotSeg = [ +big, -big ]

  CALL GetTopSeg3D( xs_3D, ray3D( 1 )%t, .TRUE. )   ! identify the top    segment above the source
  CALL GetBotSeg3D( xs_3D, ray3D( 1 )%t, .TRUE. )   ! identify the bottom segment below the source

  ! Trace the beam (note that Reflect alters the step index is)
  is = 0

  CALL Distances3D( ray3D( 1 )%x, Topx,  Botx, Topn, Botn, DistBegTop, DistBegBot )

  IF ( DistBegTop <= 0 .OR. DistBegBot <= 0 ) THEN
     Beam%Nsteps = 1
     WRITE( PRTFile, * ) 'Terminating the ray trace because the source is on or outside the boundaries'
     RETURN       ! source must be within the medium
  END IF

  Stepping: DO WHILE ( .TRUE. )
     is  = is + 1
     is1 = is + 1

     CALL Step3D( ray3D( is ), ray3D( is1 ), topRefl, botRefl )

     CALL GetTopSeg3D( ray3D( is1 )%x, ray3D( is1 )%t, .FALSE. )   ! identify the top    segment above the source
     CALL GetBotSeg3D( ray3D( is1 )%x, ray3D( is1 )%t, .FALSE. )   ! identify the bottom segment below the source

     IF ( IsegTopx == 0 .OR. IsegTopy == 0 .OR. IsegBotx == 0 .OR. IsegBoty == 0 ) THEN ! we escaped the box
        Beam%Nsteps = is
        EXIT Stepping
     END IF

     ! Reflections?
     ! Tests that ray at step is is inside, and ray at step is+1 is outside
     ! to detect only a crossing from inside to outside
     ! DistBeg is the distance at step is, which is saved
     ! DistEnd is the distance at step is+1, which needs to be calculated

     CALL Distances3D( ray3D( is1 )%x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )

     IF ( topRefl ) THEN
        IF ( STEP_DEBUGGING ) THEN
           WRITE( PRTFile, * ) 'Top reflecting'
        END IF
        IF ( atiType == 'C' ) THEN
           s1 = ( ray3D( is1 )%x( 1 ) - Topx( 1 ) ) / ( xTopSeg( 2 ) - xTopSeg( 1 ) )   ! proportional distance along segment
           s2 = ( ray3D( is1 )%x( 2 ) - Topx( 2 ) ) / ( yTopSeg( 2 ) - yTopSeg( 1 ) )   ! proportional distance along segment

           TopnInt = Top( IsegTopx,     IsegTopy     )%Noden * ( 1 - s1 ) * ( 1 - s2 ) +  &
                     Top( IsegTopx + 1, IsegTopy     )%Noden * ( s1     ) * ( 1 - s2 ) +  &
                     Top( IsegTopx + 1, IsegTopy + 1 )%Noden * ( s1     ) * ( s2     ) +  &
                     Top( IsegTopx,     IsegTopy + 1 )%Noden * ( 1 - s1 ) * ( s2     )
           z_xx = Top( IsegTopx, IsegTopy )%z_xx
           z_xy = Top( IsegTopx, IsegTopy )%z_xy
           z_yy = Top( IsegTopx, IsegTopy )%z_yy

           kappa_xx = Top( IsegBotx, IsegBoty )%kappa_xx
           kappa_xy = Top( IsegBotx, IsegBoty )%kappa_xy
           kappa_yy = Top( IsegBotx, IsegBoty )%kappa_yy
        ELSE
           TopnInt = Topn   ! normal is constant in a segment
           z_xx = 0
           z_xy = 0
           z_yy = 0

           kappa_xx = 0
           kappa_xy = 0
           kappa_yy = 0
        END IF

        CALL Reflect3D( is, Bdry%Top%HS, 'TOP', TopnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RTop, NTopPTS )
        ray3D( is + 1 )%NumTopBnc = ray3D( is )%NumTopBnc + 1
        CALL Distances3D( ray3D( is + 1 )%x, Topx,  Botx, Topn, Botn, DistEndTop, DistEndBot )

     ELSE IF ( botRefl ) THEN
        IF ( STEP_DEBUGGING ) THEN
           WRITE( PRTFile, * ) 'Bottom reflecting'
        END IF

        IF ( btyType == 'C' ) THEN
           s1 = ( ray3D( is1 )%x( 1 ) - Botx( 1 ) ) / ( xBotSeg( 2 ) - xBotSeg( 1 ) )   ! proportional distance along segment
           s2 = ( ray3D( is1 )%x( 2 ) - Botx( 2 ) ) / ( yBotSeg( 2 ) - yBotSeg( 1 ) )   ! proportional distance along segment

           BotnInt = Bot( IsegBotx,     IsegBoty     )%Noden * ( 1 - s1 ) * ( 1 - s2 ) +  &
                     Bot( IsegBotx + 1, IsegBoty     )%Noden * ( s1     ) * ( 1 - s2 ) +  &
                     Bot( IsegBotx + 1, IsegBoty + 1 )%Noden * ( s1     ) * ( s2     ) +  &
                     Bot( IsegBotx,     IsegBoty + 1 )%Noden * ( 1 - s1 ) * ( s2     )
           z_xx = Bot( IsegBotx, IsegBoty )%z_xx
           z_xy = Bot( IsegBotx, IsegBoty )%z_xy
           z_yy = Bot( IsegBotx, IsegBoty )%z_yy

           kappa_xx = Bot( IsegBotx, IsegBoty )%kappa_xx
           kappa_xy = Bot( IsegBotx, IsegBoty )%kappa_xy
           kappa_yy = Bot( IsegBotx, IsegBoty )%kappa_yy
        ELSE
           BotnInt = Botn   ! normal is constant in a segment
           z_xx = 0
           z_xy = 0
           z_yy = 0

           kappa_xx = 0
           kappa_xy = 0
           kappa_yy = 0
        END IF

        CALL Reflect3D( is, Bdry%Bot%HS, 'BOT', BotnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RBot, NBotPTS )
        ray3D( is + 1 )%NumBotBnc = ray3D( is )%NumBotBnc + 1
        CALL Distances3D( ray3D( is + 1 )%x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )

     END IF

     ! Has the ray exited the beam box, lost its energy, escaped the BTY, ATI boundaries, or exceeded storage limit?
     ! LP: See explanation for changes in bdry3DMod: GetTopSeg3D.
     term_x = ray3D( is + 1 )%x
     term_t = ray3D( is + 1 )%t
     term_minx = MAX( Bot( 1,            1 )%x( 1 ), Top( 1,            1 )%x( 1 ) )
     term_miny = MAX( Bot( 1,            1 )%x( 2 ), Top( 1,            1 )%x( 2 ) )
     term_maxx = MIN( Bot( NBTYPts( 1 ), 1 )%x( 1 ), Top( NATIPts( 1 ), 1 )%x( 1 ) )
     term_maxy = MIN( Bot( 1, NBTYPts( 2 ) )%x( 2 ), Top( 1, NATIPts( 2 ) )%x( 2 ) )
     IF ( &
          ABS( term_x( 1 ) - xs_3D( 1 ) ) >= Beam%Box%x .OR. &
          ABS( term_x( 2 ) - xs_3D( 2 ) ) >= Beam%Box%y .OR. &
          ABS( term_x( 3 )              ) >= Beam%Box%z .OR. & ! box is centered at z=0
          term_x( 1 ) < term_minx .OR. &
          term_x( 2 ) < term_miny .OR. &
          term_x( 1 ) > term_maxx .OR. &
          term_x( 2 ) > term_maxy .OR. &
          ( term_x( 1 ) == term_minx .AND. term_t( 1 ) < 0.0 ) .OR. &
          ( term_x( 2 ) == term_miny .AND. term_t( 2 ) < 0.0 ) .OR. &
          ( term_x( 1 ) == term_maxx .AND. term_t( 1 ) > 0.0 ) .OR. &
          ( term_x( 2 ) == term_maxy .AND. term_t( 2 ) > 0.0 ) .OR. &
          ray3D( is + 1 )%Amp < 0.005 .OR. &
          iSmallStepCtr > 50 ) THEN
        Beam%Nsteps = is + 1
        EXIT Stepping
     ELSE IF ( is >= MaxN - 3 ) THEN
        Beam%Nsteps = is
        WRITE( PRTFile, * ) 'Warning in TraceRay3D : Insufficient storage for ray trajectory'
        WRITE( PRTFile, * ) 'Angles are  alpha = ', alpha * RadDeg, '    beta = ', beta * RadDeg
        EXIT Stepping
     END IF

     DistBegTop = DistEndTop
     DistBegBot = DistEndBot

  END DO Stepping
END SUBROUTINE TraceRay3D

! **********************************************************************!

SUBROUTINE Distances3D( rayx, Topx, Botx, Topn, Botn, DistTop, DistBot )
!! Computes distances from ray to boundaries

  ! Formula differs from JKPS because code uses outward pointing normals
  ! Note that Topx, Botx, just need to be any node on the diagonal that divides each square into triangles
  ! In bdry3DMod, the node is selected as the one at the lowest x, y, index and that defines the triangles

  REAL (KIND=8), INTENT( IN  ) :: rayx( 3 )             ! ray coordinate
  REAL (KIND=8), INTENT( IN  ) :: Topx( 3 ), Botx( 3 )  ! top, bottom boundary coordinate for node
  REAL (KIND=8), INTENT( IN  ) :: Topn( 3 ), Botn( 3 )  ! top, bottom boundary normal
  REAL (KIND=8), INTENT( OUT ) :: DistTop, DistBot      ! distance from the ray to top, bottom boundaries
  REAL (KIND=8)                :: dTop( 3 ), dBot( 3 )

  dTop    = rayx - Topx  ! vector pointing from top    to ray
  dBot    = rayx - Botx  ! vector pointing from bottom to ray
  DistTop = -DOT_PRODUCT( Topn, dTop )
  DistBot = -DOT_PRODUCT( Botn, dBot )

!!$  write( *, * )
!!$  write( *, * ) 'Distances3D', DistBot
!!$  write( *, * ) 'rayx', rayx
!!$  write( *, * ) 'Botx', Botx
!!$  write( *, * ) 'dBot', dBot
!!$  write( *, * ) 'Normal', Botn
!!$  write( *, * )
END SUBROUTINE Distances3D

END PROGRAM BELLHOP3D
