BellhopCore Subroutine

subroutine BellhopCore()

Uses

  • proc~~bellhopcore~~UsesGraph proc~bellhopcore BellhopCore module~anglemod anglemod proc~bellhopcore->module~anglemod module~arrmod ArrMod proc~bellhopcore->module~arrmod module~sourcereceiverpositions SourceReceiverPositions proc~bellhopcore->module~sourcereceiverpositions module~writeray WriteRay proc~bellhopcore->module~writeray module~anglemod->module~sourcereceiverpositions module~fatalerror FatalError module~anglemod->module~fatalerror module~mathconstants MathConstants module~anglemod->module~mathconstants module~sortmod SortMod module~anglemod->module~sortmod module~subtabulate SubTabulate module~anglemod->module~subtabulate module~bellhopmod bellhopMod module~arrmod->module~bellhopmod module~arrmod->module~mathconstants module~sourcereceiverpositions->module~fatalerror module~monotonicmod monotonicMod module~sourcereceiverpositions->module~monotonicmod module~sourcereceiverpositions->module~sortmod module~sourcereceiverpositions->module~subtabulate module~writeray->module~bellhopmod module~sspmod sspmod module~writeray->module~sspmod module~bellhopmod->module~mathconstants module~sspmod->module~fatalerror module~sspmod->module~monotonicmod module~splinec splinec module~sspmod->module~splinec

Core subroutine to run Bellhop algorithm

Arguments

None

Calls

