steps=[12];

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=['./' group '/' simul '/' expename '/computed_ivol_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabf_file=['./' 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=['./' group '/' simul '/' expename '/computed_ivol_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabf_file=['./' 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, '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')
		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)');
		xlim([0.5 3.5]);
		ylim([-40 160]);
		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
		legend(results_model,'location','EastOutside')
		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','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=['./' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabf_file=['./' 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=['./' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
					exptendacabf_file=['./' 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, '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')

		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')
		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={};
		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=['./' group '/' simul '/' expename '/computed_bmbfl_AIS_' group '_' simul '_' expename '.nc'];
					expbmbfl_file=['./' 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=['./' group '/' simul '/' expename '/computed_bmbfl_AIS_' group '_' simul '_' expename '.nc'];
					expbmbfl_file=['./' 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')

end %}}}

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

	
	clf;
	figure(1);  set(gcf,'position',[440   377   730   420])
	load Data/Antarctica/test_emulated.mat
	load Data/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)
	legend('AWI\_PISM1\_o + NorESM','NCAR\_CISM\_o + UKESM','ULB\_fETISh32\_o + HadGEM2','ULB\_fETISh16\_o + CSIRO',...
		'AWI\_PISM1\_s + CESM2','ILTS\_PIK\_SICOPOLIS\_s + MIROC','IMAU\_IMAUICE2\_s + IPSL','JPL1\_ISSM\_s + MIROC',...
		'LSCE\_GRISLI\_s + CCSM4','NCAR\_CISM\_s + CCMS4','UCIJPL\_ISSM\_s + NorESM2','ULB\_fETISh32\_s + IPSL2',...
		'ULB\_fETHISh16\_s + MIROC','VUB\_AISMPALEO\_s + MIROC','location','NorthWest','box','off')

	export_fig('Figures/SimulEmulTestAntarctica.pdf')

end %}}}

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

	load Data/Antarctica/ant_dyn_loss_emulated

	time=[2015:2100];
	figure(1); clear figure; set(gcf,'color','w'); hold on;
	plot(-1,-1,'color',[128 128 128]/256,'linewidth',3.0)
	plot(-1,-1,'color',[178 34 34]/256,'linewidth',3.0)
	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)
	legend('Simulated','Emulated','location','NorthWest','fontsize',18,'box','off')
	export_fig('Figures/Emulation_Antarctica.pdf')

end %}}}

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

	colors = distinguishable_colors(length(model_list)+10);
	colors(19,:)=colors(end,:);
	colors(4,:)=[0,0,0];
	colors(6,:)=[240, 228, 66]/256;
	colors(8,:)=[86, 180, 233]/256;
	colors(3,:)=[0, 158, 115]/256;
	colors(1,:)=[0, 114, 178]/256;
	colors(2,:)=[213, 94, 0]/256;
	colors(11,:)=[204, 121, 167]/256;
	colors(5,:)=[256, 0, 185]/256;

	all_experiments_list={'exp05','exp01',... %open and standard NorESM
		'exp06','exp02',...
		'exp08','exp04',...
		'expA5','expA1',... 
		'expA6','expA2',... 
		'expA7','expA3',... 
		'expB6','expB1',... 
		'expB8','expB3',...
		'expB9','expB4',...
		'expB10','expB5',...
		};

	climate_model_list={'NorESM','MIROC','CCSM4','HadGEM2','CSIRO','IPSL','CNRM','UKESM','CESM2'};

	figure(2); clf;
	for isubplot=1:12,

		figure(2);
		if isubplot==1
			set(gcf,'color','w')
			set(gcf, 'position',[69 1 809 796])
		end

		if isubplot<10
			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:2;
						expename=all_experiments_list{2*(isubplot-1)+iexp};
						specifics;

						eval(['isexp=is' expename ';'])
						if isexp,
							number_models=number_models+1;
							explimnsw_file=['./' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
							exptendacabfgr_file=['./' 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
			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;
			subplotmodel(4,3,isubplot,pairoptions('tightsubplot',1,'gap',[.04 .03],'hmargin',[0.1 0.03],'vmargin',[.1 .01]));
			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])
			if isubplot>6
				xlabel('Years','fontsize',16)
			else
				xticklabels('');
			end
			if isubplot==1 | isubplot==4 | isubplot==7
				ylabel('SLR (mm)','fontsize',16)
				set(gca, 'YTick', [0:20:100]);
			else
				yticklabels('');
			end
			text(2057,92,climate_model_list{isubplot},'fontsize',16,'HorizontalAlignment', 'center')
			if isubplot==8 | isubplot==9,
				text(2008,-2,char(96+isubplot),'fontsize',20,'fontweight','b');
			else
				text(2008,-10,char(96+isubplot),'fontsize',20,'fontweight','b');
			end
			ax = gca;
			ax.FontSize = 16;
			grid on
		end

		if isubplot==11,
			subplotmodel(4,3,isubplot,pairoptions('tightsubplot',1,'gap',[.04 .03],'hmargin',[0.01 0.01],'vmargin',[.05 .01]));
			%subplotmodel(4,3,isubplot,pairoptions('tightsubplot',1,'gap',[.04 .03],'hmargin',[0.1 0.03],'vmargin',[.1 .01]));
			colorm = distinguishable_colors(186);
			colorm(186,:)=colorm(7,:);
			colorm(92,:)=colorm(9,:);
			colorm(82,:)=colorm(10,:);
			colorm(13,:)=colorm(4,:);
			colorm(20,:)=[240, 228, 66]/256;
			colorm(172,:)=[86, 180, 233]/256;
			colorm(83,:)=[0, 158, 115]/256;
			colorm(7,:)=[0, 114, 178]/256;
			colorm(2,:)=[213, 94, 0]/256;
			colorm(25,:)=[204, 121, 167]/256;
			colorm(18,:)=[256  0 185]/256;

			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;

			[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]);
		end

		if isubplot==12,
			subplot(4,3,isubplot)
			patch([0 0 0],[0 0 0],colors(1,:));
			patch([0 0 0],[0 0 0],colors(2,:));
			patch([0 0 0],[0 0 0],colors(3,:));
			patch([0 0 0],[0 0 0],colors(4,:));
			patch([0 0 0],[0 0 0],colors(5,:));
			patch([0 0 0],[0 0 0],colors(6,:));
			patch([0 0 0],[0 0 0],colors(7,:));
			patch([0 0 0],[0 0 0],colors(8,:));
			patch([0 0 0],[0 0 0],colors(9,:));
			patch([0 0 0],[0 0 0],colors(10,:));
			patch([0 0 0],[0 0 0],colors(11,:));
			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',14)
			%l.Position=[0 0 1 1];
		end
	end

	export_fig('Figures/GlaciersAllGCMs.pdf')

