%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB Script to calculate and save the 100-member ensemble likelihood 
% weights used in the Bayesian calibration in Coulon et al. (2023)
%
% 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 this script to generate 'ENSEMBLE_SCORE.mat'
%
% 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.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

fac=10; % multiplier of obs error

% Run for the different configurations (regular, without elevation
% feedback, and with hydrofrcaturing) 

Elev_list=[0 1 0];
HydroF_list=[0 0 1];

for t=1:length(Elev_list)

ElevCorr0=Elev_list(t);
HydroFrac=HydroF_list(t);

if ElevCorr0==1
    load ('HIST_ENSEMBLE_DATA_NO_ELEVATION_FEEDBACK.mat','mbcomp_ensemble','time')
    load ('HIST_ENSEMBLE_DATA_BASIN_NO_ELEVATION_FEEDBACK.mat','mbcomp_ensemble_basin')
elseif HydroFrac==1
    load ('HIST_ENSEMBLE_DATA_HYDROFRAC.mat','mbcomp_ensemble','time')
    load ('HIST_ENSEMBLE_DATA_BASIN_HYDROFRAC.mat','mbcomp_ensemble_basin')
else
    load ('HIST_ENSEMBLE_DATA.mat','mbcomp_ensemble','time')
    load ('HIST_ENSEMBLE_DATA_BASIN.mat','mbcomp_ensemble_basin')
end
load ('CALIBRATION_DATA.mat')

n=length(val); % number of obs data
obs=val(1:n);
sig=sigma(1:n);
mod=zeros(n,size(mbcomp_ensemble,3));
for ii=1:n % loop over obs data
    for tt=1:size(mbcomp_ensemble,3) % loop over ensemble members
        if ii<5
            % calculate mean net mass balance (mbcomp_ensemble_basin(:,:,13,:) of ensemble member tt 
            %  for observational period of obs data ii over WAIS
            mod(ii,tt)=mean(sum(mbcomp_ensemble_basin([1 18:23],find(time==year1(ii)):find(time==year2(ii)),13,tt),1)); 
        elseif ii<9
            %  calculate mean net mass balance (mbcomp_ensemble_basin(:,:,13,:) of ensemble member tt 
            %  for observational period of obs data ii over EAIS
            mod(ii,tt)=mean(sum(mbcomp_ensemble_basin([2:17],find(time==year1(ii)):find(time==year2(ii)),13,tt),1)); % EAIS
        else
            %  calculate mean net mass balance (mbcomp_ensemble_basin(:,:,13,:) of ensemble member tt 
            %  for observational period of obs data ii over APIS
            mod(ii,tt)=mean(sum(mbcomp_ensemble_basin([24:27],find(time==year1(ii)):find(time==year2(ii)),13,tt),1)); % APIS
        end
    end
end

sig_obs=sig; % observational error
sig_mod=fac*sig_obs; % define structural error as fac*observational error
var=(sig_obs).^2+(sig_mod).^2; % discrepancy variance

%% SCORE CALCULATION

s=zeros(1,size(mbcomp_ensemble,3));
% calculate likelihood score for each ensemble member
for i=1:size(mbcomp_ensemble,3)
    s(i)=exp(-0.5.*sum(((mod(:,i)-obs).^2)./(var)));
end

% to create a weight w, the score for each ensemble member is then
% normalised
w=s./sum(s);

w=w./mean(w); % ensures that sum(w)=N

if ElevCorr0==1
    save('ENSEMBLE_SCORE_NO_ELEVATION_FEEDBACK','w')
elseif HydroFrac==1
    save('ENSEMBLE_SCORE_HYDROFRAC','w')
else
    save('ENSEMBLE_SCORE','w')
end

end


