Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6269)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6270)
@@ -2928,5 +2928,6 @@
 	double     L[numdof];
 	double     dh1dh6[3][6];
-	double     D_scalar;
+	double     D_scalar_conduct,D_scalar_advec;
+	double     D_scalar_trans,D_scalar_artdiff;
 	double     D[3][3];
 	double     K[2][2]={0.0};
@@ -2969,10 +2970,10 @@
 		GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
 
-		D_scalar=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
-		if(dt) D_scalar=D_scalar*dt;
-
-		D[0][0]=D_scalar; D[0][1]=0; D[0][2]=0;
-		D[1][0]=0; D[1][1]=D_scalar; D[1][2]=0;
-		D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
+		D_scalar_conduct=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
+		if(dt) D_scalar_conduct=D_scalar_conduct*dt;
+
+		D[0][0]=D_scalar_conduct; D[0][1]=0; D[0][2]=0;
+		D[1][0]=0; D[1][1]=D_scalar_conduct; D[1][2]=0;
+		D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_conduct;
 
 		TripleMultiply(&B_conduct[0][0],3,numdof,1,
@@ -2990,10 +2991,10 @@
 		vz_input->GetParameterValue(&w, gauss);
 
-		D_scalar=gauss->weight*Jdet;
-		if(dt) D_scalar=D_scalar*dt;
-
-		D[0][0]=D_scalar*u;D[0][1]=0;         D[0][2]=0;
-		D[1][0]=0;         D[1][1]=D_scalar*v;D[1][2]=0;
-		D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar*w;
+		D_scalar_advec=gauss->weight*Jdet;
+		if(dt) D_scalar_advec=D_scalar_advec*dt;
+
+		D[0][0]=D_scalar_advec*u;D[0][1]=0;         D[0][2]=0;
+		D[1][0]=0;         D[1][1]=D_scalar_advec*v;D[1][2]=0;
+		D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar_advec*w;
 
 		TripleMultiply(&B_advec[0][0],3,numdof,1,
@@ -3006,9 +3007,9 @@
 		if(dt){
 			GetNodalFunctionsP1(&L[0], gauss);
-			D_scalar=gauss->weight*Jdet;
-			D_scalar=D_scalar;
+			D_scalar_trans=gauss->weight*Jdet;
+			D_scalar_trans=D_scalar_trans;
 
 			TripleMultiply(&L[0],numdof,1,0,
-						&D_scalar,1,1,0,
+						&D_scalar_trans,1,1,0,
 						&L[0],1,numdof,0,
 						&Ke_gaussian_transient[0][0],0);
@@ -3022,8 +3023,8 @@
 		if(artdiff==1){
 			/*Build K: */
-			D_scalar=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
-			if(dt) D_scalar=D_scalar*dt;
-			K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
-			K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
+			D_scalar_artdiff=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel);
+			if(dt) D_scalar_artdiff=D_scalar_artdiff*dt;
+			K[0][0]=D_scalar_artdiff*pow(u,2);       K[0][1]=D_scalar_artdiff*fabs(u)*fabs(v);
+			K[1][0]=D_scalar_artdiff*fabs(u)*fabs(v);K[1][1]=D_scalar_artdiff*pow(v,2);
 
 			GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); 
@@ -3039,5 +3040,5 @@
 
 			tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
-			//if(dt) scalar_artdiff=scalar_artdiff*dt;
+			if(dt) scalar_artdiff=scalar_artdiff*dt;
 
 			scalar_artdiff=tau_parameter*gauss->weight*Jdet;
@@ -3828,5 +3829,5 @@
 	double viscosity,temperature;
 	double tau_parameter,diameter;
-	double u,v,w,scalar_artdiff;
+	double u,v,w;
 	double scalar_def,scalar_transient;
 	double temperature_list[NUMVERTICES];
@@ -3888,8 +3889,6 @@
 
 			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]);
+
+			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
 		}
 	}