Coverage Report: sspMod.f90

Generated from GCOV analysis of Fortran source code

42.0%
Lines Executed
421 total lines
54.3%
Branches Executed
646 total branches
100.0%
Calls Executed
87 total calls
0
-
Source:sspMod.f90
0
-
Graph:sspMod.gcno
0
-
Data:sspMod.gcda
0
-
Runs:73
1
-
!! Sound speed profile module with interpolation and derivatives
2
-
3
-
MODULE sspmod
4
-
!! Sound speed profile handling with interpolation, derivatives, and environment management
5
-
6
-
! Holds SSP input by user and associated variables
7
-
8
-
! This module is very similar to the one used by the other programs in the Acoustics Toolbox
9
-
! However, it returns the SSP *and* its derivatives
10
-
11
-
! Also, a greater premium has been placed on returning this info quickly, since BELLHOP calls it at every step
12
-
! Therefore more information is pre-computed
13
-
14
-
USE FatalError
15
-
USE monotonicMod
16
-
USE splinec
17
-
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT
18
-
19
-
IMPLICIT NONE
20
-
PUBLIC
21
-
SAVE
22
-
23
-
INTEGER, PARAMETER, PRIVATE :: ENVFile = 5, PRTFile = 6
24
-
INTEGER, PARAMETER :: MaxSSP = 100001
25
-
INTEGER :: iSegr = 1, iSegx = 1, iSegy = 1, iSegz = 1
26
-
INTEGER, PRIVATE :: iz
27
-
REAL (KIND=8), PRIVATE :: Depth, W
28
-
REAL (KIND=8) :: zTemp, betaPowerLaw = 1, fT = 1D20
29
-
REAL (KIND=8) :: alphaR = 1500, betaR = 0, alphaI = 0, betaI = 0, rhoR = 1
30
-
31
-
TYPE rxyz_vector
32
-
REAL(KIND=8), ALLOCATABLE :: r(:), x(:), y(:), z(:)
33
-
END TYPE rxyz_vector
34
-
35
-
! SSP
36
-
TYPE SSPStructure
37
-
INTEGER :: NPts, Nr, Nx, Ny, Nz
38
-
REAL (KIND=8) :: z( MaxSSP ), rho( MaxSSP )
39
-
COMPLEX (KIND=8) :: c( MaxSSP ), cz( MaxSSP ), n2( MaxSSP ), n2z( MaxSSP ), cSpline( 4, MaxSSP )
40
-
COMPLEX (KIND=8) :: cCoef( 4, MaxSSP ), CSWork( 4, MaxSSP ) ! for PCHIP coefs.
41
-
REAL (KIND=8), ALLOCATABLE :: cMat( :, : ), czMat( :, : ), cMat3( :, :, : ), czMat3( :, :, : )
42
-
TYPE ( rxyz_vector ) :: Seg
43
-
CHARACTER (LEN=1) :: Type
44
-
CHARACTER (LEN=2) :: AttenUnit
45
-
END TYPE SSPStructure
46
-
47
-
TYPE( SSPStructure ) :: SSP
48
-
49
-
! *** Halfspace properties structure ***
50
-
51
-
TYPE HSInfo
52
-
REAL (KIND=8) :: alphaR, alphaI, betaR, betaI ! compressional and shear wave speeds/attenuations in user units
53
-
COMPLEX (KIND=8) :: cP, cS ! P-wave, S-wave speeds
54
-
REAL (KIND=8) :: rho, Depth ! density, depth
55
-
CHARACTER (LEN=1) :: BC ! Boundary condition type
56
-
CHARACTER (LEN=6) :: Opt
57
-
END TYPE HSInfo
58
-
59
-
TYPE BdryPt
60
-
TYPE( HSInfo ) :: HS
61
-
END TYPE BdryPt
62
-
63
-
TYPE BdryType
64
-
TYPE( BdryPt ) :: Top, Bot
65
-
END TYPE BdryType
66
-
67
-
TYPE(BdryType) :: Bdry
68
-
69
-
CONTAINS
70
-
71
104612529*
SUBROUTINE EvaluateSSP( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
72
-
!! Evaluates sound speed profile at given location
73
-
74
-
! Call the particular profil routine indicated by the SSP%Type and perform Task
75
-
! Task = 'TAB' then tabulate cp, cs, rhoT
76
-
! Task = 'INI' then initialize
77
-
78
-
REAL (KIND=8), INTENT( IN ) :: freq
79
-
REAL (KIND=8), INTENT( IN ) :: x( 2 ) ! r-z coordinate where SSP is to be evaluated
80
-
REAL (KIND=8), INTENT( IN ) :: t( 2 ) ! ray tangent; for edge cases of updating segments
81
-
CHARACTER ( LEN=3), INTENT( IN ) :: Task
82
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho
83
-
REAL (KIND=8) :: gradc_3d( 3 ), cxx, cyy, cxy, cxz, cyz
84
-
REAL (KIND=8) :: x3( 3 ), t3( 3 )
85
-
86
2959180
SELECT CASE ( SSP%Type )
87
-
CASE ( 'N' ) ! N2-linear profile option
88
90956282
CALL n2Linear( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
89
-
CASE ( 'C' ) ! C-linear profile option
90
90934579
CALL cLinear( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
91
-
CASE ( 'P' ) ! monotone PCHIP ACS profile option
92
10482283
CALL cPCHIP( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
93
-
CASE ( 'S' ) ! Cubic spline profile option
94
10718770
CALL cCubic( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
95
-
CASE ( 'Q' )
96
3173964*
CALL Quad( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
97
-
CASE ( 'H' )
98
#####
IF ( Task == 'TAB' ) THEN
99
#####
WRITE( PRTFile, * ) 'BELLHOP: EvaluateSSP: Hexahedral is not a valid SSP in 2D mode'
100
#####
CALL ERROUT( 'BELLHOP: EvaluateSSP', 'Hexahedral is not a valid SSP in 2D mode')
101
-
END IF
102
-
! this is called by BELLHOP3D only once, during READIN
103
-
! possibly the logic should be changed to call EvaluateSSP2D or 3D
104
#####
x3 = [ 0.0D0, 0.0D0, x( 2 ) ]
105
#####
t3 = [ 0.0D0, 0.0D0, t( 2 ) ]
106
#####
CALL Hexahedral( x3, t3, c, cimag, gradc_3d, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, Task )
107
-
CASE ( 'A' ) ! Analytic profile option
108
#####
CALL Analytic( x, t, c, cimag, gradc, crr, crz, czz, rho )
109
-
CASE DEFAULT
110
#####
WRITE( PRTFile, * ) 'Profile option: ', SSP%Type
111
104612529*
CALL ERROUT( 'BELLHOP: EvaluateSSP', 'Invalid profile option' )
112
-
END SELECT
113
-
114
104612528
END SUBROUTINE EvaluateSSP
115
-
116
-
! **********************************************************************!
117
-
118
#####
SUBROUTINE EvaluateSSP2D( x2D, t2D, c, cimag, gradc, crr, crz, czz, rho, xs, tradial, freq )
119
-
!! Converts cartesian gradients to polar
120
-
121
-
! Called from BELLHOP3D to get a 2D slice out of the 3D SSP
122
-
123
-
REAL (KIND=8), INTENT( IN ) :: x2D( 2 ), t2D( 2 ), xs( 3 ), tradial( 2 ), freq
124
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), czz, crz, crr, rho
125
-
REAL (KIND=8) :: x( 3 ), t( 3 ), gradc3D(3 ), cxx, cyy, cxy, cxz, cyz
126
-
127
-
! convert polar coordinate to cartesian
128
#####
x = [ xs( 1 ) + x2D( 1 ) * tradial( 1 ), xs( 2 ) + x2D( 1 ) * tradial( 2 ), x2D( 2 ) ]
129
#####
t = [ t2D( 1 ) * tradial( 1 ), t2D( 1 ) * tradial( 2 ), t2D( 2 ) ]
130
-
131
#####
CALL EvaluateSSP3D( x, t, c, cimag, gradc3D, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, 'TAB' )
132
-
133
#####
gradc( 1 ) = DOT_PRODUCT( tradial, gradc3D( 1 : 2 ) ) ! r derivative
134
#####
gradc( 2 ) = gradc3D( 3 ) ! z derivative
135
-
136
#####
crz = tradial( 1 ) * cxz + tradial( 2 ) * cyz
137
#####
crr = cxx * ( tradial( 1 ) )**2 + 2.0 * cxy * tradial( 1 ) * tradial( 2 ) + cyy * ( tradial( 2 ) )**2
138
-
139
#####
RETURN
140
-
END SUBROUTINE EvaluateSSP2D
141
-
142
-
! **********************************************************************!
143
-
144
#####
SUBROUTINE EvaluateSSP3D( x, t, c, cimag, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, Task )
145
-
146
-
! Call the particular profil routine indicated by the SSP%Type and perform Task
147
-
! Task = 'TAB' then tabulate cp, cs, rhoT
148
-
! Task = 'INI' then initialize
149
-
150
-
REAL (KIND=8), INTENT( IN ) :: freq
151
-
REAL (KIND=8), INTENT( IN ) :: x( 3 ) ! x-y-z coordinate where SSP is to be evaluated
152
-
REAL (KIND=8), INTENT( IN ) :: t( 3 ) ! ray tangent; for edge cases of updating segments
153
-
CHARACTER ( LEN=3), INTENT( IN ) :: Task
154
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 3 ), cxx, cyy, czz, cxy, cxz, cyz, rho
155
-
REAL (KIND=8) :: x_rz( 2 ), t_rz( 2 ), gradc_rz( 2 ), crr, crz
156
-
157
#####
x_rz = [ 0.0D0, x( 3 ) ] ! convert x-y-z coordinate to cylindrical coordinate
158
#####
t_rz = [ 0.0D0, t( 3 ) ]
159
-
160
#####
SELECT CASE ( SSP%Type )
161
-
CASE ( 'N' )
162
#####
CALL n2Linear( x_rz, t_rz, c, cimag, gradc_rz, crr, crz, czz, rho, freq, Task )
163
-
CASE ( 'C' )
164
#####
CALL cLinear( x_rz, t_rz, c, cimag, gradc_rz, crr, crz, czz, rho, freq, Task )
165
-
CASE ( 'S' )
166
#####
CALL cCubic( x_rz, t_rz, c, cimag, gradc_rz, crr, crz, czz, rho, freq, Task )
167
-
CASE ( 'H' )
168
#####
CALL Hexahedral( x, t, c, cimag, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, Task )
169
-
CASE ( 'A' )
170
#####
CALL Analytic3D( x, t, c, cimag, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho )
171
-
CASE DEFAULT
172
#####
WRITE( PRTFile, * ) 'Profile option: ', SSP%Type
173
#####
CALL ERROUT( 'BELLHOP3D: EvaluateSSP3D', 'Invalid profile option' )
174
-
END SELECT
175
-
176
#####
SELECT CASE ( SSP%Type )
177
-
CASE ( 'N', 'C', 'S' )
178
#####
gradc = [ 0.0D0, 0.0D0, gradc_rz( 2 ) ]
179
-
180
#####
cxx = 0.0D0
181
#####
cyy = 0.0D0
182
#####
cxy = 0.0D0
183
#####
cxz = 0.0D0
184
#####
cyz = 0.0D0
185
-
CASE ( 'H', 'A' )
186
#####
CONTINUE
187
-
CASE DEFAULT
188
#####
WRITE( ERROR_UNIT, * ) 'EvaluateSSP3D: Unknown SSP type for derivative fixup: ', SSP%Type
189
#####
CALL ERROUT( 'BELLHOP3D: EvaluateSSP3D', 'Invalid profile option' )
190
-
END SELECT
191
-
192
#####
END SUBROUTINE EvaluateSSP3D
193
-
194
-
! **********************************************************************!
195
-
196
2959180*
SUBROUTINE n2Linear( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
197
-
!! Linear interpolation for squared buoyancy frequency
198
-
199
-
! N2-linear interpolation of SSP data
200
-
201
-
REAL (KIND=8), INTENT( IN ) :: freq
202
-
REAL (KIND=8), INTENT( IN ) :: x( 2 ) ! r-z coordinate where sound speed is evaluated
203
-
REAL (KIND=8), INTENT( IN ) :: t( 2 ) ! ray tangent; for edge cases of updating segments
204
-
CHARACTER (LEN=3), INTENT( IN ) :: Task
205
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! sound speed and its derivatives
206
-
207
2959180
IF ( Task == 'INI' ) THEN ! read in SSP data
208
2
Depth = x( 2 )
209
2
CALL ReadSSP( Depth, freq )
210
-
211
9*
SSP%n2( 1 : SSP%NPts ) = 1.0 / SSP%c( 1 : SSP%NPts ) ** 2
212
-
213
-
! compute gradient, n2z
214
7*
DO iz = 2, SSP%Npts
215
15*
SSP%n2z( iz - 1 ) = ( SSP%n2( iz ) - SSP%n2( iz - 1 ) ) / &
216
22*
( SSP%z( iz ) - SSP%z( iz - 1 ) )
217
-
END DO
218
-
ELSE ! return SSP info
219
-
220
2959178
CALL UpdateDepthSegmentT( x, t )
221
-
222
2959178*
W = ( x( 2 ) - SSP%z( iSegz ) ) / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) )
223
-
224
2959178*
c = REAL( 1.0D0 / SQRT( ( 1.0D0 - W ) * SSP%n2( iSegz ) + W * SSP%n2( iSegz + 1 ) ) )
225
2959178*
cimag = AIMAG( 1.0D0 / SQRT( ( 1.0D0 - W ) * SSP%n2( iSegz ) + W * SSP%n2( iSegz + 1 ) ) )
226
-
227
8877534*
gradc = [ 0.0D0, -0.5D0 * c * c * c * REAL( SSP%n2z( iSegz ) ) ]
228
2959178
crr = 0.0d0
229
2959178
crz = 0.0d0
230
2959178
czz = 3.0d0 * gradc( 2 ) * gradc( 2 ) / C
231
-
232
2959178*
rho = ( 1.0D0 - W ) * SSP%rho( iSegz ) + W * SSP%rho( iSegz + 1 )
233
-
END IF
234
-
235
2959180
END SUBROUTINE n2Linear
236
-
237
-
! **********************************************************************!
238
-
239
87997102*
SUBROUTINE cLinear( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
240
-
!! c-linear interpolation of SSP data
241
-
242
-
REAL (KIND=8), INTENT( IN ) :: freq
243
-
REAL (KIND=8), INTENT( IN ) :: x( 2 ) ! r-z coordinate where sound speed is evaluated
244
-
REAL (KIND=8), INTENT( IN ) :: t( 2 ) ! ray tangent; for edge cases of updating segments
245
-
CHARACTER (LEN=3), INTENT( IN ) :: Task
246
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! sound speed and its derivatives
247
-
248
87997102
IF ( Task == 'INI' ) THEN ! read in SSP data
249
63
Depth = x( 2 )
250
63
CALL ReadSSP( Depth, freq )
251
-
ELSE ! return SSP info
252
-
253
87997039
CALL UpdateDepthSegmentT( x, t )
254
-
255
87997039*
c = REAL( SSP%c( iSegz ) + ( x( 2 ) - SSP%z( iSegz ) ) * SSP%cz( iSegz ) )
256
87997039*
cimag = AIMAG( SSP%c( iSegz ) + ( x( 2 ) - SSP%z( iSegz ) ) * SSP%cz( iSegz ) )
257
263991117*
gradc = [ 0.0D0, REAL( SSP%cz( iSegz ) ) ]
258
87997039
crr = 0.0d0
259
87997039
crz = 0.0d0
260
87997039
czz = 0.0d0
261
-
262
87997039*
W = ( x( 2 ) - SSP%z( iSegz ) ) / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) )
263
87997039*
rho = ( 1.0D0 - W ) * SSP%rho( iSegz ) + W * SSP%rho( iSegz + 1 )
264
-
END IF
265
-
266
87997101
END SUBROUTINE cLinear
267
-
268
-
! **********************************************************************!
269
-
270
2937477*
SUBROUTINE cPCHIP( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
271
-
!! PCHIP for interpolation of sound speed
272
-
273
-
! This implements the new monotone piecewise cubic Hermite interpolating
274
-
! polynomial (PCHIP) algorithm for the interpolation of the sound speed c.
275
-
276
-
USE pchipMod
277
-
REAL (KIND=8), INTENT( IN ) :: freq
278
-
REAL (KIND=8), INTENT( IN ) :: x( 2 ) ! r-z coordinate where sound speed is evaluated
279
-
REAL (KIND=8), INTENT( IN ) :: t( 2 ) ! ray tangent; for edge cases of updating segments
280
-
CHARACTER (LEN=3), INTENT( IN ) :: Task
281
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! sound speed and its derivatives
282
-
REAL (KIND=8) :: xt
283
-
COMPLEX (KIND=8) :: c_cmplx
284
-
285
2937477
IF ( Task == 'INI' ) THEN ! read in SSP data
286
-
287
1
Depth = x( 2 )
288
1
CALL ReadSSP( Depth, freq )
289
-
290
-
! 2 3
291
-
! compute coefficients of std cubic polynomial: c0 + c1*x + c2*x + c3*x
292
-
!
293
-
294
1
CALL PCHIP( SSP%z, SSP%c, SSP%NPts, SSP%cCoef, SSP%CSWork )
295
-
296
-
ELSE ! return SSP info
297
-
298
2937476
CALL UpdateDepthSegmentT( x, t )
299
-
300
2937476*
xt = x( 2 ) - SSP%z( iSegz )
301
2937476*
c_cmplx = SSP%cCoef( 1, iSegz ) &
302
2937476*
+ ( SSP%cCoef( 2, iSegz ) &
303
2937476*
+ ( SSP%cCoef( 3, iSegz ) &
304
11749904*
+ SSP%cCoef( 4, iSegz ) * xt ) * xt ) * xt
305
-
306
2937476
c = REAL( c_cmplx )
307
2937476
cimag = AIMAG( c_cmplx )
308
-
309
2937476*
gradc = [ 0.0D0, REAL( SSP%cCoef( 2, iSegz ) &
310
2937476*
+ ( 2.0D0 * SSP%cCoef( 3, iSegz ) &
311
14687380*
+ 3.0D0 * SSP%cCoef( 4, iSegz ) * xt ) * xt ) ]
312
-
313
2937476
crr = 0.0D0
314
2937476
crz = 0.0D0
315
2937476*
czz = REAL( 2.0D0 * SSP%cCoef( 3, iSegz ) + 6.0D0 * SSP%cCoef( 4, iSegz ) * xt )
316
-
317
2937476*
W = ( x( 2 ) - SSP%z( iSegz ) ) / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) )
318
2937476*
rho = ( 1.0D0 - W ) * SSP%rho( iSegz ) + W * SSP%rho( iSegz + 1 ) ! linear interp of density
319
-
320
-
END IF
321
-
322
2937477
END SUBROUTINE cPCHIP
323
-
324
-
! **********************************************************************!
325
-
326
7544806*
SUBROUTINE cCubic( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
327
-
!! Cubic spline interpolation for sound speed
328
-
329
-
REAL (KIND=8), INTENT( IN ) :: freq
330
-
REAL (KIND=8), INTENT( IN ) :: x( 2 ) ! r-z coordinate where sound speed is evaluated
331
-
REAL (KIND=8), INTENT( IN ) :: t( 2 ) ! ray tangent; for edge cases of updating segments
332
-
CHARACTER (LEN=3), INTENT( IN ) :: Task
333
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! sound speed and its derivatives
334
-
INTEGER :: iBCBeg, iBCEnd
335
-
REAL (KIND=8) :: hSpline
336
-
COMPLEX (KIND=8) :: c_cmplx, cz_cmplx, czz_cmplx
337
-
338
7544806
IF ( Task == 'INI' ) THEN
339
-
340
-
! *** Task 'INIT' for initialization ***
341
-
342
3
Depth = x( 2 )
343
3
CALL ReadSSP( Depth, freq )
344
-
345
17*
SSP%cSpline( 1, 1 : SSP%NPts ) = SSP%c( 1 : SSP%NPts )
346
-
347
-
! Compute spline coefs
348
3
iBCBeg = 0
349
3
iBCEnd = 0
350
3
CALL CSpline( SSP%z, SSP%cSpline( 1, 1 ), SSP%NPts, iBCBeg, iBCEnd, SSP%NPts )
351
-
ELSE
352
-
353
-
! *** Section to return SSP info ***
354
-
355
7544803
CALL UpdateDepthSegmentT( x, t )
356
-
357
7544803*
hSpline = x( 2 ) - SSP%z( iSegz )
358
-
359
-
! c = Spline( SSP%cSpline( 1, iSegz ), hSpline )
360
-
! cz = SplineX( SSP%cSpline( 1, iSegz ), hSpline )
361
-
! czz = SplineXX( SSP%cSpline( 1, iSegz ), hSpline )
362
-
363
7544803*
CALL SplineALL( SSP%cSpline( 1, iSegz ), hSpline, c_cmplx, cz_cmplx, czz_cmplx )
364
-
365
7544803
c = DBLE( c_cmplx )
366
7544803
cimag = AIMAG( c_cmplx )
367
22634409
gradc = [ 0.0D0, DBLE( cz_cmplx ) ]
368
7544803
czz = DBLE( czz_cmplx )
369
7544803
crr = 0.0d0
370
7544803
crz = 0.0d0
371
-
372
-
! linear interpolation for density
373
7544803*
W = ( x( 2 ) - SSP%z( iSegz ) ) / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) )
374
7544803*
rho = ( 1.0D0 - W ) * SSP%rho( iSegz ) + W * SSP%rho( iSegz + 1 )
375
-
END IF
376
-
377
7544806
END SUBROUTINE cCubic
378
-
379
-
! **********************************************************************!
380
-
381
3173964*
SUBROUTINE Quad( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, Task )
382
-
!! Quadrilateral interpolation for sound speed profiles
383
-
384
-
! Bilinear quadrilateral interpolation of SSP data in 2D
385
-
386
-
INTEGER, PARAMETER :: SSPFile = 40
387
-
REAL (KIND=8), INTENT( IN ) :: freq
388
-
REAL (KIND=8), INTENT( IN ) :: x( 2 ) ! r-z coordinate where sound speed is evaluated
389
-
REAL (KIND=8), INTENT( IN ) :: t( 2 ) ! ray tangent; for edge cases of updating segments
390
-
CHARACTER (LEN=3), INTENT( IN ) :: Task
391
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! sound speed and its derivatives
392
-
INTEGER :: AllocateStatus, iSegT, iz2
393
-
REAL (KIND=8) :: c1, c2, cz1, cz2, cr, cz, s1, s2, delta_r, delta_z
394
-
395
3173964
IF ( Task == 'INI' ) THEN
396
-
397
-
! *** read in SSP data ***
398
-
399
3
Depth = x( 2 )
400
3
CALL ReadSSP( Depth, freq )
401
-
402
-
! Read the 2D SSP matrix
403
3
WRITE( PRTFile, * ) '__________________________________________________________________________'
404
3
WRITE( PRTFile, * )
405
3
WRITE( PRTFile, * ) 'Using range-dependent sound speed'
406
-
407
3
READ( SSPFile, * ) SSP%Nr
408
3
WRITE( PRTFile, * ) 'Number of SSP ranges = ', SSP%Nr
409
-
410
3
IF ( SSP%Nr < 2 ) THEN
411
#####
CALL ERROUT( 'sspMod: Quad', 'You must have at least two profiles in your 2D SSP field' )
412
-
END IF
413
-
414
3*
ALLOCATE( SSP%cMat( SSP%NPts, SSP%Nr ), SSP%czMat( SSP%NPts - 1, SSP%Nr ), SSP%Seg%r( SSP%Nr ), STAT = AllocateStatus )
415
3*
IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Quad', 'Insufficient memory to store SSP' )
416
-
417
3*
READ( SSPFile, * ) SSP%Seg%r( 1 : SSP%Nr )
418
3
WRITE( PRTFile, * )
419
3
WRITE( PRTFile, * ) 'Profile ranges (km):'
420
3*
WRITE( PRTFile, FMT="( F10.2 )" ) SSP%Seg%r( 1 : SSP%Nr )
421
-
422
3*
IF ( .NOT. monotonic( SSP%Seg%r, SSP%Nr ) ) THEN
423
#####
CALL ERROUT( 'sspMod: Quad', 'The ranges in the SSP must be monotone increasing' )
424
-
END IF
425
-
426
93*
SSP%Seg%r = 1000.0 * SSP%Seg%r ! convert km to m
427
-
428
3
WRITE( PRTFile, * )
429
3
WRITE( PRTFile, * ) 'Sound speed matrix:'
430
3
WRITE( PRTFile, * ) ' Depth (m ) Soundspeed (m/s)'
431
9*
DO iz2 = 1, SSP%NPts
432
6*
READ( SSPFile, * ) SSP%cMat( iz2, : )
433
-
! WRITE( PRTFile, FMT="( 'iSegz depth = ', F10.2, ' m' )" ) SSP%z( iz2 )
434
-
! WRITE( PRTFile, FMT="( 12F10.2 )" ) SSP%cMat( iz2, : )
435
9*
WRITE( PRTFile, FMT="( 12F10.2 )" ) SSP%z( iz2 ), SSP%cMat( iz2, : )
436
-
END DO
437
-
438
3
CLOSE( SSPFile )
439
-
440
-
! calculate cz
441
93*
DO iSegt = 1, SSP%Nr
442
183*
DO iz2 = 2, SSP%NPts
443
90*
delta_z = ( SSP%z( iz2 ) - SSP%z( iz2 - 1 ) )
444
180*
SSP%czMat( iz2 - 1, iSegt ) = ( SSP%cMat( iz2, iSegt ) - SSP%cMat( iz2 - 1, iSegt ) ) / delta_z
445
-
END DO
446
-
END DO
447
-
448
3
SSP%Nz = SSP%NPts
449
3
RETURN
450
-
451
-
ELSE
452
-
453
-
! *** Section to return SSP info ***
454
-
455
3173961
CALL UpdateDepthSegmentT( x, t )
456
-
457
-
! Check that x is inside the box where the sound speed is defined
458
3173961*
IF ( x( 1 ) < SSP%Seg%r( 1 ) .OR. x( 1 ) > SSP%Seg%r( SSP%Nr ) ) THEN ! .OR. &
459
#####
WRITE( PRTFile, * ) 'ray is outside the box where the ocean soundspeed is defined'
460
#####
WRITE( PRTFile, * ) ' x = ( r, z ) = ', x
461
#####
CALL ERROUT( 'sspMod: Quad', 'ray is outside the box where the soundspeed is defined' )
462
-
END IF
463
-
464
3173961
CALL UpdateRangeSegmentT( x, t )
465
-
466
-
! for this depth, x( 2 ), get the sound speed at both ends of the segment
467
3173961*
cz1 = SSP%czMat( iSegz, iSegr )
468
3173961*
cz2 = SSP%czMat( iSegz, iSegr + 1 )
469
-
470
3173961*
s2 = x( 2 ) - SSP%z( iSegz )
471
3173961*
delta_z = SSP%z( iSegz + 1 ) - SSP%z( iSegz )
472
-
473
3173961*
c1 = SSP%cMat( iSegz, iSegr ) + s2 * cz1
474
3173961*
c2 = SSP%cMat( iSegz, iSegr + 1 ) + s2 * cz2
475
-
476
-
! s1 = proportional distance of x( 1 ) in range
477
3173961*
delta_r = ( SSP%Seg%r( iSegr + 1 ) - SSP%Seg%r( iSegr ) )
478
3173961*
s1 = ( x( 1 ) - SSP%Seg%r( iSegr ) ) / delta_r
479
3173961
s1 = MIN( s1, 1.0D0 ) ! force piecewise constant extrapolation for points outside the box
480
3173961
s1 = MAX( s1, 0.0D0 ) ! "
481
-
482
3173961
c = ( 1.0D0 - s1 ) * c1 + s1 * c2
483
-
484
-
! interpolate the attenuation !!!! This will use the wrong segment if the ssp in the envil is sampled at different depths
485
3173961
s2 = s2 / delta_z ! convert to a proportional depth in the layer
486
3173961*
cimag = AIMAG( ( 1.0D0 - s2 ) * SSP%c( Isegz ) + s2 * SSP%c( Isegz + 1 ) )
487
-
! volume attenuation is taken from the single c(z) profile
488
-
489
3173961
cz = ( 1.0D0 - s1 ) * cz1 + s1 * cz2
490
-
491
3173961
cr = ( c2 - c1 ) / delta_r
492
3173961
crz = ( cz2 - cz1 ) / delta_r
493
-
494
9521883
gradc = [ cr, cz ]
495
3173961
crr = 0.0D0
496
3173961
czz = 0.0D0
497
-
498
-
! linear interpolation for density
499
3173961*
W = ( x( 2 ) - SSP%z( iSegz ) ) / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) )
500
3173961*
rho = ( 1.0D0 - W ) * SSP%rho( iSegz ) + W * SSP%rho( iSegz + 1 )
501
-
END IF
502
-
503
-
END SUBROUTINE Quad
504
-
505
-
! **********************************************************************!
506
-
507
#####
SUBROUTINE Hexahedral( x, t, c, cimag, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, Task )
508
-
!! Trilinear hexahedral interpolation of SSP data
509
-
510
-
! assumes a rectilinear case (not the most general hexahedral)
511
-
512
-
INTEGER, PARAMETER :: SSPFile = 40
513
-
REAL (KIND=8), INTENT( IN ) :: freq
514
-
REAL (KIND=8), INTENT( IN ) :: x( 3 ) ! x-y-z coordinate where sound speed is evaluated
515
-
REAL (KIND=8), INTENT( IN ) :: t( 3 ) ! ray tangent; for edge cases of updating segments
516
-
CHARACTER (LEN =3), INTENT( IN ) :: Task
517
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 3 ), cxx, cyy, czz, cxy, cxz, cyz, rho
518
-
! sound speed and its derivatives
519
-
INTEGER :: AllocateStatus, iSegxt, iSegyt, iy2, iz2
520
-
REAL (KIND=8) :: c1, c2, c11, c12, c21, c22, cz11, cz12, cz21, cz22, cz1, cz2, &
521
-
cx, cy, cz, s1, s2, s3
522
-
523
#####
IF ( Task == 'INI' ) THEN
524
-
525
-
! *** Section to read in SSP data ***
526
-
527
-
! Read dummy SSP information from the environmental file
528
-
! This is over-ridden by the info in the SSP file
529
-
! However, cz info is used in beam selection
530
#####
Depth = x( 3 )
531
#####
CALL ReadSSP( Depth, freq )
532
-
533
-
! Read the 3D SSP matrix
534
-
535
#####
WRITE( PRTFile, * )
536
#####
WRITE( PRTFile, * ) 'Reading sound speed profile from file'
537
-
538
-
! x coordinates
539
#####
READ( SSPFile, * ) SSP%Nx
540
#####
WRITE( PRTFile, * )
541
#####
WRITE( PRTFile, * ) 'Number of points in x = ', SSP%Nx
542
#####
ALLOCATE( SSP%Seg%x( SSP%Nx ), STAT = AllocateStatus )
543
#####
IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Hexahedral', 'Insufficient memory to store SSP' )
544
#####
READ( SSPFile, * ) SSP%Seg%x
545
-
!WRITE( PRTFile, * )
546
-
!WRITE( PRTFile, * ) 'x-coordinates of SSP (km):'
547
-
!WRITE( PRTFile, FMT="( F10.2 )" ) SSP%Seg%x( 1 : SSP%Nx )
548
#####
IF ( .NOT. monotonic( SSP%Seg%x, SSP%Nx ) ) THEN
549
#####
CALL ERROUT( 'sspMod: Hexahedral', 'The x coordinates in the SSP must be monotone increasing' )
550
-
END IF
551
-
552
-
! y coordinates
553
#####
READ( SSPFile, * ) SSP%Ny
554
#####
WRITE( PRTFile, * )
555
#####
WRITE( PRTFile, * ) 'Number of points in y = ', SSP%Ny
556
#####
ALLOCATE( SSP%Seg%y( SSP%Ny ), STAT = AllocateStatus )
557
#####
IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Hexahedral', 'Insufficient memory to store SSP' )
558
#####
READ( SSPFile, * ) SSP%Seg%y
559
-
!WRITE( PRTFile, * )
560
-
!WRITE( PRTFile, * ) 'y-coordinates of SSP (km):'
561
-
!WRITE( PRTFile, FMT="( F10.2 )" ) SSP%Seg%y( 1 : SSP%Ny )
562
#####
IF ( .NOT. monotonic( SSP%Seg%y, SSP%Ny ) ) THEN
563
#####
CALL ERROUT( 'sspMod: Hexahedral', 'The y coordinates in the SSP must be monotone increasing' )
564
-
END IF
565
-
566
-
! z coordinates
567
#####
READ( SSPFile, * ) SSP%Nz
568
#####
WRITE( PRTFile, * )
569
#####
WRITE( PRTFile, * ) 'Number of points in z = ', SSP%Nz
570
#####
ALLOCATE( SSP%Seg%z( SSP%Nz ), STAT = AllocateStatus )
571
#####
IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Hexahedral', 'Insufficient memory to store SSP' )
572
#####
READ( SSPFile, * ) SSP%Seg%z
573
-
!WRITE( PRTFile, * )
574
-
!WRITE( PRTFile, * ) 'z-coordinates of SSP (km):'
575
-
!WRITE( PRTFile, FMT="( F10.2 )" ) SSP%Seg%z( 1 : SSP%Nz )
576
#####
IF ( .NOT. monotonic( SSP%Seg%z, SSP%Nz ) ) THEN
577
#####
CALL ERROUT( 'sspMod: Hexahedral', 'The z coordinates in the SSP must be monotone increasing' )
578
-
END IF
579
-
580
-
! SSP matrix should be bigger than 2x2x2
581
#####
IF ( SSP%Nx < 2 .OR. SSP%Ny < 2 .OR. SSP%Nz < 2 ) THEN
582
-
CALL ERROUT( 'sspMod: Hexahedral', &
583
#####
'You must have at least two points in x, y, z directions in your 3D SSP field' )
584
-
END IF
585
-
586
#####
IF ( SSP%Nz >= MaxSSP ) THEN
587
-
! LP: SSP%Nz / SSP%Seg%z will get assigned to SSP%NPts / SSP%z.
588
-
CALL ERROUT( 'sspMod: Hexahedral', &
589
#####
'Number of SSP points in Z exceeds limit' )
590
-
END IF
591
-
592
#####
ALLOCATE( SSP%cMat3( SSP%Nx, SSP%Ny, SSP%Nz ), SSP%czMat3( SSP%Nx, SSP%Ny, SSP%Nz - 1 ), STAT = AllocateStatus )
593
#####
IF ( AllocateStatus /= 0 ) CALL ERROUT( 'sspMod: Hexahedral', 'Insufficient memory to store SSP' )
594
-
595
#####
WRITE( PRTFile, * )
596
-
! WRITE( PRTFile, * ) 'Sound speed matrix:'
597
#####
DO iz2 = 1, SSP%Nz
598
#####
DO iy2 = 1, SSP%Ny
599
#####
READ( SSPFile, * ) SSP%cMat3( :, iy2, iz2 )
600
-
! WRITE( PRTFile, FMT="( 'z-section = ', F10.2, ' m' )" ) SSP%Seg%z( iz2 )
601
-
! WRITE( PRTFile, FMT="( 'y-section = ', F10.2, ' km' )" ) SSP%Seg%y( iy2 )
602
-
! WRITE( PRTFile, FMT="( 12F10.2 )" ) SSP%cMat3( :, iy2, iz2 )
603
-
END DO
604
-
END DO
605
-
606
#####
CLOSE( SSPFile )
607
-
608
#####
SSP%Seg%x = 1000.0 * SSP%Seg%x ! convert km to m
609
#####
SSP%Seg%y = 1000.0 * SSP%Seg%y ! convert km to m
610
-
611
-
! calculate cz
612
#####
DO iSegxt = 1, SSP%Nx
613
#####
DO iSegyt = 1, SSP%Ny
614
#####
DO iz2 = 2, SSP%Nz
615
#####
SSP%czMat3( iSegxt, iSegyt, iz2 - 1 ) = &
616
#####
( SSP%cMat3( iSegxt, iSegyt, iz2 ) - SSP%cMat3( iSegxt, iSegyt, iz2 - 1 ) ) / &
617
#####
( SSP%Seg%z( iz2 ) - SSP%Seg%z( iz2 - 1 ) )
618
-
END DO
619
-
END DO
620
-
END DO
621
-
622
-
! over-ride the SSP%z values read in from the environmental file with these new values
623
#####
SSP%NPts = SSP%Nz
624
#####
SSP%z( 1 : SSP%Nz ) = SSP%Seg%z
625
-
626
#####
RETURN
627
-
ELSE
628
-
629
-
! *** Section to return SSP info ***
630
-
631
-
! Check that x is inside the box where the sound speed is defined
632
#####
IF ( x( 1 ) < SSP%Seg%x( 1 ) .OR. x( 1 ) > SSP%Seg%x( SSP%Nx ) .OR. &
633
#####
x( 2 ) < SSP%Seg%y( 1 ) .OR. x( 2 ) > SSP%Seg%y( SSP%Ny ) ) THEN ! .OR. &
634
-
!!$ x( 3 ) < SSP%Seg%z( 1 ) .OR. x( 3 ) > SSP%Seg%z( SSP%Nz ) ) THEN
635
#####
WRITE( PRTFile, * ) 'ray is outside the box where the ocean soundspeed is defined'
636
#####
WRITE( PRTFile, * ) ' x = ( x, y, z ) = ', x
637
#####
CALL ERROUT( 'Hexahedral', 'ray is outside the box where the soundspeed is defined' )
638
-
END IF
639
-
640
#####
CALL Update3DXSegmentT( x, t )
641
#####
CALL Update3DYSegmentT( x, t )
642
#####
CALL Update3DZSegmentT( x, t )
643
-
644
-
! cz at the corners of the current rectangle
645
#####
cz11 = SSP%czMat3( iSegx, iSegy , iSegz )
646
#####
cz12 = SSP%czMat3( iSegx + 1, iSegy , iSegz )
647
#####
cz21 = SSP%czMat3( iSegx, iSegy + 1, iSegz )
648
#####
cz22 = SSP%czMat3( iSegx + 1, iSegy + 1, iSegz )
649
-
650
-
! for this depth, x( 3 ) get the sound speed at the corners of the current rectangle
651
#####
s3 = x( 3 ) - SSP%Seg%z( iSegz )
652
#####
c11 = SSP%cMat3( iSegx, iSegy , iSegz ) + s3 * cz11
653
#####
c21 = SSP%cMat3( iSegx + 1, iSegy , iSegz ) + s3 * cz12
654
#####
c12 = SSP%cMat3( iSegx, iSegy + 1, iSegz ) + s3 * cz21
655
#####
c22 = SSP%cMat3( iSegx + 1, iSegy + 1, iSegz ) + s3 * cz22
656
-
657
-
! s1 = proportional distance of x( 1 ) in x
658
#####
s1 = ( x( 1 ) - SSP%Seg%x( iSegx ) ) / ( SSP%Seg%x( iSegx + 1 ) - SSP%Seg%x( iSegx ) )
659
#####
s1 = MIN( s1, 1.0D0 ) ! force piecewise constant extrapolation for points outside the box
660
#####
s1 = MAX( s1, 0.0D0 ) ! "
661
-
662
-
! s2 = proportional distance of x( 2 ) in y
663
#####
s2 = ( x( 2 ) - SSP%Seg%y( iSegy ) ) / ( SSP%Seg%y( iSegy + 1 ) - SSP%Seg%y( iSegy ) )
664
#####
s2 = MIN( s2, 1.0D0 ) ! force piecewise constant extrapolation for points outside the box
665
#####
s2 = MAX( s2, 0.0D0 ) ! "
666
-
667
-
! interpolate the soundspeed in the x direction, at the two endpoints in y (top and bottom sides of rectangle)
668
#####
c1 = c11 + s1 * ( c21 - c11 )
669
#####
c2 = c12 + s1 * ( c22 - c12 )
670
-
!c = ( 1.0D0 - s2 ) * c1 + s2 * c2 ! interpolation in y
671
#####
cy = ( c2 - c1 ) / ( SSP%Seg%y( iSegy + 1 ) - SSP%Seg%y( iSegy ) )
672
-
673
-
! interpolate the soundspeed in the y direction, at the two endpoints in x (left and right sides of rectangle)
674
#####
c1 = c11 + s2 * ( c12 - c11 )
675
#####
c2 = c21 + s2 * ( c22 - c21 )
676
-
677
#####
c = c1 + s1 * ( c2 - c1 ) ! interpolation in x
678
-
679
-
! interpolate the attenuation !!!! This will use the wrong segment if the ssp in the envfil is sampled at different depths
680
#####
s3 = s3 / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) ) ! convert s3 to a proportional distance in the layer
681
#####
cimag = AIMAG( ( 1.0D0 - s3 ) * SSP%c( Isegz ) + s3 * SSP%c( Isegz + 1 ) )
682
-
! volume attenuation is taken from the single c(z) profile
683
-
684
#####
cx = ( c2 - c1 ) / ( SSP%Seg%x( iSegx + 1 ) - SSP%Seg%x( iSegx ) )
685
-
686
-
! same thing on cz
687
#####
cz1 = cz11 + s2 * ( cz21 - cz11 )
688
#####
cz2 = cz12 + s2 * ( cz22 - cz12 )
689
#####
cz = cz1 + s1 * ( cz2 - cz1 ) ! interpolation in z
690
-
691
-
!gradc = [ cx, cy, cz ]
692
#####
gradc( 1 ) = cx
693
#####
gradc( 2 ) = cy
694
#####
gradc( 3 ) = cz
695
-
696
#####
cxx = 0.0D0
697
#####
cyy = 0.0D0
698
#####
czz = 0.0D0
699
#####
cxy = 0.0D0
700
#####
cxz = 0.0D0
701
#####
cyz = 0.0D0
702
-
! write( *, FMT="( 'SSP', I3, 2X, 4F10.2, F10.4 )" ) layer, x, c, cz
703
-
704
-
! linear interpolation for density
705
#####
W = ( x( 3 ) - SSP%z( iSegz ) ) / ( SSP%z( iSegz + 1 ) - SSP%z( iSegz ) )
706
#####
rho = ( 1.0D0 - W ) * SSP%rho( iSegz ) + W * SSP%rho( iSegz + 1 )
707
-
END IF
708
-
709
-
END SUBROUTINE Hexahedral
710
-
711
-
! **********************************************************************!
712
-
713
#####
SUBROUTINE Analytic( x, t, c, cimag, gradc, crr, crz, czz, rho )
714
-
715
-
REAL (KIND=8), INTENT( IN ) :: x( 2 )
716
-
REAL (KIND=8), INTENT( IN ) :: t( 2 ) ! ray tangent; for edge cases of updating segments
717
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho
718
-
REAL (KIND=8) :: c0, cr, cz, DxtDz, xt
719
-
720
#####
iSegz = 1
721
#####
c0 = 1500.0
722
#####
rho = 1.0
723
-
724
-
! homogeneous halfspace was removed since BELLHOP needs to get gradc just a little below the boundaries, on ray reflection
725
-
726
-
!!$ IF ( x( 2 ) < 5000.0 ) THEN
727
#####
xt = 2.0 * ( x( 2 ) - 1300.0 ) / 1300.0
728
#####
DxtDz = 2.0 / 1300.0
729
#####
c = C0 * ( 1.0 + 0.00737*( xt - 1.0 + EXP( -xt ) ) )
730
#####
cimag = 0.
731
#####
cz = C0 * 0.00737 * ( 1.0 - EXP( -xt ) ) * DxtDz
732
#####
czz = C0 * 0.00737 * EXP( -xt ) * DxtDz ** 2
733
-
!!$ ELSE
734
-
!!$ ! Homogeneous half-space
735
-
!!$ xt = 2.0 * ( 5000.0 - 1300.0 ) / 1300.0
736
-
!!$ c = C0 * ( 1.0 + 0.00737 * ( xt - 1.0 + EXP( -xt ) ) )
737
-
!!$ cimag = 0.0
738
-
!!$ cz = 0.0
739
-
!!$ czz = 0.0
740
-
!!$ END IF
741
-
742
#####
cr = 0.0
743
#####
gradc = [ cr, cz ]
744
#####
crz = 0.0
745
#####
crr = 0.0
746
-
747
#####
RETURN
748
-
END SUBROUTINE Analytic
749
-
750
-
! **********************************************************************!
751
-
752
#####
SUBROUTINE AnalyticCosh( x, t, c, cimag, gradc, crr, crz, czz, rho )
753
-
754
-
REAL (KIND=8), INTENT( IN ) :: x( 2 )
755
-
REAL (KIND=8), INTENT( IN ) :: t( 2 ) ! ray tangent; for edge cases of updating segments
756
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 2 ), crr, crz, czz, rho
757
-
REAL (KIND=8) :: cr, cz, A, B, D, z
758
-
759
#####
iSegz = 1
760
#####
rho = 1.0
761
-
762
#####
D = 3000.0
763
#####
z = x( 2 )
764
-
765
#####
A = 1500.0
766
#####
B = 0.0003
767
#####
D = 1500.0
768
-
769
#####
c = A * COSH( B * ( Z - D ) )
770
#####
cz = B * A * SINH( B * ( Z - D ) )
771
#####
czz = B * B * A * COSH( B * ( Z - D ) )
772
-
773
#####
cimag = 0.
774
-
775
#####
cr = 0.0
776
#####
gradc = [ cr, cz ]
777
#####
crz = 0.0
778
#####
crr = 0.0
779
-
780
#####
RETURN
781
-
END SUBROUTINE AnalyticCosh
782
-
783
-
! **********************************************************************!
784
-
785
#####
SUBROUTINE Analytic3D( x, t, c, cimag, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho )
786
-
787
-
REAL (KIND=8), INTENT( IN ) :: x( 3 )
788
-
REAL (KIND=8), INTENT( OUT ) :: c, cimag, gradc( 3 )
789
-
REAL (KIND=8), INTENT( OUT ) :: cxx, cyy, czz, cxy, cxz, cyz
790
-
REAL (KIND=8), INTENT( IN ) :: t( 3 ) ! ray tangent; for edge cases of updating segments
791
-
REAL (KIND=8), INTENT( INOUT ) :: rho
792
-
793
-
REAL (KIND=8) :: c0, W, Wz, epsilon, epsilon_y
794
-
795
#####
iSegz = 1
796
#####
c0 = 1500.0
797
#####
rho = 1.0
798
-
799
-
!!$ IF ( x( 3 ) .LT. 5000.0 ) THEN
800
#####
epsilon = 0.00737 + x( 2 ) / 100000.0 * 0.003
801
#####
epsilon_y = 0.003 / 100000.0
802
-
803
#####
W = 2.0 * ( x( 3 ) - 1300.0 ) / 1300.0
804
#####
Wz = 2.0 / 1300.0
805
-
806
#####
c = c0 * ( 1.0 + epsilon * ( W - 1.0 + EXP( -W ) ) )
807
#####
cimag = 0.
808
#####
gradc( 2 ) = c0 * epsilon_y * ( W - 1.0 + EXP( -W ) )
809
#####
gradc( 3 ) = c0 * epsilon * ( 1.0 - EXP( -W ) ) * Wz
810
#####
czz = c0 * epsilon * EXP( -W ) * Wz **2
811
#####
cyz = c0 * epsilon_y * ( 1.0 - EXP( -W ) ) * Wz
812
-
!!$ ELSE ! HOMOGENEOUS HALF-SPACE
813
-
!!$ W = 2.0 * ( 5000.0 - 1300.0 ) / 1300.0
814
-
!!$ c = c0*( 1.0 + 0.00737 * ( W - 1.0 + EXP( -W ) ) )
815
-
!!$ gradc( 2 ) = 0.0
816
-
!!$ gradc( 3 ) = 0.0
817
-
!!$ czz = 0.0
818
-
!!$ cyz = 0.0
819
-
!!$ END IF
820
-
821
#####
gradc( 1 ) = 0.0
822
#####
cxx = 0.0
823
#####
cyy = 0.0
824
#####
cxz = 0.0
825
#####
cxy = 0.0
826
-
827
#####
RETURN
828
-
END SUBROUTINE Analytic3D
829
-
830
-
! **********************************************************************!
831
-
832
72
SUBROUTINE ReadSSP( Depth, freq )
833
-
!! Reads SSP data from the environmental file and convert to Nepers/m
834
-
835
-
USE AttenMod
836
-
837
-
REAL (KIND=8), INTENT(IN) :: freq, Depth
838
-
839
72
WRITE( PRTFile, * )
840
72
WRITE( PRTFile, * ) 'Sound speed profile:'
841
-
842
72
WRITE( PRTFile, "( ' z alphaR betaR rho alphaI betaI' )" )
843
72
WRITE( PRTFile, "( ' (m) (m/s) (m/s) (g/cm^3) (m/s) (m/s)', / )" )
844
-
845
72
SSP%NPts = 1
846
-
847
181*
DO iz = 1, MaxSSP
848
-
849
181*
READ( ENVFile, * ) SSP%z( iz ), alphaR, betaR, rhoR, alphaI, betaI
850
181*
WRITE( PRTFile, FMT="( F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )" ) SSP%z( iz ), alphaR, betaR, rhoR, alphaI, betaI
851
-
852
181*
SSP%c( iz ) = CRCI( SSP%z( iz ), alphaR, alphaI, freq, freq, SSP%AttenUnit, betaPowerLaw, fT )
853
180*
SSP%rho( iz ) = rhoR
854
-
855
-
! verify that the depths are monotone increasing
856
180
IF ( iz > 1 ) THEN
857
109*
IF ( SSP%z( iz ) <= SSP%z( iz - 1 ) ) THEN
858
#####
WRITE( PRTFile, * ) 'Bad depth in SSP: ', SSP%z( iz )
859
#####
CALL ERROUT( 'ReadSSP', 'The depths in the SSP must be monotone increasing' )
860
-
END IF
861
-
END IF
862
-
863
-
! compute gradient, cz
864
507*
IF ( iz > 1 ) SSP%cz( iz - 1 ) = ( SSP%c( iz ) - SSP%c( iz - 1 ) ) / &
865
436*
( SSP%z( iz ) - SSP%z( iz - 1 ) )
866
-
867
-
! Did we read the last point?
868
180*
IF ( ABS( SSP%z( iz ) - Depth ) < 100. * EPSILON( 1.0e0 ) ) THEN
869
71
SSP%Nz = SSP%NPts
870
71
IF ( SSP%NPts == 1 ) THEN
871
#####
WRITE( PRTFile, * ) '#SSP points: ', SSP%NPts
872
#####
CALL ERROUT( 'ReadSSP', 'The SSP must have at least 2 points' )
873
-
END IF
874
-
875
71
RETURN
876
-
END IF
877
-
878
109*
SSP%NPts = SSP%NPts + 1
879
-
END DO
880
-
881
-
! Fall through means too many points in the profile
882
#####
WRITE( PRTFile, * ) 'Max. #SSP points: ', MaxSSP
883
#####
CALL ERROUT( 'ReadSSP', 'Number of SSP points exceeds limit' )
884
-
885
-
END SUBROUTINE ReadSSP
886
-
887
104612457
SUBROUTINE UpdateDepthSegmentT( x, t )
888
-
REAL (KIND=8), INTENT(IN) :: x( 2 ), t( 2 )
889
-
890
-
! LP: Handles edge cases based on which direction the ray is going. If the
891
-
! ray takes a small step in the direction of t, it will remain in the same
892
-
! segment as it is now.
893
-
894
104612457
IF ( t( 2 ) >= 0.0 ) THEN
895
-
! SSP%z( iSegz ) <= x( 2 ) < SSP%z( iSegz + 1 )
896
59035938*
DO WHILE ( x( 2 ) < SSP%z( iSegz ) .AND. iSegz > 1 )
897
#####
iSegz = iSegz - 1
898
-
END DO
899
59740624*
DO WHILE ( x( 2 ) >= SSP%z( iSegz + 1 ) .AND. iSegz < SSP%Npts-1 )
900
704686
iSegz = iSegz + 1
901
-
END DO
902
-
ELSE
903
-
! SSP%z( iSegz ) < x( 2 ) <= SSP%z( iSegz + 1 )
904
45576519*
DO WHILE ( x( 2 ) > SSP%z( iSegz + 1 ) .AND. iSegz < SSP%Npts-1 )
905
#####
iSegz = iSegz + 1
906
-
END DO
907
46129314*
DO WHILE ( x( 2 ) <= SSP%z( iSegz ) .AND. iSegz > 1 )
908
552795
iSegz = iSegz - 1
909
-
END DO
910
-
END IF
911
-
912
104612457
END SUBROUTINE UpdateDepthSegmentT
913
-
914
3173961
SUBROUTINE UpdateRangeSegmentT( x, t )
915
-
REAL (KIND=8), INTENT(IN) :: x( 2 ), t( 2 )
916
-
917
-
! LP: See UpdateDepthSegmentT
918
-
919
3173961
IF ( t( 1 ) >= 0.0 ) THEN
920
-
! SSP%Seg%r( iSegr ) <= x( 1 ) < SSP%Seg%r( iSegr + 1 )
921
3150787*
DO WHILE ( x( 1 ) < SSP%Seg%r( iSegr ) .AND. iSegr > 1 )
922
#####
iSegr = iSegr - 1
923
-
END DO
924
3164340*
DO WHILE ( x( 1 ) >= SSP%Seg%r( iSegr + 1 ) .AND. iSegr < SSP%Nr-1 )
925
13553
iSegr = iSegr + 1
926
-
END DO
927
-
ELSE
928
-
! SSP%Seg%r( iSegr ) < x( 1 ) <= SSP%Seg%r( iSegr + 1 )
929
23440*
DO WHILE ( x( 1 ) > SSP%Seg%r( iSegr + 1 ) .AND. iSegr < SSP%Nr-1 )
930
266
iSegr = iSegr + 1
931
-
END DO
932
30224*
DO WHILE ( x( 1 ) <= SSP%Seg%r( iSegr ) .AND. iSegr > 1 )
933
7050
iSegr = iSegr - 1
934
-
END DO
935
-
END IF
936
-
937
3173961
END SUBROUTINE UpdateRangeSegmentT
938
-
939
#####
SUBROUTINE Update3DXSegmentT( x, t )
940
-
REAL (KIND=8), INTENT(IN) :: x( 3 ), t( 3 )
941
-
942
-
! LP: See UpdateDepthSegmentT
943
-
944
#####
IF ( t( 1 ) >= 0.0 ) THEN
945
-
! SSP%Seg%x( iSegx ) <= x( 1 ) < SSP%Seg%x( iSegx + 1 )
946
#####
DO WHILE ( x( 1 ) < SSP%Seg%x( iSegx ) .AND. iSegx > 1 )
947
#####
iSegx = iSegx - 1
948
-
END DO
949
#####
DO WHILE ( x( 1 ) >= SSP%Seg%x( iSegx + 1 ) .AND. iSegx < SSP%Nx-1 )
950
#####
iSegx = iSegx + 1
951
-
END DO
952
-
ELSE
953
-
! SSP%Seg%x( iSegx ) < x( 1 ) <= SSP%Seg%x( iSegx + 1 )
954
#####
DO WHILE ( x( 1 ) > SSP%Seg%x( iSegx + 1 ) .AND. iSegx < SSP%Nx-1 )
955
#####
iSegx = iSegx + 1
956
-
END DO
957
#####
DO WHILE ( x( 1 ) <= SSP%Seg%x( iSegx ) .AND. iSegx > 1 )
958
#####
iSegx = iSegx - 1
959
-
END DO
960
-
END IF
961
-
962
#####
END SUBROUTINE Update3DXSegmentT
963
-
964
#####
SUBROUTINE Update3DYSegmentT( x, t )
965
-
REAL (KIND=8), INTENT(IN) :: x( 3 ), t( 3 )
966
-
967
-
! LP: See UpdateDepthSegmentT
968
-
969
#####
IF ( t( 2 ) >= 0.0 ) THEN
970
-
! SSP%Seg%y( iSegy ) <= x( 2 ) < SSP%Seg%y( iSegy + 1 )
971
#####
DO WHILE ( x( 2 ) < SSP%Seg%y( iSegy ) .AND. iSegy > 1 )
972
#####
iSegy = iSegy - 1
973
-
END DO
974
#####
DO WHILE ( x( 2 ) >= SSP%Seg%y( iSegy + 1 ) .AND. iSegy < SSP%Ny-1 )
975
#####
iSegy = iSegy + 1
976
-
END DO
977
-
ELSE
978
-
! SSP%Seg%y( iSegy ) < x( 2 ) <= SSP%Seg%y( iSegy + 1 )
979
#####
DO WHILE ( x( 2 ) > SSP%Seg%y( iSegy + 1 ) .AND. iSegy < SSP%Ny-1 )
980
#####
iSegy = iSegy + 1
981
-
END DO
982
#####
DO WHILE ( x( 2 ) <= SSP%Seg%y( iSegy ) .AND. iSegy > 1 )
983
#####
iSegy = iSegy - 1
984
-
END DO
985
-
END IF
986
-
987
#####
END SUBROUTINE Update3DYSegmentT
988
-
989
#####
SUBROUTINE Update3DZSegmentT( x, t )
990
-
REAL (KIND=8), INTENT(IN) :: x( 3 ), t( 3 )
991
-
992
-
! LP: See UpdateDepthSegmentT
993
-
994
#####
IF ( t( 3 ) >= 0.0 ) THEN
995
-
! SSP%Seg%z( iSegz ) <= x( 3 ) < SSP%Seg%z( iSegz + 1 )
996
#####
DO WHILE ( x( 3 ) < SSP%Seg%z( iSegz ) .AND. iSegz > 1 )
997
#####
iSegz = iSegz - 1
998
-
END DO
999
#####
DO WHILE ( x( 3 ) >= SSP%Seg%z( iSegz + 1 ) .AND. iSegz < SSP%Nz-1 )
1000
#####
iSegz = iSegz + 1
1001
-
END DO
1002
-
ELSE
1003
-
! SSP%Seg%z( iSegz ) < x( 3 ) <= SSP%Seg%z( iSegz + 1 )
1004
#####
DO WHILE ( x( 3 ) > SSP%Seg%z( iSegz + 1 ) .AND. iSegz < SSP%Nz-1 )
1005
#####
iSegz = iSegz + 1
1006
-
END DO
1007
#####
DO WHILE ( x( 3 ) <= SSP%Seg%z( iSegz ) .AND. iSegz > 1 )
1008
#####
iSegz = iSegz - 1
1009
-
END DO
1010
-
END IF
1011
-
1012
#####
END SUBROUTINE Update3DZSegmentT
1013
-
1014
#####
END MODULE sspmod
1014
#####
END MODULE sspmod
1014
#####
END MODULE sspmod
1014
#####
END MODULE sspmod
1014
#####
END MODULE sspmod
1014
#####
END MODULE sspmod
1014
#####
END MODULE sspmod
1014
#####
END MODULE sspmod

Legend:

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