clear all

% read GMST from 60 members of NoVolc and Volc 
NoVolc=load('Data/gmst_historicalNoVolc_1990-2004_monthly.txt');
NoVolc=reshape(NoVolc(:,3:end),12,[],60); % dimensions: month,year,member
NoVolc_ensave=squeeze(mean(NoVolc,3));
%
Volc=load('Data/gmst_historicalVolc_1990-2004_monthly.txt');
Volc=reshape(Volc(:,3:end),12,[],60); 
%Volc=Volc(:,3:end);
Volc_ensave=squeeze(mean(Volc,3));
%
tvec=[(1990+0.5/12):(1/12):(2004+11.5/12)]';

% read obs data 
Obs=load('Data/GLB.Ts+dSST_nohead_1880-2013.txt');
% extract SAT data (skip first i.e. time column) and transpose 
Obs=Obs(:,2:13)'/100;
% pick 1990-2004 (dimensions: month,year)
Obs=Obs(:,1+1990-1880:1+2004-1880); 

% read CESM/NorESM historical volcanic forcing 
REARTH=6371000;
FHIS='Data/CCSM4_volcanic_1850-2008_prototype1.nc';
colmass=ncread(FHIS,'colmass');
lat=ncread(FHIS,'lat');
for n=1:length(lat)
if n==1
lats=-90;
latn=0.5*(lat(n)+lat(n+1));
elseif n==length(lat)
lats=0.5*(lat(n)+lat(n-1));
latn=90;
else
lats=0.5*(lat(n)+lat(n-1));
latn=0.5*(lat(n)+lat(n+1));
end
area(n)=2*pi*REARTH^2*(sin(latn/180*pi)-sin(lats/180*pi));
end
% computate global aerosol mass and convert units ot Tg  
forcing=area*colmass*1e-9;
% select time period Jan 1990 to Dec 2004 (data in file starts Dec 1949)
forcing=forcing(:,2+(1990-1850)*12:13+(2004-1850)*12);

% remove seasonal cycle 
scycle=squeeze(mean(NoVolc_ensave,2))-mean(mean(NoVolc_ensave));

for y=1:size(Volc,2)
Volc_ensave(:,y)=Volc_ensave(:,y)-scycle;
NoVolc_ensave(:,y)=NoVolc_ensave(:,y)-scycle;
for mem=1:60
Volc(:,y,mem)=Volc(:,y,mem)-scycle;
NoVolc(:,y,mem)=NoVolc(:,y,mem)-scycle;
end
end 

% subtract (reference to) first year of ensemble mean 
Volc=Volc-mean(Volc_ensave(1:12));
Volc_ensave=Volc_ensave-mean(Volc_ensave(1:12));
NoVolc=NoVolc-mean(NoVolc_ensave(1:12));
NoVolc_ensave=NoVolc_ensave-mean(NoVolc_ensave(1:12));
Obs=Obs-mean(Obs(1:12)); 

% compute quantiles
Volc=reshape(Volc,[],60);
NoVolc=reshape(NoVolc,[],60);
for n=1:180
Volc_95(n)=quantile(bootstrapmean(Volc(n,:),1000),0.95);
Volc_5(n)=quantile(bootstrapmean(Volc(n,:),1000),0.05);
NoVolc_95(n)=quantile(bootstrapmean(NoVolc(n,:),1000),0.95);
NoVolc_5(n)=quantile(bootstrapmean(NoVolc(n,:),1000),0.05);
end 
% find overlap 
VolcNoVolc_95=min(Volc_95,NoVolc_95);
VolcNoVolc_5=max(Volc_5,NoVolc_5);
VolcNoVolc_5(find(VolcNoVolc_5>VolcNoVolc_95))=VolcNoVolc_95(find(VolcNoVolc_5>VolcNoVolc_95));

% prepare linear model
f_linear=inline('-0.5*(t).*exp(-0.6*(t))');
tvec_linear=[1991-0.5/12:1/12:2000-0.5/12]+0.5; 
B=regress(NoVolc_ensave(:)-mean(NoVolc_ensave(:)),tvec-mean(tvec));
i1=find(abs(tvec-tvec_linear(1))<0.01);
i2=find(abs(tvec-tvec_linear(end))<0.01);
Fit=(tvec(:)-mean(tvec(:)))*B+mean(NoVolc_ensave(:));;
Fit(i1:i2)=Fit(i1:i2)+f_linear(tvec_linear-tvec_linear(1))';

