[Sea-mat] update for matlab function

[John Klinck]

Attached is  the matlab function for calculating dynamic normal modes 
(dynmodes.m) which I submitted a few years ago. The header now contains a more 
extensive explanation of what the routine does. I get a few requests per year to 
explain the routine, which was poorly documented in its original form. Please 
replace the klinck.tar.Z and klinck.zip with this file (unless you use tar or 
zip wrappers to avoid differences in end of file for different computers).

Let me know if you have any questions. 

John

John M. Klinck
Center for Coastal Physical Oceanography
Department of Ocean, Earth and Atmospheric Sciences
Crittenton Hall
Old Dominion University
Norfolk, Virginia 23529
phone: 757 683-6005, fax: 757 683-5550 email: klinck at ccpo.odu.edu
-------------- next part --------------
function [wmodes,pmodes,ce]=dynmodes(Nsq,p,nmodes)
% DYNMODES calculates ocean dynamic vertical modes
%  taking a column vector of Brunt-Vaisala values (Nsq) at
%  different pressures (p) and calculating some number of 
%  dynamic modes (nmodes). 
%  Note: The input pressures need not be uniformly spaced, 
%    and the deepest pressure is assumed to be the bottom.
%
%  The program is based on the developement in section 6.13 
%  (p. 167ff) of Atmosphere-Ocean Dynamics by Adrian Gill, 
%  Academic Press, 1982.
%
%  The variables describing flow in the ocean or atmosphere are separated
%  into vertical (depending only on height or depth) and horizontal and
%  time. For example, the vertical velocity (w) is represented as a sum of
%  different vertical structures as 
%  $$ 
%  w = \sum_{n=0}^\infty \hat h_n(z) \hat w_n (x,y,t).  
%  $$
%
%  There is a related set of vertical functions ($\hat p_n(z)$) which
%  represent the vertical structure of horizontal velocity and pressure
%  anomaly. In the program below, wmodes is a matrix containing the $\hat
%  h_n(z)$ functions and pmodes are the $\hat p_n(z)$ functions.  wmodes(:,n)
%  or pmodes(:,n) is the nth vertical structure at the depths (p) provided as
%  input.
%
%  Each set of functions is mutually orthogonal with appropriate weighting. 
%
%  The governing equation for the vertical structures (6.11.18 in Gill) is 
%  $$ 
%    d^2/dz^2 \hat h_n + (N^2(z) /c_e) \hat h_n = 0, n = 1, 2, ... ,
%  $$
%  where $N^2(z)$ is the Brunt-Vaiasala freqency squared and $c_e$ is the
%  eigenvalue, which is the phase speed of the internal gravity wave with the
%  indicated vertical structure. The boundary conditions are that $\hat h$
%  vanishes at the surface and bottom (the rigid lid assumption). The
%  solution for n=0 is a special case described in Gill which is very weakly
%  affected by the density structure.
%
%  The vertical velocity structure (wmodes) is calculated at the depths (p)
%  for which $N^2$ is provided on input. The differential equation above is
%  converted to a set of algebraic equations using second order central
%  finite differences. The linear algebraic equation is 
%  $A * wmodes + \lambda B * wmodes = 0$, where B contains $N^2$ along its
%  diagonal. The matlab routine eig(A,B) calculate the eigenvalues $\lambda$
%  which are the reciprocals of the wave speed $c_e$.
%
%  USAGE: [wmodes,pmodes,ce]=dynmodes(Nsq,p,nmodes);
%                               or
%                            dynmodes;  % to show demo 
%
%   Inputs:   Nsq = column vector of Brunt-Vaisala buoyancy frequency (s^-2)
%               p = column vector of pressure (decibars)
%          nmodes = number of vertical modes to calculate 
%  
%   Outputs:   wmodes = vertical velocity structure
%              pmodes = horizontal velocity structure
%                  ce = modal speed (m/s)
%  developed by John Klinck. July, 1999
%  send comments and corrections to klinck at ccpo.odu.edu

