figure(1), clf,

subplot(2,2,1)
[h]=plot(Time/YEAR,Ts1,'r');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Temperature (K)','Fontsize',15);
title('Tropical Ts','Fontsize',15);

subplot(2,2,2)
[h]=plot(Time/YEAR,Ts2,'b');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Temperature (K)','Fontsize',15);
title('Extra-Tropical Ts','Fontsize',15);

subplot(2,2,3)
[h]=plot(Time/YEAR,Ts,'g');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Temperature (K)','Fontsize',15);
title('Global Ts','Fontsize',15);

subplot(2,2,4)
[h]=plot(Time/YEAR,Ts1-Ts2,'c');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Temperature (K)','Fontsize',15);
title('Equator-Pole \Delta Ts','Fontsize',15);

figure(2), clf,

subplot(2,2,1)
[h]=plot(Time(1:end-1)/YEAR,Ft,'g');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Energy flux (Wm-2)','Fontsize',15);
title('Global TOA net radiation','Fontsize',15);

subplot(2,2,2)
[h]=plot(Time(1:end-1)/YEAR,Ft1,'r',Time(1:end-1)/YEAR,Ft2,'b');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Energy flux (Wm-2)','Fontsize',15);
title('TOA net radiation (Box 1 & 2)','Fontsize',15);

subplot(2,2,3)
[h]=plot(Time(1:end-1)/YEAR,Fs1,'r',Time(1:end-1)/YEAR,Fs2,'b');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Energy flux (Wm-2)','Fontsize',15);
title('Net surface heat flux(box 1 & 2)','Fontsize',15);

subplot(2,2,4)
[h]=plot(Time(1:end-1)/YEAR,Feva1,'r',Time(1:end-1)/YEAR,Feva2,'b');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Energy flux (Wm-2)','Fontsize',15);
title('Evaporative cooling (box 1 & 2)','Fontsize',15);

figure(3), clf,

subplot(2,2,1)
[h]=plot(Time(1:Nt-1)/YEAR,(Fa*pi*RADIUS^2)/PW,'m'); hold on
set(h,'Linewidth',2);
[h]=plot(Time(1:Nt-1)/YEAR,(Fo*pi*RADIUS^2)/PW,'c'); hold on
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Energy transport(PW)','Fontsize',15);
title('Heat transport (O & A)','Fontsize',15);

subplot(2,2,2)
[h]=plot(Time(1:Nt-1)/YEAR,((Fa+Fo)*pi*RADIUS^2)/PW,'k'); hold on
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Energy transport(PW)','Fontsize',15);
title('Total heat transport (PW)','Fontsize',15);

subplot(2,2,3)
[h]=plot(Time(1:Nt-1)/YEAR,PSIa/SV,'k');
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Circulation (10^9 kg/s)','Fontsize',15);
title('Atmospheric circulation strength','Fontsize',15);

figure(4), clf,

[h]=plot(Time/YEAR,CO2,'k'); hold on
set(h,'Linewidth',2);
set(gca,'Fontsize',15);
xlabel('Time (year)','Fontsize',15);
ylabel('Concentration (ppm)','Fontsize',15);
title('Imposed CO2','Fontsize',15);

return

figure(3), clf,
for n = 1:200:Nt-1
    subplot(2,3,1)
    [h]=plot(Ta1(n),To1(n),'r.'); hold on;set(h,'markersize',10); 
    [h]=plot(Ta2(n),To2(n),'b.'); set(h,'markersize',10);
    %[h]=plot((Ta1(n)+Ta2(n))/2,(To1(n)+To2(n))/2,'g.'); set(h,'markersize',10);
    plot([Ta1_end Ta1_end],[0 400],'r');plot([0 400],[To1_end To1_end],'r');
    plot([Ta2_end Ta2_end],[0 400],'b');plot([0 400],[To2_end To2_end],'b');
    set(gca,'Fontsize',10);
    axis([225 325 225 325]);
    xlabel('TA (K)','Fontsize',15);
    ylabel('TO (K)','Fontsize',15);
    pause(1/5);
    
    subplot(2,3,2)
    [h]=plot(Fa(n)*pi*(RADIUS^2)/PW,Fo(n)*pi*(RADIUS^2)/PW,'m.'); 
    hold on;set(h,'markersize',10); 
    plot([Fa_end Fa_end],[-400 400],'m');plot([-400 400],[Fo_end Fo_end],'m');
    H=5; axis([-H H -H H]);
    xlabel('HA (PW)','Fontsize',15);
    ylabel('HO (PW)','Fontsize',15);
    pause(1/5);
    
    subplot(2,3,3)
    [h]=plot(Feva1(n),Feva2(n),'c.'); 
    hold on;set(h,'markersize',10); 
    H=200; axis([0 H 0 H]);
    xlabel('LvE1 (W m-2)','Fontsize',15);
    ylabel('LvE2 (W m-2)','Fontsize',15);
    pause(1/5);
    
    subplot(2,3,4)
    [h]=plot(Ts1(n),Ts2(n),'g.'); 
    hold on;set(h,'markersize',10);
    plot([200:1:400],[200:1:400]-DTcrit_circ,'g--');
    plot([Ts1_end Ts1_end],[0 400],'g');plot([0 400],[Ts2_end Ts2_end],'g');
    axis([240 320 240 320]);
    xlabel('Ts1 (K)','Fontsize',15);
    ylabel('Ts2 (K)','Fontsize',15);
    pause(1/5);
    
    subplot(2,3,5)
    [h]=plot(Ts(n),Ft(n),'k.'); 
    hold on;set(h,'markersize',10);
    plot([Ts_end Ts_end],[-500 500],'k');plot([250 400],[Ft_end Ft_end],'k');
    axis([250 325 -200 200]);
    xlabel('Ts (K)','Fontsize',15);
    ylabel('Ft (Wm-2)','Fontsize',15);
    pause(1/5);
end