Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 642)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 643)
@@ -11,4 +11,5 @@
 int ParamEnum(void){                    return           5; }
 int RgbEnum(void){                      return           6; }
+int ResultEnum(void){                   return           7; }
 
 /*analysis types: */
@@ -38,4 +39,5 @@
 int MaterialsEnum(void){                return          44; }
 int ParametersEnum(void){               return          45; }
+int ResultsEnum(void){                  return          46; }
 
 /*Elements: */
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 642)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 643)
@@ -25,4 +25,5 @@
 int RgbEnum(void);
 int InputEnum(void);
+int ResultEnum(void);
 
 /*formulations: */
@@ -60,4 +61,5 @@
 int MaterialsEnum(void);
 int ParametersEnum(void);
+int ResultsEnum(void);
 
 /*Functions on enums: */
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 642)
+++ /issm/trunk/src/c/Makefile.am	(revision 643)
@@ -297,4 +297,6 @@
 					./objects/Node.h\
 					./objects/Node.cpp\
+					./objects/Result.h\
+					./objects/Result.cpp\
 					./objects/DakotaPlugin.h\
 					./objects/DakotaPlugin.cpp\
@@ -534,16 +536,15 @@
 					./parallel/diagnostic_core_nonlinear.cpp\
 					./parallel/thermal_core.cpp\
+					./parallel/thermal_core_nonlinear.cpp\
 					./parallel/CreateFemModel.cpp\
-					./parallel/OutputDiagnostic.cpp\
 					./parallel/WriteLockFile.cpp\
 					./parallel/control.cpp\
-					./parallel/OutputControl.cpp\
-					./parallel/OutputThermal.cpp\
-					./parallel/OutputPrognostic.cpp\
 					./parallel/objectivefunctionC.cpp\
 					./parallel/GradJCompute.cpp\
 					./parallel/qmu.cpp\
 					./parallel/SpawnCore.cpp\
-					./parallel/DakotaResponses.cpp
+					./parallel/DakotaResponses.cpp\
+					./parallel/ProcessResults.cpp\
+					./parallel/OutputResults.cpp
 
 libpISSM_a_CXXFLAGS = -fPIC -D_PARALLEL_   -D_C_
@@ -552,5 +553,5 @@
 bin_PROGRAMS = 
 else 
-bin_PROGRAMS = diagnostic.exe 
+bin_PROGRAMS = diagnostic.exe  thermal.exe
 #control.exe thermal.exe prognostic.exe
 endif
Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 642)
+++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 643)
@@ -25,6 +25,8 @@
 	int      numberofdofspernode;
 	int      dim;
-	int*     epart=NULL;
-	int*     part=NULL;
+	int*  epart=NULL;
+	int*  part=NULL;
+	double*  dpart=NULL;
+	int      elements_width;
 
 
@@ -76,4 +78,5 @@
 	param->SetInteger(model->ismacayealpattyn);
 	parameters->AddObject(param);
+
 
 	count++;
@@ -220,5 +223,4 @@
 		char*  descriptor=NULL;
 		char*  tag=NULL;
-		int    elements_width; //size of elements
 
 		descriptors=(char**)xmalloc(model->numberofresponses*sizeof(char*));
@@ -246,12 +248,4 @@
 		xfree((void**)&descriptors);
 
-		/*Width of elements: */
-		if(strcmp(model->meshtype,"2d")==0){
-			elements_width=3; //tria elements
-		}
-		else{
-			elements_width=6; //penta elements
-		}
-
 		#ifdef _PARALLEL_
 		/*partition grids in model->qmu_npart parts: */
@@ -259,14 +253,19 @@
 		if(strcmp(model->meshtype,"2d")==0){
 			ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+			elements_width=3; //tria elements
 		}
 		else{
 			ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
+			elements_width=6; //penta elements
 		}
 
 		MeshPartitionx(&epart, &part,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,model->qmu_npart);
+
+		dpart=(double*)xmalloc(model->numberofnodes*sizeof(double));
+		for(i=0;i<model->numberofnodes;i++)dpart[i]=part[i];
 
 		count++;
 		param= new Param(count,"qmu_part",DOUBLEVEC);
-		param->SetDoubleVec((double*)part,model->numberofnodes);
+		param->SetDoubleVec(dpart,model->numberofnodes);
 		parameters->AddObject(param);
 
@@ -276,4 +275,5 @@
 		xfree((void**)&epart);
 		xfree((void**)&part);
+		xfree((void**)&dpart);
 
 		#endif
