Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6267)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6268)
@@ -2918,5 +2918,5 @@
 	double     heatcapacity,thermalconductivity,dt;
 	double     scalar_artdiff;
-	double     tau_parameter,diameter,normu;
+	double     tau_parameter,diameter;
 	double     xyz_list[NUMVERTICES][3];
 	double     B[3][numdof];
@@ -3038,11 +3038,8 @@
 			GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
 
-			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));
-			}
-			else tau_parameter=diameter/(2*normu);
+			tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
+			//if(dt) scalar_artdiff=scalar_artdiff*dt;
+
 			scalar_artdiff=tau_parameter*gauss->weight*Jdet;
-
 			for(i=0;i<numdof;i++){
 				for(j=0;j<numdof;j++){
@@ -3830,5 +3827,5 @@
 	double thermalconductivity;
 	double viscosity,temperature;
-	double tau_parameter,diameter,normu;
+	double tau_parameter,diameter;
 	double u,v,w,scalar_artdiff;
 	double scalar_def,scalar_transient;
@@ -3890,11 +3887,8 @@
 			vz_input->GetParameterValue(&w, gauss);
 
-			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));
-			}
-			else tau_parameter=diameter/(2*normu);
-
+			tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
 			scalar_artdiff=tau_parameter*gauss->weight*Jdet*phi/(rho_ice*heatcapacity);
+			//if(dt) scalar_artdiff=scalar_artdiff*dt;
+
 			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
 		}
@@ -4296,4 +4290,20 @@
 	delete gauss;
 	xfree((void**)&doflist);
+}
+/*}}}*/
+/*FUNCTION Penta::GetStabilizationParameter {{{1*/
+double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
+	/*Compute stabilization parameter*/
+
+	double normu;
+	double tau_parameter;
+
+	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));
+	}
+	else tau_parameter=diameter/(2*normu);
+
+	return tau_parameter;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 6267)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 6268)
@@ -196,4 +196,5 @@
 		void	  GetSolutionFromInputsDiagnosticVert(Vec solutiong);
 		void	  GetSolutionFromInputsThermal(Vec solutiong);
+		double GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity);
 		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);