function [V,D,Rn,Cn]=dynmods(N,z,nmods,lat)
%
%       [V,D,Rn,Cn]=dynmods(N,z,nmods,lat)
%
% Calculates vertical modes [V], phase speeds [Cn], Rossby radii [Rn],
% for an arbitrary Brunt-Vaisala frequency profile, N(z).   
%  i.e. solves equation (2.2) of Flierl, 1978, Dyn.Atm.Oce.,2, 341-381.
%
%  Inputs:
%       N,z : Equal length vectors with values of N (radians/sec)
%         at depths z(meters) at fixed depth intervals.
%       nmods : Number of modes to compute.
%       lat : Latitude (degrees)
%  Outputs, sorted by mode number (1 ... nmods):  
%       V : Vertical modes (Eigenvector) matrix, one column per mode,
%               orthonormalized.
%       D : Row vector of Eigenvalues
%       Rn : Row vector of Rossby radii  (km)
%       Cn : Row vector of Rossby phase speeds (m/sec)
 
 
%       Based on an algorithm generously provided by:
% John Wilkin, CSIRO Division of Oceanography, Australia
% wilkin@flood.ml.csiro.au
%
%       Modified to solve Flierl's equation (2.2) 
%  Ref:  Flierl, 1978, Dyn.Atm.Oce.,2, 341-381
%  instead of equation (6.11.18) in Gill's text as originally 
%  designed by Wilkin.   Modification is based on adaptation of a fortran
%  algorithm by David Adamec, NASA, GSFC.  
%  Other modifications to handle irregularly [vertically]
%  spaced data for N, as well as other minor adaptations.
%       Gary S.E. Lagerloef        November 1991
%       Science Applications International Corporation
 
%       Science Applications International Corporation
%       13400B Northup Way NE, Bellevue, WA  98005  USA
%       (206-747-7152)   Omnet: g.lagerloef
 
if all([size(N)] > 1);error('N must be a vector'),end
if all([size(z)] > 1);error('z must be a vector'),end
if length(N) ~= length(z);error('N and z must be equal length'),end
if any([size(lat)] > 1);error('lat must be a scalar'),end
if any(diff(diff(z)) ~=0) | any(diff(z) == 0);
  error('Depths must be at equally spaced intervals.\n'),
  fprintf('Try Matlab "spline.m" utility.'),end
 
f = 4*pi*sin(lat*pi/180)/86400; f=abs(f)          % Coriolis
delz=z(2)-z(1); 
 
%  Estimate N at the mid-grid points:
n=length(N); X = N(1:n-1) + .5*diff(N);
 
X = X.^2;   % N squared
 
n=length(X);fprintf('Length X = %g\n',n)
X = f.^2 ./ X / delz^2;         % f^2 / N^2 / delz^2
X=X(:);  %make sure X is a column vector
 
%       Set up diagonals of matrix for Eigenvalue solution
%   Central diagonal (double difference at gridpoint n)
%    A(1,1) and A(n+1,n+1) satisfy the upper and lower boundary conditions
%    and ensure symetrical matrix so that eigenvectors will be orthogonal.
 
A = diag(-[X(1); X(1:n-1)+X(2:n); X(n)]);
 
%   Add upper and lower diagonals
%   Upper diagonal for gridpoint n+1, lower diagonal for n-1
 
 
A = A + diag(X,1) + diag(X,-1);
 
[am,an]=size(A);fprintf('Size A = %g, %g\n',am,an);
 
%       Solve for Eigenvectors, Eigenvalues 
fprintf('Solving for Eigenvectors & Eigenvalues\n\n')
[V,D]  = eig(A);        % V are Eigenvectors = vertical modes
 
% Reverse Eigenvalue sign because Matlab's eig.m returns neg values
D = -D;
D = diag(D);  D=D(:)';          % Make Eigenvalues a row vector
 
%       Invert and sort Eigenvalues
[D,I] = sort(D);
Y      = 1.0./sqrt(D);
 
%  Mode 1 eigenvalue will be near zero if barotropic
fprintf('Mode information:   %9.0f  %9.0f  %9.0f\n',1,2,3), 
 D = D(1:nmods) * 1.e6;
fprintf('Eigenvalues (km^-2):   %8.4e  %8.4e  %8.4e\n',D(1),D(2),D(3)),
 Rn  = Y(1:nmods)/1000; % Rossby radius of each mode (kilometers)
fprintf('Rossby radii (km):     %8.4e  %8.4e  %8.4e\n',Rn(1),Rn(2),Rn(3)),
 Cn = Y(1:nmods)*f;     % Phase speed of each mode (m/sec)
fprintf('Phase speeds (m/s):    %8.4e  %8.4e  %8.4g\n',Cn(1),Cn(2),Cn(3)),
V = V(:,I(1:nmods));    % vertical structure of each mode
fprintf('\n Warning:  Mode 1 is barotropic if its Eigenvalue is near zero.')
fprintf('\n In that case, the Rossby radius and phase speed for mode 1')
fprintf('\n will be in error (excessively large)')
fprintf('\n Type any key to continue\n'),pause(2)
 
%
% Ortho-normalize the modes according to Flierl's Eq (2.3)
V = V * sqrt(n);
% Change signs of Eigenvectors so that they are all positive at surface:
Vo=sign(V(1,:)); 
for nn=1:nmods; if Vo(nn)==0; Vo(nn)=1; end, end,
V = V * diag(Vo);
 
%
% plot
%
subplot(121);plot(V,-abs(z));title('Vertical Modes');
ylabel('Depth (m)'),grid;
subplot(122);plot(N,-abs(z));title('Brunt-Vaisala Frequency (N)');
ylabel('Depth (m)');grid,  subplot(111);
 
% fini