Index: /issm/trunk/src/c/objects/Result.cpp
===================================================================
--- /issm/trunk/src/c/objects/Result.cpp	(revision 643)
+++ /issm/trunk/src/c/objects/Result.cpp	(revision 643)
@@ -0,0 +1,170 @@
+/*!\file Result.cpp
+ * \brief: implementation of the Result object
+ */
+
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "stdio.h"
+#include "./Result.h"
+#include <string.h>
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./ParameterInputs.h"
+#include "../shared/shared.h"
+#include "../include/typedefs.h"
+
+Result::Result(){
+	return;
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Result::Result"
+Result::Result(int result_id,double result_time,int result_step,char* result_fieldname,Vec result_field){
+
+	id=result_id;
+	time=result_time;
+	step=result_step;
+	
+	if(!result_fieldname){
+		throw  ErrorException(__FUNCT__," NULL fieldname in constructor argument");
+	}
+	else{
+		fieldname=(char*)xmalloc((strlen(result_fieldname)+1)*sizeof(char));
+		strcpy(fieldname,result_fieldname);
+	}
+
+	field=result_field; //do not copy, as the results are large in memory size.
+	dfield=NULL;
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Result::Result"
+Result::Result(int result_id,double result_time,int result_step,char* result_fieldname,double* result_field,int result_size){
+
+	id=result_id;
+	time=result_time;
+	step=result_step;
+	
+	if(!result_fieldname){
+		throw  ErrorException(__FUNCT__," NULL fieldname in constructor argument");
+	}
+	else{
+		fieldname=(char*)xmalloc((strlen(result_fieldname)+1)*sizeof(char));
+		strcpy(fieldname,result_fieldname);
+	}
+
+	dfield=result_field; //do not copy, as the results are large in memory size.
+	size=result_size;
+	field=NULL;
+}
+
+
+Result::~Result(){
+	xfree((void**)&fieldname);
+	VecFree(&field);
+	xfree((void**)&dfield);
+}
+		
+void Result::Echo(void){
+
+	printf("Result:\n");
+	printf("   id: %i\n",id);
+	printf("   time: %g\n",time);
+	printf("   step: %i\n",step);
+	printf("   field name: %s\n",fieldname);
+	if(field){
+		printf("   field pointer %p\n",field);
+	}
+	else{
+		printf("   field pointer %p\n",dfield);
+		printf("   field size %i\n",size);
+	}
+}
+	
+#undef __FUNCT__ 
+#define __FUNCT__ "Result::Marshall"
+void  Result::Marshall(char** pmarshalled_dataset){
+
+	throw  ErrorException(__FUNCT__," not supported yet!");
+
+}
+	
+#undef __FUNCT__ 
+#define __FUNCT__ "Result::MarshallSize"
+int   Result::MarshallSize(){
+
+	throw  ErrorException(__FUNCT__," not supported yet!");
+}
+
+char* Result::GetName(void){
+	return "result";
+}
+		
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Result::Demarshall"
+void  Result::Demarshall(char** pmarshalled_dataset){
+
+	throw  ErrorException(__FUNCT__," not supported yet!");
+}
+
+int Result::Enum(void){
+
+	return ResultEnum();
+
+}
+int    Result::GetId(void){ return id; }
+
+int    Result::MyRank(void){ 
+	extern int my_rank;
+
+	return my_rank; 
+}
+
+double Result::GetTime(){
+	return time;
+}
+
+int    Result::GetStep(){
+	return step;
+}
+
+char*    Result::GetFieldName(){
+	return fieldname;
+}
+		
+void   Result::WriteData(FILE* fid){
+
+	int length;
+
+	/*First write field name :*/
+	length=(strlen(fieldname)+1)*sizeof(char);
+
+	fwrite(&length,sizeof(int),1,fid);
+	fwrite(fieldname,length,1,fid);
+
+	/*Now write time and step: */
+	fwrite(&time,sizeof(double),1,fid);
+	fwrite(&step,sizeof(int),1,fid);
+
+	/*Now write field: */
+	fwrite(&size,sizeof(int),1,fid);
+	fwrite(dfield,size*sizeof(double),1,fid);
+}
+
+
+void  Result::GetField(Vec* pfield){
+	*pfield=field;
+}
+		
+void  Result::GetField(double** pfield){
+	*pfield=dfield;
+}
+
+Object* Result::copy() {
+	return new Result(*this); 
+}
Index: /issm/trunk/src/c/objects/Result.h
===================================================================
--- /issm/trunk/src/c/objects/Result.h	(revision 643)
+++ /issm/trunk/src/c/objects/Result.h	(revision 643)
@@ -0,0 +1,50 @@
+/*!\file Result.h
+ * \brief: header file for result object
+ */
+
+#ifndef _RESULT_H_
+#define _RESULT_H_
+
+#include "stdio.h"
+#include "./Object.h"
+#include "../toolkits/toolkits.h"
+
+class Result: public Object{
+
+	private: 
+		int    id;
+		double time;
+		int    step;
+		char*  fieldname;
+		Vec    field;
+		double* dfield;
+		int     size;
+
+	public:
+
+		Result();
+		Result(int result_id,double result_time,int result_step,char* result_fieldname,Vec result_field);
+		Result(int result_id,double result_time,int result_step,char* result_fieldname,double* result_field,int result_size);
+		~Result();
+
+		void  Echo();
+		int   GetId(void); 
+		int   MyRank(void);
+		void  Marshall(char** pmarshalled_dataset);
+		int   MarshallSize();
+		char* GetName();
+		void  Demarshall(char** pmarshalled_dataset);
+		int   Enum();
+		Object* copy();
+
+		double GetTime();
+		void  GetField(Vec* pfield);
+		void  GetField(double** pfield);
+		int    GetStep();
+		void   WriteData(FILE* fid);
+		char*  GetFieldName();
+
+};
+
+#endif  /* _RESULT_H_ */
+
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 642)
+++ /issm/trunk/src/c/objects/objects.h	(revision 643)
@@ -19,4 +19,5 @@
 #include "./Beam.h"
 #include "./Spc.h"
+#include "./Result.h"
 #include "./Rgb.h"
 #include "./Icefront.h"
Index: /issm/trunk/src/c/parallel/DakotaResponses.cpp
===================================================================
--- /issm/trunk/src/c/parallel/DakotaResponses.cpp	(revision 642)
+++ /issm/trunk/src/c/parallel/DakotaResponses.cpp	(revision 643)
@@ -9,8 +9,10 @@
 #endif
 
+#include "../DataSet/DataSet.h"    
+#include "../shared/shared.h"
+
 #undef __FUNCT__ 
 #define __FUNCT__ "DakotaResponses"
-    
-void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Vec u_g,Vec p_g,int analysis_type,int sub_analysis_type){
+void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,DataSet* results,int analysis_type,int sub_analysis_type){
 
 	int i;
Index: sm/trunk/src/c/parallel/OutputControl.cpp
===================================================================
--- /issm/trunk/src/c/parallel/OutputControl.cpp	(revision 642)
+++ 	(revision )
@@ -1,68 +1,0 @@
-/*
-	OutputControl.c: output model results for control solution.
-*/
-#undef __FUNCT__ 
-#define __FUNCT__ "OutputControl"
-
-#include "../toolkits/toolkits.h"
-#include "../shared/shared.h"
-#include "../io/io.h"
-#include "../objects/objects.h"
-			
-void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets){
-
-	int i;
-	extern int my_rank;
-	
-	/* output: */
-	FILE* fid=NULL;
-
-	/* standard output: */
-	Vec     partition_shifted=NULL;
-	double* serial_partition=NULL;
-
-	double* serial_ug=NULL;
-	int     one=1;
-	int     u_g_size;
-	int     nods;
-
-	/*serialize outputs: */
-	VecDuplicate(partition,&partition_shifted);
-	VecCopy(partition,partition_shifted);
-	VecShift(partition_shifted,1.0); //matlab indexing starts at 1
-
-	VecToMPISerial(&serial_partition,partition_shifted);
-	VecGetSize(partition,&nods);
-
-	VecToMPISerial(&serial_ug,u_g);
-	VecGetSize(u_g,&u_g_size);
-	
-	/* Open output file to write raw binary data: */
-	if(my_rank==0){
-		fid=pfopen(filename,"wb");
-
-		/*Write solution type: */
-		WriteDataToDisk((void*)"control",NULL,NULL,"String",fid);
-
-		/*Write partition: */
-		WriteDataToDisk(serial_partition,&nods,&one,"Mat",fid);
-	
-		/*Write solution to disk: */
-		WriteDataToDisk(serial_ug,&u_g_size,&one,"Mat",fid);
-
-		/*Write parameter to disk: */
-		WriteDataToDisk(p_g,&u_g_size,&one,"Mat",fid);
-
-		/*Write J to disk: */
-		WriteDataToDisk(&nsteps,NULL,NULL,"Integer",fid);
-		WriteDataToDisk(J,&nsteps,&one,"Mat",fid);
-
-		/*Close file: */
-		pfclose(fid,filename);
-	}
-
-	VecFree(&partition_shifted);
-	xfree((void**)&serial_partition);
-	xfree((void**)&serial_ug);
-
-}	
Index: sm/trunk/src/c/parallel/OutputDiagnostic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/OutputDiagnostic.cpp	(revision 642)
+++ 	(revision )
@@ -1,81 +1,0 @@
-
-/*
-	OutputDiagnostic.c: output model results for diagnostic solution.
-*/
-#undef __FUNCT__ 
-#define __FUNCT__ "OutputDiagnostic"
-
-#include "../toolkits/toolkits.h"
-#include "../shared/shared.h"
-#include "../io/io.h"
-#include "../objects/objects.h"
-
-void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename){
-
-	int i;
-	extern int my_rank;
-	
-	/* output: */
-	FILE* fid=NULL;
-
-	/*intermediary: */
-	FemModel* fem_dh=NULL;
-
-	NodeSets* nodesets=NULL;
-	Vec       partition=NULL;
-	char*     analysis_type="diagnostic";
-	
-	/* standard output: */
-	Vec partition_shifted=NULL;
-	double* serial_partition=NULL;
-	double* serial_ug=NULL;
-	double* serial_pg=NULL;
-	int     u_g_size;
-	int     p_g_size;
-	int     one=1;
-	int     nods;
-
-	/*recover fem models: */
-	fem_dh=femmodels+0;
-	
-	/*Recover diagnostic horiz femmodel: */
-	partition=fem_dh->partition;
-	
-	/*serialize outputs: */
-	VecDuplicate(partition,&partition_shifted);
-	VecCopy(partition,partition_shifted);
-	VecShift(partition_shifted,1.0); //matlab indexing starts at 1
-
-	VecToMPISerial(&serial_partition,partition_shifted);
-	VecGetSize(partition,&nods);
-
-	VecToMPISerial(&serial_ug,u_g);
-	VecGetSize(u_g,&u_g_size);
-	
-	VecToMPISerial(&serial_pg,p_g);
-	VecGetSize(p_g,&p_g_size);
-
-	/* Open output file to write raw binary data: */
-	if(my_rank==0){
-		fid=pfopen(filename,"wb");
-
-		/*Write solution type: */
-		WriteDataToDisk(analysis_type,NULL,NULL,"String",fid);
-
-		/*Write partition: */
-		WriteDataToDisk(serial_partition,&nods,&one,"Mat",fid);
-		
-		/*Write solution to disk: */
-		WriteDataToDisk(serial_ug,&u_g_size,&one,"Mat",fid);
-		WriteDataToDisk(serial_pg,&p_g_size,&one,"Mat",fid);
-
-		/*Close file: */
-		pfclose(fid,filename);
-	}
-
-	/*Free ressources: */
-	VecFree(&partition_shifted);
-	xfree((void**)&serial_partition);
-	xfree((void**)&serial_ug);
-	xfree((void**)&serial_pg);
-}	
Index: sm/trunk/src/c/parallel/OutputPrognostic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/OutputPrognostic.cpp	(revision 642)
+++ 	(revision )
@@ -1,69 +1,0 @@
-
-/*
-	OutputPrognostic.c: output model results for prognostic solution.
-*/
-#undef __FUNCT__ 
-#define __FUNCT__ "OutputPrognostic"
-
-#include "../toolkits/toolkits.h"
-#include "../shared/shared.h"
-#include "../io/io.h"
-#include "../objects/objects.h"
-	
-void OutputPrognostic(Vec h_g,FemModel* femmodel,char* filename){
-
-	int i;
-	extern int my_rank;
-	
-	/* output: */
-	FILE* fid=NULL;
-
-	NodeSets* nodesets=NULL;
-	Vec       partition=NULL;
-	char*     analysis_type="prognostic";
-	
-	/* standard output: */
-	Vec partition_shifted=NULL;
-	double* serial_partition=NULL;
-	
-	double* serial_hg=NULL;
-	int     h_g_size;
-	int     one=1;
-	int     nods;
-
-	/*Recover prognostic horiz femmodel: */
-	partition=femmodel->partition;
-	
-	/*serialize outputs: */
-	VecDuplicate(partition,&partition_shifted);
-	VecCopy(partition,partition_shifted);
-	VecShift(partition_shifted,1.0); //matlab indexing starts at 1
-
-	VecToMPISerial(&serial_partition,partition_shifted);
-	VecGetSize(partition,&nods);
-
-	VecToMPISerial(&serial_hg,h_g);
-	VecGetSize(h_g,&h_g_size);
-	
-	/* Open output file to write raw binary data: */
-	if(my_rank==0){
-		fid=pfopen(filename,"wb");
-
-		/*Write solution type: */
-		WriteDataToDisk(analysis_type,NULL,NULL,"String",fid);
-
-		/*Write partition: */
-		WriteDataToDisk(serial_partition,&nods,&one,"Mat",fid);
-		
-		/*Write solution to disk: */
-		WriteDataToDisk(serial_hg,&h_g_size,&one,"Mat",fid);
-
-		/*Close file: */
-		pfclose(fid,filename);
-	}
-
-	/*Free ressources: */
-	VecFree(&partition_shifted);
-	xfree((void**)&serial_partition);
-	xfree((void**)&serial_hg);
-}	
Index: /issm/trunk/src/c/parallel/OutputResults.cpp
===================================================================
--- /issm/trunk/src/c/parallel/OutputResults.cpp	(revision 643)
+++ /issm/trunk/src/c/parallel/OutputResults.cpp	(revision 643)
@@ -0,0 +1,39 @@
+/*!\file:  OutputResults.cpp
+ * \brief: go through results dataset, and for each result, write it to disk.
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "stdio.h"
+#include "../DataSet/DataSet.h"
+#include "../io/io.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "OutputResults"
+void OutputResults(DataSet* results,char* filename){
+
+	int i;
+	Result* result=NULL;
+	FILE* fid=NULL;
+	extern int my_rank;
+
+	/* Open output file to write raw binary data: */
+	if(my_rank==0){
+		fid=pfopen(filename,"wb");
+
+		for(i=0;i<results->Size();i++){
+			result=(Result*)results->GetObjectByOffset(i);
+
+			/*write result to disk: */
+			result->WriteData(fid);
+		}
+		
+		/*Close file: */
+		pfclose(fid,filename);
+	}
+
+}
Index: sm/trunk/src/c/parallel/OutputThermal.cpp
===================================================================
--- /issm/trunk/src/c/parallel/OutputThermal.cpp	(revision 642)
+++ 	(revision )
@@ -1,112 +1,0 @@
-
-/*
-	OutputThermal.c: output model results for thermal solution.
-*/
-#undef __FUNCT__ 
-#define __FUNCT__ "OutputThermal"
-
-#include "../toolkits/toolkits.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-#include "../shared/shared.h"
-#include "../io/io.h"
-#include "../objects/objects.h"
-
-void OutputThermal(Vec* t_g,Vec* m_g, double* time,FemModel* femmodels,char* filename){
-
-	int i;
-	extern int my_rank;
-	
-	/* output: */
-	FILE* fid=NULL;
-
-	/*intermediary: */
-	FemModel* fem=NULL;
-	NodeSets* nodesets=NULL;
-	Vec       partition=NULL;
-	char*     analysis_type="thermal";
-	
-	/* standard output: */
-	Vec     partition_shifted=NULL;
-	double* serial_partition=NULL;
-
-	double** serial_tg=NULL;
-	double** serial_mg=NULL;
-
-	double ndt,dt;
-	int    nsteps;
-	int    nstepsplusone;
-	int    sub_analysis_type;
-	
-	int     one=1;
-	int     nods;
-
-	/*Recover thermal horiz femmodel: */
-	fem=femmodels+0;
-	
-	partition=fem->partition;
-	nodesets=fem->nodesets;
-
-	fem->parameters->FindParam((void*)&dt,"dt");
-	fem->parameters->FindParam((void*)&ndt,"ndt");
-	fem->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
-
-	/*serialize outputs: */
-	VecDuplicate(partition,&partition_shifted);
-	VecCopy(partition,partition_shifted);
-	VecShift(partition_shifted,1.0); //matlab indexing starts at 1
-	
-	VecToMPISerial(&serial_partition,partition_shifted);
-	VecGetSize(partition,&nods);
-
-
-	if(sub_analysis_type==SteadyAnalysisEnum()){
-		nsteps=0;
-	}
-	else{
-		nsteps=(int)(ndt/dt);
-	}
-	nstepsplusone=nsteps+1;
-
-	/*allocate t_g and m_g arrays: */
-	serial_tg=(double**)xmalloc((nsteps+1)*sizeof(double*));
-	serial_mg=(double**)xmalloc((nsteps+1)*sizeof(double*));
-
-	/*get tg_i*/
-	for(i=0;i<=nsteps;i++){
-		VecToMPISerial(&serial_tg[i],t_g[i]);
-		VecToMPISerial(&serial_mg[i],m_g[i]);
-	}
-	
-	/* Open output file to write raw binary data: */
-	if(my_rank==0){
-
-		fid=pfopen(filename,"wb");
-
-		/*Write solution type: */
-		WriteDataToDisk(analysis_type,NULL,NULL,"String",fid);
-
-		/*Write partition: */
-		WriteDataToDisk(serial_partition,&nods,&one,"Mat",fid);
-		
-		/*Write number of steps: */
-		WriteDataToDisk(&nsteps,NULL,NULL,"Integer",fid);
-
-		/*Write time steps: */
-		WriteDataToDisk(&time,&nstepsplusone,&one,"Mat",fid);
-
-		/*Write solutions to disk: */
-		for(i=0;i<=nsteps;i++){
-			WriteDataToDisk(serial_tg[i],&nods,&one,"Mat",fid);
-			WriteDataToDisk(serial_mg[i],&nods,&one,"Mat",fid);
-		}
-	
-		/*Close file: */
-		pfclose(fid,filename);
-	}
-
-	/*Free ressources: */
-	VecFree(&partition_shifted);
-	xfree((void**)&serial_partition);
-	xfree((void**)&serial_tg);
-	xfree((void**)&serial_mg);
-}	
Index: /issm/trunk/src/c/parallel/ProcessResults.cpp
===================================================================
--- /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 643)
+++ /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 643)
@@ -0,0 +1,180 @@
+/*!\file:  ProcessResults.cpp
+ * \brief: go through results dataset, and for each result, process it for easier retrieval 
+ * by the Matlab side. This usually means splitting the velocities from the g-size nodeset 
+ * to the grid set (ug->vx,vy,vz), same for pressure (p_g->pressure), etc ... It also implies 
+ * departitioning of the results.
+ * This phase is necessary prior to outputting the results on disk.
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#undef __FUNCT__ 
+#define __FUNCT__ "ProcessResults"
+
+#include "../DataSet/DataSet.h"
+#include "../objects/objects.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../shared/shared.h"
+
+void ProcessResults(DataSet** pnewresults,DataSet* results,FemModel* fems,int analysis_type){
+
+	int i,n;
+	Result* result=NULL;
+	Result* newresult=NULL;
+
+	/*output: */
+	DataSet* newresults=NULL;
+
+	/*fem models: */
+	FemModel* fem_dh=NULL;
+	FemModel* fem_dv=NULL;
+	FemModel* fem_dhu=NULL;
+	FemModel* fem_ds=NULL;
+	FemModel* fem_sl=NULL;
+
+	int ishutter;
+	int ismacayealpattyn;
+	int isstokes;
+
+	/*intermediary: */
+	Vec     u_g=NULL;
+	double* u_g_serial=NULL;
+	double* vx=NULL;
+	double* vy=NULL;
+	double* vz=NULL;
+	Vec     p_g=NULL;
+	double* p_g_serial=NULL;
+	double* pressure=NULL;
+	double* partition=NULL;
+	double  yts;
+
+	int numberofnodes;
+
+	/*Initialize new results: */
+	newresults=new DataSet(ResultsEnum());
+		
+
+	/*Recover femmodels first: */
+	if(analysis_type==DiagnosticAnalysisEnum()){
+
+		fem_dh=fems+0;
+		fem_dv=fems+1;
+		fem_ds=fems+2;
+		fem_dhu=fems+3;
+		fem_sl=fems+4;
+	
+		/*some flags needed: */
+		fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter");
+		fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
+		fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
+	}
+
+	
+
+	for(n=0;n<results->Size();n++){
+	
+		result=(Result*)results->GetObjectByOffset(n);
+		
+		if(strcmp(result->GetFieldName(),"u_g")==0){
+			/*Ok, are we dealing with velocities coming from MacAyeal, Pattyin, Hutter, on 2 dofs? or from 
+			 *Stokes on 4 dofs: */
+			result->GetField(&u_g);
+			VecToMPISerial(&u_g_serial,u_g);
+
+			if(!isstokes){
+				/*ok, 2 dofs, on number of nodes: */
+				if(ismacayealpattyn){
+					fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+					VecToMPISerial(&partition,fem_dh->partition);
+					fem_dh->parameters->FindParam((void*)&yts,"yts");
+				}
+				else{
+					fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+					VecToMPISerial(&partition,fem_dhu->partition);
+					fem_dhu->parameters->FindParam((void*)&yts,"yts");
+				}
+				vx=(double*)xmalloc(numberofnodes*sizeof(double));
+				vy=(double*)xmalloc(numberofnodes*sizeof(double));
+				vz=(double*)xmalloc(numberofnodes*sizeof(double)); 
+
+				for(i=0;i<numberofnodes;i++){
+					vx[i]=u_g_serial[2*(int)partition[i]+0]*yts;
+					vy[i]=u_g_serial[2*(int)partition[i]+1]*yts;
+					vz[i]=0;
+				}
+			}
+			else{
+				/* 4 dofs on number of nodes. discard pressure: */
+				fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+				VecToMPISerial(&partition,fem_ds->partition);
+				fem_ds->parameters->FindParam((void*)&yts,"yts");
+				vx=(double*)xmalloc(numberofnodes*sizeof(double));
+				vy=(double*)xmalloc(numberofnodes*sizeof(double));
+				vz=(double*)xmalloc(numberofnodes*sizeof(double));
+				for(i=0;i<numberofnodes;i++){
+					vx[i]=u_g_serial[4*(int)partition[i]+0]*yts;
+					vy[i]=u_g_serial[4*(int)partition[i]+1]*yts;
+					vz[i]=u_g_serial[4*(int)partition[i]+2]*yts;
+				}
+			}
+
+			/*Ok, add vx,vy and vz to newresults: */
+			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vx",vx,numberofnodes);
+			newresults->AddObject(newresult);
+
+			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vy",vy,numberofnodes);
+			newresults->AddObject(newresult);
+
+			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"vz",vz,numberofnodes);
+			newresults->AddObject(newresult);
+
+			/*do some cleanup: */
+			xfree((void**)&u_g_serial);
+			xfree((void**)&partition);
+		}
+		else if(strcmp(result->GetFieldName(),"p_g")==0){
+
+			/*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
+			result->GetField(&p_g);
+			VecToMPISerial(&p_g_serial,p_g);
+
+			if(!isstokes){
+				if(ismacayealpattyn){
+					fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+					VecToMPISerial(&partition,fem_dh->partition);
+				}
+				else{
+					fem_dhu->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+					VecToMPISerial(&partition,fem_dhu->partition);
+				}
+			}
+			else{
+				fem_ds->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+				VecToMPISerial(&partition,fem_ds->partition);
+			}
+
+			pressure=(double*)xmalloc(numberofnodes*sizeof(double));
+
+			for(i=0;i<numberofnodes;i++){
+				pressure[i]=p_g_serial[(int)partition[i]];
+			}
+			
+			/*Ok, add pressure,vy and vz to newresults: */
+			newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"pressure",pressure,numberofnodes);
+			newresults->AddObject(newresult);
+
+			/*do some cleanup: */
+			xfree((void**)&p_g_serial);
+			xfree((void**)&partition);
+		}
+	}
+
+	/*Assign output pointers:*/
+	*pnewresults=newresults;
+
+}
+
Index: /issm/trunk/src/c/parallel/SpawnCore.cpp
===================================================================
--- /issm/trunk/src/c/parallel/SpawnCore.cpp	(revision 642)
+++ /issm/trunk/src/c/parallel/SpawnCore.cpp	(revision 643)
@@ -34,4 +34,5 @@
 	double* qmu_part=NULL;
 	int     qmu_npart;
