%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB Script to reproduce the sensitivity of future Antarctic surface mass
% balance over the grounded ice sheet and the ice shelves to Antarctic (90–60°S) 
% near-surface warming from Coulon et al. (2023): Yearly Bayesian calibrated 
% median surface mass balance, snow accumulation and surface runoff (Gt yr−1)
% projected over the period 2015-2300 are compared with the mean annual near surface 
% temperature anomaly (in °C, with respect to the 1995-2014 average) projected by 
% four GCMs from the sixth phase of the Coupled Model Intercomparison Project (CMIP6) 
% over the Antarctic domain under the Shared Socioeconomic Pathways 1-2.6 and 5-8.5
%
% Author: Violaine Coulon (violaine.coulon@ulb.be)
% Date: October 2023
%
% Data Source:
% The data used in this script are available at 10.5281/zenodo.8398772
%
% MATLAB Version and Toolboxes:
% This script is designed to work with MATLAB R20XX or later and does not
% require any additional toolboxes beyond the standard MATLAB functions.
% Some functions used in the script might be available in earlier versions
% of MATLAB as well.
%
% Please make sure you have the necessary MATLAB version installed to run
% this script.
%
% Usage:
% 1. Ensure you have the necessary MATLAB version and toolboxes installed.
% 2. Download the paper and the required data from the provided sources.
% 3. Run the script 'Compute_Bayesian_Weight.m' to generate 'ENSEMBLE_SCORE.mat'
% 4. Run this script to generate the figure.
%
% Citation:
% If you use this script or the reproduced figure in your work, please
% consider citing the original paper as follows:
%
% Coulon, V., Klose, A. K., Kittel, C., Edwards, T., Turner, F., Winkelmann, R., and Pattyn, F.: 
% Disentangling the drivers of future Antarctic ice loss with a historically-calibrated ice-sheet model, 
% EGUsphere [preprint], https://doi.org/10.5194/egusphere-2023-1532, 2023.
%
% Disclaimer:
% This script is provided as-is without any warranties. The author is not
% responsible for any inaccuracies or errors in the reproduced figure.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fig=figure('Position', [100, 300, 1200, 500]);
subplot(1,2,1)
plot_SMB_med_vs_aTs('GR')
title('GROUNDED ICE')
subplot(1,2,2)
plot_SMB_med_vs_aTs('FL')
title('ICE SHELVES')
legend off

function plot_SMB_med_vs_aTs(type)

load ENSEMBLE_SCORE

scenario='SSP585';
polynomial_deg=3;

% PLOT ZERO LINES
x=0:0.1:20;
y=polyval(x.*0,x);
plot(x,y,'color',[0.5 0.5 0.5],'LineWidth',1,'LineStyle','--','HandleVisibility','off'); hold on  

CMIP6model={'MRI-ESM2-0', 'IPSL-CM6A-LR', 'CESM2-WACCM', 'UKESM1-0-LL'}; 
aTs_mean_TOT=zeros(length(CMIP6model),286);
for i=1:length(CMIP6model)
    load(['',char(CMIP6model(i)),'_',scenario,'_mean_aTs'])
    aTs_mean_TOT(i,:)=aTs_mean;
end
aTs_mean=mean(aTs_mean_TOT,1);

load(['',scenario,'_ENSEMBLE_DATA.mat'],'mbcomp_ensemble','time')
mbcomp_ensemble=mbcomp_ensemble(1:286,:,:);

% SMB 
cl=[8/255 105/255 142/255];
if isequal(type,'AIS')
    var=(squeeze(mbcomp_ensemble(:,1,:)))';
elseif isequal(type,'GR')
    var=(squeeze(mbcomp_ensemble(:,8,:)))';
elseif isequal(type,'FL')
    var=(squeeze(mbcomp_ensemble(:,1,:)-mbcomp_ensemble(:,8,:)))';
