% ro.m
%
% Calculates phase speeds, Rossby radii, and modal structures (of vertical
% displacement) for an arbitrary stratification profile specified by 
% Brunt-Vaisala frequency profile (supplied by bvfsq.m) under the 
% Bousinnesq approximation, i.e. solves equation (6.11.18) of Gill 
%
% Only the first  "keepmodes" modes are kept. 
%
% John Wilkin, CSIRO Division of Oceanography, Australia
% wilkin@flood.ml.csiro.au
% April 24, 1991
%
% Quick and dirty... hardwire everything
%
keepmodes = 3;
H    = 450;                  % maximum water depth
f    = 0.9e-4;                % Coriolis
n    = 100;                   % number of grid points in z 
delz = H/n;
z    = (1:n-1)*delz - H;
%
coef = -1/delz^2  ./ bvfsq(z);            
A    = diag(-2*ones(n-1,1),0) + diag(ones(n-2,1),1) + diag(ones(n-2,1),-1);
B    = diag(coef) * A;
[V,D]  = eig(B);
eigval = diag(D);
c      = 1.0./sqrt(eigval);
[Y,I]  = sort(-c);
%
speed  = -Y(1:keepmodes);     % phase speed of each mode
radius = speed(1:keepmodes)/f;% Rossby radius of deformation of each mode
mode   = V(:,I(1:keepmodes)); % structure of vertical displacement ech mode
%
% add in points at z=-H,0 to make the plots complete
%
zplt      = [ -H ; z'; 0];
mode   = [zeros(1,keepmodes) ; mode; zeros(1,keepmodes) ];
%
% normalize the modes
%
mode = diag( 1.0 ./ max(abs(mode)))*mode';
%
% plot
%
subplot(121);plot(mode,zplt);title('Modes of vertical displacement');
ylabel('Depth (m)');
subplot(122);plot(bvfsq(z),z);title('N^2');
