1 #====================================================================== |
1 #====================================================================== |
2 # L I B G M . P L |
2 # L I B G M . P L |
3 # doc: Sun Feb 20 14:43:47 2011 |
3 # doc: Sun Feb 20 14:43:47 2011 |
4 # dlm: Sat Apr 20 20:05:51 2013 |
4 # dlm: Mon Oct 6 09:48:22 2014 |
5 # (c) 2011 A.M. Thurnherr |
5 # (c) 2011 A.M. Thurnherr |
6 # uE-Info: 62 0 NIL 0 0 72 2 2 4 NIL ofnI |
6 # uE-Info: 21 10 NIL 0 0 70 2 2 4 NIL ofnI |
7 #====================================================================== |
7 #====================================================================== |
8 |
8 |
9 # HISTORY: |
9 # HISTORY: |
10 # Feb 20, 2011: - created |
10 # Feb 20, 2011: - created |
11 # Feb 28, 2011: - cosmetics |
11 # Feb 28, 2011: - cosmetics |
16 # - return nan for omega outside internal-wave band |
16 # - return nan for omega outside internal-wave band |
17 # Mar 29, 2012: - re-wrote using definition of B(omega) from Munk (1981) |
17 # Mar 29, 2012: - re-wrote using definition of B(omega) from Munk (1981) |
18 # Aug 23, 2012: - cosmetics? |
18 # Aug 23, 2012: - cosmetics? |
19 # Sep 7, 2012: - made N0, E0, b, jstar global |
19 # Sep 7, 2012: - made N0, E0, b, jstar global |
20 # Dec 28, 2012: - added allowance for small roundoff error to Sw() |
20 # Dec 28, 2012: - added allowance for small roundoff error to Sw() |
|
21 # Oct 6, 2014: - made omega optional in Sw() |
21 |
22 |
22 require "$ANTS/libEOS83.pl"; |
23 require "$ANTS/libEOS83.pl"; |
23 |
24 |
24 my($pi) = 3.14159265358979; |
25 my($pi) = 3.14159265358979; |
25 |
26 |
51 # level. Interestingly, the only N dependence is in m and m*. As far |
52 # level. Interestingly, the only N dependence is in m and m*. As far |
52 # as I know, little is known about its intermittency compared to horizontal |
53 # as I know, little is known about its intermittency compared to horizontal |
53 # velocity. Since w WKB-scales inversely with N, the largest signals should |
54 # velocity. Since w WKB-scales inversely with N, the largest signals should |
54 # be in the abyss where you therefore likely have the best chance of |
55 # be in the abyss where you therefore likely have the best chance of |
55 # measuring it. |
56 # measuring it. |
|
57 # |
|
58 # E. Kunze (email, Sep 19, 2013): |
|
59 # |
|
60 # S[w](omega, m) = PI*E0*b*[f*sqrt(omega^2-f^2)/omega]*[jstar/(m+mstar)^2] with |
|
61 # |
|
62 # S[w](m) = PI*E0*b*N*f*[jstar/(m+mstar)^2] |
|
63 # |
|
64 # where the nondimensional spectral energy level E0 = 6.3e-5, stratification |
|
65 # lengthscale b = 1500 m, jstar = 3, mstar = jstar*PI*N/b/N0, and N0 = |
|
66 # 5.3e-3 rad/s. |
|
67 # |
|
68 # NOTES: |
|
69 # - b=1500m is a likely typo, as Gregg & Kunze (1991) have b=1300m |
|
70 # - k_z == m |
|
71 # |
56 #====================================================================== |
72 #====================================================================== |
57 |
73 |
58 sub m($$) # vertical wavenumber as a function of mode number & stratification params |
74 sub m($$) # vertical wavenumber as a function of mode number & stratification params |
59 { |
75 { |
60 my($j,$N,$omega) = @_; |
76 my($j,$N,$omega) = @_; |
61 |
77 |
62 return defined($omega) && ($omega <= $GM_N0) |
78 return defined($omega) && ($omega <= $GM_N0) |
63 ? $pi / $GM_b * sqrt(($N**2 - $omega**2) / ($GM_N0**2 - $omega**2)) * $j |
79 ? $pi / $GM_b * sqrt(($N**2 - $omega**2) / ($GM_N0**2 - $omega**2)) * $j |
64 : $pi * $j * $N / ($GM_b * $GM_N0); # valid, except in vicinity of buoyancy turning frequency (p. 285) |
80 : $pi * $j * $N / ($GM_b * $GM_N0); # valid, except in vicinity of buoyancy turning frequency (Munk 1981, p.285) |
65 } |
81 } |
66 |
82 |
67 sub B($) # structure function (omega dependence) |
83 sub B($) # structure function (omega dependence) |
68 { # NB: f must be defined |
84 { # NB: f must be defined |
69 my($omega) = @_; |
85 my($omega) = @_; |
71 unless defined($f); |
87 unless defined($f); |
72 return 2 / $pi * $f / $omega / sqrt($omega**2 - $f**2); |
88 return 2 / $pi * $f / $omega / sqrt($omega**2 - $f**2); |
73 } |
89 } |
74 |
90 |
75 |
91 |
76 sub Sw($$$$) |
92 sub Sw(@) |
77 { |
93 { |
78 my($omega,$m,$lat,$N) = &antsFunUsage(4,'fff','<frequency[1/s]> <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]>',@_); |
94 my($omega,$m,$lat,$N) = &antsFunUsage(-3,'fff','[frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]>',@_); |
79 |
|
80 local($f) = abs(&f($lat)); |
|
81 $omega += $PRACTICALLY_ZERO if ($omega < $f); |
|
82 $omega -= $PRACTICALLY_ZERO if ($omega > $N); |
|
83 return nan if ($omega < $f || $omega > $N); |
|
84 |
|
85 my($GM_b) = 1300; #m # pycnocline lengthscale |
95 my($GM_b) = 1300; #m # pycnocline lengthscale |
86 |
96 |
87 my($mstar) = &m($GM_jstar,$N,$omega); |
97 if (defined($N)) { # Sw(omega,m) |
88 |
98 local($f) = abs(&f($lat)); |
89 return $GM_E0 * $GM_b * 2 * $f**2/$omega**2/B($omega) * $GM_jstar / ($m+$mstar)**2; |
99 $omega += $PRACTICALLY_ZERO if ($omega < $f); |
|
100 $omega -= $PRACTICALLY_ZERO if ($omega > $N); |
|
101 return nan if ($omega < $f || $omega > $N); |
|
102 my($mstar) = &m($GM_jstar,$N,$omega); |
|
103 return $GM_E0 * $GM_b * 2 * $f**2/$omega**2/B($omega) * $GM_jstar / ($m+$mstar)**2; |
|
104 } else { # Sw(m) |
|
105 $N = $lat; # shift arguments to account for missing omega |
|
106 $lat = $m; |
|
107 local($f) = abs(&f($lat)); |
|
108 $m = $omega; |
|
109 undef($omega); |
|
110 my($mstar) = &m($GM_jstar,$N); |
|
111 return $pi * $GM_E0 * $GM_b * $N * $f * $GM_jstar / ($m+$mstar)**2; |
|
112 } |
90 } |
113 } |
91 |
114 |
92 #---------------------------------------------------------------------- |
115 #---------------------------------------------------------------------- |
93 # GM76, as per Gregg and Kunze (JGR 1991) |
116 # GM76, as per Gregg and Kunze (JGR 1991) |
94 # - beta is vertical wavenumber (m above) |
117 # - beta is vertical wavenumber (m above) |