Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18871)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 18872)
@@ -42,11 +42,9 @@
 /*}}}*/
 
 /*Object constructors and destructor*/
-FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm){/*{{{*/
+FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm,bool trace){/*{{{*/
 
 	/*configuration: */
-	int* analyses=NULL;
-	int  numanalyses;
 	int  solution_type;
 	int  ierr;
 
@@ -76,11 +74,10 @@
 
 	/*Create femmodel from input files: */
 	profiler->Tag(StartInit);
-	this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,solution_type);
+	this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,solution_type,trace,NULL);
 	profiler->Tag(FinishInit);
 
 	/*Free resources */
-	xDelete<int>(analyses);
 	xDelete<char>(lockfilename);
 	xDelete<char>(binfilename);
 	xDelete<char>(outbinfilename);
@@ -89,6 +86,20 @@
 
 }
 /*}}}*/
+FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/
+
+	bool traceon=true;
+	
+	/*Store the communicator, but do not set it as a global variable, as this has already 
+	 * been done by the FemModel that called this copy constructor: */
+	this->comm=incomm;
+	this->SetStaticComm();
+
+	/*Create femmodel from input files, with trace activated: */
+	this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,solution_type,traceon,X);
+
+}
+/*}}}*/
 FemModel::~FemModel(){/*{{{*/
 
 	/*Intermediary*/
@@ -130,7 +141,7 @@
 
 }
 /*}}}*/
-void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, const int in_solution_type){/*{{{*/
+void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/
 
 	/*intermediary*/
 	int         i;
@@ -147,10 +158,9 @@
 	this->solution_type    = in_solution_type;
 	this->analysis_counter = nummodels-1;   //point to last analysis_type carried out.
 	this->results          = new Results(); //not initialized by CreateDataSets
-
 	/*Open input file on cpu 0 and create IoModel */
 	if(my_rank==0) IOMODEL = pfopen0(inputfilename ,"rb");
-	IoModel* iomodel = new IoModel(IOMODEL);
+	IoModel* iomodel = new IoModel(IOMODEL,trace,X);
 
 	/*Figure out what analyses are activated for this solution*/
 	SolutionAnalysesList(&this->analysis_type_list,&this->nummodels,iomodel,this->solution_type);
@@ -187,12 +197,15 @@
 	if(my_rank==0) pfclose(IOMODEL,inputfilename);
 	pfclose(toolkitsoptionsfid,toolkitsfilename);
 
-	/*Open output file once for all and add output file name and file descriptor to parameters*/
+	/*Open output file once for all and add output file descriptor to parameters*/
 	output_fid=pfopen(outputfilename,"wb");
+	this->parameters->SetParam(output_fid,OutputFilePointerEnum);
+	
+	/*Now save all of these file names into parameters, you never know when you might need them: */
+	this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename));
+	this->parameters->AddObject(new StringParam(RootPathEnum,rootpath));
+	this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename));
 	this->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
-	this->parameters->SetParam(output_fid,OutputFilePointerEnum);
-
-	/*Save lock file name for later: */
 	this->parameters->AddObject(new StringParam(LockFileNameEnum,lockfilename));
 
 	/*Clean up*/