proc~~bellhopcore~~CallsGraph proc~bellhopcore BellhopCore proc~errout ERROUT proc~bellhopcore->proc~errout proc~evaluatessp3d EvaluateSSP3D proc~bellhopcore->proc~evaluatessp3d proc~influence3dgeogaussiancart Influence3DGeoGaussianCart proc~bellhopcore->proc~influence3dgeogaussiancart proc~influence3dgeogaussianraycen Influence3DGeoGaussianRayCen proc~bellhopcore->proc~influence3dgeogaussianraycen proc~influence3dgeohatcart Influence3DGeoHatCart proc~bellhopcore->proc~influence3dgeohatcart proc~influence3dgeohatraycen Influence3DGeoHatRayCen proc~bellhopcore->proc~influence3dgeohatraycen proc~influencecervenycart InfluenceCervenyCart proc~bellhopcore->proc~influencecervenycart proc~influencecervenyraycen InfluenceCervenyRayCen proc~bellhopcore->proc~influencecervenyraycen proc~influencegeogaussiancart InfluenceGeoGaussianCart proc~bellhopcore->proc~influencegeogaussiancart proc~influencegeohatcart InfluenceGeoHatCart proc~bellhopcore->proc~influencegeohatcart proc~influencegeohatraycen InfluenceGeoHatRayCen proc~bellhopcore->proc~influencegeohatraycen proc~influencesgb InfluenceSGB proc~bellhopcore->proc~influencesgb proc~pickepsilon PickEpsilon proc~bellhopcore->proc~pickepsilon proc~scalepressure ScalePressure proc~bellhopcore->proc~scalepressure proc~scalepressure3d ScalePressure3D proc~bellhopcore->proc~scalepressure3d proc~traceray2d TraceRay2D proc~bellhopcore->proc~traceray2d proc~traceray3d TraceRay3D proc~bellhopcore->proc~traceray3d proc~writearrivalsascii3d WriteArrivalsASCII3D proc~bellhopcore->proc~writearrivalsascii3d proc~writearrivalsbinary3d WriteArrivalsBinary3D proc~bellhopcore->proc~writearrivalsbinary3d proc~writeray3d WriteRay3D proc~bellhopcore->proc~writeray3d sngl sngl proc~bellhopcore->sngl proc~evaluatessp3d->proc~errout proc~analytic3d Analytic3D proc~evaluatessp3d->proc~analytic3d proc~ccubic cCubic proc~evaluatessp3d->proc~ccubic proc~clinear cLinear proc~evaluatessp3d->proc~clinear proc~hexahedral Hexahedral proc~evaluatessp3d->proc~hexahedral proc~n2linear n2Linear proc~evaluatessp3d->proc~n2linear proc~applycontribution ApplyContribution proc~influence3dgeogaussiancart->proc~applycontribution proc~raynormal_unit RayNormal_unit proc~influence3dgeogaussiancart->proc~raynormal_unit proc~scalebeam ScaleBeam proc~influence3dgeogaussiancart->proc~scalebeam interface~cross_product cross_product proc~influence3dgeogaussianraycen->interface~cross_product proc~influence3dgeogaussianraycen->proc~applycontribution proc~raynormal RayNormal proc~influence3dgeogaussianraycen->proc~raynormal proc~rtoir RToIR proc~influence3dgeogaussianraycen->proc~rtoir proc~influence3dgeogaussianraycen->proc~scalebeam proc~influence3dgeohatcart->proc~applycontribution proc~influence3dgeohatcart->proc~raynormal_unit proc~influence3dgeohatcart->proc~scalebeam proc~influence3dgeohatraycen->interface~cross_product proc~influence3dgeohatraycen->proc~applycontribution proc~influence3dgeohatraycen->proc~raynormal proc~influence3dgeohatraycen->proc~rtoir proc~influence3dgeohatraycen->proc~scalebeam proc~influencecervenycart->proc~errout proc~branchcut BranchCut proc~influencecervenycart->proc~branchcut proc~evaluatessp EvaluateSSP proc~influencecervenycart->proc~evaluatessp proc~hermite Hermite proc~influencecervenycart->proc~hermite proc~influencecervenycart->proc~rtoir proc~influencecervenyraycen->proc~errout proc~influencecervenyraycen->proc~branchcut proc~influencecervenyraycen->proc~hermite proc~influencecervenyraycen->proc~rtoir proc~applycontribution~2 ApplyContribution proc~influencegeogaussiancart->proc~applycontribution~2 proc~finalphase FinalPhase proc~influencegeogaussiancart->proc~finalphase proc~incphaseifcaustic IncPhaseIfCaustic proc~influencegeogaussiancart->proc~incphaseifcaustic proc~influencegeohatcart->proc~applycontribution~2 proc~influencegeohatcart->proc~finalphase proc~influencegeohatcart->proc~incphaseifcaustic proc~influencegeohatraycen->proc~applycontribution~2 proc~influencegeohatraycen->proc~finalphase proc~influencegeohatraycen->proc~incphaseifcaustic proc~influencegeohatraycen->proc~rtoir proc~influencesgb->proc~applycontribution~2 proc~influencesgb->proc~incphaseifcaustic proc~scalepressure->sngl proc~distances3d Distances3D proc~traceray2d->proc~distances3d proc~getbotseg3d GetBotSeg3D proc~traceray2d->proc~getbotseg3d proc~gettopseg3d GetTopSeg3D proc~traceray2d->proc~gettopseg3d proc~raytooceant RayToOceanT proc~traceray2d->proc~raytooceant proc~raytooceanx RayToOceanX proc~traceray2d->proc~raytooceanx proc~reflect2d Reflect2D proc~traceray2d->proc~reflect2d proc~step2d~2 Step2D proc~traceray2d->proc~step2d~2 proc~traceray3d->proc~distances3d proc~traceray3d->proc~getbotseg3d proc~traceray3d->proc~gettopseg3d proc~reflect3d Reflect3D proc~traceray3d->proc~reflect3d proc~step3d Step3D proc~traceray3d->proc~step3d proc~writearrivalsbinary3d->sngl proc~cross_product_dble cross_product_dble interface~cross_product->proc~cross_product_dble proc~cross_product_sngl cross_product_sngl interface~cross_product->proc~cross_product_sngl proc~applycontribution->proc~writeray3d proc~applycontribution->sngl proc~addarr3d AddArr3D proc~applycontribution->proc~addarr3d proc~applycontribution~2->proc~writeray3d proc~applycontribution~2->sngl proc~addarr AddArr proc~applycontribution~2->proc~addarr proc~writeray2d WriteRay2D proc~applycontribution~2->proc~writeray2d proc~readssp ReadSSP proc~ccubic->proc~readssp proc~splineall SPLINEALL proc~ccubic->proc~splineall proc~updatedepthsegmentt UpdateDepthSegmentT proc~ccubic->proc~updatedepthsegmentt proc~clinear->proc~readssp proc~clinear->proc~updatedepthsegmentt proc~evaluatessp->proc~errout proc~evaluatessp->proc~ccubic proc~evaluatessp->proc~clinear proc~evaluatessp->proc~hexahedral proc~evaluatessp->proc~n2linear proc~analytic Analytic proc~evaluatessp->proc~analytic proc~cpchip cPCHIP proc~evaluatessp->proc~cpchip proc~quad Quad proc~evaluatessp->proc~quad proc~isatcaustic IsAtCaustic proc~finalphase->proc~isatcaustic proc~getbotseg3d->proc~errout proc~gettopseg3d->proc~errout proc~hexahedral->proc~errout interface~monotonic monotonic proc~hexahedral->interface~monotonic proc~hexahedral->proc~readssp proc~update3dxsegmentt Update3DXSegmentT proc~hexahedral->proc~update3dxsegmentt proc~update3dysegmentt Update3DYSegmentT proc~hexahedral->proc~update3dysegmentt proc~update3dzsegmentt Update3DZSegmentT proc~hexahedral->proc~update3dzsegmentt proc~incphaseifcaustic->proc~isatcaustic proc~n2linear->proc~readssp proc~n2linear->proc~updatedepthsegmentt proc~evaluatessp2d EvaluateSSP2D proc~reflect2d->proc~evaluatessp2d proc~interpolatereflectioncoefficient InterpolateReflectionCoefficient proc~reflect2d->proc~interpolatereflectioncoefficient proc~reflect3d->proc~evaluatessp3d proc~reflect3d->interface~cross_product proc~reflect3d->proc~raynormal proc~reflect3d->proc~interpolatereflectioncoefficient proc~step2d~2->proc~raytooceant proc~step2d~2->proc~raytooceanx proc~step2d~2->proc~evaluatessp2d proc~oceantorayx OceanToRayX proc~step2d~2->proc~oceantorayx proc~reducestep3d ReduceStep3D proc~step2d~2->proc~reducestep3d proc~steptobdry3d StepToBdry3D proc~step2d~2->proc~steptobdry3d proc~step3d->proc~evaluatessp3d proc~step3d->interface~cross_product proc~step3d->proc~raynormal proc~step3d->proc~raynormal_unit proc~computedeltapq ComputeDeltaPQ proc~step3d->proc~computedeltapq proc~step3d->proc~reducestep3d proc~step3d->proc~steptobdry3d proc~updateraypq UpdateRayPQ proc~step3d->proc~updateraypq proc~monotonic_dble monotonic_dble interface~monotonic->proc~monotonic_dble proc~monotonic_sngl monotonic_sngl interface~monotonic->proc~monotonic_sngl proc~addarr->sngl proc~addarr3d->sngl proc~cpchip->proc~readssp proc~cpchip->proc~updatedepthsegmentt proc~pchip PCHIP proc~cpchip->proc~pchip proc~evaluatessp2d->proc~evaluatessp3d proc~oceantorayx->proc~raytooceanx proc~quad->proc~errout proc~quad->interface~monotonic proc~quad->proc~readssp proc~quad->proc~updatedepthsegmentt proc~updaterangesegmentt UpdateRangeSegmentT proc~quad->proc~updaterangesegmentt proc~readssp->proc~errout proc~crci CRCI proc~readssp->proc~crci proc~crci->proc~errout proc~franc_garr Franc_Garr proc~crci->proc~franc_garr proc~cspline CSPLINE proc~pchip->proc~cspline proc~fprime_interior_cmplx fprime_interior_Cmplx proc~pchip->proc~fprime_interior_cmplx proc~fprime_left_end_cmplx fprime_left_end_Cmplx proc~pchip->proc~fprime_left_end_cmplx proc~fprime_right_end_cmplx fprime_right_end_Cmplx proc~pchip->proc~fprime_right_end_cmplx proc~h_del h_del proc~pchip->proc~h_del proc~fprime_interior fprime_interior proc~fprime_interior_cmplx->proc~fprime_interior proc~fprime_left_end fprime_left_end proc~fprime_left_end_cmplx->proc~fprime_left_end proc~fprime_right_end fprime_right_end proc~fprime_right_end_cmplx->proc~fprime_right_end

