% MATLAB file: Chapter_6_plot.m    template for problems M6.1, 6.2, 6.3
% Sample quiver and contour plotting for Chapter 6.
% First define the grid points on which fields are computed.
% (Distances are in units of km).
clear all                       % clear the graphics window
%close all                       % clear the workspace
xx=linspace(-3000,3000,30);     % 30 gridpoints in x   
yy=linspace( -1000,1000,10);    % 10 gridpoints in y
[x,y]=meshgrid(xx,yy);          % Sets matrix for grid system in x and y
%   *********Define the vectors to be plotted*********
k = 2*pi/6000;                  % zonal wavenumber in units of 1/ km
l = 2*pi/4000;                  % meridional wavenumber in 1/km
V = 25; %*k/(2*pi/6000);           % wave amplitude speed m/s
U = 30;                         % zonal mean wind
c = 25;                         % phase speed
cor = 1.e-4;                    % Coriolis parameter
beta= 1.67e-11;                 % df/dy
% NOTE that x and y are in km and k is in km-1 for convenience in 
%graphing the factors of 1000 are to convert to meters
%specified geopotential
phi = 55000-U*cor*(y)*1.e3+ V*cor/k*1.e+3*sin(k*x).*cos(l*(y));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dummy fields used to setup graphics
% Students to insert correct formulas for the fields below
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ug = U*ones(size(x));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following fields are given dummy values and must be calculated
% by the student
vg =  zeros(size(x));           % meridional wind
zeta = phi;                     %  vorticity 
advzeta = phi;                  % relative vorticity advection
advbeta = phi;                  % planetary vorticity advection
dzetadt = phi;                  % vorticity tendency
div=phi;                        % divergence field
uad = ug;           % divergent part of zonal ageostrophic wind
vad = vg;           % divergent part of ageostrophic wind    
divag= phi;         % diveegence of ageostrophic wind
vortag = phi;       % vorticity of ageostrophic wind
uan = ug;           % nondivergent ageostrophic wind
van = vg;           % nondivergent ageostrophic wind

%%%% Cecilia modified the script here. Check her work.
% Matlab scales vector plots so that the 
% longest arrow equals the separation between gridpoints. 
% Thus it is difficult to compare the magnitude of vector fields 
% in different figures. You can change this, but it is a pain. 
% First, you must normalize the wind fields and then scale the arrows
% appropriately. See "help quiver" if you are determined.
% It is not required.

% fields for exercise 6.1
ug=U+V*l/k*sin(k*x).*sin(l*y);  %%% the period before the * is essential
vg=V*cos(k*x).*cos(l*y);  %%% the period before the * is essential
zeta = -V*1e-3/k*(k^2+l^2)*sin(k*x).*cos(l*y);
advzeta = U*V*1e-6*(k^2+l^2)*cos(k*x).*cos(l*y);
advbeta = -beta*vg;

% fields for exercise 6.2
totadv=advzeta+advbeta; % from pg 152 unnumbered equations
dzetadt=V*c*(k^2+l^2)*1e-6*cos(k*x).*cos(l*y); % tendency of zeta 
div=(-dzetadt+totadv)/cor; % from Eq 6.18

% fields for exercise 6.3, examine if you wish
% note the exercise says that the nondivergent part of 
% the ageostrophic is proportional to beta. This statement is misleading 
% because vad has a term that is proportional to beta too. 
% The same term appears with the opposite sign in van. Without this 
% term in van, you can easily see that the divergence of (uan,van) is nonzero
% Note  divag = div as it should
uan = -beta*ug/cor.*y*1.e3;    % nondivergent part of ageostrophic wind
van = -beta*vg/cor.*y*1.e3+... 
   beta/cor*V/l*1e+3*cos(k*x).*sin(l*y);  
uad = V*k/cor*(U-c)*1e-3*sin(k*x).*cos(l*y); % divergent part of ageo wind
vad = V*l/cor*(U-c)*1e-3*cos(k*x).*sin(l*y) - ...
   beta/cor*V/l*1e+3*cos(k*x).*sin(l*y);  
divag= V/cor*1e-6*(k^2+l^2)*(U-c)*cos(k*x).*cos(l*y) - beta*vg/cor;
vortag = beta/cor*(ug-y.*zeta);

