Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16620)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16621)
@@ -3776,5 +3776,5 @@
 	IssmDouble epsvel=2.220446049250313e-16;
 	IssmDouble heatcapacity,thermalconductivity,dt;
-	IssmDouble pressure[NUMVERTICES],enthalpy[NUMVERTICES];
+	IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
 	IssmDouble latentheat,kappa;
 	IssmDouble tau_parameter,diameter;
@@ -3824,10 +3824,7 @@
 		/*Conduction: */  
 		GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
-
-		//enthalpy_input->GetInputValue(&enthalpy, gauss);
-		//pressure_input->GetInputValue(&pressure, gauss);
-		GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum);
+		GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
 		GetInputListOnVertices(&pressure[0],PressureEnum);
-		kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);  // TODO: switch to pointers
+		kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);
 
 		D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
@@ -4232,7 +4229,6 @@
 	IssmDouble rho_ice,heatcapacity;
 	IssmDouble thermalconductivity,kappa;
-	IssmDouble viscosity,pressure;
-	IssmDouble enthalpy;
-	IssmDouble enthalpy_v[NUMVERTICES], pressure_v[NUMVERTICES];
+	IssmDouble viscosity,pressure,enthalpy;
+	IssmDouble enthalpypicard[NUMVERTICES], pressure_p[NUMVERTICES];
 	IssmDouble tau_parameter,diameter;
 	IssmDouble u,v,w;
@@ -4300,9 +4296,7 @@
 			vy_input->GetInputValue(&v, gauss);
 			vz_input->GetInputValue(&w, gauss);
-			//pressure_input->GetInputValue(&pressure, gauss);
-			//enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
-			GetInputListOnVertices(&enthalpy_v[0],EnthalpyPicardEnum);
-			GetInputListOnVertices(&pressure_v[0],PressureEnum);
-			kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy_v,pressure_v);
+			GetInputListOnVertices(&enthalpypicard[0],EnthalpyPicardEnum);
+			GetInputListOnVertices(&pressure_p[0],PressureEnum);
+			kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpypicard,pressure_p);
 			kappa/=rho_ice;
 			tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa); 
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16620)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16621)
@@ -397,7 +397,7 @@
 
 	IssmDouble lambda; // fraction of cold ice
-	IssmDouble kappa,kappa_c, kappa_t;
-	IssmDouble Hc, Ht;
-	IssmDouble PIE[6], dHpmp[6];
+	IssmDouble kappa,kappa_c,kappa_t;  //enthalpy conductivities
+	IssmDouble Hc,Ht;
+	IssmDouble PIE[6],dHpmp[6];
 	int iv;
 
@@ -419,5 +419,6 @@
 	}
 	else {
-		/* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice*/
+		/* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
+		 cf Patankar 1980, pp44 */
 		kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
 		kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);