1
-
!! Provides BELLHOP program definition
4
-
!! BELLHOP beam tracing for ocean acoustics
6
-
! BELLHOP Beam tracing for ocean acoustics
8
-
! Copyright (C) 2009 Michael B. Porter
10
-
! This program is free software: you can redistribute it and/or modify
11
-
! it under the terms of the GNU General Public License as published by
12
-
! the Free Software Foundation, either version 3 of the License, or
13
-
! (at your option) any later version.
15
-
! This program is distributed in the hope that it will be useful,
16
-
! but WITHOUT ANY WARRANTY; without even the implied warranty of
17
-
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18
-
! GNU General Public License for more details.
20
-
! You should have received a copy of the GNU General Public License
21
-
! along with this program. If not, see <http://www.gnu.org/licenses/>.
23
-
! First version (1983) originally developed with Homer Bucker, Naval Ocean Systems Center
25
29
USE ReadEnvironmentBell
29
-
USE SourceReceiverPositions
37
-
LOGICAL, PARAMETER :: Inline = .FALSE.
39
-
CHARACTER ( LEN=2 ) :: AttenUnit
40
-
CHARACTER ( LEN=80 ) :: FileRoot
44
-
! get the file root for naming all input and output files
45
-
! should add some checks here ...
47
29
CALL GET_COMMAND_ARGUMENT( 1, FileRoot )
48
-
! Open the print file
49
29*
OPEN( UNIT = PRTFile, FILE = TRIM( FileRoot ) // '.prt', STATUS = 'UNKNOWN', IOSTAT = iostat )
51
-
! Read in or otherwise initialize inline all the variables used by BELLHOP
54
-
! NPts, Sigma not used by BELLHOP
55
-
Title = 'BELLHOP- Calibration case with envfil passed as parameters'
57
-
! NMedia variable is not used by BELLHOP
59
-
! *** Boundary information (type of boundary condition and, if a halfspace, then halfspace info)
62
-
Bdry%Top%HS%BC = 'V'
63
-
Bdry%Top%HS%Depth = 0
64
-
Bdry%Bot%HS%Depth = 100
65
-
Bdry%Bot%HS%Opt = 'A_'
66
-
Bdry%Bot%HS%BC = 'A'
67
-
Bdry%Bot%HS%cp = CRCI( 1D20, 1590D0, 0.5D0, freq, freq, AttenUnit, betaPowerLaw, fT ) ! compressional wave speed
68
-
Bdry%Bot%HS%cs = CRCI( 1D20, 0D0, 0D0, freq, freq, AttenUnit, betaPowerLaw, fT ) ! shear wave speed
69
-
Bdry%Bot%HS%rho = 1.2 ! density
71
-
! *** sound speed in the water column ***
73
-
SSP%Type = 'C' ! interpolation method for SSP
74
-
SSP%NPts = 2 ! number of SSP points
75
-
SSP%z( 1 : 2 ) = [ 0, 100 ]
76
-
SSP%c( 1 : 2 ) = [ 1500, 1500 ]
77
-
SSP%cz( 1 : 2 ) = [ 0, 0 ] ! user should really not have to supply this ...
79
-
! *** source and receiver positions ***
85
-
ALLOCATE( Pos%sz( Pos%NSz ), Pos%ws( Pos%NSz ), Pos%isz( Pos%NSz ) )
86
-
ALLOCATE( Pos%rz( Pos%NRz ), Pos%wr( Pos%NRz ), Pos%irz( Pos%NRz ) )
87
-
ALLOCATE( Pos%Rr( Pos%NRr ) )
90
-
!Pos%Rz = [ 0, 50, 100 ]
91
-
!Pos%Rr = 1000. * [ 1, 2, 3, 4, 5 ] ! meters !!!
92
-
Pos%Rz = [ ( jj, jj = 1, Pos%NRz ) ]
93
-
Pos%Rr = 10. * [ ( jj, jj = 1 , Pos%NRr ) ] ! meters !!!
99
-
Beam%Box%r = 5100 ! meters
101
-
Angles%Nalpha = 1789
102
-
! Angles%alpha = [ -80, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80 ] ! -89 89
103
-
Angles%alpha = ( 180. / Angles%Nalpha ) * [ ( jj, jj = 1, Angles%Nalpha ) ] - 90.
105
-
! *** altimetry ***
107
-
ALLOCATE( Top( 2 ), Stat = IAllocStat )
108
-
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', 'Insufficient memory for altimetry data' )
109
-
Top( 1 )%x = [ -sqrt( huge( Top( 1 )%x( 1 ) ) ) / 1.0d5, 0.d0 ]
110
-
Top( 2 )%x = [ sqrt( huge( Top( 1 )%x( 1 ) ) ) / 1.0d5, 0.d0 ]
112
-
CALL ComputeBdryTangentNormal( Top, 'Top' )
114
-
! *** bathymetry ***
116
-
ALLOCATE( Bot( 2 ), Stat = IAllocStat )
117
-
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', 'Insufficient memory for bathymetry data' )
118
-
Bot( 1 )%x = [ -sqrt( huge( Bot( 1 )%x( 1 ) ) ) / 1.0d5, 5000.d0 ]
119
-
Bot( 2 )%x = [ sqrt( huge( Bot( 1 )%x( 1 ) ) ) / 1.0d5, 5000.d0 ]
121
-
CALL ComputeBdryTangentNormal( Bot, 'Bot' )
123
-
ALLOCATE( RBot( 1 ), Stat = IAllocStat ) ! bottom reflection coefficient
124
-
ALLOCATE( RTop( 1 ), Stat = iAllocStat ) ! top reflection coefficient
125
-
NBotPts = 1 ! LP: missing initialization
128
-
! *** Source Beam Pattern ***
130
-
ALLOCATE( SrcBmPat( 2, 2 ), Stat = IAllocStat )
131
-
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP-ReadPat', 'Insufficient memory' )
132
-
SrcBmPat( 1, : ) = [ -180.0, 0.0 ]
133
-
SrcBmPat( 2, : ) = [ 180.0, 0.0 ]
134
-
SrcBmPat( :, 2 ) = 10 ** ( SrcBmPat( :, 2 ) / 20 ) ! convert dB to linear scale !!!
136
29
CALL ReadEnvironment( FileRoot, ThreeD )
137
14
CALL ReadATI( FileRoot, Bdry%Top%HS%Opt( 5 : 5 ), Bdry%Top%HS%Depth, PRTFile ) ! AlTImetry
138
14
CALL ReadBTY( FileRoot, Bdry%Bot%HS%Opt( 2 : 2 ), Bdry%Bot%HS%Depth, PRTFile ) ! BaThYmetry
139
14
CALL ReadReflectionCoefficient( FileRoot, Bdry%Bot%HS%Opt( 1 : 1 ), Bdry%Top%HS%Opt( 2 : 2 ), PRTFile ) ! (top and bottom)
140
14
SBPFlag = Beam%RunType( 3 : 3 )
141
14
CALL ReadPat( FileRoot, PRTFile ) ! Source Beam Pattern
142
-
! dummy bearing angles
144
14*
ALLOCATE( Pos%theta( Pos%Ntheta ), Stat = IAllocStat )
145
14*
Pos%theta( 1 ) = 0.
148
14
CALL OpenOutputFiles( FileRoot, ThreeD )
153
-
! **********************************************************************!
155
14
SUBROUTINE BellhopCore
156
-
!! Core subroutine to run Bellhop algorithm
161
-
INTEGER, PARAMETER :: ArrivalsStorage = 200000000, MinNArr = 10
162
-
INTEGER :: IBPvec( 1 ), ibp, is, iBeamWindow2, Irz1, Irec, NalphaOpt, iSeg
163
-
REAL :: Tstart, Tstop
164
-
REAL (KIND=8) :: Amp0, DalphaOpt, xs( 2 ), tdummy( 2 ), RadMax, s, &
165
-
c, cimag, gradc( 2 ), crr, crz, czz, rho
166
14
COMPLEX, ALLOCATABLE :: U( :, : )
167
-
COMPLEX (KIND=8) :: epsilon
169
14
CALL CPU_TIME( Tstart )
171
14
omega = 2.0 * pi * freq
173
14
IF ( Beam%deltas == 0.0 ) THEN
174
14
Beam%deltas = ( Bdry%Bot%HS%Depth - Bdry%Top%HS%Depth ) / 10.0 ! Automatic step size selection
175
14
WRITE( PRTFile, * )
176
14
WRITE( PRTFile, fmt = '( '' Step length, deltas = '', G11.4, '' m (automatically selected)'' )' ) Beam%deltas
179
55964*
Angles%alpha = DegRad * Angles%alpha ! convert to radians
180
14
Angles%Dalpha = 0.0
181
14
IF ( Angles%Nalpha /= 1 ) &
182
14*
Angles%Dalpha = ( Angles%alpha( Angles%Nalpha ) - Angles%alpha( 1 ) ) / ( Angles%Nalpha - 1 ) ! angular spacing between beams
184
-
! convert range-dependent geoacoustic parameters from user to program units
185
14
IF ( atiType( 2 : 2 ) == 'L' ) THEN
186
#####
DO iSeg = 1, NatiPts
187
#####
Top( iSeg )%HS%cp = CRCI( 1D20, Top( iSeg )%HS%alphaR, Top( iSeg )%HS%alphaI, freq, freq, 'W ', &
188
#####
betaPowerLaw, ft ) ! compressional wave speed
189
#####
Top( iSeg )%HS%cs = CRCI( 1D20, Top( iSeg )%HS%betaR, Top( iSeg )%HS%betaI, freq, freq, 'W ', &
190
#####
betaPowerLaw, ft ) ! shear wave speed
194
14
IF ( btyType( 2 : 2 ) == 'L' ) THEN
195
#####
DO iSeg = 1, NbtyPts
196
#####
Bot( iSeg )%HS%cp = CRCI( 1D20, Bot( iSeg )%HS%alphaR, Bot( iSeg )%HS%alphaI, freq, freq, 'W ', &
197
#####
betaPowerLaw, ft ) ! compressional wave speed
198
#####
Bot( iSeg )%HS%cs = CRCI( 1D20, Bot( iSeg )%HS%betaR, Bot( iSeg )%HS%betaI, freq, freq, 'W ', &
199
#####
betaPowerLaw, ft ) ! shear wave speed
203
#####
SELECT CASE ( Beam%RunType( 5 : 5 ) )
205
#####
NRz_per_range = 1 ! irregular grid
207
14
NRz_per_range = Pos%NRz ! rectilinear grid
210
-
! for a TL calculation, allocate space for the pressure matrix
211
4
SELECT CASE ( Beam%RunType( 1 : 1 ) )
212
-
CASE ( 'C', 'S', 'I' ) ! TL calculation
213
4*
ALLOCATE ( U( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
214
4
IF ( IAllocStat /= 0 ) &
215
4*
CALL ERROUT( 'BELLHOP', 'Insufficient memory for TL matrix: reduce Nr * NRz' )
216
-
CASE ( 'A', 'a', 'R', 'E' ) ! Arrivals calculation
217
14*
ALLOCATE ( U( 1, 1 ), Stat = IAllocStat ) ! open a dummy variable
220
-
! for an arrivals run, allocate space for arrivals matrices
221
8
SELECT CASE ( Beam%RunType( 1 : 1 ) )
223
8
MaxNArr = MAX( ArrivalsStorage / ( NRz_per_range * Pos%NRr ), MinNArr ) ! allow space for at least 10 arrivals
224
8
WRITE( PRTFile, * )
225
8
WRITE( PRTFile, * ) '( Maximum # of arrivals = ', MaxNArr, ')'
227
8*
ALLOCATE ( Arr( NRz_per_range, Pos%NRr, MaxNArr ), NArr( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
228
8
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', &
229
8*
'Insufficient memory to allocate arrivals matrix; reduce parameter ArrivalsStorage' )
232
14*
ALLOCATE ( Arr( NRz_per_range, Pos%NRr, 1 ), NArr( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
235
14
WRITE( PRTFile, * )
237
28*
SourceDepth: DO is = 1, Pos%NSz
238
42*
xs = [ 0.0, Pos%sz( is ) ] ! source coordinate
240
4
SELECT CASE ( Beam%RunType( 1 : 1 ) )
241
-
CASE ( 'C', 'S', 'I' ) ! TL calculation, zero out pressure matrix
243
-
CASE ( 'A', 'a' ) ! Arrivals calculation, zero out arrival matrix
247
-
! LP: This initialization was missing (only done once globally). In some
248
-
! cases (a source on a boundary), in conjunction with other subtle issues,
249
-
! the missing initialization caused the initial state of one ray to be
250
-
! dependent on the final state of the previous ray, which is obviously
252
-
! This initialization is here in addition to TraceRay2D because of EvaluateSSP.
256
14
tdummy = [ 1.0, 0.0 ]
257
14
CALL EvaluateSSP( xs, tdummy, c, cimag, gradc, crr, crz, czz, rho, freq, 'TAB' )
258
14
RadMax = 5 * c / freq ! 5 wavelength max radius
260
-
! Are there enough beams?
261
14*
DalphaOpt = SQRT( c / ( 6.0 * freq * Pos%Rr( Pos%NRr ) ) )
262
14*
NalphaOpt = 2 + INT( ( Angles%alpha( Angles%Nalpha ) - Angles%alpha( 1 ) ) / DalphaOpt )
264
14
IF ( Beam%RunType( 1 : 1 ) == 'C' .AND. Angles%Nalpha < NalphaOpt ) THEN
265
1
WRITE( PRTFile, * ) 'Warning in BELLHOP : Too few beams'
266
1
WRITE( PRTFile, * ) 'Nalpha should be at least = ', NalphaOpt
269
-
! Trace successive beams
271
55964*
DeclinationAngle: DO ialpha = 1, Angles%Nalpha
272
55950*
SrcDeclAngle = RadDeg * Angles%alpha( ialpha ) ! take-off declination angle in degrees
274
55964
IF ( Angles%iSingle_alpha == 0 .OR. ialpha == Angles%iSingle_alpha ) THEN ! Single beam run?
276
167850*
IBPvec = maxloc( SrcBmPat( :, 1 ), mask = SrcBmPat( :, 1 ) < SrcDeclAngle ) ! index of ray angle in beam pattern
277
55950
IBP = IBPvec( 1 )
278
55950
IBP = MAX( IBP, 1 ) ! don't go before beginning of table
279
55950
IBP = MIN( IBP, NSBPPts - 1 ) ! don't go past end of table
281
-
! linear interpolation to get amplitude
282
55950*
s = ( SrcDeclAngle - SrcBmPat( IBP, 1 ) ) / ( SrcBmPat( IBP + 1, 1 ) - SrcBmPat( IBP, 1 ) )
283
55950*
Amp0 = ( 1 - s ) * SrcBmPat( IBP, 2 ) + s * SrcBmPat( IBP + 1, 2 )
285
-
! Lloyd mirror pattern for semi-coherent option
286
55950
IF ( Beam%RunType( 1 : 1 ) == 'S' ) &
287
#####
Amp0 = Amp0 * SQRT( 2.0 ) * ABS( SIN( omega / c * xs( 2 ) * SIN( Angles%alpha( ialpha ) ) ) )
289
-
! show progress ...
290
55950
IF ( MOD( ialpha - 1, max( Angles%Nalpha / 50, 1 ) ) == 0 ) THEN
291
700
WRITE( PRTFile, FMT = "( 'Tracing beam ', I7, F10.2 )" ) ialpha, SrcDeclAngle
295
55950*
CALL TraceRay2D( xs, Angles%alpha( ialpha ), Amp0 ) ! Trace a ray
297
55950
IF ( Beam%RunType( 1 : 1 ) == 'R' ) THEN ! Write the ray trajectory to RAYFile
298
50
CALL WriteRay2D( SrcDeclAngle, Beam%Nsteps )
299
-
ELSE ! Compute the contribution to the field
301
#####
epsilon = PickEpsilon( Beam%Type( 1 : 2 ), omega, c, gradc, Angles%alpha( ialpha ), &
302
55900
Angles%Dalpha, Beam%rLoop, Beam%epsMultiplier ) ! 'optimal' beam constant
304
#####
SELECT CASE ( Beam%Type( 1 : 1 ) )
306
#####
iBeamWindow2 = Beam%iBeamWindow **2
307
#####
RadMax = 50 * c / freq ! 50 wavelength max radius
308
#####
CALL InfluenceCervenyRayCen( U, epsilon, Angles%alpha( ialpha ), iBeamWindow2, RadMax )
310
#####
iBeamWindow2 = Beam%iBeamWindow **2
311
#####
RadMax = 50 * c / freq ! 50 wavelength max radius
312
#####
CALL InfluenceCervenyCart( U, epsilon, Angles%alpha( ialpha ), iBeamWindow2, RadMax )
314
#####
CALL InfluenceGeoHatRayCen( U, Angles%alpha( ialpha ), Angles%Dalpha )
316
#####
RadMax = 50 * c / freq ! 50 wavelength max radius
317
#####
CALL InfluenceSGB( U, Angles%alpha( ialpha ), Angles%Dalpha, RadMax )
319
4600*
CALL InfluenceGeoGaussianCart( U, Angles%alpha( ialpha ), Angles%Dalpha )
321
60500*
CALL InfluenceGeoHatCart( U, Angles%alpha( ialpha ), Angles%Dalpha )
326
-
END DO DeclinationAngle
328
-
! write results to disk
330
28
SELECT CASE ( Beam%RunType( 1 : 1 ) )
331
-
CASE ( 'C', 'S', 'I' ) ! TL calculation
332
4*
CALL ScalePressure( Angles%Dalpha, ray2D( 1 )%c, Pos%Rr, U, NRz_per_range, Pos%NRr, Beam%RunType, freq )
333
4
IRec = 10 + NRz_per_range * ( is - 1 )
334
1312*
RcvrDepth: DO Irz1 = 1, NRz_per_range
336
1308*
WRITE( SHDFile, REC = IRec ) U( Irz1, 1 : Pos%NRr )
338
-
CASE ( 'A' ) ! arrivals calculation, ascii
339
8*
CALL WriteArrivalsASCII( Pos%Rr, NRz_per_range, Pos%NRr, Beam%RunType( 4 : 4 ) )
340
-
CASE ( 'a' ) ! arrivals calculation, binary
341
26*
CALL WriteArrivalsBinary( Pos%Rr, NRz_per_range, Pos%NRr, Beam%RunType( 4 : 4 ) )
347
14
CALL CPU_TIME( Tstop )
348
14
WRITE( PRTFile, "( /, ' CPU Time = ', G15.3, 's' )" ) Tstop - Tstart
351
4
SELECT CASE ( Beam%RunType( 1 : 1 ) )
352
-
CASE ( 'C', 'S', 'I' ) ! TL calculation
354
-
CASE ( 'A', 'a' ) ! arrivals calculation
356
-
CASE ( 'R', 'E' ) ! ray trace
362
14
END SUBROUTINE BellhopCore
364
-
! **********************************************************************!
366
55900*
COMPLEX (KIND=8 ) FUNCTION PickEpsilon( BeamType, omega, c, gradc, alpha, Dalpha, rLoop, EpsMultiplier )
367
-
!! Picks the optimum value for epsilon
369
-
REAL (KIND=8), INTENT( IN ) :: omega, c, gradc( 2 ) ! angular frequency, sound speed and gradient
370
-
REAL (KIND=8), INTENT( IN ) :: alpha, Dalpha ! angular spacing for ray fan
371
-
REAL (KIND=8), INTENT( IN ) :: epsMultiplier, Rloop ! multiplier to manually adjust result, loop range
372
-
CHARACTER (LEN= 2), INTENT( IN ) :: BeamType
373
-
LOGICAL, SAVE :: INIFlag = .TRUE.
374
-
REAL (KIND=8) :: HalfWidth
375
-
REAL (KIND=8) :: cz
376
-
COMPLEX (KIND=8) :: epsilonOpt
377
-
CHARACTER (LEN=40) :: TAG
379
-
! LP: BUG: Multiple codepaths do not set epsilonOpt, leads to UB
381
#####
SELECT CASE ( BeamType( 1 : 1 ) )
383
#####
TAG = 'Paraxial beams'
384
#####
SELECT CASE ( BeamType( 2 : 2 ) )
386
#####
TAG = 'Space filling beams'
387
#####
halfwidth = 2.0 / ( ( omega / c ) * Dalpha )
388
#####
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
390
#####
TAG = 'Minimum width beams'
391
#####
halfwidth = SQRT( 2.0 * c * 1000.0 * rLoop / omega )
392
#####
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
394
#####
TAG = 'WKB beams'
395
#####
halfwidth = HUGE( halfwidth )
396
#####
cz = gradc( 2 )
397
#####
IF ( cz == 0.0 ) THEN
398
#####
epsilonOpt = 1.0D10
400
#####
epsilonOpt = ( -SIN( alpha ) / COS( alpha ** 2 ) ) * c * c / cz
404
-
CASE ( 'G', 'g' ) ! LP: BUG: Missing ^
405
51300
TAG = 'Geometric hat beams'
406
51300
halfwidth = 2.0 / ( ( omega / c ) * Dalpha )
407
51300
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
410
4600
TAG = 'Geometric Gaussian beams'
411
4600
halfwidth = 2.0 / ( ( omega / c ) * Dalpha )
412
4600*
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
415
#####
CALL ERROUT( 'BELLHOP', 'Geo Gaussian beams in ray-centered coords. not implemented in BELLHOP' )
418
#####
TAG = 'Simple Gaussian beams'
419
#####
halfwidth = 2.0 / ( ( omega / c ) * Dalpha )
420
55900*
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
423
55900
PickEpsilon = EpsMultiplier * epsilonOpt
425
-
! On first call write info to prt file
426
55900
IF ( INIFlag ) THEN
427
13
WRITE( PRTFile, * )
428
13
WRITE( PRTFile, * ) TAG
429
13
WRITE( PRTFile, * ) 'halfwidth = ', halfwidth
430
13
WRITE( PRTFile, * ) 'epsilonOpt = ', epsilonOpt
431
13
WRITE( PRTFile, * ) 'EpsMult = ', EpsMultiplier
432
13
WRITE( PRTFile, * )
436
55900
END FUNCTION PickEpsilon
438
-
! **********************************************************************!
440
111900
SUBROUTINE TraceRay2D( xs, alpha, Amp0 )
441
-
!! Traces the beam corresponding to a particular take-off angle
446
-
REAL (KIND=8), INTENT( IN ) :: xs( 2 ) ! x-y coordinate of the source
447
-
REAL (KIND=8), INTENT( IN ) :: alpha, Amp0 ! initial angle, amplitude
448
-
INTEGER :: is, is1 ! index for a step along the ray
449
-
REAL (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho
450
-
REAL (KIND=8) :: dEndTop( 2 ), dEndBot( 2 ), TopnInt( 2 ), BotnInt( 2 ), ToptInt( 2 ), BottInt( 2 )
451
-
REAL (KIND=8) :: DistBegTop, DistEndTop, DistBegBot, DistEndBot ! Distances from ray beginning, end to top and bottom
452
-
REAL (KIND=8) :: sss
453
-
REAL (KIND=8) :: tinit( 2 )
454
-
LOGICAL :: topRefl, botRefl
456
-
! Initial conditions
458
-
! LP: This initialization was missing (only done once globally). In some cases
459
-
! (a source on a boundary), in conjunction with other subtle issues, the
460
-
! missing initialization caused the initial state of one ray to be dependent
461
-
! on the final state of the previous ray, which is obviously non-physical.
465
55950
iSmallStepCtr = 0
466
167850
tinit = [ COS( alpha ), SIN( alpha ) ]
467
55950
CALL EvaluateSSP( xs, tinit, c, cimag, gradc, crr, crz, czz, rho, freq, 'TAB' )
468
55950
ray2D( 1 )%c = c
469
167850
ray2D( 1 )%x = xs
470
167850
ray2D( 1 )%t = tinit / c
471
167850
ray2D( 1 )%p = [ 1.0, 0.0 ]
472
167850
ray2D( 1 )%q = [ 0.0, 1.0 ]
473
55950
ray2D( 1 )%tau = 0.0
474
55950
ray2D( 1 )%Amp = Amp0
475
55950
ray2D( 1 )%Phase = 0.0
476
55950
ray2D( 1 )%NumTopBnc = 0
477
55950
ray2D( 1 )%NumBotBnc = 0
479
-
! second component of qv is not used in geometric beam tracing
480
-
! set I.C. to 0 in hopes of saving run time
481
158650
IF ( Beam%RunType( 2 : 2 ) == 'G' ) ray2D( 1 )%q = [ 0.0, 0.0 ]
483
55950
CALL GetTopSeg( xs( 1 ), ray2D( 1 )%t( 1 ) ) ! identify the top segment above the source
484
55950
CALL GetBotSeg( xs( 1 ), ray2D( 1 )%t( 1 ) ) ! identify the bottom segment below the source
486
-
! convert range-dependent geoacoustic parameters from user to program units
487
-
! compiler is not accepting the copy of the whole structure at once ...
488
55950
IF ( atiType( 2 : 2 ) == 'L' ) THEN
489
#####
Bdry%Top%HS%cp = Top( IsegTop )%HS%cp ! grab the geoacoustic info for the new segment
490
#####
Bdry%Top%HS%cs = Top( IsegTop )%HS%cs
491
#####
Bdry%Top%HS%rho = Top( IsegTop )%HS%rho
494
55950
IF ( btyType( 2 : 2 ) == 'L' ) THEN
495
#####
Bdry%Bot%HS%cp = Bot( IsegBot )%HS%cp
496
#####
Bdry%Bot%HS%cs = Bot( IsegBot )%HS%cs
497
#####
Bdry%Bot%HS%rho = Bot( IsegBot )%HS%rho
500
-
! Trace the beam (note that Reflect alters the step index, is)
502
#####
CALL Distances2D( ray2D( 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop, dEndBot, &
503
55950*
Top( IsegTop )%n, Bot( IsegBot )%n, DistBegTop, DistBegBot )
505
55950
IF ( DistBegTop <= 0 .OR. DistBegBot <= 0 ) THEN
506
#####
Beam%Nsteps = 1
507
#####
WRITE( PRTFile, * ) 'Terminating the ray trace because the source is on or outside the boundaries'
508
#####
RETURN ! source must be within the medium
511
-
! LP: Changed from loop istep to MaxN-1, as is will hit MaxN-3 first.
512
13930471
Stepping: DO WHILE ( .TRUE. )
514
13986421
is1 = is + 1
516
27972842*
CALL Step2D( ray2D( is ), ray2D( is1 ), &
517
-
topRefl, botRefl, &
518
#####
Top( IsegTop )%x, Top( IsegTop )%n, &
519
41959263*
Bot( IsegBot )%x, Bot( IsegBot )%n )
521
-
! New altimetry segment?
522
13986421*
IF ( ray2D( is1 )%x( 1 ) < rTopSeg( 1 ) &
523
27972842*
.OR. ( ray2D( is1 )%x( 1 ) == rTopSeg( 1 ) .AND. ray2D( is1 )%t( 1 ) < 0.0 ) &
524
13986421*
.OR. ray2D( is1 )%x( 1 ) > rTopSeg( 2 ) &
525
69932105*
.OR. ( ray2D( is1 )%x( 1 ) == rTopSeg( 2 ) .AND. ray2D( is1 )%t( 1 ) >= 0.0 )) THEN
526
#####
CALL GetTopSeg( ray2D( is1 )%x( 1 ), ray2D( is1 )%t( 1 ) )
527
#####
IF ( atiType( 2 : 2 ) == 'L' ) THEN
528
#####
Bdry%Top%HS%cp = Top( IsegTop )%HS%cp ! grab the geoacoustic info for the new segment
529
#####
Bdry%Top%HS%cs = Top( IsegTop )%HS%cs
530
#####
Bdry%Top%HS%rho = Top( IsegTop )%HS%rho
534
-
! New bathymetry segment?
535
13986421*
IF ( ray2D( is1 )%x( 1 ) < rBotSeg( 1 ) &
536
27972842*
.OR. ( ray2D( is1 )%x( 1 ) == rBotSeg( 1 ) .AND. ray2D( is1 )%t( 1 ) < 0.0 ) &
537
13986421*
.OR. ray2D( is1 )%x( 1 ) > rBotSeg( 2 ) &
538
69932105*
.OR. ( ray2D( is1 )%x( 1 ) == rBotSeg( 2 ) .AND. ray2D( is1 )%t( 1 ) >= 0.0 )) THEN
539
24353*
CALL GetBotSeg( ray2D( is1 )%x( 1 ), ray2D( is1 )%t( 1 ) )
540
24353
IF ( btyType( 2 : 2 ) == 'L' ) THEN
541
#####
Bdry%Bot%HS%cp = Bot( IsegBot )%HS%cp ! grab the geoacoustic info for the new segment
542
#####
Bdry%Bot%HS%cs = Bot( IsegBot )%HS%cs
543
#####
Bdry%Bot%HS%rho = Bot( IsegBot )%HS%rho
548
-
! Tests that ray at step is is inside, and ray at step is+1 is outside
549
-
! to detect only a crossing from inside to outside
550
-
! DistBeg is the distance at step is, which is saved
551
-
! DistEnd is the distance at step is+1, which needs to be calculated
553
27972842*
CALL Distances2D( ray2D( is1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop, dEndBot, &
554
41959263*
Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )
556
13986421
IF ( topRefl ) THEN
557
-
! WRITE( PRTFile, * ) 'Top reflecting'
558
328819
IF ( atiType == 'C' ) THEN
559
#####
sss = DOT_PRODUCT( dEndTop, Top( IsegTop )%t ) / Top( IsegTop )%Len ! proportional distance along segment
560
#####
TopnInt = ( 1 - sss ) * Top( IsegTop )%Noden + sss * Top( 1 + IsegTop )%Noden
561
#####
ToptInt = ( 1 - sss ) * Top( IsegTop )%Nodet + sss * Top( 1 + IsegTop )%Nodet
563
986457*
TopnInt = Top( IsegTop )%n ! normal is constant in a segment
564
986457*
ToptInt = Top( IsegTop )%t
567
328819*
CALL Reflect2D( is, Bdry%Top%HS, 'TOP', ToptInt, TopnInt, Top( IsegTop )%kappa, RTop, NTopPTS )
568
328819*
ray2D( is + 1 )%NumTopBnc = ray2D( is )%NumTopBnc + 1
570
657638*
CALL Distances2D( ray2D( is + 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop, dEndBot, &
571
986457*
Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )
573
13657602
ELSE IF ( botRefl ) THEN
574
-
! WRITE( PRTFile, * ) 'Bottom reflecting'
575
346748
IF ( btyType == 'C' ) THEN
576
#####
sss = DOT_PRODUCT( dEndBot, Bot( IsegBot )%t ) / Bot( IsegBot )%Len ! proportional distance along segment
577
#####
BotnInt = ( 1 - sss ) * Bot( IsegBot )%Noden + sss * Bot( 1 + IsegBot )%Noden
578
#####
BottInt = ( 1 - sss ) * Bot( IsegBot )%Nodet + sss * Bot( 1 + IsegBot )%Nodet
580
1040244*
BotnInt = Bot( IsegBot )%n ! normal is constant in a segment
581
1040244*
BottInt = Bot( IsegBot )%t
584
346748*
CALL Reflect2D( is, Bdry%Bot%HS, 'BOT', BottInt, BotnInt, Bot( IsegBot )%kappa, RBot, NBotPTS )
585
346748*
ray2D( is + 1 )%NumBotBnc = ray2D( is )%NumBotBnc + 1
586
693496*
CALL Distances2D( ray2D( is + 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop, dEndBot, &
587
1040244*
Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )
591
-
! Has the ray left the box, lost its energy, escaped the boundaries, or exceeded storage limit?
592
-
! There is a test here for when two adjacent ray points are outside the boundaries
593
-
! BELLHOP3D handles this differently using a counter-limit for really small steps.
594
13986421*
IF ( ABS( ray2D( is + 1 )%x( 1 ) ) >= Beam%Box%r .OR. &
595
27972842*
ABS( ray2D( is + 1 )%x( 2 ) ) >= Beam%Box%z .OR. ray2D( is + 1 )%Amp < 0.005 .OR. &
596
55945684
( DistBegTop < 0.0 .AND. DistEndTop < 0.0 ) .OR. &
597
-
( DistBegBot < 0.0 .AND. DistEndBot < 0.0 ) ) THEN
598
-
! ray2D( is + 1 )%t( 1 ) < 0 ) THEN ! this last test kills off a backward traveling ray
599
55944
Beam%Nsteps = is + 1
601
13930477
ELSE IF ( is >= MaxN - 3 ) THEN
602
6
WRITE( PRTFile, * ) 'Warning in TraceRay2D : Insufficient storage for ray trajectory'
607
13930471
DistBegTop = DistEndTop
608
13930471
DistBegBot = DistEndBot
611
-
END SUBROUTINE TraceRay2D
613
-
! **********************************************************************!
615
14717938
SUBROUTINE Distances2D( rayx, Topx, Botx, dTop, dBot, Topn, Botn, DistTop, DistBot )
616
-
!! Calculates the distances to the boundaries
618
-
! Formula differs from JKPS because code uses outward pointing normals
620
-
REAL (KIND=8), INTENT( IN ) :: rayx( 2 ) ! ray coordinate
621
-
REAL (KIND=8), INTENT( IN ) :: Topx( 2 ), Botx( 2 ) ! top, bottom coordinate
622
-
REAL (KIND=8), INTENT( IN ) :: Topn( 2 ), Botn( 2 ) ! top, bottom normal vector (outward)
623
-
REAL (KIND=8), INTENT( OUT ) :: dTop( 2 ), dBot( 2 ) ! vector pointing from top, bottom bdry to ray
624
-
REAL (KIND=8), INTENT( OUT ) :: DistTop, DistBot ! distance (normal to bdry) from the ray to top, bottom boundary
626
44153814
dTop = rayx - Topx ! vector pointing from top to ray
627
44153814
dBot = rayx - Botx ! vector pointing from bottom to ray
628
44153814
DistTop = -DOT_PRODUCT( Topn, dTop )
629
44153814
DistBot = -DOT_PRODUCT( Botn, dBot )
631
14717938
END SUBROUTINE Distances2D
633
-
! **********************************************************************!
635
675567*
SUBROUTINE Reflect2D( is, HS, BotTop, tBdry, nBdry, kappa, RefC, Npts )
637
-
INTEGER, INTENT( IN ) :: Npts
638
-
REAL (KIND=8), INTENT( IN ) :: tBdry( 2 ), nBdry( 2 ) ! Tangent and normal to the boundary
639
-
REAL (KIND=8), INTENT( IN ) :: kappa ! Boundary curvature
640
-
CHARACTER (LEN=3), INTENT( IN ) :: BotTop ! Flag indicating bottom or top reflection
641
-
TYPE( HSInfo ), INTENT( IN ) :: HS ! half-space properties
642
-
TYPE(ReflectionCoef), INTENT( IN ) :: RefC( NPts ) ! reflection coefficient
643
-
INTEGER, INTENT( INOUT ) :: is
645
-
REAL (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! derivatives of sound speed
646
-
REAL (KIND=8) :: RM, RN, Tg, Th, rayt( 2 ), rayn( 2 ), rayt_tilde( 2 ), rayn_tilde( 2 ), cnjump, csjump ! for curvature change
647
-
REAL (KIND=8) :: ck, co, si, cco, ssi, pdelta, rddelta, sddelta, theta_bot ! for beam shift
648
-
COMPLEX (KIND=8) :: kx, kz, kzP, kzS, kzP2, kzS2, mu, f, g, y2, y4, Refl ! for tabulated reflection coef.
649
-
COMPLEX (KIND=8) :: ch, a, b, d, sb, delta, ddelta ! for beam shift
650
-
TYPE(ReflectionCoef) :: RInt
655
2026701*
Tg = DOT_PRODUCT( ray2D( is )%t, tBdry ) ! component of ray tangent, along boundary
656
2026701*
Th = DOT_PRODUCT( ray2D( is )%t, nBdry ) ! component of ray tangent, normal to boundary
658
675567*
ray2D( is1 )%NumTopBnc = ray2D( is )%NumTopBnc
659
675567*
ray2D( is1 )%NumBotBnc = ray2D( is )%NumBotBnc
660
2026701*
ray2D( is1 )%x = ray2D( is )%x
661
2026701*
ray2D( is1 )%t = ray2D( is )%t - 2.0 * Th * nBdry ! changing the ray direction
663
-
! Calculate the change in curvature
664
-
! Based on formulas given by Muller, Geoph. J. R.A.S., 79 (1984).
666
675567*
CALL EvaluateSSP( ray2D( is )%x, ray2D( is )%t, c, cimag, gradc, crr, crz, czz, rho, freq, 'TAB' ) ! just to get c
668
-
! incident unit ray tangent and normal
669
2026701*
rayt = c * ray2D( is )%t ! unit tangent to ray
670
2026701
rayn = [ -rayt( 2 ), rayt( 1 ) ] ! unit normal to ray
672
-
! reflected unit ray tangent and normal (the reflected tangent, normal system has a different orientation)
673
2026701*
rayt_tilde = c * ray2D( is1 )%t ! unit tangent to ray
674
2026701
rayn_tilde = -[ -rayt_tilde( 2 ), rayt_tilde( 1 ) ] ! unit normal to ray
676
675567
RN = 2 * kappa / c ** 2 / Th ! boundary curvature correction
678
-
! get the jumps (this could be simplified, e.g. jump in rayt is roughly 2 * Th * nbdry
679
2026701
cnjump = -DOT_PRODUCT( gradc, rayn_tilde - rayn )
680
2026701
csjump = -DOT_PRODUCT( gradc, rayt_tilde - rayt )
682
675567
IF ( BotTop == 'TOP' ) THEN
683
328819
cnjump = -cnjump ! this is because the (t,n) system of the top boundary has a different sense to the bottom boundary
687
675567
RM = Tg / Th ! this is tan( alpha ) where alpha is the angle of incidence
688
675567
RN = RN + RM * ( 2 * cnjump - RM * csjump ) / c ** 2
690
#####
SELECT CASE ( Beam%Type( 3 : 3 ) )
697
675567*
ray2D( is1 )%c = c
698
675567*
ray2D( is1 )%tau = ray2D( is )%tau
699
2026701*
ray2D( is1 )%p = ray2D( is )%p + ray2D( is )%q * RN
700
2026701*
ray2D( is1 )%q = ray2D( is )%q
702
-
! account for phase change
704
#####
SELECT CASE ( HS%BC )
705
-
CASE ( 'R' ) ! rigid
706
#####
ray2D( is1 )%Amp = ray2D( is )%Amp
707
#####
ray2D( is1 )%Phase = ray2D( is )%Phase
708
-
CASE ( 'V' ) ! vacuum
709
328312*
ray2D( is1 )%Amp = ray2D( is )%Amp
710
328312*
ray2D( is1 )%Phase = ray2D( is )%Phase + pi
711
-
CASE ( 'F' ) ! file
712
115260
RInt%theta = RadDeg * ABS( ATAN2( Th, Tg ) ) ! angle of incidence (relative to normal to bathymetry)
713
115260*
IF ( RInt%theta > 90 ) RInt%theta = 180. - RInt%theta ! reflection coefficient is symmetric about 90 degrees
714
115260*
CALL InterpolateReflectionCoefficient( RInt, RefC, Npts, PRTFile )
715
115260*
ray2D( is1 )%Amp = ray2D( is )%Amp * RInt%R
716
115260*
ray2D( is1 )%Phase = ray2D( is )%Phase + RInt%phi
717
-
CASE ( 'A', 'G' ) ! half-space
718
231995
kx = omega * Tg ! wavenumber in direction parallel to bathymetry
719
231995
kz = omega * Th ! wavenumber in direction perpendicular to bathymetry (in ocean)
721
-
! notation below is a bit misleading
722
-
! kzS, kzP is really what I called gamma in other codes, and differs by a factor of +/- i
723
231995
IF ( REAL( HS%cS ) > 0.0 ) THEN
724
#####
kzS2 = kx ** 2 - ( omega / HS%cS ) ** 2
725
#####
kzP2 = kx ** 2 - ( omega / HS%cP ) ** 2
726
#####
kzS = SQRT( kzS2 )
727
#####
kzP = SQRT( kzP2 )
728
#####
mu = HS%rho * HS%cS ** 2
730
#####
y2 = ( ( kzS2 + kx ** 2 ) ** 2 - 4.0D0 * kzS * kzP * kx ** 2 ) * mu
731
#####
y4 = kzP * ( kx ** 2 - kzS2 )
733
#####
f = omega ** 2 * y4
736
231995
kzP = SQRT( kx ** 2 - ( omega / HS%cP ) ** 2 )
738
-
! Intel and GFortran compilers return different branches of the SQRT for negative reals
739
231995
IF ( REAL( kzP ) == 0.0D0 .AND. AIMAG( kzP ) < 0.0D0 ) kzP = -kzP
744
231995
Refl = - ( rho * f - i * kz * g ) / ( rho * f + i * kz * g ) ! complex reflection coef.
746
463990*
IF ( ABS( Refl ) < 1.0E-5 ) THEN ! kill a ray that has lost its energy in reflection
747
507*
ray2D( is1 )%Amp = 0.0
748
507*
ray2D( is1 )%Phase = ray2D( is )%Phase
750
231488*
ray2D( is1 )%Amp = ABS( Refl ) * ray2D( is )%Amp
751
231488*
ray2D( is1 )%Phase = ray2D( is )%Phase + ATAN2( AIMAG( Refl ), REAL( Refl ) )
753
-
! compute beam-displacement Tindle, Eq. (14)
754
-
! needs a correction to beam-width as well ...
755
-
! IF ( REAL( kz2Sq ) < 0.0 ) THEN
756
-
! rhoW = 1.0 ! density of water
757
-
! rhoWSq = rhoW * rhoW
758
-
! rhoHSSq = rhoHS * rhoHS
759
-
! DELTA = 2 * GK * rhoW * rhoHS * ( kz1Sq - kz2Sq ) /
760
-
! &( kz1 * i * kz2 *
761
-
! &( -rhoWSq * kz2Sq + rhoHSSq * kz1Sq ) )
762
-
! RV( is + 1 ) = RV( is + 1 ) + DELTA
765
231488
if ( Beam%Type( 4 : 4 ) == 'S' ) then ! beam displacement & width change (Seongil's version)
766
#####
ch = ray2D( is )%c / conjg( HS%cP )
767
#####
co = ray2D( is )%t( 1 ) * ray2D( is )%c
768
#####
si = ray2D( is )%t( 2 ) * ray2D( is )%c
769
#####
ck = omega / ray2D( is )%c
771
#####
a = 2 * HS%rho * ( 1 - ch * ch )
772
#####
b = co * co - ch * ch
773
#####
d = HS%rho * HS%rho * si * si + b
778
#####
IF ( si /= 0.0 ) THEN
779
#####
delta = a * co / si / ( ck * sb * d ) ! Do we need an abs() on this???
784
#####
pdelta = real( delta ) / ( ray2D( is )%c / co)
785
-
ddelta = -a / ( ck*sb*d ) - a*cco / ssi / (ck*sb*d) + a*cco / (ck*b*sb*d) &
786
#####
-a*co / si / (ck*sb*d*d) * (2* HS%rho * HS%rho *si*co-2*co*si)
787
#####
rddelta = -real( ddelta )
788
#####
sddelta = rddelta / abs( rddelta )
790
-
! next 3 lines have an update by Diana McCammon to allow a sloping bottom
791
-
! I think the formulas are good, but this won't be reliable because it doesn't have the logic
792
-
! that tracks crossing into new segments after the ray displacement.
794
#####
theta_bot = datan( tBdry( 2 ) / tBdry( 1 )) ! bottom angle
795
#####
ray2D( is1 )%x( 1 ) = ray2D( is1 )%x( 1 ) + real( delta ) * dcos( theta_bot ) ! range displacement
796
#####
ray2D( is1 )%x( 2 ) = ray2D( is1 )%x( 2 ) + real( delta ) * dsin( theta_bot ) ! depth displacement
797
#####
ray2D( is1 )%tau = ray2D( is1 )%tau + pdelta ! phase change
798
#####
ray2D( is1 )%q = ray2D( is1 )%q + sddelta * rddelta * si * c * ray2D( is )%p ! beam-width change
803
#####
WRITE( PRTFile, * ) 'HS%BC = ', HS%BC
804
675567*
CALL ERROUT( 'Reflect2D', 'Unknown boundary condition type' )
807
675567
END SUBROUTINE Reflect2D
809
-
END PROGRAM BELLHOP