steps=[1];

model_list={'AWI_PISM1','DOE_MALI','ILTS_PIK_SICOPOLIS','IMAU_IMAUICE1','IMAU_IMAUICE2','JPL1_ISSM','LSCE_GRISLI','NCAR_CISM','PIK_PISM1','PIK_PISM2','UCIJPL_ISSM','ULB_FETISH32','ULB_FETISH16','UTAS_ElmerIce','VUB_AISMPALEO','VUW_PISM'};

model_list2=model_list;
for i=1:numel(model_list2)
	   model_list2{i} = strrep(model_list2{i},'_','\_');
end
model_list3=model_list; %backwards for pcolor
for i=1:numel(model_list2)
	model_list3{i} = model_list2{end-i+1};
end

org=organizer('repository','./Models','prefix','Projections.','steps',steps);

%Figure 1
if perform(org,'Plot_NonLinearity_Volume')% {{{

	model_list={'ILTS_PIK_SICOPOLIS','JPL1_ISSM','LSCE_GRISLI','NCAR_CISM','UCIJPL_ISSM'};
	colors = distinguishable_colors(20);
	experiments_list={'expC2','expC3','expC5','expC6','expC11','expC12',... %open and standard NorESM %only RCP85
		};
	%experiments_list={'expC2','expC3','expC5','expC6','expC8','expC9','expC11','expC12',... %open and standard NorESM

	results=[]; results_baseline=[];
	resultstendacabf=[]; resultstendacabf_baseline=[];

	results_model={};
	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					explim_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivol_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabf_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smb_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explim_file,'time');
					lim=ncread(explim_file,'ivol')*ice_density/(10^9*1000); %from m^3 to Gt;
					tendacabf=ncread(exptendacabf_file,'smb')*yearday_model*3600*24/(10^9*1000); %from kg/s to Gt/yr;
					results(end+1)=lim(85);
					resultstendacabf(end+1)=sum(tendacabf(1:85));
					if length(time_model)>90, error('too long '); end
					if strcmpi(expename,'expC2'),
						expename='exp01';
						results_model{end+1}=[model_list2{imodel} '\_NorESM'] ;
					elseif strcmpi(expename,'expC3'),
						expename='exp05';
						results_model{end+1}=[model_list2{imodel} '\_NorESM'] ;
					elseif strcmpi(expename,'expC5'),
						expename='exp02';
						results_model{end+1}=[model_list2{imodel} '\_MIROC'] ;
					elseif strcmpi(expename,'expC6'),
						expename='exp06';
						results_model{end+1}=[model_list2{imodel} '\_MIROC'] ;
					elseif strcmpi(expename,'expC8'),
						expename='exp03';
						results_model{end+1}=[model_list2{imodel} '\_NorESM\_2.6'] ;
					elseif strcmpi(expename,'expC9'),
						expename='exp07';
						results_model{end+1}=[model_list2{imodel} '\_NorESM\_2.6'] ;
					elseif strcmpi(expename,'expC11'),
						expename='exp04';
						results_model{end+1}=[model_list2{imodel} '\_CCSM4'] ;
					elseif strcmpi(expename,'expC12'),
						expename='exp08';
						results_model{end+1}=[model_list2{imodel} '\_CCSM'] ;
					end
					explim_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivol_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabf_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smb_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explim_file,'time');
					lim=ncread(explim_file,'ivol')*ice_density/(10^9*1000); %from m^3 to Gt;
					tendacabf=ncread(exptendacabf_file,'smb')*yearday_model*3600*24/(10^9*1000); %from kg/s to Gt/yr;
					results_baseline(end+1)=lim(85);
					resultstendacabf_baseline(end+1)=sum(tendacabf(1:85));
				end

			end %end of model
		end %end of isexp
		allresults=[results_baseline; results; results_baseline-resultstendacabf_baseline]'; 

		%Plot results
		figure(2); clf; set(gcf,'color','w'); set(gcf,'Position',[400 500 900 400]);
		hb=bar(allresults*(-1/362.5),'grouped'); hold on 
		set(gca, 'fontsize',14);
		legend('All forcings','Ocean-only ','Ocean-calculated','location','NorthEast','fontsize',14)
		%set(gca, 'XTick', 1:size(allresults,1) );
		%set(gca, 'XTicklabels', results_model ,'fontsize',8);
		%xtickangle(90)
		set(gca, 'XTickLabel', [])
		ylabel('Ice mass loss (mm SLE)','fontsize',16);
		ylim([-120 500])
		set(gca, 'YGrid', 'on', 'XGrid', 'off')
		text(-1.0,-100,'a','fontsize',18,'fontweight','b');

		export_fig('Figures/NonLinearityVolume.pdf')
		error
		h=text(1,-48,'WAIS','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		h=text(2,-48,'EAIS','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		h=text(3,-48,'Peninsula','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		for iresults=1:length(results_model),
			bi=bar(iresults,results(iresults)*(-1/361.8),'facecolor',colors(iresults,:));
			h=text(iresults,-42,results{iresults},'VerticalAlignment','middle','HorizontalAlignment','right','color',colors(iresults,:),'fontweight','b');
			set(h,'Rotation',90);
		end
		ylabel('Sea Level Contribution (mm SLE)');
		%title('NorESM');
		xlim([0.5 3.5]);
		ylim([-40 160]);
		%set(gca,'position',[0.15 0.47 0.8 0.45])
		hold on

		%Add SMB results
		resultstendacabfgr=[resultstendacabfgr_west,resultstendacabfgr_east,resultstendacabfgr_penin];
		for ib = 1:numel(hb)
			%XData property is the tick labels/group centers; XOffset is the offset
			%of each distinct group
			xData = hb(ib).XData+hb(ib).XOffset;
			plot(xData(1),resultstendacabfgr_west(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
			plot(xData(2),resultstendacabfgr_east(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
			plot(xData(3),resultstendacabfgr_penin(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
		end
		%bar(resultstendacbf*(-1/361.8),'grouped','color','k');
		legend(results_model,'location','EastOutside')
		error
		export_fig('Figures/NonLinearityVolume.pdf')
	end

end %}}}
if perform(org,'Plot_NonLinearity_VAF')% {{{

	model_list={'ILTS_PIK_SICOPOLIS','JPL1_ISSM','LSCE_GRISLI','NCAR_CISM','UCIJPL_ISSM'};
	colors = distinguishable_colors(20);
	%experiments_list={'expC2','expC3','expC5','expC6','expC8','expC9','expC11','expC12',... %open and standard NorESM
	experiments_list={'expC2','expC3','expC5','expC6','expC11','expC12',... %open and standard NorESM
		};

	results=[]; results_baseline=[];
	resultstendacabf=[]; resultstendacabf_baseline=[];

	results_model={};
	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					explim_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabf_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explim_file,'time');
					lim=ncread(explim_file,'ivaf')*ice_density/(10^9*1000); %from m^3 to Gt;
					tendacabf=ncread(exptendacabf_file,'smbgr')*yearday_model*3600*24/(10^9*1000); %from kg/s to Gt/yr;
					results(end+1)=lim(85);
					resultstendacabf(end+1)=sum(tendacabf(1:85));
					if length(time_model)>90, error('too long '); end
					if strcmpi(expename,'expC2'),
						expename='exp01';
						results_model{end+1}=[model_list2{imodel} '\_NorESM'] ;
					elseif strcmpi(expename,'expC3'),
						expename='exp05';
						results_model{end+1}=[model_list2{imodel} '\_NorESM'] ;
					elseif strcmpi(expename,'expC5'),
						expename='exp02';
						results_model{end+1}=[model_list2{imodel} '\_MIROC'] ;
					elseif strcmpi(expename,'expC6'),
						expename='exp06';
						results_model{end+1}=[model_list2{imodel} '\_MIROC'] ;
					elseif strcmpi(expename,'expC8'),
						expename='exp03';
						results_model{end+1}=[model_list2{imodel} '\_NorESM\_2.6'] ;
					elseif strcmpi(expename,'expC9'),
						expename='exp07';
						results_model{end+1}=[model_list2{imodel} '\_NorESM\_2.6'] ;
					elseif strcmpi(expename,'expC11'),
						expename='exp04';
						results_model{end+1}=[model_list2{imodel} '\_CCSM4'] ;
					elseif strcmpi(expename,'expC12'),
						expename='exp08';
						results_model{end+1}=[model_list2{imodel} '\_CCSM4'] ;
					end
					explim_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabf_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explim_file,'time');
					lim=ncread(explim_file,'ivaf')*ice_density/(10^9*1000); %from m^3 to Gt;
					tendacabf=ncread(exptendacabf_file,'smbgr')*yearday_model*3600*24/(10^9*1000); %from kg/s to Gt/yr;
					results_baseline(end+1)=lim(85);
					resultstendacabf_baseline(end+1)=sum(tendacabf(1:85));
				end

			end %end of model
		end %end of isexp
		allresults=[results_baseline; results; results_baseline-resultstendacabf_baseline]'; 

		%Plot results
		figure(2); clf; set(gcf,'color','w'); set(gcf,'Position',[400 500 900 400]);
		hb=bar(allresults*(-1/362.5),'grouped'); hold on 
		set(gca, 'fontsize',14);
		legend('All forcings','Ocean-only','Ocean-calculated','location','NorthEast','fontsize',14)
		%set(gca, 'XTick', 1:size(allresults,1) );
		%set(gca, 'XTicklabels', results_model ,'fontsize',8);
		%xtickangle(90)
		set(gca, 'XTickLabel', [])
		ylabel('Ice mass above floatation loss (mm SLE)','fontsize',16);
		set(gca, 'YGrid', 'on', 'XGrid', 'off')
		ylim([-100 100]);
		text(-1.0,-100,'b','fontsize',18,'fontweight','b');
		export_fig('Figures/NonLinearityVAF.pdf')
		error
		h=text(1,-48,'WAIS','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		h=text(2,-48,'EAIS','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		h=text(3,-48,'Peninsula','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		for iresults=1:length(results_model),
			bi=bar(iresults,results(iresults)*(-1/361.8),'facecolor',colors(iresults,:));
			h=text(iresults,-42,results{iresults},'VerticalAlignment','middle','HorizontalAlignment','right','color',colors(iresults,:),'fontweight','b');
			set(h,'Rotation',90);
		end
		ylabel('Sea Level Contribution (mm SLE)');
		%title('NorESM');
		xlim([0.5 3.5]);
		ylim([-40 160]);
		%set(gca,'position',[0.15 0.47 0.8 0.45])
		hold on

		%Add SMB results
		resultstendacabfgr=[resultstendacabfgr_west,resultstendacabfgr_east,resultstendacabfgr_penin];
		for ib = 1:numel(hb)
			%XData property is the tick labels/group centers; XOffset is the offset
			%of each distinct group
			xData = hb(ib).XData+hb(ib).XOffset;
			plot(xData(1),resultstendacabfgr_west(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
			plot(xData(2),resultstendacabfgr_east(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
			plot(xData(3),resultstendacabfgr_penin(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
		end
		%bar(resultstendacbf*(-1/361.8),'grouped','color','k');
		legend(results_model,'location','EastOutside')
		error
		export_fig('Figures/NonLinearityVAF.pdf')
	end

end %}}}
if perform(org,'Plot_NonLinearity_Melt')% {{{

	model_list={'ILTS_PIK_SICOPOLIS','JPL1_ISSM','LSCE_GRISLI','NCAR_CISM','UCIJPL_ISSM'};
	colors = distinguishable_colors(20);
	experiments_list={'expC2','expC3','expC5','expC6','expC11','expC12',... %open and standard NorESM
		};
	%experiments_list={'expC2','expC3','expC5','expC6','expC8','expC9','expC11','expC12',... %open and standard NorESM

	results=[]; results_baseline=[];
	resultsbmbfl=[]; resultsbmbfl_baseline=[];

	results_model={};
	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					%expbmbfl_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_bmbfl_AIS_' group '_' simul '_' expename '.nc'];
					expbmbfl_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_bmbfl_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(expbmbfl_file,'time');
					bmbfl=ncread(expbmbfl_file,'bmbfl')*yearday_model*3600*24/(10^9*1000); %from kg/s to Gt/yr;
					resultsbmbfl(end+1)=sum(bmbfl(1:85));
					if length(time_model)>90, error('too long '); end
					if strcmpi(expename,'expC2'),
						expename='exp01';
						results_model{end+1}=[model_list2{imodel} '\_open\_NorESM'] ;
					elseif strcmpi(expename,'expC3'),
						expename='exp05';
						results_model{end+1}=[model_list2{imodel} '\_standard\_NorESM'] ;
					elseif strcmpi(expename,'expC5'),
						expename='exp02';
						results_model{end+1}=[model_list2{imodel} '\_open\_MIROC'] ;
					elseif strcmpi(expename,'expC6'),
						expename='exp06';
						results_model{end+1}=[model_list2{imodel} '\_standard\_MIROC'] ;
					elseif strcmpi(expename,'expC8'),
						expename='exp03';
						results_model{end+1}=[model_list2{imodel} '\_open\_NorESM\_2.6'] ;
					elseif strcmpi(expename,'expC9'),
						expename='exp07';
						results_model{end+1}=[model_list2{imodel} '\_standard\_NorESM\_2.6'] ;
					elseif strcmpi(expename,'expC11'),
						expename='exp04';
						results_model{end+1}=[model_list2{imodel} '\_open\_CCSM4'] ;
					elseif strcmpi(expename,'expC12'),
						expename='exp08';
						results_model{end+1}=[model_list2{imodel} '\_standard\_CCSM4'] ;
					end
					%expbmbfl_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_bmbfl_AIS_' group '_' simul '_' expename '.nc'];
					expbmbfl_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_bmbfl_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(expbmbfl_file,'time');
					bmbfl=ncread(expbmbfl_file,'bmbfl')*yearday_model*3600*24/(10^9*1000); %from kg/s to Gt/yr;
					resultsbmbfl_baseline(end+1)=sum(bmbfl(1:85));
				end

			end %end of model
		end %end of isexp
		allresults=[resultsbmbfl; resultsbmbfl_baseline]';

		%Plot results
		figure(2); clf; set(gcf,'color','w'); set(gcf,'Position',[400 400 900 500]);
		hb=bar(allresults*(-1),'grouped'); hold on  %keep in Gt
		set(gca, 'fontsize',14);
		legend('All forcings','Ocean-only','location','NorthEast','fontsize',14);
		set(gca, 'XTick', 1:size(allresults,1) );
		set(gca, 'XTicklabels', results_model ,'fontsize',14);
		xtickangle(90)
		ylabel('Sub-ice shelf melt (Gt)','fontsize',16);
		set(gca, 'YGrid', 'on', 'XGrid', 'off')
		text(-0.8,-1.5*10^4,'c','fontsize',18,'fontweight','b');
		export_fig('Figures/NonLinearityBasalMelt.pdf')
		error
		h=text(1,-48,'WAIS','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		h=text(2,-48,'EAIS','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		h=text(3,-48,'Peninsula','VerticalAlignment','middle','HorizontalAlignment','center','color','k','fontweight','b');
		for iresults=1:length(results_model),
			bi=bar(iresults,results(iresults)*(-1/361.8),'facecolor',colors(iresults,:));
			h=text(iresults,-42,results{iresults},'VerticalAlignment','middle','HorizontalAlignment','right','color',colors(iresults,:),'fontweight','b');
			set(h,'Rotation',90);
		end
		ylabel('Sea Level Contribution (mm SLE)');
		%title('NorESM');
		xlim([0.5 3.5]);
		ylim([-40 160]);
		%set(gca,'position',[0.15 0.47 0.8 0.45])
		hold on
		text(-1,0,'c','fontsize',14,'fontweight','b');
		error
		export_fig('Figures/NonLinearityBasalMelt.pdf')

		%Add SMB results
		resultsbmbflgr=[resultsbmbflgr_west,resultsbmbflgr_east,resultsbmbflgr_penin];
		for ib = 1:numel(hb)
			%XData property is the tick labels/group centers; XOffset is the offset
			%of each distinct group
			xData = hb(ib).XData+hb(ib).XOffset;
			plot(xData(1),resultsbmbflgr_west(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
			plot(xData(2),resultsbmbflgr_east(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
			plot(xData(3),resultsbmbflgr_penin(ib)*(-1/362.5),'d','MarkerFaceColor',colors(ib,:),'MarkerSize',6,'MarkerEdgeColor','k')
		end
		%bar(resultstendacbf*(-1/361.8),'grouped','color','k');
		legend(results_model,'location','EastOutside')
		error
		export_fig('Figures/NonLinearityBasalMelt.pdf')
	end

end %}}}

%Figure 2
if perform(org,'ComparisonEmulatedSimulatedTestRuns')% {{{

	clf;
	load Data/EmulationData/Antarctica/test_emulated.mat
	load Data/EmulationData/Antarctica/ant_dyn_loss_emulated
	time=2016:2100;
	resh_test=reshape(test_values,[85,14]);
	resh_emul=reshape(test_emulated,[85,14]);
	plot(time,-resh_test,'linewidth',2)
	hold on
	plot(time,-resh_emul,':','linewidth',2)
	set(gcf,'color','w')
	xlim([2016 2100])
	grid on
	ax = gca;
   ax.FontSize = 16;
	xlabel('Time (yrs)','fontsize',18)
	ylabel({'Dynamic sea level contribution', '(mm SLE)'},'fontsize',18)

	export_fig('Figures/SimulEmulTestAntarctica.pdf')

end %}}}

%Figure 3
if perform(org,'EmulatedSimulatedResults')% {{{

	load Data/EmulationData/Antarctica/ant_dyn_loss_emulated

	time=[2015:2100];
	figure(1); clear figure; set(gcf,'color','w'); hold on;
	set(gcf,'position',[100 100 800 400]);
	plot(time,-1*[zeros(1,210); individual_results_limnsw_timeseries'.*(1-emulated_results)'],'color',[128 128 128]/256,'linewidth',1.0)
	hold on
	plot(time,-1*[zeros(1,210); individual_results_limnsw_timeseries'.*(emulated_results)'],'color',[178 34 34]/256,'linewidth',1.0)
	xlabel('Time (years)','fontsize',18)
	xlim([2015 2100])
	grid on
	ax = gca;
   ax.FontSize = 16;
	ylabel({'Dynamic sea level contribution', '(mm SLE)'},'fontsize',18)
	export_fig('Figures/Emulation_Antarctica.pdf')

end %}}}

%Figure 4
if perform(org,'Plot_GlacierReduced_NorESM')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'exp05','exp01',... %open and standard NorESM
		};

	number_models=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number_models=number_models+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
						if strcmpi(modelname,'UCIJPL_ISSM') & strcmpi(expename,'exp01') & iglacier==56,
							sensitivity_results_limnsw_glacier56(end,45)=(sensitivity_results_limnsw_glacier56(end,44)+sensitivity_results_limnsw_glacier56(end,46))/2;
						end

					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor = [colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('NorESM','fontsize',18)
	text(2008,-10,'a','fontsize',20,'fontweight','b');
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	export_fig('Figures/GlaciersNorESM.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_MIROC')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'exp06','exp02',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor = [colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('MIROC','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'b','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersMIROC.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_CCSM4')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'exp08','exp04',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor=[colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('CCSM4','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'c','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersCCSM4.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_HadGEM')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'expA5','expA1',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					modelname
					expename

					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor=[colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('HadGEM2','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'d','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersHadGEM.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_CSIRO')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'expA6','expA2',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor=[colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('CSIRO','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'e','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersCSIRO.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_IPSL')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'expA7','expA3',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor=[colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('IPSL','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'f','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersIPSL.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_CNRM')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'expB6','expB1',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor=[colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('CNRM','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'g','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersCNRM.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_UKESM')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'expB8','expB3',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor=[colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('UKESM','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'h','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersUKESM.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_CESM2')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'expB9','expB4',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor=[colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('CESM2','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'i','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersCESM2.pdf')

end %}}}
if perform(org,'Plot_GlacierReduced_CNRMESM')% {{{

	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 450]);
	colors = distinguishable_colors(length(model_list)+10);
	colors(6,:)=colors(6,:)*0.93;
	colors(19,:)=colors(end,:);
	experiments_list={'expB10','expB5',... %open and standard NorESM
		};

	number=0;
	results_model={};

	for iglacier=1:198,
		eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '=[];'])
	end
	mean_sensitivity_results_limnsw_glacier=zeros(198,86);

	if 1,
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list),
				expename=experiments_list{iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					number=number+1;
					explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					time_model=ncread(explimnsw_file,'time');
					for iglacier=1:198,
						eval(['results_limnsw_glacier' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(iglacier) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
						eval(['results_tendacabfgr_glacier' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(iglacier) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
						eval(['sensitivity_results_limnsw_glacier' num2str(iglacier) '(end+1,:)=[0;results_limnsw_glacier' num2str(iglacier) '(1:85)];'])
						eval(['sensitivity_results_tendacabfgr_glacier' num2str(iglacier) '(end+1,:)=[0;cumsum(results_tendacabfgr_glacier' num2str(iglacier) '(1:85))];'])
					end
				end

			end %end of model
		end %end of isexp
	end
	figure(1); clf
	total_ant_vector=zeros(size(sensitivity_results_limnsw_glacier1));
	for iglacier=1:198,
		eval(['mean_sensitivity_results_limnsw_glacier(iglacier,:)=-mean(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ',1)/362.5;'])
		eval(['total_ant_vector=total_ant_vector-(sensitivity_results_limnsw_glacier' num2str(iglacier) '-sensitivity_results_tendacabfgr_glacier' num2str(iglacier) ')/362.5;'])
	end
	total11=mean_sensitivity_results_limnsw_glacier(7,end)+ mean_sensitivity_results_limnsw_glacier(2,end)+ mean_sensitivity_results_limnsw_glacier(83,end)+mean_sensitivity_results_limnsw_glacier(13,end)+mean_sensitivity_results_limnsw_glacier(18,end)+mean_sensitivity_results_limnsw_glacier(20,end)+mean_sensitivity_results_limnsw_glacier(186,end)+mean_sensitivity_results_limnsw_glacier(172,end)+mean_sensitivity_results_limnsw_glacier(92,end)+mean_sensitivity_results_limnsw_glacier(82,end)+mean_sensitivity_results_limnsw_glacier(25,end);
	total11_vector=-(sensitivity_results_limnsw_glacier7-sensitivity_results_tendacabfgr_glacier7)/362.5 -(sensitivity_results_limnsw_glacier2-sensitivity_results_tendacabfgr_glacier2)/362.5 -(sensitivity_results_limnsw_glacier83-sensitivity_results_tendacabfgr_glacier83)/362.5 -(sensitivity_results_limnsw_glacier13-sensitivity_results_tendacabfgr_glacier13)/362.5 -(sensitivity_results_limnsw_glacier18-sensitivity_results_tendacabfgr_glacier18)/362.5 -(sensitivity_results_limnsw_glacier20-sensitivity_results_tendacabfgr_glacier20)/362.5 -(sensitivity_results_limnsw_glacier186-sensitivity_results_tendacabfgr_glacier186)/362.5 -(sensitivity_results_limnsw_glacier172-sensitivity_results_tendacabfgr_glacier172)/362.5 -(sensitivity_results_limnsw_glacier92-sensitivity_results_tendacabfgr_glacier92)/362.5 -(sensitivity_results_limnsw_glacier82-sensitivity_results_tendacabfgr_glacier82)/362.5 -(sensitivity_results_limnsw_glacier25-sensitivity_results_tendacabfgr_glacier25)/362.5;
	total11percent=total11/sum(mean_sensitivity_results_limnsw_glacier(:,end));
	total11percent_vector=total11_vector(:,end)./total_ant_vector(:,end)*100;
	total_ant=sum(mean_sensitivity_results_limnsw_glacier,1);
	results_11glaciers=[mean_sensitivity_results_limnsw_glacier(7,:);mean_sensitivity_results_limnsw_glacier(2,:);mean_sensitivity_results_limnsw_glacier(83,:); mean_sensitivity_results_limnsw_glacier(13,:);...
		mean_sensitivity_results_limnsw_glacier(18,:);mean_sensitivity_results_limnsw_glacier(20,:);mean_sensitivity_results_limnsw_glacier(186,:);mean_sensitivity_results_limnsw_glacier(172,:); ...
		mean_sensitivity_results_limnsw_glacier(92,:);mean_sensitivity_results_limnsw_glacier(82,:);mean_sensitivity_results_limnsw_glacier(25,:)];
	results_summary=[results_11glaciers;total_ant-sum(results_11glaciers,1)];
	years=2015:2100;
	a=area(years,results_summary');
	for iglacier=1:11,
		a(iglacier).FaceColor=[colors(iglacier,:)];
		a(iglacier).EdgeColor = [colors(iglacier,:)];
	end
	a(12).FaceColor=[255 255 255]/255;
	hold on
	plot(years,total_ant,'k','LineWidth',2)
	xlim([2015 2100])
	ylim([-5 100])
	xlabel('Years','fontsize',16)
	ylabel('SLR (mm)','fontsize',16)
	title('CNRMESM','fontsize',18)
	ax = gca;
	ax.FontSize = 16;
	set(get(gca,'title'),'Position',[2057.5 92 1.00011])
	grid on
	text(2008,-10,'j','fontsize',20,'fontweight','b');
	export_fig('Figures/GlaciersCNRMESM.pdf')

end %}}}
if perform(org,'Map11glaciers')% {{{

	%Now create inset map to show the color of each region
	colorm = distinguishable_colors(186);

	glaciers_8km=rot90(double(ncread('Data/glaciers_8km_iceonly.nc','glaciers'))); %12 glaciers with most ice loss
	pos=find(glaciers_8km~=7 &  glaciers_8km~=2 & glaciers_8km~=83 & glaciers_8km~=13 & glaciers_8km~=18 & glaciers_8km~=20 & ...
		glaciers_8km~=186 & glaciers_8km~=172 & glaciers_8km~=92 & glaciers_8km~=82 & glaciers_8km~=25);
	glaciers_8km(pos)=0;
	colorm(186,:)=colorm(7,:);
	colorm(172,:)=colorm(8,:);
	colorm(92,:)=colorm(9,:);
	colorm(83,:)=colorm(3,:);
	colorm(82,:)=colorm(10,:);
	colorm(25,:)=colorm(11,:);
	colorm(20,:)=colorm(6,:);
	colorm(18,:)=colorm(5,:);
	colorm(13,:)=colorm(4,:);
	colorm(2,:)=colorm(2,:);
	colorm(7,:)=colorm(1,:);

	figure(1);
	set(gcf,'position',[100 100 300 300])
	patch([0 0 0],[0 0 0],colorm(7,:));
	patch([0 0 0],[0 0 0],colorm(2,:));
	patch([0 0 0],[0 0 0],colorm(83,:));
	patch([0 0 0],[0 0 0],colorm(13,:));
	patch([0 0 0],[0 0 0],colorm(18,:));
	patch([0 0 0],[0 0 0],colorm(20,:));
	patch([0 0 0],[0 0 0],colorm(186,:));
	patch([0 0 0],[0 0 0],colorm(172,:));
	patch([0 0 0],[0 0 0],colorm(92,:));
	patch([0 0 0],[0 0 0],colorm(82,:));
	patch([0 0 0],[0 0 0],colorm(25,:));
	patch([0 0 0],[0 0 0],'w','FaceColor','w','EdgeColor','w');
	axis('off')
	l=legend('Thwaites','Pine Island','Totten','Getz','Kamb','MacAyeal','Institute','Moller','Leopold and Astrid Coast','Moscow','Whillans',...
		'Location','EastOutside','box','off','fontsize',18)
	l.Position=[0 0 1 1];
	set(gcf,'color','w');
	h = gcf;
	set(h,'Units','Inches');
	pos = get(h,'Position');
	set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
	print(gcf, '-dpdf', '-painters', 'Figures/AntarcticMap12GlaciersLegend.pdf');

	error
	[pos_nani pos_nanj]=find(isnan(glaciers_8km) | glaciers_8km==0);
	data_min=1; data_max=187;
	image_rgb = ind2rgb(uint16((max(data_min,min(data_max,glaciers_8km)) - data_min)*(size(colorm,1)/(data_max-data_min))),colorm);
	image_rgb(sub2ind(size(image_rgb),repmat(pos_nani,1,3),repmat(pos_nanj,1,3),repmat(1:3,size(pos_nani,1),1))) = repmat([1 1 1],size(pos_nani,1),1);
	imagesc(-3040:6080/size(glaciers_8km,2):3040,-3040:6080/size(glaciers_8km,1):3040,image_rgb); colorbar('off');
	set(gca,'fontsize',14);
	colorbar('off'); axis('equal','off'); xlim([-3040 3040]); ylim([-3040 3040]);
	expdisp('./Data/IceContour.exp','linewidth',1,'linestyle','k');
	expdisp('./Data/IceShelvesContour.exp','linewidth',1,'linestyle','k');
	xlim([-2900 2900])
	ylim([-2250 2160])
	set(gca,'position',[0.01 0.01 0.49 0.99]);
	h = gcf;
	set(h,'Units','Inches');
	pos = get(h,'Position');
	set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
	print(gcf, '-dpdf', '-painters', 'Figures/AntarcticMap12Glaciers.pdf');

end %}}}

%Figure 5
if perform(org,'Plot_GlacierSensitivity_ByClimateModel')% {{{

	colors = distinguishable_colors(12); %number of basins
	experiments_list={'exp01','exp02','exp04','exp05','exp06','exp08','expA1','expA2','expA3','expA5','expA6','expA7','expB1','expB3','expB4','expB5','expB6','expB8','expB9','expB10'};
	%experiments_list={'exp01','exp02','exp04','exp05','exp06','exp08','expA1','expA2','expA3','expA5','expA6','expA7','expB1','expB3','expB4','expB5','expB6','expB8','expB9','expB10','exp09','exp10','exp13','expD1','expD2','expD3','expD4','expD5','expD6','expD7','expD8','expD9','expD10','expD11','expD12','expD13','expD14','expD15','expD16','expD17','expD18','expD51','expD52','expD53','expD54','expD55','expD56','expD57','expD58'}; %open and standard RCP8.5 & ssp585

	glacier_list=[7 2 83 13 18 20 186 172 92 82 25];

	if 1,
	for iglacier=1:11,
		eval(['sensitivity_results_limnsw_noresm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_noresm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_noresm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_limnsw_miroc_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_miroc_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_miroc_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_limnsw_ccsm4_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_ccsm4_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_ccsm4_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_limnsw_hadgem2_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_hadgem2_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_hadgem2_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_limnsw_csiro_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_csiro_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_csiro_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_limnsw_ipsl_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_ipsl_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_ipsl_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_limnsw_cnrm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_cnrm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_cnrm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_limnsw_ukesm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_ukesm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_ukesm_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_limnsw_cesm2_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendacabfgr_cesm2_sector' num2str(iglacier) '=[];'])
		eval(['sensitivity_results_tendlibmassbffl_cesm2_sector' num2str(iglacier) '=[];'])
	end

	for imodel=1:length(model_list),
		modelname=model_list{imodel};
		for iexp=1:length(experiments_list),
			expename=experiments_list{iexp};
			specifics;

			eval(['isexp=is' expename ';'])
			if isexp,
				explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
				exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
				exptendlibmassbffl_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_bmbfl_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];

				for iglacier=1:11,
					indexglaciernetcdf=glacier_list(iglacier);
					eval(['results_tendlibmassbffl_sector' num2str(iglacier) '=ncread(exptendlibmassbffl_file,''bmbfl_glacier_' num2str(indexglaciernetcdf) ''') *yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
					eval(['results_tendacabfgr_sector' num2str(iglacier) '=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(indexglaciernetcdf) ''')*yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
					eval(['results_limnsw_sector' num2str(iglacier) '=ncread(explimnsw_file,''ivaf_glacier_' num2str(indexglaciernetcdf) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt

					if strcmpi(expename,'expB5') | strcmpi(expename,'expB10'),
						% do nothing removed from list
					elseif strcmpi(expename,'expB4') | strcmpi(expename,'expB9') | strcmpi(expename,'expD17') | strcmpi(expename,'expD18'),
						eval(['sensitivity_results_tendacabfgr_cesm2_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_cesm2_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_cesm2_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					elseif strcmpi(expename,'expB3') | strcmpi(expename,'expB8') | strcmpi(expename,'expD15') | strcmpi(expename,'expD16'),
						eval(['sensitivity_results_tendacabfgr_ukesm_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_ukesm_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_ukesm_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					elseif strcmpi(expename,'expB1') | strcmpi(expename,'expB2') | strcmpi(expename,'expB6') | strcmpi(expename,'expB7') | strcmpi(expename,'expD13') | strcmpi(expename,'expD14'),
						eval(['sensitivity_results_tendacabfgr_cnrm_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_cnrm_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_cnrm_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					elseif strcmpi(expename,'expA3') | strcmpi(expename,'expA4') | strcmpi(expename,'expA7') | strcmpi(expename,'expA8') | strcmpi(expename,'expD11') | strcmpi(expename,'expD12'),
						eval(['sensitivity_results_tendacabfgr_ipsl_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_ipsl_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_ipsl_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					elseif strcmpi(expename,'expA2') | strcmpi(expename,'expA6') | strcmpi(expename,'expD9') | strcmpi(expename,'expD10'),
						eval(['sensitivity_results_tendacabfgr_csiro_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_csiro_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_csiro_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					elseif strcmpi(expename,'expA1') | strcmpi(expename,'expA5') | strcmpi(expename,'expD7') | strcmpi(expename,'expD8'),
						eval(['sensitivity_results_tendacabfgr_hadgem2_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_hadgem2_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_hadgem2_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					elseif strcmpi(expename,'exp04') | strcmpi(expename,'exp08') | strcmpi(expename,'exp11') | strcmpi(expename,'exp12') | strcmpi(expename,'expD5') | strcmpi(expename,'expD6') | strcmpi(expename,'expD56') | strcmpi(expename,'expD57') | strcmpi(expename,'expD58'),
						eval(['sensitivity_results_tendacabfgr_ccsm4_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_ccsm4_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_ccsm4_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					elseif strcmpi(expename,'exp02') | strcmpi(expename,'exp06') | strcmpi(expename,'expD1') | strcmpi(expename,'expD2') | strcmpi(expename,'expD53') | strcmpi(expename,'expD54') | strcmpi(expename,'expD55'),
						eval(['sensitivity_results_tendacabfgr_ccsm4_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_ccsm4_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_ccsm4_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					else
						eval(['sensitivity_results_tendacabfgr_noresm_sector' num2str(iglacier) '(end+1)=sum(results_tendacabfgr_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_tendlibmassbffl_noresm_sector' num2str(iglacier) '(end+1)=sum(results_tendlibmassbffl_sector' num2str(iglacier) '(1:85));'])
						eval(['sensitivity_results_limnsw_noresm_sector' num2str(iglacier) '(end+1)=results_limnsw_sector' num2str(iglacier) '(85);'])
					end
				end
			end

		end %end of model
	end %end of isexp
	end

	for iglacier=1:11,
		eval(['melt_' num2str(iglacier) '=transp([-sensitivity_results_tendlibmassbffl_noresm_sector' num2str(iglacier) ' , -sensitivity_results_tendlibmassbffl_miroc_sector' num2str(iglacier) ' , -sensitivity_results_tendlibmassbffl_ccsm4_sector' num2str(iglacier) ' , -sensitivity_results_tendlibmassbffl_hadgem2_sector' num2str(iglacier) ' , -sensitivity_results_tendlibmassbffl_csiro_sector' num2str(iglacier) ' , -sensitivity_results_tendlibmassbffl_ipsl_sector' num2str(iglacier) ' , -sensitivity_results_tendlibmassbffl_cnrm_sector' num2str(iglacier) ' , -sensitivity_results_tendlibmassbffl_ukesm_sector' num2str(iglacier) ' , -sensitivity_results_tendlibmassbffl_cesm2_sector' num2str(iglacier) ']) ;'])
		eval(['vaf_' num2str(iglacier) '=transp([-(sensitivity_results_limnsw_noresm_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_noresm_sector' num2str(iglacier) ')/362.5 , -(sensitivity_results_limnsw_miroc_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_miroc_sector' num2str(iglacier) ')/362.5 , -(sensitivity_results_limnsw_ccsm4_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_ccsm4_sector' num2str(iglacier) ')/362.5 , -(sensitivity_results_limnsw_hadgem2_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_hadgem2_sector' num2str(iglacier) ')/362.5 , -(sensitivity_results_limnsw_csiro_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_csiro_sector' num2str(iglacier) ')/362.5 , -(sensitivity_results_limnsw_ipsl_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_ipsl_sector' num2str(iglacier) ')/362.5 , -(sensitivity_results_limnsw_cnrm_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_cnrm_sector' num2str(iglacier) ')/362.5 , -(sensitivity_results_limnsw_ukesm_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_ukesm_sector' num2str(iglacier) ')/362.5 , -(sensitivity_results_limnsw_cesm2_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_cesm2_sector' num2str(iglacier) ')/362.5]);'])
%		eval(['pos=find(melt_' num2str(iglacier) '>0) ; '])
%		eval(['melt_' num2str(iglacier) '=melt_' num2str(iglacier) '(pos) ; '])
%		eval(['vaf_' num2str(iglacier) '=vaf_' num2str(iglacier) '(pos) ; '])
		eval([' Melt_' num2str(iglacier) ' = [ melt_' num2str(iglacier) ' ones(length( melt_' num2str(iglacier) '), 1)] ;'])
		eval(['coeff_' num2str(iglacier) '= Melt_' num2str(iglacier) '\vaf_' num2str(iglacier) ' ;'])
		eval([' vaf_data_' num2str(iglacier) ' =Melt_' num2str(iglacier) ' *coeff_' num2str(iglacier) ' ;'])
		eval(['melt_calc_' num2str(iglacier) '=[-1:6]*10^4 ;'])
		eval(['p_' num2str(iglacier) '=polyfit(melt_' num2str(iglacier) ',vaf_' num2str(iglacier) ',1) ;'])
		eval([' vaf_calc_' num2str(iglacier) '=polyval(p_' num2str(iglacier) ',melt_calc_' num2str(iglacier) ') ;'])
		eval([' vaf_data_' num2str(iglacier) '=polyval(p_' num2str(iglacier) ',melt_' num2str(iglacier) ') ;'])
		eval([' Rsq_' num2str(iglacier) '=1-sum((vaf_' num2str(iglacier) '-vaf_data_' num2str(iglacier) ').^2)/sum((vaf_' num2str(iglacier) '-mean(vaf_data_' num2str(iglacier) ')).^2) ;'])
	end

	%Plot results per region
	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 500]);
	plot(-2*10^4,-20,'o','color','k'); hold on
	plot(-2*10^4,-20,'+','color','k');
	plot(-2*10^4,-20,'*','color','k');
	plot(-2*10^4,-20,'.','color','k');
	plot(-2*10^4,-20,'x','color','k');
	plot(-2*10^4,-20,'s','color','k');
	plot(-2*10^4,-20,'d','color','k');
	plot(-2*10^4,-20,'v','color','k');
	plot(-2*10^4,-20,'^','color','k');
	hold on
	for iglacier=1:11,
		eval(['plot(-sensitivity_results_tendlibmassbffl_noresm_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_noresm_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_noresm_sector' num2str(iglacier) ')/362.5,''o'',''color'',colors(iglacier,:));'])
		hold on
		eval(['plot(-sensitivity_results_tendlibmassbffl_miroc_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_miroc_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_miroc_sector' num2str(iglacier) ')/362.5,''+'',''color'',colors(iglacier,:));'])
		eval(['plot(-sensitivity_results_tendlibmassbffl_ccsm4_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_ccsm4_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_ccsm4_sector' num2str(iglacier) ')/362.5,''*'',''color'',colors(iglacier,:));'])
		eval(['plot(-sensitivity_results_tendlibmassbffl_hadgem2_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_hadgem2_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_hadgem2_sector' num2str(iglacier) ')/362.5,''.'',''color'',colors(iglacier,:));'])
		eval(['plot(-sensitivity_results_tendlibmassbffl_csiro_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_csiro_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_csiro_sector' num2str(iglacier) ')/362.5,''x'',''color'',colors(iglacier,:));'])
		eval(['plot(-sensitivity_results_tendlibmassbffl_ipsl_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_ipsl_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_ipsl_sector' num2str(iglacier) ')/362.5,''s'',''color'',colors(iglacier,:));'])
		eval(['plot(-sensitivity_results_tendlibmassbffl_cnrm_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_cnrm_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_cnrm_sector' num2str(iglacier) ')/362.5,''d'',''color'',colors(iglacier,:));'])
		eval(['plot(-sensitivity_results_tendlibmassbffl_ukesm_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_ukesm_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_ukesm_sector' num2str(iglacier) ')/362.5,''v'',''color'',colors(iglacier,:));'])
		eval(['plot(-sensitivity_results_tendlibmassbffl_cesm2_sector' num2str(iglacier) ',-(sensitivity_results_limnsw_cesm2_sector' num2str(iglacier) '-sensitivity_results_tendacabfgr_cesm2_sector' num2str(iglacier) ')/362.5,''^'',''color'',colors(iglacier,:));'])
	end
	ylabel('Dynamic mass loss (mm SLE)','fontsize',16);
	xlabel('Ice shelf melt (Gt)','fontsize',16);
	xlim([-0.1 0.6]*10^5)
	ylim([-10 60])
	%text(-1.4*10^4,-12,'a','VerticalAlignment','middle','HorizontalAlignment','right','fontsize',18,'fontweight','b');
	legend('NorESM','MIROC','CCSM4','HadGEM2','CSIRO','IPSL','CNRM','UKESM','CESM2','fontsize',16,'location','NorthWest')
	ax = gca;
   ax.FontSize = 16;

	cplotpos=get(gca,'pos');
	%compute inset position
	PosInset=[0.45, 0.37, 0.55, 0.55];
	axes('pos',PosInset);
	axis equal off

	%Now create inset map to show the color of each region
	glaciers_8km=rot90(double(ncread('Data/glaciers_8km_iceonly.nc','glaciers'))); %12 glaciers with most ice loss
	pos=find(glaciers_8km~=7 &  glaciers_8km~=2 & glaciers_8km~=83 & glaciers_8km~=13 & glaciers_8km~=18 & glaciers_8km~=20 & ...
		glaciers_8km~=186 & glaciers_8km~=172 & glaciers_8km~=92 & glaciers_8km~=82 & glaciers_8km~=25);
	glaciers_8km(pos)=0;
	colorm=distinguishable_colors(186);
	colorm(186,:)=colorm(7,:);
	colorm(172,:)=colorm(8,:);
	colorm(97,:)=colorm(12,:);
	colorm(92,:)=colorm(9,:);
	colorm(83,:)=colorm(3,:);
	colorm(82,:)=colorm(10,:);
	colorm(25,:)=colorm(11,:);
	colorm(20,:)=colorm(6,:);
	colorm(18,:)=colorm(5,:);
	colorm(13,:)=colorm(4,:);
	colorm(2,:)=colorm(2,:);
	colorm(7,:)=colorm(1,:);

	set(gcf,'color','w');
	[pos_nani pos_nanj]=find(isnan(glaciers_8km) | glaciers_8km==0);
	data_min=1; data_max=187;
	image_rgb = ind2rgb(uint16((max(data_min,min(data_max,glaciers_8km)) - data_min)*(size(colorm,1)/(data_max-data_min))),colorm);
	image_rgb(sub2ind(size(image_rgb),repmat(pos_nani,1,3),repmat(pos_nanj,1,3),repmat(1:3,size(pos_nani,1),1))) = repmat([1 1 1],size(pos_nani,1),1);
	imagesc(-3040:6080/size(glaciers_8km,2):3040,-3040:6080/size(glaciers_8km,1):3040,image_rgb); colorbar('off');
	set(gca,'fontsize',14);
	colorbar('off'); axis('equal','off'); xlim([-3040 3040]); ylim([-3040 3040]);
	expdisp('./Data/IceContour.exp','linewidth',1,'linestyle','k');
	expdisp('./Data/IceShelvesContour.exp','linewidth',1,'linestyle','k');
	colorbar('off');

	h = gcf;
	set(h,'Units','Inches');
	pos = get(h,'Position');
	set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
	error
	print(gcf, '-dpdf', '-painters', 'Figures/Figure15GlaciersByClimateModel.pdf');

end %}}}

%Figure 6
if perform(org,'VarianceRCP85NoSMBTotal_emulator_anovan')% {{{

	if 0,
		load Data/EmulationData/Antarctica/ant_dyn_loss_emulated  

		icesheetmodels=icesheetmodellong_results(:,1);
		climatemodels=climatemodel_results(:,1);
		var_ensemble=var(-individual_results_limnsw_timeseries,1);
		var_ensemble_anova=zeros(85,1);
		var_climate=zeros(85,1);
		var_ice=zeros(85,1);
		var_distrib=zeros(85,1);
		var_4=zeros(85,1);
		var_5=zeros(85,1);
		var_6=zeros(85,1);
		var_7=zeros(85,1);
		var_error=zeros(85,1);
		num_sample_distrib=100;
		for yeari=1:85,
			yeari
			distribution=repmat(1:num_sample_distrib,210,1); distribution=distribution(:);
			individual_results=-individual_results_limnsw_timeseries(:,yeari);
			individual_std=individual_results_std_timeseries(:,yeari);
			distribution_results=normrnd(repmat(individual_results,num_sample_distrib,1),repmat(individual_std,num_sample_distrib,1));
			[p,tabl]=anovan(distribution_results,{repmat(icesheetmodels,num_sample_distrib,1),repmat(climatemodels,num_sample_distrib,1),distribution},'model',2,'display','off') ;
			var_ice(yeari)=tabl{2,2}/(210*num_sample_distrib); %first line and column are names of variables and parameters
			var_climate(yeari)=tabl{3,2}/(210*num_sample_distrib);
			var_distrib(yeari)=tabl{4,2}/(210*num_sample_distrib);
			var_4(yeari)=tabl{5,2}/(210*num_sample_distrib);
			var_5(yeari)=tabl{6,2}/(210*num_sample_distrib);
			var_6(yeari)=tabl{7,2}/(210*num_sample_distrib);
			var_error(yeari)=tabl{8,2}/(210*num_sample_distrib);
			var_ensemble_anova(yeari)=tabl{end,2}/(210*num_sample_distrib);
		end
		save Data/EmulationData/Antarctica/var_ant_data_peter var_ensemble_anova var_ice var_climate var_distrib var_4 var_5 var_6 var_error
	end

	years=2016:2100;
	if 1,
		load Data/EmulationData/Antarctica/var_ant_data
		figure(3);
		set(gcf,'color','w')
		set(gcf, 'position',[100 100 650 400])
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		a=area(years,([var_climate,var_ice,var_4,var_emulator_total]./var_ensemble_anova*100));
		a(1).FaceColor=[0.8500 0.2550 0.0980];
		a(2).FaceColor=[0 0.5570 0.7410]
		a(4).FaceColor=a(3).FaceColor;
		a(3).FaceColor=[157 102 255]/255;
		legend('climate','ice','ice-climate','emulation','location','East')
		set(gca,'fontsize',14);
		xlim([2016 2100])
		ylim([0 100])
		xlabel('Years','fontsize',14)
		ylabel('% of \sigma^2_{tot}','fontsize',14)
		text(2008,0,'b','fontsize',18,'fontweight','b');
		grid on
		ax = gca;
		ax.Layer = 'top';
		error
		export_fig('Figures/RelativeVarianceAntarctica.pdf')

	end
	if 1,
		load Data/EmulationData/Antarctica/var_ant_data
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		set(gcf,'color','w')
		set(gcf, 'position',[100 100 650 400])
		plot(years,sqrt(var_ensemble_anova),'k','linewidth',2); hold off
		hold on
		plot(years,sqrt(var_climate),'linewidth',2,'color',[0.8500 0.2550 0.0980])
		plot(years,sqrt(var_ice),'linewidth',2,'color',[0 0.5570 0.7410])
		plot(years,sqrt(var_4),'linewidth',2,'color',[157 102 255]/255)
		plot(years,sqrt(var_emulator_total),'linewidth',2,'color',[0.9290 0.6940 0.1250])
		legend('\sigma total','\sigma climate','\sigma ice','\sigma ice-climate','\sigma emulation','location','NorthWest','fontsize',14,'box','off')
		xlim([2016 2100])
		ylim([0 50])
		xlabel('Years','fontsize',14)
		ylabel('\sigma (mm SLE)','fontsize',18)
		text(2010,0,'a','fontsize',18,'fontweight','b');
		set(gca,'fontsize',14);
		grid on
		ax = gca;
		ax.Layer = 'top';

		export_fig('Figures/UncertaintyAntarctica.pdf')
	end

end %}}}

%Figure 7 and 8
if perform(org,'VarianceRCP85NoSMBTotal_emulator_anovan_glaciers')% {{{

	glacier_list=[7 2 83 13 18 20 186 172 92 82 25];
	glacier_names={'Thwaites','Pine Island','Totten','Getz','Kamb','MacAyeal','Institute',...
		'Moller','Leopold and Astrid Coast','Moscow','Whillans'};
	for iglacier=11,
		close all
		indexglaciernetcdf=glacier_list(iglacier);

		if 0,
			eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/glacier_dyn_loss_emulated'])

			icesheetmodels=icesheetmodellong_results(:,1);
			climatemodels=climatemodel_results(:,1);
			var_ensemble=var(-individual_results_limnsw_timeseries,1);
			var_ensemble_anova=zeros(85,1);
			var_climate=zeros(85,1);
			var_ice=zeros(85,1);
			var_distrib=zeros(85,1);
			var_4=zeros(85,1);
			var_5=zeros(85,1);
			var_6=zeros(85,1);
			var_7=zeros(85,1);
			var_error=zeros(85,1);
			num_sample_distrib=100;
			for yeari=1:85,
				yeari
				distribution=repmat(1:num_sample_distrib,210,1); distribution=distribution(:);
				individual_results=-individual_results_limnsw_timeseries(:,yeari);
				individual_std=individual_results_std_timeseries(:,yeari);
				distribution_results=normrnd(repmat(individual_results,num_sample_distrib,1),repmat(individual_std,num_sample_distrib,1));
				[p,tabl]=anovan(distribution_results,{repmat(icesheetmodels,num_sample_distrib,1),repmat(climatemodels,num_sample_distrib,1),distribution},'model',2,'display','off') ;
				var_ice(yeari)=tabl{2,2}/(210*num_sample_distrib); %first line and column are names of variables and parameters
				var_climate(yeari)=tabl{3,2}/(210*num_sample_distrib);
				var_distrib(yeari)=tabl{4,2}/(210*num_sample_distrib);
				var_4(yeari)=tabl{5,2}/(210*num_sample_distrib);
				var_5(yeari)=tabl{6,2}/(210*num_sample_distrib);
				var_6(yeari)=tabl{7,2}/(210*num_sample_distrib);
				var_error(yeari)=tabl{8,2}/(210*num_sample_distrib);
				var_ensemble_anova(yeari)=tabl{end,2}/(210*num_sample_distrib);
			end
			error
			eval(['save Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/var_glacier_' num2str(indexglaciernetcdf) ' var_ensemble_anova var_ice var_climate var_distrib var_4 var_5 var_6 var_error'])
		end
	end

	if 0, %legend
		set(gcf, 'position',[100 100 300 400])
		set(gcf,'color','w')
		a=area(years,[var_climate,var_ice,var_4,var_emulator_total]./var_ensemble_anova*100);
		a(1).FaceColor=[0.8500 0.2550 0.0980];
		a(2).FaceColor=[0 0.5570 0.7410];
		a(4).FaceColor=a(3).FaceColor;
		a(3).FaceColor=[157 102 255]/255;
		a(1).BaseLine.Visible = 'off';
		xlim([0 1])
		ylim([0 1])
		axis off
		l=legend('climate','ice','ice-climate','emulation','fontsize',22)
		h = gcf;
		set(h,'Units','Inches');
      pos = get(h,'Position');
      set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)]);
		%print(gcf, '-dpdf', '-painters', './Figures/RelativeVariance_Legend.pdf');
	end

	years=2016:2100;
	if 1,
		eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/var_glacier_' num2str(indexglaciernetcdf) ' var_ensemble_anova var_ice var_climate var_distrib var_4 var_5 var_6 var_error'])
		figure(3);
		set(gcf,'color','w')
		set(gcf, 'position',[100 100 650 400])
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		a=area(years,[var_climate,var_ice,var_4,var_emulator_total]./var_ensemble_anova*100);
		a(1).FaceColor=[0.8500 0.2550 0.0980];
		a(2).FaceColor=[0 0.5570 0.7410];
		a(4).FaceColor=a(3).FaceColor;
		a(3).FaceColor=[157 102 255]/255;
		xlim([2016 2100])
		ylim([0 100])
		xlabel('Years','fontsize',18)
		ylabel('% of \sigma^2_{tot}','fontsize',18)
		set(gca,'fontsize',18);
		grid on
		ax = gca;
		ax.Layer = 'top';
		text(2009,0,char(96+iglacier),'fontsize',20,'fontweight','b')
		text(2055,104,glacier_names{iglacier},'fontsize',20,'HorizontalAlignment', 'center')

		error
		var_4(end)./var_ensemble_anova(end)*100
		%eval(['export_fig(''Figures/RelativeVariance_Glacier_' num2str(indexglaciernetcdf) '.pdf'')'])
	end
	if 0,
		eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/var_glacier_' num2str(indexglaciernetcdf) ''])
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		figure(4)
		set(gcf,'color','w')
		set(gcf, 'position',[100 100 650 400])
		plot(years,sqrt(var_ensemble_anova),'k','linewidth',2); hold off
		hold on
		plot(years,sqrt(var_climate),'linewidth',2,'color',[0.8500 0.2550 0.0980])
		plot(years,sqrt(var_ice),'linewidth',2,'color',[0 0.5570 0.7410])
		plot(years,sqrt(var_4),'linewidth',2,'color',[157 102 255]/255)
		plot(years,sqrt(var_emulator_total),'linewidth',2,'color',[0.9290 0.6940 0.1250])
		xlim([2016 2100])
		if iglacier==1,
			ylim([0 20])
		else
			ylim([0 6])
		end
		xlabel('Years','fontsize',18)
		ylabel('\sigma (mm SLE)','fontsize',18)
		set(gca,'fontsize',18);
		grid on
		text(2007,0,char(96+iglacier),'fontsize',18,'fontweight','b')
		if iglacier==1,
			text(2055,20.8,glacier_names{iglacier},'fontsize',20,'HorizontalAlignment', 'center')
		else
			text(2055,6.25,glacier_names{iglacier},'fontsize',20,'HorizontalAlignment', 'center')
		end
		ax = gca;
		ax.Layer = 'top';

		eval(['export_fig(''Figures/UncertaintyAntarctica_Glacier_' num2str(indexglaciernetcdf) '.pdf'')'])
		error
	end

	if 0, %legend
		set(gcf, 'position',[100 100 300 400])
		set(gcf,'color','w')
		eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/var_glacier_' num2str(indexglaciernetcdf) ''])
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		plot(years,sqrt(var_ensemble_anova),'k','linewidth',2); hold off
		hold on
		plot(years,sqrt(var_climate),'linewidth',2,'color',[0.8500 0.2550 0.0980])
		plot(years,sqrt(var_ice),'linewidth',2,'color',[0 0.5570 0.7410])
		plot(years,sqrt(var_4),'linewidth',2,'color',[157 102 255]/255)
		plot(years,sqrt(var_emulator_total),'linewidth',2,'color',[0.9290 0.6940 0.1250])
		xlim([0 1])
		ylim([0 1])
		l=legend('\sigma total','\sigma climate','\sigma ice','\sigma ice-climate','\sigma emulation',...
			'location','East','fontsize',20,'box','off')
		axis off
		h = gcf;
		set(h,'Units','Inches');
      pos = get(h,'Position');
		print(gcf, '-dpdf', '-painters', './Figures/UncertaintyAntarctica_Legend.pdf');
		error
	end

end %}}}

%Figure 9
if perform(org,'VarianceRCP85NoSMBTotalWithCtrl_emulator_anovan')% {{{

	if 0,
		load Data/EmulationData/AntarcticaWithCtrl/ant_dyn_loss_emulated_peter

		icesheetmodels=icesheetmodellong_results(:,1);
		climatemodels=climatemodel_results(:,1);
		var_ensemble=var(-individual_results_limnsw_timeseries,1);
		var_ensemble_anova=zeros(85,1);
		var_climate=zeros(85,1);
		var_ice=zeros(85,1);
		var_distrib=zeros(85,1);
		var_4=zeros(85,1);
		var_5=zeros(85,1);
		var_6=zeros(85,1);
		var_7=zeros(85,1);
		var_error=zeros(85,1);
		num_sample_distrib=100;
		for yeari=1:85,
			yeari
			distribution=repmat(1:num_sample_distrib,210,1); distribution=distribution(:);
			individual_results=-individual_results_limnsw_timeseries(:,yeari);
			individual_std=individual_results_std_timeseries(:,yeari);
			distribution_results=normrnd(repmat(individual_results,num_sample_distrib,1),repmat(individual_std,num_sample_distrib,1));
			[p,tabl]=anovan(distribution_results,{repmat(icesheetmodels,num_sample_distrib,1),repmat(climatemodels,num_sample_distrib,1),distribution},'model',2,'display','off') ;
			var_ice(yeari)=tabl{2,2}/(210*num_sample_distrib); %first line and column are names of variables and parameters
			var_climate(yeari)=tabl{3,2}/(210*num_sample_distrib);
			var_distrib(yeari)=tabl{4,2}/(210*num_sample_distrib);
			var_4(yeari)=tabl{5,2}/(210*num_sample_distrib);
			var_5(yeari)=tabl{6,2}/(210*num_sample_distrib);
			var_6(yeari)=tabl{7,2}/(210*num_sample_distrib);
			var_error(yeari)=tabl{8,2}/(210*num_sample_distrib);
			var_ensemble_anova(yeari)=tabl{end,2}/(210*num_sample_distrib);
		end
		save Data/EmulationData/AntarcticaWithCtrl/var_ant_data_peter var_ensemble_anova var_ice var_climate var_distrib var_4 var_5 var_6 var_error
	end

	years=2016:2100;
	if 0,
		load Data/EmulationData/AntarcticaWithCtrl/var_ant_data
		figure(3);
		set(gcf,'color','w')
		set(gcf, 'position',[100 100 650 400])
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		a=area(years,([var_climate,var_ice,var_4,var_emulator_total]./var_ensemble_anova*100));
		a(1).FaceColor=[0.8500 0.2550 0.0980];
		a(2).FaceColor=[0 0.5570 0.7410]
		a(4).FaceColor=a(3).FaceColor;
		a(3).FaceColor=[157 102 255]/255;
		legend('climate','ice','ice-climate','emulation','location','East','fontsize',16)
		xlim([2016 2100])
		ylim([0 100])
		xlabel('Years','fontsize',16)
		ylabel('% of \sigma^2_{tot}','fontsize',16)
		set(gca,'fontsize',16);
		text(2008,0,'b','fontsize',18,'fontweight','b');
		grid on
		ax = gca;
		ax.Layer = 'top';
		yticks(0:10:100)
		%proportion of error
		error
		export_fig('Figures/RelativeVarianceAntarcticaWithCtrl.pdf')
		sigmatomean=sqrt(var_ensemble_anova)./(-mean(individual_results_limnsw_timeseries,1)')*100;

	end
	if 1,
		load Data/EmulationData/AntarcticaWithCtrl/var_ant_data
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		set(gcf,'color','w')
		set(gcf, 'position',[100 100 650 400])
		plot(years,sqrt(var_ensemble_anova),'k','linewidth',2); hold off
		hold on
		plot(years,sqrt(var_climate),'linewidth',2,'color',[0.8500 0.2550 0.0980])
		plot(years,sqrt(var_ice),'linewidth',2,'color',[0 0.5570 0.7410])
		plot(years,sqrt(var_4),'linewidth',2,'color',[157 102 255]/255)
		plot(years,sqrt(var_emulator_total),'linewidth',2,'color',[0.9290 0.6940 0.1250])
		legend('\sigma total','\sigma climate','\sigma ice','\sigma ice-climate','\sigma emulation','location','NorthWest','fontsize',16,'box','off')
		xlim([2016 2100])
		text(2010,0,'a','fontsize',18,'fontweight','b');
		%ylim([0 50])
		xlabel('Years','fontsize',16)
		ylabel('\sigma (mm SLE)','fontsize',16)
		set(gca,'fontsize',16);
		grid on
		ax = gca;
		ax.Layer = 'top';
		error

		export_fig('Figures/UncertaintyAntarcticaWithCtrl.pdf')
	end

end %}}}

