Index: /issm/trunk-jpl/src/m/plot/plot_gridded.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_gridded.m	(revision 15398)
+++ /issm/trunk-jpl/src/m/plot/plot_gridded.m	(revision 15399)
@@ -23,4 +23,5 @@
 %Interpolating data on grid
 [x_m y_m data_grid]=InterpFromMeshToGrid(elements,x,y,data,xlim(1),ylim(2),post,post,round(diff(ylim)/post),round(diff(xlim)/post),NaN);
+data_grid_save = data_grid;
 if size(data_grid,1)<3 | size(data_grid,2)<3,
 	error('data_grid size too small in plot_gridded, check posting and units');
@@ -34,4 +35,5 @@
 
 %Process data_grid: add white in NaN and correct caxis accordingly
+[data_nani data_nanj]=find(isnan(data_grid) | data_grid==-9999);
 if exist(options,'caxis'),
 	caxis_opt=getfieldvalue(options,'caxis');
@@ -44,8 +46,4 @@
 	data_max=max(data_grid(:));
 end
-options = changefieldvalue(options,'cbYLim',[data_min data_max]);
-white   = data_min - (data_max-data_min)/(lenmap);
-options = changefieldvalue(options,'caxis',[white data_max]);
-data_grid(isnan(data_grid))=white;
 
 %Select plot area 
@@ -53,10 +51,35 @@
 
 %shading interp;
-if exist(options,'forcecolormap'),
-	image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(map)/(data_max-data_min))),map);
-	h=imagesc(xlim,ylim,image_rgb);
-else
-	h=imagesc(xlim,ylim,data_grid);
+image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(map)/(data_max-data_min))),map);
+if exist(options,'shaded'),
+	a    = -45;
+	scut = 0.2;
+	c    = 1;
+	% computes lighting from elevation gradient
+	[fx,fy] = gradient(data_grid_save,x_m,y_m);
+	fxy = -fx*sind(a) - fy*cosd(a);
+	clear fx fy % free some memory...
+	fxy(isnan(fxy)) = 0;
+
+	% computes maximum absolute gradient (median-style), normalizes, saturates and duplicates in 3-D matrix
+	r = repmat(max(min(fxy/nmedian(abs(fxy),1 - scut/100),1),-1),[1,1,3]);
+
+	% applies contrast using exponent
+	rp = (1 - abs(r)).^c;
+	image_rgb = image_rgb.*rp;
+
+	% lighter for positive gradient
+	k = find(r > 0);
+	image_rgb(k) = image_rgb(k) + (1 - rp(k));
 end
+
+% set novalues / NaN to black color
+if ~isempty(data_nani)
+	nancolor=getfieldvalue(options,'nancolor',[0 0 0]);
+	image_rgb(sub2ind(size(image_rgb),repmat(data_nani,1,3),repmat(data_nanj,1,3),repmat(1:3,size(data_nani,1),1))) = repmat(nancolor,size(data_nani,1),1);
+end
+
+%plot grid
+h=imagesc(xlim,ylim,image_rgb);
 axis xy
 
@@ -68,3 +91,7 @@
 
 %Apply options
+if ~isnan(data_min),
+	options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready
+end
+options=addfielddefault(options,'axis','xy equal off'); % default axis
 applyoptions(md,data,options);
