Coverage Report: ReadEnvironmentBell.f90

Generated from GCOV analysis of Fortran source code

56.4%
Lines Executed
303 total lines
51.5%
Branches Executed
169 total branches
100.0%
Calls Executed
262 total calls
0
-
Source:ReadEnvironmentBell.f90
0
-
Graph:ReadEnvironmentBell.gcno
0
-
Data:ReadEnvironmentBell.gcda
0
-
Runs:73
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
73*
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
73
WRITE( PRTFile, * ) 'BELLHOP/BELLHOP3D'
37
73
WRITE( PRTFile, * )
38
-
39
-
! Open the environmental file
40
73*
OPEN( UNIT = ENVFile, FILE = TRIM( FileRoot ) // '.env', STATUS = 'OLD', IOSTAT = iostat, ACTION = 'READ' )
41
73
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
73
IF ( ThreeD ) THEN
48
#####
Title( 1 : 11 ) = 'BELLHOP3D- '
49
#####
READ( ENVFile, * ) Title( 12 : 80 )
50
-
ELSE
51
73
Title( 1 : 9 ) = 'BELLHOP- '
52
73
READ( ENVFile, * ) Title( 10 : 80 )
53
-
END IF
54
-
55
73
WRITE( PRTFile, * ) Title
56
-
57
73
READ( ENVFile, * ) freq
58
73
WRITE( PRTFile, '('' frequency = '', G11.4, '' Hz'', / )' ) freq
59
-
60
73
READ( ENVFile, * ) NMedia
61
73
WRITE( PRTFile, * ) 'Dummy parameter NMedia = ', NMedia
62
73
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
73
CALL ReadTopOpt( Bdry%Top%HS%Opt, Bdry%Top%HS%BC, AttenUnit, FileRoot )
66
-
67
-
! *** Top BC ***
68
-
69
73
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
73
CALL TopBot( freq, AttenUnit, Bdry%Top%HS )
75
-
76
-
! ****** Read in ocean SSP data ******
77
-
78
72
READ( ENVFile, * ) NPts, Sigma, Bdry%Bot%HS%Depth
79
72
WRITE( PRTFile, * )
80
72
WRITE( PRTFile, FMT = "( ' Depth = ', F10.2, ' m' )" ) Bdry%Bot%HS%Depth
81
-
82
72
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
216
x = [ 0.0D0, Bdry%Bot%HS%Depth ] ! tells SSP Depth to read to
90
72
t = [ 0.0, 1.0 ]
91
72
CALL EvaluateSSP( x, t, c, cimag, gradc, crr, crz, czz, rho, freq, 'INI' )
92
-
END IF
93
-
94
71
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
71
Bdry%Bot%HS%Opt = ' ' ! initialize to blanks
99
71
READ( ENVFile, * ) Bdry%Bot%HS%Opt, Sigma
100
71
WRITE( PRTFile, * )
101
71
WRITE( PRTFile, FMT = "(33X, '( RMS roughness = ', G10.3, ' )' )" ) Sigma
102
-
103
11
SELECT CASE ( Bdry%Bot%HS%Opt( 2 : 2 ) )
104
-
CASE ( '~', '*' )
105
71*
WRITE( PRTFile, * ) ' Bathymetry file selected'
106
-
CASE( '-', '_', ' ' )
107
-
CASE DEFAULT
108
71*
CALL ERROUT( 'READIN', 'Unknown bottom option letter in second position' )
109
-
END SELECT
110
-
111
71
Bdry%Bot%HS%BC = Bdry%Bot%HS%Opt( 1 : 1 )
112
71
CALL TopBot( freq, AttenUnit, Bdry%Bot%HS )
113
-
114
-
! *** source and receiver locations ***
115
-
116
71
CALL ReadSxSy( ThreeD ) ! Read source/receiver x-y coordinates
117
-
118
71
ZMin = SNGL( Bdry%Top%HS%Depth )
119
71
ZMax = SNGL( Bdry%Bot%HS%Depth )
120
-
! CALL ReadSzRz( ZMin + 100 * SPACING( ZMin ), ZMax - 100 * SPACING( ZMax ) ) ! not sure why I had this
121
71
CALL ReadSzRz( ZMin, ZMax )
122
71
CALL ReadRcvrRanges
123
71*
IF ( ThreeD ) CALL ReadRcvrBearings
124
71
CALL ReadfreqVec( freq, Bdry%Top%HS%Opt( 6:6 ) )
125
71
CALL ReadRunType( Beam%RunType, PlotType )
126
-
127
71
Depth = Zmax - Zmin ! water depth
128
71
CALL ReadRayElevationAngles( freq, Depth, Bdry%Top%HS%Opt, Beam%RunType )
129
71*
IF ( ThreeD ) CALL ReadRayBearingAngles( freq, Bdry%Top%HS%Opt, Beam%RunType )
130
-
131
71
WRITE( PRTFile, * )
132
71
WRITE( PRTFile, * ) '__________________________________________________________________________'
133
71
WRITE( PRTFile, * )
134
-
135
-
! Limits for tracing beams
136
-
137
71
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
71
READ( ENVFile, * ) Beam%deltas, Beam%Box%z, Beam%Box%r
153
-
154
71
WRITE( PRTFile, * )
155
71
WRITE( PRTFile, fmt = '( '' Step length, deltas = '', G11.4, '' m'' )' ) Beam%deltas
156
71
WRITE( PRTFile, * )
157
71
WRITE( PRTFile, fmt = '( '' Maximum ray depth, Box%z = '', G11.4, '' m'' )' ) Beam%Box%z
158
71
WRITE( PRTFile, fmt = '( '' Maximum ray range, Box%r = '', G11.4, ''km'' )' ) Beam%Box%r
159
-
160
71
Beam%Box%r = 1000.0 * Beam%Box%r ! convert km to m
161
-
END IF
162
-
163
-
! *** Beam characteristics ***
164
71
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
71
WRITE( PRTFile, * ) 'No beam shift in effect'
171
-
END SELECT
172
-
173
71
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
60
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
#####
READ( ENVFile, * ) Beam%Type( 2 : 3 ), Beam%epsMultiplier, Beam%rLoop
203
#####
WRITE( PRTFile, * )
204
#####
WRITE( PRTFile, * )
205
#####
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
#####
WRITE( PRTFile, * ) 'Standard curvature condition'
214
-
CASE DEFAULT
215
#####
CALL ERROUT( 'READIN', 'Unknown curvature condition' )
216
-
END SELECT
217
-
218
#####
WRITE( PRTFile, * ) 'Epsilon multiplier', Beam%epsMultiplier
219
#####
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
#####
Beam%Nimage = 1
225
#####
Beam%iBeamWindow = 4
226
#####
Beam%Component = 'P'
227
#####
READ( ENVFile, * ) Beam%Nimage, Beam%iBeamWindow, Beam%Component
228
#####
WRITE( PRTFile, * )
229
#####
WRITE( PRTFile, * ) 'Number of images, Nimage = ', Beam%Nimage
230
#####
WRITE( PRTFile, * ) 'Beam windowing parameter = ', Beam%iBeamWindow
231
#####
WRITE( PRTFile, * ) 'Component = ', Beam%Component
232
-
CASE DEFAULT
233
60*
CALL ERROUT( 'READIN', 'Unknown beam type (second letter of run type)' )
234
-
END SELECT
235
-
END IF
236
-
237
71
WRITE( PRTFile, * )
238
71
CLOSE( ENVFile )
239
-
240
71
END SUBROUTINE ReadEnvironment
241
-
242
-
! **********************************************************************!
243
-
244
73*
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
73
TopOpt = ' ' ! initialize to blanks
253
73
READ( ENVFile, * ) TopOpt
254
73
WRITE( PRTFile, * )
255
-
256
73
SSP%Type = TopOpt( 1 : 1 )
257
73
BC = TopOpt( 2 : 2 )
258
73
AttenUnit = TopOpt( 3 : 4 )
259
73
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
64
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
73*
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
5*
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
73*
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
71*
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
73*
CALL ERROUT( 'READIN', 'Unknown top option letter in fourth position' )
337
-
END SELECT
338
-
339
6
SELECT CASE ( TopOpt( 5 : 5 ) )
340
-
CASE ( '~', '*' )
341
73*
WRITE( PRTFile, * ) ' Altimetry file selected'
342
-
CASE ( '-', '_', ' ' )
343
-
CASE DEFAULT
344
73*
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
73*
WRITE( PRTFile, * ) ' Development options enabled'
350
-
CASE ( ' ' )
351
-
CASE DEFAULT
352
73*
CALL ERROUT( 'READIN', 'Unknown top option letter in sixth position' )
353
-
END SELECT
354
-
355
73
END SUBROUTINE ReadTopOpt
356
-
357
-
! **********************************************************************!
358
-
359
71*
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
71
READ( ENVFile, * ) RunType
370
71
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
16
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
71*
CALL ERROUT( 'READIN', 'Unknown RunType selected' )
389
-
END SELECT
390
-
391
#####
SELECT CASE ( RunType( 2 : 2 ) )
392
-
CASE ( 'C' )
393
#####
WRITE( PRTFile, * ) 'Cartesian beams'
394
-
CASE ( 'R' )
395
#####
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
71
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
68
RunType( 4 : 4 ) = 'R'
416
71
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
70
WRITE( PRTFile, * ) 'Rectilinear receiver grid: Receivers at Rr( : ) x Rz( : )'
429
70
RunType( 5 : 5 ) = 'R'
430
141
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
71
RunType( 6 : 6 ) = '2'
440
-
END SELECT
441
-
442
71
END SUBROUTINE ReadRunType
443
-
444
-
! **********************************************************************!
445
-
446
144*
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
70
SELECT CASE ( HS%BC )
457
-
CASE ( 'V' )
458
70*
WRITE( PRTFile, * ) ' VACUUM'
459
-
CASE ( 'R' )
460
#####
WRITE( PRTFile, * ) ' Perfectly RIGID'
461
-
CASE ( 'A' )
462
72*
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
144*
CALL ERROUT( 'TopBot', 'Unknown boundary condition type' )
473
-
END SELECT
474
-
475
-
! ****** Read in BC parameters depending on particular choice ******
476
-
477
144
HS%cp = 0.0
478
144
HS%cs = 0.0
479
144
HS%rho = 0.0
480
-
481
72
SELECT CASE ( HS%BC )
482
-
CASE ( 'A' ) ! *** Half-space properties ***
483
72
zTemp = 0.0
484
72
READ( ENVFile, * ) zTemp, alphaR, betaR, rhoR, alphaI, betaI
485
-
WRITE( PRTFile, FMT = "( F10.2, 3X, 2F10.2, 3X, F6.2, 3X, 2F10.4 )" ) &
486
71
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
71
betaPowerLaw = 1.0
492
71
ft = 1000.0
493
-
494
71
HS%cp = CRCI( zTemp, alphaR, alphaI, freq, freq, AttenUnit, betaPowerLaw, ft )
495
71
HS%cs = CRCI( zTemp, betaR, betaI, freq, freq, AttenUnit, betaPowerLaw, ft )
496
-
497
71*
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
72*
CONTINUE
547
-
CASE DEFAULT
548
#####
WRITE( ERROR_UNIT, * ) 'TopBot: Unknown boundary condition type: ', HS%BC
549
144*
CALL ERROUT( 'TopBot', 'Unknown boundary condition type' )
550
-
END SELECT
551
-
552
143
END SUBROUTINE TopBot
553
-
554
-
! **********************************************************************!
555
-
556
71*
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
16
atten = 0.0
647
-
648
-
! following to set PlotType has already been done in READIN if that was used for input
649
16
SELECT CASE ( Beam%RunType( 5 : 5 ) )
650
-
CASE ( 'R' )
651
16*
PlotType = 'rectilin '
652
-
CASE ( 'I' )
653
#####
PlotType = 'irregular '
654
-
CASE DEFAULT
655
16*
PlotType = 'rectilin '
656
-
END SELECT
657
-
658
87*
CALL WriteHeader( TRIM( FileRoot ) // '.shd', Title, REAL( freq ), atten, PlotType )
659
-
END SELECT
660
-
661
71
END SUBROUTINE OpenOutputFiles
662
-
663
-
END MODULE ReadEnvironmentBell

Legend:

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