Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 15054)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 15055)
@@ -321,5 +321,4 @@
 					./modules/InputConvergencex/InputConvergencex.cpp\
 					./modules/InputConvergencex/InputConvergencex.h\
-					./analyses/convergence.cpp\
 					./analyses/ProcessArguments.cpp\
 					./analyses/ResetBoundaryConditions.cpp\
@@ -333,4 +332,5 @@
 					./solutionsequences/solutionsequence_nonlinear.cpp\
 					./solutionsequences/solutionsequence_newton.cpp\
+					./solutionsequences/convergence.cpp\
 					./classes/Options/Options.h\
 					./classes/Options/Options.cpp\
@@ -361,6 +361,5 @@
 #}}}
 #Steadystate sources  {{{
-steadystate_sources = ./analyses/steadystate_core.cpp\
-							 ./analyses/steadystateconvergence.cpp
+steadystate_sources = ./analyses/steadystate_core.cpp
 #}}}
 #Prognostic sources  {{{
@@ -436,6 +435,4 @@
 					  ./analyses/control_core.cpp\
 					  ./analyses/controltao_core.cpp\
-					  ./analyses/controlrestart.cpp\
-					  ./analyses/controlconvergence.cpp\
 					  ./analyses/objectivefunction.cpp\
 					  ./analyses/gradient_core.cpp\
Index: /issm/trunk-jpl/src/c/analyses/CMakeLists.txt
===================================================================
--- /issm/trunk-jpl/src/c/analyses/CMakeLists.txt	(revision 15054)
+++ /issm/trunk-jpl/src/c/analyses/CMakeLists.txt	(revision 15055)
@@ -23,6 +23,5 @@
 # }}}
 # STEADYSTATE_SOURCES {{{
-set(STEADYSTATE_SOURCES $ENV{ISSM_DIR}/src/c/solutions/steadystate_core.cpp
-                  $ENV{ISSM_DIR}/src/c/solutions/steadystateconvergence.cpp PARENT_SCOPE)
+set(STEADYSTATE_SOURCES $ENV{ISSM_DIR}/src/c/solutions/steadystate_core.cpp PARENT_SCOPE)
 # }}}
 # PROGNOSTIC_SOURCES {{{
Index: /issm/trunk-jpl/src/c/analyses/analyses.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/analyses.h	(revision 15054)
+++ /issm/trunk-jpl/src/c/analyses/analyses.h	(revision 15055)
@@ -39,9 +39,4 @@
 IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs);
 
-//convergence:
-void convergence(bool* pconverged, Matrix<IssmDouble>* K_ff,Vector<IssmDouble>* p_f,Vector<IssmDouble>* u_f,Vector<IssmDouble>* u_f_old,Parameters* parameters);
-bool controlconvergence(IssmDouble J,IssmDouble tol_cm);
-bool steadystateconvergence(FemModel* femmodel);
-
 //optimization
 int GradJSearch(IssmDouble* search_vector,FemModel* femmodel,int step);
@@ -50,5 +45,4 @@
 void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prootpath,int argc,char **argv);
 void WriteLockFile(char* filename);
-void controlrestart(FemModel* femmodel,IssmDouble* J);
 void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);
 COMM EnvironmentInit(int argc,char** argv);
Index: /issm/trunk-jpl/src/c/analyses/control_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/control_core.cpp	(revision 15054)
+++ /issm/trunk-jpl/src/c/analyses/control_core.cpp	(revision 15055)
@@ -10,32 +10,35 @@
 #include "../solutionsequences/solutionsequences.h"
 
