Index: /issm/trunk/src/c/DataSet/DataSet.cpp
===================================================================
--- /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 201)
+++ /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 202)
@@ -146,7 +146,12 @@
 }
 
-int  DataSet::DeleteObject(int id){
-
-	return 0;
+int  DataSet::DeleteObject(Object* object){
+	
+	vector<Object*>::iterator iterator;
+
+	iterator = find(objects.begin(), objects.end(),object);
+
+	objects.erase(iterator);
+
 }
 
@@ -293,4 +298,30 @@
 }
 
+Object*   DataSet::FindParamObject(char* name){
+
+	/*Go through a dataset, and find a Param* object 
+	 *which parameter name is "name" : */
+	
+	vector<Object*>::iterator object;
+	Param* param=NULL;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+		/*Find param type objects: */
+		if((*object)->Enum()==ParamEnum()){
+
+			/*Ok, this object is a parameter, recover it and ask which name it has: */
+			param=(Param*)(*object);
+
+			if (strcmp(param->GetParameterName(),name)==0){
+				/*Ok, this is the one! Return the object: */
+				return (*object);
+			}
+		}
+	}
+	return NULL;
+}
+
+
 void   DataSet::NodeRank(int* ranks){
 
@@ -324,7 +355,11 @@
 
 	int  dofcount=0;
+	int  dofcount1=0;
 	int* alldofcount=NULL;
+	int* alldofcount1=NULL;
 	int* borderdofs=NULL;
+	int* borderdofs1=NULL;
 	int* allborderdofs=NULL;
+	int* allborderdofs1=NULL;
 	int  i;
 	vector<Object*>::iterator object;
@@ -343,5 +378,5 @@
 			
 			/*Ok, this object is a node, ask it to distribute dofs, and update the dofcount: */
-			node->DistributeDofs(&dofcount);
+			node->DistributeDofs(&dofcount,&dofcount1);
 			
 		}
@@ -356,13 +391,20 @@
 	MPI_Bcast(alldofcount,num_procs,MPI_INT,0,MPI_COMM_WORLD);
 
+	alldofcount1=(int*)xmalloc(num_procs*sizeof(int));
+	MPI_Gather(&dofcount1,1,MPI_INT,alldofcount1,1,MPI_INT,0,MPI_COMM_WORLD);
+	MPI_Bcast(alldofcount1,num_procs,MPI_INT,0,MPI_COMM_WORLD);
+
 	/*Ok, now every cpu should start its own dof count at the end of the dofcount 
 	 * from cpu-1. : */
 	dofcount=0;
+	dofcount1=0;
 	if(my_rank==0){
 		dofcount=0;
+		dofcount1=0;
 	}
 	else{
 		for(i=0;i<my_rank;i++){
 			dofcount+=alldofcount[i];
+			dofcount1+=alldofcount1[i];
 		}
 	}
@@ -378,5 +420,5 @@
 			
 			/*Ok, this object is a node, ask it to update his dofs: */
-			node->UpdateDofs(dofcount);
+			node->UpdateDofs(dofcount,dofcount1);
 			
 		}
@@ -387,4 +429,7 @@
 	borderdofs=(int*)xcalloc(numberofnodes*numdofspernode,sizeof(int));
 	allborderdofs=(int*)xcalloc(numberofnodes*numdofspernode,sizeof(int));
