Skip to content

Commit 2480e56

Browse files
author
Bob Kopp
committed
Initial commit
0 parents  commit 2480e56

30 files changed

+2465
-0
lines changed

‎Calc_SESLConterfact.m‎

Lines changed: 224 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,224 @@
1+
function S = Calc_SESLConterfact(S,CFtemp,column)
2+
3+
% Calculate sea level with different temperature inputs, factual and counterfactual.
4+
%
5+
% S = Calc_SESLConterfact(S,CFtemp,column)
6+
%
7+
% INPUT:
8+
% - S -> Output Structure of calc_SESL_Prc.m
9+
% - CFtemp -> string, input temperature to use:
10+
% - 'HadCRUT' - HadCRUT4 temperature,
11+
% - 'CMIP5_mean' - 'historicalNat_CMIP35_tas' mean from Gareth Jones <gareth.s.jones@metoffice.gov.uk>
12+
% - 'CMIP5' - 'historicalNat_CMIP35_tas' single simu (data column given as input) from Gareth Jones <gareth.s.jones@metoffice.gov.uk>
13+
% - 'mean' - mean temperature between 500-1800 CE
14+
% - 'linrate' - linear rate temperature between 500-1800 CE
15+
% - column -> Data column of 'historicalNat_CMIP35_tas' if CFtemp=='CMIP5'
16+
% - Below set the first ['fyr'] and last ['lyr'] year of simulation and the
17+
% whether to output also the input structure or only the calculation below
18+
%
19+
% OUTPUT: For each in the input selected CFtemp, the calulated sea level
20+
% (sl), temperature (Tcf), equilibrium temperature (T01) and year
21+
% values (time) are given alongside with the general settings. If
22+
% below output_all == true some of the output of previous scripts
23+
% such as data is repeated.
24+
25+
fyr = 1900; % first forecast year
26+
lyr = 2015; % last forecastyear
27+
28+
output_all = false; %Output all (true) or only the forecast structure
29+
30+
sl = []; Tcf = []; T01 = [];
31+
if strcmp(CFtemp,'HadCRUT')
32+
Tfc = load(fullfile('Data','T_HC4_gl_2015.mat'));
33+
elseif strcmp(CFtemp, 'CMIP5_mean')
34+
Tfc = load(fullfile('Data','historicalNat_CMIP35_tas.mat'));
35+
elseif strcmp(CFtemp, 'CMIP5_rand')
36+
Tfc = load(fullfile('Data','historicalNat_CMIP35_tas.mat'));
37+
Tfc.Tfc = Tfc.T(:,column);
38+
end
39+
fprintf('\t "%s" temperature scenario:\n',CFtemp)
40+
pb = 1000; % show a progress bar, progressing every 1000 sample*Tnum
41+
if S.settings.sample*S.settings.Tnum<pb; pb=S.settings.sample*S.settings.Tnum;end;
42+
fprintf('%1.0f',zeros(1,S.settings.sample*S.settings.Tnum/pb));
43+
fprintf('\n')
44+
for i_P = 1:S.settings.sample
45+
for i_T = 1:S.settings.Tnum
46+
i = (i_P-1)*S.settings.Tnum + i_T;
47+
if mod(i,pb) == 0
48+
fprintf('%1.0f',0)
49+
end
50+
if strcmp(CFtemp,'mean')
51+
Tfc_ = S.MH.Tm500_1800(i);
52+
yrfc_ = fyr:lyr;
53+
Tfc_ = Tfc_*ones(size(yrfc_));
54+
elseif strcmp(CFtemp,'linrate')
55+
po = S.MH.Tpo500_1800(i,:);
56+
yrfc_ = fyr:lyr;
57+
Tfc_ = polyval(po,yrfc_);
58+
Tfc_ = Tfc_ - Tfc_(1) + S.MH.Tm1801_1900(i);
59+
elseif strcmp(CFtemp,'HadCRUT')
60+
Tfc_ = Tfc.T(:,2)';
61+
yrfc_ = Tfc.T(:,1)';
62+
Tfc_ = Tfc_ - mean(Tfc_(yrfc_>=1850 & yrfc_<=1900)) + S.MH.Tm1850_1900(i);
63+
Tfc_ = Tfc_(yrfc_>=fyr);
64+
yrfc_ = yrfc_(yrfc_>=fyr);
65+
elseif strcmp(CFtemp, 'CMIP5_mean')
66+
Tfc_ = Tfc.T(:,3)';
67+
yrfc_ = Tfc.T(:,1);
68+
Tfc_ = Tfc_ - mean(Tfc_(yrfc_>=1860 & yrfc_<=1900)) + S.MH.Tm1860_1900(i);
69+
Tfc_ = Tfc_(yrfc_>=fyr);
70+
yrfc_ = yrfc_(yrfc_>=fyr);
71+
elseif strcmp(CFtemp, 'CMIP5')
72+
Tfc_ = Tfc.Tfc';
73+
yrfc_ = Tfc.T(:,1);
74+
Tfc_ = Tfc_ - mean(Tfc_(yrfc_>=1860 & yrfc_<=1900)) + S.MH.Tm1860_1900(i);
75+
Tfc_ = Tfc_(yrfc_>=fyr);
76+
yrfc_ = yrfc_(yrfc_>=fyr);
77+
end
78+
79+
FC = calc_fc(S,Tfc_,i_P,i);
80+
sl = [sl;FC.sl];
81+
Tcf = [Tcf;Tfc_];
82+
T01 = [T01;FC.T01];
83+
84+
end
85+
end
86+
87+
fprintf('\n')
88+
89+
S.proj.sl = sl;
90+
S.proj.Tcf = Tcf;
91+
S.proj.T01 = T01;
92+
S.proj.time = yrfc_;
93+
94+
if lyr<=2000
95+
% Calculate likelihoods & Differences
96+
L = calc_lik(S);
97+
S.proj.Lik = L.Lik;
98+
S.proj.Lik_yrly = L.Lik_yrly;
99+
S.proj.Diff = L.Diff;
100+
S.proj.Diff_yrly = L.Diff_yrly;
101+
S.proj.Perc = L.Perc;
102+
S.proj.Perc_yrly = L.Perc_yrly;
103+
S.proj.slFactual = L.sl;
104+
end
105+
if strcmp(CFtemp, 'CMIP5')
106+
S.proj.column_CMIP5 = column;
107+
end
108+
if ~output_all
109+
S_ = S.proj;
110+
S_.settings = S.settings;
111+
S = S_;
112+
end
113+
114+
end
115+
116+
function fc = calc_fc(S,Tfc,i_P,i)
117+
118+
set = S.settings;
119+
120+
T01 = S.MH.T01_1900(i);
121+
if strcmp(S.settings.model,'TwoTau')
122+
T02 = S.MH.T02_1900(i);
123+
end
124+
125+
if strcmp(S.settings.model,'CRdecay')
126+
c = S.MH.c_1900(i);
127+
else
128+
c = S.MH.Params(i_P,3);
129+
end
130+
131+
a1 = S.MH.Params(i_P,1);
132+
a2 = S.MH.Params(i_P,2);
133+
tau1 = S.MH.Params(i_P,4);
134+
tau2 = S.MH.Params(i_P,5);
135+
136+
nyrs = max(size(Tfc));
137+
138+
for ii = 1:nyrs-1
139+
% Compute here T0(sample,year) as the tau-year memory of temp(j)
140+
T01(ii+1) = T01(ii) + (1/tau1)*(Tfc(ii) - T01(ii));
141+
if strcmp(set.model,'TwoTau')
142+
T02(ii+1) = T02(ii) + (1/tau2)*(Tfc(ii) - T02(ii));
143+
elseif strcmp(set.model,'CRdecay')
144+
c(ii+1) = c(ii) * (1-1/tau2);
145+
end
146+
end
147+
if strcmp(set.model,'TwoTau')
148+
dsl = a1*(Tfc - T01) + a2*(Tfc - T02);
149+
fc.T02 = T02;
150+
elseif strcmp(set.model,'ConstRate')
151+
dsl = c + a1*(Tfc - T01);
152+
elseif strcmp(set.model,'CRdecay')
153+
dsl = c + a1*(Tfc - T01);
154+
elseif strcmp(set.model,'CRovTau')
155+
dsl = c/tau1 + a1*(Tfc - T01);
156+
elseif strcmp(set.model,'simpel')
157+
dsl = a1*(Tfc - T01);
158+
end
159+
sl = cumsum(dsl);
160+
sl = sl-sl(1);
161+
162+
fc.sl = sl;
163+
fc.T01 = T01;
164+
end
165+
166+
function L = calc_lik(S)
167+
168+
% which SL to take as factual
169+
fact_SL = 'calib'; % 'calib': semi-emp calibration
170+
% 'prox': proxy data
171+
% 'recalc': semi-empiricaly recalculate with factual temperature
172+
%---- SL projection under counterfactual T
173+
sl = S.proj.sl;
174+
yr = S.proj.time;
175+
176+
%---- SL to calculate likelihood against: either Proxy or semi emp. calibration
177+
if strcmp(fact_SL,'calib')
178+
yrPrx = 1900:2000;
179+
slPrx = S.MH.sl_2000;
180+
elseif strcmp(fact_SL,'recalc')
181+
Tfc = S.data.temp.obs';
182+
yrPrx = S.data.temp.yr;
183+
ix = yrPrx>=yr(1) & yrPrx<=yr(end);
184+
fc = calc_fc(S,Tfc(ix),yrPrx(ix),[50,16,84]);
185+
slPrx = fc.sl;
186+
yrPrx = yrPrx(ix);
187+
end
188+
% cut proxy to same period as T scenario
189+
ix = yrPrx>=yr(1) & yrPrx<=yr(end);
190+
yrPrx = yrPrx(ix);
191+
slPrx = slPrx(:,ix);
192+
for i = 1:size(sl,1)%count sample
193+
Diff_(i,:) = -sl(i,:)+slPrx(i,:);
194+
Perc_(i,:) = (sl(i,:)./slPrx(i,:))*100;
195+
for j = 1:size(sl,2)%count years
196+
L_(i,j) = sum(sl(:,j)>=slPrx(i,j))/size(sl,1);
197+
end
198+
end
199+
200+
% index to give data the spacing of SL & T
201+
ix = ismember(yrPrx,S.MH.yr_);
202+
ix(yrPrx==yr(1))=true;
203+
ix(yrPrx==yr(end))=true;
204+
205+
Lik(:,1) = yrPrx; % the year at which the likelihoods are calculated
206+
Lik(:,2) = mean(L_,1);
207+
L.Lik_yrly = Lik;
208+
L.Lik = Lik(ix,:);
209+
210+
D = prctile(Diff_,[50 17 83 5 95]);
211+
L.Diff_yrly(:,1) = yrPrx;
212+
L.Diff_yrly(:,2:6) = D';
213+
L.Diff = L.Diff_yrly(ix,:);
214+
215+
P = prctile(Perc_,[50 17 83 5 95]);
216+
P(:,1) = 100;
217+
L.Perc_yrly(:,1) = yrPrx;
218+
L.Perc_yrly(:,2:6) = P';
219+
L.Perc = L.Perc_yrly(ix,:);
220+
221+
L.sl(:,1) = yrPrx(ix);
222+
L.sl(:,2:6) = prctile(slPrx(:,ix),[50 17 83 5 95])';
223+
224+
end

