1 #====================================================================== |
|
2 # L I B G M . P L |
|
3 # doc: Sun Feb 20 14:43:47 2011 |
|
4 # dlm: Tue Nov 18 12:42:30 2014 |
|
5 # (c) 2011 A.M. Thurnherr |
|
6 # uE-Info: 22 53 NIL 0 0 70 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # HISTORY: |
|
10 # Feb 20, 2011: - created |
|
11 # Feb 28, 2011: - cosmetics |
|
12 # Mar 28, 2012: - BUG: N had been ignored (but only affects vertical |
|
13 # wavelengths > 1000m in any was significantly |
|
14 # - changed from Munk eqn 9.23b to 9.23a, which also |
|
15 # affects only long wavelengths |
|
16 # - return nan for omega outside internal-wave band |
|
17 # Mar 29, 2012: - re-wrote using definition of B(omega) from Munk (1981) |
|
18 # Aug 23, 2012: - cosmetics? |
|
19 # Sep 7, 2012: - made N0, E0, b, jstar global |
|
20 # Dec 28, 2012: - added allowance for small roundoff error to Sw() |
|
21 # Oct 6, 2014: - made omega optional in Sw() |
|
22 # Nov 18, 2014: - made b & jstar mandatory for Sw() |
|
23 |
|
24 require "$ANTS/libEOS83.pl"; |
|
25 |
|
26 my($pi) = 3.14159265358979; |
|
27 |
|
28 #====================================================================== |
|
29 # Global Constants |
|
30 #====================================================================== |
|
31 |
|
32 $GM_N0 = 5.24e-3; # rad/s # reference stratification (from Gregg + Kunze, 1991) |
|
33 $GM_E0 = 6.3e-5; # dimensionless # spectral level (Munk 1981) |
|
34 $GM_b = 1300; # m # pycnocline e-folding scale |
|
35 $GM_jstar = 3; # dimless # peak mode number |
|
36 |
|
37 #====================================================================== |
|
38 # Vertical velocity spectral density |
|
39 # |
|
40 # Units: K.E. per frequency per wavenumber [m^2/s^2*1/s*1/m = m^3/s] |
|
41 # Version: GM79? |
|
42 # |
|
43 # E. Kunze (email, Feb 2011): The GM vertical velocity w spectrum is described by |
|
44 # |
|
45 # S[w](omega, k_z) = PI*E_0*b*{f*sqrt(omega^2-f^2)/omega}*{j*/(k_z + k_z*)^2} |
|
46 # |
|
47 # where E_0 = 6.3 x 10^-5 is the dimensionless spectral level, b = 1300 m is |
|
48 # the pycnocline lengthscale, j* = 3 the peak mode number and k_z* the |
|
49 # corresponding vertical wavenumber. The flat log-log spectrum implies w is |
|
50 # dominated by near-N frequencies (where we know very little though Yves |
|
51 # Desaubies wrote some papers back in the late 70's/early 80's about the |
|
52 # near-N peak) and low modes. The rms w = 0.6 cm/s, right near your noise |
|
53 # level. Interestingly, the only N dependence is in m and m*. As far |
|
54 # as I know, little is known about its intermittency compared to horizontal |
|
55 # velocity. Since w WKB-scales inversely with N, the largest signals should |
|
56 # be in the abyss where you therefore likely have the best chance of |
|
57 # measuring it. |
|
58 # |
|
59 # E. Kunze (email, Sep 19, 2013): |
|
60 # |
|
61 # S[w](omega, m) = PI*E0*b*[f*sqrt(omega^2-f^2)/omega]*[jstar/(m+mstar)^2] with |
|
62 # |
|
63 # S[w](m) = PI*E0*b*N*f*[jstar/(m+mstar)^2] |
|
64 # |
|
65 # where the nondimensional spectral energy level E0 = 6.3e-5, stratification |
|
66 # lengthscale b = 1500 m, jstar = 3, mstar = jstar*PI*N/b/N0, and N0 = |
|
67 # 5.3e-3 rad/s. |
|
68 # |
|
69 # NOTES: |
|
70 # - b=1500m is a likely typo, as Gregg & Kunze (1991) have b=1300m |
|
71 # - k_z == m |
|
72 # |
|
73 #====================================================================== |
|
74 |
|
75 sub m($$) # vertical wavenumber as a function of mode number & stratification params |
|
76 { |
|
77 my($j,$N,$omega) = @_; |
|
78 |
|
79 return defined($omega) && ($omega <= $GM_N0) |
|
80 ? $pi / $GM_b * sqrt(($N**2 - $omega**2) / ($GM_N0**2 - $omega**2)) * $j |
|
81 : $pi * $j * $N / ($GM_b * $GM_N0); # valid, except in vicinity of buoyancy turning frequency (Munk 1981, p.285) |
|
82 } |
|
83 |
|
84 sub B($) # structure function (omega dependence) |
|
85 { # NB: f must be defined |
|
86 my($omega) = @_; |
|
87 croak("coriolis parameter not defined\n") |
|
88 unless defined($f); |
|
89 return 2 / $pi * $f / $omega / sqrt($omega**2 - $f**2); |
|
90 } |
|
91 |
|
92 |
|
93 sub Sw(@) |
|
94 { |
|
95 my($omega,$m,$lat,$b,$jstar,$N) = |
|
96 &antsFunUsage(-5,'fff','[frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]> <b[m]> <j*>',@_); |
|
97 |
|
98 if (defined($N)) { # Sw(omega,m) |
|
99 local($f) = abs(&f($lat)); |
|
100 $omega += $PRACTICALLY_ZERO if ($omega < $f); |
|
101 $omega -= $PRACTICALLY_ZERO if ($omega > $N); |
|
102 return nan if ($omega < $f || $omega > $N); |
|
103 my($mstar) = &m($jstar,$N,$omega); |
|
104 return $GM_E0 * $b * 2 * $f**2/$omega**2/B($omega) * $jstar / ($m+$mstar)**2; |
|
105 } else { # Sw(m) |
|
106 $N = $lat; # shift arguments to account for missing omega |
|
107 $lat = $m; |
|
108 local($f) = abs(&f($lat)); |
|
109 $m = $omega; |
|
110 undef($omega); |
|
111 my($mstar) = &m($jstar,$N); |
|
112 return $pi * $GM_E0 * $b * $N * $f * $jstar / ($m+$mstar)**2; |
|
113 } |
|
114 } |
|
115 |
|
116 #---------------------------------------------------------------------- |
|
117 # GM76, as per Gregg and Kunze (JGR 1991) |
|
118 # - beta is vertical wavenumber (m above) |
|
119 #---------------------------------------------------------------------- |
|
120 |
|
121 sub Su($$) |
|
122 { |
|
123 my($beta,$N) = @_; |
|
124 |
|
125 my($beta_star) = &m($GM_jstar,$N); # A3 |
|
126 return 3*$GM_E0*$GM_b**3*$GM_N0**2 / (2*$GM_jstar*$pi) / (1+$beta/$beta_star)**2; # A2 |
|
127 } |
|
128 |
|
129 sub Su_z($$) |
|
130 { |
|
131 my($beta,$N) = &antsFunUsage(2,'ff','<vertical wavenumber[rad/m]> <N[rad/s]>',@_); |
|
132 return $beta**2 * &Su($beta,$N); |
|
133 } |
|
134 |
|
135 1; |
|