Coverage Report: ReadEnvironmentBell.f90

Generated from GCOV analysis of Fortran source code

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

Legend:

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