Coverage Report: AttenMod.f90

Generated from GCOV analysis of Fortran source code

58.5%
Lines Executed
53 total lines
45.5%
Branches Executed
22 total branches
87.5%
Calls Executed
8 total calls
0
-
Source:AttenMod.f90
0
-
Graph:AttenMod.gcno
0
-
Data:AttenMod.gcda
0
-
Runs:73
1
-
!! Attenuation module for converting sound speed and attenuation to complex values
2
-
!! Provides routines to convert sound speed and attenuation in user units to complex sound speed
3
-
MODULE AttenMod
4
-
!! Acoustic attenuation calculations including volume attenuation formulas and unit conversions
5
-
6
-
! Attenuation module
7
-
! Routines to convert a sound speed and attenuation in user units to a complex sound speed
8
-
! Includes a formula for volume attenuation
9
-
10
-
USE FatalError
11
-
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT
12
-
13
-
IMPLICIT NONE
14
-
PUBLIC
15
-
16
-
INTEGER, PRIVATE, PARAMETER :: PRTFile = 6
17
-
INTEGER, PARAMETER :: MaxBioLayers = 200
18
-
INTEGER :: iBio, NBioLayers
19
-
REAL (KIND=8) :: T = 20, Salinity = 35, pH = 8, z_bar = 0, FG
20
-
! Francois-Garrison volume attenuation; temperature, salinity, ...
21
-
22
-
TYPE bioStructure
23
-
REAL (KIND=8) :: Z1, Z2, f0, Q, a0
24
-
END TYPE bioStructure
25
-
26
-
TYPE( bioStructure ) :: bio( MaxBioLayers )
27
-
28
-
CONTAINS
29
-
!! Converts real wave speed and attenuation to complex wave speed
30
323*
FUNCTION CRCI( z, c, alpha, freq, freq0, AttenUnit, beta, fT )
31
-
32
-
! Converts real wave speed and attenuation to a single
33
-
! complex wave speed (with positive imaginary part)
34
-
35
-
! AttenUnit
36
-
! 6 Cases: N Nepers/meter
37
-
! M dB/meter (M for Meters)
38
-
! m dB/meter with a power law
39
-
! F dB/m-kHz (F for frequency dependent)
40
-
! W dB/wavelength (W for Wavelength)
41
-
! Q Q
42
-
! L Loss parameter
43
-
!
44
-
! second letter adds volume attenuation according to standard laws:
45
-
! T for Thorp
46
-
! F for Francois Garrison
47
-
! B for biological
48
-
!
49
-
! freq is the current frequency
50
-
! freq0 is the reference frequency for which the dB/meter was specified (used only for 'm')
51
-
52
-
! Returns
53
-
! c real part of sound speed
54
-
! alpha imaginary part of sound speed
55
-
56
-
USE MathConstants
57
-
58
-
REAL (KIND=8), INTENT( IN ) :: freq, freq0, alpha, c, z, beta, fT
59
-
CHARACTER (LEN=2), INTENT( IN ) :: AttenUnit
60
-
REAL (KIND=8) :: f2, omega, alphaT, Thorp, a, FG
61
-
COMPLEX (KIND=8) :: CRCI
62
-
63
323
omega = 2.0 * pi * freq
64
323
alphaT = 0.0
65
-
! Convert to Nepers/m
66
#####
SELECT CASE ( AttenUnit( 1 : 1 ) )
67
-
CASE ( 'N' )
68
#####
alphaT = alpha
69
-
CASE ( 'M' ) ! dB/m
70
#####
alphaT = alpha / 8.6858896D0
71
-
CASE ( 'm' ) ! dB/m with power law
72
#####
alphaT = alpha / 8.6858896D0
73
#####
IF ( freq < fT ) THEN ! frequency raised to the power beta
74
#####
alphaT = alphaT * ( freq / freq0 ) ** beta
75
-
ELSE ! linear in frequency
76
#####
alphaT = alphaT * ( freq / freq0 ) * ( fT / freq0 ) ** ( beta - 1 )
77
-
END IF
78
-
CASE ( 'F' ) ! dB/(m kHz)
79
283
alphaT = alpha * freq / 8685.8896D0
80
-
CASE ( 'W' ) ! dB/wavelength
81
40*
IF ( c /= 0.0 ) alphaT = alpha * freq / ( 8.6858896D0 * c )
82
-
! The following lines give f^1.25 frequency dependence
83
-
! FAC = SQRT( SQRT( freq / 50.0 ) )
84
-
! IF ( c /= 0.0 ) alphaT = FAC * alpha * freq / ( 8.6858896D0 * c )
85
-
CASE ( 'Q' ) ! Quality factor
86
#####
IF ( c * alpha /= 0.0 ) alphaT = omega / ( 2.0 * c * alpha )
87
-
CASE ( 'L' ) ! loss parameter
88
#####
IF ( c /= 0.0 ) alphaT = alpha * omega / c
89
-
CASE ( ' ' )
90
-
! default = do nothing
91
-
CASE DEFAULT
92
#####
WRITE( ERROR_UNIT, * ) 'AttenMod: CRCI: Unknown attenuation unit (1st char): ', AttenUnit( 1 : 1 )
93
323*
CALL ERROUT( 'AttenMod : CRCI', 'Unknown attenuation unit (first character)' )
94
-
END SELECT
95
-
96
-
! added volume attenuation
97
#####
SELECT CASE ( AttenUnit( 2 : 2 ) )
98
-
CASE ( 'T' )
99
#####
f2 = ( freq / 1000.0 ) ** 2
100
-
101
-
! Original formula from Thorp 1967
102
-
! Thorp = 40.0 * f2 / ( 4100.0 + f2 ) + 0.1 * f2 / ( 1.0 + f2 ) ! dB/kyard
103
-
! Thorp = Thorp / 914.4D0 ! dB / m
104
-
! Thorp = Thorp / 8.6858896D0 ! Nepers / m
105
-
106
-
! Updated formula from JKPS Eq. 1.34
107
#####
Thorp = 3.3d-3 + 0.11 * f2 / ( 1.0 + f2 ) + 44.0 * f2 / ( 4100.0 + f2 ) + 3d-4 * f2 ! dB/km
108
#####
Thorp = Thorp / 8685.8896 ! Nepers / m
109
-
110
#####
alphaT = alphaT + Thorp
111
-
CASE ( 'F' ) ! Francois-Garrison
112
10
FG = Franc_Garr( freq / 1000 ) ! dB/km
113
10
FG = FG / 8685.8896 ! Nepers / m
114
10*
alphaT = alphaT + FG
115
-
CASE ( 'B' ) ! biological attenuation per Orest Diachok
116
313*
DO iBio = 1, NBioLayers
117
#####
IF ( z >= bio( iBio )%Z1 .AND. z <= bio( iBio )%Z2 ) THEN
118
#####
a = bio( iBio )%a0 / ( ( 1.0 - bio( iBio )%f0 ** 2 / freq ** 2 ) ** 2 + 1.0 / bio( iBio )%Q ** 2 ) ! dB/km
119
#####
a = a / 8685.8896 ! Nepers / m
120
#####
alphaT = alphaT + a
121
-
END IF
122
-
END DO
123
-
CASE ( ' ' )
124
-
! default = do nothing
125
-
CASE DEFAULT
126
#####
WRITE( ERROR_UNIT, * ) 'AttenMod: CRCI: Unknown volume attenuation unit (2nd char): ', AttenUnit( 2 : 2 )
127
323*
CALL ERROUT( 'AttenMod : CRCI', 'Unknown attenuation unit (second character)' )
128
-
END SELECT
129
-
130
-
! Convert Nepers/m to equivalent imaginary sound speed
131
323
alphaT = alphaT * c * c / omega
132
323
CRCI = CMPLX( c, alphaT, KIND=8 )
133
-
134
323
IF ( alphaT > c ) THEN
135
1
WRITE( PRTFile, * ) 'Complex sound speed: ', CRCI
136
1
WRITE( PRTFile, * ) 'Usually this means you have an attenuation that is way too high'
137
1
CALL ERROUT( 'AttenMod : CRCI ', 'The complex sound speed has an imaginary part > real part' )
138
-
END IF
139
-
140
322
RETURN
141
-
END FUNCTION CRCI
142
-
143
-
! **********************************************************************!
144
-
145
-
!! Computes Francois-Garrison volume attenuation
146
10
FUNCTION Franc_Garr( f )
147
-
148
-
! Francois Garrison formulas for attenuation
149
-
! Based on a Matlab version by D. Jackson APL-UW
150
-
151
-
! mbp Feb. 2019
152
-
! Verified using F-G Table IV
153
-
154
-
! alpha = attenuation (dB/km)
155
-
! f = frequency (kHz)
156
-
! T = temperature (deg C)
157
-
! S = salinity (psu)
158
-
! pH = 7 for neutral water
159
-
! z_bar = depth (m)
160
-
161
-
! Returns
162
-
! alpha = volume attenuation in dB/km
163
-
164
-
REAL (KIND=8), INTENT( IN ) :: f
165
-
REAL (KIND=8) :: Franc_Garr
166
-
REAL (KIND=8) :: c, A1, A2, A3, P1, P2, P3, f1, f2
167
-
! LP: Bug (at least technically): Single-precision and double-precision
168
-
! literals are all mixed together here, both including values which are
169
-
! not representable as floats, thus leading to different results.
170
-
! Practically, these are approximate, experimentally-derived constants, so
171
-
! errors at the 1e-7 scale are completely not meaningful.
172
-
173
10
c = 1412 + 3.21 * T + 1.19 * Salinity + 0.0167 * z_bar
174
-
175
-
! Boric acid contribution
176
10
A1 = 8.86 / c * 10 ** ( 0.78 * pH - 5 )
177
10
P1 = 1
178
10
f1 = 2.8 * sqrt( Salinity / 35 ) * 10 ** ( 4 - 1245 / ( T + 273 ) )
179
-
180
-
! Magnesium sulfate contribution
181
10
A2 = 21.44 * Salinity / c * ( 1 + 0.025 * T )
182
10
P2 = 1 - 1.37D-4 * z_bar + 6.2D-9 * z_bar ** 2
183
10
f2 = 8.17 * 10 ** ( 8 - 1990 / ( T + 273 ) ) / ( 1 + 0.0018 * ( Salinity - 35 ) )
184
-
185
-
! Viscosity
186
10
P3 = 1 - 3.83D-5 * z_bar + 4.9D-10 * z_bar ** 2
187
10
if ( T < 20 ) THEN
188
10
A3 = 4.937D-4 - 2.59D-5 * T + 9.11D-7 * T ** 2 - 1.5D-8 * T ** 3
189
-
else
190
#####
A3 = 3.964D-4 -1.146D-5 * T + 1.45D-7 * T ** 2 - 6.5D-10 * T ** 3
191
-
end if
192
-
193
-
Franc_Garr = A1 * P1 * ( f1 * f ** 2 ) / ( f1 ** 2 + f ** 2 ) + A2 * P2 * ( f2 * f ** 2 ) / ( f2 ** 2 + f ** 2 ) + &
194
10
A3 * P3 * f ** 2
195
-
196
10
END FUNCTION Franc_Garr
197
-
198
#####
END MODULE AttenMod

Legend:

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