Called by

proc~~bellhopcore~~CalledByGraph proc~bellhopcore BellhopCore program~bellhop3d BELLHOP3D program~bellhop3d->proc~bellhopcore

Source Code

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 ) &
     Angles%Dalpha = ( Angles%alpha( Angles%Nalpha ) - Angles%alpha( 1 ) ) / ( Angles%Nalpha - 1 )  ! angular spacing between beams
  SELECT CASE ( Beam%RunType( 5 : 5 ) )
  CASE ( 'I' )
     NRz_per_range = 1         ! irregular grid
  CASE ( 'R' )
     NRz_per_range = Pos%NRz   ! rectilinear grid
  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 ) &
          CALL ERROUT( 'BELLHOP', 'Insufficient memory for TL matrix: reduce Nr * NRz'  )
  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
  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' ) &
                          Amp0 = Amp0 * SQRT( 2.0 ) * ABS( SIN( omega / c0 * xs_3D( 3 ) * SIN( Angles%alpha( ialpha ) ) ) )

                       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
                       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 )
                       ENDIF
                    ENDIF   ! 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
                    END SELECT
                 END IF
              ENDIF   ! 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 )
              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 )
           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' )                ! ray trace
     CLOSE( RAYFile )
  END SELECT

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

END SUBROUTINE BellhopCore