if nargin<1
   help(mfilename);
   nplot=3;
%      Calling the function with no input arguements runs
%      one of these two test problems. 
%      This also serves as an example driver to calculate 
%      vertical modes.
%
%      problem 1
%    solution is h = ho sin(z /ce) where ce = 1 / n pi
%     ce = 0.3183 / n = [ 0.3183 0.1591 0.1061]
%  p=0:.05:1;
%  z=-p;
%  n=length(p);
%  Nsq(1:n)=1;
%  nmodes=3;
%
%      problem 2
%    solution is h = ho sin(No z /ce) where ce = No H / n pi
%    for No=1.e-3 and H = 400, the test values are 
%     ce = 0.127324 / n = [ 0.127324, 0.063662, 0.042441]
%
  p=0:10:400;
  z=-p;
  n=length(p);
  Nsq(1:n)=1.e-6;
  nmodes=3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%
  [wmodes,pmodes,ce]=dynmodes(Nsq,p,nmodes);

  figure(1)
  plot(Nsq,z);
  title('Buoyancy Frequency Squared (s^{-2})')

  figure(2)
  plot(ce(1:nplot),'r:o');
  title(' Modal Speed (m/s)')

  figure(3)
  plot(wmodes(:,1:nplot),z);
  title('Vertical Velocity Structure')

  figure(4)
  plot(pmodes(:,1:nplot),z);
  title('Horizontal Velocity Structure')

  figure(gcf)
  return
end
  
rho0=1028;
%    convert to column vector if necessary
[m,n] = size(p);
if n == 1
   p=p';
end
[m,n] = size(Nsq);
if n == 1
   Nsq=Nsq';
   n=m;
end

%                 check for surface value
if p(1) > 0
%             add surface pressure with top Nsq value
    z(1)=0;
    z(2:n+1)=-p(1:n);
    N2(1)=Nsq(1);
    N2(2:n+1)=Nsq(1:n);
    nz=n+1;
else
    z=-p;
    N2=Nsq;
    nz=n;
end

%          calculate depths and spacing
%        spacing
dz(1:nz-1)=z(1:nz-1)-z(2:nz);
%        midpoint depth
zm(1:nz-1)=z(1:nz-1)-.5*dz(1:nz-1)'';
%        midpoint spacing
dzm=zeros(1,nz);
dzm(2:nz-1)=zm(1:nz-2)-zm(2:nz-1);
dzm(1)=dzm(2);
dzm(nz)=dzm(nz-1);

%        get dynamic modes
A = zeros(nz,nz);
B = zeros(nz,nz);
%             create matrices   
for i=2:nz-1
  A(i,i) = 1/(dz(i-1)*dzm(i))  + 1/(dz(i)*dzm(i));
  A(i,i-1) = -1/(dz(i-1)*dzm(i));
  A(i,i+1) = -1/(dz(i)*dzm(i));
end
for i=1:nz
  B(i,i)=N2(i);
end
%             set boundary conditions
A(1,1)=-1.;
A(nz,1)=-1.;

[wmodes,e] = eig(A,B);

%          extract eigenvalues
e=diag(e);
%
ind=find(imag(e)==0);
e=e(ind);
wmodes=wmodes(:,ind);
%
ind=find(e>=1.e-10);
e=e(ind);
wmodes=wmodes(:,ind);
%
[e,ind]=sort(e);
wmodes=wmodes(:,ind);

nm=length(e);
ce=1./sqrt(e);
%                   create pressure structure
pmodes=zeros(size(wmodes));

for i=1:nm
%           calculate first deriv of vertical modes
  pr=diff(wmodes(:,i));   
  pr(1:nz-1)= pr(1:nz-1)./dz(1:nz-1)';
  pr=pr*rho0*ce(i)*ce(i);
%       linearly interpolate back to original depths
  pmodes(2:nz-1,i)=.5*(pr(2:nz-1)+pr(1:nz-2));
  pmodes(1,i)=pr(1);
  pmodes(nz,i)=pr(nz-1);
end