Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 22988)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 22989)
@@ -319,5 +319,4 @@
 		virtual void       ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
 		virtual void       WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue){_error_("not implemented yet");};
-		virtual void       WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");};
 		virtual void       ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 22988)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 22989)
@@ -4066,82 +4066,55 @@
 	if(!input) _error_("Cannot calculate distance to "<<EnumToStringx(fieldenum)<<", input not found");
 	IssmDouble minvalue = input->Min();
+	if(minvalue>fieldvalue) return;
 	IssmDouble maxvalue = input->Max();
-	if(minvalue>fieldvalue || maxvalue<fieldvalue) return;
-
-	/* check #2: If only one vertex is on fieldvalue, there is no segment here */
-	IssmDouble lsf[NUMVERTICES];
-	this->GetInputListOnVertices(&lsf[0],fieldenum);
-	int nrice=0;       
-	for(int i=0;i<NUMVERTICES;i++) if(lsf[i]==fieldvalue) nrice++;
-	if(nrice==1) return;
-
-	/*2. Write segments*/
-	IssmDouble* xyz_list_zero = NULL;
+	if(maxvalue<fieldvalue) return;
+
+	/*2. Find coordinates of where levelset crosses 0*/
+	int         numiceverts;
+	IssmDouble  s[2],x[2],y[2];
+	int        *indices = NULL;
+	this->GetLevelsetIntersection(&indices, &numiceverts,&s[0],fieldenum,fieldvalue);
+	_assert_(numiceverts); 
+
+	/*3 Write coordinates*/
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
-	this->ZeroLevelsetCoordinates(&xyz_list_zero,&xyz_list[0][0],fieldenum);
-	if(xyz_list_zero){
-		IssmDouble x[2],y[2];
-		x[0] = xyz_list_zero[0*3 + 0]; x[1] = xyz_list_zero[1*3 + 0];
-		y[0] = xyz_list_zero[0*3 + 1]; y[1] = xyz_list_zero[1*3 + 1];
-		segments->AddObject(new Contour<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));
-	}
-	xDelete<IssmDouble>(xyz_list_zero);
-}
-/*}}}*/
-void       Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/
-
-	if(!this->IsZeroLevelset(MaskIceLevelsetEnum)) return;
-
-	IssmDouble* xyz_list_zero = NULL;
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
-	this->ZeroLevelsetCoordinates(&xyz_list_zero,&xyz_list[0][0], MaskIceLevelsetEnum);
-	if(xyz_list_zero){
-		IssmDouble x[2];
-		IssmDouble y[2];
-		x[0] = xyz_list_zero[0*3 + 0]; x[1] = xyz_list_zero[1*3 + 0];
-		y[0] = xyz_list_zero[0*3 + 1]; y[1] = xyz_list_zero[1*3 + 1];
-		segments->AddObject(new Contour<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));
-	}
-	xDelete<IssmDouble>(xyz_list_zero);
-
-//	IssmDouble ls[NUMVERTICES];
-//	IssmDouble  xyz_list[NUMVERTICES][3];
-//
-//	if(IsIceInElement()){
-//
-//		/*Retrieve all inputs and parameters*/
-//		GetInputListOnVertices(&ls[0],levelset_enum);
-//
-//		/*If the level set is awlays <0, there is no ice front here*/
-//		bool iszerols= false;
-//		if(IsIceInElement()){
-//			if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
-//				iszerols = true;
-//			}
-//		}
-//
-//		if(iszerols){
-//			/*OK we have one segment!*/
-//			::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-//			int count = 0;
-//			IssmDouble x[2];
-//			IssmDouble y[2];
-//
-//			for(int i=0;i<NUMVERTICES,i++){
-//				int index1 = i;
-//				int index1 = (i+1)%3;
-//				if(ls[index1]<=0 && ls[index2]>=0){
-//
-//				}
-//
-//			}
-//			Contour* segment = new Contour<IssmDouble>(segment->Size()+1,2,x,y,false);
-//		}
-//
-//	}
-//
-//	_error_("STOP");
+	int counter = 0;
+	if((numiceverts>0) && (numiceverts<NUMVERTICES)){
+		for(int i=0;i<numiceverts;i++){
+			for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
+				x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
+				y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
+				counter++;
+			}
+		}
+	}
+	else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
+		IssmDouble lsf[NUMVERTICES];
+		this->GetInputListOnVertices(&lsf[0],fieldenum);
+		for(int i=0;i<NUMVERTICES;i++){
+			if(lsf[indices[i]]==0.){
+				x[counter]=xyz_list[indices[i]][0];
+				y[counter]=xyz_list[indices[i]][1];
+				counter++;
+			}
+			if(counter==2) break;
+		}
+		if(counter==1){
+			/*We actually have only 1 vertex on levelset, write a single point as a segment*/
+			_error_("TODO");
+			counter==2;
+		}
+	}
+	else{
+		_error_("not sure what's going on here...");
+	}
+
+	/*4. Write segment*/
+	_assert_(counter==2);
+	segments->AddObject(new Contour<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));
+
+	/*Cleanup and return*/
+	xDelete<int>(indices);
 }
 /*}}}*/