+	borderdofs1=(int*)xcalloc(numberofnodes*3,sizeof(int));
+	allborderdofs1=(int*)xcalloc(numberofnodes*3,sizeof(int));
+
 	for ( object=objects.begin() ; object < objects.end(); object++ ){
 
@@ -395,9 +440,10 @@
 			
 			/*Ok, let this object show its border dofs, if is is a border dof: */
-			node->ShowBorderDofs(borderdofs);
+			node->ShowBorderDofs(borderdofs,borderdofs1);
 			
 		}
 	}
 	MPI_Allreduce ( (void*)borderdofs,(void*)allborderdofs,numberofnodes*numdofspernode,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
+	MPI_Allreduce ( (void*)borderdofs1,(void*)allborderdofs1,numberofnodes*3,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
 
 	/*Ok, now every cpu knows everyone else's border node dofs, update the border dofs accordingly: */
@@ -410,5 +456,5 @@
 			
 			/*Ok, this object is a node, ask it to update his dofs: */
-			node->UpdateBorderDofs(allborderdofs);
+			node->UpdateBorderDofs(allborderdofs,allborderdofs1);
 			
 		}
@@ -420,4 +466,8 @@
 	xfree((void**)&borderdofs);
 	xfree((void**)&allborderdofs);
+	xfree((void**)&alldofcount1);
+	xfree((void**)&borderdofs1);
+	xfree((void**)&allborderdofs1);
+
 	return;
 
@@ -433,10 +483,9 @@
 
 	/*Create partition vector: */
-	partition=NewVec(numberofnodes*numdofspernode);
-
-	/*Go through all nodes, and ask each node to plug its doflist in 
+	partition=NewVec(numberofnodes);
+
+	/*Go through all nodes, and ask each node to plug its 1D doflist in 
 	 * partition. The location where each node plugs its doflist into 
-	 * partition is determined by its (id-1)*numdofspernode (ie, serial 
-	 * organisation of the dofs).
+	 * partition is determined by its (id-1)*3 (ie, serial * organisation of the dofs).
 	 */
 	
Index: /issm/trunk/src/c/DataSet/DataSet.h
===================================================================
--- /issm/trunk/src/c/DataSet/DataSet.h	(revision 201)
+++ /issm/trunk/src/c/DataSet/DataSet.h	(revision 202)
@@ -45,4 +45,5 @@
 		int   Size();
 		int   FindParam(void* pvalue, char* name);
+		Object* FindParamObject(char* name);
 		void  NodeRank(int* ranks);
 		void  DistributeDofs(int numberofnodes,int numdofspernode);
@@ -76,4 +77,5 @@
 		void  Misfit(double* pJ, double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type);
 		void  VelocityExtrude(Vec ug,double* ug_serial);
+		int   DeleteObject(Object* object);
 
 
Index: /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 201)
+++ /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 202)
@@ -25,5 +25,7 @@
 	double* optscal=NULL;
 	double* maxiter=NULL; 
-	double* control_parameter=NULL;
+	Vec     control_parameter=NULL;
+	Vec     vx_obs=NULL;
+	Vec     vy_obs=NULL;
 
 	/*Get counter : */
@@ -98,29 +100,29 @@
 
 	/*Get vx_obs, vy_obs, and the parameter value: */
-	ModelFetchData((void**)&model->vx_obs,NULL,NULL,model_handle,"vx_obs","Matrix","Mat");
-	ModelFetchData((void**)&model->vy_obs,NULL,NULL,model_handle,"vy_obs","Matrix","Mat");
-	ModelFetchData((void**)&control_parameter,NULL,NULL,model_handle,model->control_type,"Matrix","Mat");
+	ModelFetchData((void**)&vx_obs,NULL,NULL,model_handle,"vx_obs","Vector",NULL);
+	ModelFetchData((void**)&vy_obs,NULL,NULL,model_handle,"vy_obs","Vector",NULL);
+	ModelFetchData((void**)&control_parameter,NULL,NULL,model_handle,model->control_type,"Vector",NULL);
 
-	for(i=0;i<model->numberofnodes;i++)model->vx_obs[i]=model->vx_obs[i]/model->yts;
-	for(i=0;i<model->numberofnodes;i++)model->vy_obs[i]=model->vy_obs[i]/model->yts;
+	VecScale(vx_obs,1.0/model->yts);
+	VecScale(vy_obs,1.0/model->yts);
 
 	count++;
-	param= new Param(count,"vx_obs",DOUBLEVEC);
-	param->SetDoubleVec(model->vx_obs,model->numberofnodes);
+	param= new Param(count,"vx_obs",PETSCVEC);
+	param->SetVec(vx_obs);
 	parameters->AddObject(param);
 
 	count++;
-	param= new Param(count,"vy_obs",DOUBLEVEC);
-	param->SetDoubleVec(model->vy_obs,model->numberofnodes);
+	param= new Param(count,"vy_obs",PETSCVEC);
+	param->SetVec(vy_obs);
 	parameters->AddObject(param);
 
 	count++;
-	param= new Param(count,"control_parameter",DOUBLEVEC);
-	param->SetDoubleVec(control_parameter,model->numberofnodes);
+	param= new Param(count,"control_parameter",PETSCVEC);
+	param->SetVec(control_parameter);
 	parameters->AddObject(param);
 	
-	xfree((void**)&model->vx_obs);
-	xfree((void**)&model->vy_obs);
-	xfree((void**)&control_parameter);
+	VecFree(&vx_obs);
+	VecFree(&vy_obs);
+	VecFree(&control_parameter);
 
 
Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 201)
+++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 202)
@@ -27,6 +27,7 @@
 	double* maxiter=NULL; 
 	double* control_parameter=NULL;