+	DataSet* results=NULL;
 
 	extern int my_rank;
@@ -91,5 +92,5 @@
 
 	/*Modify core inputs to reflect the dakota variables inputs: */
-	inputs->UpdateFromDakota(variables,variables_descriptors,numvariables,qmu_part,qmu_npart);
+	//inputs->UpdateFromDakota(variables,variables_descriptors,numvariables,qmu_part,qmu_npart);
 
 	/*Run the analysis core solution sequence, with the updated inputs: */
@@ -97,5 +98,5 @@
 
 		_printf_("Starting diagnostic core\n");
-		diagnostic_core(&u_g,&p_g,femmodels,inputs);
+		diagnostic_core(results,femmodels,inputs);
 
 	}
@@ -105,5 +106,5 @@
 
 	/*compute responses on cpu 0: dummy for now! */
-	DakotaResponses(responses,responses_descriptors,numresponses,u_g,p_g,analysis_type,sub_analysis_type);
+	DakotaResponses(responses,responses_descriptors,numresponses,results,analysis_type,sub_analysis_type);
 
 	/*Free ressources:*/
Index: /issm/trunk/src/c/parallel/control.cpp
===================================================================
--- /issm/trunk/src/c/parallel/control.cpp	(revision 642)
+++ /issm/trunk/src/c/parallel/control.cpp	(revision 643)
@@ -158,5 +158,5 @@
 			inputs->Add("fit",fit[n]);
 			diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
