Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16613)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16614)
@@ -3776,5 +3776,5 @@
 	IssmDouble epsvel=2.220446049250313e-16;
 	IssmDouble heatcapacity,thermalconductivity,dt;
-	IssmDouble pressure,enthalpy;
+	IssmDouble pressure[NUMVERTICES],enthalpy[NUMVERTICES];
 	IssmDouble latentheat,kappa;
 	IssmDouble tau_parameter,diameter;
@@ -3815,5 +3815,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussPenta(3,3);
+	gauss=new GaussPenta(2,2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -3825,7 +3825,9 @@
 		GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
 
-		enthalpy_input->GetInputValue(&enthalpy, gauss);
-		pressure_input->GetInputValue(&pressure, gauss);
-		kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); _assert_(kappa>0.);
+		//enthalpy_input->GetInputValue(&enthalpy, gauss);
+		//pressure_input->GetInputValue(&pressure, gauss);
+		GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum);
+		GetInputListOnVertices(&pressure[0],PressureEnum);
+		kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);  // TODO: switch to pointers
 
 		D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
@@ -4231,5 +4233,6 @@
 	IssmDouble thermalconductivity,kappa;
 	IssmDouble viscosity,pressure;
-	IssmDouble enthalpy,enthalpypicard;
+	IssmDouble enthalpy;
+	IssmDouble enthalpy_v[NUMVERTICES], pressure_v[NUMVERTICES];
 	IssmDouble tau_parameter,diameter;
 	IssmDouble u,v,w;
@@ -4297,7 +4300,9 @@
 			vy_input->GetInputValue(&v, gauss);
 			vz_input->GetInputValue(&w, gauss);
-			pressure_input->GetInputValue(&pressure, gauss);
-			enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
-			kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
+			//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);
 			kappa/=rho_ice;
 			tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa); 
@@ -5021,5 +5026,5 @@
 			else{
 				enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
-				kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy[iv],enthalpy[iv+NUMVERTICES2D], pressure[iv],pressure[iv+NUMVERTICES2D]); _assert_(kappa>0.);
+				kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy, pressure); _assert_(kappa>0.);
 				for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
 			}
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16613)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16614)
@@ -387,46 +387,50 @@
 /*FUNCTION Matpar::GetEnthalpyDiffusionParameter{{{*/
 IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
-	/*if (enthalpy<=PureIceEnthalpy(pressure))
+	if (enthalpy<PureIceEnthalpy(pressure))
 		return thermalconductivity/heatcapacity;
 	else
-		return thermalconductivity/heatcapacity*pow(10.,-5);
-	*/
-
-	IssmDouble eps=0.1*heatcapacity;
-	IssmDouble hpmp=PureIceEnthalpy(pressure);
-	IssmDouble kappa_c=thermalconductivity/heatcapacity;
-	IssmDouble kappa_t=temperateiceconductivity/heatcapacity;
-	
-	if(enthalpy<=hpmp-eps)
-		return kappa_c;
-	else if(enthalpy>=hpmp+eps)
-		return kappa_t;
+		return temperateiceconductivity/heatcapacity;
+}
+/*}}}*/
+/*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/
+IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]){
+
+	IssmDouble lambda; // fraction of cold ice
+	IssmDouble kappa,kappa_c, kappa_t;
+	IssmDouble Hc, Ht;
+	IssmDouble PIE[6], dHpmp[6];
+	int iv;
+
+	for (iv=0; iv<6; iv++){
+		PIE[iv]=PureIceEnthalpy(pressure[iv]);
+		dHpmp[iv]=enthalpy[iv]-PIE[iv];
+	}
+
+	bool allequalsign=true;
+	if(dHpmp[0]<0)
+		for(iv=1; iv<6;iv++)
+			allequalsign=(allequalsign && (dHpmp[iv]<0));
+	else
+		for(iv=1; iv<6;iv++)
+			allequalsign=(allequalsign && (dHpmp[iv]>=0));
+
+	if(allequalsign){
+		kappa=GetEnthalpyDiffusionParameter(enthalpy[0], pressure[0]);
+	}
 	else {
-		IssmDouble xi=enthalpy-hpmp;
-		IssmDouble pi=3.141592653589793238462643;
-		return kappa_c + (kappa_t-kappa_c)*((xi+eps)/(2*eps) + sin(pi*xi/eps)/(2*pi));
-	}
-
-	//return 1./(1.+exp(-(enthalpy-(hpmp))/eps))*(kappa_t-kappa_c) + kappa_c;
-}
-/*}}}*/
-/*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/
-IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy0, IssmDouble enthalpy1, IssmDouble pressure0, IssmDouble pressure1){
-	/*returns kappa depending on distribution of enthalpy over edge of element
-		lambda is the barycentric coordinate that solves H0+(H1-H0)*lambda=H_pureice.
-		it represents the fraction of the ice column which is temperate/cold like H0.
-		if lambda<0 or lambda>1, then the whole ice column is cold or temperate, respectively.
-	*/
-	IssmDouble kappa, kappa0, kappa1;
-	if (enthalpy0!=enthalpy1){
-		IssmDouble lambda=(PureIceEnthalpy(pressure0)-enthalpy0)/(enthalpy1-enthalpy0);
-		if ((lambda>=0.) && (lambda<=1.)){
-			kappa0=GetEnthalpyDiffusionParameter(enthalpy0,pressure0);
-			kappa1=GetEnthalpyDiffusionParameter(enthalpy1,pressure1);
-			kappa=lambda*kappa0+(1.-lambda)*kappa1;
+		/* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice*/
+		kappa_c=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
+		kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
+		Hc=0.; Ht=0.;
+		for(iv=0; iv<6;iv++){
+			if(enthalpy[iv]<PIE[iv])
+				Hc+=(PIE[iv]-enthalpy[iv]);
+			else
+				Ht+=(enthalpy[iv]-PIE[iv]);
 		}
-	}
-	else
-		kappa=GetEnthalpyDiffusionParameter(enthalpy0, pressure0);
+		_assert_((Hc+Ht)>0.);
+		lambda=Hc/(Hc+Ht);
+		kappa=1./(lambda/kappa_c + (1.-lambda)/kappa_t);
+	}
 	return kappa;
 }
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 16613)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 16614)
@@ -128,5 +128,5 @@
 		IssmDouble PureIceEnthalpy(IssmDouble pressure);
 		IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
-		IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy0, IssmDouble enthalpy1, IssmDouble pressure0, IssmDouble pressure1);
+		IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]);
 		IssmDouble GetLithosphereShearModulus();
 		IssmDouble GetLithosphereDensity();
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 16613)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 16614)
@@ -181,5 +181,5 @@
 	MaterialsThermalExchangeVelocityEnum,
 	MaterialsThermalconductivityEnum,
-  MaterialsTemperateiceconductivityEnum,
+  	MaterialsTemperateiceconductivityEnum,
 	MaterialsLithosphereShearModulusEnum,
 	MaterialsLithosphereDensityEnum,