end
%meanvar=mean(var/1000.*w',1);
q=[0.5];
quart=zeros(length(q),size(var,2));
for ii=1:size(var,2)
    quart(:,ii)=quantile_eval(var(:,ii)./1000,w',q);
end
meanvar=quart;
scatter(movmean(aTs_mean,5),meanvar,8,'o','MarkerEdgeColor',cl,'MarkerFaceColor',cl,'HandleVisibility','off','MarkerEdgeAlpha',0.4); hold on
SMB=meanvar;
% ACC 
cl=[124 196 204]/255;
if isequal(type,'AIS')
    var=(squeeze(mbcomp_ensemble(:,2,:)))';
elseif isequal(type,'GR')
    var=(squeeze(mbcomp_ensemble(:,9,:)))';
elseif isequal(type,'FL')
    var=(squeeze(mbcomp_ensemble(:,2,:)-mbcomp_ensemble(:,9,:)))';
end
%meanvar=mean(var/1000.*w',1);
q=[0.5];
quart=zeros(length(q),size(var,2));
for ii=1:size(var,2)
    quart(:,ii)=quantile_eval(var(:,ii)./1000,w',q);
end
meanvar=quart;
scatter(movmean(aTs_mean,5),meanvar,8,'o','MarkerEdgeColor',cl,'MarkerFaceColor',cl,'HandleVisibility','off','MarkerEdgeAlpha',0.4); hold on
ACC=meanvar;
% RUNOFF 
cl=[237 177 32]/255;
if isequal(type,'AIS')
    var=(squeeze(mbcomp_ensemble(:,4,:)))';
elseif isequal(type,'GR')
    var=(squeeze(mbcomp_ensemble(:,11,:)))';
elseif isequal(type,'FL')
    var=(squeeze(mbcomp_ensemble(:,4,:)-mbcomp_ensemble(:,11,:)))';
end
%meanvar=mean(var/1000.*w',1);
q=[0.5];
quart=zeros(length(q),size(var,2));
for ii=1:size(var,2)
    quart(:,ii)=quantile_eval(var(:,ii)./1000,w',q);
end
meanvar=quart;
scatter(movmean(aTs_mean,5),meanvar,8,'o','MarkerEdgeColor',cl,'MarkerFaceColor',cl,'HandleVisibility','off','MarkerEdgeAlpha',0.4); hold on
RUN=meanvar;

aTS=movmean(aTs_mean,5);

scenario='SSP126'; 

CMIP6model={'MRI-ESM2-0', 'IPSL-CM6A-LR', 'CESM2-WACCM', 'UKESM1-0-LL'}; 
aTs_mean_TOT=zeros(length(CMIP6model),286);
for i=1:length(CMIP6model)
    load(['',char(CMIP6model(i)),'_',scenario,'_mean_aTs'])
    aTs_mean_TOT(i,:)=aTs_mean;
end
aTs_mean=mean(aTs_mean_TOT,1);

load(['',scenario,'_ENSEMBLE_DATA.mat'],'mbcomp_ensemble','time')
mbcomp_ensemble=mbcomp_ensemble(1:286,:,:);

% SMB 
cl=[8/255 105/255 142/255];
if isequal(type,'AIS')
    var=(squeeze(mbcomp_ensemble(:,1,:)))';
elseif isequal(type,'GR')
    var=(squeeze(mbcomp_ensemble(:,8,:)))';
elseif isequal(type,'FL')
    var=(squeeze(mbcomp_ensemble(:,1,:)-mbcomp_ensemble(:,8,:)))';
end
q=[0.5];
quart=zeros(length(q),size(var,2));
for ii=1:size(var,2)
    quart(:,ii)=quantile_eval(var(:,ii)./1000,w',q);
end
meanvar=quart;
scatter(movmean(aTs_mean,5),meanvar,8,'^','MarkerEdgeColor',cl,'MarkerFaceColor','none','HandleVisibility','off','MarkerEdgeAlpha',0.4); hold on
SMB=[SMB meanvar];
% ACC 
cl=[124 196 204]/255;
if isequal(type,'AIS')
    var=(squeeze(mbcomp_ensemble(:,2,:)))';
elseif isequal(type,'GR')
    var=(squeeze(mbcomp_ensemble(:,9,:)))';
elseif isequal(type,'FL')
    var=(squeeze(mbcomp_ensemble(:,2,:)-mbcomp_ensemble(:,9,:)))';
end
q=[0.5];
quart=zeros(length(q),size(var,2));
for ii=1:size(var,2)
    quart(:,ii)=quantile_eval(var(:,ii)./1000,w',q);
end
meanvar=quart;
scatter(movmean(aTs_mean,5),meanvar,8,'^','MarkerEdgeColor',cl,'MarkerFaceColor','none','HandleVisibility','off','MarkerEdgeAlpha',0.4); hold on
ACC=[ACC meanvar];
% RUNOFF 
cl=[237 177 32]/255;
if isequal(type,'AIS')
    var=(squeeze(mbcomp_ensemble(:,4,:)))';
elseif isequal(type,'GR')
    var=(squeeze(mbcomp_ensemble(:,11,:)))';
elseif isequal(type,'FL')
    var=(squeeze(mbcomp_ensemble(:,4,:)-mbcomp_ensemble(:,11,:)))';
end
q=[0.5];
quart=zeros(length(q),size(var,2));
for ii=1:size(var,2)
    quart(:,ii)=quantile_eval(var(:,ii)./1000,w',q);
end
meanvar=quart;
scatter(movmean(aTs_mean,5),meanvar,8,'^','MarkerEdgeColor',cl,'MarkerFaceColor','none','HandleVisibility','off','MarkerEdgeAlpha',0.4); hold on
RUN=[RUN meanvar];

aTS=[aTS movmean(aTs_mean,5)];

% PLOT FITS:

% SMB
cl=[8/255 105/255 142/255];
p1=polyfit(aTS,SMB,polynomial_deg);
x=linspace(min(aTS),15,100);
y=polyval(p1,x);
plot(x,y,'color',cl,'LineWidth',1)
% ACC
cl=[124 196 204]/255;
p1=polyfit(aTS,ACC,polynomial_deg);
x=linspace(min(aTS),15,100);
y=polyval(p1,x);
plot(x,y,'color',cl,'LineWidth',1)
% RUNOFF
cl=[237 177 32]/255;
p1=polyfit(aTS,RUN,polynomial_deg);
x=linspace(min(aTS),15,100);
y=polyval(p1,x);
plot(x,y,'color',cl,'LineWidth',1)
ylabel('Mass Fluxes (10^3 Gt/yr)')
xlabel('Antarctic Near-Surface Warming (°C)')
xlim([0 15])
xticks(0:2.5:15)
grid on
set(gca,'FontSize',11)

legend({'Surface Mass Balance','Snow Accumulation','Surface Runoff'},'Location','northwest','EdgeColor','none','FontSize',12,'color','none')
end

% Evaluation of the alpha-quantile of the data set Y
% weight is a vector a weight assigned to the data
% alpha can be a vector

function quant = quantile_eval(Y,weight,alpha)

N = length(Y);

[Y_sorted,I] = sort(Y); % sort data in ascending order
weight = weight(I);
occ = countmember(Y_sorted,Y_sorted);
[Y_unique,ia] = unique(Y_sorted,'last');
pmasstilde = occ(ia)/N;
weighttilde = weight(ia);
Ftilde = cumsum(pmasstilde.*weighttilde);
Ftilde(end) = 1;
if(alpha>=Ftilde(end))
    quant = Y_unique(end);
    return
elseif(alpha<Ftilde(1))  
    quant = -inf;
    return
end 

quant = zeros(length(alpha),1);
% for ii=1:length(alpha)
%     ind = find(Ftilde>=alpha(ii),1,'first');
%     quant(ii) = Y_unique(ind);
% end
for ii=1:length(alpha)
  if(alpha(ii)>=Ftilde(end))
    quant(ii) = Y_unique(end);
  elseif(alpha(ii)<Ftilde(1))
    quant(ii) = -inf;
  else
    ind = find(Ftilde>=alpha(ii),1,'first');
    quant(ii) = Y_unique(ind);
  end
end

end

function COUNT = countmember(A,B)
% COUNTMEMBER - count members
%
%   COUNT = COUNTMEMBER(A,B) counts the number of times the elements of array A are
%   present in array B, so that C(k) equals the number of occurences of
%   A(k) in B. A may contain non-unique elements. C will have the same size as A.
%   A and B should be of the same type, and can be cell array of strings.
%

narginchk(2,2) ;
if ~isequal(class(A), class(B))
    error('Both inputs should be of the same class.') ;
end
if isempty(A) || isempty(B)
    % nothing to do
    COUNT = zeros(size(A)) ;
else    
    % which elements are unique in A,
    % also store the position to re-order later on
    [AUnique, ~, posJ] = unique(A(:)) ;
    % assign each element in B a number corresponding to the element of A
    [~, loc] = ismember(B, AUnique) ;
    % count these numbers
    N = histcounts(loc(:), 1:length(AUnique)+1) ;
    % re-order according to A, and reshape
    COUNT = reshape(N(posJ),size(A)) ;
end

end