%
% 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),CO2(n),0); %Box1 emissivity (no dim.)
     EPSA2 = calc_EPSA(Ts2(n),Ta2(n),CO2(n),0); %Box2 emissivity (no dim.)
   else
     EPSA1 = calc_EPSA(Ts1(1),Ta1(1),CO2(n),0); %Box1 emissivity (no dim.)
     EPSA2 = calc_EPSA(Ts2(1),Ta2(1),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

   % "Moist" diagnostics
   MSE1(n) = calc_MSE(Ts1(n),Ta1(n)); %moist static energy
   MSE2(n) = calc_MSE(Ts2(n),Ta2(n));
   MTSPT(n) = calc_MTSPT(Ts1(n),Ts2(n),Ta1(n),Ta2(n)); %moisture transport 
   
   % 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),MSE1(n),MSE2(n)); %atm. heat flux (W m-2)
   Fo(n) = calc_Fo(PSIo(n),Ts1(n),To2(n)); %oce. heat flux (W m-2)
   PSIres = PSIo(n)*CPO/(pi*RADIUS^2); %rescaling of PSIo (W m-2 K-1)
   
   % Time stepping Box 1 (Tropics)
   Tend_atm = (Fs1(n) + Ft1(n) - Fa(n))/HCA;
   Tend_oce1_ml = -( Fs1(n) - PSIres*(To1(n)-Ts1(n)) )/HCM; %mixed layer 
   Tend_oce1_th = PSIres*(To2(n)-To1(n))/HCO; %thermocline 
   Ta1(n+1) = Ta1(n) + Dt*Tend_atm;
   Ts1(n+1) = Ts1(n) + Dt*Tend_oce1_ml;
   To1(n+1) = To1(n) + Dt*Tend_oce1_th;
   
   % Time stepping Box 2 (Extra-Tropics)
   Tend_atm = (Fs2(n) + Ft2(n) + Fa(n))/HCA;
   Tend_oce2_ml = -( Fs2(n) - PSIres*(Ts1(n)-Ts2(n)) )/HCM; %mixed layer 
   Tend_oce2_th = PSIres*(Ts2(n)-To2(n))/HCO; %thermocline 
   Ta2(n+1) = Ta2(n) + Dt*Tend_atm;
   Ts2(n+1) = Ts2(n) + Dt*Tend_oce2_ml;
   To2(n+1) = To2(n) + Dt*Tend_oce2_th;   
   
   %Test energy conservation for Ocean
   %Sum1 = HCM*Tend_oce1_ml + Fs1(n) + HCO*Tend_oce1_th, %in Wm-2
   %Sum2 = HCM*Tend_oce2_ml + Fs2(n) + HCO*Tend_oce2_th, %in Wm-2
   %Sum3 = Fo(n), %in Wm-2
   %return
   %NB: worked fine.
   
   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
