Warning, /verification/obcs_ctrl/input_ad/VERT_FSFB2.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit bce9e67d on 2011-04-20 19:18:26 UTC
bce9e67dfb Matt*0001 function [c2, Psi, G, N2, Pmid] = VERT_FSFB2(N2,Pmid)
0002 %function [c2, Psi, G, N2, Pmid] = VERT_FSFB2(N2,Pmid)
0003 %
0004 % VERT_FSFB.m
0005 %
0006 % Gabriel A. Vecchi - May 12, 1998
0007 %%%%%%%%%%%%%%%%
0008 %
0009 % Solves the discretized wave projection problem
0010 % given the vertical profiles of Temperature, Salinity, Pressure
0011 % and the depth inteval length.
0012 %
0013 % Uses the seawater function sw_bfrq to calculate N2.
0014 %%%%%%%%%%%%%%%%
0015 %
0016 % Arguments:
0017 % T = temperature vector at same depths as salinity and pressure.
0018 % S = salinity vector at same depths as temperature and pressure.
0019 % P = pressure vector at same depths as temperature and salinity.
0020 % Dz = length of depth interval in meters.
0021 %%%%%%%%%%%%%%%%
0022 %
0023 % Returns:
0024 % c2 = vector of square of the wavespeed.
0025 % Psi = matrix of eigenvectors (horizontal velocity structure functions).
0026 % G = matrix of integral of eigenvectors (vertical velocity structure functions).
0027 % N2 = Brunt-Vaisla frequency calculated at the midpoint pressures.
0028 % Pmid = midpoint pressures.
0029 %%%%%%%%%%%%%%%%
0030
0031 % Find N2 - get a M-1 sized vector, at the equator.
0032 %[N2,crap,Pmid] = sw_bfrq(S,T,P,0);
0033
0034 for i = 1:length(N2)
0035 if N2(i) < 0
0036 N2(i) = min(abs(N2));
0037 end;
0038 end;
0039
0040 % bdc: needs equally-spaced depths!
0041 Dz= median(diff(Pmid));
0042
0043 % add a point for the surface
0044 M = length(N2)+1;
0045
0046 % Fill in D - the differential operator matrix.
0047 % Surface (repeat N2 from midpoint depth)
0048 D(1,1) = -2/N2(1);
0049 D(1,2) = 2/N2(1);
0050 % Interior
0051 for i = 2 : M-1,
0052 D(i,i-1) = 1/N2(i-1);
0053 D(i,i) = -1/N2(i-1)-1/N2(i);
0054 D(i,i+1) = 1/N2(i);
0055 end
0056 % Bottom
0057 D(M,M-1) = 2/N2(M-1);
0058 D(M,M) = -2/N2(M-1);
0059 D=-D./(Dz*Dz);
0060 %bdc: no need for A?
0061 % A = eye(M);
0062
0063 % Calculate generalized eigenvalue problem
0064 % bdc: eigs gets top M-1
0065 %[Psi,lambda] = eigs(D,[],M-1);
0066 % use eig:
0067 [Psi,lambda] = eig(D);
0068
0069 % Calculate square of the wavespeed.
0070 c2 = sum(lambda);
0071 c2=1./c2;
0072
0073 Psi = fliplr(Psi);
0074 c2 = fliplr(c2);
0075 for i=1:size(Psi,2)
0076 Psi(:,i) = Psi(:,i)/Psi(1,i);
0077 end
0078
0079 % normalize?
0080 G = INTEGRATOR(M,Dz)*Psi;
0081
0082 function [INT] = INTEGRATOR(M,Del)
0083 %function [INT] = INTEGRATOR(M,Del)
0084 %
0085 % INTEGRATOR.m
0086 %
0087 % Gabriel A. Vecchi - June 7, 1998
0088 %%%%%%%%%%%%%%%%
0089 % Generates and integration matrix.
0090 % Integrates from first point to each point.
0091 %%%%%%%%%%%%%%%%
0092
0093 INT = tril(ones(M));
0094 INT = INT - 0.5*(eye(M));
0095 INT(:,1) = INT(:,1) - 0.5;
0096 INT = INT*Del;
0097
0098
