clear
clf
load x.out
load y.out
load t.out
load v.out
load u.out
load w.out
load s.out
load s0.out
load m.out
load whov.out
a=size(whov);
n=round(a(1)/2);
w2hov=whov(n:a(1),:);
t2=t(n:a(1));
load vhov.out
v2hov=vhov(n:a(1),:);
load uhov.out
u2hov=uhov(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 Power spectra:  \n   9) u \n  10) v \n  11) w  \n \n  0) Exit \n ');
if it == 1
pcolor(x,y,s)
shading interp
hold
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')
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 == 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(['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(['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(['w (cm/s), from ',num2str(a),' to ',num2str(b)],'fontweight','bold')
elseif it == 9
	ut=u2hov(:,3);
   utf=fft(ut);
   a=abs(utf);
   amin=min(a);
   a(1)=amin;
   b=a+1-amin;
%   c=2.*log(b);
   c=log(b);
   c2=c(1:j);
   c2(1)=0.0;
   t2=tg(1:j);
   t2=t2-t2(1);
   h=plot(t2,c2);
   set(h,'linewidth',2)
   xlabel('Frequency (cycles/day)','fontweight','bold')
   ylabel('log(Power)','fontweight','bold')
   title('Zonal velocity','fontweight','bold')
   set(gca,'fontweight','bold')
   pause
   uk=u2hov(jm,:);
   ukf=fft(uk);
   a=abs(ukf);
   amin=min(a);
   a(1)=amin;
   b=a+1-amin;
%   c=2.*log(b);
   c=log(b);
   c2=c(1:im2);
   c2(1)=0.0;
   x2=xg(1:im2);
   x2=x2-2;
   h=plot(x2,c2);
   set(h,'linewidth',2)
   xlabel('Wavenumber','fontweight','bold')
   ylabel('log(Power)','fontweight','bold')
   title('Zonal velocity','fontweight','bold')
   set(gca,'fontweight','bold')
elseif it == 10
   ut=v2hov(:,3);
   utf=fft(ut);
   a=abs(utf);
   amin=min(a);
   a(1)=amin;
   b=a+1-amin;
%   c=2.*log(b);
   c=log(b);
   c2=c(1:j);
   c2(1)=0.0;
   t2=tg(1:j);
   t2=t2-t2(1);
   h=plot(t2,c2);
   set(h,'linewidth',2)
   xlabel('Frequency (cycles/day)','fontweight','bold')
   ylabel('log(Power)','fontweight','bold')
   title('Meridional velocity','fontweight','bold')
   set(gca,'fontweight','bold')
   pause
   uk=v2hov(jm,:);
   ukf=fft(uk);
   a=abs(ukf);
   amin=min(a);
   a(1)=amin;
   b=a+1-amin;
%   c=2.*log(b);
   c=log(b);
   c2=c(1:im2);
   c2(1)=0.0;
   x2=xg(1:im2);
   x2=x2-2;
   h=plot(x2,c2);
   set(h,'linewidth',2)
   xlabel('Wavenumber','fontweight','bold')
   ylabel('log(Power)','fontweight','bold')
   title('Meridional velocity','fontweight','bold')
   set(gca,'fontweight','bold')
elseif it == 11
   ut=w2hov(:,3);
   utf=fft(ut);
   a=abs(utf);
   amin=min(a);
   a(1)=amin;
   b=a+1-amin;
%   c=2.*log(b);
   c=log(b);
   c2=c(1:j);
   c2(1)=0.0;
   t2=tg(1:j);
   t2=t2-t2(1);
   h=plot(t2,c2);
   set(h,'linewidth',2)
   xlabel('Frequency (cycles/day)','fontweight','bold')
   ylabel('log(Power)','fontweight','bold')
   title('Vertical velocity','fontweight','bold')
   set(gca,'fontweight','bold')
   pause
   uk=w2hov(jm,:);
   ukf=fft(uk);
   a=abs(ukf);
   amin=min(a);
   a(1)=amin;
   b=a+1-amin;
   c=log(b);
   c2=c(1:im2);
   c2(1)=0.0;
   x2=xg(1:im2);
   x2=x2-2;
   h=plot(x2,c2);
   set(h,'linewidth',2)
   xlabel('Wavenumber','fontweight','bold')
   ylabel('log(Power)','fontweight','bold')
   title('Vertical velocity','fontweight','bold')
   set(gca,'fontweight','bold')
elseif it == 0
break
end
end
