Skip to content

Commit

Permalink
Initial commit
Browse files Browse the repository at this point in the history
  • Loading branch information
Bob Kopp committed Dec 6, 2016
0 parents commit 2480e56
Show file tree
Hide file tree
Showing 30 changed files with 2,465 additions and 0 deletions.
224 changes: 224 additions & 0 deletions Calc_SESLConterfact.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,224 @@
function S = Calc_SESLConterfact(S,CFtemp,column)

% Calculate sea level with different temperature inputs, factual and counterfactual.
%
% S = Calc_SESLConterfact(S,CFtemp,column)
%
% INPUT:
% - S -> Output Structure of calc_SESL_Prc.m
% - CFtemp -> string, input temperature to use:
% - 'HadCRUT' - HadCRUT4 temperature,
% - 'CMIP5_mean' - 'historicalNat_CMIP35_tas' mean from Gareth Jones <[email protected]>
% - 'CMIP5' - 'historicalNat_CMIP35_tas' single simu (data column given as input) from Gareth Jones <[email protected]>
% - 'mean' - mean temperature between 500-1800 CE
% - 'linrate' - linear rate temperature between 500-1800 CE
% - column -> Data column of 'historicalNat_CMIP35_tas' if CFtemp=='CMIP5'
% - Below set the first ['fyr'] and last ['lyr'] year of simulation and the
% whether to output also the input structure or only the calculation below
%
% OUTPUT: For each in the input selected CFtemp, the calulated sea level
% (sl), temperature (Tcf), equilibrium temperature (T01) and year
% values (time) are given alongside with the general settings. If
% below output_all == true some of the output of previous scripts
% such as data is repeated.

fyr = 1900; % first forecast year
lyr = 2015; % last forecastyear

output_all = false; %Output all (true) or only the forecast structure

sl = []; Tcf = []; T01 = [];
if strcmp(CFtemp,'HadCRUT')
Tfc = load(fullfile('Data','T_HC4_gl_2015.mat'));
elseif strcmp(CFtemp, 'CMIP5_mean')
Tfc = load(fullfile('Data','historicalNat_CMIP35_tas.mat'));
elseif strcmp(CFtemp, 'CMIP5_rand')
Tfc = load(fullfile('Data','historicalNat_CMIP35_tas.mat'));
Tfc.Tfc = Tfc.T(:,column);
end
fprintf('\t "%s" temperature scenario:\n',CFtemp)
pb = 1000; % show a progress bar, progressing every 1000 sample*Tnum
if S.settings.sample*S.settings.Tnum<pb; pb=S.settings.sample*S.settings.Tnum;end;
fprintf('%1.0f',zeros(1,S.settings.sample*S.settings.Tnum/pb));
fprintf('\n')
for i_P = 1:S.settings.sample
for i_T = 1:S.settings.Tnum
i = (i_P-1)*S.settings.Tnum + i_T;
if mod(i,pb) == 0
fprintf('%1.0f',0)
end
if strcmp(CFtemp,'mean')
Tfc_ = S.MH.Tm500_1800(i);
yrfc_ = fyr:lyr;
Tfc_ = Tfc_*ones(size(yrfc_));
elseif strcmp(CFtemp,'linrate')
po = S.MH.Tpo500_1800(i,:);
yrfc_ = fyr:lyr;
Tfc_ = polyval(po,yrfc_);
Tfc_ = Tfc_ - Tfc_(1) + S.MH.Tm1801_1900(i);
elseif strcmp(CFtemp,'HadCRUT')
Tfc_ = Tfc.T(:,2)';
yrfc_ = Tfc.T(:,1)';
Tfc_ = Tfc_ - mean(Tfc_(yrfc_>=1850 & yrfc_<=1900)) + S.MH.Tm1850_1900(i);
Tfc_ = Tfc_(yrfc_>=fyr);
yrfc_ = yrfc_(yrfc_>=fyr);
elseif strcmp(CFtemp, 'CMIP5_mean')
Tfc_ = Tfc.T(:,3)';
yrfc_ = Tfc.T(:,1);
Tfc_ = Tfc_ - mean(Tfc_(yrfc_>=1860 & yrfc_<=1900)) + S.MH.Tm1860_1900(i);
Tfc_ = Tfc_(yrfc_>=fyr);
yrfc_ = yrfc_(yrfc_>=fyr);
elseif strcmp(CFtemp, 'CMIP5')
Tfc_ = Tfc.Tfc';
yrfc_ = Tfc.T(:,1);
Tfc_ = Tfc_ - mean(Tfc_(yrfc_>=1860 & yrfc_<=1900)) + S.MH.Tm1860_1900(i);
Tfc_ = Tfc_(yrfc_>=fyr);
yrfc_ = yrfc_(yrfc_>=fyr);
end

