
if 1
% ------------  get data  ---------------
% 
% get dR
load('../Data/toaflux_NBF1850OC_f19_tn11_01_0001-1000.mat');
dr0=fluxtoa(end);

load ('../Data/toaflux_NBF18504XOC_f19_tn11_01_0001-0150.mat');
dr=fluxtoa-dr0;

% get dT
load('../../timeseries/Data/timeseries_NBF1850OC_f19_tn11_01_0001-1000.mat');
dt0=ts_yr.t2m_yr(end);


load('../../timeseries/Data/timeseries_NBF18504XOC_f19_tn11_01_0001-0150.mat');
dt=ts_yr.t2m_yr-dt0;

% get effective climate sensitivity
dt_eff=(dt*7.0)./(7.0-dr);



% ------ 1% CO2 ------

% get dR at year 70
load('../Data/toaflux_NBF1850OC_f19_tn11_02_1001-2000.mat');
r70_pi=mean(fluxtoa(60:80));

load ('../Data/toaflux_NBF1850RMOC_f19_tn11_01_0001-0150.mat');
r70=mean(fluxtoa(60:80));

% get dT at year 70
load('../../timeseries/Data/timeseries_NBF1850OC_f19_tn11_02_1001-2000.mat');
t70_pi=mean(ts_yr.t2m_yr(60:80));

load('../../timeseries/Data/timeseries_NBF1850RMOC_f19_tn11_01_0001-0150.mat');
t70=mean(ts_yr.t2m_yr(60:80));

TCR=t70-t70_pi;
dr70=r70-r70_pi;
TCR_eff=TCR*7.0/(7.0-dr70);
%save('climate_sensitivity.mat','dt','dr','TCR');

% --------------------------------------------------
end


% ---------- PLOT ----------------------------------
%

load climate_sensitivity.mat;
%dt(1:20)=[];dr(1:20)=[];

p = polyfit(dt,dr,1);

figure;
set(gcf,'renderer','painters');
plot(dt,dr,'o','color','c');hold on;
plot(dt_eff(111:150),dr(111:150),'o','color','r');hold on;
plot(mean(dt_eff(111:150)),mean(dr(111:150)),'o','MarkerFaceColor','k','markersize',10);hold on;

dtt=0:0.01:7;
plot(dtt,dtt*p(1)+p(2),'linewidth',2,'color','k');

xlim([0 7]);ylim([0 7]);grid on;
title({['ECS= ' sprintf('%5.2f',-p(2)/p(1)/2) ' ^{\circ}C;' ' ECS_{eff}= ' sprintf('%5.2f',mean(dt_eff(111:150))/2) ' ^{\circ}C' ], ...
       ['TCR= ' sprintf('%5.2f',TCR) ' ^{\circ}C; TCR_{eff}=' sprintf('%5.2f',TCR_eff) '^{\circ}C']},'fontsize',15);
set(gca,'fontsize',15);
xlabel('\DeltaT (^{\circ}C)');ylabel('change in TOA net radiation (W m^{-2})');

print -depsc climate_sensitivity.eps
