Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11360)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11361)
@@ -1132,6 +1132,8 @@
 /*}}}*/
 /*FUNCTION Penta::GetStabilizationParameter {{{1*/
-double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
+double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double kappa){
 	/*Compute stabilization parameter*/
+	/*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
+	/*kappa=enthalpydiffusionparameter for enthalpy model*/
 
 	double normu;
@@ -1139,6 +1141,6 @@
 
 	normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
-	if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
-		tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
+	if(normu*diameter/(3*2*kappa)<1){ 
+		tau_parameter=pow(diameter,2)/(3*2*2*kappa);
 	}
 	else tau_parameter=diameter/(2*normu);
@@ -3372,5 +3374,5 @@
 		else if(stabilization==2){
 			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
-			tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
+			tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
 
 			for(i=0;i<numdof;i++){
@@ -3486,5 +3488,5 @@
 	double     Jdet,u,v,w,um,vm,wm,vel;
 	double     h,hx,hy,hz,vx,vy,vz;
-	double     gravity,rho_ice,rho_water;
+	double     gravity,rho_ice,rho_water,kappa;
 	double     heatcapacity,thermalconductivity,dt;
 	double     tau_parameter,diameter;
@@ -3512,4 +3514,5 @@
 	heatcapacity=matpar->GetHeatCapacity();
 	thermalconductivity=matpar->GetThermalConductivity();
+	kappa=thermalconductivity/(rho_ice*heatcapacity);
 	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 	this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
@@ -3603,5 +3606,5 @@
 		else if(stabilization==2){
 			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
-			tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
+			tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
 
 			for(i=0;i<numdof;i++){
@@ -3708,6 +3711,6 @@
 	double Jdet,phi,dt;
 	double rho_ice,heatcapacity;
-	double thermalconductivity;
-	double viscosity,enthalpy;
+	double thermalconductivity,kappa;
+	double viscosity,enthalpy,pressure;
 	double tau_parameter,diameter;
 	double u,v,w;
@@ -3733,4 +3736,5 @@
 	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=NULL;
 	if (dt) enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(inputs);
@@ -3746,4 +3750,5 @@
 		GetNodalFunctionsP1(&L[0], gauss);
 
+		enthalpy_input->GetInputValue(&enthalpy, gauss);
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
 		matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
@@ -3757,5 +3762,4 @@
 		/* Build transient now */
 		if(dt){
-			enthalpy_input->GetInputValue(&enthalpy, gauss);
 			scalar_transient=enthalpy*Jdet*gauss->weight;
 			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
@@ -3768,6 +3772,8 @@
 			vy_input->GetInputValue(&v, gauss);
 			vz_input->GetInputValue(&w, gauss);
-
-			tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
+			pressure_input->GetInputValue(&pressure, gauss);
+			kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
+
+			tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
 
 			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
@@ -3961,5 +3967,5 @@
 	double Jdet,phi,dt;
 	double rho_ice,heatcapacity;
-	double thermalconductivity;
+	double thermalconductivity,kappa;
 	double viscosity,temperature;
 	double tau_parameter,diameter;
@@ -3981,4 +3987,5 @@
 	heatcapacity=matpar->GetHeatCapacity();
 	thermalconductivity=matpar->GetThermalConductivity();
+	kappa=heatcapacity/(rho_ice*thermalconductivity);
 	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 	this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
@@ -4022,5 +4029,5 @@
 			vz_input->GetInputValue(&w, gauss);
 
-			tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
+			tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
 
 			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.h	(revision 11360)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h	(revision 11361)
@@ -183,5 +183,5 @@
 		void	  GetPhi(double* phi, double*  epsilon, double viscosity);
 		void	  GetSolutionFromInputsEnthalpy(Vec solutiong);
-		double  GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity);
+		double  GetStabilizationParameter(double u, double v, double w, double diameter, double kappa);
 		void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
 		void    GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);