FC = calc_fc(S,Tfc_,i_P,i);
sl = [sl;FC.sl];
Tcf = [Tcf;Tfc_];
T01 = [T01;FC.T01];

end
end

fprintf('\n')

S.proj.sl = sl;
S.proj.Tcf = Tcf;
S.proj.T01 = T01;
S.proj.time = yrfc_;

if lyr<=2000
% Calculate likelihoods & Differences
L = calc_lik(S);
S.proj.Lik = L.Lik;
S.proj.Lik_yrly = L.Lik_yrly;
S.proj.Diff = L.Diff;
S.proj.Diff_yrly = L.Diff_yrly;
S.proj.Perc = L.Perc;
S.proj.Perc_yrly = L.Perc_yrly;
S.proj.slFactual = L.sl;
end
if strcmp(CFtemp, 'CMIP5')
S.proj.column_CMIP5 = column;
end
if ~output_all
S_ = S.proj;
S_.settings = S.settings;
S = S_;
end

end

function fc = calc_fc(S,Tfc,i_P,i)

set = S.settings;

T01 = S.MH.T01_1900(i);
if strcmp(S.settings.model,'TwoTau')
T02 = S.MH.T02_1900(i);
end

if strcmp(S.settings.model,'CRdecay')
c = S.MH.c_1900(i);
else
c = S.MH.Params(i_P,3);
end

a1 = S.MH.Params(i_P,1);
a2 = S.MH.Params(i_P,2);
tau1 = S.MH.Params(i_P,4);
tau2 = S.MH.Params(i_P,5);

nyrs = max(size(Tfc));

for ii = 1:nyrs-1
% Compute here T0(sample,year) as the tau-year memory of temp(j)
T01(ii+1) = T01(ii) + (1/tau1)*(Tfc(ii) - T01(ii));
if strcmp(set.model,'TwoTau')
T02(ii+1) = T02(ii) + (1/tau2)*(Tfc(ii) - T02(ii));
elseif strcmp(set.model,'CRdecay')
c(ii+1) = c(ii) * (1-1/tau2);
end
end
if strcmp(set.model,'TwoTau')
dsl = a1*(Tfc - T01) + a2*(Tfc - T02);
fc.T02 = T02;
elseif strcmp(set.model,'ConstRate')
dsl = c + a1*(Tfc - T01);
elseif strcmp(set.model,'CRdecay')
dsl = c + a1*(Tfc - T01);
elseif strcmp(set.model,'CRovTau')
dsl = c/tau1 + a1*(Tfc - T01);
elseif strcmp(set.model,'simpel')
dsl = a1*(Tfc - T01);
end
sl = cumsum(dsl);
sl = sl-sl(1);

fc.sl = sl;
fc.T01 = T01;
end

function L = calc_lik(S)

% which SL to take as factual
fact_SL = 'calib'; % 'calib': semi-emp calibration
% 'prox': proxy data
% 'recalc': semi-empiricaly recalculate with factual temperature
%---- SL projection under counterfactual T
sl = S.proj.sl;
yr = S.proj.time;