@@ -740,26 +753,35 @@
 			}
 		}
 		else{
-			switch(output_enum){
+			/*last chance for the output definition, if the enum is one of Outputdefinition[1-10]Enum:*/
+			if(output_enum>=Outputdefinition1Enum && output_enum <=Outputdefinition10Enum){
+				double_result = OutputDefinitionsResponsex(this,output_enum);
+				if(save_results){
+					results->AddResult(new GenericExternalResult<IssmPDouble>(results->Size()+1,output_string,reCast<IssmPDouble>(double_result),step,time));
+					continue;
+				}
+			}
+			else{
+				switch(output_enum){
 
-				/*Scalar output*/
-				case DivergenceEnum:               this->Divergencex(&double_result);               break;
-				case MaxDivergenceEnum:            this->MaxDivergencex(&double_result);            break;
-				case IceVolumeEnum:                this->IceVolumex(&double_result);                break;
-				case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break;
-				case MinVelEnum:                   this->MinVelx(&double_result);                   break;
-				case MaxVelEnum:                   this->MaxVelx(&double_result);                   break;
-				case MinVxEnum:                    this->MinVxx(&double_result);                    break;
-				case MaxVxEnum:                    this->MaxVxx(&double_result);                    break;
-				case MaxAbsVxEnum:                 this->MaxAbsVxx(&double_result);                 break;
-				case MinVyEnum:                    this->MinVyx(&double_result);                    break;
-				case MaxVyEnum:                    this->MaxVyx(&double_result);                    break;
-				case MaxAbsVyEnum:                 this->MaxAbsVyx(&double_result);                 break;
-				case MinVzEnum:                    this->MinVzx(&double_result);                    break;
-				case MaxVzEnum:                    this->MaxVzx(&double_result);                    break;
-				case MaxAbsVzEnum:                 this->MaxAbsVzx(&double_result);                 break;
-				case MassFluxEnum:                 this->MassFluxx(&double_result);                 break;
-				case TotalSmbEnum:                 this->TotalSmbx(&double_result);                 break;
+					/*Scalar output*/
+					case DivergenceEnum:               this->Divergencex(&double_result);               break;
+					case MaxDivergenceEnum:            this->MaxDivergencex(&double_result);            break;
+					case IceVolumeEnum:                this->IceVolumex(&double_result);                break;
+					case IceVolumeAboveFloatationEnum: this->IceVolumeAboveFloatationx(&double_result); break;
+					case MinVelEnum:                   this->MinVelx(&double_result);                   break;
+					case MaxVelEnum:                   this->MaxVelx(&double_result);                   break;
+					case MinVxEnum:                    this->MinVxx(&double_result);                    break;
+					case MaxVxEnum:                    this->MaxVxx(&double_result);                    break;
+					case MaxAbsVxEnum:                 this->MaxAbsVxx(&double_result);                 break;
+					case MinVyEnum:                    this->MinVyx(&double_result);                    break;
+					case MaxVyEnum:                    this->MaxVyx(&double_result);                    break;
+					case MaxAbsVyEnum:                 this->MaxAbsVyx(&double_result);                 break;
+					case MinVzEnum:                    this->MinVzx(&double_result);                    break;
+					case MaxVzEnum:                    this->MaxVzx(&double_result);                    break;
+					case MaxAbsVzEnum:                 this->MaxAbsVzx(&double_result);                 break;
+					case MassFluxEnum:                 this->MassFluxx(&double_result);                 break;
+					case TotalSmbEnum:                 this->TotalSmbx(&double_result);                 break;
 
 			   /*Scalar control output*/
 				case SurfaceAbsVelMisfitEnum:       SurfaceAbsVelMisfitx(&double_result,elements,nodes,vertices,loads,materials,parameters);        break;
@@ -777,69 +799,71 @@
 				case BalancethicknessMisfitEnum:    BalancethicknessMisfitx(&double_result);                                                        break;
 				case SurfaceAbsMisfitEnum:          SurfaceAbsMisfitx(&double_result); break;
 
-			   /*Vector */
-				default:
+				   /*Vector */
+					default:
 
-					/*Vector layout*/
-					int interpolation,nodesperelement,size;
-					int rank_interpolation=-1,rank_nodesperelement=-1;
+						/*Vector layout*/
+						int interpolation,nodesperelement,size;
+						int rank_interpolation=-1,rank_nodesperelement=-1;
 
-					/*Get interpolation (and compute input if necessary)*/
-					for(int j=0;j<elements->Size();j++){
-						Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
-						element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum);
-					}
-
-					/*Broadcast for cpus that do not have any elements*/
-					ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
-					ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
-					ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-					ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-
-					if(results_on_nodes){
-
-						/*Allocate matrices*/
-						int         nbe       = this->elements->NumberOfElements();
-						IssmDouble* values    = xNewZeroInit<IssmDouble>(nbe*nodesperelement);
-						IssmDouble* allvalues = xNew<IssmDouble>(nbe*nodesperelement);
-
-						/*Fill-in matrix*/
+						/*Get interpolation (and compute input if necessary)*/
 						for(int j=0;j<elements->Size();j++){
 							Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
-							element->ResultToPatch(values,nodesperelement,output_enum);
+							element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum);
 						}
 
-						/*Gather from all cpus*/
-						ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
-						xDelete<IssmDouble>(values);
+						/*Broadcast for cpus that do not have any elements*/
+						ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
+						ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
+						ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+						ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm());
 
-						if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time));
-						xDelete<IssmDouble>(allvalues);
+						if(results_on_nodes){
 
-					}
-					else{
+							/*Allocate matrices*/
+							int         nbe       = this->elements->NumberOfElements();
+							IssmDouble* values    = xNewZeroInit<IssmDouble>(nbe*nodesperelement);
+							IssmDouble* allvalues = xNew<IssmDouble>(nbe*nodesperelement);
 
-						/*Allocate vector depending on interpolation*/
-						switch(interpolation){
-							case P0Enum: size = this->elements->NumberOfElements(); break;
-							case P1Enum: size = this->vertices->NumberOfVertices(); break;
-							default:     _error_("Interpolation "<<EnumToStringx(interpolation)<<" not supported yet");
+							/*Fill-in matrix*/
+							for(int j=0;j<elements->Size();j++){
+								Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
+								element->ResultToPatch(values,nodesperelement,output_enum);
+							}
 
+							/*Gather from all cpus*/
+							ISSM_MPI_Allreduce((void*)values,(void*)allvalues,nbe*nodesperelement,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+							xDelete<IssmDouble>(values);
+
+							if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nbe,nodesperelement,step,time));
+							xDelete<IssmDouble>(allvalues);
+
 						}
