load output\x.out
load output\t.out
load output\p.out
sa=size(x);
ss=sa(1,2)-1;
s1=size(t);
ss1=s1(1,1);
fr=[1:ss1];
xt=x+0.5.*(x(2)-x(1));
x2=xt(1:ss);
i=1;
while i == 1
   r=input('\n Time-distance plots: \n \n 1)  Precipitation \n 2)  Omega \n 3)  v \n 4)  u \n 5)  Surface T \n 51) Surface SW \n \n Cross-sections: \n \n 6)  Omega \n 7)  M_u \n 8)  M_d \n 9)  Streamfunction \n 10) u \n 11) Relative humidity \n 12) Perturbation temperature \n 13) Angular momentum and streamfunction  \n \n Other quantities \n \n 14) Fluxes \n 15) Time-mean Surface T \n \n Spectra \n \n 16) Precipitation \n 17) u \n 18) Surface T \n 19) Surface SW \n \n 0) exit \n');
if r == 2
load output\omegat.out
b=max(max(omegat));
a=min(min(omegat));
c=omegat;
contourf(x,t,c,12)   
set(gca,'fontweight','bold','fontsize',12)
title(['Level 9 Omega (mb/hour) from ',num2str(a),' to ',num2str(b)])
xlabel('x (km)')
ylabel ('time (days)')
elseif r ==1
load output\precip.out
b=max(max(precip));
a=min(min(precip));
c=precip;
contourf(x,t,c,12)   
set(gca,'fontweight','bold','fontsize',12)
title(['Precipitation (mm/day) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('time (days)')
elseif r == 3
load output\vt.out
b=max(max(vt));
a=min(min(vt));
c2=vt(:,1:ss);
contour(x2,t,c2,12)   
set(gca,'fontweight','bold','fontsize',12)
title(['Surface v (m/s) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('time (days)')
elseif r == 4
load output\ut.out
b=max(max(ut));
a=min(min(ut));
c=ut;
contourf(x,t,c,12)   
set(gca,'fontweight','bold','fontsize',12)
title(['Surface u (m/s) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('time (days)')
elseif r == 5
load output\sst.out
b=max(max(sst));
a=min(min(sst));
c=sst;
set(gca,'fontweight','bold','fontsize',12)
contourf(x,t,c,12)   
set(gca,'fontweight','bold','fontsize',12)
title(['Surface temperature (C) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('time (days)')
elseif r == 51
load output\sw.out
b=max(max(sw));
a=min(min(sw));
c=sw;
contourf(x,t,c,12)   
set(gca,'fontweight','bold','fontsize',12)
title(['Surface absorbed SW (W/m^2) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('time (days)')
elseif r == 6
load output\omegal.out
b=max(max(omegal));
a=min(min(omegal));
contourf(x,-p,omegal,12)   
u=get(gca,'YTick');
u1=transpose(-u);
v=cellstr(num2str(u1));
set(gca, 'YTickLabel',v)
set(gca,'fontweight','bold','fontsize',12)
title(['Omega (mb/hour) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('Pressure (mb)')
elseif r == 7
load output\mul.out
b=max(max(mul));
a=min(min(mul));
contourf(x,-p,mul,12)   
u=get(gca,'YTick');
u1=transpose(-u);
v=cellstr(num2str(u1));
set(gca, 'YTickLabel',v)
set(gca,'fontweight','bold','fontsize',12)
title(['Updraft mass flux from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('Pressure (mb)')
elseif r == 8
load output\mdl.out
b=max(max(mdl));
a=min(min(mdl));
contourf(x,-p,mdl)   
u=get(gca,'YTick');
u1=transpose(-u);
v=cellstr(num2str(u1));
set(gca, 'YTickLabel',v)
set(gca,'fontweight','bold','fontsize',12)
title(['Downdraft mass flux from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('Pressure (mb)')
elseif r == 9
load output\psi.out
psi2=psi(:,1:ss);
b=max(max(psi2));
a=min(min(psi2));
contourf(x2,-p,psi2,12)   
u=get(gca,'YTick');
u1=transpose(-u);
v=cellstr(num2str(u1));
set(gca, 'YTickLabel',v)
set(gca,'fontweight','bold','fontsize',12)
title(['Streamfunction from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('Pressure (mb)')
elseif r == 10
load output\ul.out
b=max(max(ul));
a=min(min(ul));
contourf(x,-p,ul,12)   
u=get(gca,'YTick');
u1=transpose(-u);
v=cellstr(num2str(u1));
set(gca, 'YTickLabel',v)
set(gca,'fontweight','bold','fontsize',12)
title(['u (m/s) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('Pressure (mb)')
elseif r == 11
load output\rh.out
b=max(max(rh));
a=min(min(rh));
contourf(x,-p,rh)   
u=get(gca,'YTick');
u1=transpose(-u);
v=cellstr(num2str(u1));
set(gca, 'YTickLabel',v)
set(gca,'fontweight','bold','fontsize',12)
title(['Relative humidity from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('Pressure (mb)')
elseif r == 12
load output\temp.out
b=max(max(temp));
a=min(min(temp));
contourf(x,-p,temp,15)
u=get(gca,'YTick');
u1=transpose(-u);
v=cellstr(num2str(u1));
set(gca, 'YTickLabel',v)
set(gca,'fontweight','bold','fontsize',12)
title(['Perturbation temperature (C) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('Pressure (mb)')
elseif r == 13
load output\am.out
load output\psi.out
b=max(max(am));
a=min(min(am));
psi2=psi(:,1:ss);
contourf(x,-p,am,15)
hold on
contour(x2,-p,psi2,15,'k')
u=get(gca,'YTick');
u1=transpose(-u);
v=cellstr(num2str(u1));
set(gca, 'YTickLabel',v)
set(gca,'fontweight','bold','fontsize',12)
title(['Angular Momentum (m/s) from ',num2str(a),' to ',num2str(b)])
xlabel('y (km)')
ylabel ('Pressure (mb)')
hold off
elseif r == 14
load output\flux.out
aflux=flux(1:ss,2);
oflux=flux(1:ss,3);
nflux=flux(1:ss,4);
h=plot(x2,aflux,x2,oflux);   
set(h,'linewidth',3)
set(gca,'fontweight','bold','fontsize',12)
title('Energy Fluxes')
legend('Atmosphere','Ocean',0)
xlabel('y (km)','fontweight','bold','fontsize',12)
ylabel ('Flux (10^{15} W)','fontweight','bold','fontsize',12)
elseif r == 15
   load output\sst.out
   sstm=mean(sst,1);
h=plot(x,sstm);   
set(h,'linewidth',3)
set(gca,'fontweight','bold','fontsize',12)
title('Time-averaged Surface Temperature')
xlabel('y (km)','fontweight','bold','fontsize',12)
ylabel ('Surface Temperature (^{o}C)','fontweight','bold','fontsize',12)
elseif r == 16
   load output\precip.out
   var2=precip(:,2);
   pp=fft(var2);
   ppp=pp.*conj(pp)./ss1;
   h=plot(fr(4:25),ppp(4:25));
   set(h,'linewidth',3)
   set(gca,'fontweight','bold','fontsize',12)
   u=get(gca,'XTick');
   u1=transpose(1./u);
   u1=ss1.*0.01.*round(100.*u1);
   v=cellstr(num2str(u1));
   set(gca, 'XTickLabel',v);
   title('Precipitation, Sea')
   xlabel('Period (days)','fontweight','bold','fontsize',12)
   ylabel ('Power','fontweight','bold','fontsize',12)
   %
   pause
   %
   var2=precip(:,3);
   pp=fft(var2);
   ppp=pp.*conj(pp)./ss1;
   h=plot(fr(4:25),ppp(4:25));
   set(h,'linewidth',3)
   set(gca,'fontweight','bold','fontsize',12)
   set(gca, 'XTickLabel',v);
   title('Precipitation, Land')
   xlabel('Period (days)','fontweight','bold','fontsize',12)
   ylabel ('Power','fontweight','bold','fontsize',12)
elseif r == 17
   load output\ut.out
   var2=ut(:,2);
   pp=fft(var2);
   ppp=pp.*conj(pp)./ss1;
   h=plot(fr(4:25),ppp(4:25));
   set(h,'linewidth',3)
   set(gca,'fontweight','bold','fontsize',12)
   u=get(gca,'XTick');
   u1=transpose(1./u);
   u1=ss1.*0.01.*round(100.*u1);
   v=cellstr(num2str(u1));
   set(gca, 'XTickLabel',v);
   title('Zonal wind, Sea')
   xlabel('Period (days)','fontweight','bold','fontsize',12)
   ylabel ('Power','fontweight','bold','fontsize',12)
   %
   pause
   %
   var2=ut(:,3);
   pp=fft(var2);
   ppp=pp.*conj(pp)./ss1;
   h=plot(fr(4:25),ppp(4:25));
   set(h,'linewidth',3)
   set(gca,'fontweight','bold','fontsize',12)
   set(gca, 'XTickLabel',v);
   title('Zonal wind, Land')
   xlabel('Period (days)','fontweight','bold','fontsize',12)
   ylabel ('Power','fontweight','bold','fontsize',12)
elseif r == 18
   load output\sst.out
   var2=sst(:,2);
   pp=fft(var2);
   ppp=pp.*conj(pp)./ss1;
   h=plot(fr(4:25),ppp(4:25));
   set(h,'linewidth',3)
   set(gca,'fontweight','bold','fontsize',12)
   u=get(gca,'XTick');
   u1=transpose(1./u);
   u1=ss1.*0.01.*round(100.*u1);
   v=cellstr(num2str(u1));
   set(gca, 'XTickLabel',v);
   title('Surface temperature, Sea')
   xlabel('Period (days)','fontweight','bold','fontsize',12)
   ylabel ('Power','fontweight','bold','fontsize',12)
   %
   pause
   %
   var2=sst(:,3);
   pp=fft(var2);
   ppp=pp.*conj(pp)./ss1;
   h=plot(fr(4:25),ppp(4:25));
   set(h,'linewidth',3)
   set(gca,'fontweight','bold','fontsize',12)
   set(gca, 'XTickLabel',v);
   title('Surface temperature, Land')
   xlabel('Period (days)','fontweight','bold','fontsize',12)
   ylabel ('Power','fontweight','bold','fontsize',12)
elseif r == 19
   load output\sw.out
   var2=sw(:,2);
   pp=fft(var2);
   ppp=pp.*conj(pp)./ss1;
   h=plot(fr(4:25),ppp(4:25));
   set(h,'linewidth',3)
   set(gca,'fontweight','bold','fontsize',12)
   u=get(gca,'XTick');
   u1=transpose(1./u);
   u1=ss1.*0.01.*round(100.*u1);
   v=cellstr(num2str(u1));
   set(gca, 'XTickLabel',v);
   title('Surface absorbed shortwave, Sea')
   xlabel('Period (days)','fontweight','bold','fontsize',12)
   ylabel ('Power','fontweight','bold','fontsize',12)
   %
   pause
   %
   var2=sw(:,3);
   pp=fft(var2);
   ppp=pp.*conj(pp)./ss1;
   h=plot(fr(4:25),ppp(4:25));
   set(h,'linewidth',3)
   set(gca,'fontweight','bold','fontsize',12)
   set(gca, 'XTickLabel',v);
   title('Surface absorbed shortwave, Land')
   xlabel('Period (days)','fontweight','bold','fontsize',12)
   ylabel ('Power','fontweight','bold','fontsize',12)
elseif r == 0
break
end 
end
