% created 12 Oct 2016 (ingo.bethke@uni.no)

% clean start 
clear all 

% Read peak stratospheric volcanic aerosol mass loads
%   These are derived from sulphate depositions of Sigl et al. 2015 and
%   converted to stratospheric aerosol mass using the Gao et al. 2007 scaling
data=load('nature14565-s6_stratAerosolMass.txt');
loc_sigl=data(:,2); % location flag (1=tropical, 2=NH, 3=SH)
mass_sigl=data(:,3); % stratospheric volcanic aerosol mass (Tg)

% define ensemble size and time period for synthetic forcing
NMEMBER=60; % number of forcing members
YEAR1=2006;
YEARN=2099;
P=1./2500./12.; % probability for eruption event to occur in specific month

% choose and seed random number generator 
%   twister is default in python and new matlab versions
%   the below setting gives the same result as the numpy python implementation 
rand('twister',0);  
 
% create forcing futures 
for member=1:NMEMBER
    count_future=0; 
    disp(['Member ' int2str(member)]) 
    for year=YEAR1:YEARN 
        for month=1:12
            for count_sigl=1:length(mass_sigl)
                if rand<P  
                   count_future=count_future+1;
                   futureVolcForcing{member}.loc(count_future)=loc_sigl(count_sigl);
                   futureVolcForcing{member}.mass(count_future)=mass_sigl(count_sigl);
                   futureVolcForcing{member}.month(count_future)=month;
                   futureVolcForcing{member}.year(count_future)=year;
                end
            end
        end
    end
end 

% save result 
save futureVolcForcing_new.mat futureVolcForcing loc_sigl mass_sigl 

% plot 
member=35
figure(1)
clf
bar(futureVolcForcing{member}.year,futureVolcForcing{member}.mass)
set(gca,'xlim',[2006 2100],'xtick',[2010:10:2100],'ylim',[0 5]);
ylabel('Aerosol mass (Tg)')
title(['Member ' int2str(member)])
print(gcf,'-dpng','-r300',['futureVolcForcing_mem' int2str(member)])