-			OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
+			//OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
 			_printf_("%s\n","      done.");
 		}
@@ -179,5 +179,5 @@
 	
 	_printf_("%s\n","      saving final results...");
-	OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
+	//OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
 	_printf_("%s\n","      done.");
 
Index: /issm/trunk/src/c/parallel/diagnostic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/diagnostic.cpp	(revision 642)
+++ /issm/trunk/src/c/parallel/diagnostic.cpp	(revision 643)
@@ -15,4 +15,5 @@
 #endif
 
+
 int main(int argc,char* *argv){
 	
@@ -28,7 +29,9 @@
 	/*Fem models : */
 	FemModel femmodels[5];
+
+	/*Results: */
+	DataSet* results=NULL;
+	DataSet* newresults=NULL;
 	
-	Vec u_g=NULL;
-	Vec p_g=NULL;
 	ParameterInputs* inputs=NULL;
 	int waitonlock=0;
@@ -70,4 +73,5 @@
 	CreateFemModel(&femmodels[4],fid,"slope_compute","");
 
+
 	_printf_("initialize inputs:\n");
 	femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g");
@@ -76,24 +80,33 @@
 	inputs=new ParameterInputs;
 	inputs->Add("velocity",u_g_initial,3,numberofnodes);
-
+	
 	/*erase velocities: */
 	param=(Param*)femmodels[0].parameters->FindParamObject("u_g");
 	femmodels[0].parameters->DeleteObject((Object*)param);
 
+	_printf_("initialize results:\n");
+	results=new DataSet(ResultsEnum());
+
 	/*are we running the solutoin sequence, or a qmu wrapper around it? : */
 	femmodels[0].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
 	if(!qmu_analysis){
+
 		/*run diagnostic analysis: */
 		_printf_("call computational core:\n");
-		diagnostic_core(&u_g,&p_g,femmodels,inputs);
+		diagnostic_core(results,femmodels,inputs);
 
-		_printf_("write results to disk:\n");
-		OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename);
 	}
 	else{
 		/*run qmu analysis: */
 		_printf_("calling qmu analysis on diagnostic core:\n");
+		
 		qmu(qmuname,&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
 	}
+
+	_printf_("process results:\n");
+	ProcessResults(&newresults,results,&femmodels[0],DiagnosticAnalysisEnum()); delete results;
+	
+	_printf_("write results to disk:\n");
+	OutputResults(newresults,outputfilename);
 
 	_printf_("write lock file:\n");
Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/diagnostic_core.cpp	(revision 642)
+++ /issm/trunk/src/c/parallel/diagnostic_core.cpp	(revision 643)
@@ -13,5 +13,5 @@
 #include "../issm.h"
 
-void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs){
+void diagnostic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs){
 
 	extern int my_rank;
@@ -23,4 +23,7 @@
 	FemModel* fem_ds=NULL;
 	FemModel* fem_sl=NULL;
+
+	/*output: */
+	Result* result=NULL;
 
 	/*solutions: */
@@ -171,6 +174,8 @@
 	}
 	