+/*Local prototypes*/
+bool controlconvergence(IssmDouble J, IssmDouble tol_cm);
+
 void control_core(FemModel* femmodel){
 
-	int     i,n;
+	int     i;
 
 	/*parameters: */
-	int     num_controls,num_responses;
-	int     nsteps;
-	IssmDouble  tol_cm;
-	bool    cm_gradient;
-	int     dim;
-	int     solution_type;
-	bool    isstokes;
-	bool    dakota_analysis=false;
+	int        num_controls,num_responses;
+	int        nsteps;
+	IssmDouble tol_cm;
+	bool       cm_gradient;
+	int        dim;
+	int        solution_type;
+	bool       isstokes;
+	bool       dakota_analysis = false;
 
-	int*    control_type = NULL;
-	IssmDouble* responses=NULL;
-	int*    step_responses=NULL;
-	IssmDouble* maxiter=NULL;
-	IssmDouble* cm_jump=NULL;
+	int        *control_type   = NULL;
+	IssmDouble *responses      = NULL;
+	int        *step_responses = NULL;
+	IssmDouble *maxiter        = NULL;
+	IssmDouble *cm_jump        = NULL;
 
 	/*intermediary: */
-	IssmDouble  search_scalar=1;
-	OptArgs optargs;
-	OptPars optpars;
+	IssmDouble search_scalar = 1;
+	OptArgs    optargs;
+	OptPars    optpars;
 
 	/*Solution and Adjoint core pointer*/
-	void (*solutioncore)(FemModel*)=NULL;
-	void (*adjointcore)(FemModel*)=NULL;
+	void (*solutioncore)(FemModel*) = NULL;
+	void (*adjointcore)(FemModel*)  = NULL;
 
 	/*output: */
@@ -75,5 +78,5 @@
 
 	/*Start looping: */
-	for(n=0;n<nsteps;n++){
+	for(int n=0;n<nsteps;n++){
 
 		/*Display info*/
@@ -115,5 +118,5 @@
 		#ifdef _HAVE_ADOLC_
 		IssmPDouble* J_passive=xNew<IssmPDouble>(nsteps);
-		for(int i=0;i<nsteps;i++)J_passive[i]=reCast<IssmPDouble>(J[i]);
+		for(i=0;i<nsteps;i++) J_passive[i]=reCast<IssmPDouble>(J[i]);
 		femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,1,0));
 		xDelete<IssmPDouble>(J_passive);
@@ -132,2 +135,14 @@
 	xDelete<IssmDouble>(J);
 }
+bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
+
+	bool converged=false;
+
+	/*Has convergence been reached?*/
+	if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
+		converged=true;
+		if(VerboseConvergence()) _pprintString_("      Convergence criterion reached: J = " << J << " < " << tol_cm);
+	}
+
+	return converged;
+}
Index: sm/trunk-jpl/src/c/analyses/controlconvergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/controlconvergence.cpp	(revision 15054)
+++ 	(revision )
@@ -1,26 +1,0 @@
-/*!\file: controlconvergence.cpp
- * \brief: determine convergence of control_core solution
- */ 
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "./analyses.h"
-#include "../classes/classes.h"
-#include "../shared/shared.h"
-#include "../modules/modules.h"
-
-bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
-
-	bool converged=false;
-
-	/*Has convergence been reached?*/
-	if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
-		converged=true;
-		if(VerboseConvergence()) _pprintString_("      Convergence criterion reached: J = " << J << " < " << tol_cm);
-	}
-
-	return converged;
-}
Index: sm/trunk-jpl/src/c/analyses/controlrestart.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/controlrestart.cpp	(revision 15054)
+++ 	(revision )
@@ -1,43 +1,0 @@
-/*!\file: controlrestart.cpp
- * \brief: save as much as possible, to be able to restart the control_core solution
- */ 
-
-#include "./analyses.h"
-#include "../modules/modules.h"
-#include "../shared/shared.h"
-
-void controlrestart(FemModel* femmodel,IssmDouble* J){
-
-	int      num_controls;
-	int*     control_type = NULL;
-	int      nsteps;
-	bool     dakota_analysis=true;
-
-	/*retrieve output file name: */
-	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-	femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
-	femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
-	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
-
-	/*only save if we are not running qmu analysis. We certainly don't want to save control results each time we 
-	 * run on control core!: */
-	if(!dakota_analysis){
-		/*we essentially want J and the parameter: */
-		for(int i=0;i<num_controls;i++) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]);
-		#ifdef _HAVE_ADOLC_
-		IssmPDouble* J_passive=xNew<IssmPDouble>(nsteps);
-		for(int i=0;i<nsteps;i++)J_passive[i]=reCast<IssmPDouble>(J[i]);
-		femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,1,0));
-		xDelete<IssmPDouble>(J_passive);
-		#else
-		femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J,nsteps,1,1,0));
-		#endif
-		//femmodel->results->AddObject(new GenericExternalResult<char*>(femmodel->results->Size()+1,InversionControlParametersEnum,EnumToStringx(control_type),1,0));
-
-		/*write to disk: */
-		OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
-	}
-
-	/*Clean up and return*/
-	xDelete<int>(control_type);
-}
Index: sm/trunk-jpl/src/c/analyses/convergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/convergence.cpp	(revision 15054)
+++ 	(revision )
@@ -1,146 +1,0 @@
-/*!\file: convergence.cpp
- * \brief: figure out if convergence has been reached
- */ 
-
-#include "../classes/classes.h"
-#include "../modules/modules.h"
-#include "../shared/shared.h"
-
-void convergence(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,Parameters* parameters){
-
-	/*output*/
-	bool converged=false;
-
-	/*intermediary*/
-	Vector<IssmDouble>* KU=NULL;
-	Vector<IssmDouble>* KUF=NULL;
-	Vector<IssmDouble>* KUold=NULL;
-	Vector<IssmDouble>* KUoldF=NULL;
-	Vector<IssmDouble>* duf=NULL;
-	IssmDouble ndu,nduinf,nu;
-	IssmDouble nKUF;
-	IssmDouble nKUoldF;
-	IssmDouble nF;
-	IssmDouble solver_residue,res;
-
-	/*convergence options*/
-	IssmDouble eps_res;
-	IssmDouble eps_rel;
-	IssmDouble eps_abs;
-	IssmDouble yts;
-
-	if(VerboseModule()) _pprintLine_("   checking convergence");
-
-	/*If uf is NULL in input, f-set is nil, model is fully constrained, therefore converged from 
-	 * the get go: */
-	if(uf->IsEmpty()){
-		*pconverged=true;
-		return;
-	}
-
-	/*get convergence options*/
-	parameters->FindParam(&eps_res,DiagnosticRestolEnum);
-	parameters->FindParam(&eps_rel,DiagnosticReltolEnum);
-	parameters->FindParam(&eps_abs,DiagnosticAbstolEnum);
-	parameters->FindParam(&yts,ConstantsYtsEnum);
-
-	/*Display solver caracteristics*/
-	if (VerboseConvergence()){
-
-		//compute KUF = KU - F = K*U - F
-		KU=uf->Duplicate(); Kff->MatMult(uf,KU);
-		KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0);
-
-		//compute norm(KUF), norm(F) and residue
-		nKUF=KUF->Norm(NORM_TWO);
-		nF=pf->Norm(NORM_TWO);
-		solver_residue=nKUF/nF;
-		_pprintLine_("\n" << "   solver residue: norm(KU-F)/norm(F)=" << solver_residue);
-
-		//clean up
-		delete KU;
-		delete KUF;
-	}
-
-	/*Force equilibrium (Mandatory)*/
-
-	//compute K[n]U[n-1]F = K[n]U[n-1] - F
-	_assert_(uf); _assert_(Kff);
-	KUold=uf->Duplicate(); Kff->MatMult(old_uf,KUold);
-	KUoldF=KUold->Duplicate();KUold->Copy(KUoldF); KUoldF->AYPX(pf,-1.0);
-	nKUoldF=KUoldF->Norm(NORM_TWO);
-	nF=pf->Norm(NORM_TWO);
-	res=nKUoldF/nF;
-	if (xIsNan<IssmDouble>(res)){
-		_pprintLine_("norm nf = " << nF << "f and norm kuold = " << nKUoldF << "f");
-		_error_("mechanical equilibrium convergence criterion is NaN!");
-	}
-
-	//clean up
-	delete KUold;
-	delete KUoldF;
-
-	//print
-	if(res<eps_res){
-		if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   mechanical equilibrium convergence criterion"<<res*100<< " < "<<eps_res*100<<" %");
-		converged=true;
-	}
-	else{ 
-		if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   mechanical equilibrium convergence criterion"<<res*100<<" > "<<eps_res*100<<" %");
-		converged=false;
-	}
-
-	/*Relative criterion (optional)*/
-	if (!xIsNan<IssmDouble>(eps_rel) || (VerboseConvergence())){
-
-		//compute norm(du)/norm(u)
-		duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
-		ndu=duf->Norm(NORM_TWO); nu=old_uf->Norm(NORM_TWO);
-
-		if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
-
-		//clean up
-		delete duf;
-
-		//print
-		if (!xIsNan<IssmDouble>(eps_rel)){
-			if((ndu/nu)<eps_rel){
-				if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %");
-			}
-			else{ 
-				if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %");
-				converged=false;
-			}
-		}
-		else _pprintLine_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " %");
-
-	}
-
-	/*Absolute criterion (Optional) = max(du)*/
-	if (!xIsNan<IssmDouble>(eps_abs) || (VerboseConvergence())){
-
-		//compute max(du)
-		duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
-		ndu=duf->Norm(NORM_TWO); nduinf=duf->Norm(NORM_INF);
-		if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
-
-		//clean up
-		delete duf;
-
-		//print
-		if (!xIsNan<IssmDouble>(eps_abs)){
-			if ((nduinf*yts)<eps_abs){
-				if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " < " << eps_abs << " m/yr");
-			}
-			else{
-				if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " > " << eps_abs << " m/yr");
-				converged=false;
-			}
-		}
-		else  _pprintLine_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " m/yr");
-
-	}
-
-	/*assign output*/
-	*pconverged=converged;
-}
Index: /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp	(revision 15054)
+++ /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp	(revision 15055)
@@ -14,4 +14,7 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+
+/*Local prototypes*/
+bool steadystateconvergence(FemModel* femmodel);
 
 void steadystate_core(FemModel* femmodel){
@@ -91,2 +94,27 @@
 	xDelete<int>(requested_outputs);
 }
