clear
clf
%
% ismooth smooths the wavenumber-frequency spectra...integer from 1 to 10
%
ismooth=1;
%
%  Activate this line to run program from within matlab
%
%!tropics3_ocean
%
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'))
xlab=1e-4.*x;
mplot=m;
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);
m=zeros(1,31);
for j=1:31,
   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)',...
      'Add colorbar to plot','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';
contourf(x,y,mplot,20);
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';
contourf(x,y,w,20);   
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));  
contourf(x,y,hum,20);
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:ismooth,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   contourf(x2,y2,plog3,20)
%   shading interp
   %vc=[5:15];
   %[c,h]=contour(x2,y2,plog3,vc,'k');
   set(gca,'fontweight','bold')
   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:ismooth,
      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   
   contourf(x2,y2,plog3,20)
%   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:ismooth,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   contourf(x2,y2,plog3,20)
   %pcolor(x2,y2,plog3)
   %shading interp
   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:ismooth,
      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,20)
   %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:ismooth,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   contourf(x2,y2,plog3,20)
   %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:ismooth,
      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)
   contourf(x2,y2,plog3,20)
   %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:ismooth,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   contourf(x2,y2,plog3,20)
   %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:ismooth,
      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); %#ok<*SAGROW>
   end   
%   plog4=mean(plog3,2);
%   for j=1:m2,
%     plog3(j,:)=plog3(j,:)-plog4(j);
%   end   
   %contourf(x2,y2,plog3)
   contourf(x2,y2,plog3,20)
   %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:ismooth,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   contourf(x2,y2,plog3,20)
   %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:ismooth,
      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)
   contourf(x2,y2,plog3,20)
   %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:ismooth,
      dummy=plog3;
      for k=2:m2-1,
         dummy(k,:)=0.25.*(plog3(k-1,:)+plog3(k+1,:)+2.*plog3(k,:));
      end
      plog3=dummy;
   end   
   contourf(x2,y2,plog3,20)
   %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:ismooth,
      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)
   contourf(x2,y2,plog3,20)
   %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
clear M
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'))
contourf(x,y,mmov01);
M(:,1)=getframe;
contourf(x,y,mmov02);
M(:,2)=getframe;
contourf(x,y,mmov03);
M(:,3)=getframe;
contourf(x,y,mmov04);
M(:,4)=getframe;
contourf(x,y,mmov05);
M(:,5)=getframe;
contourf(x,y,mmov06);
M(:,6)=getframe;
contourf(x,y,mmov07);
M(:,7)=getframe;
contourf(x,y,mmov08);
M(:,8)=getframe;
contourf(x,y,mmov09);
M(:,9)=getframe;
contourf(x,y,mmov10);
M(:,10)=getframe;
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
clear M
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'))
contourf(x,y,umov01);
M(:,1)=getframe;
contourf(x,y,umov02);
M(:,2)=getframe;
contourf(x,y,umov03);
M(:,3)=getframe;
contourf(x,y,umov04);
M(:,4)=getframe;
contourf(x,y,umov05);
M(:,5)=getframe;
contourf(x,y,umov06);
M(:,6)=getframe;
contourf(x,y,umov07);
M(:,7)=getframe;
contourf(x,y,umov08);
M(:,8)=getframe;
contourf(x,y,umov09);
M(:,9)=getframe;
contourf(x,y,umov10);
M(:,10)=getframe;
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
clear M
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'))
contourf(x,y,vmov01);
M(:,1)=getframe;
contourf(x,y,vmov02);
M(:,2)=getframe;
contourf(x,y,vmov03);
M(:,3)=getframe;
contourf(x,y,vmov04);
M(:,4)=getframe;
contourf(x,y,vmov05);
M(:,5)=getframe;
contourf(x,y,vmov06);
M(:,6)=getframe;
contourf(x,y,vmov07);
M(:,7)=getframe;
contourf(x,y,vmov08);
M(:,8)=getframe;
contourf(x,y,vmov09);
M(:,9)=getframe;
contourf(x,y,vmov10);
M(:,10)=getframe;
contourf(x,y,vmov11);
M(:,11)=getframe;
movie(M,11)
set(gca,'fontweight','bold')
xlabel('X (km)','fontweight','bold')
ylabel('Y (km)','fontweight','bold')
elseif it == m(28)
clf
clear M
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'))
contourf(x,y,tmov01);
hold on
quiver(x,y,umov01,vmov01)
hold off
M(:,1)=getframe;
contourf(x,y,tmov02);
hold on
quiver(x,y,umov02,vmov02)
hold off
M(:,2)=getframe;
contourf(x,y,tmov03);
hold on
quiver(x,y,umov03,vmov03)
hold off
M(:,3)=getframe;
contourf(x,y,tmov04);
hold on
quiver(x,y,umov04,vmov04)
hold off
M(:,4)=getframe;
contourf(x,y,tmov05);
hold on
quiver(x,y,umov05,vmov05)
hold off
M(:,5)=getframe;
contourf(x,y,tmov06);
hold on
quiver(x,y,umov06,vmov06)
hold off
M(:,6)=getframe;
contourf(x,y,tmov07);
hold on
quiver(x,y,umov07,vmov07)
hold off
M(:,7)=getframe;
contourf(x,y,tmov08);
hold on
quiver(x,y,umov08,vmov08)
hold off
M(:,8)=getframe;
contourf(x,y,tmov09);
hold on
quiver(x,y,umov09,vmov09)
hold off
M(:,9)=getframe;
contourf(x,y,tmov10);
hold on
quiver(x,y,umov10,vmov10)
hold off
M(:,10)=getframe;
contourf(x,y,tmov11);
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)
   cbb=colorbar('vert');
   set(cbb,'fontweight','bold')
elseif it == m(30)
   fig=strcat(name,'.png');
   exportfig(gcf,fig,'format','png','color','rgb','bounds','loose','Renderer','Painters','Resolution',300);
elseif it == m(31)   
   break
end
end
