Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11261)
+++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp	(revision 11262)
@@ -3207,4 +3207,5 @@
 	int        i,j,ig,found=0;
 	double     Jdet,u,v,w,um,vm,wm;
+	double     h,hx,hy,hz,vx,vy,vz,vel;
 	double     gravity,rho_ice,rho_water;
 	double     epsvel=2.220446049250313e-16;
@@ -3214,9 +3215,6 @@
 	double     tau_parameter,diameter;
 	double     xyz_list[NUMVERTICES][3];
-	double     B[3][numdof];
-	double     Bprime[3][numdof];
 	double     B_conduct[3][numdof];
 	double     B_advec[3][numdof];
-	double     B_stab[2][numdof];
 	double     Bprime_advec[3][numdof];
 	double     L[numdof];
@@ -3225,5 +3223,5 @@
 	double     D_scalar_trans,D_scalar_stab;
 	double     D[3][3];
-	double     K[2][2]={0.0};
+	double     K[3][3]={0.0};
 	Tria*      tria=NULL;
 	GaussPenta *gauss=NULL;
@@ -3262,5 +3260,4 @@
 		/*Conduction: */  
 		/*Need to change that depending on enthalpy value -> cold or temperate ice: */  
-
 		GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
 
@@ -3281,5 +3278,4 @@
 
 		/*Advection: */
-
 		GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); 
 		GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 
@@ -3291,11 +3287,12 @@
 		vym_input->GetInputValue(&vm,gauss);
 		vzm_input->GetInputValue(&wm,gauss);
+		vx=u-um; vy=v-vm; vz=w-wm;
 
 		D_scalar_advec=gauss->weight*Jdet;
 		if(dt) D_scalar_advec=D_scalar_advec*dt;
 
-		D[0][0]=D_scalar_advec*(u-um);D[0][1]=0;                    D[0][2]=0;
-		D[1][0]=0;                    D[1][1]=D_scalar_advec*(v-vm);D[1][2]=0;
-		D[2][0]=0;                    D[2][1]=0;                    D[2][2]=D_scalar_advec*(w-wm);
+		D[0][0]=D_scalar_advec*vx;D[0][1]=0;                D[0][2]=0;
+		D[1][0]=0;                D[1][1]=D_scalar_advec*vy;D[1][2]=0;
+		D[2][0]=0;                D[2][1]=0;                D[2][2]=D_scalar_advec*vz;
 
 		TripleMultiply(&B_advec[0][0],3,numdof,1,
@@ -3305,5 +3302,4 @@
 
 		/*Transient: */
-
 		if(dt){
 			GetNodalFunctionsP1(&L[0], gauss);
@@ -3318,24 +3314,25 @@
 
 		/*Artifficial diffusivity*/
-
 		if(stabilization==1){
 			/*Build K: */
-			D_scalar_stab=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel);
+			GetElementSizes(&hx,&hy,&hz);
+			vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
+			h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
+			K[0][0]=h/(2*vel)*fabs(vx*vx);  K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz);
+			K[1][0]=h/(2*vel)*fabs(vy*vx);  K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz);
+			K[2][0]=h/(2*vel)*fabs(vz*vx);  K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz);
+			D_scalar_stab=gauss->weight*Jdet;
 			if(dt) D_scalar_stab=D_scalar_stab*dt;
-			K[0][0]=D_scalar_stab*pow(u,2);       K[0][1]=D_scalar_stab*fabs(u)*fabs(v);
-			K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2);
-			_error_("TO BE RECODED");
-
-			//GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); 
-
-			TripleMultiply(&B_stab[0][0],2,numdof,1,
-						&K[0][0],2,2,0,
-						&B_stab[0][0],2,numdof,0,
+			for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j];
+
+			GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 
+
+			TripleMultiply(&Bprime_advec[0][0],3,numdof,1,
+						&K[0][0],3,3,0,
+						&Bprime_advec[0][0],3,numdof,0,
 						&Ke->values[0],1);
 		}
 		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);
 
@@ -3456,9 +3453,6 @@
 	double     tau_parameter,diameter;
 	double     xyz_list[NUMVERTICES][3];
-	double     B[3][numdof];
-	double     Bprime[3][numdof];
 	double     B_conduct[3][numdof];
 	double     B_advec[3][numdof];
-	double     B_stab[2][numdof];
 	double     Bprime_advec[3][numdof];
 	double     L[numdof];
@@ -3537,5 +3531,4 @@
 
 		/*Transient: */
-
 		if(dt){
 			GetNodalFunctionsP1(&L[0], gauss);
@@ -3550,5 +3543,4 @@
 
 		/*Artifficial diffusivity*/
-
 		if(stabilization==1){
 			/*Build K: */