Coverage Report: ReadEnvironmentBell.f90

Generated from GCOV analysis of Fortran source code

62.7%
Lines Executed
303 total lines
52.6%
Branches Executed
173 total branches
100.0%
Calls Executed
311 total calls
0
-
Source:ReadEnvironmentBell.f90
0
-
Graph:ReadEnvironmentBell.gcno
0
-
Data:ReadEnvironmentBell.gcda
0
-
Runs:75
1
-
!! Environment input file parsing and setup for BELLHOP
2
-
3
-
MODULE ReadEnvironmentBell
4
-
!! Provides environment file reading and initialization
5
-
6
-
USE BellhopMod
7
-
USE MathConstants
8
-
USE sspmod
9
-
USE AttenMod
10
-
USE FatalError
11
-
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT
12
-
13
-
IMPLICIT NONE
14
-
PUBLIC
15
-
16
-
CONTAINS
17
-
18
75*
SUBROUTINE ReadEnvironment( FileRoot, ThreeD )
19
-
!! Reads and parses the main environment file
20
-
21
-
! Routine to read in and echo all the input data
22
-
! Note that default values of SSP, DENSITY, Attenuation will not work
23
-
24
-
USE anglemod
25
-
USE SourceReceiverPositions
26
-
27
-
REAL (KIND=8), PARAMETER :: c0 = 1500.0
28
-
LOGICAL, INTENT(IN ) :: ThreeD
29
-
CHARACTER (LEN=80), INTENT(IN ) :: FileRoot
30
-
INTEGER :: NPts, NMedia, iostat
31
-
REAL :: ZMin, ZMax
32
-
REAL (KIND=8) :: x( 2 ), t( 2 ), c, cimag, gradc( 2 ), crr, crz, czz, rho, sigma, Depth
33
-
CHARACTER (LEN= 2) :: AttenUnit
34
-
CHARACTER (LEN=10) :: PlotType
35
-
36
75
WRITE( PRTFile, * ) 'BELLHOP/BELLHOP3D'
37
75
WRITE( PRTFile, * )
38
-
39
-
! Open the environmental file
40
75*
OPEN( UNIT = ENVFile, FILE = TRIM( FileRoot ) // '.env', STATUS = 'OLD', IOSTAT = iostat, ACTION = 'READ' )
41
75
IF ( IOSTAT /= 0 ) THEN ! successful open?
42
#####
WRITE( PRTFile, * ) 'ENVFile = ', TRIM( FileRoot ) // '.env'
43
#####
CALL ERROUT( 'BELLHOP - READIN', 'Unable to open the environmental file' )
44
-
END IF
45
-
46
-
! Prepend model name to title
47
75
IF ( ThreeD ) THEN
48
#####
Title( 1 : 11 ) = 'BELLHOP3D- '
49
#####
READ( ENVFile, * ) Title( 12 : 80 )
50
-
ELSE
51
75
Title( 1 : 9 ) = 'BELLHOP- '
52
75
READ( ENVFile, * ) Title( 10 : 80 )
53
-
END IF
54
-
55
75
WRITE( PRTFile, * ) Title
56
-
57
75
READ( ENVFile, * ) freq
58
75
WRITE( PRTFile, '('' frequency = '', G11.4, '' Hz'', / )' ) freq
59
-
60
75
READ( ENVFile, * ) NMedia
61
75
WRITE( PRTFile, * ) 'Dummy parameter NMedia = ', NMedia
62
75
IF ( NMedia /= 1 ) CALL ERROUT( 'READIN', &
63
#####
'Only one medium or layer is allowed in BELLHOP; sediment layers must be handled using a reflection coefficient' )
64
-
65
75
CALL ReadTopOpt( Bdry%Top%HS%Opt, Bdry%Top%HS%BC, AttenUnit, FileRoot )
66
-
67
-
! *** Top BC ***
68
-
69
75
IF ( Bdry%Top%HS%BC == 'A' ) THEN
70
3
WRITE( PRTFile, "( //, ' z alphaR betaR rho alphaI betaI', / )" )
71
3
WRITE( PRTFile, "( ' (m) (m/s) (m/s) (g/cm^3) (m/s) (m/s)', / )" )
72
-
END IF
73
-
74
75
CALL TopBot( freq, AttenUnit, Bdry%Top%HS )
75
-
76
-
! ****** Read in ocean SSP data ******
77
-
78
74
READ( ENVFile, * ) NPts, Sigma, Bdry%Bot%HS%Depth
79
74
WRITE( PRTFile, * )
80
74
WRITE( PRTFile, FMT = "( ' Depth = ', F10.2, ' m' )" ) Bdry%Bot%HS%Depth
81
-
82
74
IF ( Bdry%Top%HS%Opt( 1 : 1 ) == 'A' ) THEN
83
#####
WRITE( PRTFile, * ) 'Analytic SSP option'
84
-
! following is hokey, should be set in Analytic routine
85
#####
SSP%NPts = 2
86
#####
SSP%z( 1 ) = 0.0
87
#####
SSP%z( 2 ) = Bdry%Bot%HS%Depth
88
-
ELSE
89
222
x = [ 0.0D0, Bdry%Bot%HS%Depth ] ! tells SSP Depth to read to
90
74
t = [ 0.0, 1.0 ]
91
74
CALL EvaluateSSP( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, 'INI' )
92
-
END IF
93
-
94
73
Bdry%Top%HS%Depth = SSP%z( 1 ) ! Depth of top boundary is taken from first SSP point
95
-
! bottom depth should perhaps be set the same way?
96
-
97
-
! *** Bottom BC ***
98
73
Bdry%Bot%HS%Opt = ' ' ! initialize to blanks
99
73
READ( ENVFile, * ) Bdry%Bot%HS%Opt, Sigma
100
73
WRITE( PRTFile, * )
101
73
WRITE( PRTFile, FMT = "(33X, '( RMS roughness = ', G10.3, ' )' )" ) Sigma
102
-
103
11
SELECT CASE ( Bdry%Bot%HS%Opt( 2 : 2 ) )
104
-
CASE ( '~', '*' )
105
73*
WRITE( PRTFile, * ) ' Bathymetry file selected'
106
-
CASE( '-', '_', ' ' )
107
-
CASE DEFAULT
108
73*
CALL ERROUT( 'READIN', 'Unknown bottom option letter in second position' )
109
-
END SELECT
110
-
111
73
Bdry%Bot%HS%BC = Bdry%Bot%HS%Opt( 1 : 1 )
112
73
CALL TopBot( freq, AttenUnit, Bdry%Bot%HS )
113
-
114
-
! *** source and receiver locations ***
115
-
116
73
CALL ReadSxSy( ThreeD ) ! Read source/receiver x-y coordinates
117
-
118
73
ZMin = SNGL( Bdry%Top%HS%Depth )
119
73
ZMax = SNGL( Bdry%Bot%HS%Depth )
120
-
! CALL ReadSzRz( ZMin + 100 * SPACING( ZMin ), ZMax - 100 * SPACING( ZMax ) ) ! not sure why I had this
121
73
CALL ReadSzRz( ZMin, ZMax )
122
73
CALL ReadRcvrRanges
123
73*
IF ( ThreeD ) CALL ReadRcvrBearings
124
73
CALL ReadfreqVec( freq, Bdry%Top%HS%Opt( 6:6 ) )
125
73
CALL ReadRunType( Beam%RunType, PlotType )
126
-
127
73
Depth = Zmax - Zmin ! water depth
128
73
CALL ReadRayElevationAngles( freq, Depth, Bdry%Top%HS%Opt, Beam%RunType )
129
73*
IF ( ThreeD ) CALL ReadRayBearingAngles( freq, Bdry%Top%HS%Opt, Beam%RunType )
130
-
131
73
WRITE( PRTFile, * )
132
73
WRITE( PRTFile, * ) '__________________________________________________________________________'
133
73
WRITE( PRTFile, * )
134
-
135
-
! Limits for tracing beams
136
-
137
73
IF ( ThreeD ) THEN
138
#####
READ( ENVFile, * ) Beam%deltas, Beam%Box%x, Beam%Box%y, Beam%Box%z
139
#####
Beam%Box%x = 1000.0 * Beam%Box%x ! convert km to m
140
#####
Beam%Box%y = 1000.0 * Beam%Box%y ! convert km to m
141
-
142
#####
IF ( Beam%deltas == 0.0 ) Beam%deltas = ( Bdry%Bot%HS%Depth - Bdry%Top%HS%Depth ) / 10.0 ! Automatic step size selection
143
-
! WRITE( PRTFile, '('' frequency = '', G11.4, '' Hz'', / )' ) freq
144
-
145
#####
WRITE( PRTFile, * )
146
#####
WRITE( PRTFile, fmt = '( '' Step length, deltas = '', G11.4, '' m'' )' ) Beam%deltas
147
#####
WRITE( PRTFile, * )
148
#####
WRITE( PRTFile, fmt = '( '' Maximum ray x-range, Box%x = '', G11.4, '' m'' )' ) Beam%Box%x
149
#####
WRITE( PRTFile, fmt = '( '' Maximum ray y-range, Box%y = '', G11.4, '' m'' )' ) Beam%Box%y
150
#####
WRITE( PRTFile, fmt = '( '' Maximum ray z-range, Box%z = '', G11.4, '' m'' )' ) Beam%Box%z
151
-
ELSE
152
73
READ( ENVFile, * ) Beam%deltas, Beam%Box%z, Beam%Box%r
153
-
154
73
WRITE( PRTFile, * )
155
73
WRITE( PRTFile, fmt = '( '' Step length, deltas = '', G11.4, '' m'' )' ) Beam%deltas
156
73
WRITE( PRTFile, * )
157
73
WRITE( PRTFile, fmt = '( '' Maximum ray depth, Box%z = '', G11.4, '' m'' )' ) Beam%Box%z
158
73
WRITE( PRTFile, fmt = '( '' Maximum ray range, Box%r = '', G11.4, ''km'' )' ) Beam%Box%r
159
-
160
73
Beam%Box%r = 1000.0 * Beam%Box%r ! convert km to m
161
-
END IF
162
-
163
-
! *** Beam characteristics ***
164
73
Beam%Type( 4 : 4 ) = Beam%RunType( 7 : 7 ) ! selects beam shift option
165
-
166
#####
SELECT CASE ( Beam%Type( 4 : 4 ) )
167
-
CASE ( 'S' )
168
#####
WRITE( PRTFile, * ) 'Beam shift in effect'
169
-
CASE DEFAULT
170
73
WRITE( PRTFile, * ) 'No beam shift in effect'
171
-
END SELECT
172
-
173
73
IF ( Beam%RunType( 1 : 1 ) /= 'R' ) THEN ! no worry about the beam type if this is a ray trace run
174
-
175
-
! Beam%Type( 1 : 1 ) is
176
-
! 'G' or '^' Geometric hat beams in Cartesian coordinates
177
-
! 'g' Geometric hat beams in ray-centered coordinates
178
-
! 'B' Geometric Gaussian beams in Cartesian coordinates
179
-
! 'b' Geometric Gaussian beams in ray-centered coordinates
180
-
! 'S' Simple Gaussian beams
181
-
! 'C' Cerveny Gaussian beams in Cartesian coordinates
182
-
! 'R' Cerveny Gaussian beams in Ray-centered coordinates
183
-
! Beam%Type( 2 : 2 ) controls the setting of the beam width
184
-
! 'F' space Filling
185
-
! 'M' minimum width
186
-
! 'W' WKB beams
187
-
! Beam%Type( 3 : 3 ) controls curvature changes on boundary reflections
188
-
! 'D' Double
189
-
! 'S' Single
190
-
! 'Z' Zero
191
-
! Beam%Type( 4 : 4 ) selects whether beam shifts are implemented on boundary reflection
192
-
! 'S' yes
193
-
! 'N' no
194
-
195
-
! Curvature change can cause overflow in grazing case
196
-
! Suppress by setting BeamType( 3 : 3 ) = 'Z'
197
-
198
62
Beam%Type( 1 : 1 ) = Beam%RunType( 2 : 2 )
199
60
SELECT CASE ( Beam%Type( 1 : 1 ) )
200
-
CASE ( 'G', 'g' , '^', 'B', 'b', 'S' ) ! geometric hat beams, geometric Gaussian beams, or simple Gaussian beams
201
-
CASE ( 'R', 'C' ) ! Cerveny Gaussian Beams; read extra lines to specify the beam options
202
2
READ( ENVFile, * ) Beam%Type( 2 : 3 ), Beam%epsMultiplier, Beam%rLoop
203
2
WRITE( PRTFile, * )
204
2
WRITE( PRTFile, * )
205
2
WRITE( PRTFile, * ) 'Type of beam = ', Beam%Type( 1 : 1 )
206
-
207
#####
SELECT CASE ( Beam%Type( 3 : 3 ) )
208
-
CASE ( 'D' )
209
#####
WRITE( PRTFile, * ) 'Curvature doubling invoked'
210
-
CASE ( 'Z' )
211
#####
WRITE( PRTFile, * ) 'Curvature zeroing invoked'
212
-
CASE ( 'S' )
213
2*
WRITE( PRTFile, * ) 'Standard curvature condition'
214
-
CASE DEFAULT
215
2*
CALL ERROUT( 'READIN', 'Unknown curvature condition' )
216
-
END SELECT
217
-
218
2
WRITE( PRTFile, * ) 'Epsilon multiplier', Beam%epsMultiplier
219
2
WRITE( PRTFile, * ) 'Range for choosing beam width', Beam%rLoop
220
-
221
-
! Images, windows
222
-
! LP: These values are not initialized if not written in the file,
223
-
! and Component is not always written in the test env files.
224
2
Beam%Nimage = 1
225
2
Beam%iBeamWindow = 4
226
2
Beam%Component = 'P'
227
2
READ( ENVFile, * ) Beam%Nimage, Beam%iBeamWindow, Beam%Component
228
2
WRITE( PRTFile, * )
229
2
WRITE( PRTFile, * ) 'Number of images, Nimage = ', Beam%Nimage
230
2
WRITE( PRTFile, * ) 'Beam windowing parameter = ', Beam%iBeamWindow
231
2*
WRITE( PRTFile, * ) 'Component = ', Beam%Component
232
-
CASE DEFAULT
233
62*
CALL ERROUT( 'READIN', 'Unknown beam type (second letter of run type)' )
234
-
END SELECT
235
-
END IF
236
-
237
73
WRITE( PRTFile, * )
238
73
CLOSE( ENVFile )
239
-
240
73
END SUBROUTINE ReadEnvironment
241
-
242
-
! **********************************************************************!
243
-
244
75*
SUBROUTINE ReadTopOpt( TopOpt, BC, AttenUnit, FileRoot )
245
-
246
-
CHARACTER (LEN= 6), INTENT( OUT ) :: TopOpt
247
-
CHARACTER (LEN= 1), INTENT( OUT ) :: BC ! Boundary condition type
248
-
CHARACTER (LEN= 2), INTENT( OUT ) :: AttenUnit
249
-
CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
250
-
INTEGER :: iostat
251
-
252
75
TopOpt = ' ' ! initialize to blanks
253
75
READ( ENVFile, * ) TopOpt
254
75
WRITE( PRTFile, * )
255
-
256
75
SSP%Type = TopOpt( 1 : 1 )
257
75
BC = TopOpt( 2 : 2 )
258
75
AttenUnit = TopOpt( 3 : 4 )
259
75
SSP%AttenUnit = AttenUnit
260
-
261
-
! SSP approximation options
262
-
263
2
SELECT CASE ( SSP%Type )
264
-
CASE ( 'N' )
265
2
WRITE( PRTFile, * ) ' N2-linear approximation to SSP'
266
-
CASE ( 'C' )
267
66
WRITE( PRTFile, * ) ' C-linear approximation to SSP'
268
-
CASE ( 'P' )
269
1
WRITE( PRTFile, * ) ' PCHIP approximation to SSP'
270
-
CASE ( 'S' )
271
3
WRITE( PRTFile, * ) ' Spline approximation to SSP'
272
-
CASE ( 'Q' )
273
3
WRITE( PRTFile, * ) ' Quad approximation to SSP'
274
3*
OPEN ( FILE = TRIM( FileRoot ) // '.ssp', UNIT = SSPFile, FORM = 'FORMATTED', STATUS = 'OLD', IOSTAT = iostat )
275
3*
IF ( IOSTAT /= 0 ) THEN ! successful open?
276
#####
WRITE( PRTFile, * ) 'SSPFile = ', TRIM( FileRoot ) // '.ssp'
277
#####
CALL ERROUT( 'BELLHOP - READIN', 'Unable to open the SSP file' )
278
-
END IF
279
-
CASE ( 'H' )
280
#####
WRITE( PRTFile, * ) ' Hexahedral approximation to SSP'
281
#####
OPEN ( FILE = TRIM( FileRoot ) // '.ssp', UNIT = SSPFile, FORM = 'FORMATTED', STATUS = 'OLD', IOSTAT = iostat )
282
#####
IF ( IOSTAT /= 0 ) THEN ! successful open?
283
#####
WRITE( PRTFile, * ) 'SSPFile = ', TRIM( FileRoot ) // '.ssp'
284
#####
CALL ERROUT( 'BELLHOP - READIN', 'Unable to open the SSP file' )
285
-
END IF
286
-
CASE ( 'A' )
287
#####
WRITE( PRTFile, * ) ' Analytic SSP option'
288
-
CASE DEFAULT
289
75*
CALL ERROUT( 'READIN', 'Unknown option for SSP approximation' )
290
-
END SELECT
291
-
292
-
! Attenuation options
293
-
294
#####
SELECT CASE ( AttenUnit( 1 : 1 ) )
295
-
CASE ( 'N' )
296
#####
WRITE( PRTFile, * ) ' Attenuation units: nepers/m'
297
-
CASE ( 'F' )
298
68*
WRITE( PRTFile, * ) ' Attenuation units: dB/mkHz'
299
-
CASE ( 'M' )
300
#####
WRITE( PRTFile, * ) ' Attenuation units: dB/m'
301
-
CASE ( 'W' )
302
7*
WRITE( PRTFile, * ) ' Attenuation units: dB/wavelength'
303
-
CASE ( 'Q' )
304
#####
WRITE( PRTFile, * ) ' Attenuation units: Q'
305
-
CASE ( 'L' )
306
#####
WRITE( PRTFile, * ) ' Attenuation units: Loss parameter'
307
-
CASE DEFAULT
308
75*
CALL ERROUT( 'READIN', 'Unknown attenuation units' )
309
-
END SELECT
310
-
311
-
! optional addition of volume attenuation using standard formulas
312
-
313
#####
SELECT CASE ( AttenUnit( 2 : 2 ) )
314
-
CASE ( 'T' )
315
#####
WRITE( PRTFile, * ) ' THORP volume attenuation added'
316
-
CASE ( 'F' )
317
2
WRITE( PRTFile, * ) ' Francois-Garrison volume attenuation added'
318
2
READ( ENVFile, * ) T, Salinity, pH, z_bar
319
-
WRITE( PRTFile, "( ' T = ', G11.4, 'degrees S = ', G11.4, ' psu pH = ', G11.4, ' z_bar = ', G11.4, ' m' )" ) &
320
2*
T, Salinity, pH, z_bar
321
-
CASE ( 'B' )
322
#####
WRITE( PRTFile, * ) ' Biological attenaution'
323
#####
READ( ENVFile, * ) NBioLayers
324
#####
WRITE( PRTFile, * ) ' Number of Bio Layers = ', NBioLayers
325
-
326
73*
DO iBio = 1, NBioLayers
327
#####
READ( ENVFile, * ) bio( iBio )%Z1, bio( iBio )%Z2, bio( iBio )%f0, bio( iBio )%Q, bio( iBio )%a0
328
#####
WRITE( PRTFile, * ) ' Top of layer = ', bio( iBio )%Z1, ' m'
329
#####
WRITE( PRTFile, * ) ' Bottom of layer = ', bio( iBio )%Z2, ' m'
330
#####
WRITE( PRTFile, * ) ' Resonance frequency = ', bio( iBio )%f0, ' Hz'
331
#####
WRITE( PRTFile, * ) ' Q = ', bio( iBio )%Q
332
#####
WRITE( PRTFile, * ) ' a0 = ', bio( iBio )%a0
333
-
END DO
334
-
CASE ( ' ' )
335
-
CASE DEFAULT
336
75*
CALL ERROUT( 'READIN', 'Unknown top option letter in fourth position' )
337
-
END SELECT
338
-
339
6
SELECT CASE ( TopOpt( 5 : 5 ) )
340
-
CASE ( '~', '*' )
341
75*
WRITE( PRTFile, * ) ' Altimetry file selected'
342
-
CASE ( '-', '_', ' ' )
343
-
CASE DEFAULT
344
75*
CALL ERROUT( 'READIN', 'Unknown top option letter in fifth position' )
345
-
END SELECT
346
-
347
1
SELECT CASE ( TopOpt( 6 : 6 ) )
348
-
CASE ( 'I' )
349
75*
WRITE( PRTFile, * ) ' Development options enabled'
350
-
CASE ( ' ' )
351
-
CASE DEFAULT
352
75*
CALL ERROUT( 'READIN', 'Unknown top option letter in sixth position' )
353
-
END SELECT
354
-
355
75
END SUBROUTINE ReadTopOpt
356
-
357
-
! **********************************************************************!
358
-
359
73*
SUBROUTINE ReadRunType( RunType, PlotType )
360
-
!! Reads and validates the run type parameters
361
-
362
-
! Read the RunType variable and echo with explanatory information to the print file
363
-
364
-
USE SourceReceiverPositions
365
-
366
-
CHARACTER (LEN= 7), INTENT( OUT ) :: RunType
367
-
CHARACTER (LEN=10), INTENT( OUT ) :: PlotType
368
-
369
73
READ( ENVFile, * ) RunType
370
73
WRITE( PRTFile, * )
371
-
372
11
SELECT CASE ( RunType( 1 : 1 ) )
373
-
CASE ( 'R' )
374
11
WRITE( PRTFile, * ) 'Ray trace run'
375
-
CASE ( 'E' )
376
18*
WRITE( PRTFile, * ) 'Eigenray trace run'
377
-
CASE ( 'I' )
378
#####
WRITE( PRTFile, * ) 'Incoherent TL calculation'
379
-
CASE ( 'S' )
380
#####
WRITE( PRTFile, * ) 'Semi-coherent TL calculation'
381
-
CASE ( 'C' )
382
18
WRITE( PRTFile, * ) 'Coherent TL calculation'
383
-
CASE ( 'A' )
384
26*
WRITE( PRTFile, * ) 'Arrivals calculation, ASCII file output'
385
-
CASE ( 'a' )
386
#####
WRITE( PRTFile, * ) 'Arrivals calculation, binary file output'
387
-
CASE DEFAULT
388
73*
CALL ERROUT( 'READIN', 'Unknown RunType selected' )
389
-
END SELECT
390
-
391
1
SELECT CASE ( RunType( 2 : 2 ) )
392
-
CASE ( 'C' )
393
1
WRITE( PRTFile, * ) 'Cartesian beams'
394
-
CASE ( 'R' )
395
1*
WRITE( PRTFile, * ) 'Ray centered beams'
396
-
CASE ( 'S' )
397
#####
WRITE( PRTFile, * ) 'Simple gaussian beams'
398
-
CASE ( 'b' )
399
#####
WRITE( PRTFile, * ) 'Geometric gaussian beams in ray-centered coordinates'
400
-
CASE ( 'B' )
401
3*
WRITE( PRTFile, * ) 'Geometric gaussian beams in Cartesian coordinates'
402
-
CASE ( 'g' )
403
#####
WRITE( PRTFile, * ) 'Geometric hat beams in ray-centered coordinates'
404
-
CASE DEFAULT
405
68
RunType( 2 : 2 ) = 'G'
406
73
WRITE( PRTFile, * ) 'Geometric hat beams in Cartesian coordinates'
407
-
END SELECT
408
-
409
2
SELECT CASE ( RunType( 4 : 4 ) )
410
-
CASE ( 'R' )
411
2
WRITE( PRTFile, * ) 'Point source (cylindrical coordinates)'
412
-
CASE ( 'X' )
413
1
WRITE( PRTFile, * ) 'Line source (Cartesian coordinates)'
414
-
CASE DEFAULT
415
70
RunType( 4 : 4 ) = 'R'
416
73
WRITE( PRTFile, * ) 'Point source (cylindrical coordinates)'
417
-
END SELECT
418
-
419
1
SELECT CASE ( RunType( 5 : 5 ) )
420
-
CASE ( 'R' )
421
1
WRITE( PRTFile, * ) 'Rectilinear receiver grid: Receivers at ( Rr( ir ), Rz( ir ) ) )'
422
1*
PlotType = 'rectilin '
423
-
CASE ( 'I' )
424
#####
WRITE( PRTFile, * ) 'Irregular grid: Receivers at Rr( : ) x Rz( : )'
425
#####
IF ( Pos%NRz /= Pos%NRr ) CALL ERROUT( 'READIN', 'Irregular grid option selected with NRz not equal to Nr' )
426
#####
PlotType = 'irregular '
427
-
CASE DEFAULT
428
72
WRITE( PRTFile, * ) 'Rectilinear receiver grid: Receivers at Rr( : ) x Rz( : )'
429
72
RunType( 5 : 5 ) = 'R'
430
145
PlotType = 'rectilin '
431
-
END SELECT
432
-
433
#####
SELECT CASE ( RunType( 6 : 6 ) )
434
-
CASE ( '2' )
435
#####
WRITE( PRTFile, * ) 'N x 2D calculation (neglects horizontal refraction)'
436
-
CASE ( '3' )
437
#####
WRITE( PRTFile, * ) '3D calculation'
438
-
CASE DEFAULT
439
73
RunType( 6 : 6 ) = '2'
440
-
END SELECT
441
-
442
73
END SUBROUTINE ReadRunType
443
-
444
-
! **********************************************************************!
445
-
446
148*
SUBROUTINE TopBot( freq, AttenUnit, HS )
447
-
!! Handles top and bottom boundary conditions
448
-
449
-
REAL (KIND=8), INTENT( IN ) :: freq ! center / nominal frequency (wideband not supported)
450
-
CHARACTER (LEN=2), INTENT( IN ) :: AttenUnit
451
-
TYPE( HSInfo ), INTENT( INOUT ) :: HS
452
-
REAL (KIND=8) :: Mz, vr, alpha2_f ! values related to grain size
453
-
454
-
! Echo to PRTFile user's choice of boundary condition
455
-
456
72
SELECT CASE ( HS%BC )
457
-
CASE ( 'V' )
458
72*
WRITE( PRTFile, * ) ' VACUUM'
459
-
CASE ( 'R' )
460
#####
WRITE( PRTFile, * ) ' Perfectly RIGID'
461
-
CASE ( 'A' )
462
74*
WRITE( PRTFile, * ) ' ACOUSTO-ELASTIC half-space'
463
-
CASE ( 'G' )
464
#####
WRITE( PRTFile, * ) ' Grain size to define half-space'
465
-
CASE ( 'F' )
466
2*
WRITE( PRTFile, * ) ' FILE used for reflection loss'
467
-
CASE ( 'W' )
468
#####
WRITE( PRTFile, * ) ' Writing an IRC file'
469
-
CASE ( 'P' )
470
#####
WRITE( PRTFile, * ) ' reading PRECALCULATED IRC'
471
-
CASE DEFAULT
472
148*
CALL ERROUT( 'TopBot', 'Unknown boundary condition type' )
473
-
END SELECT
474
-
475
-
! ****** Read in BC parameters depending on particular choice ******
476
-
477
148
HS%cp = 0.0
478
148
HS%cs = 0.0
479
148
HS%rho = 0.0
480
-
481
74
SELECT CASE ( HS%BC )
482
-
CASE ( 'A' ) ! *** Half-space properties ***
483
74
zTemp = 0.0
484
74
READ( ENVFile, * ) zTemp, alphaR, betaR, rhoR, alphaI, betaI
485
-
WRITE( PRTFile, FMT = "( F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )" ) &
486
73
zTemp, alphaR, betaR, rhoR, alphaI, betaI
487
-
488
-
! dummy parameters for a layer with a general power law for attenuation
489
-
! these are not in play because the AttenUnit for this is not allowed yet
490
-
!freq0 = freq
491
73
betaPowerLaw = 1.0
492
73
ft = 1000.0
493
-
494
73
HS%cp = CRCI( zTemp, alphaR, alphaI, freq, freq, AttenUnit, betaPowerLaw, ft )
495
73
HS%cs = CRCI( zTemp, betaR, betaI, freq, freq, AttenUnit, betaPowerLaw, ft )
496
-
497
73*
HS%rho = rhoR
498
-
CASE ( 'G' ) ! *** Grain size (formulas from UW-APL HF Handbook)
499
-
500
-
! These formulas are from the UW-APL Handbook
501
-
! The code is taken from older Matlab and is unnecesarily verbose
502
-
! vr is the sound speed ratio
503
-
! rhor is the density ratio
504
#####
READ( ENVFile, * ) zTemp, Mz
505
#####
WRITE( PRTFile, FMT = "( F10.2, 3X, F10.2 )" ) zTemp, Mz
506
-
507
#####
IF ( Mz >= -1 .AND. Mz < 1 ) THEN
508
#####
vr = 0.002709 * Mz ** 2 - 0.056452 * Mz + 1.2778
509
#####
rhor = 0.007797 * Mz ** 2 - 0.17057 * Mz + 2.3139
510
#####
ELSE IF ( Mz >= 1 .AND. Mz < 5.3 ) THEN
511
#####
vr = -0.0014881 * Mz ** 3 + 0.0213937 * Mz ** 2 - 0.1382798 * Mz + 1.3425
512
#####
rhor = -0.0165406 * Mz ** 3 + 0.2290201 * Mz ** 2 - 1.1069031 * Mz + 3.0455
513
-
ELSE
514
#####
vr = -0.0024324 * Mz + 1.0019
515
#####
rhor = -0.0012973 * Mz + 1.1565
516
-
END IF
517
-
518
#####
IF ( Mz >= -1 .AND. Mz < 0 ) THEN
519
#####
alpha2_f = 0.4556
520
#####
ELSE IF ( Mz >= 0 .AND. Mz < 2.6 ) THEN
521
#####
alpha2_f = 0.4556 + 0.0245 * Mz
522
#####
ELSE IF( Mz >= 2.6 .AND. Mz < 4.5 ) THEN
523
#####
alpha2_f = 0.1978 + 0.1245 * Mz
524
#####
ELSE IF( Mz >= 4.5 .AND. Mz < 6.0 ) THEN
525
#####
alpha2_f = 8.0399 - 2.5228 * Mz + 0.20098 * Mz ** 2
526
#####
ELSE IF( Mz >= 6.0 .AND. Mz < 9.5 ) THEN
527
#####
alpha2_f = 0.9431 - 0.2041 * Mz + 0.0117 * Mz ** 2
528
-
ELSE
529
#####
alpha2_f = 0.0601
530
-
END IF
531
-
532
-
! AttenUnit = 'L' ! loss parameter
533
-
!!! following uses a reference sound speed of 1500 ???
534
-
!!! should be sound speed in the water, just above the sediment
535
-
! the term vr / 1000 converts vr to units of m per ms
536
#####
alphaR = vr * 1500.0
537
#####
alphaI = alpha2_f * ( vr / 1000 ) * 1500.0 * log( 10.0 ) / ( 40.0 * pi ) ! loss parameter Sect. IV., Eq. (4) of handbook
538
-
539
#####
HS%cp = CRCI( zTemp, alphaR, alphaI, freq, freq, 'L ', betaPowerLaw, ft )
540
#####
HS%cs = 0.0
541
#####
HS%rho = rhoR
542
-
WRITE( PRTFile, FMT = "( 'Converted sound speed =', 2F10.2, 3X, 'density = ', F10.2, 3X, 'loss parm = ', F10.4 )" ) &
543
#####
HS%cp, rhor, alphaI
544
-
545
-
CASE ( 'V', 'R', 'F', 'W', 'P' )
546
74*
CONTINUE
547
-
CASE DEFAULT
548
#####
WRITE( ERROR_UNIT, * ) 'TopBot: Unknown boundary condition type: ', HS%BC
549
148*
CALL ERROUT( 'TopBot', 'Unknown boundary condition type' )
550
-
END SELECT
551
-
552
147
END SUBROUTINE TopBot
553
-
554
-
! **********************************************************************!
555
-
556
73*
SUBROUTINE OpenOutputFiles( FileRoot, ThreeD )
557
-
!! Opens output files based on run type
558
-
559
-
! Write appropriate header information
560
-
561
-
USE SourceReceiverPositions
562
-
USE angleMod
563
-
USE bdryMod
564
-
USE RWSHDFile
565
-
566
-
LOGICAL, INTENT( IN ) :: ThreeD
567
-
CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
568
-
REAL :: atten
569
-
CHARACTER (LEN=10) :: PlotType
570
-
571
29
SELECT CASE ( Beam%RunType( 1 : 1 ) )
572
-
CASE ( 'R', 'E' ) ! Ray trace or Eigenrays
573
29*
OPEN ( FILE = TRIM( FileRoot ) // '.ray', UNIT = RAYFile, FORM = 'FORMATTED' )
574
29
WRITE( RAYFile, * ) '''', Title( 1 : 50 ), ''''
575
29
WRITE( RAYFile, * ) freq
576
29
WRITE( RAYFile, * ) Pos%NSx, Pos%NSy, Pos%NSz
577
29
WRITE( RAYFile, * ) Angles%Nalpha, Angles%Nbeta
578
29
WRITE( RAYFile, * ) Bdry%Top%HS%Depth
579
29
WRITE( RAYFile, * ) Bdry%Bot%HS%Depth
580
-
581
58
IF ( ThreeD ) THEN
582
#####
WRITE( RAYFile, * ) '''xyz'''
583
-
ELSE
584
29
WRITE( RAYFile, * ) '''rz'''
585
-
END IF
586
-
587
-
CASE ( 'A' ) ! arrival file in ascii format
588
26*
OPEN ( FILE = TRIM( FileRoot ) // '.arr', UNIT = ARRFile, FORM = 'FORMATTED' )
589
-
590
26
IF ( ThreeD ) THEN
591
#####
WRITE( ARRFile, * ) '''3D'''
592
-
ELSE
593
26
WRITE( ARRFile, * ) '''2D'''
594
-
END IF
595
-
596
26
WRITE( ARRFile, * ) freq
597
-
598
-
! write source locations
599
-
600
26
IF ( ThreeD ) THEN
601
#####
WRITE( ARRFile, * ) Pos%NSx, Pos%Sx( 1 : Pos%NSx )
602
#####
WRITE( ARRFile, * ) Pos%NSy, Pos%Sy( 1 : Pos%NSy )
603
#####
WRITE( ARRFile, * ) Pos%NSz, Pos%Sz( 1 : Pos%NSz )
604
-
ELSE
605
26*
WRITE( ARRFile, * ) Pos%NSz, Pos%Sz( 1 : Pos%NSz )
606
-
END IF
607
-
608
-
! write receiver locations
609
-
610
26*
WRITE( ARRFile, * ) Pos%NRz, Pos%Rz( 1 : Pos%NRz )
611
26*
WRITE( ARRFile, * ) Pos%NRr, Pos%Rr( 1 : Pos%NRr )
612
26*
IF ( ThreeD ) THEN
613
#####
WRITE( ARRFile, * ) Pos%Ntheta, Pos%theta( 1 : Pos%Ntheta )
614
-
END IF
615
-
616
-
CASE ( 'a' ) ! arrival file in binary format
617
#####
OPEN ( FILE = TRIM( FileRoot ) // '.arr', UNIT = ARRFile, FORM = 'UNFORMATTED' )
618
-
619
#####
IF ( ThreeD ) THEN
620
#####
WRITE( ARRFile ) '''3D'''
621
-
ELSE
622
#####
WRITE( ARRFile ) '''2D'''
623
-
END IF
624
-
625
#####
WRITE( ARRFile ) SNGL( freq )
626
-
627
-
! write source locations
628
-
629
#####
IF ( ThreeD ) THEN
630
#####
WRITE( ARRFile ) Pos%NSx, Pos%Sx( 1 : Pos%NSx )
631
#####
WRITE( ARRFile ) Pos%NSy, Pos%Sy( 1 : Pos%NSy )
632
#####
WRITE( ARRFile ) Pos%NSz, Pos%Sz( 1 : Pos%NSz )
633
-
ELSE
634
#####
WRITE( ARRFile ) Pos%NSz, Pos%Sz( 1 : Pos%NSz )
635
-
END IF
636
-
637
-
! write receiver locations
638
-
639
#####
WRITE( ARRFile ) Pos%NRz, Pos%Rz( 1 : Pos%NRz )
640
#####
WRITE( ARRFile ) Pos%NRr, Pos%Rr( 1 : Pos%NRr )
641
#####
IF ( ThreeD ) THEN
642
#####
WRITE( ARRFile ) Pos%Ntheta, Pos%theta( 1 : Pos%Ntheta )
643
-
END IF
644
-
645
-
CASE DEFAULT
646
18
atten = 0.0
647
-
648
-
! following to set PlotType has already been done in READIN if that was used for input
649
18
SELECT CASE ( Beam%RunType( 5 : 5 ) )
650
-
CASE ( 'R' )
651
18*
PlotType = 'rectilin '
652
-
CASE ( 'I' )
653
#####
PlotType = 'irregular '
654
-
CASE DEFAULT
655
18*
PlotType = 'rectilin '
656
-
END SELECT
657
-
658
91*
CALL WriteHeader( TRIM( FileRoot ) // '.shd', Title, REAL( freq ), atten, PlotType )
659
-
END SELECT
660
-
661
73
END SUBROUTINE OpenOutputFiles
662
-
663
-
END MODULE ReadEnvironmentBell

Legend:

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