% 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).
% CC modified to illustrate planetary versus relative vort adv 
% from page 153-4

clear all                       % clear the graphics window
%close all                       % clear the workspace
%   *********Define the vectors to be plotted*********
L=3000;
k = 2*pi/L/2;                  % zonal wavenumber in units of 1/ km
l = 2*pi/L;			        % meridional wavenumber in 1/km
V = 25;                         % wave amplitude speed m/s
U = 30;                         % zonal mean wind
c = 25;                         % phase speed
t=0;
cor = 1.e-4;                    % Coriolis parameter
beta= 1.67e-11;                 % df/dy
xx=linspace(-2*L,2*L,30);     % 30 gridpoints in x   
yy=linspace( -L/4,L/4,10);    % 10 gridpoints in y
[x,y]=meshgrid(xx,yy);          % Sets matrix for grid system in x and y
% 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
kk=1e-3; % so A=V/kk*1e3 on pg 153, pg 178 has kk=k;
kk=k; % this helps scale the terms nicely
phi = 55000-U*cor*(y)*1.e3+ V*cor/kk*1.e3*sin(k*(x-c*t)).*cos(l*(y));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dummy fields used to setup graphics
% Students to insert correct formulas for the fields below
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following fields are given dummy values and must be calculated
ug = U*ones(size(x));
vg =  zeros(size(x));           % meridional wind


%%%% End of Cecilia's mods
%%%% 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/kk*sin(k*x).*sin(l*y);  %%% the period before the * is essential
vg=V*k/kk*cos(k*x).*cos(l*y);  %%% the period before the * is essential
zeta = -V/kk*1e-3*(k^2+l^2)*sin(k*(x-c*t)).*cos(l*y);
advzeta = U*V*k/kk*1e-6*(k^2+l^2)*cos(k*(x-c*t)).*cos(l*y);
advbeta = -beta*vg;
%advbeta = -beta*V*k/kk*cos(k*(x-c*t)).*cos(l*y);

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

% below kk=k is hardwired
% 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);
set(gca,'fontsize',16)
clabel(cs,h),title('geopotential height')
xlabel('x (km)'), ylabel('y (km)')
subplot(2,1,2)
[cs,h]=contour(x,y,zeta);
set(gca,'fontsize',16)
clabel(cs,h)
hold on
quiver(x,y,ug,vg),title('geostrophic wind and vorticity')
xlabel('x (km)'), ylabel('y (km)')

figure(2) ;clf

sc = max(max(max(advzeta)),-min(min(advbeta)));
sc=eval(num2str(sc,1));
con=linspace( -sc,sc,31);    % 10 gridpoints in y

subplot(2,1,1)
[cs,h]=contour(x,y,advzeta,con); caxis([-1 1]*sc);%clabel(cs,h),
set(gca,'fontsize',16)
title('relative vorticity advection')
xlabel('x (km)'), ylabel('y (km)')
subplot(2,1,2)
[cs,h]=contour(x,y,advbeta,con); caxis([-1 1]*sc);%clabel(cs,h), 
set(gca,'fontsize',16)
title('planetary vorticity advection')
xlabel('x (km)'), ylabel('y (km)');

figure(4) ;clf
subplot(2,1,1)
[cs,h]=contour(x,y,advzeta); clabel(cs,h),
set(gca,'fontsize',16)
title('relative vorticity advection')
xlabel('x (km)'), ylabel('y (km)')
subplot(2,1,2)
[cs,h]=contour(x,y,advbeta);clabel(cs,h), 
set(gca,'fontsize',16)
title('planetary vorticity advection')
xlabel('x (km)'), ylabel('y (km)');


%return
figure(3); clf
subplot(2,1,1)
[cs,h] = contour(x,y,dzetadt);
set(gca,'fontsize',16)
clabel(cs,h), title('vorticity tendency')
xlabel('x (km)'), ylabel('y (km)')
subplot(2,1,2)
[cs,h] = contour(x,y,div);
clabel(cs,h)
hold on
quiver(x,y,uad,vad)
set(gca,'fontsize',16)
title('divergent ageostrophic wind and divergence')
xlabel('x (km)'), ylabel('y (km)');
%axis([-3000 3000 -1000 1000])


%return
figure(5);clf
subplot(2,1,1)
[cs,h] = contour(x,y,divag);
clabel(cs,h), title('divergence of ageostrophic wind')
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 its vorticity')
xlabel('x (km)'), ylabel('y (km)')





