Index: /issm/trunk-jpl/examples/SlrFarrell/runme.m
===================================================================
--- /issm/trunk-jpl/examples/SlrFarrell/runme.m	(revision 22812)
+++ /issm/trunk-jpl/examples/SlrFarrell/runme.m	(revision 22813)
@@ -1,5 +1,5 @@
 
 clear all;
-steps=[2]; % 
+steps=[5]; % 
 
 if any(steps==1) % Global mesh creation {{{ 
@@ -70,4 +70,7 @@
 	% initial sea-level: 1 m RSL everywhere. 
 	md.slr.sealevel=md.mask.ocean_levelset; 
+	
+	md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
+	md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
 
 	save ./Models/SlrFarrell.Loads md; 
@@ -85,14 +88,16 @@
 	% Love numbers and reference frame: CF or CM (choose one!)  
 	nlove=10001;	% up to 10,000 degree 
-	md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[];
-	md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[];
-	md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[];
+	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)=[];
 
 	% Mask: for computational efficiency only those elements that have loads are convolved! 
-	md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded 
-	md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 
-	pos=find(md.slr.deltathickness~=0);
-	md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads 
 	md.mask.land_levelset = 1-md.mask.ocean_levelset; 
+	% fake ice load in one element!  
+	md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); % no ice 
+	md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1); % floated... 
+	pos=find(md.mesh.lat <-80);
+	md.mask.ice_levelset(pos(1))=-1; % ice yes!  
+	md.mask.groundedice_levelset(pos(1))=1; % ice grounded!  
 
 	%% IGNORE BUT DO NOT DELETE %% {{{
@@ -119,7 +124,4 @@
 	md = loadmodel('./Models/SlrFarrell.Parameterization');
 
-	% Request outputs 
-	md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'};
-	
 	% Cluster info 
 	md.cluster=generic('name',oshostname(),'np',3); 
@@ -137,12 +139,7 @@
 	md = loadmodel('./Models/SlrFarrell.Solution');
 
-	% loads and solutions. 
-	sol1 = md.slr.deltathickness*100; % WEH cm 
-	sol2 = md.results.SlrSolution.SlrUmotion*1000; % [mm] 
-	sol3 = md.results.SlrSolution.SlrNmotion*1000; % [mm] 
-	sol4 = md.results.SlrSolution.SlrEmotion*1000; % [mm] 
-	sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',...
-		'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; 
-
+	% solutions. 
+	sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m)  
+	
 	res = 1.0; % degree 
 
@@ -151,50 +148,34 @@
 	sol_grid = zeros(size(lat_grid)); 
 
-	for kk=1:4 
-		sol=eval(sprintf('sol%d',kk));
+	% Make a interpolation object
+	F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
+	F.Method = 'linear';
+	F.ExtrapolationMethod = 'linear'; 
 	
-		% if data are on elements, map those on to the vertices {{{
-		if length(sol)==md.mesh.numberofelements 
-			% map on to the vertices 
-			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 % }}}
+	% Do the interpolation to get gridded solutions... 
+	sol_grid = F(lat_grid, lon_grid);
+	sol_grid(isnan(sol_grid))=0; 
+	sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan 
 
-		% Make a interpolation object
-		F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
-		F.Method = 'linear';
-		F.ExtrapolationMethod = 'linear'; 
+	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;
+	contourf(lon_grid,lat_grid,sol_grid,[96 98 100 102 104 105]); shading flat; hold on;
+	geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); 
+	plot(long,lat,'k'); hold off; 
+	% define colormap, caxis, xlim etc {{{
+	c1=colorbar;
+	%colormap('haxby'); 
+	%caxis([96 104]); 
+	xlim([-180 180]); 
+	ylim([-90 90]); 
+	% }}} 
+	grid on; 
+	title('Relative sea-level [mm]'); 
+	set(gcf,'color','w');
 
-		% Do the interpolation to get gridded solutions... 
-		sol_grid = F(lat_grid, lon_grid);
-		sol_grid(isnan(sol_grid))=0; 
-		sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to 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; 
-		plot(long,lat,'k'); hold off; 
-		c1=colorbar;
-		colormap(jet); 
-		xlim([-180 180]); 
-		ylim([-90 90]); 
-		grid on; 
-		title(sol_name(kk)); 
-		set(gcf,'color','w');
-
-		%export_fig('Fig5.pdf'); 
-	end
+	%export_fig('Fig5.pdf'); 
 
 end % }}} 
 
-