@@ -4151,24 +4124,23 @@
 
 	/*Intermediaries*/
-	const int dim=3;
-	int numiceverts;
-	int i, n, e, counter;
-	IssmDouble s[2];
-	int* indices=NULL;
-	IssmDouble* xyz_zero=NULL;
+	int         numiceverts;
+	int         i,n,e,counter;
+	IssmDouble  s[2];
+	int        *indices  = NULL;
+	IssmDouble *xyz_zero = NULL;
 
 	this->GetLevelsetIntersection(&indices, &numiceverts, s,levelsetenum,0.);
 	
-	//TODO: check if for 2 iceverts front segment is oriented in CCW way
-	
-	if(numiceverts>0) xyz_zero=xNew<IssmDouble>(2*dim);
-	if((numiceverts>0)&&(numiceverts<NUMVERTICES)){
+	printf("numiceverts = %i\n",numiceverts);
+	if(numiceverts>0) xyz_zero=xNew<IssmDouble>(2*3);
+
+	if((numiceverts>0) && (numiceverts<NUMVERTICES)){
 		counter=0;
-		for(i=0;i<numiceverts;i++){	// iterate over ice vertices
+		for(i=0;i<numiceverts;i++){
 			for(n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
-				for(e=0;e<dim;e++){ // spatial direction
-					int ind_ice		=dim*indices[i]+e;
-					int ind_noice	=dim*indices[n]+e;
-					int ind			=dim*counter+e;
+				for(e=0;e<3;e++){ // spatial direction
+					int ind_ice   = 3 *indices[i]+e;
+					int ind_noice = 3 *indices[n]+e;
+					int ind       = 3 *counter+e;
 					xyz_zero[ind]=xyz_list[ind_ice]+s[counter]*(xyz_list[ind_noice]-xyz_list[ind_ice]);
 				}
@@ -4183,5 +4155,5 @@
 		for(i=0;i<NUMVERTICES;i++){
 			if(lsf[indices[i]]==0.){
-				for(e=0;e<dim;e++)	xyz_zero[dim*counter+e]=xyz_list[dim*indices[i]+e];
+				for(e=0;e<3;e++)	xyz_zero[3*counter+e]=xyz_list[3*indices[i]+e];
 				counter++;
 			}
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 22988)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 22989)
@@ -144,5 +144,4 @@
 		void        VerticalSegmentIndicesBase(int** pindices,int* pnumseg){_error_("not implemented yet");};
 		void			WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue);
-		void			WriteLevelsetSegment(DataSet* segments);
 		void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22988)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22989)
@@ -2030,5 +2030,5 @@
 		if(!element->IsOnBase()) continue;
 		Element* basalelement = element->SpawnBasalElement();
-		basalelement->WriteLevelsetSegment(segments);
+		basalelement->WriteFieldIsovalueSegment(segments,MaskIceLevelsetEnum,0.);
 		if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	}
Index: /issm/trunk-jpl/src/c/cores/movingfront_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/movingfront_core.cpp	(revision 22988)
+++ /issm/trunk-jpl/src/c/cores/movingfront_core.cpp	(revision 22989)
@@ -85,5 +85,4 @@
 
 	/*Reset levelset if needed*/
-	printf("-------------- file: movingfront_core.cpp line: %i\n",__LINE__); 
 	if(reinit_frequency && (step%reinit_frequency==0)){
 		if(VerboseSolution()) _printf0_("   reinitializing level set\n");