%
% Run Two-Box Model
%
for n = 1:Nt-1
    
   'Year = ...', n*Dt/YEAR
   
   % Emissivity calculations
   if WAVAFEEDBACK == 1
     EPSA1 = calc_EPSA(Ts1(n),Ta1(n),RHA1,CO2(n),0); %Box1 emissivity (no dim.)
     EPSA2 = calc_EPSA(Ts2(n),Ta2(n),RHA2,CO2(n),0); %Box2 emissivity (no dim.)
   else
     EPSA1 = calc_EPSA(Ts1(1),Ta1(1),RHA1,CO2(n),0); %Box1 emissivity (no dim.)
     EPSA2 = calc_EPSA(Ts2(1),Ta2(1),RHA2,CO2(n),0); %Box2 emissivity (no dim.)
   end
   
   % Atm. & Oce. Circulation strength
   [PSIa(n),PSIo(n)] = calc_Psi(Ts1(n),Ts2(n)); %in kg/s
   
   % Flux calculations
   [Fs1(n),Feva1(n)] = calc_Fs(Te1,Ta1(n),Ts1(n),EPSA1); %net surface heat flux (W m-2)
   Ft1(n) = calc_Ft(Te1,Ta1(n),Ts1(n),EPSA1); %net top-of-atm heat flux (W m-2)
   [Fs2(n),Feva2(n)] = calc_Fs(Te2,Ta2(n),Ts2(n),EPSA2); %net surface heat flux (W m-2)
   Ft2(n) = calc_Ft(Te2,Ta2(n),Ts2(n),EPSA2); %net top-of-atm heat flux (W m-2)
   Fa(n) = calc_Fa(PSIa(n),Ta1(n),Ta2(n)); %atm. heat flux (W m-2)
   Fo(n) = calc_Fo(PSIo(n),To1(n),To2(n)); %oce. heat flux (W m-2)
   
   % Time stepping Box 1 (Tropics)
   Tend_atm = (Fs1(n) + Ft1(n) - Fa(n))/HCA;
   Tend_oce = -(Fs1(n) + Fo(n))/HCO;
   Ta1(n+1) = Ta1(n) + Dt*Tend_atm;
   To1(n+1) = To1(n) + Dt*Tend_oce;
   Ts1(n+1) = To1(n+1) + DTO;
   
   % Time stepping Box 2 (Extra-Tropics)
   Tend_atm = (Fs2(n) + Ft2(n) + Fa(n))/HCA;
   Tend_oce = -(Fs2(n) - Fo(n))/HCO;
   Ta2(n+1) = Ta2(n) + Dt*Tend_atm;
   To2(n+1) = To2(n) + Dt*Tend_oce;   
   Ts2(n+1) = To2(n+1) + DTO;
   
   if PICTUREonline==1
   %On-line plots
   figure(7), clf,
   subplot(1,2,1),
    if n==1,
        [h]=plot(.5*(Ts1(n)+Ts2(n)),.5*(Fs1(n)+Fs2(n)),'r*');
    else
       [h]=plot(.5*(Ts1(n)+Ts2(n)),.5*(Fs1(n)+Fs2(n)),'b*');
    end
    set(h,'markersize',15);
    title(['Time = ' num2str(n*Dt/YEAR) 'yr'],'Fontsize',15);
    xlabel('Global average surface temperature (K)','Fontsize',15);
    ylabel('Global average TOA imbalance (Wm-2)','Fontsize',15);
    set(gca,'Fontsize',15);
    axis([290 293 -8 1]);
   subplot(1,2,2),
    if n==1,
       [h]=plot(Fa(n)*pi*RADIUS^2/PW,Fo(n)*pi*RADIUS^2/PW,'r*');
    else
       [h]=plot(Fa(n)*pi*RADIUS^2/PW,Fo(n)*pi*RADIUS^2/PW,'b*',[0 10],[0 10],'g-',[0 10],[0 5],'b-');
    end
    set(h,'markersize',15);
    title(['Time = ' num2str(n*Dt/YEAR) 'yr'],'Fontsize',15);
    xlabel('Atmospheric Heat Transport (PW)','Fontsize',15);
    ylabel('Oceanic Heat Transport (PW)','Fontsize',15);
    set(gca,'Fontsize',15);
    axis([1 4 1 2.5]);
   end %PICTUREonline
end