%
% matlab source code to plot geopoential height data and
% use finite differences to obtain geostrophic winds
%
%.  Data created using extract_gemgrids.f
% 
%  This example uses traditional looping over the grid indices
%  to compute finite difference estimates for
%     ug = -(1/f) d[phi]/dy
%  ug(j) = (phi(j+1)-phi(j-1))/(2*f*dy)
%
% Note: that all matrices are stored as a(j,i) where j is latitude grid 
% and i is longitude grid
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Get things set up
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all

%the data directory
datadir='/home/disk/synoptic/class/atm442/matlab_data_feb03'; 

%your matlab code directory
codedir='/home/disk/synoptic/class/atm442/matlab_code';  

date= '2003020400';
lev= [500];

latr = [20 : 1.0 : 60];               %latitude range of data grid
lonr = [220 : 1.0 : 300];             %longitude range of data grid

delx = 1.0;
dely = 1.0;
J = size(latr,2);                % number of grid points in latitude
I = size(lonr,2);                % number of grid points in longitude

latv=zeros([J I]);                 % This step defines the size of the latv, 
lonv=zeros([J I]);                 % lonv arrays and fills them
                                   % with zeros
                                   % latv and lonv are used for
                                   % plotting vectors
for i=1:I
for j=1:J
  
  latv(j,i)=latr(j);
  lonv(j,i)=lonr(i);
end
end

% Read in data

cd (datadir)                           % Go to the data location
height = load(['HGHT',int2str(lev),'mb',date,'.txt']); 
phi=height*9.81 ;                      % phi is the geopotential

cd (codedir)                           % Go back to your code directory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First Plot Geopotential Height
%
% Typical contour limits and intervals for geopotential height
% at various pressure levels:
%
% clow=8040; chigh=9960; cint=120;        % 300 mb
% clow=4680; chigh=6000; cint=60;         % 500 mb
% clow=2460; chigh=3240; cint=60;         % 700 mb
% clow=1080; chigh=1620; cint=30;         % 850 mb
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1); clf                    

mapUS                             % make background map of US
clow=4680; chigh=6000; cint=60;   % contour line limits for 500 phi
conts = [clow:cint:chigh];        % define contour lines
V = [clow:cint*2:chigh];          % label every 2th line

pcolorm(height,[200,300]);  colormap(jet);   %puts in color fill
colorbar('v')
[c,h]=contorm(latr,lonr,height,conts,'linewidth',1,'color','k');
clabelm(c,h,V);                   %puts in contours
title('geopotential height at 500 hPa')

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute and plot geostrophic zonal wind using finite difference formula.
% You will copy this model for the meridional component of the wind.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ug=zeros([J I]);                 % set size of ug matrix
ug(1:J,1:I)=NaN;                 % set all values to NaN 

dy = 6.37e6*dely*pi/180;         % grid distance in meridional direction
cor= 2*7.292e-5*sin(latr*pi/180);  % Coriolis parameter
 for i=2:I-1
   for j=2:J-1
     ug(j,i) = - (phi(j+1,i)-phi(j-1,i))/(cor(j)*dy*2);
   end 
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
% Now plot the zonal geostrophic wind
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2);                   %plot the zonal wind component
mapUS
 
clow=-40; chigh=70; cint=5;  %contour lines limits and interval
conts = [clow:cint:chigh];
V = [clow:cint*2:chigh];     % label every 2nd 
[c,h]=contorm(latr,lonr,ug,conts,'linewidth',1.5,'color','r');
clabelm(c,h,V);

stations_4feb03              % this runs the stations_4feb03 script
                             % stars will be plotted at the station
                             % locations
title('geostrophic zonal wind at 500 hPa')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute meridional component of geostrophic wind here
% Don't forget that dx changes with latitude
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute the magnitude of the geostrophic wind here and magnitude
% of the observed wind
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Contouring two fields -- Use Height and ug for now
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3);                        %just contours for now
mapUS

clow=4680; chigh=6000; cint=60;   %contour line limits for 500 phi
conts = [clow:cint:chigh];        %define contour lines
V = [clow:cint*4:chigh];          % label every 4th line

[c,h]=contorm(latr,lonr,height,conts,'linewidth',1,'color','k');
clabelm(c,h,V);                   %puts in contours

clow=-40; chigh=60; cint=5;        %contour lines limits and interval
conts = [clow:cint:chigh];
V = [clow:cint*2:chigh];           % label every 2th
[c,h]=contorm(latr,lonr,ug,conts,'--b','linewidth',1);
clabelm(c,h,V);

title('geopotential height and mag ug at 500 hPa')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Plotting two vectors with two colors
% This example uses observed wind at 500 mb and at 850 mb
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lev1=850;
lev2=500
cd (datadir)
u1=load(['UREL',int2str(lev1),'mb',date,'.txt']);
u2=load(['UREL',int2str(lev2),'mb',date,'.txt']);
v1=load(['VREL',int2str(lev1),'mb',date,'.txt']);
v2=load(['VREL',int2str(lev2),'mb',date,'.txt']);


cd (codedir)
a=0.035                            %scale factor for vectors


figure(4)
mapUS

% these quiverm calls plot every other vector, suppresses the automatic 
%scaling and plot each vector in different colors.  g=green, b=blue

quiverm(latv(1:2:J,1:2:I),lonv(1:2:J,1:2:I),a*v1(1:2:J,1:2:I),a*u1(1:2:J,1:2:I),'g',0)
quiverm(latv(1:2:J,1:2:I),lonv(1:2:J,1:2:I),a*v2(1:2:J,1:2:I),a* ...
        u2(1:2:J,1:2:I),'b',0)

stations_4feb03

title('observed wind vectors 850mb and 500mb 98111000')