% This program loads and analyses global sea level and temperature data as
% described in Rahmstorf (Science 2007). It also produces the graphs shown 
% in that paper. For further description and references, please read the
% paper.
%
% Code written September 2006, upgraded July 2007, 
% by Stefan Rahmstorf, Potsdam Institute for Climate Impact Research 
% Please report any problems to the author at stefan@pik-potsdam.de
%
% Reference:
% Rahmstorf, S., 2007: A semi-empirical approach to projecting future
% sea-level rise. Science, 315, 368-370.

clear; %clear work space - take care, this removes all data currently in matlab's memory

%----------------------------
% Parameter and options setting
%----------------------------
% Note that my paper uses the SSAtrend of Moore et al. 2005 for smoothing the data. 
% As a simple matlab alternative, fitting a polynomial is available as
% an option below. This is not recommended; one can see in figure 2 (top)
% that this gives only a poor fit, although it has little effect on the future
% projections. I highly encourage the use of SSAtrend
% or similar types of filtering to remove short-period noise. Using
% SSAtrend with embedding dimension 15 exactly reproduces my published
% results and graphs. 
% Reference: J. C. Moore, A. Grinsted, S. Jevrejeva, Eos 86, 226 (2005).

detr='n'; % set to y or n, depending on whether detrending is wanted
SSA='y'; % set to n if you don't have access to the SSAtrend smoother
ed = 15 % embedding dimension (years) of SSAtrend smoother (15 was used in the paper)
degree = 8; % degree of polynomial fit - only used if no SSAtrend available - not recommended!
bin='5'; % choose 5 (as in the paper) or 15 year bins

%----------------------------
% loading data
%----------------------------

% load the church & white sea level data
load church_13221.txt;
seayear = church_13221(:,1);
sealevel = church_13221(:,2)/10;

% load GISS global temperature data 
load 'giss_landocean.txt';
gissyear=giss_landocean(:,1); % column 1 is the year
gisstemp=(giss_landocean(:,2) .* 0.01); % Column 2 is the annual global value
 
% load the IPCC 2001 global temperature scenarios
load ipcc_scenarios.dat
ipccyear=ipcc_scenarios(:,1);
% prepare vectors for shading low and high limits
xt=[ipccyear; flipud(ipccyear)];
yt=[ipcc_scenarios(:,2); flipud(ipcc_scenarios(:,3))]; 

%----------------------------
% apply smoothing to time series
%----------------------------

if SSA=='y'
  smoothsealevel = ssatrend([seayear sealevel],ed);
  smoothgisstemp = ssatrend([gissyear gisstemp],ed);
% alternative smoothing using polynomial fit
else
  [p1 s1 m1] = polyfit(seayear,sealevel,degree);
  [smoothsealevel error1] = polyval(p1,seayear,s1,m1);
  [p2 s2 m2] = polyfit(gissyear,gisstemp,degree);
  [smoothgisstemp error2] = polyval(p2,gissyear,s2,m2);
end
g1990=(smoothsealevel(120)+smoothsealevel(121))/2; %smoothed 1990 value, subtracted below
sealevel=sealevel-g1990;
smoothsealevel=smoothsealevel-g1990;

%----------------------------
% compute and look at correlations
%----------------------------

% compute average rate of rise and average temperature over time intervals (i.e., bin data) 
rateofrise=diff(smoothsealevel); % time derivative of smoothed sea level
if bin=='5'
  for i=1:24
    mtemp(i)=mean(smoothgisstemp(5*i-2:5*i+2));
    mrate(i)=mean(rateofrise(5*i+7:5*i+11));
    myear(i)=mean(gissyear(5*i-2:5*i+2));
  end
end
if bin=='15'
  for i=1:8
    mtemp(i)=mean(smoothgisstemp(15*i-12:15*i+2));
    mrate(i)=mean(rateofrise(15*i-3:15*i+11));
    myear(i)=mean(gissyear(15*i-12:15*i+2));
  end
end

% correlation of these binned data:

% optional detrending before fitting 
% (we need to keep the mean which the matlab detrend function removes)
if detr=='y'
    mtemp = detrend(mtemp) + mean(mtemp);
    mrate = detrend(mrate) + mean (mrate);
end

[r p] = corrcoef([mtemp' mrate']);
disp('correlation coefficient: '); disp(r(1,2));  % correlation coefficient
disp('p-value: '); disp(p(1,2));  % probability of this happening by chance with uncorrelated data

% linear fit to these data:
[po s] = polyfit(mtemp,mrate,1); % s is used with "polyval" to compute error bounds
fitrate = polyval(po,[-.48 -.25 0 .25 .5],s);
slope = po(1)
tbase = -po(2)/po(1) % this is the derived pre-industrial base temperature

% compute sea level hindcast from the fit found above 
[hindrate herror] = polyval(po,smoothgisstemp,s); %hindcast rate of sea level rise
hindlevel = cumsum(hindrate); %integrate up the rate
hindlevel = hindlevel - hindlevel(111); % normalise to 1990 value

