Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11354)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11355)
@@ -3854,8 +3854,10 @@
 	double     rho_ice,heatcapacity,geothermalflux_value;
 	double     basalfriction,alpha2,vx,vy;
-	double     scalar;
+	double     scalar,enthalpy,enthalpyup;
+	double     pressure,pressureup;
 	double     basis[NUMVERTICES];
 	Friction*  friction=NULL;
 	GaussPenta* gauss=NULL;
+	GaussPenta* gaussup=NULL;
 
 	/* Geothermal flux on ice sheet base and basal friction */
@@ -3875,4 +3877,6 @@
 	Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
 	Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
+	Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
+	Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
 	Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
 
@@ -3882,25 +3886,43 @@
 	/* Start looping on the number of gauss 2d (nodes on the bedrock) */
 	gauss=new GaussPenta(0,1,2,2);
+	gaussup=new GaussPenta(3,4,5,2);
 	for(ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
+		gaussup->GaussPoint(ig);
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
 		GetNodalFunctionsP1(&basis[0], gauss);
 
-		geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
-		friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
-		
-		scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
-		if(dt) scalar=dt*scalar;
-
-		for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
+		enthalpy_input->GetInputValue(&enthalpy,gauss);
+		pressure_input->GetInputValue(&pressure,gauss);
+		if(enthalpy>matpar->PureIceEnthalpy(pressure)){
+			enthalpy_input->GetInputValue(&enthalpyup,gaussup);
+			pressure_input->GetInputValue(&pressureup,gaussup);
+			if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
+				//do nothing, don't add heatflux
+			}
+			else{
+				//need to change spcenthalpy according to Aschwanden 
+				//NEED TO UPDATE
+			}
+		}
+		else{
+			geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
+			friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
+			vx_input->GetInputValue(&vx,gauss);
+			vy_input->GetInputValue(&vy,gauss);
+			basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
+
+			scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
+			if(dt) scalar=dt*scalar;
+
+			for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
+		}
 	}
 
 	/*Clean up and return*/
 	delete gauss;
+	delete gaussup;
 	delete friction;
 	return pe;
@@ -4085,10 +4107,8 @@
 	double     rho_ice,heatcapacity,geothermalflux_value;
 	double     basalfriction,alpha2,vx,vy;
-	double     scalar,enthalpy,enthalpyup;
-	double     pressure,pressureup;
 	double     basis[NUMVERTICES];
+	double     scalar;
 	Friction*  friction=NULL;
 	GaussPenta* gauss=NULL;
-	GaussPenta* gaussup=NULL;
 
 	/* Geothermal flux on ice sheet base and basal friction */
@@ -4108,6 +4128,4 @@
 	Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
 	Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
-	Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
-	Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
 	Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
 
@@ -4117,27 +4135,11 @@
 	/* Start looping on the number of gauss 2d (nodes on the bedrock) */
 	gauss=new GaussPenta(0,1,2,2);
-	gaussup=new GaussPenta(3,4,5,2);
 	for(ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
-		gaussup->GaussPoint(ig);
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
 		GetNodalFunctionsP1(&basis[0], gauss);
 
-		enthalpy_input->GetInputValue(&enthalpy,gauss);
-		pressure_input->GetInputValue(&pressure,gauss);
-//		if(enthalpy>matpar->PureIceEnthalpy(pressure)){
-//			enthalpy_input->GetInputValue(&enthalpyup,gaussup);
-//			pressure_input->GetInputValue(&pressureup,gaussup);
-//			if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
-//				//do nothing, don't add heatflux
-//			}
-//			else{
-//				//need to change spcenthalpy according to Aschwanden 
-//				//NEED TO UPDATE
-//			}
-//		}
-//		else{
 			geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
 			friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
@@ -4150,10 +4152,8 @@
 
 			for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
-//		}
 	}
 
 	/*Clean up and return*/
 	delete gauss;
-	delete gaussup;
 	delete friction;
 	return pe;