+bool steadystateconvergence(FemModel* femmodel){
+
+	/*output: */
+	bool converged=false;
+	bool velocity_converged=false;
+	bool temperature_converged=false;
+
+	/*intermediary: */
+	int velocityenums[8]={VxEnum,VxPicardEnum,VyEnum,VyPicardEnum,VzEnum,VzPicardEnum,PressureEnum,PressurePicardEnum}; //pairs of enums (new and old) on which to carry out the converence tests
+	int temperatureenums[2]={TemperatureEnum,TemperatureOldEnum};
+	int convergencecriterion[1]={RelativeEnum}; //criterions for convergence, RelativeEnum or AbsoluteEnum
+	IssmDouble convergencecriterionvalue[1]; //value of criterion to be respected
+
+	/*retrieve parameters: */
+	femmodel->parameters->FindParam(&convergencecriterionvalue[0],SteadystateReltolEnum);
+
+	/*figure out convergence at the input level, because we don't have the solution vectors!: */
+	velocity_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&velocityenums[0],8,&convergencecriterion[0],&convergencecriterionvalue[0],1);
+	temperature_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&temperatureenums[0],2,&convergencecriterion[0],&convergencecriterionvalue[0],1);
+
+	if(velocity_converged && temperature_converged) converged=true;
+
+	/*return: */
+	return converged;
+}
Index: sm/trunk-jpl/src/c/analyses/steadystateconvergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/steadystateconvergence.cpp	(revision 15054)
+++ 	(revision )
@@ -1,39 +1,0 @@
-/*!\file: steadystateconvergence.cpp
- * \brief: determine convergence of steady state solution
- */ 
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-#include "./analyses.h"
-#include "../classes/classes.h"
-#include "../shared/shared.h"
-#include "../modules/modules.h"
-
-bool steadystateconvergence(FemModel* femmodel){
-
-	/*output: */
-	bool converged=false;
-	bool velocity_converged=false;
-	bool temperature_converged=false;
-
-	/*intermediary: */
-	int velocityenums[8]={VxEnum,VxPicardEnum,VyEnum,VyPicardEnum,VzEnum,VzPicardEnum,PressureEnum,PressurePicardEnum}; //pairs of enums (new and old) on which to carry out the converence tests
-	int temperatureenums[2]={TemperatureEnum,TemperatureOldEnum};
-	int convergencecriterion[1]={RelativeEnum}; //criterions for convergence, RelativeEnum or AbsoluteEnum 
-	IssmDouble convergencecriterionvalue[1]; //value of criterion to be respected
-
-	/*retrieve parameters: */
-	femmodel->parameters->FindParam(&convergencecriterionvalue[0],SteadystateReltolEnum);
-
-	/*figure out convergence at the input level, because we don't have the solution vectors!: */
-	velocity_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&velocityenums[0],8,&convergencecriterion[0],&convergencecriterionvalue[0],1);
-	temperature_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&temperatureenums[0],2,&convergencecriterion[0],&convergencecriterionvalue[0],1);
-
-	if(velocity_converged && temperature_converged) converged=true;
-
-	/*return: */
-	return converged;
-}
Index: /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp	(revision 15055)
+++ /issm/trunk-jpl/src/c/solutionsequences/convergence.cpp	(revision 15055)
@@ -0,0 +1,146 @@
+/*!\file: convergence.cpp
+ * \brief: figure out if convergence has been reached
+ */ 
+
+#include "../classes/classes.h"
+#include "../modules/modules.h"
+#include "../shared/shared.h"
+
+void convergence(bool* pconverged, Matrix<IssmDouble>* Kff,Vector<IssmDouble>* pf,Vector<IssmDouble>* uf,Vector<IssmDouble>* old_uf,Parameters* parameters){
+
+	/*output*/
+	bool converged=false;
+
+	/*intermediary*/
+	Vector<IssmDouble>* KU=NULL;
+	Vector<IssmDouble>* KUF=NULL;
+	Vector<IssmDouble>* KUold=NULL;
+	Vector<IssmDouble>* KUoldF=NULL;
+	Vector<IssmDouble>* duf=NULL;
+	IssmDouble ndu,nduinf,nu;
+	IssmDouble nKUF;
+	IssmDouble nKUoldF;
+	IssmDouble nF;
+	IssmDouble solver_residue,res;
+
+	/*convergence options*/
+	IssmDouble eps_res;
+	IssmDouble eps_rel;
+	IssmDouble eps_abs;
+	IssmDouble yts;
+
+	if(VerboseModule()) _pprintLine_("   checking convergence");
+
+	/*If uf is NULL in input, f-set is nil, model is fully constrained, therefore converged from 
+	 * the get go: */
+	if(uf->IsEmpty()){
+		*pconverged=true;
+		return;
+	}
+
+	/*get convergence options*/
+	parameters->FindParam(&eps_res,DiagnosticRestolEnum);
+	parameters->FindParam(&eps_rel,DiagnosticReltolEnum);
+	parameters->FindParam(&eps_abs,DiagnosticAbstolEnum);
+	parameters->FindParam(&yts,ConstantsYtsEnum);
+
+	/*Display solver caracteristics*/
+	if (VerboseConvergence()){
+
+		//compute KUF = KU - F = K*U - F
+		KU=uf->Duplicate(); Kff->MatMult(uf,KU);
+		KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0);
+
+		//compute norm(KUF), norm(F) and residue
+		nKUF=KUF->Norm(NORM_TWO);
+		nF=pf->Norm(NORM_TWO);
+		solver_residue=nKUF/nF;
+		_pprintLine_("\n" << "   solver residue: norm(KU-F)/norm(F)=" << solver_residue);
+
+		//clean up
+		delete KU;
+		delete KUF;
+	}
+
+	/*Force equilibrium (Mandatory)*/
+
+	//compute K[n]U[n-1]F = K[n]U[n-1] - F
+	_assert_(uf); _assert_(Kff);
+	KUold=uf->Duplicate(); Kff->MatMult(old_uf,KUold);
+	KUoldF=KUold->Duplicate();KUold->Copy(KUoldF); KUoldF->AYPX(pf,-1.0);
+	nKUoldF=KUoldF->Norm(NORM_TWO);
+	nF=pf->Norm(NORM_TWO);
+	res=nKUoldF/nF;
+	if (xIsNan<IssmDouble>(res)){
+		_pprintLine_("norm nf = " << nF << "f and norm kuold = " << nKUoldF << "f");
+		_error_("mechanical equilibrium convergence criterion is NaN!");
+	}
+
+	//clean up
+	delete KUold;
+	delete KUoldF;
+
+	//print
+	if(res<eps_res){
+		if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   mechanical equilibrium convergence criterion"<<res*100<< " < "<<eps_res*100<<" %");
+		converged=true;
+	}
+	else{ 
+		if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<"   mechanical equilibrium convergence criterion"<<res*100<<" > "<<eps_res*100<<" %");
+		converged=false;
+	}
+
+	/*Relative criterion (optional)*/
+	if (!xIsNan<IssmDouble>(eps_rel) || (VerboseConvergence())){
+
+		//compute norm(du)/norm(u)
+		duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
+		ndu=duf->Norm(NORM_TWO); nu=old_uf->Norm(NORM_TWO);
+
+		if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+
+		//clean up
+		delete duf;
+
+		//print
+		if (!xIsNan<IssmDouble>(eps_rel)){
+			if((ndu/nu)<eps_rel){
+				if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " < " << eps_rel*100 << " %");
+			}
+			else{ 
+				if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " > " << eps_rel*100 << " %");
+				converged=false;
+			}
+		}
+		else _pprintLine_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " %");
+
+	}
+
+	/*Absolute criterion (Optional) = max(du)*/
+	if (!xIsNan<IssmDouble>(eps_abs) || (VerboseConvergence())){
+
+		//compute max(du)
+		duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0);
+		ndu=duf->Norm(NORM_TWO); nduinf=duf->Norm(NORM_INF);
+		if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+
+		//clean up
+		delete duf;
+
+		//print
+		if (!xIsNan<IssmDouble>(eps_abs)){
+			if ((nduinf*yts)<eps_abs){
+				if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " < " << eps_abs << " m/yr");
+			}
+			else{
+				if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " > " << eps_abs << " m/yr");
+				converged=false;
+			}
+		}
+		else  _pprintLine_(setw(50) << left << "   Convergence criterion: max(du)" << nduinf*yts << " m/yr");
+
+	}
+
+	/*assign output*/
+	*pconverged=converged;
+}
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp	(revision 15054)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp	(revision 15055)
@@ -3,9 +3,9 @@
  */ 
 
