steps=[10];

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

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_list={'ILTS_PIK_SICOPOLIS','JPL1_ISSM','LSCE_GRISLI','NCAR_CISM','UCIJPL_ISSM'};

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

%'exp01','exp02','exp03','exp04','exp05','exp06','exp07','exp08','exp09','exp10','exp11','exp12','exp13',...
experiments_list={...
	'hist','ctrl_proj','ctrl_proj_open',...
	};

clustername='thwaites';
%Cluster parameters{{{
if strcmpi(clustername,'pfe');
	cluster=pfe('grouplist','s1690','numnodes',1,'cpuspernode',28,'time',2*60,'interactive',interactive,'processor','bro','queue','devel'); %time in minutes
elseif strcmpi(clustername,'thwaites'),
	cluster=generic('name','thwaites','np',20);
elseif strcmpi(clustername,'murdo'),
	cluster=generic('name','murdo','np',8);
else
	error('cluster not supported yet');
end%}}}
org=organizer('repository','./Models','prefix','Emulation.','steps',steps);
%addpath('/u/astrid-r1b/seroussi/issmjpl/proj-seroussi/CODE');

%Emulation of missing results Antarctica
if perform(org,'EmulatorYearByYear')% {{{

	%Figure out the grp
	load Data/EmulationData/Antarctica/training_ant.mat
	load Data/EmulationData/Antarctica/validation_ant.mat

	training_meltparam_binary=zeros(size(training_meltparam));
	training_meltparam_binary(1:30*85)=1;
	validation_meltparam_binary=zeros(size(validation_meltparam));
	validation_meltparam_binary(1:10*85)=1;
	std_global_training=zeros(size(training_meltparam));
	pred_global_training=zeros(size(training_meltparam));
	ycit_global_training=zeros(size(training_meltparam));
	pred_global_validation=zeros(size(validation_meltparam));
	ycit_global_validation=zeros(size(validation_meltparam));
	nugget_chosen=zeros(85,1);
	
	indexglaciernetcdf=0
	updateclimatenumber;

	for icase=1:length(training_climate_number),
		if strcmp(training_climate{icase},'NorESM'),
			training_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(training_climate{icase},'HadGEM2'),
			training_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(training_climate{icase},'CSIRO'),
			training_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(training_climate{icase},'MIROC'),
			training_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(training_climate{icase},'IPSL'),
			training_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(training_climate{icase},'CCSM4'),
			training_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(training_climate{icase},'CESM2'),
			training_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(training_climate{icase},'UKESM'),
			training_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(training_climate{icase},'CNRM-ESM'),
			training_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(training_climate{icase},'CNRM-CM6'),
			training_climate_number(icase)=meanoceantemp(10);
		else 
			training_climate{icase}
			error('model not recognized');
		end
	end

	for icase=1:length(validation_climate_number),
		if strcmp(validation_climate{icase},'NorESM'),
			validation_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(validation_climate{icase},'HadGEM2'),
			validation_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(validation_climate{icase},'CSIRO'),
			validation_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(validation_climate{icase},'MIROC'),
			validation_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(validation_climate{icase},'IPSL'),
			validation_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(validation_climate{icase},'CCSM4'),
			validation_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(validation_climate{icase},'CESM2'),
			validation_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(validation_climate{icase},'UKESM'),
			validation_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(validation_climate{icase},'CNRM-ESM'),
			validation_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(validation_climate{icase},'CNRM-CM6'),
			validation_climate_number(icase)=meanoceantemp(10);
		else 
			validation_climate{icase}
			error('model not recognized');
		end
	end

	for i=1:85,
		%i
		tab_training=table(training_ice(i:85:end),training_climate_number(i:85:end),training_meltparam_binary(i:85:end),training_values(i:85:end),'VariableNames', {'ice', 'climate','meltparam','values'});
		tab_validation=table(validation_ice(i:85:end),validation_climate_number(i:85:end),validation_meltparam_binary(i:85:end),'VariableNames', {'ice', 'climate','meltparam'});
		gpr2 = fitrgp(tab_training,'values','KernelFunction','matern32');
		eval(['gpr_' num2str(i) '=fitrgp(tab_training,''values'',''KernelFunction'',''matern32'');'])
		nugget_chosen2(i)=gpr2.Sigma;
		[ytraining,ystdtraining,ycit] = predict(gpr2,tab_training);
		[yvalidation,~,ycitvalidation] = predict(gpr2,tab_validation);
		std_global_training(i:85:end)=ystdtraining;
		pred_global_training(i:85:end)=ytraining;
		ycit_global_training(i:85:end)=ycit(:,2)-ycit(:,1);
		pred_global_validation(i:85:end)=yvalidation;
		ycit_global_validation(i:85:end)=ycitvalidation(:,2)-ycitvalidation(:,1);
		eval(['save Data/EmulationData/Antarctica/gpr_final_year_' num2str(i) ' gpr_' num2str(i) ''])

	end

	save Data/EmulationData/Antarctica/EmulationResultsYearlyEmulation pred_global_training ycit_global_training pred_global_validation ycit_global_validation

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

	experiments_list={'exp01','exp05',...
		'exp02','exp06',...
		'exp04','exp08',...
		'expA1','expA5',...
		'expA2','expA6',...
		'expA3','expA7',...
		'expB1','expB6',...
		'expB3','expB8',...
		'expB4','expB9',...
		'expB5','expB10'};

	experiments_table={'exp01','exp05';...
		'exp02','exp06';...
		'exp04','exp08';...
		'expA1','expA5';...
		'expA2','expA6';...
		'expA3','expA7';...
		'expB1','expB6';...
		'expB3','expB8';...
		'expB4','expB9';...
		'expB5','expB10'};
	
	climatemodel_list={'NorESM','NorESM',...
		'MIROC','MIROC',...
		'CCSM4','CCSM4',...
		'HadGEM2','HadGEM2',...
		'CSIRO','CSIRO',...
		'IPSL','IPSL',...
		'CNRM-ESM','CNRM-ESM',...
		'UKESM','UKESM',...
		'CESM2','CESM2',...
		'CNRM-CM6','CNRM-CM6'};

	climatemodel_results={};
	icesheetmodel_results={};
	icesheetmodellong_results={};
	meltparam_results={};
	meltparam_results_binary=[];
	emulated_results=[];
	individual_results_limnsw_timeseries=[];
	individual_results_std_timeseries=[];
	icesheetmodel_longlist={};
	
	icesheetmodel_longlist{end+1}='AWI_PISM1_open';
	nummodel=1; %initialize at 1 because AWI put by default
	existing_results=zeros(21,10);

	%Load data and gprs
   load Data/EmulationData/Antarctica/training_ant.mat
   load Data/EmulationData/Antarctica/validation_ant.mat
   load Data/EmulationData/Antarctica/test_ant.mat

   for i=1:85,
      eval(['load Data/EmulationData/Antarctica/gpr_final_year_' num2str(i) ''])
   end

   training_meltparam_binary=zeros(size(training_meltparam));
   training_meltparam_binary(1:30*85)=1;
   validation_meltparam_binary=zeros(size(validation_meltparam));
   validation_meltparam_binary(1:10*85)=1;
   test_meltparam_binary=zeros(size(test_meltparam));
   test_meltparam_binary(1:4*85)=1;

	climatecase=0;
	indexglaciernetcdf=0
	updateclimatenumber;

	for icase=1:length(training_climate_number),
		if strcmp(training_climate{icase},'NorESM'),
			training_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(training_climate{icase},'HadGEM2'),
			training_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(training_climate{icase},'CSIRO'),
			training_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(training_climate{icase},'MIROC'),
			training_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(training_climate{icase},'IPSL'),
			training_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(training_climate{icase},'CCSM4'),
			training_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(training_climate{icase},'CESM2'),
			training_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(training_climate{icase},'UKESM'),
			training_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(training_climate{icase},'CNRM-ESM'),
			training_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(training_climate{icase},'CNRM-CM6'),
			training_climate_number(icase)=meanoceantemp(10);
		else 
			training_climate{icase}
			error('model not recognized');
		end
	end

	for icase=1:length(validation_climate_number),
		if strcmp(validation_climate{icase},'NorESM'),
			validation_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(validation_climate{icase},'HadGEM2'),
			validation_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(validation_climate{icase},'CSIRO'),
			validation_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(validation_climate{icase},'MIROC'),
			validation_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(validation_climate{icase},'IPSL'),
			validation_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(validation_climate{icase},'CCSM4'),
			validation_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(validation_climate{icase},'CESM2'),
			validation_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(validation_climate{icase},'UKESM'),
			validation_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(validation_climate{icase},'CNRM-ESM'),
			validation_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(validation_climate{icase},'CNRM-CM6'),
			validation_climate_number(icase)=meanoceantemp(10);
		else 
			validation_climate{icase}
			error('model not recognized');
		end
	end
	
	for icase=1:length(test_climate_number),
		if strcmp(test_climate{icase},'NorESM'),
			test_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(test_climate{icase},'HadGEM2'),
			test_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(test_climate{icase},'CSIRO'),
			test_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(test_climate{icase},'MIROC'),
			test_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(test_climate{icase},'IPSL'),
			test_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(test_climate{icase},'CCSM4'),
			test_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(test_climate{icase},'CESM2'),
			test_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(test_climate{icase},'UKESM'),
			test_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(test_climate{icase},'CNRM-ESM'),
			test_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(test_climate{icase},'CNRM-CM6'),
			test_climate_number(icase)=meanoceantemp(10);
		else 
			test_climate{icase}
			error('model not recognized');
		end
	end

	error
	for imodel=1:length(model_list),
		modelname=model_list{imodel};
		for iexp=1:length(experiments_list)/2,
			expename=experiments_list{2*iexp-1};
			specifics;

			eval(['isexp=is' expename ';'])
			if isexp,
				if strcmp(icesheetmodel_longlist{end},['' modelname '_open'])==0,
					icesheetmodel_longlist{end+1}=['' modelname '_open'];
					nummodel=nummodel+1;
				end
				existing_results(nummodel,iexp)=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'];
				eval(['results_limnsw=ncread(explimnsw_file,''ivaf'') *ice_density/(10^9*1000);']) %from m^3 to Gt
				eval(['results_tendacabfgr=ncread(exptendacabfgr_file,''smbgr'') *yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
				individual_results_limnsw_timeseries(end+1,:)=(results_limnsw(1:85)-cumsum(results_tendacabfgr(1:85)))/362.5;
				if strcmpi(modelname,'UCIJPL_ISSM') & strcmpi(expename,'exp01'),
					individual_results_limnsw_timeseries(end,44)=(individual_results_limnsw_timeseries(end,43)+individual_results_limnsw_timeseries(end,45))/2;
				end
				individual_results_std_timeseries(end+1,:)=zeros(85,1);
				icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
				icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_open']},1,85)];
				climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp-1}},1,85)];
				emulated_results(end+1,:)=0*ones(1,85);
				meltparam_results=[meltparam_results;repmat({'open'},85,1)];
				meltparam_results_binary(end+1,:)=1*ones(1,85);
			else %no existing results so emulate them
				if isexp01==1,
					%there is some data for this model so you can emulate the missing simulations
					meltparam_tosimulate=repmat({'open'},85,1);
					meltparam_tosimulate_binary=ones(85,1);
					icesheetmodel_tosimulate=repmat({modelname},85,1);
					if strcmp(climatemodel_list{2*iexp-1},'NorESM'),
						climate_number=meanoceantemp(1);
					elseif strcmp(climatemodel_list{2*iexp-1},'HadGEM2'),
						climate_number=meanoceantemp(2);
					elseif strcmp(climatemodel_list{2*iexp-1},'CSIRO'),
						climate_number=meanoceantemp(3);
					elseif strcmp(climatemodel_list{2*iexp-1},'MIROC'),
						climate_number=meanoceantemp(4);
					elseif strcmp(climatemodel_list{2*iexp-1},'IPSL'),
						climate_number=meanoceantemp(5);
					elseif strcmp(climatemodel_list{2*iexp-1},'CCSM4'),
						climate_number=meanoceantemp(6);
					elseif strcmp(climatemodel_list{2*iexp-1},'CESM2'),
						climate_number=meanoceantemp(7);
					elseif strcmp(climatemodel_list{2*iexp-1},'UKESM'),
						climate_number=meanoceantemp(8);
					elseif strcmp(climatemodel_list{2*iexp-1},'CNRM-ESM'),
						climate_number=meanoceantemp(9);
					elseif strcmp(climatemodel_list{2*iexp-1},'CNRM-CM6'),
						climate_number=meanoceantemp(10);
					else 
						climatemodel_list{2*iexp-1}
						error('model not recognized');
					end
					climatemodel_tosimulate=repmat({climatemodel_list{2*iexp-1}},85,1);
					climatemodel_tosimulate_number=climate_number*ones(85,1);
					time_tosimulate=1:85;
					yemul=[];
					ystd_emul=[];
					for timesim=1:85,
						table_indiv=table(icesheetmodel_tosimulate(timesim),climatemodel_tosimulate_number(timesim),meltparam_tosimulate_binary(timesim),'VariableNames', {'ice', 'climate','meltparam'});
						eval(['[yemul(timesim),ystd_emul(timesim),yint_emul]= predict(gpr_' num2str(timesim) ',table_indiv);']);
					end

					%Save with the overall results
					emulated_results(end+1,:)=ones(1,85);
					individual_results_limnsw_timeseries(end+1,:)=yemul;
					individual_results_std_timeseries(end+1,:)=ystd_emul;
					icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
					icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_open']},1,85)];
					climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp-1}},1,85)];
					meltparam_results=[meltparam_results;repmat({'open'},85,1)];
					meltparam_results_binary(end+1,:)=1*ones(1,85);
				else %there is no data for this standard model so don't create one from nothing
					%do nothing
				end
			end

		end %end of model
	end %end of isexp
	for imodel=1:length(model_list),
		modelname=model_list{imodel};
		for iexp=1:length(experiments_list)/2,
			expename=experiments_list{2*iexp};
			specifics;

			eval(['isexp=is' expename ';'])
			if isexp,
				if strcmp(icesheetmodel_longlist{end},['' modelname '_standard'])==0,
					icesheetmodel_longlist{end+1}=['' modelname '_standard'];
					nummodel=nummodel+1;
				end
				existing_results(nummodel,iexp)=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'];
				eval(['results_limnsw=ncread(explimnsw_file,''ivaf'') *ice_density/(10^9*1000);']) %from m^3 to Gt
				eval(['results_tendacabfgr=ncread(exptendacabfgr_file,''smbgr'') *yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
				individual_results_limnsw_timeseries(end+1,:)=(results_limnsw(1:85)-cumsum(results_tendacabfgr(1:85)))/362.5; 
				individual_results_std_timeseries(end+1,:)=zeros(85,1);
				icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
				icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_standard']},1,85)];
				climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp}},1,85)];
				emulated_results(end+1,:)=0*ones(1,85);
				meltparam_results=[meltparam_results;repmat({'standard'},85,1)];
				meltparam_results_binary(end+1,:)=0*ones(1,85);
			else %no existing results so emulate them
				if isexp05==1,
					%there is some data for this model so you can emulate the missing simulations
					meltparam_tosimulate=repmat({'standard'},85,1);
					meltparam_tosimulate_binary=zeros(85,1);
					icesheetmodel_tosimulate=repmat({modelname},85,1);
					if strcmp(climatemodel_list{2*iexp},'NorESM'),
						climate_number=meanoceantemp(1);
					elseif strcmp(climatemodel_list{2*iexp},'HadGEM2'),
						climate_number=meanoceantemp(2);
					elseif strcmp(climatemodel_list{2*iexp},'CSIRO'),
						climate_number=meanoceantemp(3);
					elseif strcmp(climatemodel_list{2*iexp},'MIROC'),
						climate_number=meanoceantemp(4);
					elseif strcmp(climatemodel_list{2*iexp},'IPSL'),
						climate_number=meanoceantemp(5);
					elseif strcmp(climatemodel_list{2*iexp},'CCSM4'),
						climate_number=meanoceantemp(6);
					elseif strcmp(climatemodel_list{2*iexp},'CESM2'),
						climate_number=meanoceantemp(7);
					elseif strcmp(climatemodel_list{2*iexp},'UKESM'),
						climate_number=meanoceantemp(8);
					elseif strcmp(climatemodel_list{2*iexp},'CNRM-ESM'),
						climate_number=meanoceantemp(9);
					elseif strcmp(climatemodel_list{2*iexp},'CNRM-CM6'),
						climate_number=meanoceantemp(10);
					else 
						climatemodel_list{2*iexp}
						error('model not recognized');
					end
					climatemodel_tosimulate=repmat({climatemodel_list{2*iexp}},85,1);
					climatemodel_tosimulate_number=climate_number*ones(85,1);
					time_tosimulate=1:85;
					yemul=[];
					ystd_emul=[];
					for timesim=1:85,
						table_indiv=table(icesheetmodel_tosimulate(timesim),climatemodel_tosimulate_number(timesim),meltparam_tosimulate_binary(timesim),'VariableNames', {'ice', 'climate','meltparam'});
						eval(['[yemul(timesim),ystd_emul(timesim),yint_emul]= predict(gpr_' num2str(timesim) ',table_indiv);']);
					end

					%Save with the overall results
					emulated_results(end+1,:)=1*ones(1,85);
					individual_results_limnsw_timeseries(end+1,:)=yemul;
					individual_results_std_timeseries(end+1,:)=ystd_emul;
					icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
					icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_standard']},1,85)];
					climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp}},1,85)];
					meltparam_results=[meltparam_results;repmat({'standard'},85,1)];
					meltparam_results_binary(end+1,:)=0*ones(1,85);
				else %there is no data for this standard model so don't create one from nothing
					%do nothing
				end
			end

		end %end of model
	end %end of isexp

	save Data/EmulationData/Antarctica/ant_dyn_loss_emulated  individual_results_limnsw_timeseries individual_results_std_timeseries icesheetmodel_results icesheetmodellong_results climatemodel_results emulated_results meltparam_results_binary meltparam_results 

