Back to home page

MITgcm

 
 

    


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