Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 16894)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 16895)
@@ -236,5 +236,5 @@
 
 	/*Enthalpy diffusion parameter*/
-	IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element); _assert_(kappa>0.);
+	IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyEnum); _assert_(kappa>0.);
 
 	/* Start  looping on the number of gaussian points: */
@@ -304,5 +304,5 @@
 				for(int j=0;j<numnodes;j++){
 					Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
-					  ((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i])*((u-um)*dbasis[0*3+j]+(v-vm)*dbasis[1*3+j]+(w-wm)*dbasis[2*3+j]);
+					  ((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
 				}
 			}
@@ -310,5 +310,5 @@
 				for(int i=0;i<numnodes;i++){
 					for(int j=0;j<numnodes;j++){
-						Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i]);
+						Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
 					}
 				}
@@ -410,6 +410,4 @@
 	IssmDouble*    basis          = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis         = xNew<IssmDouble>(3*numnodes);
-	IssmDouble*    pressure       = xNew<IssmDouble>(numvertices);
-	IssmDouble*    enthalpypicard = xNew<IssmDouble>(numvertices);
 
 	/*Retrieve all inputs and parameters*/
@@ -425,11 +423,10 @@
 	if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
 	if(stabilization==2){
-		element->GetInputListOnVertices(enthalpypicard,EnthalpyPicardEnum);
-		element->GetInputListOnVertices(pressure,PressureEnum);
 		diameter=element->MinEdgeLength(xyz_list);
+		kappa=this->EnthalpyDiffusionParameterVolume(element,EnthalpyPicardEnum); _assert_(kappa>0.);
 	}
 
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(3);
+	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
@@ -440,5 +437,5 @@
 
 		scalar_def=phi/rho_ice*Jdet*gauss->weight;
-		if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
+		if(dt!=0.) scalar_def=scalar_def*dt;
 
 		/*TODO: add -beta*laplace T_m(p)*/
@@ -458,10 +455,10 @@
 			vy_input->GetInputValue(&v,gauss);
 			vz_input->GetInputValue(&w,gauss);
-			kappa          = element->EnthalpyDiffusionParameterVolume(numvertices,enthalpypicard,pressure) / rho_ice;
-			tau_parameter  = element->StabilizationParameter(u,v,w,diameter,kappa);
-
-			for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
-			if(reCast<bool,IssmDouble>(dt)){
-				for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
+			tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
+
+			for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
+
+			if(dt!=0.){
+				for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
 			}
 		}
@@ -471,6 +468,4 @@
 	xDelete<IssmDouble>(basis);
 	xDelete<IssmDouble>(dbasis);
-	xDelete<IssmDouble>(pressure);
-	xDelete<IssmDouble>(enthalpypicard);
 	xDelete<IssmDouble>(xyz_list);
 	delete gauss;
@@ -798,5 +793,5 @@
 	}
 }/*}}}*/
-IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element){/*{{{*/
+IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
 
 	int         iv;
@@ -813,5 +808,5 @@
 	IssmDouble* dHpmp       = xNew<IssmDouble>(numvertices);
 	element->GetInputListOnVertices(pressures,PressureEnum);
-	element->GetInputListOnVertices(enthalpies,EnthalpyEnum);
+	element->GetInputListOnVertices(enthalpies,enthalpy_enum);
 	for(iv=0;iv<numvertices;iv++){
 		PIE[iv]   = PureIceEnthalpy(element,pressures[iv]);
Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 16894)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 16895)
@@ -36,5 +36,5 @@
 		/*Intermediaries*/
 		IssmDouble EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure);
-		IssmDouble EnthalpyDiffusionParameterVolume(Element* element);
+		IssmDouble EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum);
 		IssmDouble PureIceEnthalpy(Element* element,IssmDouble pressure);
 		IssmDouble TMeltingPoint(Element* element,IssmDouble pressure);