diff --git a/Calc_SESLConterfact.m b/Calc_SESLConterfact.m new file mode 100755 index 0000000..83360a3 --- /dev/null +++ b/Calc_SESLConterfact.m @@ -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 +% - 'CMIP5' - 'historicalNat_CMIP35_tas' single simu (data column given as input) from Gareth Jones +% - '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=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 \ No newline at end of file diff --git a/Calc_SESLProjection.m b/Calc_SESLProjection.m new file mode 100755 index 0000000..bd4a435 --- /dev/null +++ b/Calc_SESLProjection.m @@ -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 diff --git a/Calc_SESL_Prc.m b/Calc_SESL_Prc.m new file mode 100755 index 0000000..8e15793 --- /dev/null +++ b/Calc_SESL_Prc.m @@ -0,0 +1,388 @@ +function P = Calc_SESL_Prc(T, num,folder) + +% Calculates percentiles of sea-level and temperature (which were previously +% stored by Calc_SL_from_Param.m) alongside with other values important for +% projections or counterfactual scenarios +% +% P = Calc_SESL_Prc(T, num,folder) +% +% INPUT: +% - T -> string, name of the temperature used for calibration +% - num -> number of files previously stored by Calc_SL_from_Param.m +% - folder -> string, name of the folder where Calc_SL_from_Param.m stored +% away files +% +% OUTPUT: +% - settings -> settings that were also in previous outputs +% - data -> original sea-level and temperature data used for calibration +% - MH -> many items related to temperature T, equilibrium +% temperature T0, the parameter c and sea level SL . Examples: +% - T01_2000: the values of T01 in the year 2000 CE +% - psl: percentiles of SL as prescribed by Prc below +% - T01_1000BC: the values of T01 in the year 1000 BC +% - Tm1801_1900: the mean value of T between 1801-1900 CE +% - SLpo1800_2000: 1st degree polynomial (polyfit(x,y,1)) of SL +% between 1800-2000 CE + + calc_Full_T0 = true; % another function at the end of the script to combine all 'num' files for full T0 sample + + Prc = [5 17 50 83 95]; % Percentiles to calculate + + S = []; % EM added + load(fullfile(folder, [T, '_1'])); + P.settings = S.settings; + P.data = S.data; + P.MH.Params = S.MH.Params; + P.MH.alpha = S.MH.alpha; + P.MH.yrT = S.MH.yrT; + P.MH.yr = S.MH.yrsl; + if strcmp(P.settings.model,'CRdecay') + P.MH.yrc = S.MH.yrc; + end + % find proper divisor of the number of sl length + DivSL = []; + for x = 1:size(S.MH.sl,2) + if mod(size(S.MH.sl,2),x)==0 + DivSL = [DivSL,x]; + end + end + Nyst = DivSL(2); % choose smallest so it won't run out of memory + yst = size(S.MH.sl,2)/Nyst; + clear S + + % Define Output + sl_ = []; + psl = []; + T01_2000 = []; + T01_1900 = []; + T01_501 = []; + T01m501_700 = []; + T01_1000BC = []; + T01_2000BC = []; + if strcmp(P.settings.model,'TwoTau') + T02_2000 = []; + T02_1900 = []; + T02_501 = []; + T02m501_700 = []; + T02_1000BC = []; + T02_2000BC = []; + end + if strcmp(P.settings.model,'CRdecay') + c_2000 = []; + c_1900 = []; + c_501 = []; + end + sl_20C = []; + Tpo500_1800 = []; + Tm1801_1900 = []; + Tm500_1800 = []; + Tm501_700 = []; + Tm2000_1800BC = []; + Tm1850_1900 = []; + Tm1860_1900 = []; + Tm1970_2000 = []; + SLpo800_1800 = []; + SLpo1800_2000 = []; + SLpo0_1700 = []; + SLpo0_400 = []; + SLpo400_800 = []; + SLpo800_1200 = []; + SLpo1200_1600 = []; + SLpo1600_1800 = []; + SLpo1800_1900 = []; + SLpo1900_2000 = []; + + fprintf('%1.0f',zeros(1,Nyst*num));fprintf('\n'); + for i = 1:Nyst + for j = 1:num + fprintf('0') + load(fullfile(folder, [T, '_', num2str(j)])); + yrs = S.data.temp.yrs; + % sl ---------------------------------------------------------- + sl_ = [sl_; S.MH.sl(:,(i-1)*yst+1:i*yst)]; +% + if i==1 + % T01(2000) ------------------------------------------------- + lyr =2000; + if sum(S.MH.yrT==lyr)==0 + T01_2000 = [T01_2000; S.MH.T01(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T01_2000 = [T01_2000;S.MH.T01(:,S.MH.yrT == lyr)]; + end + % T01(1900) ------------------------------------------------- + lyr =1900; + if sum(S.MH.yrT==lyr)==0 + T01_1900 = [T01_1900; S.MH.T01(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T01_1900 = [T01_1900;S.MH.T01(:,S.MH.yrT == lyr)]; + end + % T01(501) -------------------------------------------------- + lyr =501; + if sum(S.MH.yrT==lyr)==0 + T01_501 = [T01_501; S.MH.T01(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T01_501 = [T01_501;S.MH.T01(:,S.MH.yrT == lyr)]; + end + %T01(1000BC) ------------------------------------------------ + lyr =-1000; + if sum(S.MH.yrT==lyr)==0 + T01_1000BC = [T01_1000BC; S.MH.T01(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T01_1000BC = [T01_1000BC;S.MH.T01(:,S.MH.yrT == lyr)]; + end + %T01(2000BC) ------------------------------------------------ + lyr =-2000; + if sum(S.MH.yrT==lyr)==0 + T01_2000BC = [T01_2000BC; S.MH.T01(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T01_2000BC = [T01_2000BC;S.MH.T01(:,S.MH.yrT == lyr)]; + end + %T01(501-700) --------------------------------------------- + ix = S.MH.yrT>=501 & S.MH.yrT<=700; + T01m501_700 = [T01m501_700; mean(S.MH.T01(:,ix),2)]; + if strcmp(P.settings.model,'TwoTau') + % T02(2000) ------------------------------------------------- + lyr =2000; + if sum(S.MH.yrT==lyr)==0 + T02_2000 = [T02_2000; S.MH.T02(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T02_2000 = [T02_2000;S.MH.T02(:,S.MH.yrT == lyr)]; + end + % T02(1900) ------------------------------------------------- + lyr =1900; + if sum(S.MH.yrT==lyr)==0 + T02_1900 = [T02_1900; S.MH.T02(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T02_1900 = [T02_1900;S.MH.T02(:,S.MH.yrT == lyr)]; + end + % T02(501) -------------------------------------------------- + lyr =501; + if sum(S.MH.yrT==lyr)==0 + T02_501 = [T02_501; S.MH.T02(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T02_501 = [T02_501;S.MH.T02(:,S.MH.yrT == lyr)]; + end + %T02(1000BC) ------------------------------------------------ + lyr =-1000; + if sum(S.MH.yrT==lyr)==0 + T02_1000BC = [T02_1000BC; S.MH.T02(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T02_1000BC = [T02_1000BC;S.MH.T02(:,S.MH.yrT == lyr)]; + end + %T02(2000BC) ------------------------------------------------ + lyr =-2000; + if sum(S.MH.yrT==lyr)==0 + T02_2000BC = [T02_2000BC; S.MH.T02(:,S.MH.yrT>=lyr-ceil((yrs-1)/2) & S.MH.yrT<=lyr+floor((yrs-1)/2))]; + else + T02_2000BC = [T02_2000BC;S.MH.T02(:,S.MH.yrT == lyr)]; + end + %T02(501-700) --------------------------------------------- + ix = S.MH.yrT>=501 & S.MH.yrT<=700; + T02m501_700 = [T02m501_700; mean(S.MH.T02(:,ix),2)]; + end + if strcmp(P.settings.model,'CRdecay') + % c(2000) ----------------------------------------------- + lyr =2000; + if sum(S.MH.yrT==lyr)==0 + c_2000 = [c_2000; S.MH.c(:,S.MH.yrc>=lyr-ceil((yrs-1)/2) & S.MH.yrc<=lyr+floor((yrs-1)/2))]; + else + c_2000 = [c_2000;S.MH.c(:,S.MH.yrc == lyr)]; + end + % c(1900) ----------------------------------------------- + lyr =1900; + if sum(S.MH.yrT==lyr)==0 + c_1900 = [c_1900; S.MH.c(:,S.MH.yrc>=lyr-ceil((yrs-1)/2) & S.MH.yrc<=lyr+floor((yrs-1)/2))]; + else + c_1900 = [c_1900;S.MH.c(:,S.MH.yrc == lyr)]; + end + % c(501) ------------------------------------------------ + lyr =501; + if sum(S.MH.yrT==lyr)==0 + c_501 = [c_501; S.MH.c(:,S.MH.yrc>=lyr-ceil((yrs-1)/2) & S.MH.yrc<=lyr+floor((yrs-1)/2))]; + else + c_501 = [c_501;S.MH.c(:,S.MH.yrc == lyr)]; + end + end + % sl(1900-2000) --------------------------------------------- + sl_20C_ = S.MH.sl(:,S.MH.yrsl>=1900 & S.MH.yrsl<=2000); + sl_20C = [sl_20C;sl_20C_ - repmat(sl_20C_(:,1),1,size(sl_20C_,2))]; + % T linear rate 500-1800 ------------------------------------ + ix = S.MH.yrT>=500 & S.MH.yrT<=1800; + for k = 1: size(S.MH.temp,1) + Tpo500_1800_(k,:) = polyfit(S.MH.yrT(ix),S.MH.temp(k,ix)',1); + end + Tpo500_1800 = [Tpo500_1800;Tpo500_1800_]; clear Tpo500_1800_; + % mean T(1801-1900) ----------------------------------------- + ix = S.MH.yrT>=1801 & S.MH.yrT<=1900; + Tm1801_1900 = [Tm1801_1900; mean(S.MH.temp(:,ix),2)]; + % mean T(500-1800) ------------------------------------------ + ix = S.MH.yrT>=501 & S.MH.yrT<=1800; + Tm500_1800 = [Tm500_1800; mean(S.MH.temp(:,ix),2)]; + % mean T(500-700) ------------------------------------------- + ix = S.MH.yrT>=501 & S.MH.yrT<=700; + Tm501_700 = [Tm501_700; mean(S.MH.temp(:,ix),2)]; + % mean T(-2000--1800) --------------------------------------- + ix = S.MH.yrT>=-2000 & S.MH.yrT<=-1800; + Tm2000_1800BC = [Tm2000_1800BC; mean(S.MH.temp(:,ix),2)]; + % mean T(1850-1900) ----------------------------------------- + ix = S.MH.yrT>=1850 & S.MH.yrT<=1900; + Tm1850_1900 = [Tm1850_1900; mean(S.MH.temp(:,ix),2)]; + %mean T(1860-1900) ---------------------------------------- + ix = S.MH.yrT>=1860 & S.MH.yrT<=1900; + Tm1860_1900 = [Tm1860_1900; mean(S.MH.temp(:,ix),2)]; + % mean T(1970-2000) ----------------------------------------- + ix = S.MH.yrT>=1970 & S.MH.yrT<=2000; + Tm1970_2000 = [Tm1970_2000; mean(S.MH.temp(:,ix),2)]; + % SL linear rate 0-1700 ------------------------------------- + ix = S.MH.yrsl>=0 & S.MH.yrsl<=1700; + for k = 1: size(S.MH.temp,1) + SLpo0_1700_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo0_1700 = [SLpo0_1700;SLpo0_1700_]; clear SLpo0_1700_; + % SL linear rate 800-1800 ----------------------------------- + ix = S.MH.yrsl>=800 & S.MH.yrsl<=1800; + for k = 1: size(S.MH.temp,1) + SLpo800_1800_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo800_1800 = [SLpo800_1800;SLpo800_1800_]; clear SLpo800_1800_; + % SL linear rate 1800-2000 ---------------------------------- + ix = S.MH.yrsl>=1800 & S.MH.yrsl<=2000; + for k = 1: size(S.MH.temp,1) + SLpo1800_2000_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo1800_2000 = [SLpo1800_2000;SLpo1800_2000_]; clear SLpo1800_2000_; + % SL linear rate 0-400 -------------------------------------- + ix = S.MH.yrsl>=0 & S.MH.yrsl<=400; + for k = 1: size(S.MH.temp,1) + SLpo0_400_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo0_400 = [SLpo0_400;SLpo0_400_]; clear SLpo0_400_; + % SL linear rate 400-800 ------------------------------------ + ix = S.MH.yrsl>=400 & S.MH.yrsl<=800; + for k = 1: size(S.MH.temp,1) + SLpo400_800_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo400_800 = [SLpo400_800;SLpo400_800_]; clear SLpo400_800_; + % SL linear rate 800-1200 ----------------------------------- + ix = S.MH.yrsl>=800 & S.MH.yrsl<=1200; + for k = 1: size(S.MH.temp,1) + SLpo800_1200_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo800_1200 = [SLpo800_1200;SLpo800_1200_]; clear SLpo800_1200_; + % SL linear rate 1200-1600 ---------------------------------- + ix = S.MH.yrsl>=1200 & S.MH.yrsl<=1600; + for k = 1: size(S.MH.temp,1) + SLpo1200_1600_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo1200_1600 = [SLpo1200_1600;SLpo1200_1600_]; clear SLpo1200_1600_; + % SL linear rate 1600-1800 ---------------------------------- + ix = S.MH.yrsl>=1600 & S.MH.yrsl<=1800; + for k = 1: size(S.MH.temp,1) + SLpo1600_1800_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo1600_1800 = [SLpo1600_1800;SLpo1600_1800_]; clear SLpo1600_1800_; + % SL linear rate 1800-1900 ---------------------------------- + ix = S.MH.yrsl>=1800 & S.MH.yrsl<=1900; + for k = 1: size(S.MH.temp,1) + SLpo1800_1900_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo1800_1900 = [SLpo1800_1900;SLpo1800_1900_]; clear SLpo1800_1900_; + % SL linear rate 1900-2000 ---------------------------------- + ix = S.MH.yrsl>=1900 & S.MH.yrsl<=2000; + for k = 1: size(S.MH.temp,1) + SLpo1900_2000_(k,:) = polyfit(S.MH.yrsl(ix),S.MH.sl(k,ix),1); + end + SLpo1900_2000 = [SLpo1900_2000;SLpo1900_2000_]; clear SLpo1900_2000_; + end + clear S + end + psl_ = prctile(sl_,Prc); + if size(psl_,1) == length(Prc) + psl = [psl, prctile(sl_,Prc)]; + elseif size(psl_,2) == length(Prc); + psl = [psl, prctile(sl_,Prc)']; + end + sl(:,(i-1)*yst+1:i*yst) = sl_; + sl_ = []; + end + + P.MH.psl = psl; + P.MH.sl = sl; + P.MH.sl_20C = sl_20C; + P.MH.T01_2000 = T01_2000; + P.MH.T01_1900 =T01_1900; + if size(T01_501,2)>1 + P.MH.T01_501 =T01_501(:,2); + else + P.MH.T01_501 =T01_501; + end + P.MH.T01_1000BC = T01_1000BC; + P.MH.T01_2000BC = T01_2000BC; + P.MH.T01m501_700 = T01m501_700; + if strcmp(P.settings.model,'TwoTau') + P.MH.T02_2000 = T02_2000; + P.MH.T02_1900 =T02_1900; + if size(T02_501,2)>1 + P.MH.T02_501 =T02_501(:,2); + else + P.MH.T02_501 =T02_501; + end + P.MH.T02_1000BC = T02_1000BC; + P.MH.T02_2000BC = T02_2000BC; + P.MH.T02m501_700 = T02m501_700; + end + if strcmp(P.settings.model,'CRdecay') + P.MH.c_2000 = c_2000; + P.MH.c_1900 = c_1900; + if size(c_501,2)>1 + P.MH.c_501 = c_501(:,2); + else + P.MH.c_501 = c_501; + end + end + P.MH.Tpo500_1800 = Tpo500_1800; + P.MH.Tm1801_1900 = Tm1801_1900; + P.MH.Tm500_1800 = Tm500_1800; + P.MH.Tm501_700 = Tm501_700; + P.MH.Tm2000_1800BC = Tm2000_1800BC; + P.MH.Tm1850_1900 = Tm1850_1900; + P.MH.Tm1860_1900 = Tm1860_1900; + P.MH.Tm1970_2000 = Tm1970_2000; + P.MH.SLpo0_1700 = SLpo0_1700; + P.MH.SLpo800_1800 = SLpo800_1800; + P.MH.SLpo1800_2000 = SLpo1800_2000; + P.MH.SLpo0_400 = SLpo0_400 ; + P.MH.SLpo400_800 = SLpo400_800 ; + P.MH.SLpo800_1200 = SLpo800_1200 ; + P.MH.SLpo1200_1600 = SLpo1200_1600 ; + P.MH.SLpo1600_1800 = SLpo1600_1800 ; + P.MH.SLpo1800_1900 = SLpo1800_1900 ; + P.MH.SLpo1900_2000 = SLpo1900_2000 ; + P.MH.sl = sl; + +if calc_Full_T0 + for i = 1:num + load(fullfile(folder, [T, '_', num2str(i)])) + if i==1 + Num = size(S.MH.sl,1); + P.MH.T0 = ones(Num*num, size(S.MH.T01,2)); + if strcmp(P.settings.model,'TwoTau') + P.MH.T02 = ones(Num*num, size(S.MH.T02,2)); + end + end + S.MH.Params = []; + S.MH.alpha= []; + S.MH.temp= []; + S.MH.c= []; + S.MH.sl= []; + S.MH.yrsl= []; + S.MH.yrc= []; + S.MH.yrT= []; + P.MH.T0((i-1)*Num+1:i*Num,:) = S.MH.T01; + if strcmp(P.settings.model,'TwoTau') + P.MH.T02((i-1)*Num+1:i*Num,:) = S.MH.T02; + end + end +end +fprintf('\n') +end \ No newline at end of file diff --git a/Calc_SL_from_Param.m b/Calc_SL_from_Param.m new file mode 100755 index 0000000..6964286 --- /dev/null +++ b/Calc_SL_from_Param.m @@ -0,0 +1,128 @@ +function Calc_SL_from_Param(S,storeNum,Meth,folder) + +% Calculates time series of sea level, temperatures and c from the +% parameter set found in SESL.m +% +% Calc_SL_from_Param(S,storeNum,Meth,folder) +% +% INPUT: +% - S -> Output structure from SESL.m +% - storeNum -> storeNum time series of sea level, temperature, equilibrium +% temperature and c are calculated, numbered & saved in 'folder'. +% If there is no memory issue storeNum can be set to sample*Tnum +% - Meth -> The method used befor by SESL.m +% 'MH' for Metropolis Hastings, +% 'Slice' for Slice sampling (not working correctly) +% - folder -> folder to save samples with sizes of 'storeNum' as prescribed below +% +% OUTPUT: a number of output structures S are saved and numbered in the folder +% defined above. Each of these outputs is similar to the SESL.m output +% only that the structure S.MH or S.Slice (dependent on the input +% Method 'Meth') will be more detailed. +% Besides Params & alpha, S.MH/Slice includes time series of T01 +% (and T02 if needed), c, temp (temperature), sl (sea level) and +% the according year counts yrT, yrc, yrsl. The number of samples these +% time series include is defined by storeNum and is usually smaller +% than the full set of samples, that is why these outputs later +% need to be combined by Calc_SESL_Prc.m. This is a memory issue + +Tnum = []; % calc. Tnum random temperatures and sea-levels for each param. set and take all. + % if isempty(Tnum) take as many random temps. as prescribend in settings + +dat = S.data; +set = S.settings; +if strcmp(Meth,'SimAnn') + smpl = S.(Meth).OptParam; +else + smpl = S.(Meth).Params; +end + +pb = 100; % show a progress bar, progressing every 100 samples +if size(smpl,1) dSL/dt = a * (T(t)-T0(t)) +% TwoTau => dSL/dt = a1 * (T(t)-T0_1(t)) + a2 * (T(t)-T0_2(t)) +% ConstRate => dSL/dt = c + a * (T(t)-T0(t)) +% CRdecay => dSL/dt = c(t) + a * (T(t)-T0(t)) +% CRovTau => dSL/dt = c/tau + a * (T(t)-T0(t)) +% -> dT0(t)/dt = (T-T0)/tau +% -> dc(t)/dt = -c(t)/tau +defmodel = 'CRdecay'; + valmodel = {'ConstRate' 'CRdecay' 'CRovTau' 'TwoTau' 'simpel'}; + checkmodel = @(x) any(validatestring(x,valmodel)); + addParamValue(P,'model',defmodel,checkmodel) + +% Prior distribution for Parameters -> important for Likelihood +defall_flat_priors = false; + addParamValue(P,'all_flat_priors',defall_flat_priors,@islogical); +% adjust the prior in check_settings, according to the model used +defAdjustPrior = false; + addParamValue(P,'AdjustPrior',defAdjustPrior,@islogical); + +defStartDistr = [0.5, 0.3, 0, log10(50), log10(100) , 0]; % Start distribution for simulated annealing + addParamValue(P,'StartDistr',defStartDistr,@isnumeric); + +defa1_prior = {'uniform' 0 2}; + addParamValue(P,'a1_prior',defa1_prior,@iscell); +defa2_prior = {'uniform' 0 2}; % a2(if model=='TwoTau'), b(if useb==true) + addParamValue(P,'a2_prior',defa2_prior,@iscell); +deftau1_prior = {'uniform',log(30),log(3000)}; % log uniform distribution + addParamValue(P,'tau1_prior',deftau1_prior,@iscell); +deftau2_prior = {'uniform',log(1000),log(20000)}; % log uniform distribution + addParamValue(P,'tau2_prior',deftau2_prior,@iscell); +defc_prior = {'flat', [],[]}; % flat prior + addParamValue(P,'c_prior',defc_prior,@iscell); +defc2000_prior = {'uniform', -.2,.2}; + addParamValue(P,'c2000_prior',defc2000_prior,@iscell); +defT01st_prior = {'uniform',-.6,.6}; % T0 at the beginning of T if it is before 500CE + addParamValue(P,'T01st_prior',defT01st_prior,@iscell); +defT02nd_prior = {'normal',0,.2}; % T0(500-700) prior + addParamValue(P,'T02nd_prior',defT02nd_prior,@iscell); + +defTauLogUniform = true; % true => Tau priors are loguniform, else uniform + addParamValue(P,'TauLogUniform',defTauLogUniform,@islogical); + +% optimize T0(1) (true) or set to mean(T(T0period)) +defOptimT0 = true; + addParamValue(P,'OptimT0',defOptimT0,@islogical); + +% Jumping distribution size +defJumpDist = .005; + addParamValue(P,'JumpDist',defJumpDist,@isnumeric); + +% Number of ar1 temperatures to evaluate for likelihood +defTnum= 100; + addParamValue(P,'Tnum',defTnum,@isnumeric); + +% Sea Level data to use +defSL_dat ='GLMW-1ts'; % which SL proxy should be used? Add dummy data to folder '\Data' + valSL_dat = {'GLMW','GLMW-1amp1ts','GLMW-1ts','GLMW-Gr','GLMW-NC'}; + checkSL_dat = @(x) any(validatestring(x,valSL_dat)); + addParamValue(P,'SL_dat',defSL_dat,checkSL_dat) + +defCalibperiod = -1000:2010; % years to include for model calibration. KE11[-120 - 2000 AD], Mann Temp. [500 - 2006 AD] + addParamValue(P,'calibperiod',defCalibperiod,@isnumeric); +defPeriod = -1000:2010; % years of data to include for SL calculation + addParamValue(P,'period',defPeriod,@isnumeric); +defBaseperiod = 1400:1800; % Period for which simulated sl and data are normed + addParamValue(P,'baseperiod',defBaseperiod,@isnumeric); +defT0period =-2000:-1800; % period to initialise T0(1) + addParamValue(P,'T0period',defT0period,@isnumeric); +defT0length = 0; % if==0: T0period is used, else T0(1) = mean(T(first T0length yrs)) + addParamValue(P,'T0length',defT0length,@isnumeric); + +defUseMarT0 = true; % use Marcott temperature to calc. T0 until 'T_data' starts + addParamValue(P,'UseMarT0',defUseMarT0,@islogical); +defT0temp_level = 100; % number of years over which to level Marcott and T_data if UseMarT0==true + addParamValue(P,'T0temp_level',defT0temp_level,@isnumeric) + +defOptHo = true; % for each MC sample, calculate the optimal offset between simulated and calibration sea level via MLS + addParamValue(P,'optHo',defOptHo,@islogical); + +defT_data = 'Mann09_11'; % Choose temperature data for calibration + valT_data = {'Mann08_eiv','Mann08_cps','Mann09','Mann09_11','Marcott13_RegEM-HC3_20','PAGES2k_13'}; + checkT_data = @(x) any(validatestring(x,valT_data)); + addParamValue(P,'T_data',defT_data,checkT_data); + % 'Mann08_eiv' -> Mann et al. PNAS 2008 eiv reco (default) + % 'Mann08_cps' -> Mann et al. PNAS 2008 cps reco + % 'Mann09' -> From Scott Rutherford from Mann et al. 2009 Science + % 'Mann09_11' -> as above but 11 year averages. + % 'Marcott13_RegEM'-> Marcott et al. 2013 + HadCRUT3v (20yr mean) starting 9350BC + % 'PAGES2k13' -> Pages2K -9-2000AD + % 'Pages2k' -> Pages2K -10-1909AD, from 1910-2010AD mean(Pages2k(1850-1909AD)) + +defBurnin = 1000; % number of monte carlo samples used in the burning-in or spin-up period + addParamValue(P,'burning',defBurnin,@isnumeric); +defNumSkip =100; % only every NumSkip sample will be selected because of autocorrelation + addParamValue(P,'NumSkip',defNumSkip,@isnumeric); +defSample = 5000; % number of monte carlo samples + addParamValue(P,'sample',defSample,@isnumeric); + +defT_err = 'ar1ts'; % How should the temperature variance be treated ? + valT_err = {'default','no','ar1','ar1ts'}; + checkT_err = @(x) any(validatestring(x,valT_err)); + addParamValue(P,'T_err',defT_err,checkT_err); + % 'default' -> T + random noise as in KE11 + % 'no' -> Don't add uncertainty + % 'ar1' -> T as AR(1) process with sig as 'default' + % 'ar1ts' -> AR(1) Parameter timescale == exp(-abs(t2-t1)/timescale) +deftau_ar1 = 10; % if Terr==ar1ts: timescale of ar1ts uncertainty above + % if Terr==ar1: AR(1) parameter + addParamValue(P,'tau_ar1',deftau_ar1,@isnumeric); +defRyrs = 10; % If T_err==default: Number of years for which to add same noise on temp. in calc_sl + addParamValue(P,'ryrs',defRyrs,@isnumeric); + % 10 -> Mann et al. default (Kemp et al. 2011) + % 20 -> Marcott et al. time series are 20yr means + % 30 -> PAGES2k time series are 30yr means +defTerrSc = 1; % scaling of temp. uncertainty + addParamValue(P,'TerrSc',defTerrSc,@isnumeric); + +defNormProb = false; % take into account the normfactor 1/sqrt(2*pi*|C|) + addParamValue(P,'NormProb',defNormProb,@islogical); + +defUseCov = true; % if existing, use covariance matrix of SL data to calc likelihood + addParamValue(P,'useCov',defUseCov,@islogical); +defCovShrink = 0; % Sea level covariance matrix (C) shrikage parameter f. C gets 'shrinked' to C' while C0 = diag(diag(C)): C' = f*C0 + (1-f)*C + addParamValue(P,'CovShrink',defCovShrink,@isnumeric); +defCovTau = 100; % if ~isnan(CovTau): take the elementwise product of the covariance and a tapering function exp(-delta(t)/CovTau) where CovTau is a time scale + addParamValue(P,'CovTau',defCovTau,@isnumeric); +defNoNegCov = true; % remove negative Cov Mat entries and replace by 0 + addParamValue(P,'NoNegCov',defNoNegCov,@islogical); + +defFac1 = 1; + addParamValue(P,'fac1',defFac1,@isnumeric); + % adjustment factor for autocorrelated NC sea level data as + % described in KE11 "Bayesian updating to estimate the model + % parameters" (default = 10). Gets multiplied to SL uncertainty + +parse(P,in{:}); + +set = P.Results; +end + diff --git a/LoadData_SESL.m b/LoadData_SESL.m new file mode 100755 index 0000000..fc5621e --- /dev/null +++ b/LoadData_SESL.m @@ -0,0 +1,78 @@ +function d = LoadData_SESL(set) + +% Load input temperature and sea-level data. +% +% d = LoadData_SESL(set) +% +% INPUT: general settings from DefineSettings_SESL.m +% +% OUTPUT: data structure + + + folder = 'Data'; + + %%%%%%%%%%%%%%%%%%%%%%% SEA LEVEL in cm %%%%%%%%%%%%%%%%%%%%%%% + + SL_dat = set.SL_dat; + + load(fullfile(folder,[SL_dat, '.mat'])); + + d.sea.proxy.yr = sl(:,1)'; + d.sea.proxy.sl = sl(:,2)'/10; + d.sea.proxy.slerr = sl(:,3)'/10; + + C = C/100; + C = C + eye(length(C))*eps; + + if set.useCov + % Scale Cov Mat and remove negative elements + if ~isnan(set.CovTau) && exist('C','var')==1 + Csc = exp(-abs(bsxfun(@minus,d.sea.proxy.yr,d.sea.proxy.yr')/set.CovTau)); + d.sea.proxy.C = C.*Csc; + else + d.sea.proxy.C = set.CovShrink*diag(diag(C))+ (1-set.CovShrink)*C; + end + if set.NoNegCov && exist('C','var')==1 + ix = d.sea.proxy.C<0; + d.sea.proxy.C(ix) = 0; + end + end + + % offset proxy SL data + offset = mean(d.sea.proxy.sl(d.sea.proxy.yr>=set.baseperiod(1) & d.sea.proxy.yr<=set.baseperiod(end))); % offset to baseperiod + d.sea.proxy.sl = d.sea.proxy.sl - offset; + + + %%%%%%%%%%%%%%%%%%%%%%% TEMPERATURE %%%%%%%%%%%%%%%%%%%%%%% + + temperature = load(fullfile(folder, set.T_data)); % Load Temp + d.temp.T = temperature.T; + + dyr = diff(d.temp.T(:,1)); + if dyr(1)==mean(dyr) + d.temp.yrs = dyr(1); + else + error('Temperature not recorded in constant steps') + end + d.temp.T(:,3) = d.temp.T(:,3)*set.TerrSc; + + if set.UseMarT0 %Use Marcott RegEM temperature to spin up T0 + temperature = load(fullfile(folder, 'Marcott13_RegEM-HC3_20')); + + d.temp.T0yrs = 20; % year steps + T = temperature.T; + + % bring Marcott temp. (for T0 calculation) to one level whith the + % other temperature's first 'T0temp_level' years of data + ix = T(:,1)>=d.temp.T(1,1) & T(:,1)<=d.temp.T(1,1)+set.T0temp_level; + offs = mean(d.temp.T(d.temp.T(:,1)<=d.temp.T(1,1)+set.T0temp_level,2)); + T(:,2) = T(:,2) - mean(T(ix,2)) + offs; + + % merge Mar and other temperature + ix = T(:,1) SESL.m (Main SESL script to determine semi-empirical parameter distributions) +% -> Calc_SL_from_Param.m (Calculates time series of sea level, temperatures and c from the parameter set found in SESL.m) +% -> Calc_SESL_Prc.m (Calculate percentiles of sea-level, temperature & other important values) +% -> Calc_SESLProjection.m (Calculate sea-level projections under different RCP scenarios) +% -> Calc_SESLConterfact.m (Calculate sea level with different temperature inputs, factual and counterfactual) + + +start = NaN; % if NaN, starting parameter distribution will be set by simulated annealing, + % else a parameterset should be given in the following order [a1, a2/b, c, log(tau1), log(tau2), T0(1)] + +% MAIN SETTINGS + % Model formulation: + % simpel => dSL/dt = a * (T(t)-T0(t)) + % TwoTau => dSL/dt = a1 * (T(t)-T0_1(t)) + a2 * (T(t)-T0_2(t)) + % ConstRate => dSL/dt = c + a * (T(t)-T0(t)) + % CRdecay => dSL/dt = c(t) + a * (T(t)-T0(t)) + % CRovTau => dSL/dt = c/tau + a * (T(t)-T0(t)) + % -> dT0(t)/dt = (T-T0)/tau + % -> dc(t)/dt = -c(t)/tau + model = 'CRdecay'; % 'ConstRate' 'CRdecay' 'CRovTau' 'TwoTau' 'simpel' + % Sea-level data to load has the name 'SL_dat'_'SL_dat2' + SL_dat = 'GLMW-1ts'; % right now only 'GLMW-1ts'. 'GLMW','GLMW-1amp1ts','GLMW-1ts','GLMW-Gr','GLMW-NC' + T_data = 'Mann09_11'; % 'Mann08_eiv','Mann08_cps','Mann09','Mann09_11','Marcott13_RegEM-HC3_20','PAGES2k_13' + +% MH SETTINGS + JumpDist = 0.005; % Jumping distribution size, default = 0.005 + burning = 1000; % number of monte carlo samples used in the burning-in or spin-up period, default = 1000 + NumSkip = 500; % Thinning: only every NumSkip sample will be selected, default = 500 + sample = 1000; % number of monte carlo samples + NormProb = false; % Use normalized probabilities -> should make no difference + +% TEMPERATURE SETTINGS + Tnum = 100; % Number of ar1 temperatures to evaluate for likelihood, default = 100 + T_err = 'ar1ts'; % How should the temperature variance be treated ? + % 'default' -> T + random noise as in Kemp et al. 2011 PNAS + % 'no' -> Don't add uncertainty + % 'ar1' -> T as AR(1) process with sig as 'default' + % 'ar1ts' -> AR(1) Parameter timescale == exp(-abs(t2-t1)/timescale) + tau_ar1 = 10; % if Terr==ar1ts: timescale of ar1ts uncertainty above (default==10) + % if Terr==ar1: AR(1) parameter + ryrs = 10; % If T_err==default: Number of years for which to add same size noise on temp. in calc_sl + TerrSc = 1; % scaling of temperature uncertainty + +% SEA-LEVEL SETTINGS + UseCov = true; % if existing, use covariance matrix of SL data to calc likelihood + CovShrink = 0; % Sea level covariance matrix (C) shrikage parameter f. C gets shrinked to C' while C0 = diag(diag(C)): C' = f*C0 + (1-f)*C + CovTau = 100; % if ~isnan(CovTau): take the elementwise product of the covariance and a tapering function exp(-delta(t)/CovTau) where CovTau is a time scale (default==100) + NoNegCov = true; % remove negative Cov Mat entries and replace by 0 + fac1 = 1; % multiplier to SL uncertainty to adjust for autocorrelation ('Kemp et al. PNAS 2011' use 10) + +% MODEL SETTINGS + optHo = true; % for each MC sample, calculate the optimal offset between simulated and calibration sea level via MLS + OptimT0 = true; % optimize T0(1) (true) or set to mean(T(T0period)) + UseMarT0 = true; % use Marcott temperature to calc. T0 until 'T_data' starts + T0temp_level = 100; % number of years over which to level Marcott and T_data if UseMarT0==true + calibperiod = 500:2010; % CE years to include for model calibration. + period = -1000:2010; % CE years to include for SL calculation. + baseperiod = 1400:1800; % % Period for which simulated sl and data are normed. Obsolete if OptimT0==true. + T0period = -2000:-1800; % period to initialise T0(1) + T0length = 0; % if==0: T0period is used, else T0(1) = mean(T(first 'T0length' years)) + prior + +% PRIOR SETTINGS + StartDistr = [0.5, 0.3, 0, log10(50), log10(100) , 0.5]; % starting distribution for simulated annealing: [a1, a2/b, c, log(tau1), log(tau2), T0(1)] + all_flat_priors = false; % true, false + AdjustPrior = false; % % adjust the prior in check_settings, according to the model used:true, false + a1_prior = {'uniform' 0 2}; + a2_prior = {'uniform' 0 2}; + tau1_prior = {'uniform',log(30),log(3000)}; + tau2_prior = {'uniform',log(30),log(20000)}; + c_prior = {'flat', [],[]}; + c2000_prior = {'uniform', -.2,.2}; + T01st_prior = {'uniform',-.6,.6}; + T02nd_prior = {'normal',0,.2}; + TauLogUniform = true; % true => Tau priors are loguniform, else uniform + +% SETTINGS FOR CALCULATING TIME SERIES FROM POSTERIOR PARAMETERS + storeNum = 5000; % storeNum sl, T0, c & temperature curves are calculated from + % the posterior parameter set & saved to the folder 'Simu\' + % default = 5000, maximum = sample*Tnum + +% SL PROJECTION SETTINGS + rcps = {'RCP3PD' 'RCP45' 'RCP85'}; % which RCPs to take - {'RCP3PD' 'RCP45' 'RCP85'} + dig = []; % number of digits to use in order to save memory, else [] + magicc_realizations = 600; % Number of MAGICC realizations to use for each RCP; max & default == 600 + +% COUNTERFACTUAL/FACTUAL SL SETTINGS + Tcf = {'mean' 'linrate' 'HadCRUT'}; % which Temperature input to use:' + % 'HadCRUT' - HadCRUT4 temperature, + % 'CMIP5_mean' - 'historicalNat_CMIP35_tas' mean from Gareth Jones + % 'CMIP5' - 'historicalNat_CMIP35_tas' single simu (data column given as input) from Gareth Jones + % 'mean' - mean temperature between 500-1800 CE + % 'linrate' - linear rate temperature between 500-1800 CE + col = []; % Data columns of 'historicalNat_CMIP35_tas' to use, if Tcf=='CMIP5' + +% RUN MAIN SCRIPT SESL + fprintf('\n RUN MAIN SCRIPT SESL \n') + SL = SESL(start,... + 'model', model,... + 'SL_dat', SL_dat,... + 'T_data', T_data,... + 'StartDistr', StartDistr,... + 'all_flat_priors', all_flat_priors,... + 'AdjustPrior', AdjustPrior,... + 'a1_prior', a1_prior,... + 'a2_prior', a2_prior,... + 'tau1_prior', tau1_prior,... + 'tau2_prior', tau2_prior,... + 'c_prior', c_prior,... + 'c2000_prior', c2000_prior,... + 'T01st_prior', T01st_prior,... + 'T02nd_prior', T02nd_prior,... + 'TauLogUniform', TauLogUniform,... + 'optHo', optHo,... + 'OptimT0', OptimT0,... + 'UseMarT0', UseMarT0,... + 'T0temp_level', T0temp_level,... + 'calibperiod', calibperiod,... + 'period', period,... + 'baseperiod', baseperiod,... + 'T0period', T0period,... + 'T0length', T0length,... + 'JumpDist', JumpDist,... + 'burning', burning,... + 'NumSkip', NumSkip,... + 'sample', sample,... + 'NormProb', NormProb,... + 'Tnum', Tnum,... + 'T_err', T_err,... + 'tau_ar1', tau_ar1,... + 'ryrs', ryrs,... + 'TerrSc', TerrSc,... + 'UseCov', UseCov,... + 'CovShrink', CovShrink,... + 'CovTau', CovTau,... + 'NoNegCov', NoNegCov,... + 'fac1', fac1); + save(fullfile('Out', [T_data '_' SL_dat]), 'SL'); + +% CALCULATE TIME SERIES FROM PARAMETERS + num = ceil((SL.settings.sample*SL.settings.Tnum)/storeNum); % number of file to be stored + fprintf('\n CALCULATE TIME SERIES FROM PARAMETERS \n'); + fprintf('\t %d time series file(s) will be generated\n',num); + % storeNum sl, T0, c & temperature curves are calculated from the posterior + % parameter set & saved to the folder 'Simu\' + Calc_SL_from_Param(SL,storeNum,'MH','Simu'); + clear SL + +% CALCULATE PERCENTILES OF SL etc. + fprintf('\n CALCULATE PERCENTILES OF SL etc. \n'); + P = Calc_SESL_Prc(T_data, num,'Simu'); + save(fullfile('Out', [T_data, '_', SL_dat, '_prc']),'P'); + +% CALCULATE SL PROJECTIONS + fprintf('\n CALCULATE SL PROJECTIONS \n'); + SLpr = Calc_SESLProjection(P,rcps,dig,magicc_realizations); + save(fullfile('Out', [T_data, '_', SL_dat, '_proj']),'SLpr'); + clear SLpr + +% CALCULATE SL COUNTERFACTUALS + fprintf('\n CALCULATE SL COUNTERFACTUALS \n'); + for i = 1: size(Tcf,2) + if strcmp(cell2mat(Tcf(i)),'CMIP5') + for ii = 1:length(col) + SLcf.([cell2mat(Tcf(i)) '_' num2str(col(ii))]) = Calc_SESLConterfact(P,cell2mat(Tcf(i)),col); + end + else + SLcf.(cell2mat(Tcf(i))) = Calc_SESLConterfact(P,cell2mat(Tcf(i)),col); + end + end + save(fullfile('Out', [T_data, '_', SL_dat, '_cf']),'SLcf'); + clear SLcf +end diff --git a/SESL.m b/SESL.m new file mode 100755 index 0000000..90040b3 --- /dev/null +++ b/SESL.m @@ -0,0 +1,620 @@ +function out = SESL(P0,varargin) + +% Main semi-empirical sea-level (SESL) script to determine semi-empirical +% parameter distributions +% +% out = SESL(P0,varargin) +% +% INPUT: +% - P0 -> startin parameter set [a1, a2, c, tau1, tau2, To(1)]. +% If P0=NaN it will be asigned via simualted annealing below. +% Parameters are scaled automatically to have the same order of magnitude +% - varargin -> settings as defined by DefineSettings_SESL +% +% OUTPUT: out. +% - setting -> all the setting from above plus some internal changes and additions +% - data -> all data loaded, temperature & sea level +% - StartParam -> the parameter set from which the simulated annealing starts +% - lowerBound -> lower bound for simulated annealing +% - upperBound -> upper bound for simulated annealing +% - SimAnn -> the output of the simulated annealing if applied +% - MH -> Montecarlo-Hastings output if applied +% - MH.Params: selected parameters in columns [a1, a2, c, tau1, tau2, To(1)] +% - MH.alpha: the MH mixing +% - Slice -> the output of Slice sampling if applied +% - TimeElapsed -> time needed for calculation in seconds + + calc_lik_only = false; %if true: Calc_SL_SLice_MH2_ only calculates the likelihood + % of the last parameter draw from another run with calc_lik_only==false + + if ~calc_lik_only + tic; + set = DefineSettings_SESL(varargin); + dat = LoadData_SESL(set); + set = checksettings(dat,set); % Check that period, calibperiod, baseperiod + % & T0period are not in conflict with temperature & SL data-set lengths. + % Check conflicts of ryrs & data.temp.yrs and data.temp.yrs & tau_lim(1). + % Make prior pdfs. + + % truncate Mar temp for T0 calculation + if set.UseMarT0 + ix = dat.temp.T0temp(:,1)>=set.T0period(1); + dat.temp.T0temp = dat.temp.T0temp(ix,:); + dat.temp.T0burnin = sum(dat.temp.T0temp(:,1) 0 + Lik = zeros(100,1); + for iLik=1:100 + Lik(iLik) = target_distr([P0(1:3) 10^P0(4) 10^P0(5) P0(6)],set,dat,1); + end + if mean(Lik)==0; error('Likelihood of starting parameter set is zero, please find another one.');end + + % Start the Metropolis-Hastings or Slice Sampling. + if set.sample>0 + %----------------------------- MH sampling -------------------------------- + out.MH = MH_sampler(set,dat,P0); + %--------------------------- Slice Sampling ------------------------------- + % out.Slice = Slice_sampler(set,dat,P0); % It never correctly worked if Terr~='no' + end + out.TimeElapsed = toc; + else + out = target_distr(P0.MH.Params(end,:),P0.settings,P0.data,0); + end +end % end main function + +function [P0, lb, ub, P0offset] = def_Start_Lim(set,dat) + +% Define the starting parameter distribution for simulated annealing, +% dependent on the semi-empirical model used. +% +% P0 = [a1, a2, c, log10(tau1), log10(tau2), T0(0)] + P0 = set.StartDistr; + if strcmp(set.model,'TwoTau') + %lower and upper bounds for parameters + lb = [0, 0, -1, log10(dat.temp.yrs), log10(50), -1]; + ub = [10, 10, 1, log10(10000), log10(3000) , 1]; + elseif strcmp(set.model,'ConstRate') + %lower and upper bounds for parameters + lb = [0, 0, -1, log10(dat.temp.yrs), log10(dat.temp.yrs), -1]; + ub = [10, 10, 1, log10(10000), log10(300) , 1]; + elseif strcmp(set.model,'CRdecay') % from Bob + %lower and upper bounds for parameters + lb = [0, 0, -1, log10(dat.temp.yrs), log10(50), -1]; + ub = [10, 10, 1, log10(5000), log10(20000) , 1]; + elseif strcmp(set.model,'CRovTau') + %lower and upper bounds for parameters + lb = [0, 0, -1, log10(dat.temp.yrs), log10(dat.temp.yrs), -1]; + ub = [10, 10, 1, log10(10000), log10(300) , 1]; + elseif strcmp(set.model,'simpel') + %lower and upper bounds for parameters + lb = [0, 0, -5, log10(dat.temp.yrs), log10(dat.temp.yrs), -1]; + ub = [10, 10, 5, log10(100000), log10(1000) , 1]; + end + P0offset = min(lb)*2; +end + +function SimAnn = SimulatedAnnealing(P0,P0offset,lb,ub,set,dat) + +% Start the matlab simulated annealing function simulannealbnd to find an +% maximum likelihood parameter set to start the semi-empirical model from. + + fprintf('\nStart simulated annealing...') + % -----------Probability function (-Likelihood), bounds & start -------- + set_ = set; + + ProbFunction = @(X) target_distr([X(1:3) 10^X(4) 10^X(5) X(6)],set_,dat,1); + + %--------------------- simulated annealing options ------------------------ + options = []; + options = saoptimset(options,'Display','iter','DisplayInterval',400); + options = saoptimset(options,'InitialTemperature',500,'TemperatureFcn',@temperaturefast,'MaxFunEval',8000,'TolFun',1e-3); + + % Find fval=min(1/Likelihood) and the corresponding parameters x0 with simulated annealing + [x0,fval,exitFlag,output] = simulannealbnd(@(x) ProbFunction(10.^(x)+P0offset),log10(P0-P0offset),log10(lb-P0offset),log10(ub-P0offset),options); + + x0 = 10.^(x0)+P0offset; + + fprintf('The number of iterations was : %d\n', output.iterations); + fprintf('The number of function evaluations was : %d\n', output.funccount); + prob = -fval; + + fprintf('The maximum likelihood found was : %g\n', prob); + + fprintf('Optimal simulannealbnd Parameters: a1 = %g \n',x0(1)) + if strcmp(set.model,'TwoTau') + fprintf(' a2 = %g \n',x0(2));end + if strcmp(set.model,'ConstRate') || strcmp(set.model,'CRdecay') + fprintf(' c = %g \n',x0(3));end + if strcmp(set.model,'CRovTau') + fprintf(' c = %g \n',x0(3)/10^x0(4));end + fprintf(' tau1 = %g \n',10^x0(4)) + if strcmp(set.model,'TwoTau') || strcmp(set.model,'CRdecay') + fprintf(' tau2 = %g \n',10^x0(5));end + if set.OptimT0 + fprintf(' T0(1) = %g \n',x0(6)) + end + + SimAnn.StartParam = P0; + SimAnn.OptParam = x0; + SimAnn.exitFlag = exitFlag; + SimAnn.output = output; + SimAnn.MaxLik = prob; +end + +function MH = MH_sampler(set,dat,x0) + +% Start the slightly adapted Matlab Metropolis-Hastings algorithm mhsample, +% now called mhsample2 and added as a function below. +% Here the final parameter set gets calculated. + + fprintf('Start MH Sampling: \n# Burnin: %g \n# Samples: %g \n# Thinning: %g \n',set.burning,set.sample, set.NumSkip) + + x0offset = -2; % to avoid taking the log of negative numbers + % Probability function which outputs Likelihood + ProbFunction = @(x) target_distr([x(1:3) 10^x(4) 10^x(5) x(6)],set,dat,0); + + pr = set.JumpDist; + proprnd = @(x) x + randn*pr; + + samplerset = struct('pdf', @(xx) ProbFunction(10.^(xx)+x0offset), ... + 'proprnd',proprnd, ... + 'burnin',set.burning, ... + 'thin',set.NumSkip, ... + 'symmetric',1); + + [smpl, accept] = mhsample2(log10(x0-x0offset),set.sample,samplerset); + + smpl = 10.^smpl+x0offset; + smpl(:,4) = 10.^smpl(:,4); + smpl(:,5) = 10.^smpl(:,5); + MH.Params = smpl; + MH.alpha = accept; +end + +function Slice = Slice_sampler(set,dat,x0) + +% A theoretical alternative to the MH-sampler is the slice samplig method +% which did by now not work in this context. + + fprintf('Start Slice Sampling: \n# Burnin: %g \n# Samples: %g \n# Thinning: %g \n',set.burning,set.sample, set.NumSkip) + % Probability function which outputs Likelihood + ProbFunction = @(x) target_distr(x,set,dat,0); + w = .005;% width -> default 10 + + [rnd, neval] = slicesample2(x0,set.sample,'pdf',ProbFunction,'width',w,'burnin',set.burning,'thin',set.NumSkip); + + Slice.Params = rnd; + Slice.neval = neval; +end + +function prob = target_distr(param,set,dat,negLik) + +% Here the target distribution for the above chosen sampling method gets +% evaluated by calculating a semi-empirical sea-level (calc_sl.m) and +% comparing it to the input sea-level data (calc_prob). + + % Add uncertainty to the temperature input + temp = calc_temp(set,dat.temp); + % Calculate the equilibrium temperature + [T01, T02] = calc_T0(set,temp,dat.temp,param); + % Calculate the semi-empirical sea-level and a few parameters + pri = calc_sl(set, dat.temp,temp,T01,T02,param); + if strcmp(set.model,'CRdecay') && sum(pri.yr==2000) == 1 + c2000 = pri.c(pri.yr==2000,:); + else + c2000 = pri.c; % c is fixed + end + T0500_700 =mean(pri.T01(pri.yr>=500 & pri.yr<=700,:)); + T500_700 =mean(pri.T(pri.yr>=500 & pri.yr<=700,:)); + % Compare sl data input to semi-emp. sl -> Gaussian likelihood + prob = calc_prob(pri,param,T0500_700,T500_700,c2000,set,dat); + if negLik % simulated annealing is looking for a minimum + prob = -mean(prob); + else + prob = mean(prob); + end + +end % end calc_prop + +function printit(settings) + +% Print a few of the settings to the command window + + fprintf('\n The MODEL used is: ') + if strcmp(settings.model,'simpel') + fprintf('dSL/dt = a * (T(t)-T0(t)) \n') + elseif strcmp(settings.model,'ConstRate') + fprintf('dSL/dt = c + a * (T(t)-T0(t)) \n') + elseif strcmp(settings.model,'CRdecay') + fprintf('dSL/dt = c(t) + a * (T(t)-T0(t)) \n') + fprintf(' dc(t)/dt = c(t)/tau_ \n') + elseif strcmp(settings.model,'CRovTau') + fprintf('dSL/dt = c/tau + a * (T(t)-T0(t)) \n') + elseif strcmp(settings.model,'TwoTau') + fprintf('dSL/dt = a1 * (T(t)-T0_1(t)) + a2 * (T(t)-T0_2(t)) \n') + end + fprintf(' dT0(t)/dt=(T-T0)/tau \n') + fprintf(' -> Calib. period: %4d - %4d \n',settings.calibperiod(1),settings.calibperiod(end)) + fprintf(' -> Calc. period : %4d - %4d \n',settings.period(1),settings.period(end)) + fprintf(' -> Base period : %4d - %4d \n',settings.baseperiod(1),settings.baseperiod(end)) + fprintf(' -> To(1) period : %4d - %4d \n',settings.T0period(1),settings.T0period(end)) + + if ~settings.OptimT0 + fprintf(' T0(1) is not optimized T0(1) = mean(T(T0period))\n') + end + + fprintf('-> "%s" temperature + "%s" error \n',settings.T_data,settings.T_err); + if settings.TerrSc>1 + fprintf(' -> T unc. scaled by %1.0f !!!\n',settings.TerrSc); + end + if strcmp(settings.T_err,'default') + fprintf(' added every %2d years \n',settings.ryrs); + elseif strcmp(settings.T_err,'ar1ts') + fprintf(' with ar1 timescale = %2d years \n',settings.tau_ar1); + end + + fprintf('-> "%s" SL proxy \n',settings.SL_dat); + + if ~settings.optHo + fprintf('-> H_0 is NOT optimal \n') + end + + if settings.useCov + fprintf('-> COV MAT for MC update\n'); + if settings.CovShrink ~=0 + fprintf(' Cov Mat scaling f: %4.0f (C_ = f*C0 + (1-f)*C) \n',settings.CovShrink); + end + if ~isnan(settings.CovTau) + fprintf(' Cov Mat scaling tau: %4.0f (C_ = C.*exp(-delta(t)/CovTau)) \n',settings.CovTau); + end + if settings.NoNegCov + fprintf(' Negative terms will be REMOVED from the Cov Mat \n') + end + end +end + +function set = checksettings(dat,set) + +% check that period, calibperiod, baseperiod & T0period are not in +% conflict with temperature & SL data. Make Parameter prior PDFs. + + % CHECK LENGTH OF PERIOD + if set.UseMarT0 + Tpath = 'T0temp'; + else + Tpath = 'T'; + end + if set.period(1)round(dat.temp.T(end,1)+ceil((dat.temp.yrs-1)/2)) + period = period(1):round(dat.temp.T(end,1)+ceil((dat.temp.yrs-1)/2)); + else + period = period(1):set.period(end); + end + if ~(isequal(period,set.period)) + fprintf('\n!!!WARNING!!! "period" changed from %4.0f-%4.0f CE to %4.0f-%4.0f CE',set.period(1),set.period(end),period(1),period(end)) + set.period = period; + end + + % CHECK LENGTH OF T0-PERIOD + if set.T0length==0 + if set.T0period(1)round(dat.temp.(Tpath)(end,1)+ceil((dat.temp.yrs-1)/2)) + T0period = T0period(1):round(dat.temp.(Tpath)(end,1)+ceil((dat.temp.yrs-1)/2)); + elseif set.T0period(end)min([round(dat.temp.(Tpath)(end,1)+ceil((dat.temp.yrs-1)/2)),round(dat.sea.proxy.yr(end))]); + calibperiod = calibperiod(1):min([round(dat.temp.(Tpath)(end,1)+ceil((dat.temp.yrs-1)/2)),round(dat.sea.proxy.yr(end))]); + else + calibperiod = calibperiod(1):set.calibperiod(end); + end + if ~(isequal(calibperiod,set.calibperiod)) + fprintf('\n!!!WARNING!!! "calibperiod" changed from %4.0f-%4.0f CE to %4.0f-%4.0f CE',set.calibperiod(1),set.calibperiod(end),calibperiod(1),calibperiod(end)) + set.calibperiod = calibperiod; + end + + %CHECK LENGTH OF BASEPERIOD + if set.baseperiod(1)round(dat.sea.proxy.yr(end)) + baseperiod = baseperiod(1):round(dat.sea.proxy.yr(end)); + else + baseperiod = baseperiod(1):set.baseperiod(end); + end + if ~(isequal(baseperiod,set.baseperiod)) + fprintf('\n!!!WARNING!!! "baseperiod" changed from %4.0f-%4.0f CE to %4.0f-%4.0f CE',set.baseperiod(1),set.baseperiod(end),baseperiod(1),baseperiod(end)) + set.baseperiod = baseperiod; + end + + %CHECK Tnum & T_err + if (strcmp(set.T_err,'no') || strcmp(set.T_err,'default')) && set.Tnum>1 + fprintf('\n!!!WARNING!!! "Tnum" changed from %4.0f to %4.0f',set.Tnum,1) + set.Tnum = 1; + end + + % ryrs should not be smaller than temperature time steps + if ~(isequal(set.ryrs/dat.temp.yrs,round(set.ryrs/dat.temp.yrs))) + if set.ryrs=set.ryrs/2 + ryrs = set.ryrs+mod(set.ryrs,dat.temp.yrs); + else + ryrs = set.ryrs-mod(set.ryrs,dat.temp.yrs); + end + end + fprintf('\n!!!WARNING!!! "ryrs" changed from %4.0f to %4.0f',set.ryrs,ryrs) + set.ryrs = ryrs; + end + + % Depending on the temp. uncertainty setting, tau_ar1 needs to be >= or <= 1. + % One time it is the AR(1) timescale the other time it represents the AR(1) parameter + if strcmp(set.T_err,'ar1ts') && set.tau_ar1<1 + error('For the temperature uncertainty setting ar1ts, tau_ar1 needs to be >=1. Please adjust!') + elseif strcmp(set.T_err,'ar1') && set.tau_ar1>1 + error('For the temperature uncertainty setting ar1, tau_ar1 needs to be <=1. Please adjust!') + end + + % Adjust param prior according to model + if set.AdjustPrior + fprintf('\n!!!WARNING!!! The PRIORS get ADJUSTED to the model used.') + if strcmp(set.model,'TwoTau') + set.a1_prior = {'uniform' 0 10}; + set.a2_prior = {'uniform' 0 10}; + if set.TauLogUniform + set.tau1_prior = {'uniform',log(1000),log(10000)}; % log uniform distribution + set.tau2_prior = {'uniform',log(dat.temp.yrs),log(3000)}; % log uniform distribution + else + set.tau1_prior = {'uniform',301,10000}; % uniform distribution + set.tau2_prior = {'uniform',dat.temp.yrs,300}; % uniform distribution + end + set.c_prior = {'flat', [],[]}; + elseif strcmp(set.model,'ConstRate') + set.a1_prior = {'uniform' 0 10}; + set.a2_prior = {'flat' [] []}; + if set.TauLogUniform + set.tau1_prior = {'uniform',log(dat.temp.yrs),log(10000)}; % log uniform distribution + set.tau2_prior = {'flat',[],[]}; % log uniform distributio + else + set.tau1_prior = {'uniform',dat.temp.yrs,10000}; % uniform distribution + set.tau2_prior = {'flat',[],[]}; + end + set.c_prior = {'uniform',-1,1}; + elseif strcmp(set.model,'CRdecay') + set.a1_prior = {'uniform' 0 2}; + set.a2_prior = {'flat' [] []}; + if set.TauLogUniform + set.tau1_prior = {'uniform',log(30),log(3000)}; % log uniform distribution + set.tau2_prior = {'uniform',log(1000),log(20000)}; % log uniform distribution + else + set.tau1_prior = {'uniform',30,3000}; % uniform distribution + set.tau2_prior = {'uniform',1000,10000}; % uniform distribution + end + elseif strcmp(set.model,'CRovTau') + set.a1_prior = {'uniform' 0 10}; + set.a2_prior = {'flat' [] []}; + if set.TauLogUniform + set.tau1_prior = {'uniform',log(dat.temp.yrs),log(10000)}; % log uniform distribution + set.tau2_prior = {'flat',[],[]}; + else + set.tau1_prior = {'uniform',dat.temp.yrs,10000}; % uniform distribution + set.tau2_prior = {'flat',[],[]}; + end + set.c_prior = {'uniform', -1, 1}; + elseif strcmp(set.model,'simpel') + set.a1_prior = {'uniform' 0 10}; + set.a2_prior = {'flat' [] []}; + if set.TauLogUniform + set.tau1_prior = {'uniform',log(dat.temp.yrs),log(10000)}; % log uniform distribution + set.tau2_prior = {'flat',[],[]}; % log uniform distribution + else + set.tau1_prior = {'uniform',dat.temp.yrs,10000}; % uniform distribution + set.tau2_prior = {'flat',[],[]}; + end + set.c_prior = {'flat',[],[]}; + end + set.T0_prior = {'normal',0,.2}; + end + + % Check that adjustment time scale is not smaller than timesteps + if set.TauLogUniform + if cell2mat(set.tau1_prior(2))0 && mod(i,thin)==0; % burnin and thin + if mod(i,100)==0 + fprintf('%1.0f ',i) + if mod(i,1000)==0 + fprintf('\n') + end + end + smpl(i/thin,:,:) = x0; + end +end +fprintf('\n') +% Accept rate can be used to optimize the choice of scale parameters in +% random walk MH sampler. See for example Roberts, Gelman and Gilks (1997). +accept = accept/(nsamples*thin+burnin); + +% Move the replicates dimension to the end to make samples easier to +% manipulate. +smpl = permute(smpl,[1 3 2]); + +end +function y = mylog(x) +% my log function is to define to avoid the warnings. +y = -Inf(size(x)); +y(x>0) = log(x(x>0)); +end diff --git a/calc_T0.m b/calc_T0.m new file mode 100755 index 0000000..8108d60 --- /dev/null +++ b/calc_T0.m @@ -0,0 +1,101 @@ +function [T01, T02] = calc_T0(set,temp,tempr,param) +% T0(t) is calculated via T0 = TAU * T_ with +% [g 1 0 0 ...] +% [g^2 g 1 0 ...] +% TAU = [g^3 g^2 g 1 ...] +% [ . . . . ...] +% [ . . . . ...] +% g = 1 - 1/tau +% T_ = [T0(1), T(2)/tau, T(3)/tau, ...] +% if n is the size of the steps of our time series, then +% tau = 1 / (1 -g_original^n) with g_original = 1 - 1/tau_original +% +% [T01 T02] = calc_T0(set,temp,tempr,param) +% +% INPUT: +% - set -> general settings +% - temp -> temperature draw, dependent on the temperature error model +% - tempr -> original temperature input +% - param -> the parameter set to test [a1, a2, c, tau1, tau2, To(1)] +% +% OUTPUT: +% - T01 -> equilibrium temperature 1 +% - T02 -> equilibrium temperature 2, only if strcmp(model,'TwoTau')==1 + + if set.OptimT0 + T0_rnd = param(6); + else + T0_rnd = 0; + end + tau1 = param(4); + tau2 = param(5); + + nyrs = size(temp,1); + + % Calc T01 with long adjustment time scale, either with burnin in (UseMarT0 + % -> two different year steps) or without + if set.UseMarT0 + % number of years over which obs is avereaged (year steps of time series) + yrs1 = tempr.T0yrs; + yrs2 = tempr.yrs; + tau1_1 = tau1/yrs1; + tau1_2 = tau1/yrs2; + G1 = 1-1/tau1_1; + G2 = 1-1/tau1_2; + + G1_ = (ones(1,nyrs)*G1).^[0:nyrs-1]; + G1_M = toeplitz(G1_,[G1_(1) zeros(1,size(G1_,2)-1)]); + G2_ = (ones(1,nyrs)*G2).^[0:nyrs-1]; + G2_M = toeplitz(G2_,[G2_(1) zeros(1,size(G2_,2)-1)]); + + G_M1 = [G1_M(:,1:tempr.T0burnin) G2_M(:,tempr.T0burnin+1:end)]; + temp_1 = [temp(1:tempr.T0burnin,:)./tau1_1;temp(tempr.T0burnin+1:end,:)./tau1_2]; + temp_1(1,:) = mean(temp(tempr.T0temp(:,1)<=set.T0period(end),:),1)+ T0_rnd; + else + yrs = tempr.yrs; + tau1 = tau1/yrs; + G = 1-1/tau1; + G_ = (ones(1,nyrs)*G).^[0:nyrs-1]; + G_M1 = toeplitz(G_,[G_(1) zeros(1,size(G_,2)-1)]); + temp_1 = temp/tau1; + temp_1(1,:) = mean(temp(tempr.T(:,1)<=set.T0period(end),:),1)+ T0_rnd; + end + + % Calc T02 with short time scale + if strcmp(set.model,'TwoTau') + if set.UseMarT0 + tau2_1 = tau2/yrs1; + tau2_2 = tau2/yrs2; + G1 = 1-1/tau2_1; + G2 = 1-1/tau2_2; + + G1_ = (ones(1,nyrs)*G1).^[0:nyrs-1]; + G1_M = toeplitz(G1_,[G1_(1) zeros(1,size(G1_,2)-1)]); + G2_ = (ones(1,nyrs)*G2).^[0:nyrs-1]; + G2_M = toeplitz(G2_,[G2_(1) zeros(1,size(G2_,2)-1)]); + + G_M2 = [G1_M(:,1:tempr.T0burnin) G2_M(:,tempr.T0burnin+1:end)]; + temp_2 = [temp(1:tempr.T0burnin,:)./tau2_1;temp(tempr.T0burnin+1:end,:)./tau2_2]; + temp_2(1,:) = mean(temp(tempr.T0temp(:,1)<=set.T0period(end),:),1); + else + tau2 = tau2/yrs; + G = 1-1/tau2; + G_ = (ones(1,nyrs)*G).^[0:nyrs-1]; + G_M2 = toeplitz(G_,[G_(1) zeros(1,size(G_,2)-1)]); + temp_2 = temp/tau2; + temp_2(1,:) = mean(temp(tempr.T(:,1)<=set.T0period(end),:),1); + end + end + + % Calc T0 for each of the Tnum temperature realizations + for i = 1:set.Tnum + T01(:,i) = G_M1 * temp_1(:,i); + if strcmp(set.model,'TwoTau') + T02(:,i) = G_M2 * temp_2(:,i); + end + end + if ~strcmp(set.model,'TwoTau') + T02 = nan(size(T01)); + end + +end diff --git a/calc_prob.m b/calc_prob.m new file mode 100755 index 0000000..146920a --- /dev/null +++ b/calc_prob.m @@ -0,0 +1,110 @@ +function prob = calc_prob(S,param,T0500_700,T500_700,c2000,set,dat) + +% Calculate the Gaussian probability of a certain parameter set. +% +% prob = calc_prob(S,param,T0500_700,T500_700,c2000,set,dat) +% +% INPUT: +% - S -> output of calc_sl.m +% - param -> the parameter set to test [a1, a2, c, tau1, tau2, To(1)] +% - T0500_700 -> mean equilibrium temperature between 500-700 CE +% - T500_700 -> mean temperature between 500-700 CE +% - c2000 -> c in the year 2000 CE +% - set -> general settings +% - dat -> input data, temperature and sea-level +% +% OUTPUT: prob is a gaussian probability as defined in Kopp et al. 2016 PNAS + + % Proxy SL data + P = dat.sea.proxy; + + % calculate optimal Ho between first & last calibration year + firstyear_calib = set.calibperiod(1); + lastyear_calib = set.calibperiod(end); + if set.optHo + Ho_offs = find_optim_Ho(set,P,S); + else + offset = mean(S.sl(S.yr>=set.baseperiod(1) & S.yr<=set.baseperiod(end),:)); % offset to baseperiod + offset = repmat(offset,size(S.sl,1),1); + S.sl = S.sl - offset; + end + + % calc residual 'sq_' of simulated and proxy sea level + S.sl = S.sl(S.yr>=firstyear_calib & S.yr<=lastyear_calib,:); + S.yr = S.yr(S.yr>=firstyear_calib & S.yr<=lastyear_calib); + + ix1 = ismember(round(P.yr),S.yr); + P.yr = round(P.yr(ix1)); + P.sl = P.sl(ix1); + P.sl = repmat(P.sl',1,size(S.sl,2)); + + ix2 = ismember(S.yr,P.yr); + S.sl = S.sl(ix2,:); + if set.optHo + Ho_offs = repmat(Ho_offs,size(S.sl,1),1); + sq_ = (S.sl - P.sl + Ho_offs); + else + sq_ = (S.sl - P.sl); + end + + % truncate SL covariance matrix to the size of residuals + if set.useCov + C = P.C(ix1,ix1); + else + P.slerr = P.slerr(ix1); + C = diag((P.slerr.^2).*set.fac1); + end + + % calc likelihood of priors + if ~set.all_flat_priors + + %param = [a1, a2, c, tau1, tau2, To(1)] + lik(1) = log(set.lik_a1(param(1))); + lik(2) = log(set.lik_a2(param(2))); + lik(3) = log(set.lik_c(param(3))); + if param(4)>0 + if set.TauLogUniform + lik(4) = log(set.lik_tau1(log(param(4)))); % loguniform + else + lik(4) = log(set.lik_tau1(param(4))); % uniform + end + else + lik(4) = log(0); + end + if param(5)>0 + if set.TauLogUniform + lik(5) = log(set.lik_tau2(log(param(5)))); % loguniform + else + lik(5) = log(set.lik_tau2(param(5))); %uniform + end + else + lik(5) = log(0); + end + if set.OptimT0 + lik(6) = log(set.lik_T01st(param(6))); + end + + else + lik = 0; + end + + w = nan(1,size(sq_,2)); + for i = 1:size(sq_,2) + if ~set.all_flat_priors + lik(7) = log(set.lik_T02nd(T0500_700(i),T500_700(i))); + if strcmp(set.model,'CRdecay') + lik(8) = log(set.lik_c2000(c2000(i))); + end + end + sq = sq_(:,i); + w(i) = (-.5 * sq')*(C\sq); + if set.NormProb + w(i) = w(i) - .5*log(2*pi*det(C)); + end + w(i) = w(i) + sum(lik); + w(i) = exp(w(i)); + end + + prob = w; + +end diff --git a/calc_sl.m b/calc_sl.m new file mode 100755 index 0000000..5249e22 --- /dev/null +++ b/calc_sl.m @@ -0,0 +1,64 @@ +function pri2 = calc_sl(set,tempr,temp,T01,T02,param) + +% Calculates a sea-level time-series, dependent on the model used +% +% pri2 = calc_sl(set,tempr,temp,T01,T02,param) +% +% INPUT: +% - set -> general settings +% - tempr -> original temperature input +% - temp -> temperature draw, dependent on the temperature error model +% - T01 -> equilibrium temperature 1 +% - T02 -> equilibrium temperature 2, only if strcmp(model,'TwoTau')==1 +% - param -> the parameter set to test [a1, a2, c, tau1, tau2, To(1)] +% +% OUTPUT: +% - sl -> sea level; +% - dsl -> diff(sl); +% - yr -> year; +% - c -> c; +% - T -> temperature; +% - T01 -> equilibrium temperature 1 +% - T02 -> equilibrium temperature 2, only if strcmp(model,'TwoTau')==1 + + + a1 = param(1); + a2 = param(2); + c = param(3); + tau1 = param(4); + tau2 = param(5); + + [year, temp, T01, T02] = resize_T(set,tempr,temp,T01,T02); % make a step-like yearly time series + + % Calc dSL/dt = c + a1*(temp - T01) + a2*(temp - T02); + if strcmp(set.model,'TwoTau') + dsea = a1*(temp - T01) + a2*(temp - T02); + pri2.T02 = T02; + elseif strcmp(set.model,'ConstRate') + dsea = c + a1*(temp - T01); + elseif strcmp(set.model,'CRdecay') + % calc. c decaying exponentially on timescale tau2 + g = 1 - 1/tau2; + G = (ones(size(temp,1),1)*g).^([0:size(temp,1)-1]'); + c = c*G; + c = repmat(c,1,set.Tnum); + dsea = c + a1*(temp - T01); + elseif strcmp(set.model,'CRovTau') + dsea = c/tau1 + a1*(temp - T01); + elseif strcmp(set.model,'simpel') + dsea = a1*(temp - T01); + end + + sea = cumsum(dsea); + +if real(sea)~=sea + disp('sea not real') +end + pri2.sl = sea; + pri2.dsl = dsea; + pri2.yr = year; + pri2.c = c; + pri2.T = temp; + pri2.T01 = T01; + +end % end calc_sl diff --git a/calc_temp.m b/calc_temp.m new file mode 100755 index 0000000..da8dce7 --- /dev/null +++ b/calc_temp.m @@ -0,0 +1,63 @@ +function temp = calc_temp(set,tempr) + +% Represent the temperature uncertainty +% +% temp = calc_temp(set,tempr) +% +% INPUT: +% - set -> general settings +% - tempr -> original temperature input +% +% OUTPUT: +% - temp -> temperature draws with uncertainty depending on the temperature error model + + if strcmp(set.T_err,'ar1') + ar1 = set.tau_ar1; + elseif strcmp(set.T_err,'ar1ts') + ar1 = set.tau_ar1; + else + ar1 = []; + end + + if strcmp(set.T_err,'no') + if set.UseMarT0 %additionally to T_data, Marcott is used to calc. T0 + temp = tempr.T0temp(:,2)'; + else + temp = tempr.T(:,2)'; + end + elseif strcmp(set.T_err,'default') + if set.UseMarT0 %additionally to T_data, Marcott is used to calc. T0 + T = tempr.T0temp; + else + T = tempr.T; + end + % add noise, sampled from a normal distribution with sig=Terr, to T every ryrs + ryrs = set.ryrs/tempr.yrs; % Number of years for which to add same noise + nyrs = length(T); + Terr = T(:,3); + r = randn(round(nyrs/ryrs + 1), 1); + + mat = r(1:max(size(Terr(1:ryrs:end))),1) * ones(1,ryrs); + tempE = reshape(mat',numel(mat),1); + + tempE = Terr.* tempE(1:nyrs); + + temp = (T(:,2) + tempE)'; + + elseif strcmp(set.T_err(1:2),'ar') + if set.UseMarT0 %additionally to T_data, Marcott is used to calc. T0 + T = tempr.T0temp; + else + T = tempr.T; + end + if strcmp(set.T_err,'ar1'); + Cov_ar1 = bsxfun(@times,T(:,3),T(:,3)').* ((ar1*ones(length(T))).^abs(bsxfun(@minus,T(:,1),T(:,1)'))); + elseif strcmp(set.T_err,'ar1ts'); + Cov_ar1 = bsxfun(@times,T(:,3),T(:,3)').* exp(-abs(bsxfun(@minus,T(:,1),T(:,1)')/ar1)); + end + temp = mvnrnd(T(:,2)',Cov_ar1,set.Tnum); + end + + temp = temp'; + +end diff --git a/find_optim_Ho.m b/find_optim_Ho.m new file mode 100755 index 0000000..c6ff126 --- /dev/null +++ b/find_optim_Ho.m @@ -0,0 +1,49 @@ +function Ho_offs = find_optim_Ho(set,P,S) + +% Find optimal offset between simulated and proxy sea level in a least square sense. +% +% Ho_offs = find_optim_Ho(set,P,S) +% +% INPUT: +% - set -> general settings +% - P -> sea-level data input +% - S -> output of calc_sl, i.e. semi-empirical sea-level simulation +% +% OUTPUT: Ho_offs is the optimal (in a least square sense) offset of semi-empirical +% and input sea-level data to be subtracted from the latter one. + + firstyear_calib = set.calibperiod(1); + lastyear_calib = set.calibperiod(end); + + % reduce Proxy timeseries to calibration time frame + ix = P.yr>=firstyear_calib & P.yr<=lastyear_calib; + Psea = P.sl(ix); + Pyear = P.yr(ix); + Pslerr = P.slerr(ix).^2; + + % bring SL_simu to size of proxy (same years) + ix_ = ismember(S.yr,round(Pyear)); + sea = S.sl(ix_,:); + year = S.yr(ix_); + + % double check size of proxy and andjust to SL_simu if necessary + ix__ = ismember(round(Pyear),year); + Psea = Psea(ix__); + Pslerr = Pslerr(ix__); + + % Define error structure + if exist('C','var') == 1 && set.useCov + Pcov = P.C(ix,ix); %Pcov = diag(diag(C(ix,ix))); + Pcov = Pcov(ix__,ix__); + else + Pcov = diag(Pslerr); + end + if ~(size(Pslerr,2)==size(Psea,2)) || ~(size(Pslerr,2)==size(sea,1)) || ~(size(Psea,2)==size(sea,1)) || ~(size(Pslerr,2)==size(Pcov,2)) || ~(size(sea,1)==size(Pcov,2)) || ~(size(Psea,2)==size(Pcov,2)) + error('check size') + end + + % Calc opt. offset + Psea = repmat(Psea',1,size(sea,2)); + Ho_offs = lscov(ones(size(Pslerr')),Psea-sea,Pcov); + +end diff --git a/readme.md b/readme.md new file mode 100644 index 0000000..8761177 --- /dev/null +++ b/readme.md @@ -0,0 +1,66 @@ +# SESL: Semi-Empirical Sea Level + +README file last updated by Robert Kopp, robert-dot-kopp-at-rutgers-dot-edu, Mon Dec 05 10:26:00 EST 2016 + +## Citation and Acknowledgements + +Version 1.0 of this code was released on 5 December 2016 to accompany + + Kopp, R. E., A. C. Kemp, K. Bittermann, B. P. Horton, + J. P. Donnelly, W. R. Gehrels, C. C. Hay, J. X. Mitrovica, + E. D. Morrow, and S. Rahmstorf (2016). Temperature-driven + global sea-level variability in the Common Era. Proceedings + of the National Academy of Sciences. doi: 10.1073/pnas.1517056113. + +Please cite this source when using this code. + +This code with developed by Klaus Bittermann, with assistance from Eric Morrow. + +## Folders + +- Data: holds all input data +- Simu: holds the output of 'Calc_SL_from_Param.m' +- Out: holds all other output + +## Main Code + +- Run_SESL.m: Runs everything and gives the opportunity to define all settings + +- SESL.m: Main script with MH-algorithm which outputs the posterior parameter set +- DefineSettings_SESL.m: Defines all settings for SESL.m +- LoadData_SESL.m: Loads the input data for SESL.m + +- Calc_SL_from_Param.m: Calculates sl, T0, c & temperature time series from the posterior parameter set (SESL.m output) and saves them. The number of files to split the output to can be set in order to save memory. + +- Calc_SESL_Prc.m: Calculates percentiles of sea level and other parameters necessary for projections. To do this it loads the files stored by Calc_SL_from_Param.m above. + +- Calc_SESLProjection.m: Calculates projections with RCP temperature inputs from the output of Calc_SESL_Prc.m + +- Calc_SESLConterfact.m: Calculates sea level for counterfactual temperature inputs from the output of Calc_SESL_Prc.m + +Other external functions called by SESL.m & Calc_SL_from_Param.m: + +- calc_sl +- calc_temp +- calc_T0 +- calc_prob +- find_optim_Ho +- resize_T + + +---- + + Copyright (C) 2016 by Klaus Bittermann + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . diff --git a/resize_T.m b/resize_T.m new file mode 100755 index 0000000..bc0d61e --- /dev/null +++ b/resize_T.m @@ -0,0 +1,88 @@ +function [year, temp, T01, T02] = resize_T(set,tempr,temp, T01, T02) + +% If the input temperature is not given in yearly values, here it will, +% alongside with the equilibrium temperatures, be resized to a step-like +% yearly time series +% +% [year temp T01 T02] = resize_T(set,tempr,temp, T01, T02) +% +% INPUT: +% - set -> general settings +% - tempr -> original temperature input +% - temp -> temperature draw, dependent on the temperature error model +% - T01 -> equilibrium temperature 1 +% - T02 -> equilibrium temperature 2, only if strcmp(model,'TwoTau')==1 +% +% OUTPUT: +% - year -> year; +% - temp -> step-like yearly temperature; +% - T01 -> step-like yearly equilibrium temperature 1 +% - T02 -> step-like yearly equilibrium temperature 2, only if strcmp(model,'TwoTau')==1 + + if set.UseMarT0 + Tpath = 'T0temp'; + else + Tpath = 'T'; + end + + % Turn temp & T0 to steplike yearly data and truncate to desired length + fyr = max([tempr.T(1,1) min([set.period(1) set.calibperiod(1)])]); + lyr = max([set.period(end) set.calibperiod(end)]); + + year = repmat(tempr.T(tempr.T(:,1)>=fyr & tempr.T(:,1)<=lyr,1),tempr.yrs,1); + [year, ix] = sort(year); + year1 = [year(1)-floor((tempr.yrs-1)/2) : year(end)+ceil((tempr.yrs-1)/2)];%tempr.yr; + + temp1 = temp(tempr.(Tpath)(:,1)>=max([tempr.T(1,1) fyr]),:); + temp1 = repmat(temp1,tempr.yrs,1); + temp1 = temp1(ix,:); + + T01_1 = T01(tempr.(Tpath)(:,1)>=max([tempr.T(1,1) fyr]),:); + T01_1 = repmat(T01_1,tempr.yrs,1); + T01_1 = T01_1(ix,:); + if strcmp(set.model,'TwoTau') + T02_1 = T02(tempr.(Tpath)(:,1)>=max([tempr.T(1,1) fyr]),:); + T02_1 = repmat(T02_1,tempr.yrs,1); + T02_1 = T02_1(ix,:); + end + + % Turn temp & T0 to steplike yearly data and truncate to desired length + fyr_ = max([tempr.(Tpath)(1,1),min([set.period(1) set.calibperiod(1)])]); + lyr_ = fyr-1; + + if fyr~=fyr_ && set.UseMarT0 % if Marcott temperature is used to spin up T0 + year = repmat(tempr.T0temp(tempr.T0temp(:,1)>=fyr_ & tempr.T0temp(:,1)<=lyr_,1),tempr.T0yrs,1); + [year, ix] = sort(year); + year2 = [year(1)-floor((tempr.T0yrs-1)/2) : year(end)+ceil((tempr.T0yrs-1)/2)];%tempr.yr; + + temp2 = temp(tempr.(Tpath)(:,1)>=fyr_& tempr.(Tpath)(:,1)<=lyr_,:); + temp2 = repmat(temp2,tempr.T0yrs,1); + temp2 = temp2(ix,:); + + T01_2 = T01(tempr.(Tpath)(:,1)>=fyr_& tempr.(Tpath)(:,1)<=lyr_,:); + T01_2 = repmat(T01_2,tempr.T0yrs,1); + T01_2 = T01_2(ix,:); + if strcmp(set.model,'TwoTau') + T02_2 = T02(tempr.(Tpath)(:,1)>=fyr_& tempr.(Tpath)(:,1)<=lyr_,:); + T02_2 = repmat(T02_2,tempr.T0yrs,1); + T02_2 = T02_2(ix,:); + end + + temp = [temp2(year2