end %}}}

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

	colorm = distinguishable_colors(186);
	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,:)=[0,0,0];
	colorm(2,:)=colorm(2,:);
	colorm(7,:)=colorm(1,:);

	colorm(20,:)=[240, 228, 66]/256;
	colorm(172,:)=[86, 180, 233]/256;
	colorm(83,:)=[0, 158, 115]/256;
	colorm(7,:)=[0, 114, 178]/256;
	colorm(2,:)=[213, 94, 0]/256;
	colorm(25,:)=[204, 121, 167]/256;
	colorm(18,:)=[256, 0, 185]/256;

	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=['./' group '/' simul '/' expename '/computed_ivaf_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
				exptendacabfgr_file=['./' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
				exptendlibmassbffl_file=['./' 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
	colors=colorm;
	figure(1); set(gcf,'color','w'); set(gcf,'Position',[400 500 800 500]);
	%subplot(1,2,1)
	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');
	legend('NorESM','MIROC','CCSM4','HadGEM2','CSIRO','IPSL','CNRM','UKESM','CESM2','fontsize',16,'location','NorthWest')
	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=colors;
	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');

	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');

	p1=patch([0 0 0],[0 0 0],colorm(7,:));
	p2=patch([0 0 0],[0 0 0],colorm(2,:));
	p3=patch([0 0 0],[0 0 0],colorm(83,:));
	p4=patch([0 0 0],[0 0 0],colorm(13,:));
	p5=patch([0 0 0],[0 0 0],colorm(18,:));
	p6=patch([0 0 0],[0 0 0],colorm(20,:));
	p7=patch([0 0 0],[0 0 0],colorm(186,:));
	p8=patch([0 0 0],[0 0 0],colorm(172,:));
	p9=patch([0 0 0],[0 0 0],colorm(92,:));
	p10=patch([0 0 0],[0 0 0],colorm(82,:));
	p11=patch([0 0 0],[0 0 0],colorm(25,:));
	p12=patch([0 0 0],[0 0 0],'w','FaceColor','w','EdgeColor','w');

	l=legend([p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11],'Thwaites','Pine Island','Totten',...
		'Getz','Kamb','MacAyeal','Institute','Moller','Leopold and Astrid Coast','Moscow','Whillans',...
		'position',[0.6 0.11515 0.33 0.337],'box','off','fontsize',12)
	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/Figure5.pdf');

end %}}}

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

	if 0,
		load Data/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/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 0,
		load Data/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';
		export_fig('Figures/RelativeVarianceAntarctica.pdf')

	end
	if 1,
		load Data/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
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'};
	close all
	for iglacier=1:11
		indexglaciernetcdf=glacier_list(iglacier);

		if 0,
			eval(['load Data/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
			eval(['save Data/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

	years=2016:2100;
	if 1,
		eval(['load Data/Glacier_' num2str(indexglaciernetcdf) '/var_glacier_' num2str(indexglaciernetcdf) ''])
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		if iglacier==1
			set(gcf,'color','w')
			set(gcf, 'position',[100 100 650 400])
		end
		subplotmodel(4,3,iglacier,pairoptions('tightsubplot',1,'gap',[.04 .03],'hmargin',[0.1 0.03],'vmargin',[.1 .01]));
		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
		if iglacier>8,
			xlabel('Years','fontsize',18)
		else
			xticklabels('');
		end
		if iglacier==1 | iglacier==4 | iglacier ==7 | iglacier==10,
			ylabel('\sigma (mm SLE)','fontsize',18)
		end
		set(gca,'fontsize',18);
		grid on
		if iglacier==1,
			text(2055,18.2,glacier_names{iglacier},'fontsize',20,'HorizontalAlignment', 'center')
			text(2020,16,char(96+iglacier),'fontsize',18,'fontweight','b')
		else
			text(2055,5.50,glacier_names{iglacier},'fontsize',20,'HorizontalAlignment', 'center')
			text(2020,4.7,char(96+iglacier),'fontsize',18,'fontweight','b')
		end
		ax = gca;
		ax.Layer = 'top';

		if iglacier==11, %legend
			subplotmodel(4,3,12,pairoptions('tightsubplot',1,'gap',[.04 .03],'hmargin',[0.1 0.01],'vmargin',[.1 .01]));
			eval(['load Data/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','South','fontsize',20,'box','off')
			axis off
			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/Figure7.pdf');
		end
	end

end

end %}}}

%Figure 8
if perform(org,'RelativeUncertaintyRCP85NoSMBTotal_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'};
	close all
	for iglacier=1:11
		indexglaciernetcdf=glacier_list(iglacier);

	years=2016:2100;
	if 1,
		if iglacier==1
			set(gcf,'color','w')
			set(gcf, 'position',[100 100 650 400])
		end
		eval(['load Data/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);
		if iglacier==1
			set(gcf,'color','w')
			set(gcf, 'position',[100 100 650 400])
		end
		var_emulator_total=var_distrib+var_5+var_6+var_error;
		subplotmodel(4,3,iglacier,pairoptions('tightsubplot',1,'gap',[.04 .03],'hmargin',[0.1 0.03],'vmargin',[.1 .01]));
		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]);
		if iglacier>8;
			xlabel('Years','fontsize',18); 
		else
			xticklabels('');
		end
		if rem(iglacier,3)==1;
			ylabel('% of \sigma^2_{tot}','fontsize',18);
		else
			yticklabels('');
		end
		set(gca,'fontsize',18);
		grid on
		ax = gca;
		ax.Layer = 'top';
		text(2020,10,char(96+iglacier),'fontsize',20,'fontweight','b')
		text(2055,100,glacier_names{iglacier},'fontsize',20,'HorizontalAlignment', 'center','VerticalAlignment','top')

		var_4(end)./var_ensemble_anova(end)*100
		if iglacier==11, %legend
			subplotmodel(4,3,12,pairoptions('tightsubplot',1,'gap',[.04 .03],'hmargin',[0.1 0.01],'vmargin',[.1 .01]));
			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,'location','South')
			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/Figure8.pdf');
		end
	end

end

end %}}}

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

	if 0,
		load Data/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/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/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)

		export_fig('Figures/RelativeVarianceAntarcticaWithCtrl.pdf')
		sigmatomean=sqrt(var_ensemble_anova)./(-mean(individual_results_limnsw_timeseries,1)')*100;

	end
	if 1,
		load Data/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');
		xlabel('Years','fontsize',16)
		ylabel('\sigma (mm SLE)','fontsize',16)
		set(gca,'fontsize',16);
		grid on
		ax = gca;
		ax.Layer = 'top';

		export_fig('Figures/UncertaintyAntarcticaWithCtrl.pdf')
	end

end %}}}

