obcs_ctrl/input_ad/VERT_FSFB2.m

function [c2, Psi, G, N2, Pmid] = VERT_FSFB2(N2,Pmid)
%function [c2, Psi, G, N2, Pmid] = VERT_FSFB2(N2,Pmid)
%
%    VERT_FSFB.m
%
%     Gabriel A. Vecchi - May 12, 1998
%%%%%%%%%%%%%%%%
% 
%    Solves the discretized wave projection problem
%    given the vertical profiles of Temperature, Salinity, Pressure 
%    and the depth inteval length.
%
%  Uses the seawater function sw_bfrq to calculate N2.
%%%%%%%%%%%%%%%%
%
%    Arguments:
%       T = temperature vector at same depths as salinity and pressure.
%       S = salinity vector at same depths as temperature and pressure.
%       P = pressure vector at same depths as temperature and salinity.
%       Dz = length of depth interval in meters.
%%%%%%%%%%%%%%%%
%
%    Returns:
%       c2 = vector of square of the wavespeed.
%       Psi = matrix of eigenvectors (horizontal velocity structure functions).
%       G  =  matrix of integral of eigenvectors (vertical velocity structure functions).
%       N2 = Brunt-Vaisla frequency calculated at the midpoint pressures.
%       Pmid = midpoint pressures.
%%%%%%%%%%%%%%%%

%  Find N2 - get a M-1 sized vector, at the equator.
%[N2,crap,Pmid] = sw_bfrq(S,T,P,0);

for i = 1:length(N2)
   if N2(i) < 0
      N2(i) = min(abs(N2));
   end;
end;

% bdc: needs equally-spaced depths!
Dz= median(diff(Pmid));

% add a point for the surface
M = length(N2)+1;

%  Fill in D - the differential operator matrix.
%  Surface (repeat N2 from midpoint depth)
D(1,1) = -2/N2(1);
D(1,2) = 2/N2(1);
%  Interior
for i = 2 : M-1,
D(i,i-1) = 1/N2(i-1);
D(i,i) = -1/N2(i-1)-1/N2(i);
D(i,i+1) = 1/N2(i);
end
%  Bottom
D(M,M-1) = 2/N2(M-1);
D(M,M) = -2/N2(M-1);
D=-D./(Dz*Dz);
%bdc: no need for A?
% A = eye(M);

% Calculate generalized eigenvalue problem
% bdc: eigs gets top M-1
%[Psi,lambda] = eigs(D,[],M-1);
% use eig:
[Psi,lambda] = eig(D);

% Calculate square of the wavespeed.
c2 = sum(lambda);
c2=1./c2;

Psi = fliplr(Psi);
c2 = fliplr(c2);
for i=1:size(Psi,2)
  Psi(:,i) = Psi(:,i)/Psi(1,i);
end

% normalize?
G = INTEGRATOR(M,Dz)*Psi;

function [INT] = INTEGRATOR(M,Del)
%function [INT] = INTEGRATOR(M,Del)
%
%    INTEGRATOR.m
%
%     Gabriel A. Vecchi - June 7, 1998
%%%%%%%%%%%%%%%%
%    Generates and integration matrix.
%    Integrates from first point to each point.
%%%%%%%%%%%%%%%%

INT = tril(ones(M));
INT = INT - 0.5*(eye(M));
INT(:,1) = INT(:,1) - 0.5;
INT = INT*Del;


1	mmazloff	1.1	function [c2, Psi, G, N2, Pmid] = VERT_FSFB2(N2,Pmid)
2			%function [c2, Psi, G, N2, Pmid] = VERT_FSFB2(N2,Pmid)
3			%
4			% VERT_FSFB.m
5			%
6			% Gabriel A. Vecchi - May 12, 1998
7			%%%%%%%%%%%%%%%%
8			%
9			% Solves the discretized wave projection problem
10			% given the vertical profiles of Temperature, Salinity, Pressure
11			% and the depth inteval length.
12			%
13			% Uses the seawater function sw_bfrq to calculate N2.
14			%%%%%%%%%%%%%%%%
15			%
16			% Arguments:
17			% T = temperature vector at same depths as salinity and pressure.
18			% S = salinity vector at same depths as temperature and pressure.
19			% P = pressure vector at same depths as temperature and salinity.
20			% Dz = length of depth interval in meters.
21			%%%%%%%%%%%%%%%%
22			%
23			% Returns:
24			% c2 = vector of square of the wavespeed.
25			% Psi = matrix of eigenvectors (horizontal velocity structure functions).
26			% G = matrix of integral of eigenvectors (vertical velocity structure functions).
27			% N2 = Brunt-Vaisla frequency calculated at the midpoint pressures.
28			% Pmid = midpoint pressures.
29			%%%%%%%%%%%%%%%%
30
31			% Find N2 - get a M-1 sized vector, at the equator.
32			%[N2,crap,Pmid] = sw_bfrq(S,T,P,0);
33
34			for i = 1:length(N2)
35			if N2(i) < 0
36			N2(i) = min(abs(N2));
37			end;
38			end;
39
40			% bdc: needs equally-spaced depths!
41			Dz= median(diff(Pmid));
42
43			% add a point for the surface
44			M = length(N2)+1;
45
46			% Fill in D - the differential operator matrix.
47			% Surface (repeat N2 from midpoint depth)
48			D(1,1) = -2/N2(1);
49			D(1,2) = 2/N2(1);
50			% Interior
51			for i = 2 : M-1,
52			D(i,i-1) = 1/N2(i-1);
53			D(i,i) = -1/N2(i-1)-1/N2(i);
54			D(i,i+1) = 1/N2(i);
55			end
56			% Bottom
57			D(M,M-1) = 2/N2(M-1);
58			D(M,M) = -2/N2(M-1);
59			D=-D./(Dz*Dz);
60			%bdc: no need for A?
61			% A = eye(M);
62
63			% Calculate generalized eigenvalue problem
64			% bdc: eigs gets top M-1
65			%[Psi,lambda] = eigs(D,[],M-1);
66			% use eig:
67			[Psi,lambda] = eig(D);
68
69			% Calculate square of the wavespeed.
70			c2 = sum(lambda);
71			c2=1./c2;
72
73			Psi = fliplr(Psi);
74			c2 = fliplr(c2);
75			for i=1:size(Psi,2)
76			Psi(:,i) = Psi(:,i)/Psi(1,i);
77			end
78
79			% normalize?
80			G = INTEGRATOR(M,Dz)*Psi;
81
82			function [INT] = INTEGRATOR(M,Del)
83			%function [INT] = INTEGRATOR(M,Del)
84			%
85			% INTEGRATOR.m
86			%
87			% Gabriel A. Vecchi - June 7, 1998
88			%%%%%%%%%%%%%%%%
89			% Generates and integration matrix.
90			% Integrates from first point to each point.
91			%%%%%%%%%%%%%%%%
92
93			INT = tril(ones(M));
94			INT = INT - 0.5*(eye(M));
95			INT(:,1) = INT(:,1) - 0.5;
96			INT = INT*Del;
97
98