-	/*Assign output pointers: */
-	*pug=ug;
-	*ppg=pg;
+	/*Plug results into output dataset: */
+	result=new Result(results->Size()+1,0,1,"u_g",ug);
+	results->AddObject(result);
+	result=new Result(results->Size()+1,0,1,"p_g",pg);
+	results->AddObject(result);
 }
Index: /issm/trunk/src/c/parallel/parallel.h
===================================================================
--- /issm/trunk/src/c/parallel/parallel.h	(revision 642)
+++ /issm/trunk/src/c/parallel/parallel.h	(revision 643)
@@ -11,8 +11,11 @@
 Vec GradJCompute(ParameterInputs* inputs,FemModel* femmodel,double* u_g_obs);
 
-void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs);
+void diagnostic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
+
+void thermal_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
+void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
+
 void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
 void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int  analysis_type,int sub_analysis_type);
-void thermal_core(Vec* pt_g,double* pmelting_offset,FemModel* femmodel,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
 
 //int GradJOrth(WorkspaceParams* workspaceparams);
@@ -28,8 +31,5 @@
 
 //int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams);
-void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename);
-void OutputThermal(Vec* t_g,Vec* m_g, double* time,FemModel* femmodels,char* filename);
-void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets);
-void OutputPrognostic(Vec h_g,FemModel* femmodel,char* filename);
+void OutputResults(DataSet* results,char* filename);
 void WriteLockFile(char* filename);
 
@@ -38,5 +38,6 @@
 void qmu(const char* dakota_input_file,FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
 void SpawnCore(double* responses,double* variables,char** variable_descriptors,int numvariables, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
-void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Vec u_g,Vec p_g,int analysis_type,int sub_analysis_type);
+void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,DataSet* results,int analysis_type,int sub_analysis_type);
+void ProcessResults(DataSet** pnewresults,DataSet* results,FemModel* fems,int analysis_type);
 
 #endif
Index: /issm/trunk/src/c/parallel/prognostic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/prognostic.cpp	(revision 642)
+++ /issm/trunk/src/c/parallel/prognostic.cpp	(revision 643)
@@ -91,5 +91,5 @@
 
 	_printf_("write results to disk:\n");
-	OutputPrognostic(h_g,&fem,outputfilename);
+	//OutputPrognostic(h_g,&fem,outputfilename);
 
 	_printf_("write lock file:\n");
Index: /issm/trunk/src/c/parallel/thermal.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermal.cpp	(revision 642)
+++ /issm/trunk/src/c/parallel/thermal.cpp	(revision 643)
@@ -32,21 +32,15 @@
 	double* u_g=NULL;
 	double* p_g=NULL;
-	double* time=NULL;
-	double  dt;
-	double  ndt;
-	int     nsteps;
-	int     debug=0;
 
-	/*solution vectors: */
-	Vec* t_g=NULL;
-	Vec* m_g=NULL;
-
+	/*Results: */
+	DataSet* results=NULL;
+	DataSet* newresults=NULL;
+	
 	ParameterInputs* inputs=NULL;
 	Param*           param=NULL;
+	double  dt;
 
 	int    waitonlock=0;
-	int    sub_analysis_type;
-	double melting_offset;
-	
+		
 	MODULEBOOT();
 
@@ -77,9 +71,6 @@
 	femmodels[0].parameters->FindParam((void*)&p_g,"p_g");
 	femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
 	femmodels[0].parameters->FindParam((void*)&dt,"dt");
-	femmodels[0].parameters->FindParam((void*)&ndt,"ndt");
-	femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
-	femmodels[0].parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
-	femmodels[0].parameters->FindParam((void*)&debug,"debug");
 
 	inputs=new ParameterInputs;
@@ -99,58 +90,13 @@
 	femmodels[1].parameters->DeleteObject((Object*)param);
 
-	if(sub_analysis_type==SteadyAnalysisEnum()){
+	
+	_printf_("call computational core:\n");
+	thermal_core(results,femmodels,inputs);
 
-		time=(double*)xmalloc(sizeof(double));
-		time[0]=0;
-
-		/*allocate t_g and m_g arrays: */
-		t_g=(Vec*)xmalloc(sizeof(Vec));
-		m_g=(Vec*)xmalloc(sizeof(Vec));
-
-		if(debug)_printf_("computing temperatures:\n");
-		thermal_core(&t_g[0],&melting_offset,&femmodels[0],inputs,ThermalAnalysisEnum(),SteadyAnalysisEnum());
-		inputs->Add("temperature",t_g[0],1,numberofnodes);
-		inputs->Add("melting_offset",melting_offset);
-		
-		if(debug)_printf_("computing melting:\n");
-		diagnostic_core_linear(&m_g[0],&femmodels[1],inputs,MeltingAnalysisEnum(),SteadyAnalysisEnum());
-	}
-	else{
-		
-		nsteps=(int)(ndt/dt);
-		time=(double*)xmalloc((nsteps+1)*sizeof(double));
-		time[0]=0;
-
-		/*allocate t_g and m_g arrays: */
-		t_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec));
-		m_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec));
-
-		//initialize temperature and melting
-		femmodels[0].parameters->FindParam((void*)&t_g[0],"t_g");
-		femmodels[1].parameters->FindParam((void*)&m_g[0],"m_g");
-
-		//erase temperature and melting embedded in parameters
-		param=(Param*)femmodels[0].parameters->FindParamObject("t_g");
-		femmodels[0].parameters->DeleteObject((Object*)param);
-		param=(Param*)femmodels[1].parameters->FindParamObject("m_g");
-		femmodels[1].parameters->DeleteObject((Object*)param);
-
-		for(i=0;i<nsteps;i++){
-			if(debug)_printf_("time step: %i/%i\n",i+1,nsteps);
-			time[i+1]=(i+1)*dt;
-			
-			if(debug)_printf_("computing temperatures:\n");
-			inputs->Add("temperature",t_g[i],1,numberofnodes);
-			thermal_core(&t_g[i+1],&melting_offset,&femmodels[0],inputs,ThermalAnalysisEnum(),TransientAnalysisEnum());
-			
-			if(debug)_printf_("computing melting:\n");
-			inputs->Add("temperature",t_g[i+1],1,numberofnodes);
-			inputs->Add("melting_offset",melting_offset);
-			diagnostic_core_linear(&m_g[i+1],&femmodels[1],inputs,MeltingAnalysisEnum(),TransientAnalysisEnum());
-		}
-	}
+	_printf_("process results:\n");
+	ProcessResults(&newresults,results,&femmodels[0],ThermalAnalysisEnum()); delete results;
 
 	_printf_("write results to disk:\n");
-	OutputThermal(t_g,m_g,time,&femmodels[0],outputfilename);
+	OutputResults(newresults,outputfilename);
 
 	_printf_("write lock file:\n");
Index: /issm/trunk/src/c/parallel/thermal_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermal_core.cpp	(revision 642)
+++ /issm/trunk/src/c/parallel/thermal_core.cpp	(revision 643)
@@ -8,121 +8,119 @@
 #include "../toolkits/toolkits.h"
 #include "../objects/objects.h"
+#include "../shared/shared.h"
 #include "../EnumDefinitions/EnumDefinitions.h"
+#include "./parallel.h"
 #include "../issm.h"
 
