Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16386)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16387)
@@ -2182,4 +2182,23 @@
 }
 /*}}}*/
+/*FUNCTION Tria::EdgeOnBedIndex{{{*/
+int Tria::EdgeOnBedIndex(void){
+
+	IssmDouble values[NUMVERTICES];
+	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
+
+	for(int i=0;i<3;i++){
+		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
+			return i;
+		}
+	}
+
+	_printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
+	_error_("Could not find 2 vertices on bed");
+}
+/*}}}*/
 /*FUNCTION Tria::IsFloating {{{*/
 bool   Tria::IsFloating(){
@@ -2408,4 +2427,48 @@
 
 	/*clean-up*/
+	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::ResetCoordinateSystem{{{*/
+void  Tria::ResetCoordinateSystem(void){
+
+	int        approximation;
+	int        numindices;
+	int       *indices = NULL;
+	IssmDouble slope;
+	IssmDouble xz_plane[6];
+
+	/*For FS only: we want the CS to be tangential to the bedrock*/
+	inputs->GetInputValue(&approximation,ApproximationEnum);
+	if(IsFloating() || !HasEdgeOnBed() || approximation!=FSApproximationEnum) return;
+
+	/*Get number of nodes for velocity only and base*/
+	int index = this->EdgeOnBedIndex();
+	NodeOnEdgeIndices(&numindices,&indices,index,this->VelocityInterpolation());
+
+	/*Get inputs*/
+	Input* slope_input=inputs->GetInput(BedSlopeXEnum); _assert_(slope_input);
+
+	/*Loop over basal nodes and update their CS*/
+	GaussTria* gauss = new GaussTria();
+	for(int i=0;i<numindices;i++){//FIXME
+
+
+		gauss->GaussNode(this->VelocityInterpolation(),indices[i]);
+
+		slope_input->GetInputValue(&slope,gauss);
+
+		IssmDouble theta = atan(slope);
+
+		/*New X axis                  New Z axis*/
+		xz_plane[0]=cos(theta);       xz_plane[3]=0.;  
+		xz_plane[1]=sin(theta);       xz_plane[4]=0.;  
+		xz_plane[2]=0.;               xz_plane[5]=1.;          
+
+		XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]);
+	}
+
+	/*cleanup*/
+	xDelete<int>(indices);
 	delete gauss;
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16386)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16387)
@@ -87,4 +87,5 @@
 		bool        HasEdgeOnBed();
 		void        EdgeOnBedIndices(int* pindex1,int* pindex);
+		int         EdgeOnBedIndex();
 		bool        IsFloating();
 		bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
@@ -106,5 +107,5 @@
 		void        PatchFill(int* pcount, Patch* patch);
 		void        PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
-		void        ResetCoordinateSystem(void){_error_("not implemented yet");};
+		void        ResetCoordinateSystem(void);
 		void	      SmbGradients();
 		IssmDouble  SurfaceArea(void);
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 16386)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 16387)
@@ -1038,2 +1038,62 @@
 }
 /*}}}*/
+/*FUNCTION TriaRef::NodeOnEdgeIndices{{{*/
+void TriaRef::NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement){
+
+	/*Output*/
+	int  numindices;
+	int* indices = NULL;
+
+	switch(finiteelement){
+		case P1Enum: case P1DGEnum: case P1bubbleEnum: case P1bubblecondensedEnum:
+			numindices = 2;
+			indices    = xNew<int>(numindices);
+			switch(index){
+				case 0:
+					indices[0] = 1;
+					indices[1] = 2;
+					break;
+				case 1:
+					indices[0] = 2;
+					indices[1] = 0;
+					break;
+				case 2:
+					indices[0] = 0;
+					indices[1] = 1;
+					break;
+				default:
+					_error_("Edge index provided ("<<index<<") is not between 0 and 2");
+			}
+			break;
+		case P2Enum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			switch(index){
+				case 0:
+					indices[0] = 1;
+					indices[1] = 2;
+					indices[2] = 5;
+					break;
+				case 1:
+					indices[0] = 2;
+					indices[1] = 0;
+					indices[2] = 3;
+					break;
+				case 2:
+					indices[0] = 0;
+					indices[1] = 1;
+					indices[2] = 4;
+					break;
+				default:
+					_error_("Edge index provided ("<<index<<") is not between 0 and 2");
+			}
+			break;
+		default:
+			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+
+	/*Assign output pointer*/
+	*pnumindices = numindices;
+	*pindices    = indices;
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 16386)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 16387)
@@ -55,4 +55,5 @@
 		void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
 
+		void NodeOnEdgeIndices(int* pnumindices,int** pindices,int index,int finiteelement);
 		int  NumberofNodes(void);
 		int  NumberofNodes(int finiteelement);