Index: /issm/trunk/src/c/objects/Matice.cpp
===================================================================
--- /issm/trunk/src/c/objects/Matice.cpp	(revision 88)
+++ /issm/trunk/src/c/objects/Matice.cpp	(revision 89)
@@ -240,2 +240,67 @@
 	*pviscosity2=viscosity2;
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "MatIce::GetViscosity3d"
+void  Matice::GetViscosity3d(double* pviscosity3d, double* epsilon){
+
+	/*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
+	 *
+	 *                                 2*B
+	 * viscosity3d= -------------------------------------------------------------------
+	 *     2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
+	 *
+	 *     where mu is the viscotiy, B the flow law parameter , (u,v) the velocity 
+	 *     vector, and n the flow law exponent.
+	 *
+	 * If epsilon is NULL, it means this is the first time Emg is being run, and we 
+	 * return g, initial viscosity.
+	 */
+	
+	/*output: */
+	double viscosity3d;
+
+	/*input strain rate: */
+	double exx,eyy,exy,exz,eyz;
+
+	/*Intermediary value A and exponent e: */
+	double A,e;
+
+	if (n==1){
+		/*Viscous behaviour! viscosity3d=B: */
+		viscosity3d=B;
+	}
+	else{
+		if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
+				(epsilon[3]==0) && (epsilon[4]==0)){
+			viscosity3d=pow(10,14);
+		}
+		else{
+
+			/*Retrive strain rate components: */
+			exx=*(epsilon+0);
+			eyy=*(epsilon+1);
+			exy=*(epsilon+2);
+			exz=*(epsilon+3);
+			eyz=*(epsilon+4);
+
+			/*Build viscosity: viscosity3d=2*B/(2*A^e) */
+			A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy;
+			if(A==0){
+				/*Maxiviscosity3dm viscosity for 0 shear areas: */
+				viscosity3d=4.5*pow(10,17);
+			}
+			else{
+				e=(n-1)/2/n;
+			
+				viscosity3d=2*B/(2*pow(A,e));
+			}
+		}
+	}
+	#ifdef _DEBUG_
+	_printf_("Viscosity %lf\n",viscosity3d);
+	#endif
+
+	/*Assign output pointers:*/
+	*pviscosity3d=viscosity3d;
+}