-void thermal_core(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
+void thermal_core(DataSet* results,FemModel* fems, ParameterInputs* inputs){
 
-	/*solution : */
-	Vec tg=NULL; 
-	Vec tf=NULL; 
+	extern int my_rank;
+	int i;
+
+	/*fem models: */
+	FemModel* fem_t=NULL;
+	FemModel* fem_m=NULL;
+
+	/*output: */
+	Result* result=NULL;
+
+	/*solutions vectors: */
+	Vec* t_g=NULL;
+	Vec* m_g=NULL;
+	double* time=NULL;
+
+	/*flags: */
+	int debug=0;
+	int numberofdofspernode;
+	int numberofnodes;
+	int nsteps;
+	double  dt;
+	double  ndt;
+
+	int    sub_analysis_type;
 	double melting_offset;
+	
+	Param*           param=NULL;
 
-	/*intermediary: */
-	Mat Kgg_nopenalty=NULL;
-	Mat Kgg=NULL;
-	Mat Kff=NULL;
-	Mat Kfs=NULL;
-	Vec pg_nopenalty=NULL;
-	Vec pg=NULL;
-	Vec pf=NULL;
+	/*recover fem models: */
+	fem_t=fems+0;
+	fem_m=fems+1;
 
-	int converged;
-	int constraints_converged;
-	int num_unstable_constraints;
-	int count;
-	int numberofnodes;
-	int min_thermal_constraints;
+	//first recover parameters common to all solutions
+	fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	fem_t->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
+	fem_t->parameters->FindParam((void*)&debug,"debug");
+	fem_t->parameters->FindParam((void*)&dt,"dt");
+	fem_t->parameters->FindParam((void*)&ndt,"ndt");
 
-	/*parameters:*/
-	int kflag,pflag,connectivity,numberofdofspernode;
-	char* solver_string=NULL;
-	int debug=0;
-	int lowmem=0;
+	if(sub_analysis_type==SteadyAnalysisEnum()){
 
-	/*Recover parameters: */
-	kflag=1; pflag=1;
+		time=(double*)xmalloc(sizeof(double));
+		time[0]=0;
 
-	fem->parameters->FindParam((void*)&connectivity,"connectivity");
-	fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
-	fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
-	fem->parameters->FindParam((void*)&solver_string,"solverstring");
-	fem->parameters->FindParam((void*)&debug,"debug");
-	fem->parameters->FindParam((void*)&lowmem,"lowmem");
-	fem->parameters->FindParam((void*)&min_thermal_constraints,"min_thermal_constraints");
+		/*allocate t_g and m_g arrays: */
+		t_g=(Vec*)xmalloc(sizeof(Vec));
+		m_g=(Vec*)xmalloc(sizeof(Vec));
 
-	count=1;
-	converged=0;
-	for(;;){
+		if(debug)_printf_("computing temperatures:\n");
+		thermal_core_nonlinear(&t_g[0],&melting_offset,fem_t,inputs,ThermalAnalysisEnum(),SteadyAnalysisEnum());
+		inputs->Add("temperature",t_g[0],1,numberofnodes);
+		inputs->Add("melting_offset",melting_offset);
+		
+		if(debug)_printf_("computing melting:\n");
+		diagnostic_core_linear(&m_g[0],fem_m,inputs,MeltingAnalysisEnum(),SteadyAnalysisEnum());
+	}
+	else{
+		
+		nsteps=(int)(ndt/dt);
+		time=(double*)xmalloc((nsteps+1)*sizeof(double));
+		time[0]=0;
 
-		if(debug)_printf_("%s\n","starting direct shooting method");
+		/*allocate t_g and m_g arrays: */
+		t_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec));
+		m_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec));
 
-		/*Update parameters: */
-		UpdateFromInputsx(fem->elements,fem->nodes,fem->loads, fem->materials,inputs);
+		//initialize temperature and melting
+		fem_t->parameters->FindParam((void*)&t_g[0],"t_g");
+		fem_m->parameters->FindParam((void*)&m_g[0],"m_g");
 
-		//*Generate system matrices
-		if (!lowmem){
+		//erase temperature and melting embedded in parameters
+		param=(Param*)fem_t->parameters->FindParamObject("t_g");
+		fem_t->parameters->DeleteObject((Object*)param);
+		param=(Param*)fem_m->parameters->FindParamObject("m_g");
+		fem_m->parameters->DeleteObject((Object*)param);
 
-			/*Compute Kgg_nopenalty and pg_nopenalty once for all: */
-			if (count==1){
-				SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
-			}
+		for(i=0;i<nsteps;i++){
+			if(debug)_printf_("time step: %i/%i\n",i+1,nsteps);
+			time[i+1]=(i+1)*dt;
+			
+			if(debug)_printf_("computing temperatures:\n");
+			inputs->Add("temperature",t_g[i],1,numberofnodes);
+			thermal_core_nonlinear(&t_g[i+1],&melting_offset,fem_t,inputs,ThermalAnalysisEnum(),TransientAnalysisEnum());
+			
+			if(debug)_printf_("computing melting:\n");
+			inputs->Add("temperature",t_g[i+1],1,numberofnodes);
+			inputs->Add("melting_offset",melting_offset);
+			diagnostic_core_linear(&m_g[i+1],fem_m,inputs,MeltingAnalysisEnum(),TransientAnalysisEnum());
+		}
+	}
+	
+	/*Plug results into output dataset: */
+	if(sub_analysis_type==SteadyAnalysisEnum()){
+		result=new Result(results->Size()+1,0,1,"t_g",t_g[0]);
+		results->AddObject(result);
+		
+		result=new Result(results->Size()+1,0,1,"m_g",m_g[0]);
+		results->AddObject(result);
+	}
+	else{
+		for(i=0;i<nsteps+1;i++){
+			result=new Result(results->Size()+1,time[i],i+1,"t_g",t_g[i]);
+			results->AddObject(result);
 
-			/*Copy K_gg_nopenalty into Kgg, same for pg: */
-			Kgg=(Mat)xmalloc(sizeof(Mat));
-			MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg);
-			pg=(Vec)xmalloc(sizeof(Vec));
-			VecDuplicate(pg_nopenalty,&pg);VecCopy(pg_nopenalty,pg);
-
-			//apply penalties each time
-			PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); 
+			result=new Result(results->Size()+1,time[i],i+1,"m_g",m_g[i]);
+			results->AddObject(result);
 		}
-		else{
-			SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
-			//apply penalties
-			PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); 
-		}
-
-		/*!Reduce matrix from g to f size:*/
-		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);
-
-		/*Free ressources: */
-		MatFree(&Kgg);
-	
-		if (debug) _printf_("   reducing load from g to f set\n");
-		/*!Reduce load from g to f size: */
-		Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);
-
-		//no need for pg and Kfs anymore 
-		VecFree(&pg); 
-		MatFree(&Kfs);
-
-		/*Solve: */
-		if(debug)_printf_("%s\n","solving");
-		Solverx(&tf, Kff, pf, tf, solver_string);
-	
-		//no need for Kff and pf anymore
-		MatFree(&Kff);VecFree(&pf);
-
-		if (debug) _printf_("   merging solution from f to g set\n");
-		//Merge back to g set
-		Mergesolutionfromftogx(&tg, tf,fem->Gmn,fem->ys,fem->nodesets);
-
-		//Deal with penalty loads
-		if (debug) _printf_("   penalty constraints\n");
-		inputs->Add("temperature",tg,numberofdofspernode,numberofnodes);
-		
-		PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->loads,fem->materials,inputs,analysis_type,sub_analysis_type); 
-
-		if (!converged){
-			if(debug)_printf_("%s%i\n","   #unstable constraints = ",num_unstable_constraints);
-			if (num_unstable_constraints <= min_thermal_constraints)converged=1;
-		}
-		count++;
-		
-		if(converged==1)break;
 	}
 
-	/*Assign output pointers: */
-	*ptg=tg;
-	*pmelting_offset=melting_offset;
 }
