BellhopCore Subroutine

SUBROUTINE BellhopCore
  !! Core subroutine to run Bellhop algorithm

  USE ArrMod
  USE AttenMod

  INTEGER, PARAMETER   :: ArrivalsStorage = 200000000, MinNArr = 10
  INTEGER              :: IBPvec( 1 ), ibp, is, iBeamWindow2, Irz1, Irec, NalphaOpt, iSeg
  REAL                 :: Tstart, Tstop
  REAL        (KIND=8) :: Amp0, DalphaOpt, xs( 2 ), tdummy( 2 ), RadMax, s, &
                          c, cimag, gradc( 2 ), crr, crz, czz, rho
  COMPLEX, ALLOCATABLE :: U( :, : )
  COMPLEX     (KIND=8) :: epsilon

  CALL CPU_TIME( Tstart )

  omega = 2.0 * pi * freq

  IF ( Beam%deltas == 0.0 ) THEN
     Beam%deltas = ( Bdry%Bot%HS%Depth - Bdry%Top%HS%Depth ) / 10.0   ! Automatic step size selection
     WRITE( PRTFile, * )
     WRITE( PRTFile, fmt = '(  '' Step length,       deltas = '', G11.4, '' m (automatically selected)'' )' ) Beam%deltas
  END IF

  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

  ! convert range-dependent geoacoustic parameters from user to program units
  IF ( atiType( 2 : 2 ) == 'L' ) THEN
     DO iSeg = 1, NatiPts
        Top( iSeg )%HS%cp = CRCI( 1D20, Top( iSeg )%HS%alphaR, Top( iSeg )%HS%alphaI, freq, freq, 'W ', &
             betaPowerLaw, ft )   ! compressional wave speed
        Top( iSeg )%HS%cs = CRCI( 1D20, Top( iSeg )%HS%betaR,  Top( iSeg )%HS%betaI,  freq, freq, 'W ', &
             betaPowerLaw, ft )   ! shear         wave speed
     END DO
  END IF

  IF ( btyType( 2 : 2 ) == 'L' ) THEN
     DO iSeg = 1, NbtyPts
        Bot( iSeg )%HS%cp = CRCI( 1D20, Bot( iSeg )%HS%alphaR, Bot( iSeg )%HS%alphaI, freq, freq, 'W ', &
             betaPowerLaw, ft )   ! compressional wave speed
        Bot( iSeg )%HS%cs = CRCI( 1D20, Bot( iSeg )%HS%betaR,  Bot( iSeg )%HS%betaI,  freq, freq, 'W ', &
             betaPowerLaw, ft )   ! shear         wave speed
     END DO
  END IF

  SELECT CASE ( Beam%RunType( 5 : 5 ) )
  CASE ( 'I' )
     NRz_per_range = 1         ! irregular grid
  CASE DEFAULT
     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 ( U( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
     IF ( IAllocStat /= 0 ) &
          CALL ERROUT( 'BELLHOP', 'Insufficient memory for TL matrix: reduce Nr * NRz'  )
  CASE ( 'A', 'a', 'R', 'E' )   ! Arrivals calculation
     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 / ( NRz_per_range * Pos%NRr ), MinNArr )   ! allow space for at least 10 arrivals
     WRITE( PRTFile, * )
     WRITE( PRTFile, * ) '( Maximum # of arrivals = ', MaxNArr, ')'

     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' )
  CASE DEFAULT
     MaxNArr = 1
     ALLOCATE ( Arr( NRz_per_range, Pos%NRr, 1 ), NArr( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
  END SELECT

  WRITE( PRTFile, * )

  SourceDepth: DO is = 1, Pos%NSz
     xs = [ 0.0, Pos%sz( is ) ]   ! source coordinate

     SELECT CASE ( Beam%RunType( 1 : 1 ) )
     CASE ( 'C', 'S', 'I' ) ! TL calculation, zero out pressure matrix
        U = 0.0
     CASE ( 'A', 'a' )      ! Arrivals calculation, zero out arrival matrix
        NArr = 0
     END SELECT

     ! LP: This initialization was missing (only done once globally). In some
     ! cases (a source on a boundary), in conjunction with other subtle issues,
     ! the missing initialization caused the initial state of one ray to be
     ! dependent on the final state of the previous ray, which is obviously
     ! non-physical.
     ! This initialization is here in addition to TraceRay2D because of EvaluateSSP.
     iSegz = 1
     iSegr = 1

     tdummy = [ 1.0, 0.0 ]
     CALL EvaluateSSP( xs, tdummy, c, cimag, gradc, crr, crz, czz, rho, freq, 'TAB' )
     RadMax = 5 * c / freq  ! 5 wavelength max radius

     ! Are there enough beams?
     DalphaOpt = SQRT( c / ( 6.0 * freq * Pos%Rr( Pos%NRr ) ) )
     NalphaOpt = 2 + INT( ( Angles%alpha( Angles%Nalpha ) - Angles%alpha( 1 ) ) / DalphaOpt )

     IF ( Beam%RunType( 1 : 1 ) == 'C' .AND. Angles%Nalpha < NalphaOpt ) THEN
        WRITE( PRTFile, * ) 'Warning in BELLHOP : Too few beams'
        WRITE( PRTFile, * ) 'Nalpha should be at least = ', NalphaOpt
     ENDIF

     ! Trace successive beams

     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?

           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 / c * xs( 2 ) * SIN( Angles%alpha( ialpha ) ) ) )

           ! show progress ...
           IF ( MOD( ialpha - 1, max( Angles%Nalpha / 50, 1 ) ) == 0 ) THEN
              WRITE( PRTFile, FMT = "( 'Tracing beam ', I7, F10.2 )" ) ialpha, SrcDeclAngle
              FLUSH( PRTFile )
           END IF

           CALL TraceRay2D( xs, Angles%alpha( ialpha ), Amp0 )   ! Trace a ray

           IF ( Beam%RunType( 1 : 1 ) == 'R' ) THEN   ! Write the ray trajectory to RAYFile
              CALL WriteRay2D( SrcDeclAngle, Beam%Nsteps )
           ELSE                                       ! Compute the contribution to the field

              epsilon = PickEpsilon( Beam%Type( 1 : 2 ), omega, c, gradc, Angles%alpha( ialpha ), &
                   Angles%Dalpha, Beam%rLoop, Beam%epsMultiplier ) ! 'optimal' beam constant

              SELECT CASE ( Beam%Type( 1 : 1 ) )
              CASE ( 'R' )
                 iBeamWindow2 = Beam%iBeamWindow **2
                 RadMax       = 50 * c / freq  ! 50 wavelength max radius
                 CALL InfluenceCervenyRayCen(   U, epsilon, Angles%alpha( ialpha ), iBeamWindow2, RadMax )
              CASE ( 'C' )
                 iBeamWindow2 = Beam%iBeamWindow **2
                 RadMax       = 50 * c / freq  ! 50 wavelength max radius
                 CALL InfluenceCervenyCart(     U, epsilon, Angles%alpha( ialpha ), iBeamWindow2, RadMax )
              CASE ( 'g' )
                 CALL InfluenceGeoHatRayCen(    U, Angles%alpha( ialpha ), Angles%Dalpha )
              CASE ( 'S' )
                 RadMax       = 50 * c / freq  ! 50 wavelength max radius
                 CALL InfluenceSGB(             U, Angles%alpha( ialpha ), Angles%Dalpha, RadMax )
              CASE ( 'B' )
                 CALL InfluenceGeoGaussianCart( U, Angles%alpha( ialpha ), Angles%Dalpha )
             CASE DEFAULT
                 CALL InfluenceGeoHatCart(      U, Angles%alpha( ialpha ), Angles%Dalpha )
              END SELECT

           END IF
        END IF
     END DO DeclinationAngle

     ! write results to disk

     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 )
        IRec = 10 + NRz_per_range * ( is - 1 )
        RcvrDepth: DO Irz1 = 1, NRz_per_range
           IRec = IRec + 1
           WRITE( SHDFile, REC = IRec ) U( Irz1, 1 : Pos%NRr )
        END DO RcvrDepth
     CASE ( 'A' )             ! arrivals calculation, ascii
        CALL WriteArrivalsASCII(  Pos%Rr, NRz_per_range, Pos%NRr, Beam%RunType( 4 : 4 ) )
     CASE ( 'a' )             ! arrivals calculation, binary
        CALL WriteArrivalsBinary( Pos%Rr, NRz_per_range, Pos%NRr, Beam%RunType( 4 : 4 ) )
     END SELECT

  END DO SourceDepth

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

  ! 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
     CLOSE( RAYFile )
  END SELECT

  CLOSE( PRTFile )

END SUBROUTINE BellhopCore