dynmodes

 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.

  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 J. Klinck. July, 1999
  send comments and corrections to klinck@ccpo.odu.edu

0001 function [wmodes,pmodes,ce]=dynmodes(Nsq,p,nmodes)
0002 % DYNMODES calculates ocean dynamic vertical modes
0003 %  taking a column vector of Brunt-Vaisala values (Nsq) at
0004 %  different pressures (p) and calculating some number of
0005 %  dynamic modes (nmodes).
0006 %  Note: The input pressures need not be uniformly spaced,
0007 %    and the deepest pressure is assumed to be the bottom.
0008 %
0009 %  USAGE: [wmodes,pmodes,ce]=dynmodes(Nsq,p,nmodes);
0010 %                               or
0011 %                            dynmodes;  % to show demo
0012 %
0013 %     Inputs:     Nsq = column vector of Brunt-Vaisala buoyancy frequency (s^-2)
0014 %                  p = column vector of pressure (decibars)
0015 %           nmodes = number of vertical modes to calculate
0016 %
0017 %       Outputs:   wmodes = vertical velocity structure
0018 %                  pmodes = horizontal velocity structure
0019 %                      ce = modal speed (m/s)
0020 %  developed by J. Klinck. July, 1999
0021 %  send comments and corrections to klinck@ccpo.odu.edu
0022 
0023 if nargin<1
0024    help(mfilename);
0025    nplot=3;
0026 %    test problems
0027 %      problem 1
0028 %    solution is h = ho sin(z /ce) where ce = 1 / n pi
0029 %     ce = 0.3183 / n = [ 0.3183 0.1591 0.1061]
0030 %p=0:.05:1;
0031 %z=-p;
0032 %n=length(p);
0033 %Nsq(1:n)=1;
0034 %
0035 %      problem 2
0036 %    solution is h = ho sin(No z /ce) where ce = No H / n pi
0037 %    for No=1.e-3 and H = 400, the test values are
0038 %     ce = 0.127324 / n = [ 0.127324, 0.063662, 0.042441]
0039 %
0040     p=0:10:400;
0041     z=-p;
0042     n=length(p);
0043     Nsq(1:n)=1.e-6;
0044 
0045     nmodes=3;
0046 
0047     [wmodes,pmodes,ce]=dynmodes(Nsq,p,nmodes);
0048 
0049     figure(1)
0050     plot(Nsq,z);
0051     title('Buoyancy Frequency Squared (s^{-2})')
0052 
0053     figure(2)
0054     plot(ce(1:nplot),'r:o');
0055     title(' Modal Speed (m/s)')
0056 
0057     figure(3)
0058     plot(wmodes(:,1:nplot),z);
0059     title('Vertical Velocity Structure')
0060 
0061     figure(4)
0062     plot(pmodes(:,1:nplot),z);
0063     title('Horizontal Velocity Structure')
0064 
0065         figure(gcf)
0066         return
0067 end
0068   
0069 rho0=1028;
0070 
0071 %    convert to column vector if necessary
0072 [m,n] = size(p);
0073 if n == 1
0074    p=p';
0075 end
0076 [m,n] = size(Nsq);
0077 if n == 1
0078    Nsq=Nsq';
0079    n=m;
0080 end
0081 
0082 %                 check for surface value
0083 if p(1) > 0
0084 %             add surface pressure with top Nsq value
0085     z(1)=0;
0086     z(2:n+1)=-p(1:n);
0087     N2(1)=Nsq(1);
0088     N2(2:n+1)=Nsq(1:n);
0089     nz=n+1;
0090 else
0091     z=-p;
0092     N2=Nsq;
0093     nz=n;
0094 end
0095 
0096 %          calculate depths and spacing
0097 %        spacing
0098 dz(1:nz-1)=z(1:nz-1)-z(2:nz);
0099 %        midpoint depth
0100 zm(1:nz-1)=z(1:nz-1)-.5*dz(1:nz-1)'';
0101 %        midpoint spacing
0102 dzm=zeros(1,nz);
0103 dzm(2:nz-1)=zm(1:nz-2)-zm(2:nz-1);
0104 dzm(1)=dzm(2);
0105 dzm(nz)=dzm(nz-1);
0106 
0107 %        get dynamic modes
0108 A = zeros(nz,nz);
0109 B = zeros(nz,nz);
0110 %             create matrices
0111 for i=2:nz-1
0112   A(i,i) = 1/(dz(i-1)*dzm(i))  + 1/(dz(i)*dzm(i));
0113   A(i,i-1) = -1/(dz(i-1)*dzm(i));
0114   A(i,i+1) = -1/(dz(i)*dzm(i));
0115 end
0116 for i=1:nz
0117   B(i,i)=N2(i);
0118 end
0119 %             set boundary conditions
0120 A(1,1)=-1.;
0121 A(nz,1)=-1.;
0122 
0123 [wmodes,e] = eig(A,B);
0124 
0125 %          extract eigenvalues
0126 e=diag(e);
0127 %
0128 ind=find(imag(e)==0);
0129 e=e(ind);
0130 wmodes=wmodes(:,ind);
0131 %
0132 ind=find(e>=1.e-10);
0133 e=e(ind);
0134 wmodes=wmodes(:,ind);
0135 %
0136 [e,ind]=sort(e);
0137 wmodes=wmodes(:,ind);
0138 
0139 nm=length(e);
0140 ce=1./sqrt(e);
0141 %                   create pressure structure
0142 pmodes=zeros(size(wmodes));
0143 
0144 for i=1:nm
0145 %           calculate first deriv of vertical modes
0146   pr=diff(wmodes(:,i));   
0147   pr(1:nz-1)= pr(1:nz-1)./dz(1:nz-1)';
0148   pr=pr*rho0*ce(i)*ce(i);
0149 %       linearly interpolate back to original depths
0150   pmodes(2:nz-1,i)=.5*(pr(2:nz-1)+pr(1:nz-2));
0151   pmodes(1,i)=pr(1);
0152   pmodes(nz,i)=pr(nz-1);
0153 end