% compute sea level projections from the fit found above
for i=[2 3 6:11]
   [predrate perror] = polyval(po,ipcc_scenarios(:,i)+smoothgisstemp(111),s); %forecast rate of sea level rise
   pred = cumsum(5*predrate); %integrate up the rate; time interval here is 5 years
   predmax = cumsum(5*(predrate+perror));
   predmin = cumsum(5*(predrate-perror));
   predlevel(:,i) = pred - pred(1); % normalise for 1990
   predlevelmax(:,i) = predmax-predmax(1);
   predlevelmin(:,i) = predmin-predmin(1);
end

%-----------------------
% Prepare graphs
%-----------------------
ratio = [1.5 1 1]; % aspect ratio of the plots

% scatter plot of actual vs modelled rate of sea level rise
figure(1); clf; hold on; box on;
p1=plot(mtemp,mrate*10,'or');
set(p1,'markersize',5,'markerfacecolor','r');
plot([-.48 -.25 0 .25 .5],fitrate*10,'--');
ylabel('Rate of Sea Level Change (mm/year)')
xlabel('Warming above 1951-1980 mean')
set(gca,'FontSize',8,'plotboxaspectratio',ratio)
set(gcf,'paperunits','centimeters','paperpos',[4 15.0 12.5 7])

% plot of time series of actual and modelled rate of rise
figure(2); clf; hold on; box on;
%plot error band of modelled rate
xh=[gissyear; flipud(gissyear)];
yh=[(hindrate+herror)*10; flipud((hindrate-herror)*10)];
f2=fill(xh,yh,[.8 .8 1],'marker','.','markersize',0.01,'edgecolor',[.8 .8 1]);

% in case of SSAtrend smoothing, the observed rate (rateofrise) is smoothed, 
% since taking the derivative makes it look noisy (unlike the polynomial fit)
if SSA=='y'
  smoothrate = ssatrend([seayear(1:131) rateofrise],ed);
else
  smoothrate = rateofrise;
end

plot(gissyear,hindrate*10,seayear(1:131),smoothrate*10,'r');

axis([1865 2010 0 4.3])
ylabel('Rate of Change (mm/yr)')
xlabel('Year')
set(gca,'FontSize',8,'plotboxaspectratio',ratio,'Position',[0.2 0.65 .6 0.3])
set(gcf,'paperunits','centimeters','paperpos',[4 5.0 15 18])

% plot time series of actual and modelled sea level
width=.3; %linewidth
axb=axes('Position',[.2 0.35 .6 0.3]); hold on; box on;

c5 = line(seayear,sealevel,'color',[.9 .6 .6],'marker','square','markerfacecolor',[.9 .6 .6],'markersize',2,'linewidth',.1);
c7 = line(gissyear,hindlevel,'Color','b','LineWidth',width);
c6 = line(seayear,smoothsealevel,'Color','r','LineWidth',width); 

axis([1865 2010 -17 7])
xlabel('Year')
ylabel('Sea Level (cm)')
set(gca,'FontSize',8,'plotboxaspectratio',ratio);

% plot with future projections
figure(3); clf; hold on; box on;

%plot ipcc range as grey band
xs=[ipccyear; flipud(ipccyear)];
ys=[predlevel(:,2); flipud(predlevel(:,3))];
f1=fill(xs,ys,[.8 .8 .8],'marker','.','markersize',0.01,'edgecolor',[.8 .8 .8]); %IPCC band
plot([1000 2100],[0,0],'linewidth',.5,'color',[.5 .5 .5]);
plot([1990 1990],[-100,300],'linewidth',.5,'color',[.5 .5 .5]);

% plot projected sea level for ipcc marker scenarios
a1fi=plot(ipccyear,predlevel(:,6),'--');
a1b=plot(ipccyear,predlevel(:,7),'--');
a2=plot(ipccyear,predlevel(:,9),'--');
b1=plot(ipccyear,predlevel(:,10),'--');
b2=plot(ipccyear,predlevel(:,11),'--');
a1t=plot(ipccyear,predlevel(:,8),'--');
max=plot(ipccyear,predlevelmax(:,2),'--');
min=plot(ipccyear,predlevelmin(:,3),'--');

% add observed and hindcast sea level
c5 = line(seayear,sealevel,'color',[.9 .6 .6],'marker','square','markerfacecolor',[.9 .6 .6],'markersize',2,'linewidth',.1);
c6 = line(seayear,smoothsealevel,'Color','r','LineWidth',width); 
c7 = line(gissyear,hindlevel,'Color','b','LineWidth',width);

% set colors etc
width=0.4;
set(max,'linewidth',width,'color',[0.5 0.5 0.5])
set(min,'linewidth',width,'color',[0.5 0.5 0.5])
set(a1fi,'linewidth',width,'color',[0 1 1])
set(a1b,'linewidth',width,'color',[1 0 1])
set(a1t,'linewidth',width,'color','b')
set(a2,'linewidth',width,'color','r')
set(b1,'linewidth',width,'color',[1 1 0])
set(b2,'linewidth',width,'color','g')

axis([1870 2100 -20 140])
set(gca,'FontSize',8,'plotboxaspectratio',ratio)
xlabel('Year')
ylabel('Sea Level Change (cm)')
set(gcf,'paperunits','centimeters','paperpos',[4 15.0 12.5 7])
