clear
clf
file=input('Enter output file (default = output)','s');
if isempty(file)
   file='output';
end   
load(strcat(file,'/x.out'))
load(strcat(file,'/y.out'))
load(strcat(file,'/t.out'))
load(strcat(file,'/v.out'))
load(strcat(file,'/u.out'))
load(strcat(file,'/w.out'))
load(strcat(file,'/s.out'))
load(strcat(file,'/s0.out'))
load(strcat(file,'/sm.out'))
load(strcat(file,'/m.out'))
load(strcat(file,'/uhov.out'))
load(strcat(file,'/uhova.out'))
load(strcat(file,'/vhov.out'))
load(strcat(file,'/vhova.out'))
load(strcat(file,'/whov.out'))
load(strcat(file,'/whova.out'))
a=size(whov);
n=round(a(1)/2);
w2hov=whov(n:a(1),:);
w2hova=whova(n:a(1),:);
t2=t(n:a(1));
v2hov=vhov(n:a(1),:);
v2hova=vhova(n:a(1),:);
u2hov=uhov(n:a(1),:);
u2hova=uhova(n:a(1),:);
xs=size(x);
xb=x(2:xs(2)-1);
xs=size(xb);
xg=round(xb*xs(2)/xb(xs(2)));
im=round(xs(2));
im2=round(im/2);
xs2=xs/2;
ts=size(t2);
jm=round(ts(1));
ts2=ts(1)/2;
j=round(ts2(1));
ftmax=ts/t2(ts(1));
%ftmax=t*ts(1)./t(ts(1));
tg=4.*t2*ftmax(1)/t2(jm);
for j=1:30,
   m(j)=j;
end   
i=1;
while i == 1
   it=menu('Select from the following:','s, u and v (final time)',...
      'Mass flux (final time)','vertical velocity (final time)','zonal velocity (final time',...
      's_0 (final time)','relative humidity (final time)',...
      'Symmetric zonal velocity (hovmuller)','Asymmetric zonal velocity (hovmuller)',...
      'Symmetric meridional velocity (hovmuller)','Asymmetric meridional velocity (hovmuller)',...
      'Symmetric vertical velocity (hovmuller)','Asymmetric vertical velocity (hovmuller)',...
      'Symmetric zonal velocity (wavenumber-frequency spectrum)','Symmetric zonal velocity - background (wavenumber-frequency spectrum)',...
      'Symmetric meridional velocity (wavenumber-frequency spectrum)','Symmetric meridional velocity - background (wavenumber-frequency spectrum)',...
      'Symmetric vertical velocity (wavenumber-frequency spectrum)','Symmetric vertical velocity - background (wavenumber-frequency spectrum)',...
      'Asymmetric zonal velocity (wavenumber-frequency spectrum)','Asymmetric zonal velocity - background (wavenumber-frequency spectrum)',...
      'Asymmetric meridional velocity (wavenumber-frequency spectrum)','Asymmetric meridional velocity - background (wavenumber-frequency spectrum)',...
      'Asymmetric vertical velocity (wavenumber-frequency spectrum)','Asymmetric vertical velocity - background (wavenumber-frequency spectrum)',...
      'Mass flux (movie)','Zonal velocity (movie)','Merdional velocity (movie)','Entropy and surface wind (movie)',...
      'Make .png plot','Exit');