%Figure Summary
if perform(org,'Plot_SummaryUncertainty')% {{{

	glacier_names={'Thwaites','Pine Island','Totten','Getz','Kamb','MacAyeal','Institute',...
		'Moller','Leopold and Astrid Coast','Moscow','Whillans'};
	glacier_uncertainty = [18.7 3.7 3.5 1.8 1.9 1.6 5.7 2.7 1.7 2.8 1.3];
	glacier_climate = [4 5 25 11 37 49 29 32 9 18 53];
	glacier_ice = [49 56 36 37 23 14 26 24 49 39 14];
	glacier_iceclimate = [39 22 28 26 28 26 36 33 28 29 20];
	glacier_emulator = [8 17 11 26 12 11 9 11 14 14 13];

	set(gcf,'color','w')
	h = barh([glacier_climate; glacier_ice; glacier_iceclimate; glacier_emulator]'.*glacier_uncertainty'/100,'stacked')
	h(1).FaceColor=[0.8500 0.2550 0.0980];
	h(2).FaceColor=[0 0.5570 0.7410];
	h(4).FaceColor=h(3).FaceColor;
	h(3).FaceColor=[157 102 255]/255;
	set(gca, 'YTick', 1:size(glacier_names,2) );
	set(gca, 'YTicklabels', glacier_names,'fontsize',14);
	grid on
	l= legend('climate','ice','ice-climate','emulator','fontsize',16,'location','East','box','off')
	title(l, 'Sources of uncertainty')

	title('2100 sea level uncertainty from glaciers contributing most to increased sea level','fontsize',16)
   xlabel('Uncertainty in dynamic contribution to sea level (mm)')
	export_fig('./Figures/SummaryFigure.pdf')

end %}}}

