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
73
USE ReadEnvironmentBell
29
-
USE SourceReceiverPositions
34
-
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT
38
-
LOGICAL, PARAMETER :: Inline = .FALSE.
40
-
CHARACTER ( LEN=2 ) :: AttenUnit
41
-
CHARACTER ( LEN=80 ) :: FileRoot
45
-
! get the file root for naming all input and output files
46
-
! should add some checks here ...
48
73
CALL GET_COMMAND_ARGUMENT( 1, FileRoot )
49
-
! Open the print file
50
73*
OPEN( UNIT = PRTFile, FILE = TRIM( FileRoot ) // '.prt', STATUS = 'UNKNOWN', IOSTAT = iostat )
52
-
! Read in or otherwise initialize inline all the variables used by BELLHOP
55
-
! NPts, Sigma not used by BELLHOP
56
-
Title = 'BELLHOP- Calibration case with envfil passed as parameters'
58
-
! NMedia variable is not used by BELLHOP
60
-
! *** Boundary information (type of boundary condition and, if a halfspace, then halfspace info)
63
-
Bdry%Top%HS%BC = 'V'
64
-
Bdry%Top%HS%Depth = 0
65
-
Bdry%Bot%HS%Depth = 100
66
-
Bdry%Bot%HS%Opt = 'A_'
67
-
Bdry%Bot%HS%BC = 'A'
68
-
Bdry%Bot%HS%cp = CRCI( 1D20, 1590D0, 0.5D0, freq, freq, AttenUnit, betaPowerLaw, fT ) ! compressional wave speed
69
-
Bdry%Bot%HS%cs = CRCI( 1D20, 0D0, 0D0, freq, freq, AttenUnit, betaPowerLaw, fT ) ! shear wave speed
70
-
Bdry%Bot%HS%rho = 1.2 ! density
72
-
! *** sound speed in the water column ***
74
-
SSP%Type = 'C' ! interpolation method for SSP
75
-
SSP%NPts = 2 ! number of SSP points
76
-
SSP%z( 1 : 2 ) = [ 0, 100 ]
77
-
SSP%c( 1 : 2 ) = [ 1500, 1500 ]
78
-
SSP%cz( 1 : 2 ) = [ 0, 0 ] ! user should really not have to supply this ...
80
-
! *** source and receiver positions ***
86
-
ALLOCATE( Pos%sz( Pos%NSz ), Pos%ws( Pos%NSz ), Pos%isz( Pos%NSz ) )
87
-
ALLOCATE( Pos%rz( Pos%NRz ), Pos%wr( Pos%NRz ), Pos%irz( Pos%NRz ) )
88
-
ALLOCATE( Pos%Rr( Pos%NRr ) )
91
-
!Pos%Rz = [ 0, 50, 100 ]
92
-
!Pos%Rr = 1000. * [ 1, 2, 3, 4, 5 ] ! meters !!!
93
-
Pos%Rz = [ ( jj, jj = 1, Pos%NRz ) ]
94
-
Pos%Rr = 10. * [ ( jj, jj = 1 , Pos%NRr ) ] ! meters !!!
100
-
Beam%Box%r = 5100 ! meters
102
-
Angles%Nalpha = 1789
103
-
! Angles%alpha = [ -80, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80 ] ! -89 89
104
-
Angles%alpha = ( 180. / Angles%Nalpha ) * [ ( jj, jj = 1, Angles%Nalpha ) ] - 90.
106
-
! *** altimetry ***
108
-
ALLOCATE( Top( 2 ), Stat = IAllocStat )
109
-
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', 'Insufficient memory for altimetry data' )
110
-
Top( 1 )%x = [ -sqrt( huge( Top( 1 )%x( 1 ) ) ) / 1.0d5, 0.d0 ]
111
-
Top( 2 )%x = [ sqrt( huge( Top( 1 )%x( 1 ) ) ) / 1.0d5, 0.d0 ]
113
-
CALL ComputeBdryTangentNormal( Top, 'Top' )
115
-
! *** bathymetry ***
117
-
ALLOCATE( Bot( 2 ), Stat = IAllocStat )
118
-
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', 'Insufficient memory for bathymetry data' )
119
-
Bot( 1 )%x = [ -sqrt( huge( Bot( 1 )%x( 1 ) ) ) / 1.0d5, 5000.d0 ]
120
-
Bot( 2 )%x = [ sqrt( huge( Bot( 1 )%x( 1 ) ) ) / 1.0d5, 5000.d0 ]
122
-
CALL ComputeBdryTangentNormal( Bot, 'Bot' )
124
-
ALLOCATE( RBot( 1 ), Stat = IAllocStat ) ! bottom reflection coefficient
125
-
ALLOCATE( RTop( 1 ), Stat = iAllocStat ) ! top reflection coefficient
126
-
NBotPts = 1 ! LP: missing initialization
129
-
! *** Source Beam Pattern ***
131
-
ALLOCATE( SrcBmPat( 2, 2 ), Stat = IAllocStat )
132
-
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP-ReadPat', 'Insufficient memory' )
133
-
SrcBmPat( 1, : ) = [ -180.0, 0.0 ]
134
-
SrcBmPat( 2, : ) = [ 180.0, 0.0 ]
135
-
SrcBmPat( :, 2 ) = 10 ** ( SrcBmPat( :, 2 ) / 20 ) ! convert dB to linear scale !!!
137
73
CALL ReadEnvironment( FileRoot, ThreeD )
138
71
CALL ReadATI( FileRoot, Bdry%Top%HS%Opt( 5 : 5 ), Bdry%Top%HS%Depth, PRTFile ) ! AlTImetry
139
71
CALL ReadBTY( FileRoot, Bdry%Bot%HS%Opt( 2 : 2 ), Bdry%Bot%HS%Depth, PRTFile ) ! BaThYmetry
140
71
CALL ReadReflectionCoefficient( FileRoot, Bdry%Bot%HS%Opt( 1 : 1 ), Bdry%Top%HS%Opt( 2 : 2 ), PRTFile ) ! (top and bottom)
141
71
SBPFlag = Beam%RunType( 3 : 3 )
142
71
CALL ReadPat( FileRoot, PRTFile ) ! Source Beam Pattern
143
-
! dummy bearing angles
145
71*
ALLOCATE( Pos%theta( Pos%Ntheta ), Stat = IAllocStat )
146
71*
Pos%theta( 1 ) = 0.
149
71
CALL OpenOutputFiles( FileRoot, ThreeD )
154
-
! **********************************************************************!
156
71
SUBROUTINE BellhopCore
157
-
!! Core subroutine to run Bellhop algorithm
162
-
INTEGER, PARAMETER :: ArrivalsStorage = 200000000, MinNArr = 10
163
-
INTEGER :: IBPvec( 1 ), ibp, is, iBeamWindow2, Irz1, Irec, NalphaOpt, iSeg
164
-
REAL :: Tstart, Tstop
165
-
REAL (KIND=8) :: Amp0, DalphaOpt, xs( 2 ), tdummy( 2 ), RadMax, s, &
166
-
c, cimag, gradc( 2 ), crr, crz, czz, rho
167
71
COMPLEX, ALLOCATABLE :: U( :, : )
168
-
COMPLEX (KIND=8) :: epsilon
170
71
CALL CPU_TIME( Tstart )
172
71
omega = 2.0 * pi * freq
174
71
IF ( Beam%deltas == 0.0 ) THEN
175
71
Beam%deltas = ( Bdry%Bot%HS%Depth - Bdry%Top%HS%Depth ) / 10.0 ! Automatic step size selection
176
71
WRITE( PRTFile, * )
177
71
WRITE( PRTFile, fmt = '( '' Step length, deltas = '', G11.4, '' m (automatically selected)'' )' ) Beam%deltas
180
270416*
Angles%alpha = DegRad * Angles%alpha ! convert to radians
181
71
Angles%Dalpha = 0.0
182
71
IF ( Angles%Nalpha /= 1 ) THEN
183
71*
Angles%Dalpha = ( Angles%alpha( Angles%Nalpha ) - Angles%alpha( 1 ) ) / ( Angles%Nalpha - 1 )
184
-
! angular spacing between beams
187
-
! convert range-dependent geoacoustic parameters from user to program units
188
71
IF ( atiType( 2 : 2 ) == 'L' ) THEN
189
#####
DO iSeg = 1, NatiPts
190
#####
Top( iSeg )%HS%cp = CRCI( 1D20, Top( iSeg )%HS%alphaR, Top( iSeg )%HS%alphaI, freq, freq, 'W ', &
191
#####
betaPowerLaw, ft ) ! compressional wave speed
192
#####
Top( iSeg )%HS%cs = CRCI( 1D20, Top( iSeg )%HS%betaR, Top( iSeg )%HS%betaI, freq, freq, 'W ', &
193
#####
betaPowerLaw, ft ) ! shear wave speed
197
71
IF ( btyType( 2 : 2 ) == 'L' ) THEN
198
#####
DO iSeg = 1, NbtyPts
199
#####
Bot( iSeg )%HS%cp = CRCI( 1D20, Bot( iSeg )%HS%alphaR, Bot( iSeg )%HS%alphaI, freq, freq, 'W ', &
200
#####
betaPowerLaw, ft ) ! compressional wave speed
201
#####
Bot( iSeg )%HS%cs = CRCI( 1D20, Bot( iSeg )%HS%betaR, Bot( iSeg )%HS%betaI, freq, freq, 'W ', &
202
#####
betaPowerLaw, ft ) ! shear wave speed
206
#####
SELECT CASE ( Beam%RunType( 5 : 5 ) )
208
#####
NRz_per_range = 1 ! irregular grid
210
71
NRz_per_range = Pos%NRz ! rectilinear grid
213
-
! for a TL calculation, allocate space for the pressure matrix
214
16
SELECT CASE ( Beam%RunType( 1 : 1 ) )
215
-
CASE ( 'C', 'S', 'I' ) ! TL calculation
216
16*
ALLOCATE ( U( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
217
16
IF ( IAllocStat /= 0 ) THEN
218
#####
CALL ERROUT( 'BELLHOP', 'Insufficient memory for TL matrix: reduce Nr * NRz' )
220
-
CASE ( 'A', 'a', 'R', 'E' ) ! Arrivals calculation
221
55*
ALLOCATE ( U( 1, 1 ), Stat = IAllocStat ) ! open a dummy variable
223
#####
WRITE( ERROR_UNIT, * ) 'BELLHOP: Unknown RunType for allocation: ', Beam%RunType( 1 : 1 )
224
71*
CALL ERROUT( 'BELLHOP', 'Unknown RunType in allocation' )
227
-
! for an arrivals run, allocate space for arrivals matrices
228
26
SELECT CASE ( Beam%RunType( 1 : 1 ) )
230
26
MaxNArr = MAX( ArrivalsStorage / ( NRz_per_range * Pos%NRr ), MinNArr ) ! allow space for at least 10 arrivals
231
26
WRITE( PRTFile, * )
232
26
WRITE( PRTFile, * ) '( Maximum # of arrivals = ', MaxNArr, ')'
234
26*
ALLOCATE ( Arr( NRz_per_range, Pos%NRr, MaxNArr ), NArr( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
235
26
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', &
236
26*
'Insufficient memory to allocate arrivals matrix; reduce parameter ArrivalsStorage' )
239
71*
ALLOCATE ( Arr( NRz_per_range, Pos%NRr, 1 ), NArr( NRz_per_range, Pos%NRr ), Stat = IAllocStat )
242
71
WRITE( PRTFile, * )
244
142*
SourceDepth: DO is = 1, Pos%NSz
245
213*
xs = [ 0.0, Pos%sz( is ) ] ! source coordinate
247
16
SELECT CASE ( Beam%RunType( 1 : 1 ) )
248
-
CASE ( 'C', 'S', 'I' ) ! TL calculation, zero out pressure matrix
250
-
CASE ( 'A', 'a' ) ! Arrivals calculation, zero out arrival matrix
252
-
CASE ( 'R', 'E' ) ! Ray trace or eigenrays
255
#####
WRITE( ERROR_UNIT, * ) 'BELLHOP: Unknown RunType for initialization: ', Beam%RunType( 1 : 1 )
256
71*
CALL ERROUT( 'BELLHOP', 'Unknown RunType in initialization' )
259
-
! LP: This initialization was missing (only done once globally). In some
260
-
! cases (a source on a boundary), in conjunction with other subtle issues,
261
-
! the missing initialization caused the initial state of one ray to be
262
-
! dependent on the final state of the previous ray, which is obviously
264
-
! This initialization is here in addition to TraceRay2D because of EvaluateSSP.
268
71
tdummy = [ 1.0, 0.0 ]
269
71
CALL EvaluateSSP( xs, tdummy, c, cimag, gradc, crr, crz, czz, rho, freq, 'TAB' )
270
71
RadMax = 5 * c / freq ! 5 wavelength max radius
272
-
! Are there enough beams?
273
71*
DalphaOpt = SQRT( c / ( 6.0 * freq * Pos%Rr( Pos%NRr ) ) )
274
71*
NalphaOpt = 2 + INT( ( Angles%alpha( Angles%Nalpha ) - Angles%alpha( 1 ) ) / DalphaOpt )
276
71
IF ( Beam%RunType( 1 : 1 ) == 'C' .AND. Angles%Nalpha < NalphaOpt ) THEN
277
4
WRITE( PRTFile, * ) 'Warning in BELLHOP : Too few beams'
278
4
WRITE( PRTFile, * ) 'Nalpha should be at least = ', NalphaOpt
281
-
! Trace successive beams
283
270413*
DeclinationAngle: DO ialpha = 1, Angles%Nalpha
284
270342*
SrcDeclAngle = RadDeg * Angles%alpha( ialpha ) ! take-off declination angle in degrees
286
270413
IF ( Angles%iSingle_alpha == 0 .OR. ialpha == Angles%iSingle_alpha ) THEN ! Single beam run?
288
823614*
IBPvec = maxloc( SrcBmPat( :, 1 ), mask = SrcBmPat( :, 1 ) < SrcDeclAngle ) ! index of ray angle in beam pattern
289
270338
IBP = IBPvec( 1 )
290
270338
IBP = MAX( IBP, 1 ) ! don't go before beginning of table
291
270338
IBP = MIN( IBP, NSBPPts - 1 ) ! don't go past end of table
293
-
! linear interpolation to get amplitude
294
270338*
s = ( SrcDeclAngle - SrcBmPat( IBP, 1 ) ) / ( SrcBmPat( IBP + 1, 1 ) - SrcBmPat( IBP, 1 ) )
295
270338*
Amp0 = ( 1 - s ) * SrcBmPat( IBP, 2 ) + s * SrcBmPat( IBP + 1, 2 )
297
-
! Lloyd mirror pattern for semi-coherent option
298
270338
IF ( Beam%RunType( 1 : 1 ) == 'S' ) THEN
299
#####
Amp0 = Amp0 * SQRT( 2.0 ) * ABS( SIN( omega / c * xs( 2 ) * SIN( Angles%alpha( ialpha ) ) ) )
302
-
! show progress ...
303
270338
IF ( MOD( ialpha - 1, max( Angles%Nalpha / 50, 1 ) ) == 0 ) THEN
304
3509
WRITE( PRTFile, FMT = "( 'Tracing beam ', I7, F10.2 )" ) ialpha, SrcDeclAngle
305
3509
FLUSH( PRTFile )
308
270338*
CALL TraceRay2D( xs, Angles%alpha( ialpha ), Amp0 ) ! Trace a ray
310
270338
IF ( Beam%RunType( 1 : 1 ) == 'R' ) THEN ! Write the ray trajectory to RAYFile
311
601
CALL WriteRay2D( SrcDeclAngle, Beam%Nsteps )
312
-
ELSE ! Compute the contribution to the field
314
#####
epsilon = PickEpsilon( Beam%Type( 1 : 2 ), omega, c, gradc, Angles%alpha( ialpha ), &
315
269737
Angles%Dalpha, Beam%rLoop, Beam%epsMultiplier ) ! 'optimal' beam constant
317
#####
SELECT CASE ( Beam%Type( 1 : 1 ) )
319
#####
iBeamWindow2 = Beam%iBeamWindow **2
320
#####
RadMax = 50 * c / freq ! 50 wavelength max radius
321
#####
CALL InfluenceCervenyRayCen( U, epsilon, Angles%alpha( ialpha ), iBeamWindow2, RadMax )
323
#####
iBeamWindow2 = Beam%iBeamWindow **2
324
#####
RadMax = 50 * c / freq ! 50 wavelength max radius
325
#####
CALL InfluenceCervenyCart( U, epsilon, Angles%alpha( ialpha ), iBeamWindow2, RadMax )
327
#####
CALL InfluenceGeoHatRayCen( U, Angles%alpha( ialpha ), Angles%Dalpha )
329
#####
RadMax = 50 * c / freq ! 50 wavelength max radius
330
#####
CALL InfluenceSGB( U, Angles%alpha( ialpha ), Angles%Dalpha, RadMax )
332
9200*
CALL InfluenceGeoGaussianCart( U, Angles%alpha( ialpha ), Angles%Dalpha )
334
278937*
CALL InfluenceGeoHatCart( U, Angles%alpha( ialpha ), Angles%Dalpha )
339
-
END DO DeclinationAngle
341
-
! write results to disk
343
142
SELECT CASE ( Beam%RunType( 1 : 1 ) )
344
-
CASE ( 'C', 'S', 'I' ) ! TL calculation
345
16*
CALL ScalePressure( Angles%Dalpha, ray2D( 1 )%c, Pos%Rr, U, NRz_per_range, Pos%NRr, Beam%RunType, freq )
346
16
IRec = 10 + NRz_per_range * ( is - 1 )
347
2637*
RcvrDepth: DO Irz1 = 1, NRz_per_range
349
2621*
WRITE( SHDFile, REC = IRec ) U( Irz1, 1 : Pos%NRr )
351
-
CASE ( 'A' ) ! arrivals calculation, ascii
352
26*
CALL WriteArrivalsASCII( Pos%Rr, NRz_per_range, Pos%NRr, Beam%RunType( 4 : 4 ) )
353
-
CASE ( 'a' ) ! arrivals calculation, binary
354
#####
CALL WriteArrivalsBinary( Pos%Rr, NRz_per_range, Pos%NRr, Beam%RunType( 4 : 4 ) )
355
-
CASE ( 'R', 'E' ) ! ray trace or eigenrays
358
#####
WRITE( ERROR_UNIT, * ) 'BELLHOP: Unknown RunType when writing results: ', Beam%RunType( 1 : 1 )
359
113*
CALL ERROUT( 'BELLHOP', 'Unknown RunType when writing results' )
365
71
CALL CPU_TIME( Tstop )
366
71
WRITE( PRTFile, "( /, ' CPU Time = ', G15.3, 's' )" ) Tstop - Tstart
369
16
SELECT CASE ( Beam%RunType( 1 : 1 ) )
370
-
CASE ( 'C', 'S', 'I' ) ! TL calculation
372
-
CASE ( 'A', 'a' ) ! arrivals calculation
374
-
CASE ( 'R', 'E' ) ! ray trace
377
#####
WRITE( ERROR_UNIT, * ) 'BELLHOP: Unknown RunType when closing files: ', Beam%RunType( 1 : 1 )
378
71*
CALL ERROUT( 'BELLHOP', 'Unknown RunType when closing files' )
383
71
END SUBROUTINE BellhopCore
385
-
! **********************************************************************!
387
269737*
COMPLEX (KIND=8 ) FUNCTION PickEpsilon( BeamType, omega, c, gradc, alpha, Dalpha, rLoop, EpsMultiplier )
388
-
!! Picks the optimum value for epsilon
390
-
REAL (KIND=8), INTENT( IN ) :: omega, c, gradc( 2 ) ! angular frequency, sound speed and gradient
391
-
REAL (KIND=8), INTENT( IN ) :: alpha, Dalpha ! angular spacing for ray fan
392
-
REAL (KIND=8), INTENT( IN ) :: epsMultiplier, Rloop ! multiplier to manually adjust result, loop range
393
-
CHARACTER (LEN= 2), INTENT( IN ) :: BeamType
394
-
LOGICAL, SAVE :: INIFlag = .TRUE.
395
-
REAL (KIND=8) :: HalfWidth
396
-
REAL (KIND=8) :: cz
397
-
COMPLEX (KIND=8) :: epsilonOpt
398
-
CHARACTER (LEN=40) :: TAG
400
-
! LP: BUG: Multiple codepaths do not set epsilonOpt, leads to UB
402
#####
SELECT CASE ( BeamType( 1 : 1 ) )
404
#####
TAG = 'Paraxial beams'
405
#####
SELECT CASE ( BeamType( 2 : 2 ) )
407
#####
TAG = 'Space filling beams'
408
#####
halfwidth = 2.0 / ( ( omega / c ) * Dalpha )
409
#####
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
411
#####
TAG = 'Minimum width beams'
412
#####
halfwidth = SQRT( 2.0 * c * 1000.0 * rLoop / omega )
413
#####
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
415
#####
TAG = 'WKB beams'
416
#####
halfwidth = HUGE( halfwidth )
417
#####
cz = gradc( 2 )
418
#####
IF ( cz == 0.0 ) THEN
419
#####
epsilonOpt = 1.0D10
421
#####
epsilonOpt = ( -SIN( alpha ) / COS( alpha ** 2 ) ) * c * c / cz
424
#####
WRITE( ERROR_UNIT, * ) 'BELLHOP: Unknown beam width option: ', BeamType( 2 : 2 )
425
#####
CALL ERROUT( 'BELLHOP', 'Unknown beam width option' )
428
-
CASE ( 'G', 'g', '^' ) ! LP: BUG: Missing ^
429
260537
TAG = 'Geometric hat beams'
430
260537
halfwidth = 2.0 / ( ( omega / c ) * Dalpha )
431
260537
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
434
9200
TAG = 'Geometric Gaussian beams'
435
9200
halfwidth = 2.0 / ( ( omega / c ) * Dalpha )
436
9200*
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
439
#####
CALL ERROUT( 'BELLHOP', 'Geo Gaussian beams in ray-centered coords. not implemented in BELLHOP' )
442
#####
TAG = 'Simple Gaussian beams'
443
#####
halfwidth = 2.0 / ( ( omega / c ) * Dalpha )
444
#####
epsilonOpt = i * 0.5 * omega * halfwidth ** 2
446
#####
WRITE( ERROR_UNIT, * ) 'BELLHOP: Unknown beam type: ', BeamType( 1 : 1 )
447
269737*
CALL ERROUT( 'BELLHOP', 'Unknown beam type' )
450
269737
PickEpsilon = EpsMultiplier * epsilonOpt
452
-
! On first call write info to prt file
453
269737
IF ( INIFlag ) THEN
454
60
WRITE( PRTFile, * )
455
60
WRITE( PRTFile, * ) TAG
456
60
WRITE( PRTFile, * ) 'halfwidth = ', halfwidth
457
60
WRITE( PRTFile, * ) 'epsilonOpt = ', epsilonOpt
458
60
WRITE( PRTFile, * ) 'EpsMult = ', EpsMultiplier
459
60
WRITE( PRTFile, * )
463
269737
END FUNCTION PickEpsilon
465
-
! **********************************************************************!
467
540676
SUBROUTINE TraceRay2D( xs, alpha, Amp0 )
468
-
!! Traces the beam corresponding to a particular take-off angle
473
-
REAL (KIND=8), INTENT( IN ) :: xs( 2 ) ! x-y coordinate of the source
474
-
REAL (KIND=8), INTENT( IN ) :: alpha, Amp0 ! initial angle, amplitude
475
-
INTEGER :: is, is1 ! index for a step along the ray
476
-
REAL (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho
477
-
REAL (KIND=8) :: dEndTop( 2 ), dEndBot( 2 ), TopnInt( 2 ), BotnInt( 2 ), ToptInt( 2 ), BottInt( 2 )
478
-
REAL (KIND=8) :: DistBegTop, DistEndTop, DistBegBot, DistEndBot ! Distances from ray beginning, end to top and bottom
479
-
REAL (KIND=8) :: sss
480
-
REAL (KIND=8) :: tinit( 2 )
481
-
LOGICAL :: topRefl, botRefl
483
-
! Initial conditions
485
-
! LP: This initialization was missing (only done once globally). In some cases
486
-
! (a source on a boundary), in conjunction with other subtle issues, the
487
-
! missing initialization caused the initial state of one ray to be dependent
488
-
! on the final state of the previous ray, which is obviously non-physical.
492
270338
iSmallStepCtr = 0
493
811014
tinit = [ COS( alpha ), SIN( alpha ) ]
494
270338
CALL EvaluateSSP( xs, tinit, c, cimag, gradc, crr, crz, czz, rho, freq, 'TAB' )
495
270338
ray2D( 1 )%c = c
496
811014
ray2D( 1 )%x = xs
497
811014
ray2D( 1 )%t = tinit / c
498
811014
ray2D( 1 )%p = [ 1.0, 0.0 ]
499
811014
ray2D( 1 )%q = [ 0.0, 1.0 ]
500
270338
ray2D( 1 )%tau = 0.0
501
270338
ray2D( 1 )%Amp = Amp0
502
270338
ray2D( 1 )%Phase = 0.0
503
270338
ray2D( 1 )%NumTopBnc = 0
504
270338
ray2D( 1 )%NumBotBnc = 0
506
-
! second component of qv is not used in geometric beam tracing
507
-
! set I.C. to 0 in hopes of saving run time
508
792470
IF ( Beam%RunType( 2 : 2 ) == 'G' ) ray2D( 1 )%q = [ 0.0, 0.0 ]
510
270338
CALL GetTopSeg( xs( 1 ), ray2D( 1 )%t( 1 ) ) ! identify the top segment above the source
511
270338
CALL GetBotSeg( xs( 1 ), ray2D( 1 )%t( 1 ) ) ! identify the bottom segment below the source
513
-
! convert range-dependent geoacoustic parameters from user to program units
514
-
! compiler is not accepting the copy of the whole structure at once ...
515
270338
IF ( atiType( 2 : 2 ) == 'L' ) THEN
516
#####
Bdry%Top%HS%cp = Top( IsegTop )%HS%cp ! grab the geoacoustic info for the new segment
517
#####
Bdry%Top%HS%cs = Top( IsegTop )%HS%cs
518
#####
Bdry%Top%HS%rho = Top( IsegTop )%HS%rho
521
270338
IF ( btyType( 2 : 2 ) == 'L' ) THEN
522
#####
Bdry%Bot%HS%cp = Bot( IsegBot )%HS%cp
523
#####
Bdry%Bot%HS%cs = Bot( IsegBot )%HS%cs
524
#####
Bdry%Bot%HS%rho = Bot( IsegBot )%HS%rho
527
-
! Trace the beam (note that Reflect alters the step index, is)
529
#####
CALL Distances2D( ray2D( 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop, dEndBot, &
530
270338*
Top( IsegTop )%n, Bot( IsegBot )%n, DistBegTop, DistBegBot )
532
270338
IF ( DistBegTop <= 0 .OR. DistBegBot <= 0 ) THEN
533
#####
Beam%Nsteps = 1
534
#####
WRITE( PRTFile, * ) 'Terminating the ray trace because the source is on or outside the boundaries'
535
#####
RETURN ! source must be within the medium
538
-
! LP: Changed from loop istep to MaxN-1, as is will hit MaxN-3 first.
539
34092810
Stepping: DO WHILE ( .TRUE. )
541
34363148
is1 = is + 1
543
68726296*
CALL Step2D( ray2D( is ), ray2D( is1 ), &
544
-
topRefl, botRefl, &
545
#####
Top( IsegTop )%x, Top( IsegTop )%n, &
546
103089444*
Bot( IsegBot )%x, Bot( IsegBot )%n )
548
-
! New altimetry segment?
549
34363148*
IF ( ray2D( is1 )%x( 1 ) < rTopSeg( 1 ) &
550
68726296*
.OR. ( ray2D( is1 )%x( 1 ) == rTopSeg( 1 ) .AND. ray2D( is1 )%t( 1 ) < 0.0 ) &
551
34363148*
.OR. ray2D( is1 )%x( 1 ) > rTopSeg( 2 ) &
552
171815740*
.OR. ( ray2D( is1 )%x( 1 ) == rTopSeg( 2 ) .AND. ray2D( is1 )%t( 1 ) >= 0.0 )) THEN
553
949800*
CALL GetTopSeg( ray2D( is1 )%x( 1 ), ray2D( is1 )%t( 1 ) )
554
949800
IF ( atiType( 2 : 2 ) == 'L' ) THEN
555
#####
Bdry%Top%HS%cp = Top( IsegTop )%HS%cp ! grab the geoacoustic info for the new segment
556
#####
Bdry%Top%HS%cs = Top( IsegTop )%HS%cs
557
#####
Bdry%Top%HS%rho = Top( IsegTop )%HS%rho
561
-
! New bathymetry segment?
562
34363148*
IF ( ray2D( is1 )%x( 1 ) < rBotSeg( 1 ) &
563
68726296*
.OR. ( ray2D( is1 )%x( 1 ) == rBotSeg( 1 ) .AND. ray2D( is1 )%t( 1 ) < 0.0 ) &
564
34363148*
.OR. ray2D( is1 )%x( 1 ) > rBotSeg( 2 ) &
565
171815740*
.OR. ( ray2D( is1 )%x( 1 ) == rBotSeg( 2 ) .AND. ray2D( is1 )%t( 1 ) >= 0.0 )) THEN
566
94533*
CALL GetBotSeg( ray2D( is1 )%x( 1 ), ray2D( is1 )%t( 1 ) )
567
94533
IF ( btyType( 2 : 2 ) == 'L' ) THEN
568
#####
Bdry%Bot%HS%cp = Bot( IsegBot )%HS%cp ! grab the geoacoustic info for the new segment
569
#####
Bdry%Bot%HS%cs = Bot( IsegBot )%HS%cs
570
#####
Bdry%Bot%HS%rho = Bot( IsegBot )%HS%rho
575
-
! Tests that ray at step is is inside, and ray at step is+1 is outside
576
-
! to detect only a crossing from inside to outside
577
-
! DistBeg is the distance at step is, which is saved
578
-
! DistEnd is the distance at step is+1, which needs to be calculated
580
68726296*
CALL Distances2D( ray2D( is1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop, dEndBot, &
581
103089444*
Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )
583
34363148
IF ( topRefl ) THEN
584
-
! WRITE( PRTFile, * ) 'Top reflecting'
585
568025
IF ( atiType == 'C' ) THEN
586
#####
sss = DOT_PRODUCT( dEndTop, Top( IsegTop )%t ) / Top( IsegTop )%Len ! proportional distance along segment
587
#####
TopnInt = ( 1 - sss ) * Top( IsegTop )%Noden + sss * Top( 1 + IsegTop )%Noden
588
#####
ToptInt = ( 1 - sss ) * Top( IsegTop )%Nodet + sss * Top( 1 + IsegTop )%Nodet
590
1704075*
TopnInt = Top( IsegTop )%n ! normal is constant in a segment
591
1704075*
ToptInt = Top( IsegTop )%t
594
568025*
CALL Reflect2D( is, Bdry%Top%HS, 'TOP', ToptInt, TopnInt, Top( IsegTop )%kappa, RTop, NTopPTS )
595
568025*
ray2D( is + 1 )%NumTopBnc = ray2D( is )%NumTopBnc + 1
597
1136050*
CALL Distances2D( ray2D( is + 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop, dEndBot, &
598
1704075*
Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )
600
33795123
ELSE IF ( botRefl ) THEN
601
-
! WRITE( PRTFile, * ) 'Bottom reflecting'
602
684579
IF ( btyType == 'C' ) THEN
603
#####
sss = DOT_PRODUCT( dEndBot, Bot( IsegBot )%t ) / Bot( IsegBot )%Len ! proportional distance along segment
604
#####
BotnInt = ( 1 - sss ) * Bot( IsegBot )%Noden + sss * Bot( 1 + IsegBot )%Noden
605
#####
BottInt = ( 1 - sss ) * Bot( IsegBot )%Nodet + sss * Bot( 1 + IsegBot )%Nodet
607
2053737*
BotnInt = Bot( IsegBot )%n ! normal is constant in a segment
608
2053737*
BottInt = Bot( IsegBot )%t
611
684579*
CALL Reflect2D( is, Bdry%Bot%HS, 'BOT', BottInt, BotnInt, Bot( IsegBot )%kappa, RBot, NBotPTS )
612
684579*
ray2D( is + 1 )%NumBotBnc = ray2D( is )%NumBotBnc + 1
613
1369158*
CALL Distances2D( ray2D( is + 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop, dEndBot, &
614
2053737*
Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )
618
-
! Has the ray left the box, lost its energy, escaped the boundaries, or exceeded storage limit?
619
-
! There is a test here for when two adjacent ray points are outside the boundaries
620
-
! BELLHOP3D handles this differently using a counter-limit for really small steps.
621
34363148*
IF ( ABS( ray2D( is + 1 )%x( 1 ) ) >= Beam%Box%r .OR. &
622
68726296*
ABS( ray2D( is + 1 )%x( 2 ) ) >= Beam%Box%z .OR. ray2D( is + 1 )%Amp < 0.005 .OR. &
623
137452592
( DistBegTop < 0.0 .AND. DistEndTop < 0.0 ) .OR. &
624
-
( DistBegBot < 0.0 .AND. DistEndBot < 0.0 ) ) THEN
625
-
! ray2D( is + 1 )%t( 1 ) < 0 ) THEN ! this last test kills off a backward traveling ray
626
270332
Beam%Nsteps = is + 1
628
34092816
ELSE IF ( is >= MaxN - 3 ) THEN
629
6
WRITE( PRTFile, * ) 'Warning in TraceRay2D : Insufficient storage for ray trajectory'
634
34092810
DistBegTop = DistEndTop
635
34092810
DistBegBot = DistEndBot
638
-
END SUBROUTINE TraceRay2D
640
-
! **********************************************************************!
642
35886090
SUBROUTINE Distances2D( rayx, Topx, Botx, dTop, dBot, Topn, Botn, DistTop, DistBot )
643
-
!! Calculates the distances to the boundaries
645
-
! Formula differs from JKPS because code uses outward pointing normals
647
-
REAL (KIND=8), INTENT( IN ) :: rayx( 2 ) ! ray coordinate
648
-
REAL (KIND=8), INTENT( IN ) :: Topx( 2 ), Botx( 2 ) ! top, bottom coordinate
649
-
REAL (KIND=8), INTENT( IN ) :: Topn( 2 ), Botn( 2 ) ! top, bottom normal vector (outward)
650
-
REAL (KIND=8), INTENT( OUT ) :: dTop( 2 ), dBot( 2 ) ! vector pointing from top, bottom bdry to ray
651
-
REAL (KIND=8), INTENT( OUT ) :: DistTop, DistBot ! distance (normal to bdry) from the ray to top, bottom boundary
653
107658270
dTop = rayx - Topx ! vector pointing from top to ray
654
107658270
dBot = rayx - Botx ! vector pointing from bottom to ray
655
107658270
DistTop = -DOT_PRODUCT( Topn, dTop )
656
107658270
DistBot = -DOT_PRODUCT( Botn, dBot )
658
35886090
END SUBROUTINE Distances2D
660
-
! **********************************************************************!
662
1252604*
SUBROUTINE Reflect2D( is, HS, BotTop, tBdry, nBdry, kappa, RefC, Npts )
664
-
INTEGER, INTENT( IN ) :: Npts
665
-
REAL (KIND=8), INTENT( IN ) :: tBdry( 2 ), nBdry( 2 ) ! Tangent and normal to the boundary
666
-
REAL (KIND=8), INTENT( IN ) :: kappa ! Boundary curvature
667
-
CHARACTER (LEN=3), INTENT( IN ) :: BotTop ! Flag indicating bottom or top reflection
668
-
TYPE( HSInfo ), INTENT( IN ) :: HS ! half-space properties
669
-
TYPE(ReflectionCoef), INTENT( IN ) :: RefC( NPts ) ! reflection coefficient
670
-
INTEGER, INTENT( INOUT ) :: is
672
-
REAL (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! derivatives of sound speed
673
-
REAL (KIND=8) :: RM, RN, Tg, Th, rayt( 2 ), rayn( 2 ), rayt_tilde( 2 ), rayn_tilde( 2 )
674
-
REAL (KIND=8) :: cnjump, csjump ! for curvature change
675
-
REAL (KIND=8) :: ck, co, si, cco, ssi, pdelta, rddelta, sddelta, theta_bot ! for beam shift
676
-
COMPLEX (KIND=8) :: kx, kz, kzP, kzS, kzP2, kzS2, mu, f, g, y2, y4, Refl ! for tabulated reflection coef.
677
-
COMPLEX (KIND=8) :: ch, a, b, d, sb, delta, ddelta ! for beam shift
678
-
TYPE(ReflectionCoef) :: RInt
683
3757812*
Tg = DOT_PRODUCT( ray2D( is )%t, tBdry ) ! component of ray tangent, along boundary
684
3757812*
Th = DOT_PRODUCT( ray2D( is )%t, nBdry ) ! component of ray tangent, normal to boundary
686
1252604*
ray2D( is1 )%NumTopBnc = ray2D( is )%NumTopBnc
687
1252604*
ray2D( is1 )%NumBotBnc = ray2D( is )%NumBotBnc
688
3757812*
ray2D( is1 )%x = ray2D( is )%x
689
3757812*
ray2D( is1 )%t = ray2D( is )%t - 2.0 * Th * nBdry ! changing the ray direction
691
-
! Calculate the change in curvature
692
-
! Based on formulas given by Muller, Geoph. J. R.A.S., 79 (1984).
694
1252604*
CALL EvaluateSSP( ray2D( is )%x, ray2D( is )%t, c, cimag, gradc, crr, crz, czz, rho, freq, 'TAB' ) ! just to get c
696
-
! incident unit ray tangent and normal
697
3757812*
rayt = c * ray2D( is )%t ! unit tangent to ray
698
3757812
rayn = [ -rayt( 2 ), rayt( 1 ) ] ! unit normal to ray
700
-
! reflected unit ray tangent and normal (the reflected tangent, normal system has a different orientation)
701
3757812*
rayt_tilde = c * ray2D( is1 )%t ! unit tangent to ray
702
3757812
rayn_tilde = -[ -rayt_tilde( 2 ), rayt_tilde( 1 ) ] ! unit normal to ray
704
1252604
RN = 2 * kappa / c ** 2 / Th ! boundary curvature correction
706
-
! get the jumps (this could be simplified, e.g. jump in rayt is roughly 2 * Th * nbdry
707
3757812
cnjump = -DOT_PRODUCT( gradc, rayn_tilde - rayn )
708
3757812
csjump = -DOT_PRODUCT( gradc, rayt_tilde - rayt )
710
1252604
IF ( BotTop == 'TOP' ) THEN
711
568025
cnjump = -cnjump ! this is because the (t,n) system of the top boundary has a different sense to the bottom boundary
715
1252604
RM = Tg / Th ! this is tan( alpha ) where alpha is the angle of incidence
716
1252604
RN = RN + RM * ( 2 * cnjump - RM * csjump ) / c ** 2
718
#####
SELECT CASE ( Beam%Type( 3 : 3 ) )
726
#####
WRITE( ERROR_UNIT, * ) 'BELLHOP: Unknown curvature option: ', Beam%Type( 3 : 3 )
727
1252604*
CALL ERROUT( 'BELLHOP', 'Unknown curvature option' )
730
1252604*
ray2D( is1 )%c = c
731
1252604*
ray2D( is1 )%tau = ray2D( is )%tau
732
3757812*
ray2D( is1 )%p = ray2D( is )%p + ray2D( is )%q * RN
733
3757812*
ray2D( is1 )%q = ray2D( is )%q
735
-
! account for phase change
737
#####
SELECT CASE ( HS%BC )
738
-
CASE ( 'R' ) ! rigid
739
#####
ray2D( is1 )%Amp = ray2D( is )%Amp
740
#####
ray2D( is1 )%Phase = ray2D( is )%Phase
741
-
CASE ( 'V' ) ! vacuum
742
563397*
ray2D( is1 )%Amp = ray2D( is )%Amp
743
563397*
ray2D( is1 )%Phase = ray2D( is )%Phase + pi
744
-
CASE ( 'F' ) ! file
745
143339
RInt%theta = RadDeg * ABS( ATAN2( Th, Tg ) ) ! angle of incidence (relative to normal to bathymetry)
746
143339*
IF ( RInt%theta > 90 ) RInt%theta = 180. - RInt%theta ! reflection coefficient is symmetric about 90 degrees
747
143339*
CALL InterpolateReflectionCoefficient( RInt, RefC, Npts, PRTFile )
748
143339*
ray2D( is1 )%Amp = ray2D( is )%Amp * RInt%R
749
143339*
ray2D( is1 )%Phase = ray2D( is )%Phase + RInt%phi
750
-
CASE ( 'A', 'G' ) ! half-space
751
545868
kx = omega * Tg ! wavenumber in direction parallel to bathymetry
752
545868
kz = omega * Th ! wavenumber in direction perpendicular to bathymetry (in ocean)
754
-
! notation below is a bit misleading
755
-
! kzS, kzP is really what I called gamma in other codes, and differs by a factor of +/- i
756
545868
IF ( REAL( HS%cS ) > 0.0 ) THEN
757
#####
kzS2 = kx ** 2 - ( omega / HS%cS ) ** 2
758
#####
kzP2 = kx ** 2 - ( omega / HS%cP ) ** 2
759
#####
kzS = SQRT( kzS2 )
760
#####
kzP = SQRT( kzP2 )
761
#####
mu = HS%rho * HS%cS ** 2
763
#####
y2 = ( ( kzS2 + kx ** 2 ) ** 2 - 4.0D0 * kzS * kzP * kx ** 2 ) * mu
764
#####
y4 = kzP * ( kx ** 2 - kzS2 )
766
#####
f = omega ** 2 * y4
769
545868
kzP = SQRT( kx ** 2 - ( omega / HS%cP ) ** 2 )
771
-
! Intel and GFortran compilers return different branches of the SQRT for negative reals
772
545868
IF ( REAL( kzP ) == 0.0D0 .AND. AIMAG( kzP ) < 0.0D0 ) kzP = -kzP
777
545868
Refl = - ( rho * f - i * kz * g ) / ( rho * f + i * kz * g ) ! complex reflection coef.
779
1091736*
IF ( ABS( Refl ) < 1.0E-5 ) THEN ! kill a ray that has lost its energy in reflection
780
968*
ray2D( is1 )%Amp = 0.0
781
968*
ray2D( is1 )%Phase = ray2D( is )%Phase
783
544900*
ray2D( is1 )%Amp = ABS( Refl ) * ray2D( is )%Amp
784
544900*
ray2D( is1 )%Phase = ray2D( is )%Phase + ATAN2( AIMAG( Refl ), REAL( Refl ) )
786
-
! compute beam-displacement Tindle, Eq. (14)
787
-
! needs a correction to beam-width as well ...
788
-
! IF ( REAL( kz2Sq ) < 0.0 ) THEN
789
-
! rhoW = 1.0 ! density of water
790
-
! rhoWSq = rhoW * rhoW
791
-
! rhoHSSq = rhoHS * rhoHS
792
-
! DELTA = 2 * GK * rhoW * rhoHS * ( kz1Sq - kz2Sq ) /
793
-
! &( kz1 * i * kz2 *
794
-
! &( -rhoWSq * kz2Sq + rhoHSSq * kz1Sq ) )
795
-
! RV( is + 1 ) = RV( is + 1 ) + DELTA
798
544900
if ( Beam%Type( 4 : 4 ) == 'S' ) then ! beam displacement & width change (Seongil's version)
799
#####
ch = ray2D( is )%c / conjg( HS%cP )
800
#####
co = ray2D( is )%t( 1 ) * ray2D( is )%c
801
#####
si = ray2D( is )%t( 2 ) * ray2D( is )%c
802
#####
ck = omega / ray2D( is )%c
804
#####
a = 2 * HS%rho * ( 1 - ch * ch )
805
#####
b = co * co - ch * ch
806
#####
d = HS%rho * HS%rho * si * si + b
811
#####
IF ( si /= 0.0 ) THEN
812
#####
delta = a * co / si / ( ck * sb * d ) ! Do we need an abs() on this???
817
#####
pdelta = real( delta ) / ( ray2D( is )%c / co)
818
-
ddelta = -a / ( ck*sb*d ) - a*cco / ssi / (ck*sb*d) + a*cco / (ck*b*sb*d) &
819
#####
-a*co / si / (ck*sb*d*d) * (2* HS%rho * HS%rho *si*co-2*co*si)
820
#####
rddelta = -real( ddelta )
821
#####
sddelta = rddelta / abs( rddelta )
823
-
! next 3 lines have an update by Diana McCammon to allow a sloping bottom
824
-
! I think the formulas are good, but this won't be reliable because it doesn't have the logic
825
-
! that tracks crossing into new segments after the ray displacement.
827
#####
theta_bot = ATAN2( tBdry( 2 ) , tBdry( 1 )) ! bottom angle
828
#####
ray2D( is1 )%x( 1 ) = ray2D( is1 )%x( 1 ) + real( delta ) * COS( theta_bot ) ! range displacement
829
#####
ray2D( is1 )%x( 2 ) = ray2D( is1 )%x( 2 ) + real( delta ) * SIN( theta_bot ) ! depth displacement
830
#####
ray2D( is1 )%tau = ray2D( is1 )%tau + pdelta ! phase change
831
#####
ray2D( is1 )%q = ray2D( is1 )%q + sddelta * rddelta * si * c * ray2D( is )%p ! beam-width change
836
#####
WRITE( PRTFile, * ) 'HS%BC = ', HS%BC
837
1252604*
CALL ERROUT( 'Reflect2D', 'Unknown boundary condition type' )
840
1252604
END SUBROUTINE Reflect2D
842
-
END PROGRAM BELLHOP