-	
-	
+	Vec     vx=NULL;
+	Vec     vy=NULL;
+	Vec     vz=NULL;
 	
 	/*Initialize dataset: */
@@ -168,5 +169,4 @@
 	param->SetInteger(numberofdofspernode);
 	parameters->AddObject(param);
-
 
 	/*Deal with control parameters: */
@@ -175,4 +175,33 @@
 	}
 
+	/*Diagnostic: */
+	/*Get vx and vy: */
+	ModelFetchData((void**)&vx,NULL,NULL,model_handle,"vx","Vector",NULL);
+	ModelFetchData((void**)&vy,NULL,NULL,model_handle,"vy","Vector",NULL);
+	ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Vector",NULL);
+
+	if(vx)VecScale(vx,1.0/model->yts);
+	if(vy)VecScale(vy,1.0/model->yts);
+	if(vz)VecScale(vz,1.0/model->yts);
+
+	count++;
+	param= new Param(count,"vx",PETSCVEC);
+	param->SetVec(vx);
+	parameters->AddObject(param);
+
+	count++;
+	param= new Param(count,"vy",PETSCVEC);
+	param->SetVec(vy);
+	parameters->AddObject(param);
+
+	count++;
+	param= new Param(count,"vz",PETSCVEC);
+	param->SetVec(vz);
+	parameters->AddObject(param);
+
+	VecFree(&vx);
+	VecFree(&vy);
+	VecFree(&vz);
+
 	/*All our datasets are already ordered by ids. Set presort flag so that later on, when sorting is requested on these 
 	 * datasets, it will not be redone: */
Index: /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp
===================================================================
--- /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp	(revision 201)
+++ /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp	(revision 202)
@@ -19,36 +19,115 @@
 	double* partition=NULL;
 	Param*  param=NULL;
+
+	/*diagnostic: */
+	Vec     vec_vx=NULL;
+	Vec     vec_vy=NULL;
+	Vec     vec_vz=NULL;
+
+	double* vx=NULL;
+	double* vy=NULL;
+	double* vz=NULL;
+	
+	double* u_g=NULL;
+	
+	/*control: */
+	Vec     vec_vx_obs=NULL;
+	Vec     vec_vy_obs=NULL;
+	Vec     vec_control_parameter=NULL;
+
 	double* vx_obs=NULL;
 	double* vy_obs=NULL;
 	double* control_parameter=NULL;
+
+	double* u_g_obs=NULL;
+	double* p_g=NULL;
+	
+
+	/*diverse: */
 	int     numberofnodes;
 	int     analysis_type;
 	int     count;
 
-	/*new parameters: */
-	double* u_g_obs=NULL;
-	double* p_g=NULL;
-
 	parameters->FindParam((void*)&analysis_type,"analysis_type");
 	count=parameters->Size();
