clear all

dir1 = ...
'/Users/skrueger/ Teaching/ TEACHING 2008-09/Meteo 6150/QCOM code/gfortran_4_2_3 copy/'

fn1 = 'th.dat'

ny = 162
nz = 42
nt = 90

np = ny * nz

fn = [dir1 fn1 ]; 
fid1=fopen(fn);

aviobj = avifile('cp1a.avi','fps',6);

% Define grid

% k at T-levels is related to z by z = (K-1.5)*DZ
% j at T-points is related to y by y = (J-1.5)*DY

H = 10000.
YL = 40000.

DZ = H / ( nz - 2 )
DY = YL / ( ny - 2 )

J = 1:ny;
K = 1:nz;

y_km = (J-1.5)*DY / 1.e3;
z_km = (K-1.5)*DZ / 1.e3;

h1 = 4

cmax = 0-1.e-9;
cmin = -6;
cin = 0.5;
ci = [cmin : cin : cmax];

for t = 1:nt

th = fread(fid1,np,'float64','ieee-le');

th = reshape(th,ny,nz);

%pause(0.01)

figure(1)

% shaded contour plot

h=contourf(y_km,z_km,th',ci);

caxis([cmin-1.e-9 cmax+1.e-9])

colorbar('EastOutside')

axis([0,YL/1000,0,H/1000])

grid on

set(gca,'PlotBoxAspectRatio',[h1 1 1]) 

xlabel('Distance (km)')
ylabel('Height (km)')
title('Potential Temperature Perturbation (K)')

frame = getframe(gcf);
aviobj = addframe(aviobj,frame);

end

fclose(fid1)

viobj = close(aviobj);