%---- SL to calculate likelihood against: either Proxy or semi emp. calibration
if strcmp(fact_SL,'calib')
yrPrx = 1900:2000;
slPrx = S.MH.sl_2000;
elseif strcmp(fact_SL,'recalc')
Tfc = S.data.temp.obs';
yrPrx = S.data.temp.yr;
ix = yrPrx>=yr(1) & yrPrx<=yr(end);
fc = calc_fc(S,Tfc(ix),yrPrx(ix),[50,16,84]);
slPrx = fc.sl;
yrPrx = yrPrx(ix);
end
% cut proxy to same period as T scenario
ix = yrPrx>=yr(1) & yrPrx<=yr(end);
yrPrx = yrPrx(ix);
slPrx = slPrx(:,ix);
for i = 1:size(sl,1)%count sample
Diff_(i,:) = -sl(i,:)+slPrx(i,:);
Perc_(i,:) = (sl(i,:)./slPrx(i,:))*100;
for j = 1:size(sl,2)%count years
L_(i,j) = sum(sl(:,j)>=slPrx(i,j))/size(sl,1);
end
end

% index to give data the spacing of SL & T
ix = ismember(yrPrx,S.MH.yr_);
ix(yrPrx==yr(1))=true;
ix(yrPrx==yr(end))=true;

Lik(:,1) = yrPrx; % the year at which the likelihoods are calculated
Lik(:,2) = mean(L_,1);
L.Lik_yrly = Lik;
L.Lik = Lik(ix,:);

D = prctile(Diff_,[50 17 83 5 95]);
L.Diff_yrly(:,1) = yrPrx;
L.Diff_yrly(:,2:6) = D';
L.Diff = L.Diff_yrly(ix,:);

P = prctile(Perc_,[50 17 83 5 95]);
P(:,1) = 100;
L.Perc_yrly(:,1) = yrPrx;
L.Perc_yrly(:,2:6) = P';
L.Perc = L.Perc_yrly(ix,:);

L.sl(:,1) = yrPrx(ix);
L.sl(:,2:6) = prctile(slPrx(:,ix),[50 17 83 5 95])';

end
152 changes: 152 additions & 0 deletions Calc_SESLProjection.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
function O = Calc_SESLProjection(S,rcps,dig, NumTmag)

% Calculate sea-level projections under different RCP scenarios.
%
% O = Calc_SESLProjection(S,rcps,dig, NumTmag)
%
% INPUT:
% - S -> Output of Calc_SESL_Prc.m
% - rcps -> Cell array filled with strings. RCPs op use: {'RCP3PD', 'RCP45', 'RCP85'}
% - dig -> number of digits to use, else []
% - NumTmag -> Number of MAGICC RCP realizations to take max&default == 600
%
% OUTPUT: Sea-level time-series for each selected RCP (Psl) and year dates (Pslyr).
% The columns in Psl give the percentiles Prc as defined below

usesingle = true; % to save memory

timeprojected = [2000 2100];

Prc = [5 17 50 83 95]; % Percentiles (68.27=1sigma)

folder = 'Data';

if usesingle
SFnc = @(x)single(x);
else
SFnc = @(x)x;
end

if isempty(NumTmag)
NumTmag=600;
end
Tnum = S.settings.Tnum;
sample = S.settings.sample;

