clear all 

VARS={'FSNT','FLNT'}; 
MAXMEM=9;
MONLEN=[31 28 31 30 31 30 31 31 30 31 30 31];
w=MONLEN/sum(MONLEN);


load Data/atmstats_NFPTAEROL32_f19_f19_default_0001-0020.mat 
for VAR1=VARS 
  Off.(char(VAR1))=atmstats.(char(VAR1)); 
end 

load Data/atmstats_NFPTAEROL32_f19_f19_CMIP6volc_2_0001-0020.mat 
for VAR1=VARS 
  Fix.(char(VAR1))=atmstats.(char(VAR1)); 
end 

for VAR1=VARS
  CM5.(char(VAR1))=zeros(12,20);
  CM6.(char(VAR1))=zeros(12,20);
  for mem=1:MAXMEM 
    load(sprintf('Data/atmstats_NFPTAEROL32_f19_f19_PinatuboCMIP5_%i_1990-2009.mat',mem)); 
    CM5.(char(VAR1))=CM5.(char(VAR1))+atmstats.(char(VAR1))/MAXMEM;
    load(sprintf('Data/atmstats_NFPTAEROL32_f19_f19_PinatuboCMIP6_%i_1990-2009.mat',mem)); 
    CM6.(char(VAR1))=CM6.(char(VAR1))+atmstats.(char(VAR1))/MAXMEM;
  end 
end  

% remove seasonal cycle of default
for VAR1=VARS
  clim.(char(VAR1))=mean(Off.(char(VAR1)),2); 
  for m=1:12
    Off.(char(VAR1))(m,:)=Off.(char(VAR1))(m,:)-clim.(char(VAR1))(m);
    Fix.(char(VAR1))(m,:)=Fix.(char(VAR1))(m,:)-clim.(char(VAR1))(m);
    CM5.(char(VAR1))(m,:)=CM5.(char(VAR1))(m,:)-clim.(char(VAR1))(m);
    CM6.(char(VAR1))(m,:)=CM6.(char(VAR1))(m,:)-clim.(char(VAR1))(m);
  end
end 

% plot 
figure(1)
clf 
orient tall
tvec=[(1990+0.5/12):(1/12):(2009+11.5/12)];
subplot(2,1,1)
hold on 
h(1)=plot([1990 2009],[0 0],'k','linewidth',1.5);
h(2)=plot([1990 2009],[mean(w*Fix.FSNT) mean(w*Fix.FSNT)],'k:','linewidth',1.5);
h(3)=plot(tvec,CM5.FSNT(:),'b','linewidth',1.5);
h(4)=plot(tvec,CM6.FSNT(:),'r','linewidth',1.5);
hold off 
set(gca,'xlim',[1991 1995],'xtick',[1991:2000],'ylim',[-5 1],'fontsize',12,'fontweight','bold')
ylabel('TOA net-shortwave (W/m2)','fontsize',12,'fontweight','bold')
legend({'Volc-zero = 0 W/m2','Volc-clim(1850-2000) = -0.35 W/m2','Volc-trans CMIP5','Volc-trans CMIP6'},'location','southeast','fontsize',12,'fontweight','bold');
title('Monthly TOA radiation (relative to NO-VOLC climatology)','fontsize',16,'fontweight','bold')
subplot(2,1,2)
hold on
h(1)=plot([1990 2009],[0 0],'k','linewidth',1.5);
h(2)=plot([1990 2009],[mean(w*(Fix.FSNT-Fix.FLNT)) mean(w*(Fix.FSNT-Fix.FLNT))],'k:','linewidth',1.5);
h(3)=plot(tvec,CM5.FSNT(:)-CM5.FLNT(:),'b','linewidth',1.5);
plot(tvec,CM6.FSNT(:)-CM6.FLNT(:),'r','linewidth',1.5);
hold off
set(gca,'xlim',[1991 1995],'xtick',[1991:2000],'ylim',[-5 1],'fontsize',12,'fontweight','bold')
ylabel('TOA net (W/m2)','fontsize',12,'fontweight','bold')
legend({'Volc-zero = 0 W/m2','Volc-clim(1850-2000) = -0.23 W/m2','Volc-trans CMIP5','Volc-trans CMIP6'},'location','southeast','fontsize',12,'fontweight','bold');
print(gcf,'-dpng','CMIP6_Pinatubo_test.png')  