‎Calc_SESLProjection.m‎

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,152 @@
1+
function O = Calc_SESLProjection(S,rcps,dig, NumTmag)
2+
3+
% Calculate sea-level projections under different RCP scenarios.
4+
%
5+
% O = Calc_SESLProjection(S,rcps,dig, NumTmag)
6+
%
7+
% INPUT:
8+
% - S -> Output of Calc_SESL_Prc.m
9+
% - rcps -> Cell array filled with strings. RCPs op use: {'RCP3PD', 'RCP45', 'RCP85'}
10+
% - dig -> number of digits to use, else []
11+
% - NumTmag -> Number of MAGICC RCP realizations to take max&default == 600
12+
%
13+
% OUTPUT: Sea-level time-series for each selected RCP (Psl) and year dates (Pslyr).
14+
% The columns in Psl give the percentiles Prc as defined below
15+
16+
usesingle = true; % to save memory
17+
18+
timeprojected = [2000 2100];
19+
20+
Prc = [5 17 50 83 95]; % Percentiles (68.27=1sigma)
21+
22+
folder = 'Data';
23+
24+
if usesingle
25+
SFnc = @(x)single(x);
26+
else
27+
SFnc = @(x)x;
28+
end
29+
30+
if isempty(NumTmag)
31+
NumTmag=600;
32+
end
33+
Tnum = S.settings.Tnum;
34+
sample = S.settings.sample;
35+
36+
if ~isempty(dig)
37+
digits(dig)
38+
end
39+
a = SFnc(S.MH.Params(:,1));
40+
a = reshape(repmat(a',Tnum,1),Tnum*sample,1);
41+
a = repmat(a,NumTmag,1);
42+
if strcmp(S.settings.model,'CRdecay')
43+
c_2000 = SFnc(S.MH.c_2000);
44+
c_2000 = repmat(c_2000,NumTmag,1);
45+
46+
tau_c = SFnc(S.MH.Params(:,5));
47+
tau_c = reshape(repmat(tau_c',Tnum,1),Tnum*sample,1);
48+
tau_c = repmat(tau_c,NumTmag,1);
49+
elseif strcmp(S.settings.model,'simpel')
50+
c = SFnc(zeros(size(S.MH.Params(:,3))));
51+
c = reshape(repmat(c',Tnum,1),Tnum*sample,1);
52+
c = repmat(c,NumTmag,1);
53+
elseif strcmp(S.settings.model,'CRovTau')
54+
c = SFnc(S.MH.Params(:,3)./S.MH.Params(:,4));
55+
c = reshape(repmat(c',Tnum,1),Tnum*sample,1);
56+
c = repmat(c,NumTmag,1);
57+
elseif strcmp(S.settings.model,'ConstRate')
58+
c = SFnc(S.MH.Params(:,3));
59+
c = reshape(repmat(c',Tnum,1),Tnum*sample,1);
60+
c = repmat(c,NumTmag,1);
61+
end
62+
tau = SFnc(S.MH.Params(:,4));
63+
tau = reshape(repmat(tau',Tnum,1),Tnum*sample,1);
64+
tau = repmat(tau,NumTmag,1);
65+
66+
T0_2000 = SFnc(S.MH.T01_2000);
67+
T0_2000 = repmat(T0_2000,NumTmag,1);
68+
69+
if strcmp(S.settings.model,'TwoTau')
70+
a2 = SFnc(S.MH.Params(:,2));
71+
a2 = reshape(repmat(a2',Tnum,1),Tnum*sample,1);
72+
a2 = repmat(a2,NumTmag,1);
73+
T02_2000 = SFnc(S.MH.T02_2000);
74+
T02_2000 = repmat(T02_2000,NumTmag,1);
75+
end
76+
77+
Toff1 = SFnc(S.MH.Tm1970_2000);
78+
Toff1 = repmat(Toff1,NumTmag,1);
79+
settings = S.settings;
80+
clear S
81+
82+
for j = 1:length(rcps) % Count RCPs
83+
84+
sl = SFnc(zeros(size(a)));
85+
86+
if strcmp(settings.model,'CRdecay')
87+
c = c_2000;
88+
end
89+
T0 = T0_2000;
90+
if strcmp(settings.model,'TwoTau')
91+
T02 = T0_2000;
92+
end
93+
% load MAGICC realizations of RCPS: TIME & DATA
94+
if usesingle
95+
load(fullfile(folder, [cell2mat(rcps(j)), '_single'])); % load rcp
96+
else
97+
load(fullfile(folder, cell2mat(rcps(j)))); % load rcp
98+
end
99+
100+
DATA = DATA(:,1:NumTmag);
101+
102+
Toff2 = repmat(mean(DATA(TIME>=1970 & TIME<=2000,:)),size(DATA,1),1);
103+
104+
DATA = DATA-Toff2;
105+
106+
DATA = DATA(TIME>=timeprojected(1) & TIME<=timeprojected(end),:);
107+
108+
clear TIME T_90prc colheaders textdata Toff2
109+
110+
fprintf('\t %1.0f years of %1s are now calculated: \n',size(DATA,1),cell2mat(rcps(j)));
111+
112+
for i_yr = 1:size(DATA,1); % count years to be projected
113+
fprintf('%1.0f \t',i_yr);
114+
115+
Tfc = DATA(i_yr,:);
116+
Tfc = reshape(repmat(Tfc,Tnum*sample,1),NumTmag*Tnum*sample,1);
117+
Tfc = Tfc + Toff1;
118+
119+
if i_yr>1 % so it will be relative to 2100
120+
if ~strcmp(settings.model,'TwoTau')
121+
sl = sl + a .* (Tfc-T0) + c;
122+
else
123+
sl = sl + a .* (Tfc-T0) + a2 .* (Tfc-T02);
124+
end
125+
end
126+
127+
T0 = T0 + (1./tau).*(Tfc-T0);
128+
if strcmp(settings.model,'TwoTau')
129+
T02 = T02 + (1./tau).*(Tfc-T0);
130+
end
131+
if strcmp(settings.model,'CRdecay')
132+
c = c .* (ones(size(c))-1./tau_c);
133+
end
134+
clear Tfc
135+
136+
O.(cell2mat(rcps(j))).Psl(i_yr,:) = prctile(sl,Prc);
137+
O.(cell2mat(rcps(j))).Pslyr=timeprojected(1):timeprojected(end);
138+
end
139+
fprintf('\n')
140+
end
141+
end
142+
143+
144+
145+
% sample Tnum 600
146+
% P1 T1 TR1
147+
% P1 T2 TR1
148+
% P1 T3 TR1
149+
%
150+
% P2 T1 TR1
151+
% P2 T2 TR1
152+
% P2 T3 TR1

0 commit comments

Comments
 (0)