Index: /issm/trunk/src/c/parallel/thermal_core_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermal_core_nonlinear.cpp	(revision 643)
+++ /issm/trunk/src/c/parallel/thermal_core_nonlinear.cpp	(revision 643)
@@ -0,0 +1,106 @@
+/*!\file: thermal_core_nonlinear.cpp
+ * \brief: core of the thermal solution 
+ */ 
+
+#undef __FUNCT__ 
+#define __FUNCT__ "thermal_core_nonlinear"
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../issm.h"
+
+void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
+
+	/*solution : */
+	Vec tg=NULL; 
+	Vec tf=NULL; 
+	double melting_offset;
+
+	/*intermediary: */
+	Mat Kgg=NULL;
+	Mat Kff=NULL;
+	Mat Kfs=NULL;
+	Vec pg=NULL;
+	Vec pf=NULL;
+
+	int converged;
+	int constraints_converged;
+	int num_unstable_constraints;
+	int count;
+	int numberofnodes;
+	int min_thermal_constraints;
+
+	/*parameters:*/
+	int kflag,pflag,connectivity,numberofdofspernode;
+	char* solver_string=NULL;
+	int debug=0;
+
+	/*Recover parameters: */
+	kflag=1; pflag=1;
+
+	fem->parameters->FindParam((void*)&connectivity,"connectivity");
+	fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
+	fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+	fem->parameters->FindParam((void*)&solver_string,"solverstring");
+	fem->parameters->FindParam((void*)&debug,"debug");
+	fem->parameters->FindParam((void*)&min_thermal_constraints,"min_thermal_constraints");
+
+	count=1;
+	converged=0;
+	for(;;){
+
+		if(debug)_printf_("%s\n","starting direct shooting method");
+
+		/*Update parameters: */
+		UpdateFromInputsx(fem->elements,fem->nodes,fem->loads, fem->materials,inputs);
+
+		//*Generate system matrices
+		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
+		//apply penalties
+		PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); 
+
+		/*!Reduce matrix from g to f size:*/
+		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);
+
+		/*Free ressources: */
+		MatFree(&Kgg);
+	
+		if (debug) _printf_("   reducing load from g to f set\n");
+		/*!Reduce load from g to f size: */
+		Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);
+
+		//no need for pg and Kfs anymore 
+		VecFree(&pg); 
+		MatFree(&Kfs);
+
+		/*Solve: */
+		if(debug)_printf_("%s\n","solving");
+		Solverx(&tf, Kff, pf, tf, solver_string);
+	
+		//no need for Kff and pf anymore
+		MatFree(&Kff);VecFree(&pf);
+
+		if (debug) _printf_("   merging solution from f to g set\n");
+		//Merge back to g set
+		Mergesolutionfromftogx(&tg, tf,fem->Gmn,fem->ys,fem->nodesets);
+
+		//Deal with penalty loads
+		if (debug) _printf_("   penalty constraints\n");
+		inputs->Add("temperature",tg,numberofdofspernode,numberofnodes);
+		
+		PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->loads,fem->materials,inputs,analysis_type,sub_analysis_type); 
+
+		if (!converged){
+			if(debug)_printf_("%s%i\n","   #unstable constraints = ",num_unstable_constraints);
+			if (num_unstable_constraints <= min_thermal_constraints)converged=1;
+		}
+		count++;
+		
+		if(converged==1)break;
+	}
+
+	/*Assign output pointers: */
+	*ptg=tg;
+	*pmelting_offset=melting_offset;
+}
Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 642)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 643)
@@ -276,6 +276,8 @@
 	md.qmu_analysis=0;
 	md.iresp=1;
-
 	md.npart=0;
+
+	%results
+	md.results=NaN;
 
 	%Ice solver string
Index: /issm/trunk/src/m/classes/public/ReadData.m
===================================================================
--- /issm/trunk/src/m/classes/public/ReadData.m	(revision 642)
+++ /issm/trunk/src/m/classes/public/ReadData.m	(revision 643)
@@ -1,3 +1,3 @@
-function field=ReadData(fid)
+function result=ReadData(fid)
 %READDATA - ...
 %
@@ -7,23 +7,20 @@
 
 %read field
-[type,count]=fread(fid,1,'double');
+[length,count]=fread(fid,1,'int');
 
 if count==0,
-	field=NaN;
+	result=struct([]);
 else
-	if type==0,
-		%string
-		stringsize=fread(fid,1,'double');
-		field=char(fread(fid,stringsize,'char'));
-		field=field(1:end-1)';
-	elseif type==1,
-		%matrix
-		M=fread(fid,1,'double');
-		N=fread(fid,1,'double');
-		field=fread(fid,[M,N],'double');
-	elseif ((type==2) || (type==3)),
-		field=fread(fid,1,'double');
-	else 
-		error('ReadData error message: data type not supported yet!');
-	end
+	fieldname=fread(fid,length,'char');
+	fieldname=fieldname(1:end-1)';
+	fieldname=char(fieldname);
+	time=fread(fid,1,'double');
+	step=fread(fid,1,'int');
+	ssize=fread(fid,1,'int');
+	field=fread(fid,ssize,'double');
+
+	result.fieldname=fieldname;
+	result.time=time;
+	result.step=step;
+	result.field=field;
 end
Index: /issm/trunk/src/m/classes/public/loadresultsfromdisk.m
===================================================================
--- /issm/trunk/src/m/classes/public/loadresultsfromdisk.m	(revision 642)
+++ /issm/trunk/src/m/classes/public/loadresultsfromdisk.m	(revision 643)
@@ -13,113 +13,5 @@
 
 
-results=parseresultsfromdisk(filename);
-
-%First get solution type
-analysis_type=results{1};
-
-if ~strcmpi(analysis_type,md.analysis_type),
-	error(['loadresultsfromdisk error message: trying  to load results from disk for analysis_type: ',md.analysis_type,' when results are from analysis_type: ',analysis_type]);
-end
-
-%Get part
-part=results{2};
-
-%now to specialized reading
-if strcmpi(analysis_type,'diagnostic'),
-
-	%Get u_g
-	u_g=results{3};
-	p_g=results{4};
-
-	if strcmpi(md.type,'2d'),
-		gsize=md.numberofgrids*2;
-		%Used to recover velocities
-		indx=1:2:gsize;
-		indy=2:2:gsize;
-		indx=indx(part);
-		indy=indy(part);
-
-		%Recover velocity
-		md.vx=u_g(indx)*md.yts;
-		md.vy=u_g(indy)*md.yts;
-		md.vel=sqrt(md.vx.^2+md.vy.^2);
-		md.pressure=p_g(part);
-	else
-		%Used to recover velocities
-		gsize=length(u_g);
-		offset=gsize/md.numberofgrids;
-		indx=1:offset:gsize;
-		indy=2:offset:gsize;
-		indz=3:offset:gsize;
-		indx=indx(part);
-		indy=indy(part);
-		indz=indz(part);
-
-		%Recover velocity
-		md.vx=u_g(indx)*md.yts;
-		md.vy=u_g(indy)*md.yts;
-		md.vz=u_g(indz)*md.yts;
-		md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
-		md.pressure=p_g(part);
-	end
-
-elseif strcmpi(analysis_type,'control'),
-
-	%Get u_g
-	u_g=results{3};
-	gsize=length(u_g);
-
-	%Used to recover velocities
-	indx=1:2:gsize;
-	indy=2:2:gsize;
-	indx=indx(part);
-	indy=indy(part);
-
-	%Recover velocity
-	md.cont_vx=u_g(indx)*md.yts;
-	md.cont_vy=u_g(indy)*md.yts;
-	md.cont_vel=sqrt(md.cont_vx.^2+md.cont_vy.^2);
-	
-	%recover parameter 
-	cont_parameter=results{4};
-	cont_parameter=cont_parameter(indx);
-	md.cont_parameter=cont_parameter;
-	
-	%read J
-	if(md.nsteps~=results{5}),
-		error('output from control method incompatible with model');
-	end
-	md.cont_J=results{6};
-
-elseif strcmpi(analysis_type,'thermal'),
-
-	%read results
-	time=results{4};
-	t_g=results{5};
-	m_g=results{6};
-
-	if strcmpi(md.sub_analysis_type,'steady');
-
-		%Recover temperature and melting
-		md.temperature=t_g(part);
-		md.melting=m_g(part)*md.yts;
-
-	else
-
-		%Recover temperature and melting
-		error('not implemented yet')
-
-	end
-
-elseif strcmpi(analysis_type,'prognostic'),
-
-	%read results
-	h_g=results{3};
-	md.new_thickness=h_g(part);
-
-else
-	error(['loadresultsfromdisk error message: unknow solution type ',analysis_type]);
-end
-
+md.results=parseresultsfromdisk(filename);
 
 %Check result is consistent
Index: /issm/trunk/src/m/classes/public/parseresultsfromdisk.m
===================================================================
--- /issm/trunk/src/m/classes/public/parseresultsfromdisk.m	(revision 642)
+++ /issm/trunk/src/m/classes/public/parseresultsfromdisk.m	(revision 643)
@@ -11,5 +11,5 @@
 end
 
-results={};
+results=struct();
 
 %Read fields until the end of the file.
@@ -17,8 +17,10 @@
 
 	result=ReadData(fid);
-	if isnan(result),
+	if isempty(result),
 		break;
 	else
-		results{end+1}=result;
+		eval(['results(' num2str(result.step) ').' result.fieldname '=result.field;']);
+		eval(['results(' num2str(result.step) ').time=result.time;']);
+		eval(['results(' num2str(result.step) ').step=result.step;']);
 	end
 end
