Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
function [n2, e] = bm_bvfrq( S, T, P, LAT, DP)
%*************************************************************
% ***** brunt-vaisala freq *****
%
% uses 1980 equation of state
%
% in :
% P pressure decibars
% T temperature deg celsius (ipts-68)
% S salinity psu (pss-78)
% LAT latitude
% DP intervalle de pression (db) sur lequel est calcule n2
%
% out :
% n2 bouyancy freq cph
% e stability 1/m
%
% r. millard feb. 1991
%
%*************************************************************
%------
% BEGIN
%------
% change call list xlat to glat & pass gravity m/s^2
% after formulation of n. p. fofonoff & breck owen's
%
% Calcul de la gravit
% note that sw_g expects height as argument
Z = sw_dpth(P,LAT);
gp = sw_g(LAT,-Z);
g2 = gp .* gp;
g2 = g2 * 1.e-4;
% compute least squares estimate of specific volume anamoly gradient
[m,n] = size(P);
n2 = zeros(size(P));
% Calcul du nombre d'observation sur lequel est calcule n2
dp = P(2)-P(1);
nobs = DP / dp;
% Debut du profil
for i = 1:(nobs/2)
j = 1:i+nobs/2;
[n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end
% Fin du profil
for i = m-(nobs/2):m
j = i-nobs/2:m;
[n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end
% Profil
for i = (nobs/2+1):m-(nobs/2+1)
j = i-nobs/2:i+nobs/2;
[n2(i), e(i)] = cal_bvfrq( S, T, P, g2(i), gp(i), nobs, j);
end
% Corps du calcul
function [n2, e] = cal_bvfrq( S, T, P, g2, gp, nobs, j)
% default gravity & rad/sec to cph conversion
%
radsec = 572.9578;
pav = mean(P(j));
data = sw_svan( S(j), sw_ptmp( S(j), T(j), P(j), pav), pav);
cxy = sum( data .* (P(j)-pav) );
cy = sum( data );
cxx = sum( (P(j)-pav).*(P(j)-pav) );
sigma = sw_pden(S(j), sw_ptmp( S(j), T(j), P(j), pav), P(j), pav);
a0 = cxy/cxx;
v350p = 1./sigma - data;
vbar = v350p(end) + cy/(nobs+1);
dvdp = a0;
e = -g2*dvdp/(vbar)^2;
n2 = radsec*sqrt(abs(e))*sign(e);
% define stability parameter units (1/m)
e = e/gp;
return