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	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