Index: /issm/trunk-jpl/src/m/exp/contourlevelzero.m
===================================================================
--- /issm/trunk-jpl/src/m/exp/contourlevelzero.m	(revision 22930)
+++ /issm/trunk-jpl/src/m/exp/contourlevelzero.m	(revision 22931)
@@ -16,4 +16,9 @@
 x=md.mesh.x;
 y=md.mesh.y;
+if isfield(md.mesh,'z'),
+	z=md.mesh.z;
+else
+	z=zeros(md.mesh.numberofvertices,1);
+end
 index=md.mesh.elements;
 
@@ -82,4 +87,7 @@
 y1=zeros(numelems,1);
 y2=zeros(numelems,1);
+z1=zeros(numelems,1);
+z2=zeros(numelems,1);
+
 edge_l=zeros(numelems,2);
 
@@ -96,4 +104,7 @@
 		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));
@@ -105,4 +116,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));
@@ -114,4 +128,7 @@
 		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));
@@ -130,5 +147,6 @@
 	%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
@@ -136,4 +154,5 @@
 	x1(1)=[]; x2(1)=[];
 	y1(1)=[]; y2(1)=[];
+	z1(1)=[]; z2(1)=[];
 
 	[ro1,co1]=find(edge_l==e1);
@@ -142,5 +161,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:
@@ -148,5 +167,5 @@
 
 		else
-			xc=[x1(ro1);xc]; yc=[y1(ro1);yc];
+			xc=[x1(ro1);xc]; yc=[y1(ro1);yc];zc=[z1(ro1);zc];
 
 			%next edge:
@@ -158,4 +177,5 @@
 		x1(ro1)=[]; x2(ro1)=[];
 		y1(ro1)=[]; y2(ro1)=[];
+		z1(ro1)=[]; z2(ro1)=[];
 
 		%next connection
@@ -169,10 +189,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:
@@ -184,4 +204,5 @@
 		x1(ro2)=[]; x2(ro2)=[];
 		y1(ro2)=[]; y2(ro2)=[];
+		z1(ro2)=[]; z2(ro2)=[];
 
 		%next connection
@@ -192,4 +213,5 @@
 	contours(end+1).x=xc;
 	contours(end).y=yc;
+	contours(end).z=zc;
 	contours(end).name='';
 	contours(end).nods=length(xc);