end %}}}

%Emulation of missing results Antarctica with ctrl
if perform(org,'WithCtrl_EmulatorYearByYear')% {{{

	%Figure out the grp
	load Data/EmulationData/Antarctica/training_ant_withctrl.mat
	load Data/EmulationData/Antarctica/validation_ant_withctrl.mat

	training_meltparam_binary=zeros(size(training_meltparam));
	training_meltparam_binary(1:30*85)=1;
	validation_meltparam_binary=zeros(size(validation_meltparam));
	validation_meltparam_binary(1:10*85)=1;
	std_global_training=zeros(size(training_meltparam));
	pred_global_training=zeros(size(training_meltparam));
	ycit_global_training=zeros(size(training_meltparam));
	pred_global_validation=zeros(size(validation_meltparam));
	ycit_global_validation=zeros(size(validation_meltparam));
	nugget_chosen=zeros(85,1);
	
	indexglaciernetcdf=0
	updateclimatenumber;

	for icase=1:length(training_climate_number),
		if strcmp(training_climate{icase},'NorESM'),
			training_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(training_climate{icase},'HadGEM2'),
			training_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(training_climate{icase},'CSIRO'),
			training_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(training_climate{icase},'MIROC'),
			training_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(training_climate{icase},'IPSL'),
			training_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(training_climate{icase},'CCSM4'),
			training_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(training_climate{icase},'CESM2'),
			training_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(training_climate{icase},'UKESM'),
			training_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(training_climate{icase},'CNRM-ESM'),
			training_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(training_climate{icase},'CNRM-CM6'),
			training_climate_number(icase)=meanoceantemp(10);
		else 
			training_climate{icase}
			error('model not recognized');
		end
	end

	for icase=1:length(validation_climate_number),
		if strcmp(validation_climate{icase},'NorESM'),
			validation_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(validation_climate{icase},'HadGEM2'),
			validation_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(validation_climate{icase},'CSIRO'),
			validation_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(validation_climate{icase},'MIROC'),
			validation_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(validation_climate{icase},'IPSL'),
			validation_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(validation_climate{icase},'CCSM4'),
			validation_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(validation_climate{icase},'CESM2'),
			validation_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(validation_climate{icase},'UKESM'),
			validation_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(validation_climate{icase},'CNRM-ESM'),
			validation_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(validation_climate{icase},'CNRM-CM6'),
			validation_climate_number(icase)=meanoceantemp(10);
		else 
			validation_climate{icase}
			error('model not recognized');
		end
	end

	for i=1:85,
		tab_training=table(training_ice(i:85:end),training_climate_number(i:85:end),training_meltparam_binary(i:85:end),training_values(i:85:end),'VariableNames', {'ice', 'climate','meltparam','values'});
		tab_validation=table(validation_ice(i:85:end),validation_climate_number(i:85:end),validation_meltparam_binary(i:85:end),'VariableNames', {'ice', 'climate','meltparam'});
		gpr2 = fitrgp(tab_training,'values','KernelFunction','matern32');
		eval(['gpr_' num2str(i) '=fitrgp(tab_training,''values'',''KernelFunction'',''matern32'');'])
		nugget_chosen2(i)=gpr2.Sigma;
		[ytraining,ystdtraining,ycit] = predict(gpr2,tab_training);
		[yvalidation,~,ycitvalidation] = predict(gpr2,tab_validation);
		std_global_training(i:85:end)=ystdtraining;
		pred_global_training(i:85:end)=ytraining;
		ycit_global_training(i:85:end)=ycit(:,2)-ycit(:,1);
		pred_global_validation(i:85:end)=yvalidation;
		ycit_global_validation(i:85:end)=ycitvalidation(:,2)-ycitvalidation(:,1);
		eval(['save Data/EmulationData/AntarcticaWithCtrl/gpr_final_year_' num2str(i) ' gpr_' num2str(i) ''])

	end

	save Data/EmulationData/Antarctica/EmulationResultsYearlyEmulationWithCtrl pred_global_training ycit_global_training pred_global_validation ycit_global_validation

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

	experiments_list={'exp01','exp05',...
		'exp02','exp06',...
		'exp04','exp08',...
		'expA1','expA5',...
		'expA2','expA6',...
		'expA3','expA7',...
		'expB1','expB6',...
		'expB3','expB8',...
		'expB4','expB9',...
		'expB5','expB10'};

	experiments_table={'exp01','exp05';...
		'exp02','exp06';...
		'exp04','exp08';...
		'expA1','expA5';...
		'expA2','expA6';...
		'expA3','expA7';...
		'expB1','expB6';...
		'expB3','expB8';...
		'expB4','expB9';...
		'expB5','expB10'};
	
	climatemodel_list={'NorESM','NorESM',...
		'MIROC','MIROC',...
		'CCSM4','CCSM4',...
		'HadGEM2','HadGEM2',...
		'CSIRO','CSIRO',...
		'IPSL','IPSL',...
		'CNRM-ESM','CNRM-ESM',...
		'UKESM','UKESM',...
		'CESM2','CESM2',...
		'CNRM-CM6','CNRM-CM6'};

	climatemodel_results={};
	icesheetmodel_results={};
	icesheetmodellong_results={};
	meltparam_results={};
	meltparam_results_binary=[];
	emulated_results=[];
	individual_results_limnsw_timeseries=[];
	individual_results_std_timeseries=[];
	icesheetmodel_longlist={};
	
	icesheetmodel_longlist{end+1}='AWI_PISM1_open';
	nummodel=1; %initialize at 1 because AWI put by default
	existing_results=zeros(21,10);

	%Load data and gprs
   load Data/EmulationData/Antarctica/training_ant_withctrl.mat
   load Data/EmulationData/Antarctica/validation_ant_withctrl.mat
   load Data/EmulationData/Antarctica/test_ant_withctrl.mat

   for i=1:85,
      eval(['load Data/EmulationData/AntarcticaWithCtrl/gpr_final_year_' num2str(i) ''])
   end

   training_meltparam_binary=zeros(size(training_meltparam));
   training_meltparam_binary(1:30*85)=1;
   validation_meltparam_binary=zeros(size(validation_meltparam));
   validation_meltparam_binary(1:10*85)=1;
   test_meltparam_binary=zeros(size(test_meltparam));
   test_meltparam_binary(1:4*85)=1;

	climatecase=0;
	indexglaciernetcdf=0
	updateclimatenumber;

	for icase=1:length(training_climate_number),
		if strcmp(training_climate{icase},'NorESM'),
			training_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(training_climate{icase},'HadGEM2'),
			training_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(training_climate{icase},'CSIRO'),
			training_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(training_climate{icase},'MIROC'),
			training_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(training_climate{icase},'IPSL'),
			training_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(training_climate{icase},'CCSM4'),
			training_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(training_climate{icase},'CESM2'),
			training_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(training_climate{icase},'UKESM'),
			training_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(training_climate{icase},'CNRM-ESM'),
			training_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(training_climate{icase},'CNRM-CM6'),
			training_climate_number(icase)=meanoceantemp(10);
		else 
			training_climate{icase}
			error('model not recognized');
		end
	end

	for icase=1:length(validation_climate_number),
		if strcmp(validation_climate{icase},'NorESM'),
			validation_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(validation_climate{icase},'HadGEM2'),
			validation_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(validation_climate{icase},'CSIRO'),
			validation_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(validation_climate{icase},'MIROC'),
			validation_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(validation_climate{icase},'IPSL'),
			validation_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(validation_climate{icase},'CCSM4'),
			validation_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(validation_climate{icase},'CESM2'),
			validation_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(validation_climate{icase},'UKESM'),
			validation_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(validation_climate{icase},'CNRM-ESM'),
			validation_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(validation_climate{icase},'CNRM-CM6'),
			validation_climate_number(icase)=meanoceantemp(10);
		else 
			validation_climate{icase}
			error('model not recognized');
		end
	end
	
	for icase=1:length(test_climate_number),
		if strcmp(test_climate{icase},'NorESM'),
			test_climate_number(icase)=meanoceantemp(1);
		elseif strcmp(test_climate{icase},'HadGEM2'),
			test_climate_number(icase)=meanoceantemp(2);
		elseif strcmp(test_climate{icase},'CSIRO'),
			test_climate_number(icase)=meanoceantemp(3);
		elseif strcmp(test_climate{icase},'MIROC'),
			test_climate_number(icase)=meanoceantemp(4);
		elseif strcmp(test_climate{icase},'IPSL'),
			test_climate_number(icase)=meanoceantemp(5);
		elseif strcmp(test_climate{icase},'CCSM4'),
			test_climate_number(icase)=meanoceantemp(6);
		elseif strcmp(test_climate{icase},'CESM2'),
			test_climate_number(icase)=meanoceantemp(7);
		elseif strcmp(test_climate{icase},'UKESM'),
			test_climate_number(icase)=meanoceantemp(8);
		elseif strcmp(test_climate{icase},'CNRM-ESM'),
			test_climate_number(icase)=meanoceantemp(9);
		elseif strcmp(test_climate{icase},'CNRM-CM6'),
			test_climate_number(icase)=meanoceantemp(10);
		else 
			test_climate{icase}
			error('model not recognized');
		end
	end

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

			eval(['isexp=is' expename ';'])
			if isexp,
				if strcmp(icesheetmodel_longlist{end},['' modelname '_open'])==0,
					icesheetmodel_longlist{end+1}=['' modelname '_open'];
					nummodel=nummodel+1;
				end
				existing_results(nummodel,iexp)=1;
				explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_AIS_' group '_' simul '_' expename '.nc'];
				exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
				eval(['results_limnsw=ncread(explimnsw_file,''ivaf'') *ice_density/(10^9*1000);']) %from m^3 to Gt
				eval(['results_tendacabfgr=ncread(exptendacabfgr_file,''smbgr'') *yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
				individual_results_limnsw_timeseries(end+1,:)=(results_limnsw(1:85)-cumsum(results_tendacabfgr(1:85)))/362.5;
				if strcmpi(modelname,'UCIJPL_ISSM') & strcmpi(expename,'exp01'),
					individual_results_limnsw_timeseries(end,44)=(individual_results_limnsw_timeseries(end,43)+individual_results_limnsw_timeseries(end,45))/2;
				end
				individual_results_std_timeseries(end+1,:)=zeros(85,1);
				icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
				icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_open']},1,85)];
				climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp-1}},1,85)];
				emulated_results(end+1,:)=0*ones(1,85);
				meltparam_results=[meltparam_results;repmat({'open'},85,1)];
				meltparam_results_binary(end+1,:)=1*ones(1,85);
			else %no existing results so emulate them
				if isexp01==1,
					%there is some data for this model so you can emulate the missing simulations
					meltparam_tosimulate=repmat({'open'},85,1);
					meltparam_tosimulate_binary=ones(85,1);
					icesheetmodel_tosimulate=repmat({modelname},85,1);
					if strcmp(climatemodel_list{2*iexp-1},'NorESM'),
						climate_number=meanoceantemp(1);
					elseif strcmp(climatemodel_list{2*iexp-1},'HadGEM2'),
						climate_number=meanoceantemp(2);
					elseif strcmp(climatemodel_list{2*iexp-1},'CSIRO'),
						climate_number=meanoceantemp(3);
					elseif strcmp(climatemodel_list{2*iexp-1},'MIROC'),
						climate_number=meanoceantemp(4);
					elseif strcmp(climatemodel_list{2*iexp-1},'IPSL'),
						climate_number=meanoceantemp(5);
					elseif strcmp(climatemodel_list{2*iexp-1},'CCSM4'),
						climate_number=meanoceantemp(6);
					elseif strcmp(climatemodel_list{2*iexp-1},'CESM2'),
						climate_number=meanoceantemp(7);
					elseif strcmp(climatemodel_list{2*iexp-1},'UKESM'),
						climate_number=meanoceantemp(8);
					elseif strcmp(climatemodel_list{2*iexp-1},'CNRM-ESM'),
						climate_number=meanoceantemp(9);
					elseif strcmp(climatemodel_list{2*iexp-1},'CNRM-CM6'),
						climate_number=meanoceantemp(10);
					else 
						climatemodel_list{2*iexp-1}
						error('model not recognized');
					end
					climatemodel_tosimulate=repmat({climatemodel_list{2*iexp-1}},85,1);
					climatemodel_tosimulate_number=climate_number*ones(85,1);
					time_tosimulate=1:85;
					yemul=[];
					ystd_emul=[];
					for timesim=1:85,
						table_indiv=table(icesheetmodel_tosimulate(timesim),climatemodel_tosimulate_number(timesim),meltparam_tosimulate_binary(timesim),'VariableNames', {'ice', 'climate','meltparam'});
						eval(['[yemul(timesim),ystd_emul(timesim),yint_emul]= predict(gpr_' num2str(timesim) ',table_indiv);']);
					end

					%Save with the overall results
					emulated_results(end+1,:)=ones(1,85);
					individual_results_limnsw_timeseries(end+1,:)=yemul;
					individual_results_std_timeseries(end+1,:)=ystd_emul;
					icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
					icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_open']},1,85)];
					climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp-1}},1,85)];
					meltparam_results=[meltparam_results;repmat({'open'},85,1)];
					meltparam_results_binary(end+1,:)=1*ones(1,85);
				else %there is no data for this standard model so don't create one from nothing
					%do nothing
				end
			end

		end %end of model
	end %end of isexp
	for imodel=1:length(model_list),
		modelname=model_list{imodel};
		for iexp=1:length(experiments_list)/2,
			expename=experiments_list{2*iexp};
			specifics;

			eval(['isexp=is' expename ';'])
			if isexp,
				if strcmp(icesheetmodel_longlist{end},['' modelname '_standard'])==0,
					icesheetmodel_longlist{end+1}=['' modelname '_standard'];
					nummodel=nummodel+1;
				end
				existing_results(nummodel,iexp)=1;
				explimnsw_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_ivaf_AIS_' group '_' simul '_' expename '.nc'];
				exptendacabfgr_file=['ComputedOutputs/' group '/' simul '/' expename '/computed_smbgr_minus_ctrl_proj_AIS_' group '_' simul '_' expename '.nc'];
				eval(['results_limnsw=ncread(explimnsw_file,''ivaf'') *ice_density/(10^9*1000);']) %from m^3 to Gt
				eval(['results_tendacabfgr=ncread(exptendacabfgr_file,''smbgr'') *yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
				individual_results_limnsw_timeseries(end+1,:)=(results_limnsw(1:85)-cumsum(results_tendacabfgr(1:85)))/362.5; 
				individual_results_std_timeseries(end+1,:)=zeros(85,1);
				icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
				icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_standard']},1,85)];
				climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp}},1,85)];
				emulated_results(end+1,:)=0*ones(1,85);
				meltparam_results=[meltparam_results;repmat({'standard'},85,1)];
				meltparam_results_binary(end+1,:)=0*ones(1,85);
			else %no existing results so emulate them
				if isexp05==1,
					%there is some data for this model so you can emulate the missing simulations
					meltparam_tosimulate=repmat({'standard'},85,1);
					meltparam_tosimulate_binary=zeros(85,1);
					icesheetmodel_tosimulate=repmat({modelname},85,1);
					if strcmp(climatemodel_list{2*iexp},'NorESM'),
						climate_number=meanoceantemp(1);
					elseif strcmp(climatemodel_list{2*iexp},'HadGEM2'),
						climate_number=meanoceantemp(2);
					elseif strcmp(climatemodel_list{2*iexp},'CSIRO'),
						climate_number=meanoceantemp(3);
					elseif strcmp(climatemodel_list{2*iexp},'MIROC'),
						climate_number=meanoceantemp(4);
					elseif strcmp(climatemodel_list{2*iexp},'IPSL'),
						climate_number=meanoceantemp(5);
					elseif strcmp(climatemodel_list{2*iexp},'CCSM4'),
						climate_number=meanoceantemp(6);
					elseif strcmp(climatemodel_list{2*iexp},'CESM2'),
						climate_number=meanoceantemp(7);
					elseif strcmp(climatemodel_list{2*iexp},'UKESM'),
						climate_number=meanoceantemp(8);
					elseif strcmp(climatemodel_list{2*iexp},'CNRM-ESM'),
						climate_number=meanoceantemp(9);
					elseif strcmp(climatemodel_list{2*iexp},'CNRM-CM6'),
						climate_number=meanoceantemp(10);
					else 
						climatemodel_list{2*iexp}
						error('model not recognized');
					end
					climatemodel_tosimulate=repmat({climatemodel_list{2*iexp}},85,1);
					climatemodel_tosimulate_number=climate_number*ones(85,1);
					time_tosimulate=1:85;
					yemul=[];
					ystd_emul=[];
					for timesim=1:85,
						table_indiv=table(icesheetmodel_tosimulate(timesim),climatemodel_tosimulate_number(timesim),meltparam_tosimulate_binary(timesim),'VariableNames', {'ice', 'climate','meltparam'});
						eval(['[yemul(timesim),ystd_emul(timesim),yint_emul]= predict(gpr_' num2str(timesim) ',table_indiv);']);
					end

					%Save with the overall results
					emulated_results(end+1,:)=1*ones(1,85);
					individual_results_limnsw_timeseries(end+1,:)=yemul;
					individual_results_std_timeseries(end+1,:)=ystd_emul;
					icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
					icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_standard']},1,85)];
					climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp}},1,85)];
					meltparam_results=[meltparam_results;repmat({'standard'},85,1)];
					meltparam_results_binary(end+1,:)=0*ones(1,85);
				else %there is no data for this standard model so don't create one from nothing
					%do nothing
				end
			end

		end %end of model
	end %end of isexp

	save Data/EmulationData/AntarcticaWithCtrl/ant_dyn_loss_emulated  individual_results_limnsw_timeseries individual_results_std_timeseries icesheetmodel_results icesheetmodellong_results climatemodel_results emulated_results meltparam_results_binary meltparam_results 