if ~isempty(dig)
digits(dig)
end
a = SFnc(S.MH.Params(:,1));
a = reshape(repmat(a',Tnum,1),Tnum*sample,1);
a = repmat(a,NumTmag,1);
if strcmp(S.settings.model,'CRdecay')
c_2000 = SFnc(S.MH.c_2000);
c_2000 = repmat(c_2000,NumTmag,1);

tau_c = SFnc(S.MH.Params(:,5));
tau_c = reshape(repmat(tau_c',Tnum,1),Tnum*sample,1);
tau_c = repmat(tau_c,NumTmag,1);
elseif strcmp(S.settings.model,'simpel')
c = SFnc(zeros(size(S.MH.Params(:,3))));
c = reshape(repmat(c',Tnum,1),Tnum*sample,1);
c = repmat(c,NumTmag,1);
elseif strcmp(S.settings.model,'CRovTau')
c = SFnc(S.MH.Params(:,3)./S.MH.Params(:,4));
c = reshape(repmat(c',Tnum,1),Tnum*sample,1);
c = repmat(c,NumTmag,1);
elseif strcmp(S.settings.model,'ConstRate')
c = SFnc(S.MH.Params(:,3));
c = reshape(repmat(c',Tnum,1),Tnum*sample,1);
c = repmat(c,NumTmag,1);
end
tau = SFnc(S.MH.Params(:,4));
tau = reshape(repmat(tau',Tnum,1),Tnum*sample,1);
tau = repmat(tau,NumTmag,1);

T0_2000 = SFnc(S.MH.T01_2000);
T0_2000 = repmat(T0_2000,NumTmag,1);

if strcmp(S.settings.model,'TwoTau')
a2 = SFnc(S.MH.Params(:,2));
a2 = reshape(repmat(a2',Tnum,1),Tnum*sample,1);
a2 = repmat(a2,NumTmag,1);
T02_2000 = SFnc(S.MH.T02_2000);
T02_2000 = repmat(T02_2000,NumTmag,1);
end

Toff1 = SFnc(S.MH.Tm1970_2000);
Toff1 = repmat(Toff1,NumTmag,1);
settings = S.settings;
clear S

for j = 1:length(rcps) % Count RCPs

sl = SFnc(zeros(size(a)));

if strcmp(settings.model,'CRdecay')
c = c_2000;
end
T0 = T0_2000;
if strcmp(settings.model,'TwoTau')
T02 = T0_2000;
end
% load MAGICC realizations of RCPS: TIME & DATA
if usesingle
load(fullfile(folder, [cell2mat(rcps(j)), '_single'])); % load rcp
else
load(fullfile(folder, cell2mat(rcps(j)))); % load rcp
end

DATA = DATA(:,1:NumTmag);

Toff2 = repmat(mean(DATA(TIME>=1970 & TIME<=2000,:)),size(DATA,1),1);

DATA = DATA-Toff2;

DATA = DATA(TIME>=timeprojected(1) & TIME<=timeprojected(end),:);

clear TIME T_90prc colheaders textdata Toff2

fprintf('\t %1.0f years of %1s are now calculated: \n',size(DATA,1),cell2mat(rcps(j)));

for i_yr = 1:size(DATA,1); % count years to be projected
fprintf('%1.0f \t',i_yr);

Tfc = DATA(i_yr,:);
Tfc = reshape(repmat(Tfc,Tnum*sample,1),NumTmag*Tnum*sample,1);
Tfc = Tfc + Toff1;

if i_yr>1 % so it will be relative to 2100
if ~strcmp(settings.model,'TwoTau')
sl = sl + a .* (Tfc-T0) + c;
else
sl = sl + a .* (Tfc-T0) + a2 .* (Tfc-T02);
end
end

T0 = T0 + (1./tau).*(Tfc-T0);
if strcmp(settings.model,'TwoTau')
T02 = T02 + (1./tau).*(Tfc-T0);
end
if strcmp(settings.model,'CRdecay')
c = c .* (ones(size(c))-1./tau_c);
end
clear Tfc

O.(cell2mat(rcps(j))).Psl(i_yr,:) = prctile(sl,Prc);
O.(cell2mat(rcps(j))).Pslyr=timeprojected(1):timeprojected(end);
end
fprintf('\n')
end
end



% sample Tnum 600
% P1 T1 TR1
% P1 T2 TR1
% P1 T3 TR1
%
% P2 T1 TR1
% P2 T2 TR1
% P2 T3 TR1
Loading

0 comments on commit 2480e56

Please sign in to comment.