Index: sm/trunk/src/m/solutions/dakota/qmu.m
===================================================================
--- /issm/trunk/src/m/solutions/dakota/qmu.m	(revision 642)
+++ 	(revision )
@@ -1,183 +1,0 @@
-function md=qmu(md,package,varargin)
-%INPUT function md=qmu(md,package)
-%Deal with coupled ISSM or Cielo/ Dakota runs, to do sensitivity analyses.
-
-global ISSM_DIR;
-
-% qmudir =['qmu_' datestr(now,'yyyymmdd_HHMMSS')];
-qmudir ='qmu';
-%  qmufile can not be changed unless cielo_ice_script.sh is also changed
-qmufile='qmu';
-ivar   =1;
-iresp  =1;
-imethod=1;
-iparams=1;
-runmpi =false;
-
-%  process any extra input arguments
-
-for i=1:2:nargin-2
-    switch varargin{i}
-        case 'qmudir'
-            qmudir =varargin{i+1};
-            disp(sprintf('qmudir =%s',qmudir));
-        case 'qmufile'
-            qmufile=varargin{i+1};
-            disp(sprintf('qmufile=%s',qmufile));
-        case 'ivar'
-            ivar   =varargin{i+1};
-            disp(sprintf('ivar   =%d',ivara));
-        case 'iresp'
-            iresp  =varargin{i+1};
-            disp(sprintf('iresp  =%d',iresp));
-        case 'imethod'
-            imethod=varargin{i+1};
-            disp(sprintf('imethod=%d',imethod));
-        case 'iparams'
-            iparams=varargin{i+1};
-            disp(sprintf('iparams=%d',iparams));
-        case 'overwrite'
-            overwrite=varargin{i+1};
-            disp(sprintf('overwrite=%s',overwrite));
-        case 'outfiles'
-            outfiles =varargin{i+1};
-            disp(sprintf('outfiles =%s',outfiles));
-        case 'rstfile'
-            rstfile  =varargin{i+1};
-            disp(sprintf('rstfile  =%s',rstfile));
-        case 'rundakota'
-            rundakota=varargin{i+1};
-            disp(sprintf('rundakota=%s',rundakota));
-        case 'runmpi'
-            runmpi =varargin{i+1};
-            disp(sprintf('runmpi =%d',runmpi));
-    end
-end
-
-%first create temporary directory in which we will work
-if exist(qmudir,'dir')
-    if ~exist('overwrite','var')
-        overwrite=input(['Overwrite existing ''' qmudir ''' directory? Y/N [N]: '], 's');
-    end
-    if strncmpi(overwrite,'y',1)
-        system(['rm -rf ' qmudir]);
-%    else
-%        error('Existing ''%s'' directory not overwritten.',qmudir);
-    end
-end
-mkdir(qmudir)
-cd(qmudir)
-system('cp $ISSM_DIR/startup.m .');
-
-%save our model in qmu so that it can be repeatedly used by Dakota.
-save('Qmu.model','md')
-
-%create m and in files for dakota
-if (~isfield(md.qmu_params(iparams),'direct') || ...
-    ~md.qmu_params(iparams).direct) && ...
-   (~isfield(md.qmu_params(iparams),'analysis_driver') || ...
-    isempty(md.qmu_params(iparams).analysis_driver))
-    md.qmu_params(iparams).analysis_driver=[ISSM_DIR '/src/m/solutions/dakota/cielo_ice_script.sh'];
-end
-dakota_in_data(md.qmu_method(imethod),md.variables(ivar),md.responses(iresp),md.qmu_params(iparams),qmufile,package,md);
-
-%  check for existence of results.out files to use
-if exist('results.out.1','file') || exist('results.out.zip','file')
-    if ~exist('outfiles','var')
-        outfiles=input(['Use existing ''results.out'' files? Y/N [N]: '], 's');
-    end
-    if ~strncmpi(outfiles,'y',1)
-        system('rm -f results.out.[1-9]*');
-    else
-        if exist('results.out.zip','file') && ~exist('results.out.1','file')
-            display('Inflating ''results.out.zip'' file.');
-            system('unzip -q results.out.zip');
-        end
-    end
-end
-
-%  check for existence of dakota.rst file to use
-rstflag='';
-if (~exist('outfiles','var') || ~strncmpi(outfiles,'y',1)) && ...
-    exist('dakota.rst','file')
-    if ~exist('rstfile','var')
-        rstfiles=input(['Use existing ''dakota.rst'' file? Y/N [N]: '], 's');
-    end
-    if strncmpi(rstfiles,'y',1)
-        system('rm -f results.out.[1-9]*');
-        system('dakota_restart_util print dakota.rst | grep completed');
-        rstflag=' -read_restart dakota.rst';
-    end
-end
-
-%call dakota
-if ~exist('rundakota','var')
-    rundakota=input(['Run Dakota analysis ''' qmufile '''? Y/N [N]: '], 's');
-end
-if ~strncmpi(rundakota,'y',1)
-    cd ..
-    return
-end
-
-if ~runmpi
-    system(['dakota -i ' qmufile '.in -o ' qmufile '.out -e ' qmufile '.err' rstflag]);
-else
-%  use 'mpd --ncpus=8 &' to initialize mpi and 'mpdringtest' to verify.
-%  exporting MPIRUN_NPROCS sets mpi in dakota.
-%    system('mpd --ncpus=8 &');
-%    system('mpdringtest');
-    system(['export MPIRUN_NPROCS=8;mpirun -np 4 dakota -i ' qmufile '.in -o ' qmufile '.out -e ' qmufile '.err' rstflag]);
-end
-
-%  check to see if dakota returned errors in the err file
-if exist([qmufile '.err'],'file')
-   fide=fopen([qmufile '.err'],'r');
-   fline=fgetl(fide);
-   if ischar(fline)
-       while ischar(fline)
-           disp(sprintf('%s',fline));
-           fline=fgetl(fide);
-       end
-       status=fclose(fide);
-       cd ../
-       error(['Dakota returned error in ''' qmufile '.err'' file.  ''' qmudir ''' directory retained.'])
-    end
-    status=fclose(fide);
-else
-   cd ../
-   error(['Dakota did not generate ''' qmufile '.err'' file.  ''' qmudir ''' directory retained.'])
-end
-
-%parse inputs and results from dakota
-[method,dvar,dresp_in]=dakota_in_parse([qmufile '.in']);
-md.dakotaresults.method   =method;
-md.dakotaresults.dvar     =dvar;
-md.dakotaresults.dresp_in =dresp_in;
-
-[method,dresp_out,scm,pcm,srcm,prcm]=dakota_out_parse([qmufile '.out']);
-md.dakotaresults.dresp_out=dresp_out;
-md.dakotaresults.scm      =scm;
-md.dakotaresults.pcm      =pcm;
-md.dakotaresults.srcm     =srcm;
-md.dakotaresults.prcm     =prcm;
-
-if exist('dakota_tabular.dat','file')
-    [method,dresp_dat                  ]=dakota_out_parse('dakota_tabular.dat');
-    md.dakotaresults.dresp_dat=dresp_dat;
-end
-
-%save input and output files into model
-%md.dakotain =readfile([qmufile '.in']);
-%md.dakotaout=readfile([qmufile '.out']);
-%if exist('dakota_tabular.dat','file')
-%	md.dakotadat=readfile('dakota_tabular.dat');
-%end
-
-%  move all the individual function evalutations into zip files
-system('zip -mq params.in.zip params.in.[1-9]*');
-system('zip -mq results.out.zip results.out.[1-9]*');
-system('zip -mq matlab.out.zip matlab*.out.[1-9]*');
-
-%get out of local directory and erase
-cd ../
-% system(['rm -rf ' qmudir]);
Index: /issm/trunk/src/m/utils/Nightly/testsgetanalysis.m
===================================================================
--- /issm/trunk/src/m/utils/Nightly/testsgetanalysis.m	(revision 642)
+++ /issm/trunk/src/m/utils/Nightly/testsgetanalysis.m	(revision 643)
@@ -15,9 +15,9 @@
 if strcmpi(string,'diagnostic'),
 	analysis_type='diagnostic';
-	sub_analysis_type='';
+	sub_analysis_type='none';
 
 elseif strcmpi(string,'prognostic'),
 	analysis_type='prognostic';
-	sub_analysis_type='';
+	sub_analysis_type='none';
 
 elseif strcmpi(string,'thermalsteady'),
@@ -31,5 +31,5 @@
 elseif strcmpi(string,'transient'),
 	analysis_type='transient';
-	sub_analysis_type='';
+	sub_analysis_type='none';
 
 else