Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 17602)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_la_theta.cpp	(revision 17603)
@@ -11,11 +11,12 @@
 
 	/*intermediary: */
-	Matrix<IssmDouble>*  Kff = NULL;
-	Matrix<IssmDouble>*  Kfs = NULL;
-	Vector<IssmDouble>*  ug  = NULL;
-	Vector<IssmDouble>*  uf  = NULL;
-	Vector<IssmDouble>*  pf  = NULL;
-	Vector<IssmDouble>*  df  = NULL;
-	Vector<IssmDouble>*  ys  = NULL;
+	Matrix<IssmDouble>*  Kff    = NULL;
+	Matrix<IssmDouble>*  Kfs    = NULL;
+	Vector<IssmDouble>*  ug_old = NULL;
+	Vector<IssmDouble>*  ug     = NULL;
+	Vector<IssmDouble>*  uf     = NULL;
+	Vector<IssmDouble>*  pf     = NULL;
+	Vector<IssmDouble>*  df     = NULL;
+	Vector<IssmDouble>*  ys     = NULL;
 	int  configuration_type;
 
@@ -31,8 +32,13 @@
 
 	/*Convergence criterion*/
-	int count = 0;
-	femmodel->parameters->SetParam(reCast<IssmDouble>(10.),AugmentedLagrangianREnum);
+	int        count   = 0;
+	IssmDouble eps_rel = .001;
+	femmodel->parameters->SetParam(0.6,AugmentedLagrangianREnum);
+	GetSolutionFromInputsx(&ug,femmodel);
 
 	while(true){
+
+		/*save pointer to old velocity*/
+		delete ug_old;ug_old=ug;
 
 		/*Solve KU=F*/
@@ -51,12 +57,23 @@
 		analysis->InputUpdateFromSolutionFSXTH(femmodel->elements,femmodel->parameters);
 
+		/*Check for convergence*/
+		Vector<IssmDouble>* dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
+		IssmDouble ndu=dug->Norm(NORM_TWO);   delete dug;
+		IssmDouble nu =ug_old->Norm(NORM_TWO);
+		if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+		if((ndu/nu)<eps_rel){
+			if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %\n");
+			break;
+		}
+		else{ 
+			if(VerboseConvergence()) _printf0_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %\n");
+		}
+
 		count++;
-		if(count>5) break;
+		if(count>20) break;
 	}
 
-	/*Check for convergence*/
-	//_error_("STOP");
-
 	delete ug;  
+	delete ug_old;  
 	delete analysis;
 }
