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);
i=1;
while i == 1
it=input('Plot choices: \n \n Final Time: \n   1) s, u & v \n   2) Mass flux \n   3) w \n   4) u   \n   5) s0 \n \n Hovmueller:  \n   6) Zonal velocity \n   7) Meridional velocity \n   8) Vertical velocity \n \n Wavenumber-Frequency spectra:  \n   9) u \n  10) v \n  11) w \n \n Asymmetric Wavenumber-Frequency spectra:  \n  12) u \n  13) v \n  14) w    \n \n  0) Exit \n ');
if it == 1
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 == 2
clf
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 == 3
clf   
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 == 4
clf
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 == 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')
elseif it == 55
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')
elseif it == 6
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')
pause
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')
elseif it == 7
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')
pause
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')
elseif it == 8
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')
pause
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')
elseif it == 9
   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 Zonal Wind at Equator','fontweight','bold')
   pause
   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 Zonal Wind at Equator - background','fontweight','bold')
elseif it == 10
   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 Meridional Wind at Equator','fontweight','bold')
   pause
   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 Meridional Wind at Equator - background','fontweight','bold')
elseif it == 11
   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
   set(gca,'fontweight','bold')
   xlabel('Zonal Wavenumber','fontweight','bold')
	ylabel('Frequency (day^{-1})','fontweight','bold')
   title('Log Power of Vertical Wind at Equator','fontweight','bold')
   pause
   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 Vertical Wind at Equator - background','fontweight','bold') 
elseif it == 12
   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')
   pause
   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') 
elseif it == 13
   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')
   pause
   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') 
elseif it == 14
   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
   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')
   pause
   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') 
elseif it == 0
break
end
end
