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

% read sulphate depositions from SI of Sigl et al. 2015, Nature  
data=xlsread('nature14565-s6.xlsx',2);

% set missing data to zero
data(isnan(data))=0;

% extract fields
location_sigl15=data(:,2); % location flag (1=tropical, 2=NH, 3=SH)
deposition_north=data(:,3); % Greenland
deposition_south=data(:,5); % Antarctica

% convert to stratospheric volcanic aerosol mass (Tg) using 
% scaling relation from Gao et al. 2007, J. Geophys. Res. 
for n=1:length(location_sigl15)
  if location_sigl15(n)==2
    mass_sigl15(n)=0.57*deposition_north(n);
  else
    mass_sigl15(n)=deposition_north(n)+deposition_south(n);
  end
end

% 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 to occur in specific month

% choose and seed random number generator
%   The twister algorithm 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;
  for y=YEAR1:YEARN
    for m=1:12
      for count_sigl15=1:length(mass_sigl15)
        if rand<P
          count_future=count_future+1;
          year{member}(count_future)=y;
          month{member}(count_future)=m;
          location{member}(count_future)=location_sigl15(count_sigl15);
          mass{member}(count_future)=mass_sigl15(count_sigl15);
        end
      end
    end
  end
end

% save result
save futureVolcForcing_test.mat year month location mass 


