Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16921)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16922)
@@ -2271,5 +2271,5 @@
 		this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
 
-		element->ViscosityFS(&viscosity,xyz_list,gauss,vx_input,vy_input,vz_input);
+		element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
 		for(i=0;i<epssize;i++)     D[i*bsize+i] = + 2.*viscosity*gauss->weight*Jdet;
 		for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 16921)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 16922)
@@ -109,5 +109,5 @@
 }/*}}}*/
 
-void Element::ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+void Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
 
 	/*Intermediaries*/
@@ -116,5 +116,5 @@
 	IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];            */
 
-	if(vz_input){
+	if(dim==3){
 		/*3D*/
 		this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
@@ -207,33 +207,21 @@
 	*pviscosity=viscosity;
 }/*}}}*/
+void Element::ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];    */
+
+	this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+	material->GetViscosity2d(&viscosity, &epsilon[0]);
+
+	/*Assign output pointer*/
+	*pviscosity=viscosity;
+}/*}}}*/
+void Element::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon);
+}/*}}}*/
 void Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
 	this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
-}/*}}}*/
-void Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-	/*Compute the 3d Blatter/HOStrain Rate (5 components):
-	 *
-	 * epsilon=[exx eyy exy exz eyz]
-	 *
-	 * with exz=1/2 du/dz
-	 *      eyz=1/2 dv/dz
-	 *
-	 * the contribution of vz is neglected
-	 */
-
-	int i;
-	IssmDouble epsilonvx[5];
-	IssmDouble epsilonvy[5];
-
-	/*Check that both inputs have been found*/
-	if (!vx_input || !vy_input){
-		_error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n");
-	}
-
-	/*Get strain rate assuming that epsilon has been allocated*/
-	vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
-	vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
-
-	/*Sum all contributions*/
-	for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
 }/*}}}*/
 void Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
@@ -263,18 +251,30 @@
 	for(int i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
 }/*}}}*/
-void Element::ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble viscosity;
-	IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];    */
-
-	this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-	material->GetViscosity2d(&viscosity, &epsilon[0]);
-
-	/*Assign output pointer*/
-	*pviscosity=viscosity;
-}/*}}}*/
-void Element::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon);
+void Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+	/*Compute the 3d Blatter/HOStrain Rate (5 components):
+	 *
+	 * epsilon=[exx eyy exy exz eyz]
+	 *
+	 * with exz=1/2 du/dz
+	 *      eyz=1/2 dv/dz
+	 *
+	 * the contribution of vz is neglected
+	 */
+
+	int i;
+	IssmDouble epsilonvx[5];
+	IssmDouble epsilonvy[5];
+
+	/*Check that both inputs have been found*/
+	if (!vx_input || !vy_input){
+		_error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n");
+	}
+
+	/*Get strain rate assuming that epsilon has been allocated*/
+	vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
+	vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
+
+	/*Sum all contributions*/
+	for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
 }/*}}}*/
 void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16921)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16922)
@@ -66,5 +66,5 @@
 		void       StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
 		IssmDouble TMeltingPoint(IssmDouble pressure);
-		void       ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
+		void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
 		void       ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
 		void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);