+		
+	parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+
+	/*First serialize partition vector: */
+	VecToMPISerial(&partition,part);
+	
+	
+	if (   (analysis_type==ControlAnalysisEnum()) || 
+	       (analysis_type==DiagnosticHorizAnalysisEnum()) || 
+	       (analysis_type==DiagnosticVertAnalysisEnum()) || 
+	       (analysis_type==DiagnosticStokesAnalysisEnum()) 
+	   ){
+
+		parameters->FindParam((void*)&vec_vx,"vx");
+		parameters->FindParam((void*)&vec_vy,"vy");
+		parameters->FindParam((void*)&vec_vz,"vz");
+
+		if(vec_vx)VecToMPISerial(&vx,vec_vx);
+		if(vec_vy)VecToMPISerial(&vy,vec_vy);
+		if(vec_vz)VecToMPISerial(&vz,vec_vz);
+			
+		u_g=(double*)xcalloc(numberofnodes*3,sizeof(double));
+
+		if(vx){
+			for(i=0;i<numberofnodes;i++){
+				u_g[(int)(3*partition[i]+0)]=vx[i];  
+			}
+		}
+
+		if(vy){
+			for(i=0;i<numberofnodes;i++){
+				u_g[(int)(3*partition[i]+1)]=vy[i];  
+			}
+		}
+
+		if(vz){
+			for(i=0;i<numberofnodes;i++){
+				u_g[(int)(3*partition[i]+2)]=vz[i];  
+			}
+		}
+
+
+		/*Now, create new parameters: */
+		count++;
+		param= new Param(count,"u_g",DOUBLEVEC);
+		param->SetDoubleVec(u_g,3*numberofnodes);
+		parameters->AddObject(param);
+
+		/*erase old parameters: */
+		param=(Param*)parameters->FindParamObject("vx");
+		parameters->DeleteObject((Object*)param);
+
+		param=(Param*)parameters->FindParamObject("vy");
+		parameters->DeleteObject((Object*)param);
+
+		param=(Param*)parameters->FindParamObject("vz");
+		parameters->DeleteObject((Object*)param);
+
+	}
+
 
 	if(analysis_type==ControlAnalysisEnum()){
 
-		/*First serialize partition vector: */
-		VecToMPISerial(&partition,part);
-		
-		parameters->FindParam((void*)&vx_obs,"vx_obs");
-		parameters->FindParam((void*)&vy_obs,"vy_obs");
-		parameters->FindParam((void*)&control_parameter,"control_parameter");
-		parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+		parameters->FindParam((void*)&vec_vx_obs,"vx_obs");
+		parameters->FindParam((void*)&vec_vy_obs,"vy_obs");
+		parameters->FindParam((void*)&vec_control_parameter,"control_parameter");
 
 		/*Now, from vx_obs and vy_obs, build u_g_obs, correctly partitioned: */
+		VecToMPISerial(&vx_obs,vec_vx_obs);
+		VecToMPISerial(&vy_obs,vec_vy_obs);
+		VecToMPISerial(&control_parameter,vec_control_parameter);
+
 		u_g_obs=(double*)xcalloc(numberofnodes*2,sizeof(double));
 		p_g=(double*)xcalloc(numberofnodes*2,sizeof(double));
 
 		for(i=0;i<numberofnodes;i++){
-			u_g_obs[(int)(partition[2*i+0])]=vx_obs[i];  
-			u_g_obs[(int)(partition[2*i+1])]=vy_obs[i];  
-			p_g[(int)(partition[2*i+0])]=control_parameter[i];
+			u_g_obs[(int)(2*partition[i]+0)]=vx_obs[i];  
+			u_g_obs[(int)(2*partition[i]+1)]=vy_obs[i];  
+			p_g[(int)(2*partition[i]+0)]=control_parameter[i];
 		}
 
@@ -61,13 +140,37 @@
 		count++;
 		param= new Param(count,"p_g",DOUBLEVEC);
-		param->SetDoubleVec(p_g,2*numberofnodes);
+		param->SetDoubleVec(p_g,numberofnodes);
 		parameters->AddObject(param);
+
+		/*erase old parameter: */
+		param=(Param*)parameters->FindParamObject("vx_obs");
+		parameters->DeleteObject((Object*)param);
+
+		param=(Param*)parameters->FindParamObject("vy_obs");
+		parameters->DeleteObject((Object*)param);
+
+		param=(Param*)parameters->FindParamObject("control_parameter");
+		parameters->DeleteObject((Object*)param);
+
 	}
 
 	xfree((void**)&partition);
+	
+	xfree((void**)&vx);
+	xfree((void**)&vy);
+	xfree((void**)&vz);
+	VecFree(&vec_vx);
+	VecFree(&vec_vy);
+	VecFree(&vec_vz);
+	xfree((void**)&u_g);
+
 	xfree((void**)&vx_obs);
 	xfree((void**)&vy_obs);
+	VecFree(&vec_vx_obs);
+	VecFree(&vec_vy_obs);
+	VecFree(&vec_control_parameter);
 	xfree((void**)&control_parameter);
 	xfree((void**)&u_g_obs);
 	xfree((void**)&p_g); 
+
 }
Index: /issm/trunk/src/c/objects/Node.cpp
===================================================================
--- /issm/trunk/src/c/objects/Node.cpp	(revision 201)
+++ /issm/trunk/src/c/objects/Node.cpp	(revision 202)
@@ -73,4 +73,6 @@
 		printf("%i|",doflist[i]);
 	}
+	printf("   doflist1:|");
+	printf("%i|",doflist1[1]);
 	printf("\n");
 
@@ -109,4 +111,5 @@
 	memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
 	memcpy(marshalled_dataset,&doflist,sizeof(doflist));marshalled_dataset+=sizeof(doflist);
+	memcpy(marshalled_dataset,&doflist1,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1);
 	memcpy(marshalled_dataset,&mset,sizeof(mset));marshalled_dataset+=sizeof(mset);
 	memcpy(marshalled_dataset,&nset,sizeof(nset));marshalled_dataset+=sizeof(nset);
@@ -131,4 +134,5 @@
 		sizeof(onsurface)+
 		sizeof(doflist)+
+		sizeof(doflist1)+
 		sizeof(mset)+
 		sizeof(nset)+
@@ -164,4 +168,5 @@
 	memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
 	memcpy(&doflist,marshalled_dataset,sizeof(doflist));marshalled_dataset+=sizeof(doflist);
