Coverage Report: sspMod.f90

Generated from GCOV analysis of Fortran source code

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

Legend:

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