Coverage Report: bdryMod.f90

Generated from GCOV analysis of Fortran source code

78.9%
Lines Executed
175 total lines
56.7%
Branches Executed
736 total branches
100.0%
Calls Executed
93 total calls
0
-
Source:bdryMod.f90
0
-
Graph:bdryMod.gcno
0
-
Data:bdryMod.gcda
0
-
Runs:73
1
-
!! Boundary conditions module for altimetry and bathymetry data
2
-
3
-
MODULE bdrymod
4
-
!! Boundary handling for altimetry (top) and bathymetry (bottom) with interpolation capabilities
5
-
6
-
! Loads altimetry (top bdry) and bathymetry (bottom bdry) data
7
-
8
-
!USE norms
9
-
USE monotonicMod
10
-
USE FatalError
11
-
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT
12
-
13
-
IMPLICIT NONE
14
-
PUBLIC
15
-
SAVE
16
-
17
-
INTEGER, PARAMETER :: ATIFile = 40, BTYFile = 41, Number_to_Echo = 21
18
-
INTEGER :: IsegTop, IsegBot ! indices that point to the current active segment
19
-
INTEGER, PROTECTED :: NATIPts = 2, NBTYPts = 2
20
-
INTEGER :: ii, IOStat, IAllocStat, iSmallStepCtr = 0
21
-
22
-
REAL (KIND=8) :: rTopseg( 2 ), rBotseg( 2 ) ! range intervals defining the current active segment
23
-
CHARACTER (LEN=2) :: atiType= 'LS', btyType = 'LS'
24
-
25
-
! Halfspace properties
26
-
TYPE HSInfo2
27
-
REAL (KIND=8) :: alphaR, alphaI, betaR, betaI ! compressional and shear wave speeds/attenuations in user units
28
-
COMPLEX (KIND=8) :: cP, cS ! P-wave, S-wave speeds
29
-
REAL (KIND=8) :: rho, Depth ! density, depth
30
-
CHARACTER (LEN=1) :: BC ! Boundary condition type
31
-
CHARACTER (LEN=6) :: Opt
32
-
END TYPE HSInfo2
33
-
34
-
TYPE BdryPt
35
-
REAL (KIND=8) :: x( 2 ), t( 2 ), n( 2 ) ! coordinate, tangent, and outward normal for a segment
36
-
REAL (KIND=8) :: Nodet( 2 ), Noden( 2 ) ! tangent and normal at the node, if the curvilinear option is used
37
-
REAL (KIND=8) :: Len, Kappa ! length and curvature of a segment
38
-
REAL (KIND=8) :: Dx, Dxx, Dss ! first, second derivatives wrt depth; s is along tangent
39
-
TYPE( HSInfo2 ) :: HS
40
-
END TYPE BdryPt
41
-
42
-
TYPE(BdryPt), ALLOCATABLE :: Top( : ), Bot( : )
43
-
44
-
CONTAINS
45
-
46
71*
SUBROUTINE ReadATI( FileRoot, TopATI, DepthT, PRTFile )
47
-
!! Reads in the top altimetry
48
-
49
-
INTEGER, INTENT( IN ) :: PRTFile
50
-
CHARACTER (LEN= 1), INTENT( IN ) :: TopATI
51
-
REAL (KIND=8), INTENT( IN ) :: DepthT
52
-
! REAL (KIND=8), ALLOCATABLE :: phi( : ) ! LP: Removed as seems to have been moved to ComputeBdryTangentNormal
53
-
CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
54
-
55
6
SELECT CASE ( TopATI )
56
-
CASE ( '~', '*' )
57
6
WRITE( PRTFile, * ) '__________________________________________________________________________'
58
6
WRITE( PRTFile, * )
59
6
WRITE( PRTFile, * ) 'Using top-altimetry file'
60
-
61
6*
OPEN( UNIT = ATIFile, FILE = TRIM( FileRoot ) // '.ati', STATUS = 'OLD', IOSTAT = IOStat, ACTION = 'READ' )
62
6
IF ( IOsTAT /= 0 ) THEN
63
#####
WRITE( PRTFile, * ) 'ATIFile = ', TRIM( FileRoot ) // '.ati'
64
#####
CALL ERROUT( 'ReadATI', 'Unable to open altimetry file' )
65
-
END IF
66
-
67
6
READ( ATIFile, * ) atiType
68
1
AltiType: SELECT CASE ( atiType( 1 : 1 ) )
69
-
CASE ( 'C' )
70
1
WRITE( PRTFile, * ) 'Curvilinear Interpolation'
71
-
CASE ( 'L' )
72
5*
WRITE( PRTFile, * ) 'Piecewise linear interpolation'
73
-
CASE DEFAULT
74
6*
CALL ERROUT( 'ReadATI', 'Unknown option for selecting altimetry interpolation' )
75
-
END SELECT AltiType
76
-
77
6
READ( ATIFile, * ) NatiPts
78
6
WRITE( PRTFile, * ) 'Number of altimetry points = ', NatiPts
79
6
NatiPts = NatiPts + 2 ! we'll be extending the altimetry to infinity to the left and right
80
-
81
6*
ALLOCATE( Top( NatiPts ), Stat = IAllocStat )
82
-
! ALLOCATE( phi( NatiPts ), Stat = IAllocStat ) ! LP: Removed as seems to have been moved to ComputeBdryTangentNormal
83
6
IF ( IAllocStat /= 0 ) THEN
84
#####
CALL ERROUT( 'BELLHOP:ReadATI', 'Insufficient memory for altimetry data: reduce # ati points' )
85
-
END IF
86
-
87
6
WRITE( PRTFile, * )
88
6
WRITE( PRTFile, * ) ' Range (km) Depth (m)'
89
-
90
6011*
atiPt: DO ii = 2, NatiPts - 1
91
-
92
6005
SELECT CASE ( atiType( 2 : 2 ) )
93
-
CASE ( 'S', '' )
94
18015*
READ( ATIFile, * ) Top( ii )%x
95
-
! LP: This condition was previously ii == NatiPts,
96
-
! which will never be satisfied due to the loop bounds
97
6125*
IF ( ii < Number_to_Echo .OR. ii == NatiPts - 1 ) THEN ! echo some values
98
360*
WRITE( PRTFile, FMT = "(2G11.3)" ) Top( ii )%x
99
-
END IF
100
-
CASE ( 'L' )
101
#####
READ( ATIFile, * ) Top( ii )%x, Top( ii )%HS%alphaR, Top( ii )%HS%betaR, Top( ii )%HS%rho, &
102
#####
Top( ii )%HS%alphaI, Top( ii )%HS%betaI
103
-
! LP: Same change as above
104
#####
IF ( ii < Number_to_Echo .OR. ii == NatiPts - 1 ) THEN ! echo some values
105
#####
WRITE( PRTFile, FMT = "(7G11.3)" ) Top( ii )%x, Top( ii )%HS%alphaR, Top( ii )%HS%betaR, Top( ii )%HS%rho, &
106
#####
Top( ii )%HS%alphaI, Top( ii )%HS%betaI
107
-
END IF
108
-
CASE DEFAULT
109
6005*
CALL ERROUT( 'ReadATI', 'Unknown option for selecting altimetry option' )
110
-
END SELECT
111
-
112
6011*
IF ( Top( ii )%x( 2 ) < DepthT ) THEN
113
#####
CALL ERROUT( 'BELLHOP:ReadATI', 'Altimetry rises above highest point in the sound speed profile' )
114
-
END IF
115
-
END DO atiPt
116
-
117
6
CLOSE( ATIFile )
118
-
119
6023*
Top( : )%x( 1 ) = 1000.0 * Top( : )%x( 1 ) ! Convert ranges in km to m
120
-
121
-
CASE DEFAULT ! no altimetry given, use SSP depth for flat top
122
65*
ALLOCATE( Top( 2 ), Stat = IAllocStat )
123
65*
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', 'Insufficient memory for altimetry data' )
124
195*
Top( 1 )%x = [ -sqrt( huge( Top( 1 )%x( 1 ) ) ) / 1.0d5, DepthT ]
125
266*
Top( 2 )%x = [ sqrt( huge( Top( 1 )%x( 1 ) ) ) / 1.0d5, DepthT ]
126
-
END SELECT
127
-
128
71
CALL ComputeBdryTangentNormal( Top, 'Top' )
129
-
130
6218*
IF ( .NOT. monotonic( Top%x( 1 ), NAtiPts ) ) THEN
131
#####
CALL ERROUT( 'BELLHOP:ReadATI', 'Altimetry ranges are not monotonically increasing' )
132
-
END IF
133
-
134
71
END SUBROUTINE ReadATI
135
-
136
-
! **********************************************************************!
137
-
138
71*
SUBROUTINE ReadBTY( FileRoot, BotBTY, DepthB, PRTFile )
139
-
!! Reads in the bottom bathymetry
140
-
141
-
INTEGER, INTENT( IN ) :: PRTFile
142
-
CHARACTER (LEN= 1), INTENT( IN ) :: BotBTY
143
-
REAL (KIND=8), INTENT( IN ) :: DepthB
144
-
CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
145
-
146
11
SELECT CASE ( BotBTY )
147
-
CASE ( '~', '*' )
148
11
WRITE( PRTFile, * ) '__________________________________________________________________________'
149
11
WRITE( PRTFile, * )
150
11
WRITE( PRTFile, * ) 'Using bottom-bathymetry file'
151
-
152
11*
OPEN( UNIT = BTYFile, FILE = TRIM( FileRoot ) // '.bty', STATUS = 'OLD', IOSTAT = IOStat, ACTION = 'READ' )
153
11
IF ( IOsTAT /= 0 ) THEN
154
#####
WRITE( PRTFile, * ) 'BTYFile = ', TRIM( FileRoot ) // '.bty'
155
#####
CALL ERROUT( 'ReadBTY', 'Unable to open bathymetry file' )
156
-
END IF
157
-
158
11
READ( BTYFile, * ) btyType
159
-
160
1
BathyType: SELECT CASE ( btyType( 1 : 1 ) )
161
-
CASE ( 'C' )
162
1
WRITE( PRTFile, * ) 'Curvilinear Interpolation'
163
-
CASE ( 'L' )
164
10*
WRITE( PRTFile, * ) 'Piecewise linear interpolation'
165
-
CASE DEFAULT
166
11*
CALL ERROUT( 'ReadBTY', 'Unknown option for selecting bathymetry interpolation' )
167
-
END SELECT BathyType
168
-
169
11
READ( BTYFile, * ) NbtyPts
170
11
WRITE( PRTFile, * ) 'Number of bathymetry points = ', NbtyPts
171
-
172
11
NbtyPts = NbtyPts + 2 ! we'll be extending the bathymetry to infinity on both sides
173
11*
ALLOCATE( Bot( NbtyPts ), Stat = IAllocStat )
174
11
IF ( IAllocStat /= 0 ) THEN
175
#####
CALL ERROUT( 'BELLHOP:ReadBTY', 'Insufficient memory for bathymetry data: reduce # bty points' )
176
-
END IF
177
-
178
11
WRITE( PRTFile, * )
179
11
BathyTypeB: SELECT CASE ( btyType( 2 : 2 ) )
180
-
CASE ( 'S', '' )
181
11
WRITE( PRTFile, * ) 'Short format (bathymetry only)'
182
11*
WRITE( PRTFile, * ) ' Range (km) Depth (m)'
183
-
CASE ( 'L' )
184
#####
WRITE( PRTFile, * ) 'Long format (bathymetry and geoacoustics)'
185
#####
WRITE( PRTFile, "( ' Range (km) Depth (m) alphaR (m/s) betaR rho (g/cm^3) alphaI betaI', / )" )
186
-
CASE DEFAULT
187
11*
CALL ERROUT( 'ReadBTY', 'Unknown option for selecting bathymetry interpolation' )
188
-
END SELECT BathyTypeB
189
-
190
1078*
btyPt: DO ii = 2, NbtyPts - 1
191
-
192
1067
SELECT CASE ( btyType( 2 : 2 ) )
193
-
CASE ( 'S', '' ) ! short format
194
3201*
READ( BTYFile, * ) Bot( ii )%x
195
-
! LP: This condition was previously ii == NbtyPts,
196
-
! which will never be satisfied due to the loop bounds
197
1144*
IF ( ii < Number_to_Echo .OR. ii == NbtyPts - 1 ) THEN ! echo some values
198
231*
WRITE( PRTFile, FMT = "(2G11.3)" ) Bot( ii )%x
199
-
END IF
200
-
CASE ( 'L' ) ! long format
201
#####
READ( BTYFile, * ) Bot( ii )%x, Bot( ii )%HS%alphaR, Bot( ii )%HS%betaR, Bot( ii )%HS%rho, &
202
#####
Bot( ii )%HS%alphaI, Bot( ii )%HS%betaI
203
-
! LP: Same change as above
204
#####
IF ( ii < Number_to_Echo .OR. ii == NbtyPts - 1 ) THEN ! echo some values
205
#####
WRITE( PRTFile, FMT="( F10.2, F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )" ) &
206
#####
Bot( ii )%x, Bot( ii )%HS%alphaR, Bot( ii )%HS%betaR, Bot( ii )%HS%rho, &
207
#####
Bot( ii )%HS%alphaI, Bot( ii )%HS%betaI
208
-
END IF
209
-
CASE DEFAULT
210
1067*
CALL ERROUT( 'ReadBTY', 'Unknown option for selecting bathymetry option' )
211
-
END SELECT
212
-
213
1078*
IF ( Bot( ii )%x( 2 ) > DepthB ) THEN
214
#####
CALL ERROUT( 'BELLHOP:ReadBTY', 'Bathymetry drops below lowest point in the sound speed profile' )
215
-
END IF
216
-
217
-
END DO btypt
218
-
219
11
CLOSE( BTYFile )
220
-
221
1100*
Bot( : )%x( 1 ) = 1000.0 * Bot( : )%x( 1 ) ! Convert ranges in km to m
222
-
223
-
CASE DEFAULT ! no bathymetry given, use SSP depth for flat bottom
224
60*
ALLOCATE( Bot( 2 ), Stat = IAllocStat )
225
60*
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', 'Insufficient memory for bathymetry data' )
226
180*
Bot( 1 )%x = [ -sqrt( huge( Bot( 1 )%x( 1 ) ) ) / 1.0d5, DepthB ]
227
251*
Bot( 2 )%x = [ sqrt( huge( Bot( 1 )%x( 1 ) ) ) / 1.0d5, DepthB ]
228
-
END SELECT
229
-
230
71
CALL ComputeBdryTangentNormal( Bot, 'Bot' )
231
-
232
1280*
IF ( .NOT. monotonic( Bot%x( 1 ), NBtyPts ) ) THEN
233
#####
CALL ERROUT( 'BELLHOP:ReadBTY', 'Bathymetry ranges are not monotonically increasing' )
234
-
END IF
235
-
236
71
END SUBROUTINE ReadBTY
237
-
238
-
! **********************************************************************!
239
-
240
142*
SUBROUTINE ComputeBdryTangentNormal( Bdry, BotTop )
241
-
242
-
! Does some pre-processing on the boundary points to pre-compute segment
243
-
! lengths (%Len),
244
-
! tangents (%t, %nodet),
245
-
! normals (%n, %noden), and
246
-
! curvatures (%kappa)
247
-
!
248
-
! The boundary is also extended with a constant depth to infinity to cover cases where the ray
249
-
! exits the domain defined by the user
250
-
251
-
INTEGER, SAVE :: NPts = 0
252
142
REAL (KIND=8), ALLOCATABLE :: phi( : )
253
-
REAL (KIND=8) :: sss
254
-
TYPE(BdryPt), INTENT( INOUT ) :: Bdry( : )
255
-
CHARACTER (LEN=3), INTENT( IN ) :: BotTop ! Flag indicating bottom or top reflection
256
-
CHARACTER (LEN=2), SAVE :: CurvilinearFlag = '-'
257
-
258
71
SELECT CASE ( BotTop )
259
-
CASE ( 'Bot' )
260
71
NPts = NbtyPts
261
71
CurvilinearFlag = btyType
262
-
CASE ( 'Top' )
263
71
NPts = NatiPts
264
71*
CurvilinearFlag = atiType
265
-
CASE DEFAULT
266
#####
WRITE( ERROR_UNIT, * ) 'ComputeBdryTangentNormal: Unknown BotTop flag: ', BotTop
267
142*
CALL ERROUT( 'ComputeBdryTangentNormal', 'Unknown boundary flag' )
268
-
END SELECT
269
-
270
-
! extend the bathymetry to +/- infinity in a piecewise constant fashion
271
-
272
142*
Bdry( 1 )%x( 1 ) = -sqrt( huge( Bdry( 1 )%x( 1 ) ) ) / 1.0d5
273
142*
Bdry( 1 )%x( 2 ) = Bdry( 2 )%x( 2 )
274
142*
Bdry( 1 )%HS = Bdry( 2 )%HS
275
142*
Bdry( NPts )%x( 1 ) = +sqrt( huge( Bdry( 1 )%x( 1 ) ) ) / 1.0d5
276
142*
Bdry( NPts )%x( 2 ) = Bdry( NPts - 1 )%x( 2 )
277
142*
Bdry( NPts )%HS = Bdry( NPts - 1 )%HS
278
-
279
-
! compute tangent and outward-pointing normal to each bottom segment
280
-
! tBdry( 1, : ) = xBdry( 1, 2:NPts ) - xBdry( 1, 1:NPts - 1 )
281
-
! tBdry( 2, : ) = xBdry( 2, 2:NPts ) - xBdry( 2, 1:NPts - 1 )
282
-
! above caused compiler problems
283
-
284
7356*
BoundaryPt: DO ii = 1, NPts - 1
285
21642*
Bdry( ii )%t = Bdry( ii + 1 )%x - Bdry( ii )%x
286
7214*
Bdry( ii )%Dx = Bdry( ii )%t( 2 ) / Bdry( ii )%t( 1 ) ! first derivative
287
-
! write( *, * ) 'Dx, t', Bdry( ii )%Dx, Bdry( ii )%x, 1 / ( Bdry( ii )%x( 2 ) / 500 )
288
-
289
-
! normalize the tangent vector
290
21642*
Bdry( ii )%Len = NORM2( Bdry( ii )%t )
291
21642*
Bdry( ii )%t = Bdry( ii )%t / Bdry( ii )%Len
292
-
293
142
SELECT CASE ( BotTop )
294
-
CASE ( 'Bot' )
295
1138*
Bdry( ii )%n( 1 ) = -Bdry( ii )%t( 2 )
296
1138*
Bdry( ii )%n( 2 ) = +Bdry( ii )%t( 1 )
297
-
CASE ( 'Top' )
298
6076*
Bdry( ii )%n( 1 ) = +Bdry( ii )%t( 2 )
299
6076*
Bdry( ii )%n( 2 ) = -Bdry( ii )%t( 1 )
300
-
CASE DEFAULT
301
#####
WRITE( ERROR_UNIT, * ) 'ComputeBdryTangentNormal: Unknown BotTop flag in normal calc: ', BotTop
302
7214*
CALL ERROUT( 'ComputeBdryTangentNormal', 'Unknown boundary flag' )
303
-
END SELECT
304
-
305
-
END DO BoundaryPt
306
-
307
142
IF ( CurvilinearFlag( 1 : 1 ) == 'C' ) THEN
308
-
! curvilinear option: compute tangent and normal at node by averaging normals on adjacent segments
309
-
! averaging two centered differences is equivalent to forming a single centered difference of two steps ...
310
2002*
DO ii = 2, NPts - 1
311
2000*
sss = Bdry( ii - 1 )%Len / ( Bdry( ii - 1 )%Len + Bdry( ii )%Len )
312
2000
sss = 0.5 ! LP: BUG? Line above is overwritten.
313
6002*
Bdry( ii )%Nodet = ( 1.0 - sss ) * Bdry( ii - 1 )%t + sss * Bdry( ii )%t
314
-
END DO
315
-
316
6*
Bdry( 1 )%Nodet = [ 1.0, 0.0 ] ! tangent left-end node
317
6*
Bdry( NPts )%Nodet = [ 1.0, 0.0 ] ! tangent right-end node
318
-
319
1
SELECT CASE ( BotTop )
320
-
CASE ( 'Bot' )
321
1003*
Bdry( : )%Noden( 1 ) = -Bdry( : )%Nodet( 2 )
322
1003*
Bdry( : )%Noden( 2 ) = +Bdry( : )%Nodet( 1 )
323
-
CASE ( 'Top' )
324
1003*
Bdry( : )%Noden( 1 ) = +Bdry( : )%Nodet( 2 )
325
1003*
Bdry( : )%Noden( 2 ) = -Bdry( : )%Nodet( 1 )
326
-
CASE DEFAULT
327
#####
WRITE( ERROR_UNIT, * ) 'ComputeBdryTangentNormal: Unknown BotTop flag in node normal calc: ', BotTop
328
2*
CALL ERROUT( 'ComputeBdryTangentNormal', 'Unknown boundary flag' )
329
-
END SELECT
330
-
331
-
! compute curvature in each segment
332
-
! LP: This allocation is not necessary, could just have two variables for
333
-
! current and next phi. Operating on the whole array can trigger compiler
334
-
! SIMD parallelism (AVX-512 etc.), but this is unlikely to happen for
335
-
! atan2, and this is in one-time setup code anyway.
336
2*
ALLOCATE( phi( NPts ), Stat = IAllocStat )
337
2006*
phi = atan2( Bdry( : )%Nodet( 2 ), Bdry( : )%Nodet( 1 ) ) ! this is the angle at each node
338
-
339
2004*
DO ii = 1, NPts - 1
340
2002*
Bdry( ii )%kappa = ( phi( ii + 1 ) - phi( ii ) ) / Bdry( ii )%Len ! this is curvature = dphi/ds
341
6006*
Bdry( ii )%Dxx = ( Bdry( ii + 1 )%Dx - Bdry( ii )%Dx ) / & ! second derivative
342
6006*
( Bdry( ii + 1 )%x( 1 ) - Bdry( ii )%x( 1 ) )
343
2002*
Bdry( ii )%Dss = Bdry( ii )%Dxx * Bdry( ii )%t( 1 ) ** 3 ! derivative in direction of tangent
344
-
!write( *, * ) 'kappa, Dss, Dxx', Bdry( ii )%kappa, Bdry( ii )%Dss, Bdry( ii )%Dxx, &
345
-
! 1 / ( ( 8 / 1000 ** 2 ) * ABS( Bdry( ii )%x( 2 ) ) ** 3 ), Bdry( ii )%x( 2 )
346
-
! -1 / ( 4 * ( Bdry( ii )%x( 2 ) ) ** 3 / 1000000 ), Bdry( ii )%x( 2 )
347
-
348
2004*
Bdry( ii )%kappa = Bdry( ii )%Dss !over-ride kappa !!!!!
349
-
END DO
350
-
351
2*
DEALLOCATE( phi ) ! LP: was missing deallocation
352
-
ELSE
353
5492*
Bdry%kappa = 0
354
-
END IF
355
-
356
142*
END SUBROUTINE ComputeBdryTangentNormal
357
-
358
-
! **********************************************************************!
359
-
360
1220138
SUBROUTINE GetTopSeg( r, t )
361
-
362
-
! Get the Top segment info (index and range interval) for range, r
363
-
! LP: t: range component of ray tangent. Endpoints of segments are handled
364
-
! so that if the ray moves slightly along its current direction, it will
365
-
! remain in the same segment.
366
-
367
-
INTEGER, PARAMETER :: PRTFile = 6
368
-
INTEGER :: IsegTopT( 1 )
369
-
REAL (KIND=8), INTENT( IN ) :: r, t
370
-
371
1220138
IF ( t > 0.0 ) THEN
372
1918284108*
IsegTopT = MAXLOC( Top( : )%x( 1 ), Top( : )%x( 1 ) <= r )
373
-
ELSE
374
136918*
IsegTopT = MAXLOC( Top( : )%x( 1 ), Top( : )%x( 1 ) < r )
375
-
END IF
376
-
377
1220138
IF ( IsegTopT( 1 ) > 0 .AND. IsegTopT( 1 ) < NatiPts ) THEN ! IsegTop MUST LIE IN [ 1, NatiPts-1 ]
378
1220138
IsegTop = IsegTopT( 1 )
379
3660414*
rTopSeg = [ Top( IsegTop )%x( 1 ), Top( IsegTop + 1 )%x( 1 ) ] ! segment limits in range
380
-
ELSE
381
#####
WRITE( PRTFile, * ) 'r = ', r
382
#####
WRITE( PRTFile, * ) 'rLeft = ', Top( 1 )%x( 1 )
383
#####
WRITE( PRTFile, * ) 'rRight = ', Top( NatiPts )%x( 1 )
384
#####
CALL ERROUT( 'GetTopSeg', 'Top altimetry undefined above the ray' )
385
-
END IF
386
-
387
1220138
END SUBROUTINE GetTopSeg
388
-
389
-
! **********************************************************************!
390
-
391
364871
SUBROUTINE GetBotSeg( r, t )
392
-
393
-
! Get the Bottom segment info (index and range interval) for range, r
394
-
! LP: t: range component of ray tangent. Endpoints of segments are handled
395
-
! so that if the ray moves slightly along its current direction, it will
396
-
! remain in the same segment.
397
-
398
-
INTEGER, PARAMETER :: PRTFile = 6
399
-
INTEGER :: IsegBotT( 1 )
400
-
REAL (KIND=8), INTENT( IN ) :: r, t
401
-
402
364871
IF ( t > 0.0 ) THEN
403
136791051*
IsegBotT = MAXLOC( Bot( : )%x( 1 ), Bot( : )%x( 1 ) <= r )
404
-
ELSE
405
451784*
IsegBotT = MAXLOC( Bot( : )%x( 1 ), Bot( : )%x( 1 ) < r )
406
-
END IF
407
-
408
364871
IF ( IsegBotT( 1 ) > 0 .AND. IsegBotT( 1 ) < NbtyPts ) THEN ! IsegBot MUST LIE IN [ 1, NbtyPts-1 ]
409
364871
IsegBot = IsegBotT( 1 )
410
1094613*
rBotSeg = [ Bot( IsegBot )%x( 1 ), Bot( IsegBot + 1 )%x( 1 ) ] ! segment limits in range
411
-
ELSE
412
#####
WRITE( PRTFile, * ) 'r = ', r
413
#####
WRITE( PRTFile, * ) 'rLeft = ', Bot( 1 )%x( 1 )
414
#####
WRITE( PRTFile, * ) 'rRight = ', Bot( NbtyPts )%x( 1 )
415
#####
CALL ERROUT( 'GetBotSeg', 'Bottom bathymetry undefined below the source' )
416
-
END IF
417
-
418
364871
END SUBROUTINE GetBotSeg
419
-
420
#####
END MODULE bdrymod
420
#####
END MODULE bdrymod
420
#####
END MODULE bdrymod

Legend:

Execution Count: Number of times each line was executed. ##### indicates never executed, - indicates non-executable.