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))=w*atmstats.(char(VAR1)); 
end 

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

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

% plot 
figure(1)
clf 
REFS=mean(Off.FSNT);
REFN=mean(Off.FSNT-Off.FSNT-Off.FLNT);
hold on 
plot([1990 2009],[0 0],'b');
plot([1990 2009],[mean(Fix.FSNT) mean(Fix.FSNT)]-REFS,'b:');
plot(1990:2009,CM5.FSNT-REFS,'b');
plot(1990:2009,CM6.FSNT-REFS,'r');
hold off 
set(gca,'xlim',[1990 2000])

