Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15895)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15896)
@@ -5236,5 +5236,5 @@
 	/* Intermediaries */
 	bool        isenthalpy;
-	int         i,j,ig,analysis_type,numindices,numindicesup;
+	int         i,j,analysis_type;
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
@@ -5243,4 +5243,5 @@
 	IssmDouble  normal_base[3], d1enthalpy[3];
 	IssmDouble  kappa;
+    IssmDouble  meltrate[3], watercolumn[3];
 	IssmDouble  enthalpy, enthalpyup;
 	IssmDouble  pressure, pressureup;
@@ -5250,6 +5251,4 @@
 	IssmDouble  vx,vy,vz;
 	IssmDouble  dt;
-	int        *indices   = NULL;
-	int        *indicesup = NULL;
 	Friction   *friction  = NULL;
 
@@ -5277,32 +5276,25 @@
 	friction=new Friction("3d",inputs,matpar,analysis_type);
 
-	/*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); 
-    
-	IssmDouble  meltrate[numindices], watercolumn[numindices];
-
 	/*Ok, get meltrates now from basal conditions*/
 	GaussPenta* gauss=new GaussPenta();
 	GaussPenta* gaussup=new GaussPenta();
-	for(ig=0;ig<numindices;ig++){
-		gauss->GaussVertex(this->element_type,indices[ig]);
-		gaussup->GaussNode(this->element_type,indicesup[ig]);
+	for(int iv=0;iv<3;iv++){
+		gauss->GaussVertex(iv);
+		gaussup->GaussNode(iv+3);
 
         // TODO: make sure that no node is computed twice (insert mask)
 
-		watercolumn_input->GetInputValue(&watercolumn[indices[ig]], gauss);
+		watercolumn_input->GetInputValue(&watercolumn[iv], gauss);
 		enthalpy_input->GetInputValue(&enthalpy, gauss);
 		pressure_input->GetInputValue(&pressure, gauss);
 
 		/*Calculate basal meltrate*/
-		if((watercolumn[indices[ig]]>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
+		if((watercolumn[iv]>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
 			enthalpy=matpar->PureIceEnthalpy(pressure);
 		}
 		else if(enthalpy<matpar->PureIceEnthalpy(pressure)){
-			meltrate[indices[ig]]=0.;   //TODO: set zero meltrate and watercolumn in model
-			watercolumn[indices[ig]]=0.;
-			return;
+			meltrate[iv]=0.;   
+			watercolumn[iv]=0.;
+			continue;
 		}
 
@@ -5333,12 +5325,12 @@
 		vz_input->GetInputValue(&vz,gauss);
 		basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
-		meltrate[indices[ig]]=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
+		meltrate[iv]=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
 
 		/*Update water column*/
 		this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 		if(reCast<bool,IssmDouble>(dt))
-            watercolumn[indices[ig]]+=dt*meltrate[indices[ig]];
+            watercolumn[iv]+=dt*meltrate[iv];
 		else
-            watercolumn[indices[ig]]=meltrate[indices[ig]];
+            watercolumn[iv]=meltrate[iv];
 	}  
     /*update meltrate and watercolumn*/
@@ -5347,6 +5339,4 @@
 
 	/*Clean up and return*/
-	xDelete<int>(indices);
-	xDelete<int>(indicesup);
 	delete gauss;
 	delete gaussup;
@@ -5359,5 +5349,4 @@
 
     /*Intermediaries*/
-    int ig;
     bool isenthalpy;
     IssmDouble waterfraction_array[NUMVERTICES], temperature[NUMVERTICES];
@@ -5376,17 +5365,17 @@
     this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     latentheat=matpar->GetLatentHeat();
-    gauss=new GaussPenta(2,3);
-    for(ig=0;ig<NUMVERTICES;ig++){ 
-        gauss->GaussVertex(ig);
+    gauss=new GaussPenta();
+    for(int iv=0;iv<NUMVERTICES;iv++){ 
+        gauss->GaussVertex(iv);
         /*TODO: Check whether Vertex has been drained already*/
-        enthalpy_input->GetInputValue(&enthalpy[ig], gauss);
-        pressure_input->GetInputValue(&pressure[ig], gauss);
-        matpar->EnthalpyToThermal(&temperature[ig], &waterfraction[ig], enthalpy[ig],pressure[ig]); 
+        enthalpy_input->GetInputValue(&enthalpy[iv], gauss);
+        pressure_input->GetInputValue(&pressure[iv], gauss);
+        matpar->EnthalpyToThermal(&temperature[iv], &waterfraction[iv], enthalpy[iv],pressure[iv]); 
     
         /*drain water fraction*/
-        waterfraction[ig]-=dt*DrainageFunctionWaterfraction(waterfraction[ig]);
-        if(waterfraction[ig]<0) waterfraction[ig]=0.;
+        waterfraction[iv]-=dt*DrainageFunctionWaterfraction(waterfraction[iv]);
+        if(waterfraction[iv]<0.) waterfraction[iv]=0.;
         /*update enthalpy*/
-        matpar->ThermalToEnthalpy(&enthalpy[ig], temperature[ig], waterfraction[ig], pressure[ig]);        
+        matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);        
     }
     /*feed updated results back into model*/