Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15871)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15872)
@@ -5026,5 +5026,5 @@
 
 	/*Fetch indices of basal nodes for this finite element*/
-	BasalNodeIndices(&numindices,&indices,this->VelocityInterpolation());
+	BasalNodeIndices(&numindices,&indices,this->element_type);
 
 	/*Get parameters and inputs: */
@@ -5062,5 +5062,5 @@
 	/* Intermediaries */
 	bool        isenthalpy;
-	int         i,j,analysis_type;
+	int         i,j,analysis_type,numindices,numindicesup;
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
@@ -5077,7 +5077,7 @@
 	IssmDouble  vx,vy,vz;
 	IssmDouble  dt;
-	Friction   *friction                 = NULL;
-	GaussPenta *gauss                    = NULL;
-	GaussPenta *gaussup                  = NULL;
+	int        *indices   = NULL;
+	int        *indicesup = NULL;
+	Friction   *friction  = NULL;
 
 	/* Only compute melt rates at the base of grounded ice*/
@@ -5104,10 +5104,15 @@
 	friction=new Friction("3d",inputs,matpar,analysis_type);
 
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussPenta(0,1,2,2);
-	gaussup=new GaussPenta(3,4,5,2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-		gaussup->GaussPoint(ig);
+	/*Fetch indices of basal and surface nodes for this finite element*/
+	BasalNodeIndices(&numindices,&indices,this->element_type);
+	SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type);
+	_assert_(numindices==numindicesup); 
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	GaussPenta* gauss=new GaussPenta();
+	GaussPenta* gaussup=new GaussPenta();
+	for(i=0;i<numindices;i++){
+		gauss->GaussNode(this->element_type,indices[i]);
+		gaussup->GaussNode(this->element_type,indicesup[i]);
 
 		watercolumn_input->GetInputValue(&watercolumn, gauss);
@@ -5161,4 +5166,10 @@
 		// TODO: feed meltrate & watercolumn back to model
 	}  
+
+	/*Clean up and return*/
+	xDelete<int>(indices);
+	xDelete<int>(indicesup);
+	delete gauss;
+	delete gaussup;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15871)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15872)
@@ -2118,2 +2118,60 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::SurfaceNodeIndices{{{*/
+void PentaRef::SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement){
+
+	/*Output*/
+	int  numindices;
+	int* indices = NULL;
+
+	switch(finiteelement){
+		case P1Enum: case P1DGEnum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			break;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			break;
+		case P2xP1Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			indices[3] = 9;
+			indices[4] = 10;
+			indices[5] = 11;
+			break;
+		case P1xP2Enum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			return;
+		case P2Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 3;
+			indices[1] = 4;
+			indices[2] = 5;
+			indices[3] = 12;
+			indices[4] = 13;
+			indices[5] = 14;
+			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/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15871)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15872)
@@ -70,4 +70,5 @@
 
 		void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
+		void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
 		int  NumberofNodes(void);
 		int  NumberofNodesVelocity(void);