
clear all;
addpath('../Data','../Functions');

steps=[1]; % [1:8]

if any(steps==1) 
	disp('   Step 1: Global mesh creation');

	resolution=300;			% [km] 
	radius = 6.371012*10^3;	% [km] 

	md=model; 
	md.mask=maskpsl(); 
	md.mesh=gmshplanet('radius',radius,'resolution',resolution);

	md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 

	save ./Models/SlrGRACE_Mesh md;
	
	plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');

end 

if any(steps==2) 
	disp('   Step 2: Global mesh creation: refined coastlines');

	numrefine=1;
	resolution=150*1e3;			% inital resolution [m] 
	radius = 6.371012*10^6;		% mean radius of Earth, m
	mindistance_coast=150*1e3;	% coastal resolution [m] 
	mindistance_source=75*1e3; % source resolution [m] 
	mindistance_land=300*1e3;	% resolution on the continents [m]
	maxdistance=600*1e3;			% max element size (on mid-oceans) [m]

	md=model; 
	md.mask=maskpsl(); 
	md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] 

	for i=1:numrefine,

		md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 

		distance=zeros(md.mesh.numberofvertices,1);

		pos=find(~md.mask.ocean_levelset);	coaste.lat=md.mesh.lat(pos);	coaste.long=md.mesh.long(pos);  
		pos=find(md.mask.ocean_levelset);	coasto.lat=md.mesh.lat(pos);	coasto.long=md.mesh.long(pos);  

		for j=1:md.mesh.numberofvertices
			phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; 
			if md.mask.ocean_levelset(j),
				phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; 
				deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
				d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
			else
				phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; 
				deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
				d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
			end
			distance(j)=min(d);
		end
		pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
		
		pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); 
		distance(pos2)=mindistance_land; 

		domain = 'Greenland'; %'Antarctica'; %'HMA'; %'Alaska'; %'Glaciers'
		mask = domain_mask(md.mesh.lat,md.mesh.long,domain);  
		distance(mask>0)=mindistance_source;

		dist=min(maxdistance,distance); 
		md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist);

	end

	md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 

	save ./Models/SlrGRACE_Mesh_refined md;
	
	plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');

end 

if any(steps==3) 
	disp('   Step 3: Define loads in meters of ice height equivalent');
	md = loadmodel('./Models/SlrGRACE_Mesh_refined');

	year_month = 2007+15/365;
	time_range = [year_month year_month]; 
	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); 

	index = md.mesh.elements; 
	lat=90-md.mesh.lat; 
	lon=md.mesh.long; 
	lon(lon<0)=180+(180+lon(lon<0)); 
	
	ax_0=lat(index(:,1)); ay_0=lon(index(:,1)); 
	bx_0=lat(index(:,2)); by_0=lon(index(:,2)); 
	cx_0=lat(index(:,3)); cy_0=lon(index(:,3)); 
	
	for ii=1:md.mesh.numberofelements
		if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180)
			if ay_0(ii)==0
				ay_0(ii)=360;
			end 
			if by_0(ii)==0
				by_0(ii)=360; 
			end 
			if cy_0(ii)==0 
				cy_0(ii)=360; 
			end
		end 
	end
	
	ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2; 
	by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2; 
	cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2; 
	
	ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2; 
	by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2; 
	cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2; 
	
	lat_element=(ax_0+bx_0+cx_0)/3; 
	lon_element=(ay_0+by_0+cy_0)/3;

	lat_element_0 = 90-lat_element;		lon_element_0 = lon_element;
	lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180;
	
	domain = 'Greenland'; %'HMA'; %'Alaska'; %'Antarctica'; %'Glaciers'; 
	mask = domain_mask(lat_element_0,lon_element_0,domain);  

	md.slr.deltathickness = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice; 

	save ./Models/SlrGRACE_Loads md; 

	plotmodel(md,'data',md.slr.deltathickness,'view',[90 90],'title','Ice height equivalent [m]'); 

end 

if any(steps==4) 
	disp('   Step 4: Parameterization');
	md = loadmodel('./Models/SlrGRACE_Loads');

	nlove=10001;
	md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[];
	md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[];
	md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[];

	md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); 
	md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 
	pos=find(md.slr.deltathickness~=0);
	md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; 
	md.mask.land_levelset = 1-md.mask.ocean_levelset; 

	md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
	md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
	md.slr.ocean_area_scaling = 1;

	md.slr.spcthickness=NaN*zeros(md.mesh.numberofvertices,1);
	md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
	md.slr.Ugia=zeros(md.mesh.numberofvertices,1);

	md.slr.geodetic=1; 

	di=md.materials.rho_ice/md.materials.rho_water;
	md.geometry.thickness=ones(md.mesh.numberofvertices,1);
	md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
	md.geometry.base=md.geometry.surface-md.geometry.thickness;
	md.geometry.bed=md.geometry.base;
	
	md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
	md.materials.rheology_B=paterson(md.initialization.temperature);
	md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
	
	md.miscellaneous.name='SlrGRACE';
	
	save ./Models/SlrGRACE_Parameterization md; 