-						Vector<IssmDouble> *vector_result = new Vector<IssmDouble>(size);
+						else{
 
-						/*Fill in vector*/
-						for(int j=0;j<elements->Size();j++){
-							Element* element=(Element*)elements->GetObjectByOffset(j);
-							element->ResultToVector(vector_result,output_enum);
+							/*Allocate vector depending on interpolation*/
+							switch(interpolation){
+								case P0Enum: size = this->elements->NumberOfElements(); break;
+								case P1Enum: size = this->vertices->NumberOfVertices(); break;
+								default:     _error_("Interpolation "<<EnumToStringx(interpolation)<<" not supported yet");
+
+							}
+							Vector<IssmDouble> *vector_result = new Vector<IssmDouble>(size);
+
+							/*Fill in vector*/
+							for(int j=0;j<elements->Size();j++){
+								Element* element=(Element*)elements->GetObjectByOffset(j);
+								element->ResultToVector(vector_result,output_enum);
+							}
+							vector_result->Assemble();
+
+							if(save_results)results->AddResult(new GenericExternalResult<Vector<IssmDouble>*>(results->Size()+1,output_enum,vector_result,step,time));
 						}
-						vector_result->Assemble();
+						isvec = true;
+						break;
+				}
+			}
 
-						if(save_results)results->AddResult(new GenericExternalResult<Vector<IssmDouble>*>(results->Size()+1,output_enum,vector_result,step,time));
-					}
-					isvec = true;
-					break;
-			}
 		}
 
 		/*Add result to Results*/
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 18871)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 18872)
@@ -45,13 +45,14 @@
 		ISSM_MPI_Comm        comm;                  //communicator for this particular model
 
 		/*constructors, destructors: */
-		FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init);
+		FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init,bool trace=false);
+		FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X);
 		~FemModel();
 
 		/*Methods:*/
 		void Echo();
 		FemModel* copy();
-		void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type);
+		void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,bool trace,IssmPDouble* X=NULL);
 		void SolutionAnalysesList(int** panalyses,int* pnumanalyses,IoModel* iomodel,int solutiontype);
 		void CleanUp(void);
 		void Solve(void);
