%BRENT LEITHAUSER  29 June 2004
 % the variables are 
    %sat(365 days, 73 latitudes, 144 longitudes) surface air temperature, degrees Kelvin
    %z1000  dynamic height in meters near the Earth's surface
    %      at 1000 HPa (hectoPascals ~ millibars)
    %      this is essentially the same as the atmospheric pressure
    %      variation at a fixed altitude
    %z300   dynamic height at 300 HPa, or the top of the troposphere where
    %       the jetstream is strongest
    %z30    dyn height at 30 HPa in the stratosphere.
    %    For all these the units are meters (slightly adjusted height in
    %    meters above mean sea-level.
    % The data is on a grid of 73 latitude points (-90- to 90 latitude)
    % by 144 longitude points (0-360 East longitude).
    % time resolution is 1 per day
    % 
    %  Note, plots of the height of a surface of constant pressure, Z...,
    %  are almost the same info as plots of pressure, P,  at a fixed alitude
    %  This is because hydrostatic pressure balance holds in the vertical
    %  momentum equation, so that the pressure at a given point is simply
    %  equal to the weight of the air overhead.  This eventually connects Z
    %  fields with P fields, through dP/dz = - ro g  where ro is the
    %  density and g is gravitational acceleration.
    
    
 load dh-2002.mat
 plot(lat,squeeze(sat(1,:,:)))   % plot the 1 January surface air temp against
                                    % latitude, for all longitudes
 contour(lon,lat,squeeze(sat(1,:,:)),25)  % make a contour map of same
 pcolor(lon,lat,squeeze(sat(1,:,:))), shading interp, colorbar
            %  shaded contour map with colorbar
 hold on    % retain current plot
 contour(lon,lat,squeeze(sat(1,:,:)),25,'k')  % overlay contours on color map
                                    % using 25 black contours 
 figure       % new figure
 contour(lon,lat,squeeze(z300(1,:,:)),25)  % contour map of isobars (hence flow lines
                                        % at the top of the troposphere
 plot(lat,squeeze(z300(1,:,:)));grid    % to give an idea of how high this level is
                                        % in meters (between 8.4 km and 9.8
                                        % km
 mesh(lon,lat,squeeze(z300(1,:,:)))      % wire-frame plot, which you can rotate
                                        % by clicking the right-most button
                                        % on the top of the figure frame,
                                        % and then mousing the cursor over
                                        % the figure
                                        
 %Try to get an idea of the shapes of these pressure fields from south to
 %north, as a function of season.
           % To get an idea of the vertical structure of the atmosphere try
 figure
 plot(lat,squeeze(z1000(1,:,:)))
 hold on
 plot(lat,squeeze(z300(1,:,:)))
 plot(lat,squeeze(z30(1,:,:))) 
 %  and then work on ways to show these constant pressure surfaces more
 %  clearly.  
 %  Try to explore the map-view of the atmospheric circulation, and even
 %  the polar vortex in the winter stratosphere (the z30 field).
 %      Time averages of the circulation are useful, like the wintertime
 %      average flow at the jetstream level:
 z300winterav= mean(z300(1:100,:,:));
 contour(lon,lat,squeeze(z300winterav),25)
 
 
                     % ANIMATE NORTH_SOUTH CUTS OF Z
 itmax = 10
 for it = 1:itmax
     zz=squeeze(z30(it,:,1));
     plot(lat,zz)
     axis([-90 90 2.2e4 2.45e4])
     grid
     drawnow
 end
 % Smoother if capture frames with 'movie' statement:                                    
 for it = 1:itmax
     zz=squeeze(z30(it,:,1));
     plot(lat,zz)
     axis([-90 90 2.08e4 2.45e4])
     grid
      day=it;
      text(0,2.17e4,['2002 Z30 ' num2str(fix(day))]);
      xlabel('lat'); ylabel('meters')
      drawnow
     M(it)=getframe;
 end
 % play it
     movie(M,3,10)
 % Notice the seasonal shift of mass north and south: in winter
 % radiation cools the upper atmosphere particularly, and the dense air 'slumps'
 % downward, with air converging on the Pole above, diverging away from the
 % Pole below. Angular momentum conservation makes the converging air spin
 % cyclonically (a low pressure center develops as the 'top of the
 % atmosphere' is sucked down ). Thus, ironically, cooling the polar
 % regions of the Earth makes the pressure decrease!
 %
 % ANIMATE ALL 3 LEVELS: LOOK FOR HIGH AND LOW PRESSURE 'BELTS'
   itmax=365
  for it = 1:itmax
     subplot(311)
     zz=squeeze(z30(it,:,1));
     plot(lat,zz)
     axis([-90 90 2.08e4 2.45e4])
     grid
     day=it;
     text(0,2.17e4,['2002 Z30 ' num2str(fix(day))]);
     xlabel('lat'); ylabel('meters')
   subplot(312)
     zz=squeeze(z300(it,:,1));
     plot(lat,zz)
     axis([-90 90 8000 10000])
     grid
   subplot(313)
     zz=squeeze(z1000(it,:,1));
     plot(lat,zz)
     axis([-90 90 -500 500])
     grid
      drawnow
   
 end
 
 
  m_proj('satellite', 'lon',300,'lat', 65)
  m_grid('xaxislocation','middle');
  m_coast('color','b','linewidth',1);