Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp	(revision 13063)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp	(revision 13064)
@@ -8182,5 +8182,4 @@
 	 * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
 
-	_error_("Not supported yet");
 	int        i;
 	IssmDouble a,c,d,z,s,viscosity;
@@ -8200,5 +8199,5 @@
 	surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
 	PentaRef::GetInputValue(&z,&z_list[0],gauss);
-	tau_perp = matpar->GetRhoIce() * matpar->GetG() * (s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
+	tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
 
 	/* Get eps_b*/
@@ -8208,5 +8207,5 @@
 	eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
 	if(eps_b==0.){
-		*pviscosity=2.5*pow(10.,17.);
+		*pviscosity = 2.5e+17;
 		return;
 	}
@@ -8216,28 +8215,17 @@
 	A=matice->GetA();
 
-	/*Solve for tau_par (http://en.wikipedia.org/wiki/Cubic_function)*/
-	a=A;
-	c=A*tau_perp*tau_perp;
-	d=-eps_b;
-	tau_par = -1./(3.*a) * pow(0.5*(27.*a*a*d + sqrt(27.*a*a*d*27.*a*a*d -4.*pow(-3.*a*c,3.))),1./3.);
-
-	//printf("=======================================\n");
-	//printf("tau_par = %g (%g=0?)\n",tau_par,a*tau_par*tau_par*tau_par+c*tau_par+d);
-	//printf("a=%g c=%g d=%g\n",a,c,d);
-	//IssmDouble coeff[4];
-	//coeff[0]=d;
-	//coeff[1]=c;
-	//coeff[2]=0.;
-	//coeff[3]=a;
-	//int numroots;
-	//IssmDouble roots[3];
-	//cubic(coeff,roots,&numroots);
-	//tau_par=roots[0];
-	//for(i=0;i<numroots;i++) printf(" %g ",roots[i]);
+	/*Solve for tau_perp*/
+	int        numroots;
+	IssmDouble roots[3];
+	a = A;
+	c = A *tau_perp*tau_perp;
+	d = - eps_b;
+	cubic(a,0.,c,d,roots,&numroots);
+	tau_par=roots[0];
 	//printf(" (%g =0?)\n",a*roots[0]*roots[0]*roots[0]+c*roots[0]+d);
-	
 
 	/*Viscosity*/
-	viscosity = 1./(2.*A)*pow(tau_par*tau_par + tau_perp*tau_perp ,-2.);
+	viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
+	//printf("denom = %g (%g,%g)\n",tau_par*tau_par + tau_perp*tau_perp,tau_par,tau_perp);
 	_assert_(!isnan(viscosity));