+#include "./solutionsequences.h"
 #include "../toolkits/toolkits.h"
 #include "../classes/classes.h"
 #include "../shared/shared.h"
 #include "../modules/modules.h"
-#include "../analyses/analyses.h"
 
 void solutionsequence_newton(FemModel* femmodel){
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp	(revision 15054)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp	(revision 15055)
@@ -3,9 +3,9 @@
  */ 
 
+#include "./solutionsequences.h"
 #include "../toolkits/toolkits.h"
 #include "../classes/classes.h"
 #include "../shared/shared.h"
 #include "../modules/modules.h"
-#include "../analyses/analyses.h"
 
 void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads){
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp	(revision 15054)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp	(revision 15055)
@@ -3,9 +3,9 @@
  */ 
 
+#include "./solutionsequences.h"
 #include "../toolkits/toolkits.h"
 #include "../classes/classes.h"
 #include "../shared/shared.h"
 #include "../modules/modules.h"
-#include "../analyses/analyses.h"
 
 void solutionsequence_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads){
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h	(revision 15054)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h	(revision 15055)
@@ -6,6 +6,9 @@
 #define _SOLUTION_SEQUENCES_H_
 
-struct OptArgs;
 class FemModel;
+class Parameters;
+template <class doubletype> class Matrix;
+template <class doubletype> class Vector;
+#include "../shared/Numerics/types.h"
 
 void solutionsequence_thermal_nonlinear(FemModel* femmodel);
@@ -17,3 +20,6 @@
 void solutionsequence_adjoint_linear(FemModel* femmodel);
 
+/*convergence*/
+void convergence(bool* pconverged, Matrix<IssmDouble>* K_ff,Vector<IssmDouble>* p_f,Vector<IssmDouble>* u_f,Vector<IssmDouble>* u_f_old,Parameters* parameters);
+
 #endif