Index: /issm/trunk-jpl/src/m/plot/plot_contour.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_contour.m	(revision 21455)
+++ /issm/trunk-jpl/src/m/plot/plot_contour.m	(revision 21456)
@@ -14,5 +14,5 @@
 
 %check is2d
-if ~is2d,
+if ~is2d & ~isplanet,
 	error('plot_contour error message: contour not supported for 3d meshes, project on a layer');
 end
@@ -115,4 +115,7 @@
 	y1=zeros(numelems,1);
 	y2=zeros(numelems,1);
+	z1=zeros(numelems,1);
+	z2=zeros(numelems,1);
+
 	edge_l=zeros(numelems,2);
 
@@ -129,4 +132,6 @@
 			y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
 			y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
+			z1(j)=z(Seg1(poselem(j),1))+weight1*(z(Seg1(poselem(j),2))-z(Seg1(poselem(j),1)));
+			z2(j)=z(Seg2(poselem(j),1))+weight2*(z(Seg2(poselem(j),2))-z(Seg2(poselem(j),1)));
 			edge_l(j,1)=Seg1_num(poselem(j));
 			edge_l(j,2)=Seg2_num(poselem(j));
@@ -138,4 +143,7 @@
 			y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
 			y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
+			z1(j)=z(Seg1(poselem(j),1))+weight1*(z(Seg1(poselem(j),2))-z(Seg1(poselem(j),1)));
+			z2(j)=z(Seg3(poselem(j),1))+weight3*(z(Seg3(poselem(j),2))-z(Seg3(poselem(j),1)));
+
 			edge_l(j,1)=Seg1_num(poselem(j));
 			edge_l(j,2)=Seg3_num(poselem(j));
@@ -147,4 +155,6 @@
 			y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
 			y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
+			z1(j)=z(Seg2(poselem(j),1))+weight2*(z(Seg2(poselem(j),2))-z(Seg2(poselem(j),1)));
+			z2(j)=z(Seg3(poselem(j),1))+weight3*(z(Seg3(poselem(j),2))-z(Seg3(poselem(j),1)));
 			edge_l(j,1)=Seg2_num(poselem(j));
 			edge_l(j,2)=Seg3_num(poselem(j));
@@ -161,5 +171,5 @@
 		%take the right edge of the second segment and connect it to the next segments if any
 		e1=edge_l(1,1);   e2=edge_l(1,2);
-		xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)];
+		xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)]; zc=[z1(1);z2(1)]; 
 
 		%erase the lines corresponding to this edge
@@ -167,4 +177,5 @@
 		x1(1)=[]; x2(1)=[];
 		y1(1)=[]; y2(1)=[];
+		z1(1)=[]; z2(1)=[];
 
 		[ro1,co1]=find(edge_l==e1);
@@ -173,5 +184,5 @@
 
 			if co1==1,
-				xc=[x2(ro1);xc]; yc=[y2(ro1);yc];
+				xc=[x2(ro1);xc]; yc=[y2(ro1);yc]; zc=[z2(ro1);zc];
 
 				%next edge:
@@ -179,5 +190,5 @@
 
 			else
-				xc=[x1(ro1);xc]; yc=[y1(ro1);yc];
+				xc=[x1(ro1);xc]; yc=[y1(ro1);yc]; zc=[z1(ro1);zc];
 
 				%next edge:
@@ -189,4 +200,5 @@
 			x1(ro1)=[]; x2(ro1)=[];
 			y1(ro1)=[]; y2(ro1)=[];
+			z1(ro1)=[]; z2(ro1)=[];
 
 			%next connection
@@ -200,10 +212,10 @@
 
 			if co2==1,
-				xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];
+				xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)]; zc=[zc;z2(ro2)];
 
 				%next edge:
 				e2=edge_l(ro2,2);
 			else
-				xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)];
+				xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)]; zc=[zc;z1(ro2)];
 
 				%next edge:
@@ -215,4 +227,5 @@
 			x1(ro2)=[]; x2(ro2)=[];
 			y1(ro2)=[]; y2(ro2)=[];
+			z1(ro2)=[]; z2(ro2)=[];
 
 			%next connection
@@ -221,13 +234,22 @@
 
 		%we now have one subcontour ready to be plotted
-		zc=level*ones(length(xc)+1,1);
 		if getfieldvalue(options,'contouronly',0),
-			h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'Zdata',zc,'Cdata',zc,'facecolor','none','edgecolor','flat','linewidth',linewidth)];
+			if isplanet,
+				h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'Zdata',[zc;NaN],'facecolor','none','linewidth',linewidth)];
+			else
+				h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'Zdata',zc,'Cdata',zc,'facecolor','none','edgecolor','flat','linewidth',linewidth)];
+			end
 			hold on      
 		else
 			dist = 5000;
-			if (max(xc)-min(xc)+max(yc)-min(yc))<dist, continue; end
-			h=patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'facecolor','none','edgecolor',color,'linewidth',linewidth);
-			c = horzcat([level, xc'; length(xc), yc']);
+			if isplanet,
+				if (max(xc)-min(xc)+max(yc)-min(yc)+max(zc)-min(zc))<dist, continue; end
+				h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'Zdata',[zc;NaN],'facecolor','none','linewidth',linewidth)];
+				c = horzcat([level, xc'; length(xc), yc'; length(xc), zc']);
+			else
+				if (max(xc)-min(xc)+max(yc)-min(yc))<dist, continue; end
+				h=patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'facecolor','none','edgecolor',color,'linewidth',linewidth);
+				c = horzcat([level, xc'; length(xc), yc']);
+			end
 			%clabel(c,h,'FontSize',10,'labelspacing',20000,'color',color);
 			hold on
@@ -247,5 +269,5 @@
 	else
 		%clabel(c,h,'color',color,'FontSize',10,'labelspacing',20000);
-		clabel(c,h,'FontSize',10,'labelspacing',20000);
-	end
-end
+		%clabel(c,h,'FontSize',10,'labelspacing',20000);
+	end
+end
