Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 20682)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 20683)
@@ -224,6 +224,6 @@
 
 	/*Intermediaries*/
-	IssmDouble vx,vy,vz;
-	IssmDouble dvx[3],dvy[3],dvz[3];
+	IssmDouble vx,vy,vz,vmag;
+	IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3];
 	IssmDouble epso,epsprime; 
 	int         dim;
@@ -241,4 +241,10 @@
 	int numvertices = this->GetNumberOfVertices();
 	IssmDouble* lambdas = xNew<IssmDouble>(numvertices);
+	IssmDouble* vorticityx = xNew<IssmDouble>(numvertices);
+	IssmDouble* vorticityy = xNew<IssmDouble>(numvertices);
+	IssmDouble* vorticityz = xNew<IssmDouble>(numvertices);
+	IssmDouble* omega_corrx = xNew<IssmDouble>(numvertices);
+	IssmDouble* omega_corry = xNew<IssmDouble>(numvertices);
+	IssmDouble* omega_corrz = xNew<IssmDouble>(numvertices);
 
 	/* Start looping on the number of vertices: */
@@ -264,11 +270,36 @@
 			dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
 		}
-		EarlStrainrateQuantities(&epso,&epsprime,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+		/*Calculate velocity magnitude and its derivative*/
+		vmag = sqrt(vx*vx+vy*vy+vz*vz);
+		if(vmag<1e-12){
+			vmag=1e-12;
+			dvmag[0]=0;
+			dvmag[1]=0;
+			dvmag[2]=0;
+		}
+		else{
+			dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
+			dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
+			dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
+		}
+		EarlStrainrateQuantities(&epso,&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
 		lambdas[iv]=EarlLambdaS(epso,epsprime);
+		vorticityx[iv]=dvz[1]-dvy[2];
+		vorticityy[iv]=dvx[2]-dvz[0];
+		vorticityz[iv]=dvy[0]-dvx[1];
+		omega_corrx[iv] = -vorticityx[iv] + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag;
+		omega_corry[iv] = -vorticityy[iv] + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag;
+		omega_corrz[iv] = -vorticityz[iv] + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag;
 	}
 
 	/*Add Stress tensor components into inputs*/
 	this->AddInput(LambdaSEnum,lambdas,P1Enum);
-
+	this->AddInput(Outputdefinition21Enum,vorticityx,P1Enum);
+	this->AddInput(Outputdefinition22Enum,vorticityy,P1Enum);
+	this->AddInput(Outputdefinition23Enum,vorticityz,P1Enum);
+	this->AddInput(Outputdefinition31Enum,omega_corrx,P1Enum);
+	this->AddInput(Outputdefinition32Enum,omega_corry,P1Enum);
+	this->AddInput(Outputdefinition33Enum,omega_corrz,P1Enum);
+	
 	/*Clean up and return*/
 	delete gauss;
Index: /issm/trunk-jpl/src/c/classes/Materials/Matearl.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matearl.cpp	(revision 20682)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matearl.cpp	(revision 20683)
@@ -409,6 +409,20 @@
 	IssmDouble E,lambdas;
 	IssmDouble epso,epsprime_norm;
-	
-	EarlStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
+	IssmDouble vmag,dvmag[3];
+
+	/*Calculate velocity magnitude and its derivative*/
+	vmag = sqrt(vx*vx+vy*vy+vz*vz);
+	if(vmag<1e-12){
+		dvmag[0]=0;
+		dvmag[1]=0;
+		dvmag[2]=0;
+	}
+	else{
+		dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
+		dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
+		dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
+	}
+
+	EarlStrainrateQuantities(&epso,&epsprime_norm,vx,vy,vz,vmag,dvx,dvy,dvz,&dvmag[0]);
 	lambdas=EarlLambdaS(epso,epsprime_norm);
 
Index: /issm/trunk-jpl/src/c/shared/Elements/EarlComponents.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Elements/EarlComponents.cpp	(revision 20682)
+++ /issm/trunk-jpl/src/c/shared/Elements/EarlComponents.cpp	(revision 20683)
@@ -2,36 +2,20 @@
 #include "../Numerics/types.h"
 #include "../Exceptions/exceptions.h"
-
-void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
+#include "./elements.h"
+void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag){/*{{{*/
 
 	/*Intermediaries*/
-	IssmDouble vorticity[3],vorticity_norm;
+	IssmDouble omega[3],omega_norm;
 	IssmDouble nrsp[3],nrsp_norm;
 	IssmDouble eps[3][3],epso;
 	IssmDouble epsprime[3],epsprime_norm;
 
-	/*Create vorticity vector*/
-	_assert_(dvx && dvy && dvz);
-	vorticity[0] =  dvz[1] - dvy[2];
-	vorticity[1] =  dvx[2] - dvz[0];
-	vorticity[2] =  dvy[0] - dvx[1];
-
-	/*Normalize*/
-	vorticity_norm = sqrt(vorticity[0]*vorticity[0] + vorticity[1]*vorticity[1] + vorticity[2]*vorticity[2]);
-	if(vorticity_norm==0){
-		vorticity[0] = 0.;
-		vorticity[1] = 0.;
-		vorticity[2] = 1.;
-	}
-	else{
-		vorticity[0] =vorticity[0]/vorticity_norm;
-		vorticity[1] =vorticity[1]/vorticity_norm;
-		vorticity[2] =vorticity[2]/vorticity_norm;
-	}
+	/*Get omega, correction for rigid body rotation*/
+	EarlOmega(&omega[0],vx,vy,vz,vmag,dvx,dvy,dvz,dvmag);
 
 	/*Non-rotating shear plane*/
-	nrsp[0] =  vy*vorticity[2] - vz*vorticity[1];
-	nrsp[1] =  vz*vorticity[0] - vx*vorticity[2];
-	nrsp[2] =  vx*vorticity[1] - vy*vorticity[0];
+	nrsp[0] =  vy*omega[2] - vz*omega[1];
+	nrsp[1] =  vz*omega[0] - vx*omega[2];
+	nrsp[2] =  vx*omega[1] - vy*omega[0];
 
 	/*Normalize*/
@@ -81,5 +65,5 @@
 		for(int j=0;j<3;j++){
 			for(int k=0;k<3;k++){
-				epsprime[j] += -nrsp[i]*eps[i][k]*vorticity[k]*vorticity[j];
+				epsprime[j] += -nrsp[i]*eps[i][k]*omega[k]*omega[j];
 			}
 		}
@@ -88,8 +72,48 @@
 	/*Get norm of epsprime*/
 	epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
-
+	
 	/*Assign output pointers*/
 	*pepso=epso;
 	*pepsprime_norm=epsprime_norm;
+}/*}}}*/
+void EarlOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble omega_norm;
+
+	/*Create vorticity vector*/
+	_assert_(dvx && dvy && dvz && dvmag);
+	if(vmag<1e-12)vmag=1e-12;
+
+	/*Create vorticity vector, corrected for rigid body rotation
+	 * \overline{\omega} =\omega - \Omega
+	 *                   =-\nabla\times{\bf v} - 2*\nabla U\times{\bf v}/U;
+	 * check the magnitude of the second term -- if it is small, then the two
+	 * vorticities (omega and first term in omega) are approx. equal
+	 *
+	 * */
+	omega[0] = -(dvz[1] - dvy[2]) + 2*(dvmag[1]*vz - dvmag[2]*vy)/vmag;
+	omega[1] = -(dvx[2] - dvz[0]) + 2*(dvmag[2]*vx - dvmag[0]*vz)/vmag;
+	omega[2] = -(dvy[0] - dvx[1]) + 2*(dvmag[0]*vy - dvmag[1]*vx)/vmag;
+
+	/*Take out vorticity component aligned with the velocity*/
+	IssmDouble wdotv = vx/vmag*omega[0] + vy/vmag*omega[1] + vz/vmag*omega[2];
+	omega[0] = omega[0] - wdotv*vx/vmag;
+	omega[1] = omega[1] - wdotv*vy/vmag;
+	omega[2] = omega[2] - wdotv*vz/vmag;
+
+	/*Normalize*/
+	omega_norm = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
+	if(omega_norm==0){
+		omega[0] = 0.;
+		omega[1] = 0.;
+		omega[2] = 1.;
+	}
+	else{
+		omega[0] =omega[0]/omega_norm;
+		omega[1] =omega[1]/omega_norm;
+		omega[2] =omega[2]/omega_norm;
+	}
+
 }/*}}}*/
 IssmDouble EarlLambdaS(IssmDouble epso, IssmDouble epsprime_norm){/*{{{*/
Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Elements/elements.h	(revision 20682)
+++ /issm/trunk-jpl/src/c/shared/Elements/elements.h	(revision 20683)
@@ -16,5 +16,6 @@
 // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n);
 IssmDouble EarlLambdaS(IssmDouble epso, IssmDouble epsprime_norm);
-void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz);
+void EarlOmega(IssmDouble* omega,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz, IssmDouble* dvmag);
+void EarlStrainrateQuantities(IssmDouble *pepso, IssmDouble *pepsprime_norm, IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble vmag,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz,IssmDouble* dvmag);
 IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
 				 IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, IssmDouble signorm,