end 

if any(steps==5) 
	disp('   Step 5: Solve Slr solver');
	md = loadmodel('./Models/SlrGRACE_Parameterization');

	md.slr.rigid=1; 
	md.slr.elastic=1; 
	md.slr.rotation=1;

	md.cluster=generic('name',oshostname(),'np',3); 
	md.verbose=verbose('111111111');

	md=solve(md,'Slr');

	save ./Models/SlrGRACE_Solution md; 

end 

if any(steps==6) 
	disp('   Step 6: Plot solutions');
	md = loadmodel('./Models/SlrGRACE_Solution');

	sol1 = md.slr.deltathickness*100;							% [cm] 
	sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] 
	
	sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; 
	fig_name={'Fig_dH.pdf','Fig_slr.pdf'}; 

	res = 1.0; 

	[lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
	sol_grid = zeros(size(lat_grid)); 

	for kk=1:2
		sol=eval(sprintf('sol%d',kk));
	
		if length(sol)==md.mesh.numberofelements 
			for jj=1:md.mesh.numberofelements
				ii=(jj-1)*3;
				pp(ii+1:ii+3)=md.mesh.elements(jj,:);
			end
			for jj=1:md.mesh.numberofvertices
				pos=ceil(find(pp==jj)/3); 
				temp(jj)=mean(sol(pos)); 
			end
			sol=temp'; 
		end 

		F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
		F.Method = 'linear';
		F.ExtrapolationMethod = 'linear'; 
		
		sol_grid = F(lat_grid, lon_grid);
		sol_grid(isnan(sol_grid))=0; 
		sol_grid(lat_grid>85 & sol_grid==0) = NaN; 

		set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
		figure1=figure('Position', [100, 100, 1000, 500]); 
		gcf; load coast; cla; 
		pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
		if (kk==1)
			geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); 
		else
			geoshow(lat,long,'DisplayType','polygon','FaceColor',[0.75 0.75 0.75]); 
		end 
		plot(long,lat,'k'); hold off; 
		c1=colorbar;
		colormap('haxby'); 
		caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]); 
		xlim([-180 180]); 
		ylim([-90 90]); 
		grid on; 
		title(sol_name(kk)); 
		set(gcf,'color','w');
		%export_fig(fig_name{kk}); 
	end
end 

if any(steps==7) 
	disp('   Step 7: Transient run');
	md = loadmodel('./Models/SlrGRACE_Parameterization');

	disp('Projecting  loads onto the mesh...');
	time_range = [2005+15/365 2014+350/365]; 
	%time_range = 2007 + [15 350]/365; 
	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); 

	[ne,nt]=size(water_load); 
   md.slr.deltathickness = zeros(ne+1,nt);
	md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; 
	md.slr.deltathickness(ne+1,:)=[1:nt]; % times 
	
	md.transient.issmb=0;
	md.transient.ismasstransport=0; 
	md.transient.isstressbalance=0; 
	md.transient.isthermal=0;
	md.transient.isslr=1;
	
	md.timestepping.start_time=-1;
   md.timestepping.final_time=nt; 
	md.timestepping.time_step=1; 
   md.timestepping.interp_forcings=0;
	md.settings.output_frequency=1;

	md.cluster=generic('name',oshostname(),'np',3); 
	md.verbose=verbose('111111111');

	md=solve(md,'tr');

	save ./Models/SlrGRACE_Transient md; 

end 

if any(steps==8) 
	disp('   Step 8: Plot transient');
	md = loadmodel('./Models/SlrGRACE_Transient');

	time = md.slr.deltathickness(end,:); 

	tide_x = 37+29/60;  % degree N
	tide_y = 126+20/60; % degree E

	for tt=1:length(time)
		gmsl(tt) = md.results.TransientSolution(tt).SealevelRSLEustatic*1000; 
	   sol = (md.results.TransientSolution(tt+1).Sealevel-md.results.TransientSolution(tt).Sealevel)*1000; 
		slr_incheon(tt) = griddata(md.mesh.lat,md.mesh.long,sol,tide_x,tide_y);
	end

	plot(2005+time/12,gmsl-gmsl(1),'-*','linewidth',3); grid on; hold on; 
	plot(2005+time/12,slr_incheon-slr_incheon(1),'-*r','linewidth',3); hold off; 
	xlabel('Year'); 
	ylabel('GMSL [mm]');
	set(gcf,'color','w');
	legend('GMSL','Incheon'); 

end 

