1
-
!! Boundary conditions module for altimetry and bathymetry data
4
-
!! Boundary handling for altimetry (top) and bathymetry (bottom) with interpolation capabilities
6
-
! Loads altimetry (top bdry) and bathymetry (bottom bdry) data
16
-
INTEGER, PARAMETER :: ATIFile = 40, BTYFile = 41, Number_to_Echo = 21
17
-
INTEGER :: IsegTop, IsegBot ! indices that point to the current active segment
18
-
INTEGER, PROTECTED :: NATIPts = 2, NBTYPts = 2
19
-
INTEGER :: ii, IOStat, IAllocStat, iSmallStepCtr = 0
21
-
REAL (KIND=8) :: rTopseg( 2 ), rBotseg( 2 ) ! range intervals defining the current active segment
22
-
CHARACTER (LEN=2) :: atiType= 'LS', btyType = 'LS'
24
-
! Halfspace properties
26
-
REAL (KIND=8) :: alphaR, alphaI, betaR, betaI ! compressional and shear wave speeds/attenuations in user units
27
-
COMPLEX (KIND=8) :: cP, cS ! P-wave, S-wave speeds
28
-
REAL (KIND=8) :: rho, Depth ! density, depth
29
-
CHARACTER (LEN=1) :: BC ! Boundary condition type
30
-
CHARACTER (LEN=6) :: Opt
34
-
REAL (KIND=8) :: x( 2 ), t( 2 ), n( 2 ) ! coordinate, tangent, and outward normal for a segment
35
-
REAL (KIND=8) :: Nodet( 2 ), Noden( 2 ) ! tangent and normal at the node, if the curvilinear option is used
36
-
REAL (KIND=8) :: Len, Kappa ! length and curvature of a segment
37
-
REAL (KIND=8) :: Dx, Dxx, Dss ! first, second derivatives wrt depth; s is along tangent
38
-
TYPE( HSInfo2 ) :: HS
41
-
TYPE(BdryPt), ALLOCATABLE :: Top( : ), Bot( : )
45
14*
SUBROUTINE ReadATI( FileRoot, TopATI, DepthT, PRTFile )
46
-
!! Reads in the top altimetry
48
-
INTEGER, INTENT( IN ) :: PRTFile
49
-
CHARACTER (LEN= 1), INTENT( IN ) :: TopATI
50
-
REAL (KIND=8), INTENT( IN ) :: DepthT
51
-
! REAL (KIND=8), ALLOCATABLE :: phi( : ) ! LP: Removed as seems to have been moved to ComputeBdryTangentNormal
52
-
CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
54
#####
SELECT CASE ( TopATI )
56
#####
WRITE( PRTFile, * ) '__________________________________________________________________________'
57
#####
WRITE( PRTFile, * )
58
#####
WRITE( PRTFile, * ) 'Using top-altimetry file'
60
#####
OPEN( UNIT = ATIFile, FILE = TRIM( FileRoot ) // '.ati', STATUS = 'OLD', IOSTAT = IOStat, ACTION = 'READ' )
61
#####
IF ( IOsTAT /= 0 ) THEN
62
#####
WRITE( PRTFile, * ) 'ATIFile = ', TRIM( FileRoot ) // '.ati'
63
#####
CALL ERROUT( 'ReadATI', 'Unable to open altimetry file' )
66
#####
READ( ATIFile, * ) atiType
67
#####
AltiType: SELECT CASE ( atiType( 1 : 1 ) )
69
#####
WRITE( PRTFile, * ) 'Curvilinear Interpolation'
71
#####
WRITE( PRTFile, * ) 'Piecewise linear interpolation'
73
#####
CALL ERROUT( 'ReadATI', 'Unknown option for selecting altimetry interpolation' )
76
#####
READ( ATIFile, * ) NatiPts
77
#####
WRITE( PRTFile, * ) 'Number of altimetry points = ', NatiPts
78
#####
NatiPts = NatiPts + 2 ! we'll be extending the altimetry to infinity to the left and right
80
#####
ALLOCATE( Top( NatiPts ), Stat = IAllocStat )
81
-
! ALLOCATE( phi( NatiPts ), Stat = IAllocStat ) ! LP: Removed as seems to have been moved to ComputeBdryTangentNormal
82
#####
IF ( IAllocStat /= 0 ) &
83
#####
CALL ERROUT( 'BELLHOP:ReadATI', 'Insufficient memory for altimetry data: reduce # ati points' )
85
#####
WRITE( PRTFile, * )
86
#####
WRITE( PRTFile, * ) ' Range (km) Depth (m)'
88
#####
atiPt: DO ii = 2, NatiPts - 1
90
#####
SELECT CASE ( atiType( 2 : 2 ) )
92
#####
READ( ATIFile, * ) Top( ii )%x
93
-
! LP: This condition was previously ii == NatiPts,
94
-
! which will never be satisfied due to the loop bounds
95
#####
IF ( ii < Number_to_Echo .OR. ii == NatiPts - 1 ) THEN ! echo some values
96
#####
WRITE( PRTFile, FMT = "(2G11.3)" ) Top( ii )%x
99
#####
READ( ATIFile, * ) Top( ii )%x, Top( ii )%HS%alphaR, Top( ii )%HS%betaR, Top( ii )%HS%rho, &
100
#####
Top( ii )%HS%alphaI, Top( ii )%HS%betaI
101
-
! LP: Same change as above
102
#####
IF ( ii < Number_to_Echo .OR. ii == NatiPts - 1 ) THEN ! echo some values
103
#####
WRITE( PRTFile, FMT = "(7G11.3)" ) Top( ii )%x, Top( ii )%HS%alphaR, Top( ii )%HS%betaR, Top( ii )%HS%rho, &
104
#####
Top( ii )%HS%alphaI, Top( ii )%HS%betaI
107
#####
CALL ERROUT( 'ReadATI', 'Unknown option for selecting altimetry option' )
110
#####
IF ( Top( ii )%x( 2 ) < DepthT ) THEN
111
#####
CALL ERROUT( 'BELLHOP:ReadATI', 'Altimetry rises above highest point in the sound speed profile' )
115
#####
CLOSE( ATIFile )
117
#####
Top( : )%x( 1 ) = 1000.0 * Top( : )%x( 1 ) ! Convert ranges in km to m
119
-
CASE DEFAULT ! no altimetry given, use SSP depth for flat top
120
14*
ALLOCATE( Top( 2 ), Stat = IAllocStat )
121
14*
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', 'Insufficient memory for altimetry data' )
122
42*
Top( 1 )%x = [ -sqrt( huge( Top( 1 )%x( 1 ) ) ) / 1.0d5, DepthT ]
123
56*
Top( 2 )%x = [ sqrt( huge( Top( 1 )%x( 1 ) ) ) / 1.0d5, DepthT ]
126
14
CALL ComputeBdryTangentNormal( Top, 'Top' )
128
42*
IF ( .NOT. monotonic( Top%x( 1 ), NAtiPts ) ) THEN
129
#####
CALL ERROUT( 'BELLHOP:ReadATI', 'Altimetry ranges are not monotonically increasing' )
132
14
END SUBROUTINE ReadATI
134
-
! **********************************************************************!
136
14*
SUBROUTINE ReadBTY( FileRoot, BotBTY, DepthB, PRTFile )
137
-
!! Reads in the bottom bathymetry
139
-
INTEGER, INTENT( IN ) :: PRTFile
140
-
CHARACTER (LEN= 1), INTENT( IN ) :: BotBTY
141
-
REAL (KIND=8), INTENT( IN ) :: DepthB
142
-
CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
144
4
SELECT CASE ( BotBTY )
146
4
WRITE( PRTFile, * ) '__________________________________________________________________________'
147
4
WRITE( PRTFile, * )
148
4
WRITE( PRTFile, * ) 'Using bottom-bathymetry file'
150
4*
OPEN( UNIT = BTYFile, FILE = TRIM( FileRoot ) // '.bty', STATUS = 'OLD', IOSTAT = IOStat, ACTION = 'READ' )
151
4
IF ( IOsTAT /= 0 ) THEN
152
#####
WRITE( PRTFile, * ) 'BTYFile = ', TRIM( FileRoot ) // '.bty'
153
#####
CALL ERROUT( 'ReadBTY', 'Unable to open bathymetry file' )
156
4
READ( BTYFile, * ) btyType
158
#####
BathyType: SELECT CASE ( btyType( 1 : 1 ) )
160
#####
WRITE( PRTFile, * ) 'Curvilinear Interpolation'
162
4*
WRITE( PRTFile, * ) 'Piecewise linear interpolation'
164
4*
CALL ERROUT( 'ReadBTY', 'Unknown option for selecting bathymetry interpolation' )
165
-
END SELECT BathyType
167
4
READ( BTYFile, * ) NbtyPts
168
4
WRITE( PRTFile, * ) 'Number of bathymetry points = ', NbtyPts
170
4
NbtyPts = NbtyPts + 2 ! we'll be extending the bathymetry to infinity on both sides
171
4*
ALLOCATE( Bot( NbtyPts ), Stat = IAllocStat )
172
4
IF ( IAllocStat /= 0 ) &
173
#####
CALL ERROUT( 'BELLHOP:ReadBTY', 'Insufficient memory for bathymetry data: reduce # bty points' )
175
4
WRITE( PRTFile, * )
176
4
BathyTypeB: SELECT CASE ( btyType( 2 : 2 ) )
178
4
WRITE( PRTFile, * ) 'Short format (bathymetry only)'
179
4*
WRITE( PRTFile, * ) ' Range (km) Depth (m)'
181
#####
WRITE( PRTFile, * ) 'Long format (bathymetry and geoacoustics)'
182
#####
WRITE( PRTFile, "( ' Range (km) Depth (m) alphaR (m/s) betaR rho (g/cm^3) alphaI betaI', / )" )
184
4*
CALL ERROUT( 'ReadBTY', 'Unknown option for selecting bathymetry interpolation' )
185
-
END SELECT BathyTypeB
187
44*
btyPt: DO ii = 2, NbtyPts - 1
189
40
SELECT CASE ( btyType( 2 : 2 ) )
190
-
CASE ( 'S', '' ) ! short format
191
120*
READ( BTYFile, * ) Bot( ii )%x
192
-
! LP: This condition was previously ii == NbtyPts,
193
-
! which will never be satisfied due to the loop bounds
194
70*
IF ( ii < Number_to_Echo .OR. ii == NbtyPts - 1 ) THEN ! echo some values
195
90*
WRITE( PRTFile, FMT = "(2G11.3)" ) Bot( ii )%x
197
-
CASE ( 'L' ) ! long format
198
#####
READ( BTYFile, * ) Bot( ii )%x, Bot( ii )%HS%alphaR, Bot( ii )%HS%betaR, Bot( ii )%HS%rho, &
199
#####
Bot( ii )%HS%alphaI, Bot( ii )%HS%betaI
200
-
! LP: Same change as above
201
#####
IF ( ii < Number_to_Echo .OR. ii == NbtyPts - 1 ) THEN ! echo some values
202
#####
WRITE( PRTFile, FMT="( F10.2, F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )" ) &
203
#####
Bot( ii )%x, Bot( ii )%HS%alphaR, Bot( ii )%HS%betaR, Bot( ii )%HS%rho, &
204
#####
Bot( ii )%HS%alphaI, Bot( ii )%HS%betaI
207
40*
CALL ERROUT( 'ReadBTY', 'Unknown option for selecting bathymetry option' )
210
44*
IF ( Bot( ii )%x( 2 ) > DepthB ) THEN
211
#####
CALL ERROUT( 'BELLHOP:ReadBTY', 'Bathymetry drops below lowest point in the sound speed profile' )
218
52*
Bot( : )%x( 1 ) = 1000.0 * Bot( : )%x( 1 ) ! Convert ranges in km to m
220
-
CASE DEFAULT ! no bathymetry given, use SSP depth for flat bottom
221
10*
ALLOCATE( Bot( 2 ), Stat = IAllocStat )
222
10*
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', 'Insufficient memory for bathymetry data' )
223
30*
Bot( 1 )%x = [ -sqrt( huge( Bot( 1 )%x( 1 ) ) ) / 1.0d5, DepthB ]
224
44*
Bot( 2 )%x = [ sqrt( huge( Bot( 1 )%x( 1 ) ) ) / 1.0d5, DepthB ]
227
14
CALL ComputeBdryTangentNormal( Bot, 'Bot' )
229
82*
IF ( .NOT. monotonic( Bot%x( 1 ), NBtyPts ) ) THEN
230
#####
CALL ERROUT( 'BELLHOP:ReadBTY', 'Bathymetry ranges are not monotonically increasing' )
233
14
END SUBROUTINE ReadBTY
235
-
! **********************************************************************!
237
28*
SUBROUTINE ComputeBdryTangentNormal( Bdry, BotTop )
239
-
! Does some pre-processing on the boundary points to pre-compute segment
241
-
! tangents (%t, %nodet),
242
-
! normals (%n, %noden), and
243
-
! curvatures (%kappa)
245
-
! The boundary is also extended with a constant depth to infinity to cover cases where the ray
246
-
! exits the domain defined by the user
248
-
INTEGER, SAVE :: NPts = 0
249
28
REAL (KIND=8), ALLOCATABLE :: phi( : )
250
-
REAL (KIND=8) :: sss
251
-
TYPE(BdryPt), INTENT( INOUT ) :: Bdry( : )
252
-
CHARACTER (LEN=3), INTENT( IN ) :: BotTop ! Flag indicating bottom or top reflection
253
-
CHARACTER (LEN=2), SAVE :: CurvilinearFlag = '-'
255
14
SELECT CASE ( BotTop )
258
14
CurvilinearFlag = btyType
261
28
CurvilinearFlag = atiType
264
-
! extend the bathymetry to +/- infinity in a piecewise constant fashion
266
28*
Bdry( 1 )%x( 1 ) = -sqrt( huge( Bdry( 1 )%x( 1 ) ) ) / 1.0d5
267
28*
Bdry( 1 )%x( 2 ) = Bdry( 2 )%x( 2 )
268
28*
Bdry( 1 )%HS = Bdry( 2 )%HS
269
28*
Bdry( NPts )%x( 1 ) = +sqrt( huge( Bdry( 1 )%x( 1 ) ) ) / 1.0d5
270
28*
Bdry( NPts )%x( 2 ) = Bdry( NPts - 1 )%x( 2 )
271
28*
Bdry( NPts )%HS = Bdry( NPts - 1 )%HS
273
-
! compute tangent and outward-pointing normal to each bottom segment
274
-
! tBdry( 1, : ) = xBdry( 1, 2:NPts ) - xBdry( 1, 1:NPts - 1 )
275
-
! tBdry( 2, : ) = xBdry( 2, 2:NPts ) - xBdry( 2, 1:NPts - 1 )
276
-
! above caused compiler problems
278
96*
BoundaryPt: DO ii = 1, NPts - 1
279
204*
Bdry( ii )%t = Bdry( ii + 1 )%x - Bdry( ii )%x
280
68*
Bdry( ii )%Dx = Bdry( ii )%t( 2 ) / Bdry( ii )%t( 1 ) ! first derivative
281
-
! write( *, * ) 'Dx, t', Bdry( ii )%Dx, Bdry( ii )%x, 1 / ( Bdry( ii )%x( 2 ) / 500 )
283
-
! normalize the tangent vector
284
204*
Bdry( ii )%Len = NORM2( Bdry( ii )%t )
285
204*
Bdry( ii )%t = Bdry( ii )%t / Bdry( ii )%Len
287
28
SELECT CASE ( BotTop )
289
54*
Bdry( ii )%n( 1 ) = -Bdry( ii )%t( 2 )
290
54*
Bdry( ii )%n( 2 ) = +Bdry( ii )%t( 1 )
292
14*
Bdry( ii )%n( 1 ) = +Bdry( ii )%t( 2 )
293
82*
Bdry( ii )%n( 2 ) = -Bdry( ii )%t( 1 )
298
28
IF ( CurvilinearFlag( 1 : 1 ) == 'C' ) THEN ! curvilinear option: compute tangent and normal at node by averaging normals on adjacent segments
299
-
! averaging two centered differences is equivalent to forming a single centered difference of two steps ...
300
#####
DO ii = 2, NPts - 1
301
#####
sss = Bdry( ii - 1 )%Len / ( Bdry( ii - 1 )%Len + Bdry( ii )%Len )
302
#####
sss = 0.5 ! LP: BUG? Line above is overwritten.
303
#####
Bdry( ii )%Nodet = ( 1.0 - sss ) * Bdry( ii - 1 )%t + sss * Bdry( ii )%t
306
#####
Bdry( 1 )%Nodet = [ 1.0, 0.0 ] ! tangent left-end node
307
#####
Bdry( NPts )%Nodet = [ 1.0, 0.0 ] ! tangent right-end node
309
#####
SELECT CASE ( BotTop )
311
#####
Bdry( : )%Noden( 1 ) = -Bdry( : )%Nodet( 2 )
312
#####
Bdry( : )%Noden( 2 ) = +Bdry( : )%Nodet( 1 )
314
#####
Bdry( : )%Noden( 1 ) = +Bdry( : )%Nodet( 2 )
315
#####
Bdry( : )%Noden( 2 ) = -Bdry( : )%Nodet( 1 )
318
-
! compute curvature in each segment
319
-
! LP: This allocation is not necessary, could just have two variables for
320
-
! current and next phi. Operating on the whole array can trigger compiler
321
-
! SIMD parallelism (AVX-512 etc.), but this is unlikely to happen for
322
-
! atan2, and this is in one-time setup code anyway.
323
#####
ALLOCATE( phi( NPts ), Stat = IAllocStat )
324
#####
phi = atan2( Bdry( : )%Nodet( 2 ), Bdry( : )%Nodet( 1 ) ) ! this is the angle at each node
326
#####
DO ii = 1, NPts - 1
327
#####
Bdry( ii )%kappa = ( phi( ii + 1 ) - phi( ii ) ) / Bdry( ii )%Len ! this is curvature = dphi/ds
328
#####
Bdry( ii )%Dxx = ( Bdry( ii + 1 )%Dx - Bdry( ii )%Dx ) / & ! second derivative
329
#####
( Bdry( ii + 1 )%x( 1 ) - Bdry( ii )%x( 1 ) )
330
#####
Bdry( ii )%Dss = Bdry( ii )%Dxx * Bdry( ii )%t( 1 ) ** 3 ! derivative in direction of tangent
331
-
!write( *, * ) 'kappa, Dss, Dxx', Bdry( ii )%kappa, Bdry( ii )%Dss, Bdry( ii )%Dxx, &
332
-
! 1 / ( ( 8 / 1000 ** 2 ) * ABS( Bdry( ii )%x( 2 ) ) ** 3 ), Bdry( ii )%x( 2 )
333
-
! -1 / ( 4 * ( Bdry( ii )%x( 2 ) ) ** 3 / 1000000 ), Bdry( ii )%x( 2 )
335
#####
Bdry( ii )%kappa = Bdry( ii )%Dss !over-ride kappa !!!!!
338
#####
DEALLOCATE( phi ) ! LP: was missing deallocation
343
28*
END SUBROUTINE ComputeBdryTangentNormal
345
-
! **********************************************************************!
347
55950
SUBROUTINE GetTopSeg( r, t )
349
-
! Get the Top segment info (index and range interval) for range, r
350
-
! LP: t: range component of ray tangent. Endpoints of segments are handled
351
-
! so that if the ray moves slightly along its current direction, it will
352
-
! remain in the same segment.
354
-
INTEGER, PARAMETER :: PRTFile = 6
355
-
INTEGER :: IsegTopT( 1 )
356
-
REAL (KIND=8), INTENT( IN ) :: r, t
358
55950
IF ( t > 0.0 ) THEN
359
279750*
IsegTopT = MAXLOC( Top( : )%x( 1 ), Top( : )%x( 1 ) <= r )
361
#####
IsegTopT = MAXLOC( Top( : )%x( 1 ), Top( : )%x( 1 ) < r )
364
55950
IF ( IsegTopT( 1 ) > 0 .AND. IsegTopT( 1 ) < NatiPts ) THEN ! IsegTop MUST LIE IN [ 1, NatiPts-1 ]
365
55950
IsegTop = IsegTopT( 1 )
366
167850*
rTopSeg = [ Top( IsegTop )%x( 1 ), Top( IsegTop + 1 )%x( 1 ) ] ! segment limits in range
368
#####
WRITE( PRTFile, * ) 'r = ', r
369
#####
WRITE( PRTFile, * ) 'rLeft = ', Top( 1 )%x( 1 )
370
#####
WRITE( PRTFile, * ) 'rRight = ', Top( NatiPts )%x( 1 )
371
#####
CALL ERROUT( 'GetTopSeg', 'Top altimetry undefined above the ray' )
374
55950
END SUBROUTINE GetTopSeg
376
-
! **********************************************************************!
378
80303
SUBROUTINE GetBotSeg( r, t )
380
-
! Get the Bottom segment info (index and range interval) for range, r
381
-
! LP: t: range component of ray tangent. Endpoints of segments are handled
382
-
! so that if the ray moves slightly along its current direction, it will
383
-
! remain in the same segment.
385
-
INTEGER, PARAMETER :: PRTFile = 6
386
-
INTEGER :: IsegBotT( 1 )
387
-
REAL (KIND=8), INTENT( IN ) :: r, t
389
80303
IF ( t > 0.0 ) THEN
390
1104045*
IsegBotT = MAXLOC( Bot( : )%x( 1 ), Bot( : )%x( 1 ) <= r )
392
382482*
IsegBotT = MAXLOC( Bot( : )%x( 1 ), Bot( : )%x( 1 ) < r )
395
80303
IF ( IsegBotT( 1 ) > 0 .AND. IsegBotT( 1 ) < NbtyPts ) THEN ! IsegBot MUST LIE IN [ 1, NbtyPts-1 ]
396
80303
IsegBot = IsegBotT( 1 )
397
240909*
rBotSeg = [ Bot( IsegBot )%x( 1 ), Bot( IsegBot + 1 )%x( 1 ) ] ! segment limits in range
399
#####
WRITE( PRTFile, * ) 'r = ', r
400
#####
WRITE( PRTFile, * ) 'rLeft = ', Bot( 1 )%x( 1 )
401
#####
WRITE( PRTFile, * ) 'rRight = ', Bot( NbtyPts )%x( 1 )
402
#####
CALL ERROUT( 'GetBotSeg', 'Bottom bathymetry undefined below the source' )
405
80303
END SUBROUTINE GetBotSeg
407
#####
END MODULE bdrymod
407
#####
END MODULE bdrymod
407
#####
END MODULE bdrymod