Index: /issm/trunk-jpl/examples/EsaGRACE/runme.m
===================================================================
--- /issm/trunk-jpl/examples/EsaGRACE/runme.m	(revision 22802)
+++ /issm/trunk-jpl/examples/EsaGRACE/runme.m	(revision 22803)
@@ -1,5 +1,5 @@
 
 clear all;
-steps=[5]; % [1:6]; 
+steps=[0]; % [1:5]; 
 
 if any(steps==0) % Download GRACE land_mass data {{{
@@ -42,10 +42,10 @@
 	md = loadmodel('./Models/EsaGRACE.Mesh');
 
+	% define time interval for analysis 
 	year_month = 2007+15/365;
-	time_min=year_month; 
-	time_max=year_month; 
-
+	time_range = [year_month year_month]; 
+	
 	% map GRACE water load on to the mesh for the seleted month(s) 
-	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_min,time_max); 
+	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); 
 	
 	md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent 
@@ -117,8 +117,11 @@
 	md = loadmodel('./Models/EsaGRACE.Solution');
 
-	% solutions. 
-	ur = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
-	un = md.results.EsaSolution.EsaNmotion*1000; % [mm] 
-	ue = md.results.EsaSolution.EsaEmotion*1000; % [mm] 
+	% loads and solutions. 
+	sol1 = md.esa.deltathickness*100; % WEH cm 
+	sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
+	sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] 
+	sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] 
+	sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',...
+		'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; 
 
 	res = 1.0; % degree 
@@ -128,44 +131,49 @@
 	sol_grid = zeros(size(lat_grid)); 
 
-	sol = ue; 
+	for kk=1:4 
+		sol=eval(sprintf('sol%d',kk));
+	
+		% 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 % }}}
 
-	% Make a interpolation object
-	F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
-	F.Method = 'linear';
-	F.ExtrapolationMethod = 'linear'; 
+		% Make a interpolation object
+		F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
+		F.Method = 'linear';
+		F.ExtrapolationMethod = 'linear'; 
 
-	% Do the interpolation to get gridded solutions... 
-	sol_grid = F(lat_grid, lon_grid);
-	sol_grid(isnan(sol_grid))=0; 
+		% 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 polar unphysical 0s to Nan 
-	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; 
+		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');
 
-	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; 
-	%caxis([-0.3 0.3])
-	hold on
-	plot(long,lat,'k');
-	%geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); 	
-	hold off; 
-	c1=colorbar;
-	colormap(jet); 
-	xlim([-180 180]); 
-	ylim([-90 90]); 
-	grid on; 
-	title(['Average change in relative sea-level [mm/yr]']);
-	set(gcf,'color','w');
-
-	%export_fig('Fig5.pdf'); 
+		%export_fig('Fig5.pdf'); 
+	end
 
 end % }}} 
 
-if any(steps==6) % {{{ Transient 
-	disp('   Step 6: Transient run');
-
-end % }}} 
-