% plot
figure(1)
clf
c2=[0.4 0.6 1];c1=[1 0.6 0.4];c12=(c1+c2)/2;
set(gcf,'color',[1 1 1],'InvertHardcopy','off')
subplot(3,1,1)
%%%
% plot GMST response 
hold on
patch([tvec;flipud(tvec)],[NoVolc_95';flipud(NoVolc_5')],c1,'edgecolor','none');
patch([tvec;flipud(tvec)],[Volc_95';flipud(Volc_5')],c2,'edgecolor','none');
patch([tvec;flipud(tvec)],[VolcNoVolc_95';flipud(VolcNoVolc_5')],c12,'edgecolor','none');
plot(tvec,NoVolc_ensave(:),'color',[0.8 0 0],'linewidth',1.5);
plot(tvec,Volc_ensave(:),'color',[0 0 0.8],'linewidth',1.5);
plot(tvec,Obs(:),'k','linewidth',1.5);
 plot(tvec,Fit,'g','linewidth',1.5);
hold off
set(gca,'linewidth',1,'xtick',[1990:1:2005],'xlim',[1990 2005],'ylim',[-0.5 0.5],'ytick',-0.5:0.1:0.5);
set(gca,'fontsize',12,'fontweight','bold');
ylabel('Temperature anomaly (K)','fontsize',12,'fontweight','bold')
text(2005,0.5,'c','fontsize',14,'fontname','arial','fontweight','bold','horiz','right','vert','top')
pos=get(gca,'position');
set(gca,'position',[pos(1) pos(2)-0.60 pos(3) pos(4)*1.4]);
%%%
% plot fitted GMST response in separate panel
axes('position',[pos(1) pos(2)-0.22 pos(3) pos(4)*1.3],'color','none','yaxislocation','left');
hold on
plot([1990 2005],[0 0],':k','color',[1 1 1]*0,'linewidth',1.)
plot(tvec_linear,f_linear(tvec_linear-tvec_linear(1)),'color','g','linewidth',1.5)
hold off
set(gca,'linewidth',1,'xtick',[1990:1:2005],'xlim',[1990 2005],'ylim',[-0.5 0.5],'ytick',-0.5:0.1:0.5);
set(gca,'fontsize',12,'fontweight','bold');
ylabel('Fitted response (K)','fontsize',12,'fontweight','bold')
text(2005,0.2,'b','fontsize',14,'fontname','arial','fontweight','bold','horiz','right','vert','middle')
text(1995,0.5,'-0.5*(t).*exp(-0.6*(t))','fontsize',12,'fontname','arial','fontweight','bold','horiz','left','vert','top')
%%%
% plot volcanic aerosol forcing
axes('position',[pos(1) pos(2)+0.13 pos(3) pos(4)*0.7],'color','none','yaxislocation','left');
hold on
plot(tvec,forcing,'color',[1 1 1]*0,'linewidth',1.5)
ytick=[0:5:25];
dtick=0.2;
xref=1990;
plot([xref xref],[ytick(1) ytick(end)],'-k','linewidth',1.5);
for n=1:length(ytick)
  plot([xref+dtick xref],[ytick(n) ytick(n)],'-k','linewidth',1.5)
  text(xref-0.2,ytick(n),sprintf('%2.0f',ytick(n)),'horiz','right','fontweight','bold','fontsize',12);
end
text(xref-1.04,mean(ytick),'Load (Tg)','horiz','center','vert','bottom','rotation',90,'fontweight','bold','fontsize',12);
text(2005,25,'a','fontsize',14,'fontname','arial','fontweight','bold','horiz','right','vert','middle')

set(gca,'xlim',[1990 2005],'xtick',[1990:1:2005],'ylim',[0 30],'visible','off','ytick',[],'linewidth',1)
set(gca,'fontsize',12,'fontweight','bold');

% save plot
print(gcf,'-depsc2','response_fit_GMST.eps')


