Index: /issm/trunk-jpl/src/m/plot/colormaps/paraview.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/colormaps/paraview.m	(revision 24078)
+++ /issm/trunk-jpl/src/m/plot/colormaps/paraview.m	(revision 24079)
@@ -40,3 +40,3 @@
 	n = 256;
 end
-y = interp1(1:l,J(:,2:4)/255,linspace(1,l,n),'*linear');
+cmap = interp1(1:l,J(:,2:4)/255,linspace(1,l,n),'*linear');
Index: /issm/trunk-jpl/src/m/plot/plotchannels.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plotchannels.m	(revision 24079)
+++ /issm/trunk-jpl/src/m/plot/plotchannels.m	(revision 24079)
@@ -0,0 +1,104 @@
+function plotchannels(md,channelarea,varargin)
+%PLOTCHANNELS - plot GlaDS channel area
+%
+%   Usage:
+%      plotchannels(md,channelarea,options)
+%
+%   List of supported options
+%      - 'min' minimum channel area displayed (default is max(channelarea))
+%      - 'max' maximum channel area displayed (default is min(channelarea))
+%      - 'colormap' colormap used (default is 'gray')
+%      - 'linewidth' linewidth (default is 2)
+%
+%   Example:
+%      plotchannels(md,md.results.TransientSolution(end).ChannelArea,'min',15,'max',36);
+
+%Process options
+options = pairoptions(varargin{:});
+
+%define level
+level = channelarea;
+
+%Some processing
+Min = getfieldvalue(options,'min',min(level));
+Max = getfieldvalue(options,'max',max(level));
+
+%Create colormap
+options=addfielddefault(options,'colormap',flipud(gray()));
+palette = getcolormap(options);
+numcolors=size(palette,1);
+levels=round_ice(linspace(Min,Max,numcolors+1),2);
+
+colorind=zeros(length(level),1);
+for i=1:numcolors
+	pos=find((level>=levels(i)) & (level<=levels(i+1)) );
+	colorind(pos)=i;
+end
+colorind(find(level>levels(end)))=numcolors;
+
+%Reconstruct edges
+if 0,
+	edges = [md.mesh.edges(:,1:2) md.mesh.edges(:,1)]';
+else
+	tic
+	%Maximum number of edges
+	maxnbf = 3*md.mesh.numberofelements;
+	%Initialize intermediaries
+	edges = zeros(maxnbf,2);
+	%Chaining algorithm
+	head_minv = -1*ones(md.mesh.numberofvertices,1);
+	next_face = zeros(maxnbf,1);
+	%Initialize number of faces
+	nbf = 0;
+	for i=1:md.mesh.numberofelements
+		for j=1:3
+         %Get the two indices of the edge number j of the ith triangle
+         v1 = md.mesh.elements(i,j);
+			if(j==3)
+				v2 = md.mesh.elements(i,1);
+			else
+				v2 = md.mesh.elements(i,j+1);
+			end
+			%sort
+         if(v2<v1)
+            v3=v2; v2=v1; v1=v3;
+			end
+         %This edge a priori has not been processed yet
+         exists = false;
+			%Go through all processed faces connected to v1 and check whether we have seen this edge yet
+			e=head_minv(v1);
+			while(e~=-1)
+            if(edges(e,2)==v2)
+               exists = true;
+               break;
+				end
+				e=next_face(e);
+			end
+			%If this edge is new, add it to the lists
+         if(~exists)
+            %Update edges
+            edges(nbf+1,1) = v1;
+            edges(nbf+1,2) = v2;
+            %Update chain
+            next_face(nbf+1) = head_minv(v1);
+            head_minv(v1)    = nbf+1;
+            %Increase number of faces
+            nbf=nbf+1;
+			end
+		end
+	end
+	edges = edges(1:nbf,:);
+	toc
+end
+
+%Change edges formatting so that plot looks ok
+edges = [edges(:,1:2) edges(:,1)]';
+
+%Loop over all levels and plot
+for i=1:numcolors
+	pos=find(colorind==i);
+	hprime=plot(md.mesh.x(edges(:,pos)),md.mesh.y(edges(:,pos)),...
+		'-','Color',palette(i,:),'LineWidth',2);
+	if i==1; hold on; end
+end
+hold off