%%%% End of Cecilia's mods


% phi/g is contoured with contour lines labelled
% figure 1 is template for contour plots on isobaric surfaces
figure(1); clf
subplot(2,1,1)
[cs,h]=contour(x,y,phi/9.8);
clabel(cs,h),title('geopotential height')
xlabel('x (km)'), ylabel('y (km)')
subplot(2,1,2)
[cs,h]=contour(x,y,zeta);
clabel(cs,h)
hold on
quiver(x,y,ug,vg),title('geostrophic wind and vorticity')
xlabel('x (km)'), ylabel('y (km)')

figure(2) ; clf
subplot(2,1,1)
[cs,h]=contour(x,y,advzeta);
clabel(cs,h),title('relative vorticity advection')
xlabel('x (km)'), ylabel('y (km)')
subplot(2,1,2)
[cs,h]=contour(x,y,advbeta);
clabel(cs,h), title('planetary vorticity advection')
xlabel('x (km)'), ylabel('y (km)');
%axis([-3000 3000 -1000 1000])

figure(3); clf
subplot(2,1,1)
[cs,h] = contour(x,y,dzetadt);
clabel(cs,h), title('vorticity tendency')
xlabel('x (km)'), ylabel('y (km)')
subplot(2,1,2)
[cs,h] = contour(x,y,totadv);
clabel(cs,h), title('absolute vorticity advection')
xlabel('x (km)'), ylabel('y (km)')

figure(4); clf
subplot(2,1,1)
[cs,h] = contour(x,y,div);
clabel(cs,h)
title('divergence of the wind')
xlabel('x (km)'), ylabel('y (km)');
%axis([-3000 3000 -1000 1000])


figure(5); clf
subplot(2,1,1)
[cs,h] = contour(x,y,divag);
hold on
quiver(x,y,uad,vad)
clabel(cs,h), title('divergent ageostrophic wind with its divergence')
xlabel('x (km)'), ylabel('y (km)')
subplot(2,1,2)
[cs,h] = contour(x,y,vortag);
clabel(cs,h)
xlabel('x (km)'), ylabel('y (km)')
hold on
quiver(x,y,uan,van)
title('nondivergent ageostrophic wind and ageostrophic vorticity')
xlabel('x (km)'), ylabel('y (km)')

return

% Check that geo wind and nondivergent ageo wind are nondivergent
% using centered finite difference.
% These are not accurate enough to be terribly convincing as the divergence 
% of the geostrophic wind is bigger than the divergence of the 
% ageostropic wind. The inaccuracy is due to the finite difference.
% The analytic calculations in figs 1-5 are very accurate
dx=xx(2)-xx(1); dy=yy(2)-yy(1); dx=dx*1e3; dy=dy*1e3;
tstgeo=NaN*ones(10,30); tstnageo=tstgeo; tstageo=tstgeo;
for i=2:29
    for j=2:9
        tstnageo(j,i) = (van(j+1,i)-van(j-1,i))/(dy*2) + ...
	   + (uan(j,i+1) - uan(j,i-1))/(dx*2);
        tstageo(j,i) = (vad(j+1,i)-vad(j-1,i))/(dy*2) + ...
	   + (uad(j,i+1) - uad(j,i-1))/(dx*2);
        tstgeo(j,i) = (vg(j+1,i)-vg(j-1,i))/(dy*2) + ...
	   + (ug(j,i+1) - ug(j,i-1))/(dx*2);
    end 
end

figure(6); clf
subplot(3,1,1)
[cs,h] = contour(x,y,tstgeo);
clabel(cs,h)
title('divergence of geostrophic wind (should be zero)')
xlabel('x (km)'), ylabel('y (km)');
%axis([-3000 3000 -1000 1000])
subplot(3,1,2)
[cs,h] = contour(x,y,tstnageo);
clabel(cs,h)
title('divergence of nondivergent ageo wind (should be zero)')
xlabel('x (km)'), ylabel('y (km)');
%axis([-3000 3000 -1000 1000])
subplot(3,1,3)
[cs,h] = contour(x,y,tstageo);
clabel(cs,h)
title('divergence of divergent ageo wind (should be nonzero)')
xlabel('x (km)'), ylabel('y (km)');
%axis([-3000 3000 -1000 1000])