Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17649)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17650)
@@ -7,5 +7,5 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 10
+//#define FSANALYTICAL 2
 
 /*Model processing*/
@@ -3209,9 +3209,10 @@
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
 
-	int         i,meshtype,dim;
+	int         i,meshtype,dim,fe_FS;
 	IssmDouble  x_coord,y_coord,z_coord;
 	IssmDouble  Jdet,forcex,forcey,forcez;
 	IssmDouble *xyz_list = NULL;
 
+	element->FindParam(&fe_FS,FlowequationFeFSEnum);
 	element->FindParam(&meshtype,MeshTypeEnum);
 	switch(meshtype){
@@ -3271,4 +3272,11 @@
 	xDelete<IssmDouble>(vbasis);
 	xDelete<IssmDouble>(xyz_list);
+	if(fe_FS==XTaylorHoodEnum){
+		ElementVector* pe2=CreatePVectorFSViscousXTH(element);
+		ElementVector* pe3 = new ElementVector(pe,pe2);
+		delete pe;
+		delete pe2;
+		return pe3;
+	}
 	return pe;
 }/*}}}*/
@@ -3552,5 +3560,5 @@
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(Dstar);
-xDelete<IssmDouble>(d);
+	xDelete<IssmDouble>(d);
 	xDelete<IssmDouble>(D);
 	xDelete<IssmDouble>(tau);
@@ -4382,5 +4390,9 @@
 							+2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old));
 			}
+			/*Initial guess cannot be 0 otherwise log(0)  - inf*/
+			if(dnorm==0.) dnorm=1.;
 			NewtonSolveDnorm(&dnorm,coef1,coef2,coef3,n,dnorm);
+			_assert_(dnorm>=0.);
+			_assert_(!xIsNan<IssmDouble>(dnorm));
 
 			/*Create Ke*/
@@ -4414,4 +4426,7 @@
 			Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
 			Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
+			for(int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_xx[i]));
+			for(int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_yy[i]));
+			for(int i=0;i<3;i++) _assert_(!xIsNan<IssmDouble>(d_xx[i]));
 			element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
 			element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
@@ -4420,10 +4435,10 @@
 		else{
 			_assert_(tnumnodes==4);
-			Matrix3x3Solve(&d_xx[0],Ke,pe_xx);
-			Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
-			Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
-			Matrix3x3Solve(&d_zz[0],Ke,pe_zz);
-			Matrix3x3Solve(&d_xz[0],Ke,pe_xz);
-			Matrix3x3Solve(&d_yz[0],Ke,pe_yz);
+			Matrix4x4Solve(&d_xx[0],Ke,pe_xx);
+			Matrix4x4Solve(&d_yy[0],Ke,pe_yy);
+			Matrix4x4Solve(&d_xy[0],Ke,pe_xy);
+			Matrix4x4Solve(&d_zz[0],Ke,pe_zz);
+			Matrix4x4Solve(&d_xz[0],Ke,pe_xz);
+			Matrix4x4Solve(&d_yz[0],Ke,pe_yz);
 			element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
 			element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
Index: /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp	(revision 17649)
+++ /issm/trunk-jpl/src/c/shared/Numerics/NewtonSolveDnorm.cpp	(revision 17650)
@@ -15,5 +15,5 @@
 	/*trivial solution*/
 	if(c3==0.){
-		*pdnorm =0.;
+		*pdnorm = 0.;
 		return 0;
 	}
@@ -26,4 +26,5 @@
 
 	/*Initial guess*/
+	_assert_(dnorm>0.); 
 	IssmDouble y1 = log10(dnorm);
 
@@ -41,6 +42,9 @@
 		}
 
-		if(counter>50) _error_("Could not converge");
+		if(counter>50) break;
 	}
+
+	/*Avoid extremely large values that indicate non convergence*/
+	if(y2>50.) y2 = 50;
 
 	/*Assign output pointer*/
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 17649)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 17650)
@@ -19,4 +19,5 @@
 	Vector<IssmDouble>*  df     = NULL;
 	Vector<IssmDouble>*  ys     = NULL;
+	IssmDouble           eps_rel;
 	int  configuration_type,max_nonlinear_iterations;
 
@@ -26,5 +27,7 @@
 	/*Recover parameters: */
 	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
 	femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
+	femmodel->parameters->SetParam(.6,AugmentedLagrangianREnum);
 
 	/*Update constraints and initialize d and tau if necessary*/
@@ -33,10 +36,9 @@
 
 	/*Convergence criterion*/
-	int        count   = 0;
-	IssmDouble eps_rel = .001;
-	femmodel->parameters->SetParam(.6,AugmentedLagrangianREnum);
+	int  count = 0;
 	GetSolutionFromInputsx(&ug,femmodel);
 
 	while(true){
+		count++;
 
 		/*save pointer to old velocity*/
@@ -47,20 +49,9 @@
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf, Kfs, ys); delete Kfs;
-
-	//	pf->Echo();
-	//	IssmDouble* temp=Kff->ToSerial();
-	//	int m,n;
-	//	Kff->GetSize(&m,&n);
-	//	printarray(temp,m,n);
-	//	xDelete<IssmDouble>(temp);
-
 		Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); 
-		//uf->Echo();
 		delete Kff; delete pf; delete df;
-
 		Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
 
 		/*Update solution*/
-		_printf_("ug norm = " << ug->Norm(NORM_TWO) << "\n");
 		InputUpdateFromSolutionx(femmodel,ug); 
 
@@ -81,7 +72,11 @@
 		}
 
-		count++;
-		if(count>=max_nonlinear_iterations) break;
+		if(count>=max_nonlinear_iterations){
+			_printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
+			break;
+		}
 	}
+
+	if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count-1 << "\n");
 
 	delete ug;