% MATLAB file:  eady_model_1.m
% Contour plots for eady model for maximum baroclinic instability.
% Shows results for most unstable mode for beta = 0 but nonzero m.
% This version includes quiver plot of the Q vector for linearized model.
clear all 
%close all
cor = 2*7.292e-5*sin(pi/4);     % Coriolis parameter
bv2 = 1.e-4;                    % buoyancy frequency squared
H = 10000;                      % tropopause height in meters
%Lx =  (2*pi/2.399)*(sqrt(bv2)*H/cor);
sigma = 2.0e-6;                 % static stability parameter
Ld2 = bv2*H^2/cor^2 ;           % radius of deformation squared
Ld = sqrt(Ld2) ; 

k = [0.1:0.1:3] ;   % zonal wavenumber nondimensionalized by Ld ;

m = k ;

% m = pi/3.e6;		            % meridional wavenumber in 1/m

K2 = k.^2 + m.^2 ;
K = sqrt(K2) ; % K is miu in class notes

R = 287;                        % gas constant
Shear = 30/H;                   % shear in sec^-1
%  define the grid points on which fields are computed:
% Define coefficients for the model
A = 2.e7;                       % streamfunction amplitude

Co1 = Shear*H./2.*(1+sqrt(1-4.*cosh(K)./(K.*sinh(K))+4./K2));
Co2 = Shear*H./2.*(1-sqrt(1-4.*cosh(K)./(K.*sinh(K))+4./K2));

for ii=1:length(k)
    if imag(Co1(ii))>imag(Co2(ii))
        Cog(ii) = Co1(ii) ;
        Cod(ii) = Co2(ii) ;
    else
        Cog(ii) = Co2(ii) ;
        Cod(ii) = Co1(ii) ;
    end
    if imag(Co1(ii))==imag(Co2(ii))
        if real(Co1(ii))>real(Co2(ii))
            Cog(ii) = Co1(ii) ;
            Cod(ii) = Co2(ii) ;
        else
            Cog(ii) = Co2(ii) ;
            Cod(ii) = Co1(ii) ;
        end
    end
end
%% plot phase speed
figure
subplot(211)
plot(k,real(Cog),'-ob')
hold on
plot(k,real(Cod),'-b')
grid
xlabel('wavenumber')
title('phase speed')
subplot(212)
% add growth rate plot

%% choose ko and growing or decaying mode for vertical structure
io = find(imag(Cog).*k==max(imag(Cog).*k)) ;

%io=18; just beyond short wave cutoff

ko = k(io) ; 
Lx = Ld/ko ;
Ko = K(io) ; 
Co = Cog(io) ;
mo = ko ;

xx = linspace(0,10.0,31);     % 20 gridpoints in x   
zz = linspace(0,H,21);          % 10 gridpoints in z
[x,z] = meshgrid(xx,zz/H);         % matrix for grid system in x and z

B = A*Ko*(Co+Shear*H)./(Ko*Co-Shear*H) ;

% eddy streamfunction at t = 0
psi = (A*exp(Ko*z)+B*exp(-Ko*z)).*exp(i*ko*x);

% eddy vertical velocity
w = cor/bv2*ko*i*((Co-Shear*z).*(Ko*A*exp(Ko*z)-B*Ko*exp(-Ko*z))...
    +Shear*(A*exp(Ko*z)+B*exp(-Ko*z))).*exp(i*ko*x);
% eddy temperature
T = cor*H/R*(Ko*A*exp(Ko*z)-Ko*B*exp(-Ko*z)).*exp(i*ko*x);  
psimax = max(max(real(psi)));
V = [0:psimax/4:psimax];
V1 = [-psimax:psimax/4:0];
omax = 1.1*max(max(real(w)));
V2 = [0:omax/4:omax];
V3 = [-omax:omax/4:0];
Tmax = 1.1*max(max(real(T)));
V4 = [0:Tmax/4:Tmax];
V5 = [-Tmax:Tmax/4:0];


%% plotting of meridional distribution at top and bottom
figure
subplot(3,1,1)
hold on
contour(x,z,real(psi),V,'k','linewidth',1.5)      
contour(x,z,real(psi),V1,'-k')   % dashed lines for negative contours
contour(x,z,real(psi),[0,0],'-k','linewidth',2.5)   % thick zero line
title('streamfunction (m^2/s)'), ylabel('z  (km)')
grid on
