Coverage Report: ArrMod.f90

Generated from GCOV analysis of Fortran source code

27.1%
Lines Executed
155 total lines
52.0%
Branches Executed
404 total branches
100.0%
Calls Executed
16 total calls
0
-
Source:ArrMod.f90
0
-
Graph:ArrMod.gcno
0
-
Data:ArrMod.gcda
0
-
Runs:29
1
-
!! Acoustic arrivals module for computing and storing arrival information
2
-
3
-
MODULE ArrMod
4
-
!! Management of acoustic arrival data including storage, sorting, and output formatting
5
-
6
-
USE MathConstants
7
-
USE BellhopMod
8
-
9
-
IMPLICIT NONE
10
-
PUBLIC
11
-
12
-
! Variables for arrival information
13
-
REAL, PARAMETER :: PhaseTol = 0.05 ! arrivals with essentially the same phase are grouped into one
14
-
INTEGER :: MaxNArr
15
-
INTEGER, ALLOCATABLE :: NArr( :, : ), NArr3D( :, :, : )
16
-
REAL (KIND=4) :: factor = 1.0
17
-
18
-
TYPE Arrival
19
-
INTEGER :: NTopBnc, NBotBnc
20
-
REAL :: SrcDeclAngle, SrcAzimAngle, RcvrDeclAngle, RcvrAzimAngle, A, Phase
21
-
COMPLEX :: delay
22
-
END TYPE Arrival
23
-
24
-
TYPE(Arrival), ALLOCATABLE :: Arr( :, :, : ), Arr3D( :, :, :, : )
25
-
26
-
CONTAINS
27
-
28
962
SUBROUTINE AddArr( omega, id, ir, Amp, Phase, delay, SrcDeclAngle, RcvrDeclAngle, NumTopBnc, NumBotBnc )
29
-
!! Adds an arrival to the arrival data structure
30
-
31
-
! Adds the amplitude and delay for an ARRival into a matrix of same.
32
-
! Extra logic included to keep only the strongest arrivals.
33
-
34
-
INTEGER, INTENT( IN ) :: NumTopBnc, NumBotBnc, id, ir
35
-
REAL ( KIND = 8 ), INTENT( IN ) :: omega, Amp, Phase, SrcDeclAngle, RcvrDeclAngle
36
-
COMPLEX ( KIND = 8 ), INTENT( IN ) :: delay
37
-
LOGICAL :: NewRay
38
-
INTEGER :: iArr( 1 ), Nt
39
-
REAL :: AmpTot, w1, w2
40
-
41
962*
Nt = NArr( id, ir ) ! # of arrivals
42
962
NewRay = .TRUE.
43
-
44
-
! Is this the second step of a pair (on the same ray)?
45
-
! If so, we want to combine the arrivals to conserve space.
46
-
! (test this by seeing if the arrival time is close to the previous one)
47
-
! (also need that the phase is about the same to make sure surface and direct paths are not joined)
48
-
! LP: BUG: This only checks the last arrival, whereas the first step of the
49
-
! pair could have been placed in previous slots. See readme.
50
-
51
962
IF ( Nt >= 1 ) THEN
52
1904*
IF ( omega * ABS( delay - Arr( id, ir, Nt )%delay ) < PhaseTol .AND. &
53
1409*
ABS( Arr( id, ir, Nt )%phase - Phase ) < PhaseTol ) NewRay = .FALSE.
54
-
END IF
55
-
56
962
IF ( NewRay ) THEN
57
505
IF ( Nt >= MaxNArr ) THEN ! space not available to add an arrival?
58
#####
iArr = MINLOC( Arr( id, ir, : )%A ) ! replace weakest arrival
59
#####
IF ( Amp > Arr( id, ir, iArr( 1 ) )%A ) THEN
60
#####
Arr( id, ir, iArr( 1 ) )%A = SNGL( Amp ) ! amplitude
61
#####
Arr( id, ir, iArr( 1 ) )%Phase = SNGL( Phase ) ! phase
62
#####
Arr( id, ir, iArr( 1 ) )%delay = CMPLX( delay ) ! delay time
63
#####
Arr( id, ir, iArr( 1 ) )%SrcDeclAngle = SNGL( SrcDeclAngle ) ! launch angle from source
64
#####
Arr( id, ir, iArr( 1 ) )%RcvrDeclAngle = SNGL( RcvrDeclAngle ) ! angle ray reaches receiver
65
#####
Arr( id, ir, iArr( 1 ) )%NTopBnc = NumTopBnc ! Number of top bounces
66
#####
Arr( id, ir, iArr( 1 ) )%NBotBnc = NumBotBnc ! " bottom
67
-
ENDIF
68
-
ELSE
69
505*
NArr( id, ir ) = Nt + 1 ! # of arrivals
70
505*
Arr( id, ir, Nt + 1 )%A = SNGL( Amp ) ! amplitude
71
505*
Arr( id, ir, Nt + 1 )%Phase = SNGL( Phase ) ! phase
72
505*
Arr( id, ir, Nt + 1 )%delay = CMPLX( delay ) ! delay time
73
505*
Arr( id, ir, Nt + 1 )%SrcDeclAngle = SNGL( SrcDeclAngle ) ! launch angle from source
74
505*
Arr( id, ir, Nt + 1 )%RcvrDeclAngle = SNGL( RcvrDeclAngle ) ! angle ray reaches receiver
75
505*
Arr( id, ir, Nt + 1 )%NTopBnc = NumTopBnc ! Number of top bounces
76
505*
Arr( id, ir, Nt + 1 )%NBotBnc = NumBotBnc ! " bottom
77
-
ENDIF
78
-
ELSE ! not a new ray
79
-
!PhaseArr( id, ir, Nt ) = PhaseArr( id, ir, Nt )
80
-
81
-
! calculate weightings of old ray information vs. new, based on amplitude of the arrival
82
457*
AmpTot = Arr( id, ir, Nt )%A + SNGL( Amp )
83
457*
w1 = Arr( id, ir, Nt )%A / AmpTot
84
457
w2 = REAL( Amp ) / AmpTot
85
-
86
457*
Arr( id, ir, Nt )%delay = w1 * Arr( id, ir, Nt )%delay + w2 * CMPLX( delay ) ! weighted sum
87
457*
Arr( id, ir, Nt )%A = AmpTot
88
457*
Arr( id, ir, Nt )%SrcDeclAngle = w1 * Arr( id, ir, Nt )%SrcDeclAngle + w2 * SNGL( SrcDeclAngle )
89
457*
Arr( id, ir, Nt )%RcvrDeclAngle = w1 * Arr( id, ir, Nt )%RcvrDeclAngle + w2 * SNGL( RcvrDeclAngle )
90
-
ENDIF
91
-
92
962
RETURN
93
-
END SUBROUTINE AddArr
94
-
95
-
! **********************************************************************!
96
-
97
8*
SUBROUTINE WriteArrivalsASCII( r, Nrd, Nr, SourceType )
98
-
!! Writes arrival data in ASCII format
99
-
100
-
! Writes the arrival data (Amplitude, delay for each eigenray)
101
-
! ASCII output file
102
-
103
-
INTEGER, INTENT( IN ) :: Nrd, Nr
104
-
REAL, INTENT( IN ) :: r( Nr )
105
-
CHARACTER (LEN=1), INTENT( IN ) :: SourceType
106
-
INTEGER :: ir, id, iArr
107
-
108
34*
WRITE( ARRFile, * ) MAXVAL( NArr( 1 : Nrd, 1 : Nr ) )
109
-
110
16*
DO id = 1, Nrd
111
29*
DO ir = 1, Nr
112
13
IF ( SourceType == 'X' ) THEN ! line source
113
#####
factor = 4.0 * SQRT( pi )
114
-
ELSE ! point source
115
13*
IF ( r ( ir ) == 0 ) THEN
116
#####
factor = 1e5 ! avoid /0 at origin
117
-
ELSE
118
13*
factor = 1. / SQRT( ABS( r( ir ) ) ) ! cyl. spreading [CF] ABS() used to eliminate NaN for -ve range values
119
-
END IF
120
-
END IF
121
-
122
13*
WRITE( ARRFile, * ) NArr( id, ir )
123
526*
DO iArr = 1, NArr( id, ir )
124
-
! You can compress the output file a lot by putting in an explicit format statement here ...
125
-
! However, you'll need to make sure you keep adequate precision
126
-
WRITE( ARRFile, * ) &
127
505*
factor * Arr( id, ir, iArr )%A, &
128
505*
SNGL( RadDeg ) * Arr( id, ir, iArr )%Phase, &
129
505*
REAL( Arr( id, ir, iArr )%delay ), &
130
505*
AIMAG( Arr( id, ir, iArr )%delay ), &
131
505*
Arr( id, ir, iArr )%SrcDeclAngle, &
132
505*
Arr( id, ir, iArr )%RcvrDeclAngle, &
133
505*
Arr( id, ir, iArr )%NTopBnc, &
134
1023
Arr( id, ir, iArr )%NBotBnc
135
-
END DO ! next arrival
136
-
END DO ! next receiver depth
137
-
END DO ! next receiver range
138
-
139
8
RETURN
140
-
END SUBROUTINE WriteArrivalsASCII
141
-
142
-
! **********************************************************************!
143
-
144
-
145
#####
SUBROUTINE WriteArrivalsBinary( r, Nrd, Nr, SourceType )
146
-
!! Writes arrival data in binary format
147
-
148
-
! Writes the arrival data (amplitude, delay for each eigenray)
149
-
! Binary output file
150
-
151
-
INTEGER, INTENT( IN ) :: Nrd, Nr
152
-
REAL, INTENT( IN ) :: r( Nr )
153
-
CHARACTER (LEN=1), INTENT( IN ) :: SourceType
154
-
INTEGER :: ir, id, iArr
155
-
156
#####
WRITE( ARRFile ) MAXVAL( NArr( 1 : Nrd, 1 : Nr ) )
157
-
158
#####
DO id = 1, Nrd
159
#####
DO ir = 1, Nr
160
#####
IF ( SourceType == 'X' ) THEN ! line source
161
#####
factor = 4.0 * SQRT( pi )
162
-
ELSE ! point source
163
#####
IF ( r ( ir ) == 0 ) THEN
164
#####
factor = 1e5 ! avoid /0 at origin
165
-
ELSE
166
#####
factor = 1. / SQRT( r( ir ) ) ! cyl. spreading
167
-
END IF
168
-
END IF
169
-
170
#####
WRITE( ARRFile ) NArr( id, ir )
171
-
172
#####
DO iArr = 1, NArr( id, ir )
173
-
! integers written out as reals below for fast reading in Matlab
174
-
WRITE( ARRFile ) &
175
#####
factor * Arr( id, ir, iArr )%A, &
176
#####
SNGL( RadDeg * Arr( id, ir, iArr )%Phase ), &
177
#####
Arr( id, ir, iArr )%delay, &
178
#####
Arr( id, ir, iArr )%SrcDeclAngle, &
179
#####
Arr( id, ir, iArr )%RcvrDeclAngle, &
180
#####
REAL( Arr( id, ir, iArr )%NTopBnc ), &
181
#####
REAL( Arr( id, ir, iArr )%NBotBnc )
182
-
183
-
END DO ! next arrival
184
-
END DO ! next receiver depth
185
-
END DO ! next receiver range
186
-
187
#####
RETURN
188
-
END SUBROUTINE WriteArrivalsBinary
189
-
190
-
! **********************************************************************!
191
-
192
#####
SUBROUTINE AddArr3D( omega, itheta, id, ir, Amp, Phase, delay, SrcDeclAngle, &
193
-
SrcAzimAngle, RcvrDeclAngle, RcvrAzimAngle, NumTopBnc, NumBotBnc )
194
-
!! Adds the amplitude and delay for an ARRival into a matrix of same.
195
-
196
-
! Extra logic included to keep only the strongest arrivals.
197
-
198
-
INTEGER, INTENT( IN ) :: itheta, id, ir
199
-
INTEGER, INTENT( IN ) :: NumTopBnc, NumBotBnc
200
-
REAL ( KIND = 8 ), INTENT( IN ) :: omega, Amp, Phase, SrcDeclAngle, SrcAzimAngle, RcvrDeclAngle, RcvrAzimAngle
201
-
202
-
COMPLEX ( KIND = 8 ), INTENT( IN ) :: delay
203
-
LOGICAL :: NewRay
204
-
INTEGER :: iArr( 1 ), Nt
205
-
REAL :: AmpTot, w1, w2
206
-
207
#####
Nt = NArr3D( itheta, id, ir ) ! # of arrivals
208
#####
NewRay = .TRUE.
209
-
210
-
! Is this the second bracketing ray of a pair?
211
-
! If so, we want to combine the arrivals to conserve space.
212
-
! (test this by seeing if the arrival time is close to the previous one)
213
-
! (also need that the phase is about the same to make sure surface and direct paths are not joined)
214
-
! LP: BUG: This only checks the last arrival, whereas the first step of the
215
-
! pair could have been placed in previous slots. See readme.
216
-
217
#####
IF ( Nt >= 1 ) THEN
218
#####
IF ( omega * ABS( delay - Arr3D( itheta, id, ir, Nt )%delay ) < PhaseTol .AND. &
219
#####
ABS( Arr3D( itheta, id, ir, Nt )%phase - Phase ) < PhaseTol ) NewRay = .FALSE.
220
-
END IF
221
-
222
#####
IF ( NewRay ) THEN
223
#####
IF ( Nt >= MaxNArr ) THEN ! space available to add an arrival?
224
#####
iArr = MINLOC( Arr3D( itheta, id, ir, : )%A ) ! no: replace weakest arrival
225
#####
IF ( Amp > Arr3D( itheta, id, ir, iArr( 1 ) )%A ) THEN
226
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%A = SNGL( Amp ) ! amplitude
227
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%Phase = SNGL( Phase ) ! phase
228
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%delay = CMPLX( delay ) ! delay time
229
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%SrcDeclAngle = SNGL( SrcDeclAngle ) ! angle
230
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%SrcAzimAngle = SNGL( SrcAzimAngle ) ! angle
231
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%RcvrDeclAngle = SNGL( RcvrDeclAngle ) ! angle
232
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%RcvrAzimAngle = SNGL( RcvrAzimAngle ) ! angle
233
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%NTopBnc = NumTopBnc ! Number of top bounces
234
#####
Arr3D( itheta, id, ir, iArr( 1 ) )%NBotBnc = NumBotBnc ! " bottom
235
-
ENDIF
236
-
ELSE
237
#####
NArr3D( itheta, id, ir ) = Nt + 1 ! # of arrivals
238
#####
Arr3D( itheta, id, ir, Nt + 1 )%A = SNGL( Amp ) ! amplitude
239
#####
Arr3D( itheta, id, ir, Nt + 1 )%Phase = SNGL( Phase ) ! phase
240
#####
Arr3D( itheta, id, ir, Nt + 1 )%delay = CMPLX( delay ) ! delay time
241
#####
Arr3D( itheta, id, ir, Nt + 1 )%SrcDeclAngle = SNGL( SrcDeclAngle ) ! angle
242
#####
Arr3D( itheta, id, ir, Nt + 1 )%SrcAzimAngle = SNGL( SrcAzimAngle ) ! angle
243
#####
Arr3D( itheta, id, ir, Nt + 1 )%RcvrDeclAngle = SNGL( RcvrDeclAngle ) ! angle
244
#####
Arr3D( itheta, id, ir, Nt + 1 )%RcvrAzimAngle = SNGL( RcvrAzimAngle ) ! angle
245
#####
Arr3D( itheta, id, ir, Nt + 1 )%NTopBnc = NumTopBnc ! Number of top bounces
246
#####
Arr3D( itheta, id, ir, Nt + 1 )%NBotBnc = NumBotBnc ! " bottom
247
-
ENDIF
248
-
ELSE ! not a new ray
249
-
!PhaseArr( id, ir, Nt ) = PhaseArr( id, ir, Nt )
250
-
! calculate weightings of old ray information vs. new, based on amplitude of the arrival
251
#####
AmpTot = Arr3D( itheta, id, ir, Nt )%A + SNGL( Amp )
252
#####
w1 = Arr3D( itheta, id, ir, Nt )%A / AmpTot
253
#####
w2 = REAL( Amp ) / AmpTot
254
-
255
#####
Arr3D( itheta, id, ir, Nt )%delay = w1 * Arr3D( itheta, id, ir, Nt )%delay + w2 * CMPLX( delay ) ! weighted sum
256
#####
Arr3D( itheta, id, ir, Nt )%A = AmpTot
257
#####
Arr3D( itheta, id, ir, Nt )%SrcDeclAngle = w1 * Arr3D( itheta, id, ir, Nt )%SrcDeclAngle + w2 * SNGL( SrcDeclAngle )
258
#####
Arr3D( itheta, id, ir, Nt )%SrcAzimAngle = w1 * Arr3D( itheta, id, ir, Nt )%SrcAzimAngle + w2 * SNGL( SrcAzimAngle )
259
#####
Arr3D( itheta, id, ir, Nt )%RcvrDeclAngle = w1 * Arr3D( itheta, id, ir, Nt )%RcvrDeclAngle + w2 * SNGL( RcvrDeclAngle )
260
#####
Arr3D( itheta, id, ir, Nt )%RcvrAzimAngle = w1 * Arr3D( itheta, id, ir, Nt )%RcvrAzimAngle + w2 * SNGL( RcvrAzimAngle )
261
-
ENDIF
262
-
263
#####
RETURN
264
-
END SUBROUTINE AddArr3D
265
-
266
-
! **********************************************************************!
267
-
268
#####
SUBROUTINE WriteArrivalsASCII3D( r, Ntheta, Nrd, Nr )
269
-
!! Writes the arrival data (Amplitude, delay for each eigenray); ASCII output file
270
-
271
-
INTEGER, INTENT( IN ) :: Ntheta, Nrd, Nr
272
-
REAL, INTENT( IN ) :: r( Nr )
273
-
INTEGER :: itheta, ir, id, iArr
274
-
275
#####
WRITE( ARRFile, * ) MAXVAL( NArr3D( 1 : Ntheta, 1 : Nrd, 1 : Nr ) )
276
-
277
#####
DO itheta = 1, Ntheta
278
#####
DO id = 1, Nrd
279
#####
DO ir = 1, Nr
280
#####
IF ( Beam%RunType( 6 : 6 ) == '2' ) THEN ! 2D run?
281
#####
IF ( r ( ir ) == 0 ) THEN
282
#####
factor = 1e5 ! avoid /0 at origin
283
-
ELSE
284
#####
factor = 1. / SQRT( r( ir ) ) ! cyl. spreading
285
-
END IF
286
-
END IF
287
-
288
#####
WRITE( ARRFile, * ) NArr3D( itheta, id, ir )
289
-
290
#####
DO iArr = 1, NArr3D( itheta, id, ir )
291
-
! you can compress the output file a lot by putting in an explicit format statement here ...
292
-
! However, you'll need to make sure you keep adequate precision
293
-
WRITE( ARRFile, * ) &
294
#####
factor * Arr3D( itheta, id, ir, iArr )%A, &
295
#####
RadDeg * Arr3D( itheta, id, ir, iArr )%Phase, &
296
#####
REAL( Arr3D( itheta, id, ir, iArr )%delay ), &
297
#####
AIMAG( Arr3D( itheta, id, ir, iArr )%delay ), &
298
#####
Arr3D( itheta, id, ir, iArr )%SrcDeclAngle, &
299
#####
Arr3D( itheta, id, ir, iArr )%SrcAzimAngle, &
300
#####
Arr3D( itheta, id, ir, iArr )%RcvrDeclAngle, &
301
#####
Arr3D( itheta, id, ir, iArr )%RcvrAzimAngle, &
302
#####
Arr3D( itheta, id, ir, iArr )%NTopBnc, &
303
#####
Arr3D( itheta, id, ir, iArr )%NBotBnc
304
-
END DO ! next arrival
305
-
END DO ! next receiver depth
306
-
END DO ! next receiver range
307
-
END DO ! next receiver angle
308
-
309
#####
RETURN
310
-
END SUBROUTINE WriteArrivalsASCII3D
311
-
312
-
! **********************************************************************!
313
-
314
#####
SUBROUTINE WriteArrivalsBinary3D( r, Ntheta, Nrd, Nr )
315
-
!! Writes the arrival data (amplitude, delay for each eigenray); Binary output file
316
-
317
-
INTEGER, INTENT( IN ) :: Ntheta, Nrd, Nr
318
-
REAL, INTENT( IN ) :: r( Nr )
319
-
INTEGER :: itheta, ir, id, iArr
320
-
321
#####
WRITE( ARRFile ) MAXVAL( NArr3D( 1 : Ntheta, 1 : Nrd, 1 : Nr ) )
322
-
323
#####
DO itheta = 1, Ntheta
324
#####
DO id = 1, Nrd
325
#####
DO ir = 1, Nr
326
#####
IF ( Beam%RunType( 6 : 6 ) == '2' ) THEN ! 2D run?
327
#####
IF ( r ( ir ) == 0 ) THEN
328
#####
factor = 1e5 ! avoid /0 at origin
329
-
ELSE
330
#####
factor = 1. / SQRT( r( ir ) ) ! cyl. spreading
331
-
END IF
332
-
END IF
333
-
334
#####
WRITE( ARRFile ) NArr3D( itheta, id, ir )
335
-
336
#####
DO iArr = 1, NArr3D( itheta, id, ir )
337
-
! integers written out as reals below for fast reading in Matlab
338
-
WRITE( ARRFile ) &
339
#####
factor * Arr3D( itheta, id, ir, iArr )%A, &
340
#####
SNGL( RadDeg * Arr3D( itheta, id, ir, iArr )%Phase ), &
341
#####
Arr3D( itheta, id, ir, iArr )%delay, &
342
#####
Arr3D( itheta, id, ir, iArr )%SrcDeclAngle, &
343
#####
Arr3D( itheta, id, ir, iArr )%SrcAzimAngle, &
344
#####
Arr3D( itheta, id, ir, iArr )%RcvrDeclAngle, &
345
#####
Arr3D( itheta, id, ir, iArr )%RcvrAzimAngle, &
346
#####
REAL( Arr3D( itheta, id, ir, iArr )%NTopBnc ), &
347
#####
REAL( Arr3D( itheta, id, ir, iArr )%NBotBnc )
348
-
END DO ! next arrival
349
-
END DO ! next receiver depth
350
-
END DO ! next receiver range
351
-
END DO ! next receiver angle
352
-
353
#####
RETURN
354
-
END SUBROUTINE WriteArrivalsBinary3D
355
-
356
#####
END MODULE ArrMod

Legend:

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