Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp
===================================================================
--- /issm/trunk/src/c/parallel/objectivefunctionC.cpp	(revision 1860)
+++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp	(revision 1861)
@@ -46,8 +46,8 @@
 	int     numberofnodes;
 		
-	/*thermalstatic: */
-	int     thermalstatic=0;
+	/*steadystate: */
+	int     steadystate=0;
 	int     isstokes=0;
-	DataSet* results_thermalstatic=NULL;
+	DataSet* results_steadystate=NULL;
 	int dofs01[2]={0,1};
 
@@ -72,5 +72,5 @@
 	femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
 	femmodel->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
-	femmodel->parameters->FindParam((void*)&thermalstatic,"thermalstatic");
+	femmodel->parameters->FindParam((void*)&steadystate,"steadystate");
 	femmodel->parameters->FindParam((void*)&isstokes,"isstokes");
 
@@ -88,5 +88,5 @@
 	inputs->Add(control_type,param_g_copy,1,numberofnodes);
 
-	if(!thermalstatic){
+	if(!steadystate){
 		//Run diagnostic with updated parameters.
 		diagnostic_core_nonlinear(&u_g,NULL,NULL,NULL,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type);
@@ -94,11 +94,11 @@
 	}
 	else{
-		//Run full thermalstatic solution with updated parameters.
-		results_thermalstatic=new DataSet(ResultsEnum());
-		thermalstatic_core(results_thermalstatic,model,inputs);
+		//Run full steadystate solution with updated parameters.
+		results_steadystate=new DataSet(ResultsEnum());
+		steadystate_core(results_steadystate,model,inputs);
 		
 		//get u_g 
-		results_thermalstatic->FindResult(&u_g,"u_g");
-		delete results_thermalstatic;
+		results_steadystate->FindResult(&u_g,"u_g");
+		delete results_steadystate;
 
 		//extract the correct number of dofs (3 or 4)
Index: /issm/trunk/src/c/parallel/parallel.h
===================================================================
--- /issm/trunk/src/c/parallel/parallel.h	(revision 1860)
+++ /issm/trunk/src/c/parallel/parallel.h	(revision 1861)
@@ -18,5 +18,5 @@
 void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
 
-void thermalstatic_core(DataSet* results,Model* model, ParameterInputs* inputs);
+void steadystate_core(DataSet* results,Model* model, ParameterInputs* inputs);
 
 void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, DataSet* loads, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
Index: /issm/trunk/src/c/parallel/steadystate.cpp
===================================================================
--- /issm/trunk/src/c/parallel/steadystate.cpp	(revision 1861)
+++ /issm/trunk/src/c/parallel/steadystate.cpp	(revision 1861)
@@ -0,0 +1,147 @@
+/*!\file:  steadystate.cpp
+ * \brief: steadystate solution
+ */ 
+
+#include "../issm.h"
+#include "./parallel.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "steadystate"
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+
+int main(int argc,char* *argv){
+	
+	/*I/O: */
+	FILE* fid=NULL;
+	char* inputfilename=NULL;
+	char* outputfilename=NULL;
+	char* lockname=NULL;
+	int   numberofnodes;
+	int   qmu_analysis=0;
+
+	/*Model: */
+	Model* model=NULL;
+
+	/*Results: */
+	DataSet* results=NULL;
+	Result* result=NULL;
+	
+	ParameterInputs* inputs=NULL;
+	int waitonlock=0;
+	
+	double* u_g_initial=NULL;
+	double* p_g_initial=NULL;
+	double  dt;
+	Param*  param=NULL;
+
+	MODULEBOOT();
+
+	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
+	throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!");
+	#endif
+
+	PetscInitialize(&argc,&argv,(char *)0,"");  
+
+	/*Size and rank: */
+	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
+	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
+
+	_printf_("recover , input file name and output file name:\n");
+	inputfilename=argv[2];
+	outputfilename=argv[3];
+	lockname=argv[4];
+
+	/*Open handle to data on disk: */
+	fid=pfopen(inputfilename,"rb");
+
+	/*Initialize model structure: */
+	model=new Model();
+
+	_printf_("read and create finite element model:\n");
+	_printf_("\n   reading diagnostic horiz model data:\n");
+	CreateFemModel(model->DiagnosticHorizontal(),fid,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+	_printf_("\n   reading diagnostic vert model data:\n");
+	CreateFemModel(model->DiagnosticVertical(),fid,DiagnosticAnalysisEnum(),VertAnalysisEnum());
+	_printf_("\n   reading diagnostic stokes model data:\n");
+	CreateFemModel(model->DiagnosticStokes(),fid,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
+	_printf_("\n   reading diagnostic hutter model data:\n");
+	CreateFemModel(model->DiagnosticHutter(),fid,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
+	_printf_("\n   reading surface and bed slope computation model data:\n");
+	CreateFemModel(model->Slope(),fid,SlopeComputeAnalysisEnum(),NoneAnalysisEnum());
+	_printf_("\n   read and create thermal finite element model:\n");
+	CreateFemModel(model->Thermal(),fid,ThermalAnalysisEnum(),SteadyAnalysisEnum());
+	_printf_("\n   read and create melting finite element model:\n");
+	CreateFemModel(model->Melting(),fid,MeltingAnalysisEnum(),SteadyAnalysisEnum());
+
+	_printf_("initialize inputs:\n");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&u_g_initial,"u_g");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	model->Thermal()->parameters->FindParam((void*)&dt,"dt");
+	model->Thermal()->parameters->FindParam((void*)&p_g_initial,"p_g");
+
+	inputs=new ParameterInputs;
+	inputs->Add("velocity",u_g_initial,3,numberofnodes);
+	inputs->Add("pressure",p_g_initial,1,numberofnodes);
+	inputs->Add("dt",dt);
+	
+	/*erase velocities: */
+	param=(Param*)model->DiagnosticHorizontal()->parameters->FindParamObject("u_g");
+	model->DiagnosticHorizontal()->parameters->DeleteObject((Object*)param);
+
+	_printf_("initialize results:\n");
+	results=new DataSet(ResultsEnum());
+
+	/*are we running the solution sequence, or a qmu wrapper around it? : */
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
+	if(!qmu_analysis){
+
+		/*run diagnostic analysis: */
+		_printf_("call computational core:\n");
+		steadystate_core(results,model,inputs);
+
+	}
+	else{
+		/*run qmu analysis: */
+		_printf_("calling qmu analysis on steadystate core:\n");
+
+		#ifdef _HAVE_DAKOTA_ 
+		Qmux(model,inputs,ThermalstaticAnalysisEnum(),NoneAnalysisEnum());
+	 	#else
+		throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!");
+		#endif
+	}
+
+	/*Add analysis_type to results: */
+	result=new Result(results->Size()+1,0,1,"analysis_type","steadystate");
+	results->AddObject(result);
+
+	_printf_("process results:\n");
+	ProcessResults(&results,model,ThermalstaticAnalysisEnum());
+
+	_printf_("write results to disk:\n");
+	OutputResults(results,outputfilename);
+
+	_printf_("write lock file:\n");
+	model->DiagnosticHorizontal()->parameters->FindParam((void*)&waitonlock,"waitonlock");
+	if (waitonlock){
+		WriteLockFile(lockname);
+	}
+	_printf_("closing MPI and Petsc\n");
+	PetscFinalize(); 
+	
+
+	/*end module: */
+	MODULEEND();
+	
+	/*Free ressources */
+	xfree((void**)&u_g_initial);
+	xfree((void**)&p_g_initial);
+
+	return 0; //unix success return;
+}
Index: /issm/trunk/src/c/parallel/steadystate_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/steadystate_core.cpp	(revision 1861)
+++ /issm/trunk/src/c/parallel/steadystate_core.cpp	(revision 1861)
@@ -0,0 +1,146 @@
+/*!\file: steadystate_core.cpp
+ * \brief: core of the steadystate solution 
+ */ 
+
+#undef __FUNCT__ 
+#define __FUNCT__ "steadystate_core"
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./parallel.h"
+#include "../issm.h"
+
+void steadystate_core(DataSet* results,Model* model, ParameterInputs* inputs){
+
+	extern int my_rank;
+
+	/*fem models: */
+	FemModel* fem_dh=NULL;
+	FemModel* fem_dv=NULL;
+	FemModel* fem_dhu=NULL;
+	FemModel* fem_ds=NULL;
+	FemModel* fem_sl=NULL;
+	FemModel* fem_t=NULL;
+	FemModel* fem_m=NULL;
+
+	/*output: */
+	Result* result=NULL;
+	DataSet* results_thermal=NULL;
+	DataSet* results_diagnostic=NULL;
+
+	/*solutions: */
+	Vec u_g=NULL;
+	Vec old_u_g=NULL;
+	Vec t_g=NULL;
+	Vec t_g_average=NULL;
+	Vec old_t_g=NULL;
+	Vec p_g=NULL;
+	Vec m_g=NULL;
+	Vec du_g=NULL;
+	Vec dt_g=NULL;
+	double ndu,nu;
+	double ndt,nt;
+	double eps_rel;
+
+	/*flags: */
+	int debug=0;
+	int isstokes=0;
+	int numberofnodes;
+	int ndof;
+	int converged;
+	int step;
+
+	/*recover fem models: */
+	fem_dh=model->DiagnosticHorizontal();
+	fem_dv=model->DiagnosticVertical();
+	fem_ds=model->DiagnosticStokes();
+	fem_dhu=model->DiagnosticHutter();
+	fem_sl=model->Slope();
+	fem_t=model->Thermal();
+	fem_m=model->Melting();
+
+	//first recover parameters common to all solutions
+	fem_dh->parameters->FindParam((void*)&debug,"debug");debug=1;
+	fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	fem_dh->parameters->FindParam((void*)&eps_rel,"eps_rel");
+	fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
+
+	//initialize: 
+	converged=0;
+	step=1;
+
+	if (isstokes)ndof=4;
+	else ndof=3;
+
+	for(;;){
+	
+		if(debug)_printf_("%s%i\n","computing temperature and velocity for step: ",step);
+
+		//first compute temperature at steady state.
+		if (step>1){
+			inputs->Add("velocity",u_g,ndof,numberofnodes);
+		}
+		results_thermal=new DataSet(ResultsEnum()); 
+		thermal_core(results_thermal,model,inputs);
+	
+		//get t_g and m_g;
+		VecFree(&t_g);results_thermal->FindResult(&t_g,"t_g");
+		VecFree(&m_g);results_thermal->FindResult(&m_g,"m_g");
+		delete results_thermal;
+
+		//Add temperature to inputs.
+		//compute depth averaged temperature and add to inputs
+		VecDuplicatePatch(&t_g_average,t_g); 
+		FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->loads, fem_t->materials,"temperature");
+		inputs->Add("temperature_average",t_g_average,1,numberofnodes);
+		inputs->Add("temperature",t_g,1,numberofnodes);
+		VecFree(&t_g_average); //not needed anymore
+
+		//now compute diagnostic velocity using the steady state temperature.
+		results_diagnostic=new DataSet(ResultsEnum());
+		diagnostic_core(results_diagnostic,model, inputs);
+
+		//get p_g and u_g
+		VecFree(&u_g);results_diagnostic->FindResult(&u_g,"u_g");
+		VecFree(&p_g);results_diagnostic->FindResult(&p_g,"p_g");
+		delete results_diagnostic;
+
+		//convergence? 
+		if(step>1){
+			VecFree(&du_g);VecDuplicatePatch(&du_g,old_u_g);VecAYPX(du_g,-1.0,u_g);
+			VecNorm(du_g,NORM_2,&ndu); VecNorm(old_u_g,NORM_2,&nu);
+
+			VecFree(&dt_g);VecDuplicatePatch(&dt_g,old_t_g); VecAYPX(dt_g,-1.0,t_g);
+			VecNorm(dt_g,NORM_2,&ndu); VecNorm(old_t_g,NORM_2,&nu);
+
+					
+			if (debug) _printf_("%-60s%g\n                                     %s%g\n                                     %s%g%s\n",
+					  "      relative convergence criterion: velocity -> norm(du)/norm(u)=   ",ndu/nu*100," temperature -> norm(dt)/norm(t)=",ndt/nt*100," eps_rel:                        ",eps_rel*100," %");
+		
+			if ((ndu/nu<=eps_rel)  && (ndt/nt<=eps_rel)) converged=1;
+			else converged=0;
+		}
+		else{
+			converged=0;
+		}
+
+		VecFree(&old_u_g);VecDuplicatePatch(&old_u_g,u_g);
+		VecFree(&old_t_g);VecDuplicatePatch(&old_t_g,t_g);
+
+		step++;
+		if (converged)break;
+	}
+
+
+	/*Plug results into output dataset: */
+	result=new Result(results->Size()+1,0,1,"u_g",u_g);
+	results->AddObject(result);
+	result=new Result(results->Size()+1,0,1,"p_g",p_g);
+	results->AddObject(result);
+	result=new Result(results->Size()+1,0,1,"t_g",t_g);
+	results->AddObject(result);
+	result=new Result(results->Size()+1,0,1,"m_g",m_g);
+	results->AddObject(result);
+}
Index: sm/trunk/src/c/parallel/thermalstatic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermalstatic.cpp	(revision 1860)
+++ 	(revision )
@@ -1,147 +1,0 @@
-/*!\file:  thermalstatic.cpp
- * \brief: thermalstatic solution
- */ 
-
-#include "../issm.h"
-#include "./parallel.h"
-
-#undef __FUNCT__ 
-#define __FUNCT__ "thermalstatic"
-
-#ifdef HAVE_CONFIG_H
-	#include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-
-int main(int argc,char* *argv){
-	
-	/*I/O: */
-	FILE* fid=NULL;
-	char* inputfilename=NULL;
-	char* outputfilename=NULL;
-	char* lockname=NULL;
-	int   numberofnodes;
-	int   qmu_analysis=0;
-
-	/*Model: */
-	Model* model=NULL;
-
-	/*Results: */
-	DataSet* results=NULL;
-	Result* result=NULL;
-	
-	ParameterInputs* inputs=NULL;
-	int waitonlock=0;
-	
-	double* u_g_initial=NULL;
-	double* p_g_initial=NULL;
-	double  dt;
-	Param*  param=NULL;
-
-	MODULEBOOT();
-
-	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
-	throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!");
-	#endif
-
-	PetscInitialize(&argc,&argv,(char *)0,"");  
-
-	/*Size and rank: */
-	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
-	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
-
-	_printf_("recover , input file name and output file name:\n");
-	inputfilename=argv[2];
-	outputfilename=argv[3];
-	lockname=argv[4];
-
-	/*Open handle to data on disk: */
-	fid=pfopen(inputfilename,"rb");
-
-	/*Initialize model structure: */
-	model=new Model();
-
-	_printf_("read and create finite element model:\n");
-	_printf_("\n   reading diagnostic horiz model data:\n");
-	CreateFemModel(model->DiagnosticHorizontal(),fid,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
-	_printf_("\n   reading diagnostic vert model data:\n");
-	CreateFemModel(model->DiagnosticVertical(),fid,DiagnosticAnalysisEnum(),VertAnalysisEnum());
-	_printf_("\n   reading diagnostic stokes model data:\n");
-	CreateFemModel(model->DiagnosticStokes(),fid,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
-	_printf_("\n   reading diagnostic hutter model data:\n");
-	CreateFemModel(model->DiagnosticHutter(),fid,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
-	_printf_("\n   reading surface and bed slope computation model data:\n");
-	CreateFemModel(model->Slope(),fid,SlopeComputeAnalysisEnum(),NoneAnalysisEnum());
-	_printf_("\n   read and create thermal finite element model:\n");
-	CreateFemModel(model->Thermal(),fid,ThermalAnalysisEnum(),SteadyAnalysisEnum());
-	_printf_("\n   read and create melting finite element model:\n");
-	CreateFemModel(model->Melting(),fid,MeltingAnalysisEnum(),SteadyAnalysisEnum());
-
-	_printf_("initialize inputs:\n");
-	model->DiagnosticHorizontal()->parameters->FindParam((void*)&u_g_initial,"u_g");
-	model->DiagnosticHorizontal()->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
-	model->Thermal()->parameters->FindParam((void*)&dt,"dt");
-	model->Thermal()->parameters->FindParam((void*)&p_g_initial,"p_g");
-
-	inputs=new ParameterInputs;
-	inputs->Add("velocity",u_g_initial,3,numberofnodes);
-	inputs->Add("pressure",p_g_initial,1,numberofnodes);
-	inputs->Add("dt",dt);
-	
-	/*erase velocities: */
-	param=(Param*)model->DiagnosticHorizontal()->parameters->FindParamObject("u_g");
-	model->DiagnosticHorizontal()->parameters->DeleteObject((Object*)param);
-
-	_printf_("initialize results:\n");
-	results=new DataSet(ResultsEnum());
-
-	/*are we running the solution sequence, or a qmu wrapper around it? : */
-	model->DiagnosticHorizontal()->parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
-	if(!qmu_analysis){
-
-		/*run diagnostic analysis: */
-		_printf_("call computational core:\n");
-		thermalstatic_core(results,model,inputs);
-
-	}
-	else{
-		/*run qmu analysis: */
-		_printf_("calling qmu analysis on thermalstatic core:\n");
-
-		#ifdef _HAVE_DAKOTA_ 
-		Qmux(model,inputs,ThermalstaticAnalysisEnum(),NoneAnalysisEnum());
-	 	#else
-		throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!");
-		#endif
-	}
-
-	/*Add analysis_type to results: */
-	result=new Result(results->Size()+1,0,1,"analysis_type","thermalstatic");
-	results->AddObject(result);
-
-	_printf_("process results:\n");
-	ProcessResults(&results,model,ThermalstaticAnalysisEnum());
-
-	_printf_("write results to disk:\n");
-	OutputResults(results,outputfilename);
-
-	_printf_("write lock file:\n");
-	model->DiagnosticHorizontal()->parameters->FindParam((void*)&waitonlock,"waitonlock");
-	if (waitonlock){
-		WriteLockFile(lockname);
-	}
-	_printf_("closing MPI and Petsc\n");
-	PetscFinalize(); 
-	
-
-	/*end module: */
-	MODULEEND();
-	
-	/*Free ressources */
-	xfree((void**)&u_g_initial);
-	xfree((void**)&p_g_initial);
-
-	return 0; //unix success return;
-}
Index: sm/trunk/src/c/parallel/thermalstatic_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermalstatic_core.cpp	(revision 1860)
+++ 	(revision )
@@ -1,146 +1,0 @@
-/*!\file: thermalstatic_core.cpp
- * \brief: core of the thermalstatic solution 
- */ 
-
-#undef __FUNCT__ 
-#define __FUNCT__ "thermalstatic_core"
-
-#include "../toolkits/toolkits.h"
-#include "../objects/objects.h"
-#include "../shared/shared.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-#include "./parallel.h"
-#include "../issm.h"
-
-void thermalstatic_core(DataSet* results,Model* model, ParameterInputs* inputs){
-
-	extern int my_rank;
-
-	/*fem models: */
-	FemModel* fem_dh=NULL;
-	FemModel* fem_dv=NULL;
-	FemModel* fem_dhu=NULL;
-	FemModel* fem_ds=NULL;
-	FemModel* fem_sl=NULL;
-	FemModel* fem_t=NULL;
-	FemModel* fem_m=NULL;
-
-	/*output: */
-	Result* result=NULL;
-	DataSet* results_thermal=NULL;
-	DataSet* results_diagnostic=NULL;
-
-	/*solutions: */
-	Vec u_g=NULL;
-	Vec old_u_g=NULL;
-	Vec t_g=NULL;
-	Vec t_g_average=NULL;
-	Vec old_t_g=NULL;
-	Vec p_g=NULL;
-	Vec m_g=NULL;
-	Vec du_g=NULL;
-	Vec dt_g=NULL;
-	double ndu,nu;
-	double ndt,nt;
-	double eps_rel;
-
-	/*flags: */
-	int debug=0;
-	int isstokes=0;
-	int numberofnodes;
-	int ndof;
-	int converged;
-	int step;
-
-	/*recover fem models: */
-	fem_dh=model->DiagnosticHorizontal();
-	fem_dv=model->DiagnosticVertical();
-	fem_ds=model->DiagnosticStokes();
-	fem_dhu=model->DiagnosticHutter();
-	fem_sl=model->Slope();
-	fem_t=model->Thermal();
-	fem_m=model->Melting();
-
-	//first recover parameters common to all solutions
-	fem_dh->parameters->FindParam((void*)&debug,"debug");debug=1;
-	fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
-	fem_dh->parameters->FindParam((void*)&eps_rel,"eps_rel");
-	fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
-
-	//initialize: 
-	converged=0;
-	step=1;
-
-	if (isstokes)ndof=4;
-	else ndof=3;
-
-	for(;;){
-	
-		if(debug)_printf_("%s%i\n","computing temperature and velocity for step: ",step);
-
-		//first compute temperature at steady state.
-		if (step>1){
-			inputs->Add("velocity",u_g,ndof,numberofnodes);
-		}
-		results_thermal=new DataSet(ResultsEnum()); 
-		thermal_core(results_thermal,model,inputs);
-	
-		//get t_g and m_g;
-		VecFree(&t_g);results_thermal->FindResult(&t_g,"t_g");
-		VecFree(&m_g);results_thermal->FindResult(&m_g,"m_g");
-		delete results_thermal;
-
-		//Add temperature to inputs.
-		//compute depth averaged temperature and add to inputs
-		VecDuplicatePatch(&t_g_average,t_g); 
-		FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->loads, fem_t->materials,"temperature");
-		inputs->Add("temperature_average",t_g_average,1,numberofnodes);
-		inputs->Add("temperature",t_g,1,numberofnodes);
-		VecFree(&t_g_average); //not needed anymore
-
-		//now compute diagnostic velocity using the steady state temperature.
-		results_diagnostic=new DataSet(ResultsEnum());
-		diagnostic_core(results_diagnostic,model, inputs);
-
-		//get p_g and u_g
-		VecFree(&u_g);results_diagnostic->FindResult(&u_g,"u_g");
-		VecFree(&p_g);results_diagnostic->FindResult(&p_g,"p_g");
-		delete results_diagnostic;
-
-		//convergence? 
-		if(step>1){
-			VecFree(&du_g);VecDuplicatePatch(&du_g,old_u_g);VecAYPX(du_g,-1.0,u_g);
-			VecNorm(du_g,NORM_2,&ndu); VecNorm(old_u_g,NORM_2,&nu);
-
-			VecFree(&dt_g);VecDuplicatePatch(&dt_g,old_t_g); VecAYPX(dt_g,-1.0,t_g);
-			VecNorm(dt_g,NORM_2,&ndu); VecNorm(old_t_g,NORM_2,&nu);
-
-					
-			if (debug) _printf_("%-60s%g\n                                     %s%g\n                                     %s%g%s\n",
-					  "      relative convergence criterion: velocity -> norm(du)/norm(u)=   ",ndu/nu*100," temperature -> norm(dt)/norm(t)=",ndt/nt*100," eps_rel:                        ",eps_rel*100," %");
-		
-			if ((ndu/nu<=eps_rel)  && (ndt/nt<=eps_rel)) converged=1;
-			else converged=0;
-		}
-		else{
-			converged=0;
-		}
-
-		VecFree(&old_u_g);VecDuplicatePatch(&old_u_g,u_g);
-		VecFree(&old_t_g);VecDuplicatePatch(&old_t_g,t_g);
-
-		step++;
-		if (converged)break;
-	}
-
-
-	/*Plug results into output dataset: */
-	result=new Result(results->Size()+1,0,1,"u_g",u_g);
-	results->AddObject(result);
-	result=new Result(results->Size()+1,0,1,"p_g",p_g);
-	results->AddObject(result);
-	result=new Result(results->Size()+1,0,1,"t_g",t_g);
-	results->AddObject(result);
-	result=new Result(results->Size()+1,0,1,"m_g",m_g);
-	results->AddObject(result);
-}