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: Sun Apr 1 11:29:53 2012 |
4 # dlm: Fri Sep 14 12:35:40 2012 |
5 # (c) 2011 A.M. Thurnherr |
5 # (c) 2011 A.M. Thurnherr |
6 # uE-Info: 53 1 NIL 0 0 72 2 2 4 NIL ofnI |
6 # uE-Info: 61 0 NIL 0 0 72 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 |
13 # wavelengths > 1000m in any was significantly |
13 # wavelengths > 1000m in any was significantly |
14 # - changed from Munk eqn 9.23b to 9.23a, which also |
14 # - changed from Munk eqn 9.23b to 9.23a, which also |
15 # affects only long wavelengths |
15 # affects only long wavelengths |
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? |
|
19 # Sep 7, 2012: - made N0, E0, b, jstar global |
18 |
20 |
19 require "$ANTS/libEOS83.pl"; |
21 require "$ANTS/libEOS83.pl"; |
20 |
22 |
21 my($pi) = 3.14159265358979; |
23 my($pi) = 3.14159265358979; |
|
24 |
|
25 #====================================================================== |
|
26 # Global Constants |
|
27 #====================================================================== |
|
28 |
|
29 $GM_N0 = 5.24e-3; # rad/s # reference stratification (from Gregg + Kunze, 1991) |
|
30 $GM_E0 = 6.3e-5; # dimensionless # spectral level (Munk 1981) |
|
31 $GM_b = 1300; # m # pycnocline e-folding scale |
|
32 $GM_jstar = 3; # dimless # peak mode number |
22 |
33 |
23 #====================================================================== |
34 #====================================================================== |
24 # Vertical velocity spectral density |
35 # Vertical velocity spectral density |
25 # |
36 # |
26 # Units: K.E. per frequency per wavenumber [m^2/s^2*1/s*1/m = m^3/s] |
37 # Units: K.E. per frequency per wavenumber [m^2/s^2*1/s*1/m = m^3/s] |
45 |
56 |
46 sub m($$) # vertical wavenumber as a function of mode number & stratification params |
57 sub m($$) # vertical wavenumber as a function of mode number & stratification params |
47 { |
58 { |
48 my($j,$N,$omega) = @_; |
59 my($j,$N,$omega) = @_; |
49 |
60 |
50 my($b) = 1300; #m # stratification e-folding scale (Munk 81) |
|
51 my($N0) = 5.2e-3; #rad/s # extrapolated to surface value (Munk 81) |
|
52 |
|
53 # print(STDERR "omega = $omega, N = $N\n"); |
|
54 return defined($omega) |
61 return defined($omega) |
55 ? $pi / $b * sqrt(($N**2 - $omega**2) / ($N0**2 - $omega**2)) * $j |
62 ? $pi / $GM_b * sqrt(($N**2 - $omega**2) / ($GM_N0**2 - $omega**2)) * $j |
56 : $pi * $j * $N / ($b * $N0); # valid, except in vicinity of buoyancy turning frequency (p. 285) |
63 : $pi * $j * $N / ($GM_b * $GM_N0); # valid, except in vicinity of buoyancy turning frequency (p. 285) |
57 } |
64 } |
58 |
65 |
59 sub B($) # structure function (omega dependence) |
66 sub B($) # structure function (omega dependence) |
60 { # NB: f must be defined |
67 { # NB: f must be defined |
61 my($omega) = @_; |
68 my($omega) = @_; |
65 } |
72 } |
66 |
73 |
67 |
74 |
68 sub Sw($$$$) |
75 sub Sw($$$$) |
69 { |
76 { |
70 my($omega,$m,$lat,$N) = &antsFunUsage(4,'fff','<frequency[1/s]> <vertical wavenumber[1/m]> <lat[deg]> <N[rad/s]>',@_); |
77 my($omega,$m,$lat,$N) = &antsFunUsage(4,'fff','<frequency[1/s]> <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]>',@_); |
71 |
78 |
72 local($f) = abs(&f($lat)); |
79 local($f) = abs(&f($lat)); |
73 return nan if ($omega < $f || $omega > $N); |
80 return nan if ($omega < $f || $omega > $N); |
74 |
81 |
75 my($E0) = 6.3e-5; # dimensionless spectral level |
82 my($GM_b) = 1300; #m # pycnocline lengthscale |
76 my($j_star) = 3; # peak mode number |
|
77 my($b) = 1300; #m # pycnocline lengthscale |
|
78 |
83 |
79 my($mstar) = &m($j_star,$N,$omega); |
84 my($mstar) = &m($GM_jstar,$N,$omega); |
80 |
85 |
81 return $E0 * $b * 2 * $f**2/$omega**2/B($omega) * $j_star / ($m+$mstar)**2; |
86 return $GM_E0 * $GM_b * 2 * $f**2/$omega**2/B($omega) * $GM_jstar / ($m+$mstar)**2; |
|
87 } |
|
88 |
|
89 #---------------------------------------------------------------------- |
|
90 # GM76, as per Gregg and Kunze (JGR 1991) |
|
91 # - beta is vertical wavenumber (m above) |
|
92 #---------------------------------------------------------------------- |
|
93 |
|
94 sub Su($$) |
|
95 { |
|
96 my($beta,$N) = @_; |
|
97 |
|
98 my($beta_star) = &m($GM_jstar,$N); # A3 |
|
99 return 3*$GM_E0*$GM_b**3*$GM_N0**2 / (2*$GM_jstar*$pi) / (1+$beta/$beta_star)**2; # A2 |
|
100 } |
|
101 |
|
102 sub Su_z($$) |
|
103 { |
|
104 my($beta,$N) = &antsFunUsage(2,'ff','<vertical wavenumber[rad/m]> <N[rad/s]>',@_); |
|
105 return $beta**2 * &Su($beta,$N); |
82 } |
106 } |
83 |
107 |
84 1; |
108 1; |