%
% Diagnostics on Two-Box Model
%
% Global average
  Ts = (Ts1+Ts2)/2;
  Ft = (Ft1+Ft2)/2;
  
% Heat transport
  Ha_PW = (Fa*pi*RADIUS^2)/PW;
  Ho_PW = (Fo*pi*RADIUS^2)/PW;

% Hydrological cycle
  EVAP1 = Feva1 * (pi*RADIUS^2)/LV; % evaporation in kg/s
  PRCP1 = EVAP1 - MTSPT; %precip in kg/s
  PRCP1_mmday = 1000*PRCP1*DAY/(RHO*pi*RADIUS^2);
  EVAP2 = Feva2 * (pi*RADIUS^2)/LV; %evaporation in kg/s
  PRCP2 = EVAP2 + MTSPT; %precip in kg/s
  PRCP2_mmday = 1000*PRCP2*DAY/(RHO*pi*RADIUS^2);
  Q1 = ( MSE1 - CPA*(Ta1(1:end-1)+Ts1(1:end-1))/2 )/LV; %low level specific humidity (kg/kg)
  Q2 = ( MSE2 - CPA*(Ta2(1:end-1)+Ts2(1:end-1))/2 )/LV; %low level specific humidity
  
% Heat content
  HeatContent1_ml = HCM * Ts1;
  HeatContent2_ml = HCM * Ts2;
  HeatContent1_th = HCO * To1;
  HeatContent2_th = HCO * To2;
  HeatContent1_tot = HCM * Ts1 + HCO * To1;
  HeatContent2_tot = HCM * Ts2 + HCO * To2;
  HeatContent_tot = HeatContent1_tot + HeatContent2_tot;
  
% Final time solution
Ta1_end = Ta1(end);
Ta2_end = Ta2(end);
To1_end = To1(end),
To2_end = To2(end),
Ft1_end = Ft1(end);
Ft2_end = Ft2(end);
Ft_end = (Ft1_end+Ft2_end)/2,
Ha_end = Fa(end)*pi*(RADIUS^2)/PW;
Ho_end = Fo(end)*pi*(RADIUS^2)/PW;
Ts1_end = Ts1(end),
Ts2_end = Ts2(end),
Ts_end = (Ts1_end+Ts2_end)/2,
EPSA1_end = EPSA1,
EPSA2_end = EPSA2,

% Fraction of days experiencing convection
Frac_ConvDays1 = 1 - length(find(Feva1==0))/length(Feva1),
Frac_ConvDays2 = 1 - length(find(Feva2==0))/length(Feva2),

