Index: /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp	(revision 20649)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matbafl.cpp	(revision 20650)
@@ -331,9 +331,10 @@
 	}
 	else{
+		dvx[2] = 0.;
 		vy = 0.;
 		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
 	}
 	vz = 0.;
-	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
+	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
 
 	/*Get enhancement factors and ko*/
@@ -371,9 +372,11 @@
 	}
 	else{
+		dvx[1] = 0.;
+		dvx[2] = 0.;
 		vy = 0.;
 		dvy[0] = 0.; dvy[1] = 0.; dvy[2] = 0.;
 	}
 	vz = 0.;
-	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
+	dvz[0] = 0.; dvz[1] = 0.; dvz[2] = -dvx[0]-dvy[1];
 
 	/*Get enhancement factors and ko*/
@@ -382,45 +385,4 @@
 	Input* ko_input = element->inputs->GetInput(MaterialsRheologyKobarEnum); _assert_(ko_input);
 	ec_input->GetInputValue(&Ec,gauss);
-	es_input->GetInputValue(&Es,gauss);
-	ko_input->GetInputValue(&ko,gauss);
-
-	/*Compute viscosity*/
-	*pviscosity=GetViscosityGeneral(ko,Ec,Es,vx,vy,vz,&dvx[0],&dvy[0],&dvz[0]);
-}/*}}}*/
-void  Matbafl::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
-void  Matbafl::ResetHooks(){/*{{{*/
-
-	this->element=NULL;
-
-	/*Get Element type*/
-	this->helement->reset();
-
-}
-/*}}}*/
-IssmDouble Matbafl::GetViscosityGeneral(IssmDouble ko,IssmDouble Ec, IssmDouble Es,IssmDouble vx,IssmDouble vy,IssmDouble vz,IssmDouble* dvx,IssmDouble* dvy,IssmDouble* dvz){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble viscosity;
-	IssmDouble vorticity[3],vorticity_norm;
-	IssmDouble nrsp[3],nrsp_norm;
-	IssmDouble eps[3][3],epso;
-	IssmDouble epsprime[3],epsprime_norm;
-	IssmDouble E,lambdas;
-
-	/*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;
@@ -472,5 +434,5 @@
 		for(int j=0;j<3;j++){
 			for(int k=0;k<3;k++){
-				epsprime[j] += nrsp[i]*eps[i][k]*nrsp[k]*nrsp[j];
+				epsprime[j] += -nrsp[i]*eps[i][k]*nrsp[k]*nrsp[j];
 			}
 		}
@@ -480,5 +442,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]*vorticity[k]*vorticity[j];
 			}
 		}
@@ -490,5 +452,5 @@
 	/*Compute lambda_s*/
 	_assert_(epso>0.);
-	lambdas = sqrt(2*epsprime_norm/(3*epso));
+	lambdas = sqrt(2*epsprime_norm*epsprime_norm/(3*epso*epso));
 
 	/*Get total enhancement factor E(lambdas)*/