if it == m(1)
name='theta_e';     
pcolor(x,y,s)
shading interp
hold on
quiver(x,y,u,v)
xlabel('x (km)','fontweight','bold')
ylabel('y (km)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(s));
a=min(min(s));
title(['/theta_e, from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
hold off
elseif it == m(2)
   clf
   name='Mass_flux';
pcolor(x,y,m)
shading interp
xlabel('x (km)','fontweight','bold')
ylabel('y (km)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(m));
a=min(min(m));
title(['Mass flux (cm/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
elseif it == m(3)
   clf   
   name='w';
pcolor(x,y,w)
shading interp
xlabel('x (km)','fontweight','bold')
ylabel('y (km)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(w));
a=min(min(w));
title(['w (cm/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
elseif it == m(4)
   clf
   name='u';
pcolor(x,y,u)
shading interp
xlabel('x (km)','fontweight','bold')
ylabel('y (km)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(u));
a=min(min(u));
title(['u (m/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
elseif it == m(5)
clf
pcolor(x,y,s0)
shading interp
xlabel('x (km)','fontweight','bold')
ylabel('y (km)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(s0));
a=min(min(s0));
title(['/theta_{e0} (K), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
name='theta_e0';
elseif it == m(6)
clf
hum=100.*(1.-11.*log(s./sm));  
pcolor(x,y,hum)
shading interp
xlabel('x (km)','fontweight','bold')
ylabel('y (km)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(hum));
a=min(min(hum));
title(['Relative Humidity (%), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
name='humidity';
elseif it == m(7)
clf
contourf(x,t2,u2hov)
xlabel('x (km)','fontweight','bold')
ylabel('Time (days)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(u2hov));
a=min(min(u2hov));
title(['Symmetric u (m/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
name='u_sym_hov';
elseif it == m(8)
clf
contourf(x,t2,u2hova)
xlabel('x (km)','fontweight','bold')
ylabel('Time (days)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(u2hova));
a=min(min(u2hova));
title(['Asymmetric u (m/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
name='u_asym_hov';
elseif it == m(9)
clf
contourf(x,t2,v2hov)
xlabel('x (km)','fontweight','bold')
ylabel('Time (days)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(v2hov));
a=min(min(v2hov));
title(['Symmetric v (m/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
name='v_sym_hov';
elseif it == m(10)
clf
contourf(x,t2,v2hova)
xlabel('x (km)','fontweight','bold')
ylabel('Time (days)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(v2hova));
a=min(min(v2hova));
title(['Asymmetric v (m/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
name='v_asym_hov';
elseif it == m(11)
clf
contourf(x,t2,w2hov)
xlabel('x (km)','fontweight','bold')
ylabel('Time (days)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(w2hov));
a=min(min(w2hov));
title(['Symmetric w (cm/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
name='w_sym_hov';
elseif it == m(12)
clf
contourf(x,t2,w2hova)
xlabel('x (km)','fontweight','bold')
ylabel('Time (days)','fontweight','bold')
set(gca,'fontweight','bold')
b=max(max(w2hova));
a=min(min(w2hova));
title(['Asymmetric w (cm/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
name='w_asym_hov';
elseif it == m(13)
   fu=fft2(uhov);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
   ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Symmetric Zonal Wind at Equator','fontweight','bold')
   name='u_sym_power';
elseif it == m(14)
   fu=fft2(uhov);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   plog4=mean(plog3,1);
   for j=1:n-1,
     plog3(:,j)=plog3(:,j)-plog4(j);
   end   
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Symmetric Zonal Wind at Equator - background','fontweight','bold')
   name='u_sym_power-back';
elseif it == m(15)
   fu=fft2(vhov);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Symmetric Meridional Wind at Equator','fontweight','bold')
   name='v_sym_power';
elseif it == m(16)
   fu=fft2(vhov);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   plog4=mean(plog3,1);
   for j=1:n-1,
     plog3(:,j)=plog3(:,j)-plog4(j);
   end   
   %plog4=mean(plog3,2);
   %for j=1:m2,
   %  plog3(j,:)=plog3(j,:)-plog4(j);
   %end   
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Symmetric Meridional Wind at Equator - background','fontweight','bold')
   name='v_sym_power-back';
elseif it == m(17)
   fu=fft2(whov);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   pcolor(x2,y2,plog3)
   shading interp
   %plog3=smooth2(plog3);
   %plog3=smooth2(plog3);
   %contourf(x2,y2,plog3,20)
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Symmetric Vertical Wind at Equator','fontweight','bold')
   name='w_sym_power';
elseif it == m(18)
   fu=fft2(whov);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   plog4=mean(plog3,1);
   for j=1:n-1,
     plog3(:,j)=plog3(:,j)-plog4(j);
   end   
%   plog4=mean(plog3,2);
%   for j=1:m2,
%     plog3(j,:)=plog3(j,:)-plog4(j);
%   end   
   %contourf(x2,y2,plog3)
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Symmetric Vertical Wind at Equator - background','fontweight','bold') 
   name='w_sym_power-back';
elseif it == m(19)
   fu=fft2(uhova);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Asymmetric Zonal Wind at Equator','fontweight','bold')
   name='u_asym_power';
elseif it == m(20)
   fu=fft2(uhova);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   plog4=mean(plog3,1);
   for j=1:n-1,
     plog3(:,j)=plog3(:,j)-plog4(j);
   end   
%   plog4=mean(plog3,2);
%   for j=1:m2,
%     plog3(j,:)=plog3(j,:)-plog4(j);
%   end   
   %contourf(x2,y2,plog3)
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Asymmetric Zonal Wind at Equator - background','fontweight','bold') 
   name='u_asym_power-back';
elseif it == m(21)
   fu=fft2(vhova);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Asymmetric Meridional Wind at Equator','fontweight','bold')
   name='v_asym_power';
elseif it == m(22)
   fu=fft2(vhova);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   plog4=mean(plog3,1);
   for j=1:n-1,
     plog3(:,j)=plog3(:,j)-plog4(j);
   end   
%   plog4=mean(plog3,2);
%   for j=1:m2,
%     plog3(j,:)=plog3(j,:)-plog4(j);
%   end   
   %contourf(x2,y2,plog3)
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Asymmetric Meridional Wind at Equator - background','fontweight','bold') 
   name='v_asym_power-back';
elseif it == m(23)
   fu=fft2(whova);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   pcolor(x2,y2,plog3)
   shading interp
   %plog3=smooth2(plog3);
   %plog3=smooth2(plog3);
   %contourf(x2,y2,plog3,20)
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Asymmetric Vertical Wind at Equator','fontweight','bold')
   name='w_asym_power';
elseif it == m(24)
   fu=fft2(whova);
   power=abs(fu).^2;
   power(1,1)=0;
   plog=log(power+1);
   [ma,n]=size(plog);
   m1=2.*floor(ma/2);
   n1=2.*floor(n/2);
   m2=m1/2;
   n2=n1/2;
   plog2(1:m1,n1:-1:1)=plog(1:m1,1:n1);
   plog3(1:m2,n2+1:n1)=plog2(1:m2,1:n2);
   plog3(1:m2,1:n2)=plog2(1:m2,n2+1:n1);
   x2=1:n1;
   x2=x2-n2;
   y2=0.25.*(1:m2)./(m2.*(t2(2)-t2(1)));
   x2=transpose(x2);
   y2=transpose(y2);
   for j=1:10,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   plog4=mean(plog3,1);
   for j=1:n-1,
     plog3(:,j)=plog3(:,j)-plog4(j);
   end   
%   plog4=mean(plog3,2);
%   for j=1:m2,
%     plog3(j,:)=plog3(j,:)-plog4(j);
%   end   
   %contourf(x2,y2,plog3)
   pcolor(x2,y2,plog3)
   shading interp
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Asymmetric Vertical Wind at Equator - background','fontweight','bold') 
   name='w_asym_power-back';
elseif it == m(25)
clf
load(strcat(file,'/x.out'))
load(strcat(file,'/y.out'))
load(strcat(file,'/mmov01'))
load(strcat(file,'/mmov02'))
load(strcat(file,'/mmov03'))
load(strcat(file,'/mmov04'))
load(strcat(file,'/mmov05'))
load(strcat(file,'/mmov06'))
load(strcat(file,'/mmov07'))
load(strcat(file,'/mmov08'))
load(strcat(file,'/mmov09'))
load(strcat(file,'/mmov10'))
load(strcat(file,'/mmov11'))
h=contourf(x,y,mmov01);
M(:,1)=getframe;
h=contourf(x,y,mmov02);
M(:,2)=getframe;
h=contourf(x,y,mmov03);
M(:,3)=getframe;
h=contourf(x,y,mmov04);
M(:,4)=getframe;
h=contourf(x,y,mmov05);
M(:,5)=getframe;
h=contourf(x,y,mmov06);
M(:,6)=getframe;
h=contourf(x,y,mmov07);
M(:,7)=getframe;
h=contourf(x,y,mmov08);
M(:,8)=getframe;
h=contourf(x,y,mmov09);
M(:,9)=getframe;
h=contourf(x,y,mmov10);
M(:,10)=getframe;
h=contourf(x,y,mmov11);
M(:,11)=getframe;
movie(M,5)
set(gca,'fontweight','bold')
xlabel('X (km)','fontweight','bold')
ylabel('Y (km)','fontweight','bold')
elseif it == m(26)
clf
load(strcat(file,'/x.out'))
load(strcat(file,'/y.out'))
load(strcat(file,'/umov01'))
load(strcat(file,'/umov02'))
load(strcat(file,'/umov03'))
load(strcat(file,'/umov04'))
load(strcat(file,'/umov05'))
load(strcat(file,'/umov06'))
load(strcat(file,'/umov07'))
load(strcat(file,'/umov08'))
load(strcat(file,'/umov09'))
load(strcat(file,'/umov10'))
load(strcat(file,'/umov11'))
h=contourf(x,y,umov01);
M(:,1)=getframe;
h=contourf(x,y,umov02);
M(:,2)=getframe;
h=contourf(x,y,umov03);
M(:,3)=getframe;
h=contourf(x,y,umov04);
M(:,4)=getframe;
h=contourf(x,y,umov05);
M(:,5)=getframe;
h=contourf(x,y,umov06);
M(:,6)=getframe;
h=contourf(x,y,umov07);
M(:,7)=getframe;
h=contourf(x,y,umov08);
M(:,8)=getframe;
h=contourf(x,y,umov09);
M(:,9)=getframe;
h=contourf(x,y,umov10);
M(:,10)=getframe;
h=contourf(x,y,umov11);
M(:,11)=getframe;
movie(M,5)
set(gca,'fontweight','bold')
xlabel('X (km)','fontweight','bold')
ylabel('Y (km)','fontweight','bold')
elseif it == m(27)
clf
load(strcat(file,'/x.out'))
load(strcat(file,'/y.out'))
load(strcat(file,'/vmov01'))
load(strcat(file,'/vmov02'))
load(strcat(file,'/vmov03'))
load(strcat(file,'/vmov04'))
load(strcat(file,'/vmov05'))
load(strcat(file,'/vmov06'))
load(strcat(file,'/vmov07'))
load(strcat(file,'/vmov08'))
load(strcat(file,'/vmov09'))
load(strcat(file,'/vmov10'))
load(strcat(file,'/vmov11'))
h=contourf(x,y,vmov01);
M(:,1)=getframe;
h=contourf(x,y,vmov02);
M(:,2)=getframe;
h=contourf(x,y,vmov03);
M(:,3)=getframe;
h=contourf(x,y,vmov04);
M(:,4)=getframe;
h=contourf(x,y,vmov05);
M(:,5)=getframe;
h=contourf(x,y,vmov06);
M(:,6)=getframe;
h=contourf(x,y,vmov07);
M(:,7)=getframe;
h=contourf(x,y,vmov08);
M(:,8)=getframe;
h=contourf(x,y,vmov09);
M(:,9)=getframe;
h=contourf(x,y,vmov10);
M(:,10)=getframe;
h=contourf(x,y,vmov11);
M(:,5)=getframe;
movie(M,11)
set(gca,'fontweight','bold')
xlabel('X (km)','fontweight','bold')
ylabel('Y (km)','fontweight','bold')
elseif it == m(28)
clf
load(strcat(file,'/x.out'))
load(strcat(file,'/y.out'))
load(strcat(file,'/umov01'))
load(strcat(file,'/umov02'))
load(strcat(file,'/umov03'))
load(strcat(file,'/umov04'))
load(strcat(file,'/umov05'))
load(strcat(file,'/umov06'))
load(strcat(file,'/umov07'))
load(strcat(file,'/umov08'))
load(strcat(file,'/umov09'))
load(strcat(file,'/umov10'))
load(strcat(file,'/umov11'))
load(strcat(file,'/vmov01'))
load(strcat(file,'/vmov02'))
load(strcat(file,'/vmov03'))
load(strcat(file,'/vmov04'))
load(strcat(file,'/vmov05'))
load(strcat(file,'/vmov06'))
load(strcat(file,'/vmov07'))
load(strcat(file,'/vmov08'))
load(strcat(file,'/vmov09'))
load(strcat(file,'/vmov10'))
load(strcat(file,'/vmov11'))
load(strcat(file,'/tmov01'))
load(strcat(file,'/tmov02'))
load(strcat(file,'/tmov03'))
load(strcat(file,'/tmov04'))
load(strcat(file,'/tmov05'))
load(strcat(file,'/tmov06'))
load(strcat(file,'/tmov07'))
load(strcat(file,'/tmov08'))
load(strcat(file,'/tmov09'))
load(strcat(file,'/tmov10'))
load(strcat(file,'/tmov11'))
h=pcolor(x,y,tmov01);
shading interp
hold on
quiver(x,y,umov01,vmov01)
hold off
M(:,1)=getframe;
h=pcolor(x,y,tmov02);
shading interp
hold on
quiver(x,y,umov02,vmov02)
hold off
M(:,2)=getframe;
h=pcolor(x,y,tmov03);
shading interp
hold on
quiver(x,y,umov03,vmov03)
hold off
M(:,3)=getframe;
h=pcolor(x,y,tmov04);
shading interp
hold on
quiver(x,y,umov04,vmov04)
hold off
M(:,4)=getframe;
h=pcolor(x,y,tmov05);
shading interp
hold on
quiver(x,y,umov05,vmov05)
hold off
M(:,5)=getframe;
h=pcolor(x,y,tmov06);
shading interp
hold on
quiver(x,y,umov06,vmov06)
hold off
M(:,6)=getframe;
h=pcolor(x,y,tmov07);
shading interp
hold on
quiver(x,y,umov07,vmov07)
hold off
M(:,7)=getframe;
h=pcolor(x,y,tmov08);
shading interp
hold on
quiver(x,y,umov08,vmov08)
hold off
M(:,8)=getframe;
h=pcolor(x,y,tmov09);
shading interp
hold on
quiver(x,y,umov09,vmov09)
hold off
M(:,9)=getframe;
h=pcolor(x,y,tmov10);
shading interp
hold on
quiver(x,y,umov10,vmov10)
hold off
M(:,10)=getframe;
h=pcolor(x,y,tmov11);
shading interp
hold on
quiver(x,y,umov11,vmov11)
hold off
M(:,11)=getframe;
movie(M,5,3)
set(gca,'fontweight','bold')
xlabel('X (km)','fontweight','bold')
ylabel('Y (km)','fontweight','bold')
elseif it == m(29)
   fig=strcat(name,'.png');
   exportfig(gcf,fig,'format','png','color','rgb','bounds','loose','Renderer','Painters','Resolution',300);
elseif it == m(30)   
   break
end
end