+	memcpy(&doflist1,marshalled_dataset,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1);
 	memcpy(&mset,marshalled_dataset,sizeof(mset));marshalled_dataset+=sizeof(mset);
 	memcpy(&nset,marshalled_dataset,sizeof(nset));marshalled_dataset+=sizeof(nset);
@@ -193,11 +198,13 @@
 }
 
-void  Node::DistributeDofs(int* pdofcount){ 
+void  Node::DistributeDofs(int* pdofcount,int* pdofcount1){
 
 	int i;
 	extern int my_rank;
 	int dofcount;
+	int dofcount1;
 
 	dofcount=*pdofcount;
+	dofcount1=*pdofcount1;
 	
 	if(clone){
@@ -212,11 +219,15 @@
 	dofcount+=numberofdofs;
 
+	doflist1[0]=dofcount1;
+	dofcount1+=1;
+
 	/*Assign output pointers: */
 	*pdofcount=dofcount;
+	*pdofcount1=dofcount1;
 
 	return; 
 }
 
-void  Node::UpdateDofs(int dofcount){
+void  Node::UpdateDofs(int dofcount,int dofcount1){
 	
 	int i;
@@ -232,9 +243,10 @@
 		doflist[i]+=dofcount;
 	}
+	doflist1[0]+=dofcount1;
 
 	return; 
 }
 
-void  Node::ShowBorderDofs(int* borderdofs){
+void  Node::ShowBorderDofs(int* borderdofs,int* borderdofs1){
 
 	int j;
@@ -252,9 +264,10 @@
 		*(borderdofs+numberofdofs*(id-1)+j)=doflist[j];
 	}
-
-	return;
-}
-
-void  Node::UpdateBorderDofs(int* allborderdofs){ 
+	*(borderdofs1+(id-1)+0)=doflist1[0];
+
+	return;
+}
+
+void  Node::UpdateBorderDofs(int* allborderdofs,int* allborderdofs1){ 
 
 	int j;
@@ -272,28 +285,20 @@
 	for(j=0;j<numberofdofs;j++){
 		doflist[j]=*(allborderdofs+numberofdofs*(id-1)+j);
-	}	
+	}
+	doflist1[0]=*(allborderdofs1+(id-1)+0);
 	return; 
 }
 void  Node::CreatePartition(Vec partition){ 
 
-	int      i;
-	int*     idxm=NULL;
-	double*  values=NULL;
-
-	idxm=(int*)xmalloc(numberofdofs*sizeof(int));
-	values=(double*)xmalloc(numberofdofs*sizeof(double));
-
-	for(i=0;i<numberofdofs;i++)idxm[i]=(id-1)*numberofdofs+i;
-	for(i=0;i<numberofdofs;i++)values[i]=(double)doflist[i];
-
-	VecSetValues(partition,numberofdofs,idxm,values,INSERT_VALUES);
-
-	xfree((void**)&idxm);
-	xfree((void**)&values);
-
-	return;
-}
-void  Node::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
-		
+	int      idxm;
+	double   value;
+
+	idxm=(id-1);
+	value=(double)doflist1[0];
+
+	VecSetValues(partition,1,&idxm,&value,INSERT_VALUES);
+
+	return;
+}
 		
 void  Node::SetClone(int* minranks){
Index: /issm/trunk/src/c/objects/Node.h
===================================================================
--- /issm/trunk/src/c/objects/Node.h	(revision 201)
+++ /issm/trunk/src/c/objects/Node.h	(revision 202)
@@ -36,5 +36,6 @@
 
 		/*data that is post processed : */
-		int doflist[MAXDOFSPERNODE];
+		int doflist[MAXDOFSPERNODE]; //dof list on which we solve
+		int doflist1[1]; //1d dof list to recover input parameters
 
 	public:
@@ -52,10 +53,9 @@
 		int   GetId(void); 
 		int   MyRank(void);
-		void  DistributeDofs(int* pdofcount);
-		void  UpdateDofs(int dofcount);
-		void  ShowBorderDofs(int* borderdofs);
-		void  UpdateBorderDofs(int* allborderdofs);
+		void  DistributeDofs(int* pdofcount,int* pdofcount1);
+		void  UpdateDofs(int dofcount,int dofcount1);
+		void  ShowBorderDofs(int* borderdofs,int* borderdofs1);
+		void  UpdateBorderDofs(int* allborderdofs,int* allborderdofs1);
 		void  CreatePartition(Vec partition);
-		void  DistributeNumDofs(int* numdofspernode,int analysis_type);
 		void  SetClone(int* minranks);
 		int   GetNumberOfDofs();
Index: /issm/trunk/src/c/objects/Param.cpp
===================================================================
--- /issm/trunk/src/c/objects/Param.cpp	(revision 201)
+++ /issm/trunk/src/c/objects/Param.cpp	(revision 202)
@@ -37,5 +37,35 @@
 
 Param::~Param(){
-	return;
+
+	switch(type){
+	
+		case STRING:
+			break;
+	
+		case INTEGER:
+			break;
+	
+		case DOUBLE:
+			break;
+		
+		case DOUBLEVEC:
+			xfree((void**)&doublevec);
+			break;
+	
+		case DOUBLEMAT:
+			xfree((void**)&doublemat);
+			break;
+
+		case PETSCVEC:
+			VecFree(&vec);
+			break;
+
+		case  PETSCMAT:
+			MatFree(&mat);
+			break;
+
+		default:
+			throw ErrorException(__FUNCT__,exprintf("%s%i","unknow parameter type ",type));
+	}
 }
 
@@ -139,18 +169,32 @@
 
 		case PETSCVEC:
-			VecGetSize(vec,&M);
-			VecToMPISerial(&serial_vec,vec);
-			memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-			memcpy(marshalled_dataset,serial_vec,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
-			xfree((void**)&serial_vec);
+			if(vec){
+				VecGetSize(vec,&M);
+				VecToMPISerial(&serial_vec,vec);
+				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+				memcpy(marshalled_dataset,serial_vec,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
+				xfree((void**)&serial_vec);
+			}
+			else{
+				M=0;
+				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+			}
 			break;
 
 		case PETSCMAT:
-			MatGetSize(mat,&M,&N);
-			MatToSerial(&serial_mat,mat);
-			memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
-			memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
-			memcpy(marshalled_dataset,serial_mat,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
-			xfree((void**)&serial_mat);
+			if(mat){
+				MatGetSize(mat,&M,&N);
+				MatToSerial(&serial_mat,mat);
+				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+				memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
+				memcpy(marshalled_dataset,serial_mat,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
+				xfree((void**)&serial_mat);
+			}
+			else{
+				M=0;
+				N=0;
+				memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
+				memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
+			}
 			break;
 		
@@ -263,17 +307,22 @@
 		case PETSCVEC:
 			memcpy(&M,marshalled_dataset,sizeof(M));marshalled_dataset+=sizeof(M);
-			serial_vec=(double*)xmalloc(M*sizeof(double));
-			memcpy(serial_vec,marshalled_dataset,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
-
-			vec=NewVec(M);
-			idxm=(int*)xmalloc(M*sizeof(int));
-			for(i=0;i<M;i++)idxm[i]=i;
-			VecSetValues(vec,M,idxm,serial_vec,INSERT_VALUES);
-
-			VecAssemblyBegin(vec);
-			VecAssemblyEnd(vec);
-			
-			xfree((void**)&serial_vec);
-			xfree((void**)&idxm);
+			if(M){
+				serial_vec=(double*)xmalloc(M*sizeof(double));
+				memcpy(serial_vec,marshalled_dataset,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
+
+				vec=NewVec(M);
+				idxm=(int*)xmalloc(M*sizeof(int));
+				for(i=0;i<M;i++)idxm[i]=i;
+				VecSetValues(vec,M,idxm,serial_vec,INSERT_VALUES);
+
+				VecAssemblyBegin(vec);
+				VecAssemblyEnd(vec);
+				
+				xfree((void**)&serial_vec);
+				xfree((void**)&idxm);
+			}
+			else{
+				vec=NULL;
+			}
 			break;
 
@@ -281,18 +330,23 @@
 			memcpy(&M,marshalled_dataset,sizeof(M));marshalled_dataset+=sizeof(M);
 			memcpy(&N,marshalled_dataset,sizeof(N));marshalled_dataset+=sizeof(N);
-			serial_mat=(double*)xmalloc(M*N*sizeof(double));
-			memcpy(serial_mat,marshalled_dataset,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
-			mat=NewMat(M,N,&sparsity,NULL,NULL);
-			idxm=(int*)xmalloc(M*sizeof(int));
-			idxn=(int*)xmalloc(N*sizeof(int));
-			for(i=0;i<M;i++)idxm[i]=i;
-			for(i=0;i<N;i++)idxn[i]=i;
-			MatSetValues(mat,M,idxm,N,idxn,serial_mat,INSERT_VALUES); 
-			MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 
-			MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
-
-			xfree((void**)&serial_mat);
-			xfree((void**)&idxm);
-			xfree((void**)&idxn);
+			if(M!=0 && N!=0){
+				serial_mat=(double*)xmalloc(M*N*sizeof(double));
+				memcpy(serial_mat,marshalled_dataset,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
+				mat=NewMat(M,N,&sparsity,NULL,NULL);
+				idxm=(int*)xmalloc(M*sizeof(int));
+				idxn=(int*)xmalloc(N*sizeof(int));
+				for(i=0;i<M;i++)idxm[i]=i;
+				for(i=0;i<N;i++)idxn[i]=i;
+				MatSetValues(mat,M,idxm,N,idxn,serial_mat,INSERT_VALUES); 
+				MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 
+				MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
+
+				xfree((void**)&serial_mat);
+				xfree((void**)&idxm);
+				xfree((void**)&idxn);
+			}
+			else{
+				mat=NULL;
+			}
 			break;
 
@@ -350,6 +404,8 @@
 		Vec* pvec=(Vec*)pvalue;
 		Vec  outvec=NULL;
-		VecDuplicate(vec,&outvec);
-		VecCopy(vec,outvec);
+		if(vec){
+			VecDuplicate(vec,&outvec);
+			VecCopy(vec,outvec);
+		}
 		*pvec=outvec;
 	}
@@ -357,5 +413,7 @@
 		Mat* pmat=(Mat*)pvalue;
 		Mat  outmat=NULL;
-		MatDuplicate(mat,MAT_COPY_VALUES,&outmat);
+		if(mat){
+			MatDuplicate(mat,MAT_COPY_VALUES,&outmat);
+		}
 		*pmat=outmat;
 	}
@@ -424,2 +482,20 @@
 
 }
+
+#undef __FUNCT__
+#define __FUNCT__ "SetVec"
+void  Param::SetVec(Vec value){
+
+	if(value){
+		VecDuplicate(value,&vec);
+		VecCopy(value,vec);
+		VecGetSize(value,&M);
+		N=1;
+	}
+	else{
+		vec=NULL;
+		M=0;
+		N=0;
+	}
+
+}
Index: /issm/trunk/src/c/objects/Param.h
===================================================================
--- /issm/trunk/src/c/objects/Param.h	(revision 201)
+++ /issm/trunk/src/c/objects/Param.h	(revision 202)
@@ -46,4 +46,5 @@
 		void  SetDouble(double value);
 		void  SetDoubleVec(double* value,int size);
+		void  SetVec(Vec value);
 		void  SetInteger(int value);
 		void  SetString(char* value);
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 201)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 202)
@@ -663,5 +663,11 @@
 		GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
 
-		DL_scalar=alpha2*gauss_weight*Jdet;
+		if (velocity_param){
+			DL_scalar=alpha2*gauss_weight*Jdet;
+		}
+		else{
+			DL_scalar=0;
+		}
+			
 		for (i=0;i<2;i++){
 			DL[i][i]=DL_scalar;
@@ -916,5 +922,4 @@
 	B_param=ParameterInputsRecover(inputs,"B");
 	if(temperature_average_param){
-		temperature_average-0;
 		for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
 		temperature_average=temperature_average/3.0;
@@ -2122,5 +2127,4 @@
 					&Ke_gg_gaussian[0][0],0);
 
-		
 		/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
 		for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
@@ -2352,6 +2356,6 @@
 		/*Get bed slope: */
 		GetParameterDerivativeValue(&slope[0], &b[0],&xyz_list[0][0], gauss_l1l2l3);
-		dbdx=-slope[0];
-		dbdy=-slope[1];
+		dbdx=slope[0];
+		dbdy=slope[1];
 
 		/* Get Jacobian determinant: */
Index: /issm/trunk/src/m/classes/public/setelementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/setelementstype.m	(revision 201)
+++ /issm/trunk/src/m/classes/public/setelementstype.m	(revision 202)
@@ -119,3 +119,9 @@
 deadgrids(find(md.gridonbed))=0;                      %grids from elements on bed are not dead
 md.deadgrids=deadgrids;
+
+%figure out solution types
+md.ishutter=any(md.elements_type(:,1)==hutterenum());
+md.ismacayealpattyn=any(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum() );
+md.isstokes=any(md.elements_type(:,2));
+
 end
Index: /issm/trunk/src/m/solutions/cielo/CreateFemModel.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/CreateFemModel.m	(revision 201)
+++ /issm/trunk/src/m/solutions/cielo/CreateFemModel.m	(revision 202)
@@ -13,5 +13,5 @@
 	disp('   generating degrees of freedom...');
 	[m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,parameters);
-	
+
 	disp('   generating single point constraints...');
 	[m.nodes,m.yg]=SpcNodes(m.nodes,m.constraints);
@@ -31,5 +31,5 @@
 	disp('   configuring element and loads...');
 	[m.elements,m.loads,m.nodes] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials);
-	
+
 	disp('   processing parameters...');
 	parameters= ProcessParams(parameters,m.part);
Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 201)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 202)
@@ -18,18 +18,15 @@
 	md.dof=m_dh.nodesets.fsize; %biggest dof number
 
+	%initialize inputs
+	if strcmp(md.analysis_type,'diagnostic_stokes'),
+		inputs.velocity=m_dh.parameters.u_g;
+	else
+		inputs.velocity=zeros(2*m_dh.parameters.numberofnodes,1);
+		inputs.velocity(1:2:end)=m_dh.parameters.u_g(1:3:end);
+		inputs.velocity(2:2:end)=m_dh.parameters.u_g(2:3:end);
+	end
+
 	%Get horizontal solution. 
 	disp(sprintf('\n%s',['computing horizontal velocities...']));
-
-	%plug existing velocity in inputs
-	if ~isnan(md.vx) & ~isnan(md.vy),
-		indx=m_dh.part(1:2:end);
-		indy=m_dh.part(2:2:end);
-		velocity=zeros(m_dh.nodesets.gsize,1);
-		velocity(indx)=md.vx/md.yts;
-		velocity(indy)=md.vy/md.yts;
-		inputs=struct('velocity',velocity);
-	else
-		inputs={};
-	end
 
 	u_g=diagnostic_core_nonlinear(m_dh,inputs);
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m	(revision 201)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m	(revision 202)
@@ -11,4 +11,5 @@
 	count=1;
 	converged=0;
+
 	[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
 	if velocity_is_present
@@ -18,5 +19,4 @@
 	end
 	soln(count).u_f=[];
-
 
 	if m.parameters.debug,
@@ -74,5 +74,8 @@
 	%   Figure out if convergence is reached.
 		if(count>=3 | velocity_is_present),
-			dug=soln(count).u_g-soln(count-1).u_g; nduinf=norm(dug,inf)*m.parameters.yts; ndu=norm(dug,2); nu=norm(soln(count-1).u_g,2); 
+			dug=soln(count).u_g-soln(count-1).u_g; 
+			nduinf=norm(dug,inf)*m.parameters.yts; 
+			ndu=norm(dug,2); 
+			nu=norm(soln(count-1).u_g,2); 
 
 			%Relative criterion
Index: sm/trunk/src/m/solutions/ice/DiagnosticSolutionType.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/DiagnosticSolutionType.m	(revision 201)
+++ 	(revision )
@@ -1,15 +1,0 @@
-function [ishutter,ismacayealpattyn,isstokes]=DiagnosticSolutionType(elements_type)
-%DIAGNOSTICSOLUTIONTYPE - figure out which element are present
-%
-%   This routine goes through the types of element and check if there is any
-%   Hutter, MacAyeal, Pattyn or Stokes element, and return a list of flag
-%
-%   Usage:
-%      [ishutter,ismacayealpattyn,isstokes]=DiagnosticSolutionType(elements_type)
-%
-%   See also: SETELEMENTSTYPE
-
-ishutter=any(elements_type(:,1)==hutterenum());
-ismacayealpattyn=any(elements_type(:,1)==macayealenum() | elements_type(:,1)==pattynenum() );
-isstokes=any(elements_type(:,2));
-
Index: /issm/trunk/src/m/solutions/ice/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/diagnostic.m	(revision 201)
+++ /issm/trunk/src/m/solutions/ice/diagnostic.m	(revision 202)
@@ -17,6 +17,9 @@
 %Create fem structure (input of diagnostic3d)
 fem=struct();
+
 %Figure out which type of elements are present
-[fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type);
+fem.ishutter=ishutter;
+fem.ismacayealpattyn=ismacayealpattyn;
+fem.isstokes=isstokes;
 
 if strcmpi(md.type,'2d'),
Index: /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp
===================================================================
--- /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp	(revision 201)
+++ /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp	(revision 202)
@@ -36,9 +36,10 @@
 	/*Create constraints: */
 	CreateConstraints(&constraints,model,MODEL);
-	
+
 	/*Create loads: */
 	CreateLoads(&loads,model,MODEL);
-	
+
 	/*Create parameters: */
+
 	CreateParameters(&parameters,model,MODEL);