end %}}}

%Emulation of missing results by glaciers
if perform(org,'EmulatorYearByYear_glaciers')% {{{

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

	for iglacier=1,

		indexglaciernetcdf=glacier_list(iglacier);

		eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/training_glacier_' num2str(indexglaciernetcdf) '.mat'])
		eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/validation_glacier_' num2str(indexglaciernetcdf) '.mat'])

		updateclimatenumber;

		training_meltparam_binary=zeros(size(training_meltparam));
		training_meltparam_binary(1:30*85)=1;
		validation_meltparam_binary=zeros(size(validation_meltparam));
		validation_meltparam_binary(1:10*85)=1;
		std_global_training=zeros(size(training_meltparam));
		pred_global_training=zeros(size(training_meltparam));
		ycit_global_training=zeros(size(training_meltparam));
		pred_global_validation=zeros(size(validation_meltparam));
		ycit_global_validation=zeros(size(validation_meltparam));
		nugget_chosen=zeros(85,1);

		for icase=1:length(training_climate_number),
			if strcmp(training_climate{icase},'NorESM'),
				training_climate_number(icase)=meanoceantemp(1);
			elseif strcmp(training_climate{icase},'HadGEM2'),
				training_climate_number(icase)=meanoceantemp(2);
			elseif strcmp(training_climate{icase},'CSIRO'),
				training_climate_number(icase)=meanoceantemp(3);
			elseif strcmp(training_climate{icase},'MIROC'),
				training_climate_number(icase)=meanoceantemp(4);
			elseif strcmp(training_climate{icase},'IPSL'),
				training_climate_number(icase)=meanoceantemp(5);
			elseif strcmp(training_climate{icase},'CCSM4'),
				training_climate_number(icase)=meanoceantemp(6);
			elseif strcmp(training_climate{icase},'CESM2'),
				training_climate_number(icase)=meanoceantemp(7);
			elseif strcmp(training_climate{icase},'UKESM'),
				training_climate_number(icase)=meanoceantemp(8);
			elseif strcmp(training_climate{icase},'CNRM-ESM'),
				training_climate_number(icase)=meanoceantemp(9);
			elseif strcmp(training_climate{icase},'CNRM-CM6'),
				training_climate_number(icase)=meanoceantemp(10);
			else 
				training_climate{icase}
				error('model not recognized');
			end
		end

		for icase=1:length(validation_climate_number),
			if strcmp(validation_climate{icase},'NorESM'),
				validation_climate_number(icase)=meanoceantemp(1);
			elseif strcmp(validation_climate{icase},'HadGEM2'),
				validation_climate_number(icase)=meanoceantemp(2);
			elseif strcmp(validation_climate{icase},'CSIRO'),
				validation_climate_number(icase)=meanoceantemp(3);
			elseif strcmp(validation_climate{icase},'MIROC'),
				validation_climate_number(icase)=meanoceantemp(4);
			elseif strcmp(validation_climate{icase},'IPSL'),
				validation_climate_number(icase)=meanoceantemp(5);
			elseif strcmp(validation_climate{icase},'CCSM4'),
				validation_climate_number(icase)=meanoceantemp(6);
			elseif strcmp(validation_climate{icase},'CESM2'),
				validation_climate_number(icase)=meanoceantemp(7);
			elseif strcmp(validation_climate{icase},'UKESM'),
				validation_climate_number(icase)=meanoceantemp(8);
			elseif strcmp(validation_climate{icase},'CNRM-ESM'),
				validation_climate_number(icase)=meanoceantemp(9);
			elseif strcmp(validation_climate{icase},'CNRM-CM6'),
				validation_climate_number(icase)=meanoceantemp(10);
			else 
				validation_climate{icase}
				error('model not recognized');
			end
		end

		for i=1:85,
			i
			tab_training=table(training_ice(i:85:end),training_climate_number(i:85:end),training_meltparam_binary(i:85:end),training_values(i:85:end),'VariableNames', {'ice', 'climate','meltparam','values'});
			tab_validation=table(validation_ice(i:85:end),validation_climate_number(i:85:end),validation_meltparam_binary(i:85:end),'VariableNames', {'ice', 'climate','meltparam'});
			nugget=[0.001, 0.002, 0.005, 0.01, 0.02, 0.03, 0.05, 0.07, 0.1, 0.15, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0];
			rmse_training=[]; rmse_validation=[];
			%Thwaites glacier (#7): exponential fit
			%Pine Island glacier (#2):  exponential
			%Totten Glacier (#83):  exponential
			%Getz Land (#13):  exponential
			%Kamb (#18): matern32
			%MacAyeal (#20): matern32
			%Institute (#186) matern32
			%Moller (#172):  matern32
			%Leopold and Astrid Land (#92):  exponential
			%Moscow University (#82):  exponential
			%Whillans (#25): matern32
			gpr2=fitrgp(tab_training,'values','KernelFunction','matern32');
			eval(['gpr_' num2str(i) '=fitrgp(tab_training,''values'',''KernelFunction'',''matern32'');'])
			nugget_chosen2(i)=gpr2.Sigma;
			[ytraining,ystdtraining,ycit] = predict(gpr2,tab_training);
			[yvalidation,~,ycitvalidation] = predict(gpr2,tab_validation);
			std_global_training(i:85:end)=ystdtraining;
			pred_global_training(i:85:end)=ytraining;
			ycit_global_training(i:85:end)=ycit(:,2)-ycit(:,1);
			pred_global_validation(i:85:end)=yvalidation;
			ycit_global_validation(i:85:end)=ycitvalidation(:,2)-ycitvalidation(:,1);
			eval(['gpr_' num2str(i) '=gpr2;'])

			eval(['save Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/gpr_final_year_' num2str(i) ' gpr_' num2str(i) ''])

		end

		sqrt(mean((pred_global_training-training_values).^2))
		sqrt(mean((pred_global_validation-validation_values).^2))
		sqrt(mean(ycit_global_training.^2))
		sqrt(mean(ycit_global_validation.^2))

	end

	eval(['save Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/EmulationResultsYearlyEmulation pred_global_training ycit_global_training pred_global_validation ycit_global_validation'])

	plot(validation_values); hold on; plot(pred_global_validation); plot(pred_global_validation-ycit_global_validation,'k--'); plot(pred_global_validation+ycit_global_validation,'k--')

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

	experiments_list={'exp01','exp05',...
		'exp02','exp06',...
		'exp04','exp08',...
		'expA1','expA5',...
		'expA2','expA6',...
		'expA3','expA7',...
		'expB1','expB6',...
		'expB3','expB8',...
		'expB4','expB9',...
		'expB5','expB10'};

	experiments_table={'exp01','exp05';...
		'exp02','exp06';...
		'exp04','exp08';...
		'expA1','expA5';...
		'expA2','expA6';...
		'expA3','expA7';...
		'expB1','expB6';...
		'expB3','expB8';...
		'expB4','expB9';...
		'expB5','expB10'};
	
	climatemodel_list={'NorESM','NorESM',...
		'MIROC','MIROC',...
		'CCSM4','CCSM4',...
		'HadGEM2','HadGEM2',...
		'CSIRO','CSIRO',...
		'IPSL','IPSL',...
		'CNRM-ESM','CNRM-ESM',...
		'UKESM','UKESM',...
		'CESM2','CESM2',...
		'CNRM-CM6','CNRM-CM6'};

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

	for iglacier=6, %length(glacier_list),
		indexglaciernetcdf=glacier_list(iglacier);

		climatemodel_results={};
		icesheetmodel_results={};
		icesheetmodellong_results={};
		meltparam_results={};
		meltparam_results_binary=[];
		emulated_results=[];
		individual_results_limnsw_timeseries=[];
		individual_results_std_timeseries=[];
		icesheetmodel_longlist={};
		icesheetmodel_longlist{end+1}='AWI_PISM1_open';
		nummodel=1; %initialize at 1 because AWI put by default
		existing_results=zeros(21,10);

		%Load data and gprs
		eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/training_glacier_' num2str(indexglaciernetcdf) '.mat'])
		eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/validation_glacier_' num2str(indexglaciernetcdf) '.mat'])
		eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/test_glacier_' num2str(indexglaciernetcdf) '.mat'])

		for i=1:85,
			eval(['load Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/gpr_final_year_' num2str(i) ''])
		end

		training_meltparam_binary=zeros(size(training_meltparam));
		training_meltparam_binary(1:30*85)=1;
		validation_meltparam_binary=zeros(size(validation_meltparam));
		validation_meltparam_binary(1:10*85)=1;
		test_meltparam_binary=zeros(size(test_meltparam));
		test_meltparam_binary(1:4*85)=1;

		climatecase=2;
		updateclimatenumber;

		for icase=1:length(training_climate_number),
			if strcmp(training_climate{icase},'NorESM'),
				training_climate_number(icase)=meanoceantemp(1);
			elseif strcmp(training_climate{icase},'HadGEM2'),
				training_climate_number(icase)=meanoceantemp(2);
			elseif strcmp(training_climate{icase},'CSIRO'),
				training_climate_number(icase)=meanoceantemp(3);
			elseif strcmp(training_climate{icase},'MIROC'),
				training_climate_number(icase)=meanoceantemp(4);
			elseif strcmp(training_climate{icase},'IPSL'),
				training_climate_number(icase)=meanoceantemp(5);
			elseif strcmp(training_climate{icase},'CCSM4'),
				training_climate_number(icase)=meanoceantemp(6);
			elseif strcmp(training_climate{icase},'CESM2'),
				training_climate_number(icase)=meanoceantemp(7);
			elseif strcmp(training_climate{icase},'UKESM'),
				training_climate_number(icase)=meanoceantemp(8);
			elseif strcmp(training_climate{icase},'CNRM-ESM'),
				training_climate_number(icase)=meanoceantemp(9);
			elseif strcmp(training_climate{icase},'CNRM-CM6'),
				training_climate_number(icase)=meanoceantemp(10);
			else 
				training_climate{icase}
				error('model not recognized');
			end
		end

		for icase=1:length(validation_climate_number),
			if strcmp(validation_climate{icase},'NorESM'),
				validation_climate_number(icase)=meanoceantemp(1);
			elseif strcmp(validation_climate{icase},'HadGEM2'),
				validation_climate_number(icase)=meanoceantemp(2);
			elseif strcmp(validation_climate{icase},'CSIRO'),
				validation_climate_number(icase)=meanoceantemp(3);
			elseif strcmp(validation_climate{icase},'MIROC'),
				validation_climate_number(icase)=meanoceantemp(4);
			elseif strcmp(validation_climate{icase},'IPSL'),
				validation_climate_number(icase)=meanoceantemp(5);
			elseif strcmp(validation_climate{icase},'CCSM4'),
				validation_climate_number(icase)=meanoceantemp(6);
			elseif strcmp(validation_climate{icase},'CESM2'),
				validation_climate_number(icase)=meanoceantemp(7);
			elseif strcmp(validation_climate{icase},'UKESM'),
				validation_climate_number(icase)=meanoceantemp(8);
			elseif strcmp(validation_climate{icase},'CNRM-ESM'),
				validation_climate_number(icase)=meanoceantemp(9);
			elseif strcmp(validation_climate{icase},'CNRM-CM6'),
				validation_climate_number(icase)=meanoceantemp(10);
			else 
				validation_climate{icase}
				error('model not recognized');
			end
		end

		for icase=1:length(test_climate_number),
			if strcmp(test_climate{icase},'NorESM'),
				test_climate_number(icase)=meanoceantemp(1);
			elseif strcmp(test_climate{icase},'HadGEM2'),
				test_climate_number(icase)=meanoceantemp(2);
			elseif strcmp(test_climate{icase},'CSIRO'),
				test_climate_number(icase)=meanoceantemp(3);
			elseif strcmp(test_climate{icase},'MIROC'),
				test_climate_number(icase)=meanoceantemp(4);
			elseif strcmp(test_climate{icase},'IPSL'),
				test_climate_number(icase)=meanoceantemp(5);
			elseif strcmp(test_climate{icase},'CCSM4'),
				test_climate_number(icase)=meanoceantemp(6);
			elseif strcmp(test_climate{icase},'CESM2'),
				test_climate_number(icase)=meanoceantemp(7);
			elseif strcmp(test_climate{icase},'UKESM'),
				test_climate_number(icase)=meanoceantemp(8);
			elseif strcmp(test_climate{icase},'CNRM-ESM'),
				test_climate_number(icase)=meanoceantemp(9);
			elseif strcmp(test_climate{icase},'CNRM-CM6'),
				test_climate_number(icase)=meanoceantemp(10);
			else 
				test_climate{icase}
				error('model not recognized');
			end
		end

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

				eval(['isexp=is' expename ';'])
				if isexp,
					if strcmp(icesheetmodel_longlist{end},['' modelname '_open'])==0,
						icesheetmodel_longlist{end+1}=['' modelname '_open'];
						nummodel=nummodel+1;
					end
					existing_results(nummodel,iexp)=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'];
					eval(['results_limnsw=ncread(explimnsw_file,''ivaf_glacier_' num2str(indexglaciernetcdf) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
					eval(['results_tendacabfgr=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(indexglaciernetcdf) ''') *yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
					individual_results_limnsw_timeseries(end+1,:)=(results_limnsw(1:85)-cumsum(results_tendacabfgr(1:85))-(results_limnsw(1)-results_tendacabfgr(1)))/362.5;
					if strcmpi(modelname,'UCIJPL_ISSM') & strcmpi(expename,'exp01'),
						individual_results_limnsw_timeseries(end,44)=(individual_results_limnsw_timeseries(end,43)+individual_results_limnsw_timeseries(end,45))/2;
					end
					individual_results_std_timeseries(end+1,:)=zeros(85,1);
					icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
					icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_open']},1,85)];
					climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp-1}},1,85)];
					emulated_results(end+1,:)=0*ones(1,85);
					meltparam_results=[meltparam_results;repmat({'open'},85,1)];
					meltparam_results_binary(end+1,:)=1*ones(1,85);
				else %no existing results so emulate them
					if isexp01==1,
						%there is some data for this model so you can emulate the missing simulations
						meltparam_tosimulate=repmat({'open'},85,1);
						meltparam_tosimulate_binary=ones(85,1);
						icesheetmodel_tosimulate=repmat({modelname},85,1);
						if strcmp(climatemodel_list{2*iexp-1},'NorESM'),
							climate_number=meanoceantemp(1);
						elseif strcmp(climatemodel_list{2*iexp-1},'HadGEM2'),
							climate_number=meanoceantemp(2);
						elseif strcmp(climatemodel_list{2*iexp-1},'CSIRO'),
							climate_number=meanoceantemp(3);
						elseif strcmp(climatemodel_list{2*iexp-1},'MIROC'),
							climate_number=meanoceantemp(4);
						elseif strcmp(climatemodel_list{2*iexp-1},'IPSL'),
							climate_number=meanoceantemp(5);
						elseif strcmp(climatemodel_list{2*iexp-1},'CCSM4'),
							climate_number=meanoceantemp(6);
						elseif strcmp(climatemodel_list{2*iexp-1},'CESM2'),
							climate_number=meanoceantemp(7);
						elseif strcmp(climatemodel_list{2*iexp-1},'UKESM'),
							climate_number=meanoceantemp(8);
						elseif strcmp(climatemodel_list{2*iexp-1},'CNRM-ESM'),
							climate_number=meanoceantemp(9);
						elseif strcmp(climatemodel_list{2*iexp-1},'CNRM-CM6'),
							climate_number=meanoceantemp(10);
						else 
							climatemodel_list{2*iexp-1}
							error('model not recognized');
						end
						climatemodel_tosimulate=repmat({climatemodel_list{2*iexp-1}},85,1);
						climatemodel_tosimulate_number=climate_number*ones(85,1);
						time_tosimulate=1:85;
						yemul=[];
						ystdemul=[];
						for timesim=1:85,
							table_indiv=table(icesheetmodel_tosimulate(timesim),climatemodel_tosimulate_number(timesim),meltparam_tosimulate_binary(timesim),'VariableNames', {'ice', 'climate','meltparam'});
							eval(['[yemul(timesim),ystd_emul(timesim),yint_emul]= predict(gpr_' num2str(timesim) ',table_indiv);']);
						end

						%Save with the overall results
						emulated_results(end+1,:)=ones(1,85);
						individual_results_limnsw_timeseries(end+1,:)=yemul;
						individual_results_std_timeseries(end+1,:)=ystd_emul;
						icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
						icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_open']},1,85)];
						climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp-1}},1,85)];
						meltparam_results=[meltparam_results;repmat({'open'},85,1)];
						meltparam_results_binary(end+1,:)=1*ones(1,85);
					else %there is no data for this standard model so don't create one from nothing
						%do nothing
					end
				end

			end %end of model
		end %end of isexp
		for imodel=1:length(model_list),
			modelname=model_list{imodel};
			for iexp=1:length(experiments_list)/2,
				expename=experiments_list{2*iexp};
				specifics;

				eval(['isexp=is' expename ';'])
				if isexp,
					if strcmp(icesheetmodel_longlist{end},['' modelname '_standard'])==0,
						icesheetmodel_longlist{end+1}=['' modelname '_standard'];
						nummodel=nummodel+1;
					end
					existing_results(nummodel,iexp)=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'];
					eval(['results_limnsw=ncread(explimnsw_file,''ivaf_glacier_' num2str(indexglaciernetcdf) ''') *ice_density/(10^9*1000);']) %from m^3 to Gt
					eval(['results_tendacabfgr=ncread(exptendacabfgr_file,''smbgr_glacier_' num2str(indexglaciernetcdf) ''') *yearday_model*3600*24/(10^9*1000);']) %from kg/s to Gt/yr
					individual_results_limnsw_timeseries(end+1,:)=(results_limnsw(1:85)-cumsum(results_tendacabfgr(1:85))-(results_limnsw(1)-results_tendacabfgr(1)))/362.5;
					individual_results_std_timeseries(end+1,:)=zeros(85,1);
					icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
					icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_standard']},1,85)];
					climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp}},1,85)];
					emulated_results(end+1,:)=0*ones(1,85);
					meltparam_results=[meltparam_results;repmat({'standard'},85,1)];
					meltparam_results_binary(end+1,:)=0*ones(1,85);
				else %no existing results so emulate them
					if isexp05==1,
						%there is some data for this model so you can emulate the missing simulations
						meltparam_tosimulate=repmat({'standard'},85,1);
						meltparam_tosimulate_binary=zeros(85,1);
						icesheetmodel_tosimulate=repmat({modelname},85,1);
						if strcmp(climatemodel_list{2*iexp},'NorESM'),
							climate_number=meanoceantemp(1);
						elseif strcmp(climatemodel_list{2*iexp},'HadGEM2'),
							climate_number=meanoceantemp(2);
						elseif strcmp(climatemodel_list{2*iexp},'CSIRO'),
							climate_number=meanoceantemp(3);
						elseif strcmp(climatemodel_list{2*iexp},'MIROC'),
							climate_number=meanoceantemp(4);
						elseif strcmp(climatemodel_list{2*iexp},'IPSL'),
							climate_number=meanoceantemp(5);
						elseif strcmp(climatemodel_list{2*iexp},'CCSM4'),
							climate_number=meanoceantemp(6);
						elseif strcmp(climatemodel_list{2*iexp},'CESM2'),
							climate_number=meanoceantemp(7);
						elseif strcmp(climatemodel_list{2*iexp},'UKESM'),
							climate_number=meanoceantemp(8);
						elseif strcmp(climatemodel_list{2*iexp},'CNRM-ESM'),
							climate_number=meanoceantemp(9);
						elseif strcmp(climatemodel_list{2*iexp},'CNRM-CM6'),
							climate_number=meanoceantemp(10);
						else 
							climatemodel_list{2*iexp}
							error('model not recognized');
						end
						climatemodel_tosimulate=repmat({climatemodel_list{2*iexp}},85,1);
						climatemodel_tosimulate_number=climate_number*ones(85,1);
						time_tosimulate=1:85;
						yemul=[];
						ystd_emul=[];
						for timesim=1:85,
							table_indiv=table(icesheetmodel_tosimulate(timesim),climatemodel_tosimulate_number(timesim),meltparam_tosimulate_binary(timesim),'VariableNames', {'ice', 'climate','meltparam'});
							eval(['[yemul(timesim),ystd_emul(timesim),yint_emul]= predict(gpr_' num2str(timesim) ',table_indiv);']);
						end

						%Save with the overall results
						emulated_results(end+1,:)=1*ones(1,85);
						individual_results_limnsw_timeseries(end+1,:)=yemul;
						individual_results_std_timeseries(end+1,:)=ystd_emul;
						icesheetmodel_results=[icesheetmodel_results;repmat({['' modelname '']},1,85)];
						icesheetmodellong_results=[icesheetmodellong_results;repmat({['' modelname '_standard']},1,85)];
						climatemodel_results=[climatemodel_results;repmat({climatemodel_list{2*iexp}},1,85)];
						meltparam_results=[meltparam_results;repmat({'standard'},85,1)];
						meltparam_results_binary(end+1,:)=0*ones(1,85);
					else %there is no data for this standard model so don't create one from nothing
						%do nothing
					end
				end

			end %end of model
		end %end of isexp

		eval(['save Data/EmulationData/Glacier_' num2str(indexglaciernetcdf) '/glacier_dyn_loss_emulated  individual_results_limnsw_timeseries individual_results_std_timeseries icesheetmodel_results icesheetmodellong_results climatemodel_results emulated_results meltparam_results_binary meltparam_results'])

	end
end %}}}
