Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 377)
@@ -154,4 +154,10 @@
 	param= new Param(count,"min_thermal_constraints",INTEGER);
 	param->SetInteger(model->min_thermal_constraints);
+	parameters->AddObject(param);
+
+	/*reconditioning_number: */
+	count++;
+	param= new Param(count,"reconditioning_number",DOUBLE);
+	param->SetDouble(model->reconditioning_number);
 	parameters->AddObject(param);
 
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 377)
@@ -99,4 +99,5 @@
 	int penta_thermal_steadystate;
 	double penta_viscosity_overshoot;
+	double penta_reconditioning_number;
 
 	/*matpar constructor input: */
@@ -413,5 +414,5 @@
 					penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
 					penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
-					penta_thermal_steadystate,penta_viscosity_overshoot); 
+					penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number); 
 
 			/*Add penta element to elements dataset: */
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp	(revision 377)
@@ -16,6 +16,17 @@
 void	CreateConstraintsDiagnosticStokes(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
 
+	int i;
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
 
-	DataSet* constraints = NULL;
+	/*spc intermediary data: */
+	int spc_sid;
+	int spc_node;
+	int spc_dof;
+	double spc_value;
+	int count;
+	
+	double* gridonstokes=NULL;
+
 
 	/*Create constraints: */
@@ -24,8 +35,62 @@
 	/*Now, is the flag ishutter on? otherwise, do nothing: */
 	if (!model->isstokes)goto cleanup_and_return;
+
+	/*Fetch data: */
+	ModelFetchData((void**)&gridonstokes,NULL,NULL,model_handle,"gridonstokes","Matrix","Mat");
+
+	count=0;
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<model->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((model->my_grids[i]==1)){
+	#endif
+
+		if ((int)!gridonstokes[i]){
 	
-	cleanup_and_return:	
+			/*This grid will see its vx,vy and vz dofs spc'd to pattyn velocities: */
+			spc_sid=count;
+			spc_node=i+1;
+			spc_dof=1; //we enforce first x translation degree of freedom
+			spc_value=0;
+
+			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
+			constraints->AddObject(spc);
+			count++;
+
+			spc_sid=count;
+			spc_node=i+1;
+			spc_dof=2; //we enforce first y translation degree of freedom
+			spc_value=0;
+			
+			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
+			constraints->AddObject(spc);
+			count++;
+			
+			spc_sid=count;
+			spc_node=i+1;
+			spc_dof=3; //we enforce first y translation degree of freedom
+			spc_value=0;
+			
+			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
+			constraints->AddObject(spc);
+			count++;
+
+		}
+
+	#ifdef _PARALLEL_
+	} //if((my_grids[i]==1))
+	#endif
+	}
+
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	constraints->Presort();
+
+	cleanup_and_return:
+	/*Free data: */
+	xfree((void**)&gridonstokes);
 	
 	/*Assign output pointer: */
 	*pconstraints=constraints;
-}
+}	
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 377)
@@ -18,4 +18,12 @@
 
 
+	
+	/*output: int* epart, int* my_grids, double* my_bordergrids*/
+
+
+	int i,j,k,n;
+	extern int my_rank;
+	extern int num_procs;
+
 	/*DataSets: */
 	DataSet*    elements  = NULL;
@@ -23,4 +31,97 @@
 	DataSet*    materials = NULL;
 	
+	/*Objects: */
+	Node*       node   = NULL;
+	Penta*      penta = NULL;
+	Matice*     matice  = NULL;
+	Matpar*     matpar  = NULL;
+
+	int         analysis_type;
+	
+	/*output: */
+	int* epart=NULL; //element partitioning.
+	int* npart=NULL; //node partitioning.
+	int* my_grids=NULL;
+	double* my_bordergrids=NULL;
+
+
+	/*intermediary: */
+	int elements_width; //size of elements
+	double B_avg;
+			
+	/*matice constructor input: */
+	int    matice_mid;
+	double matice_B;
+	double matice_n;
+	
+	/*penta constructor input: */
+
+	int penta_id;
+	int penta_mid;
+	int penta_mparid;
+	int penta_g[6];
+	double penta_h[6];
+	double penta_s[6];
+	double penta_b[6];
+	double penta_k[6];
+	int penta_friction_type;
+	double penta_p;
+	double penta_q;
+	int penta_shelf;
+	int penta_onbed;
+	int penta_onsurface;
+	double penta_meanvel;/*!scaling ratio for velocities*/
+	double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
+	int penta_collapse;
+	double penta_melting[6];
+	double penta_accumulation[6];
+	double penta_geothermalflux[6];
+	int penta_artdiff;
+	int penta_thermal_steadystate;
+	double penta_viscosity_overshoot;
+	double penta_reconditioning_number;
+
+	/*matpar constructor input: */
+	int	matpar_mid;
+	double matpar_rho_ice; 
+	double matpar_rho_water;
+	double matpar_heatcapacity;
+	double matpar_thermalconductivity;
+	double matpar_latentheat;
+	double matpar_beta;
+	double matpar_meltingpoint;
+	double matpar_mixed_layer_capacity;
+	double matpar_thermal_exchange_velocity;
+	double matpar_g;
+
+	/* node constructor input: */
+	int node_id;
+	int node_partitionborder=0;
+	double node_x[3];
+	int node_onbed;
+	int node_onsurface;
+	int node_upper_node_id;
+	int node_numdofs;
+
+		
+	#ifdef _PARALLEL_
+	/*Metis partitioning: */
+	int  range;
+	Vec  gridborder;
+	int  my_numgrids;
+	int* all_numgrids=NULL;
+	int  gridcount;
+	int  count;
+	#endif
+	int  first_grid_index;
+
+	/*Rifts:*/
+	int* riftsnumpenaltypairs;
+	double** riftspenaltypairs;
+	int* riftsfill;
+	double* riftsfriction;
+	double* riftpenaltypairs=NULL;
+	int el1,el2;
+	 
 	/*First create the elements, nodes and material properties: */
 	elements  = new DataSet(ElementsEnum());
@@ -31,4 +132,387 @@
 	if (!model->isstokes)goto cleanup_and_return;
 
+	/*Get analysis_type: */
+	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+
+	/*Width of elements: */
+	if(strcmp(model->meshtype,"2d")==0){
+		elements_width=3; //tria elements
+	}
+	else{
+		elements_width=6; //penta elements
+	}
+
+
+	#ifdef _PARALLEL_
+	/*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
+	if(strcmp(model->meshtype,"2d")==0){
+		/*load elements: */
+		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+	}
+	else{
+		/*load elements2d: */
+		ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
+	}
+
+
+	MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
+
+	/*Free elements and elements2d: */
+	xfree((void**)&model->elements);
+	xfree((void**)&model->elements2d);
+
+
+	/*Deal with rifts, they have to be included into one partition only, not several: */
+	FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts);
+	
+	for(i=0;i<model->numrifts;i++){
+		riftpenaltypairs=model->riftspenaltypairs[i];
+		for(j=0;j<model->riftsnumpenaltypairs[i];j++){
+			el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing
+			el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing
+			epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
+		}
+	}
+	/*Free rifts: */
+	xfree((void**)&riftsnumpenaltypairs); 
+	for(i=0;i<model->numrifts;i++){
+		double* temp=riftspenaltypairs[i];
+		xfree((void**)&temp);
+	}
+	xfree((void**)&riftspenaltypairs);
+	xfree((void**)&riftsfill);
+	xfree((void**)&riftsfriction);
+
+	/*Used later on: */
+	my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
+	#endif
+
+
+	/*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
+
+	/*2d mesh: */
+	if (strcmp(model->meshtype,"2d")==0){
+
+		throw ErrorException(__FUNCT__," stokes elements only supported in 3d!");
+
+	}
+	else{ //	if (strcmp(meshtype,"2d")==0)
+
+		/*Fetch data needed: */
+		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+		ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
+		ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
+		ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
+		ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
+		ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
+		ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
+		ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
+		ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
+		ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
+		ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
+		ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
+		ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
+		ModelFetchData((void**)&model->accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat");
+		ModelFetchData((void**)&model->melting,NULL,NULL,model_handle,"melting","Matrix","Mat");
+		
+		for (i=0;i<model->numberofelements;i++){
+		#ifdef _PARALLEL_
+		/*We are using our element partition to decide which elements will be created on this node: */
+		if(my_rank==epart[i]){
+		#endif
+
+			
+			/*name and id: */
+			penta_id=i+1; //matlab indexing.
+			penta_mid=i+1; //refers to the corresponding material property card
+			penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
+
+			/*vertices,thickness,surface,bed and drag: */
+			for(j=0;j<6;j++){
+				penta_g[j]=(int)*(model->elements+elements_width*i+j);
+				penta_h[j]=*(model->thickness+    ((int)*(model->elements+elements_width*i+j)-1)); 
+				penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1)); 
+				penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1)); 
+				penta_k[j]=*(model->drag+        ((int)*(model->elements+elements_width*i+j)-1)); 
+				penta_melting[j]=*(model->melting+        ((int)*(model->elements+elements_width*i+j)-1));
+				penta_accumulation[j]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+j)-1));
+			}
+
+			/*basal drag:*/
+			penta_friction_type=(int)model->drag_type;
+	
+			penta_p=model->p[i];
+			penta_q=model->q[i];
+
+			/*diverse: */
+			penta_shelf=(int)*(model->elementoniceshelf+i);
+			penta_onbed=(int)*(model->elementonbed+i);
+			penta_onsurface=(int)*(model->elementonsurface+i);
+			penta_meanvel=model->meanvel;
+			penta_epsvel=model->epsvel;
+			
+			/*viscosity_overshoot*/
+			penta_viscosity_overshoot=model->viscosity_overshoot;
+
+			/*reconditioning_number: */
+			penta_reconditioning_number=model->reconditioning_number;
+
+			
+			if (*(model->elements_type+2*i+1)==StokesEnum()){
+			
+				/*Create Penta using its constructor:*/
+				penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
+						penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
+						penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
+						penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number); 
+
+				/*Add penta element to elements dataset: */
+				elements->AddObject(penta);
+			}
+	
+
+			/*Deal with material:*/
+			matice_mid=i+1; //same as the material id from the geom2 elements.
+			/*Average B over 6 element grids: */
+			B_avg=0;
+			for(j=0;j<6;j++){
+				B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
+			}
+			B_avg=B_avg/6;
+			matice_B= B_avg;
+			matice_n=(double)*(model->n+i);
+			
+			if (*(model->elements_type+2*i+1)==StokesEnum()){
+	
+				/*Create matice using its constructor:*/
+				matice= new Matice(matice_mid,matice_B,matice_n);
+		
+				/*Add matice element to materials dataset: */
+				materials->AddObject(matice);
+			}
+			
+			#ifdef _PARALLEL_
+			/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
+			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
+			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
+			 will hold which grids belong to this partition*/
+			my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
+			my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
+			my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
+			my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
+			my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
+			my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
+			#endif
+
+		#ifdef _PARALLEL_
+		}//if(my_rank==epart[i])
+		#endif
+
+		}//for (i=0;i<numberofelements;i++)
+
+		/*Free data: */
+		xfree((void**)&model->elements);
+		xfree((void**)&model->thickness);
+		xfree((void**)&model->surface);
+		xfree((void**)&model->bed);
+		xfree((void**)&model->drag);
+		xfree((void**)&model->p);
+		xfree((void**)&model->q);
+		xfree((void**)&model->elementoniceshelf);
+		xfree((void**)&model->elementonbed);
+		xfree((void**)&model->elementonsurface);
+		xfree((void**)&model->elements_type);
+		xfree((void**)&model->n);
+		xfree((void**)&model->B);
+
+	} //if (strcmp(meshtype,"2d")==0)
+
+	/*Add one constant material property to materials: */
+	matpar_mid=model->numberofelements+1; //put it at the end of the materials
+	matpar_g=model->g; 
+	matpar_rho_ice=model->rho_ice; 
+	matpar_rho_water=model->rho_water; 
+	matpar_thermalconductivity=model->thermalconductivity; 
+	matpar_heatcapacity=model->heatcapacity; 
+	matpar_latentheat=model->latentheat; 
+	matpar_beta=model->beta; 
+	matpar_meltingpoint=model->meltingpoint; 
+	matpar_mixed_layer_capacity=model->mixed_layer_capacity; 
+	matpar_thermal_exchange_velocity=model->thermal_exchange_velocity; 
+
+	/*Create matpar object using its constructor: */
+	matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
+			matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
+			matpar_thermal_exchange_velocity,matpar_g);
+		
+	/*Add to materials datset: */
+	materials->AddObject(matpar);
+	
+	#ifdef _PARALLEL_
+		/*From the element partitioning, we can determine which grids are on the inside of this cpu's 
+		 *element partition, and which are on its border with other nodes:*/
+		gridborder=NewVec(model->numberofnodes);
+
+		for (i=0;i<model->numberofnodes;i++){
+			if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
+		}
+		VecAssemblyBegin(gridborder);
+		VecAssemblyEnd(gridborder);
+
+		#ifdef _DEBUG_
+		VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
+		#endif
+		
+		VecToMPISerial(&my_bordergrids,gridborder);
+
+		#ifdef _DEBUG_
+		if(my_rank==0){
+			for (i=0;i<model->numberofnodes;i++){
+				printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
+			}
+		}
+		#endif
+	#endif
+
+	/*Partition penalties in 3d: */
+	if(strcmp(model->meshtype,"3d")==0){
+	
+		/*Get penalties: */
+		ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
+
+		if(model->numpenalties){
+
+			model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
+			#ifdef _SERIAL_
+			for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=1;
+			#else
+			for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;
+
+			for(i=0;i<model->numpenalties;i++){
+				first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
+				if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
+					/*All grids that are being penalised belong to this node's internal grid partition.:*/
+					model->penaltypartitioning[i]=1;
+				}
+				if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
+					model->penaltypartitioning[i]=0;
+				}
+			}
+			#endif
+		}
+
+		/*Free penalties: */
+		xfree((void**)&model->penalties);
+	}
+
+	/*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 
+	 We can therefore determine  which grids are internal to this node's partition 
+	 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 
+	 that, go and create the grids*/
+
+	/*Create nodes from x,y,z, as well as the spc values on those grids: */
+		
+	/*First fetch data: */
+	if (strcmp(model->meshtype,"3d")==0){
+		ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
+		ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
+	}
+	ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
+	ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
+	ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonstokes,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->borderstokes,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+
+
+	/*Get number of dofs per node: */
+	DistributeNumDofs(&node_numdofs,analysis_type);
+
+	for (i=0;i<model->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((my_grids[i]==1)){
+	#endif
+
+		node_id=i+1; //matlab indexing
+			
+		
+		
+		#ifdef _PARALLEL_
+		if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
+			node_partitionborder=1;
+		}
+		else{
+			node_partitionborder=0;
+		}
+		#else
+			node_partitionborder=0;
+		#endif
+
+
+		node_x[0]=model->x[i];
+		node_x[1]=model->y[i];
+		node_x[2]=model->z[i];
+
+		
+		node_onbed=(int)model->gridonbed[i];
+		node_onsurface=(int)model->gridonsurface[i];	
+		if (strcmp(model->meshtype,"3d")==0){
+			if (isnan(model->uppernodes[i])){
+				node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
+			}
+			else{
+				node_upper_node_id=(int)model->uppernodes[i];
+			}
+		}
+		else{
+			/*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
+			node_upper_node_id=node_id;
+		}
+
+		/*Create node using its constructor: */
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+
+		/*set single point constraints.: */
+		/*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */
+		if (model->gridonstokes[i]){
+			for(k=1;k<=node_numdofs;k++){
+				node->FreezeDof(k);
+			}
+		}
+		if (model->borderstokes[i]){ 
+			//freeze everything except pressure
+			node->FreezeDof(1);
+			node->FreezeDof(2);
+			node->FreezeDof(3);
+		}
+
+
+		/*Add node to nodes dataset: */
+		nodes->AddObject(node);
+
+	#ifdef _PARALLEL_
+	} //if((my_grids[i]==1))
+	#endif
+	}
+
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	elements->Presort();
+	nodes->Presort();
+	materials->Presort();
+
+	/*Clean fetched data: */
+	xfree((void**)&model->deadgrids);
+	xfree((void**)&model->x);
+	xfree((void**)&model->y);
+	xfree((void**)&model->z);
+	xfree((void**)&model->gridonbed);
+	xfree((void**)&model->gridonsurface);
+	xfree((void**)&model->uppernodes);
+	xfree((void**)&model->gridonstokes);
+	xfree((void**)&model->borderstokes);
+		
 	cleanup_and_return:
 
@@ -37,3 +521,16 @@
 	*pnodes=nodes;
 	*pmaterials=materials;
+
+	/*Keep partitioning information into model*/
+	model->epart=epart;
+	model->my_grids=my_grids;
+	model->my_bordergrids=my_bordergrids;
+
+	/*Free ressources:*/
+	#ifdef _PARALLEL_
+	xfree((void**)&all_numgrids);
+	xfree((void**)&npart);
+	VecFree(&gridborder);
+	#endif
+
 }
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 377)
@@ -16,6 +16,57 @@
 void	CreateLoadsDiagnosticStokes(DataSet** ploads, Model* model,ConstDataHandle model_handle){
 
+	int i,j,counter;
+	int element;
+
+	extern int my_rank;
+	extern int num_procs;
+	
 	DataSet*    loads    = NULL;
-
+	Icefront*   icefront = NULL;
+	Penpair*    penpair  = NULL;
+
+	int segment_width;
+	int element_type;
+	int i1,i2,i3,i4;
+	int i0,range;
+	int grid1,grid2;
+
+	/*rift penpair: */
+	double  normal[2];
+	double  length;
+	double  friction;
+	int     fill;
+		
+	/*icefront intermediary data: */
+	char icefront_type[ICEFRONTSTRING];
+	int icefront_element_type;
+	int	icefront_sid;
+	int icefront_eid;
+	int icefront_mparid;
+	int	icefront_node_ids[MAX_ICEFRONT_GRIDS];
+	double icefront_h[MAX_ICEFRONT_GRIDS];
+	double	icefront_b[MAX_ICEFRONT_GRIDS];
+
+	/*penpair intermediary data: */
+	int penpair_id;
+	int penpair_node_ids[2];
+	double penpair_penalty_offset;
+	int penpair_numdofs;
+	int penpair_dof;
+	int penpair_penalty_lock;
+	int penpair_element_ids[2];
+	double penpair_friction;
+	int penpair_fill;
+	double penpair_normal[2];
+	double penpair_length;
+
+	/*Rifts:*/
+	int* riftsnumpenaltypairs=NULL;
+	double** riftspenaltypairs=NULL;
+	int* riftsfill=NULL;
+	int* riftsfriction=NULL;
+	double* riftpenaltypairs=NULL;
+	int el1,el2;
+	
 	/*Create loads: */
 	loads   = new DataSet(LoadsEnum());
@@ -24,6 +75,207 @@
 	if (!model->isstokes)goto cleanup_and_return;
 
+	/*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids
+	 * referenced by a certain load must belong to the cluster node): */
+	ModelFetchData((void**)&model->segmentonneumann_diag,&model->numberofsegs_diag,NULL,model_handle,"segmentonneumann_diag","Matrix","Mat");
+	ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
+	ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
+	ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
+
+	/*First load data:*/
+	for (i=0;i<model->numberofsegs_diag;i++){
+		
+		if (strcmp(model->meshtype,"2d")==0){
+			segment_width=3;
+			element_type=TriaEnum();
+		}
+		else{
+			segment_width=5;
+			element_type=PentaEnum();
+		}
+
+
+		element=(int)(*(model->segmentonneumann_diag+segment_width*i+segment_width-1)-1); //element is in the last column
+
+		#ifdef _PARALLEL_
+		if (model->epart[element]!=my_rank){
+			/*This load does not belong to this cluster node, as it references an element 
+			 *that does not belong to this node's partition. Drop this 'i':*/
+			continue;
+		}
+		#endif
+	
+		icefront_mparid=model->numberofelements+1; //matlab indexing
+		icefront_sid=i+1; //matlab indexing
+		icefront_eid=(int)*(model->segmentonneumann_diag+segment_width*i+segment_width-1); //matlab indexing
+		icefront_element_type=element_type;
+
+		i1=(int)*(model->segmentonneumann_diag+segment_width*i+0);
+		i2=(int)*(model->segmentonneumann_diag+segment_width*i+1);
+			
+		icefront_node_ids[0]=i1;
+		icefront_node_ids[1]=i2;
+		
+		icefront_h[0]=model->thickness[i1-1];
+		icefront_h[1]=model->thickness[i2-1];
+
+		icefront_b[0]=model->bed[i1-1];
+		icefront_b[1]=model->bed[i2-1];
+		
+		if ((int)*(model->elements_type+2*element+0)==MacAyealEnum()){ //this is a collapsed 3d element, icefront will be 2d
+			strcpy(icefront_type,"segment");
+		}
+		else if ((int)*(model->elements_type+2*element+0)==PattynEnum()){ //this is a real 3d element, icefront will be 3d.
+			strcpy(icefront_type,"quad");
+			i3=(int)*(model->segmentonneumann_diag+segment_width*i+2);
+			i4=(int)*(model->segmentonneumann_diag+segment_width*i+3);
+			icefront_node_ids[2]=i3;
+			icefront_node_ids[3]=i4;
+			
+			icefront_h[2]=model->thickness[i3-1];
+			icefront_h[3]=model->thickness[i4-1];
+
+			icefront_b[2]=model->bed[i3-1];
+			icefront_b[3]=model->bed[i4-1];
+		}
+		else{
+			throw ErrorException(__FUNCT__," element type not supported yet");
+		}
+
+		icefront = new Icefront(icefront_type,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b);
+		
+		loads->AddObject(icefront);
+
+	}
+	/*Free data: */
+	xfree((void**)&model->segmentonneumann_diag);
+	xfree((void**)&model->elements_type);
+	xfree((void**)&model->thickness);
+	xfree((void**)&model->bed);
+
+		
+	/*Create penalties loads for linking of collapsed formulation grids with non collapsed grids: */
+
+	/*First fetch penalties: */
+	if (strcmp(model->meshtype,"3d")==0){
+		ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
+	}
+
+	counter=0;
+	if(model->numpenalties){
+
+		/*First deal with internal grids: */
+		for (i=0;i<model->numpenalties;i++){
+			if (model->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU
+
+				for(j=0;j<model->numlayers-1;j++){
+
+					/*We are pairing grids along a vertical profile.*/
+					grid1=(int)*(model->penalties+model->numlayers*i+j);
+					grid2=(int)*(model->penalties+model->numlayers*i+j+1);
+
+					penpair_id=counter+1; //matlab indexing
+					penpair_dof=1;
+					penpair_node_ids[0]=grid1;
+					penpair_node_ids[1]=grid2;
+					penpair_penalty_offset=model->penalty_offset;
+					penpair_numdofs=1;
+
+					penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
+								penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
+					loads->AddObject(penpair);
+
+
+					counter++;
+
+					penpair_id=counter+1; //matlab indexing
+					penpair_dof=2;
+					penpair_node_ids[0]=grid1;
+					penpair_node_ids[1]=grid2;
+					penpair_penalty_offset=model->penalty_offset;
+					penpair_numdofs=1;
+
+					penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
+								penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
+					loads->AddObject(penpair);
+
+
+					counter++;
+
+				} //for ( i=0; i<numpenalties; i++) 
+			}
+		}
+	}
+
+	/*Free data: */
+	xfree((void**)&model->penalties);
+	
+	/*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */
+	/*First fetch rifts: */
+	if(model->numrifts){
+
+		for(i=0;i<model->numrifts;i++){
+			riftpenaltypairs=model->riftspenaltypairs[i];
+			for(j=0;j<model->riftsnumpenaltypairs[i];j++){
+
+				el1=(int)*(riftpenaltypairs+7*j+2);
+				#ifdef _PARALLEL_
+				if (model->epart[el1-1]!=my_rank){
+					/*This load does not belong to this cluster node, as it references an element 
+					 *that does not belong to this node's partition. Drop this 'j':*/
+					continue;
+				}
+				#endif
+
+				/*Ok, retrieve all the data needed to add a penalty between the two grids: */
+				el2=(int)*(riftpenaltypairs+7*j+3); 
+				grid1=(int)*(riftpenaltypairs+7*j+0); 
+				grid2=(int)*(riftpenaltypairs+7*j+1);
+				normal[0]=*(riftpenaltypairs+7*j+4);
+				normal[1]=*(riftpenaltypairs+7*j+5);
+				length=*(riftpenaltypairs+7*j+6);
+				friction=model->riftsfriction[i];
+				fill=model->riftsfill[i];
+
+				penpair_id=counter+1; //matlab indexing
+				penpair_node_ids[0]=grid1;
+				penpair_node_ids[1]=grid2;
+				penpair_element_ids[0]=el1;
+				penpair_element_ids[1]=el2;
+				penpair_penalty_offset=model->penalty_offset;
+				penpair_penalty_lock=model->penalty_lock;
+				penpair_numdofs=2;
+				penpair_normal[0]=normal[0];
+				penpair_normal[1]=normal[1];
+				penpair_length=length;
+				penpair_friction=friction;
+				penpair_fill=fill;
+
+				penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
+						penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
+				
+				loads->AddObject(penpair);
+
+				counter++;
+			}
+		}
+	}
+
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	loads->Presort();
+
 	cleanup_and_return:
-	
+
+	/*Free ressources:*/
+	xfree((void**)&riftsnumpenaltypairs); 
+	for(i=0;i<model->numrifts;i++){
+		double* temp=riftspenaltypairs[i];
+		xfree((void**)&temp);
+	}
+	xfree((void**)&riftspenaltypairs);
+	xfree((void**)&riftsfill);
+	xfree((void**)&riftsfriction);
+
+
 	/*Assign output pointer: */
 	*ploads=loads;
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 377)
@@ -79,4 +79,5 @@
 	int penta_thermal_steadystate;
 	double penta_viscosity_overshoot;
+	double penta_reconditioning_number;
 
 	/*matpar constructor input: */
@@ -238,5 +239,5 @@
 				penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
 				penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
-				penta_thermal_steadystate,penta_viscosity_overshoot); 
+				penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number); 
 
 		/*Add penta element to elements dataset: */
Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Model.cpp	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/Model.cpp	(revision 377)
@@ -61,4 +61,6 @@
 	model->gridonbed=NULL;
 	model->gridonsurface=NULL;
+	model->gridonstokes=NULL;
+	model->borderstokes=NULL;
 	model->thickness=NULL;
 	model->surface=NULL;
@@ -124,4 +126,5 @@
 	model->yts=0;
 	model->viscosity_overshoot=0;
+	model->reconditioning_number=0;
 	model->waitonlock=0;
 
@@ -195,4 +198,6 @@
 	xfree((void**)&model->gridonbed);
 	xfree((void**)&model->gridonsurface);
+	xfree((void**)&model->gridonstokes);
+	xfree((void**)&model->borderstokes);
 	xfree((void**)&model->thickness);
 	xfree((void**)&model->surface);
@@ -328,4 +333,5 @@
 	ModelFetchData((void**)&model->solverstring,NULL,NULL,model_handle,"solverstring","String",NULL);
 	ModelFetchData((void**)&model->viscosity_overshoot,NULL,NULL,model_handle,"viscosity_overshoot","Scalar",NULL);
+	ModelFetchData((void**)&model->reconditioning_number,NULL,NULL,model_handle,"stokesreconditioning","Scalar",NULL);
 	ModelFetchData((void**)&model->waitonlock,NULL,NULL,model_handle,"waitonlock","Integer",NULL);
 
Index: /issm/trunk/src/c/ModelProcessorx/Model.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 377)
@@ -56,4 +56,6 @@
 	double* gridonbed;
 	double* gridonsurface;
+	double* gridonstokes;
+	double* borderstokes;
 	double* thickness;
 	double* surface;
@@ -119,4 +121,5 @@
 	double  yts;
 	double  viscosity_overshoot;
+	double  reconditioning_number;
 	int     waitonlock;
 
Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 376)
+++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 377)
@@ -90,4 +90,5 @@
 	int penta_thermal_steadystate;
 	double penta_viscosity_overshoot;
+	double penta_reconditioning_number;
 
 	/* node constructor input: */
@@ -286,5 +287,5 @@
 					penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
 					penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
-					penta_thermal_steadystate,penta_viscosity_overshoot); 
+					penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number); 
 
 			/*Add penta element to elements dataset: */
@@ -446,5 +447,5 @@
 		if (strcmp(model->meshtype,"3d")==0){
 			/*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
-			if (model->gridonbed[i]==0){
+			if (model->deadgrids[i]){
 				for(k=1;k<=node_numdofs;k++){
 					node->FreezeDof(k);
Index: /issm/trunk/src/c/objects/OptArgs.h
===================================================================
--- /issm/trunk/src/c/objects/OptArgs.h	(revision 376)
+++ /issm/trunk/src/c/objects/OptArgs.h	(revision 377)
@@ -14,8 +14,10 @@
 	char* function_name;
 	mxArray* m;
+	mxArray* inputs;
 	mxArray* p_g;
 	mxArray* u_g_obs;
 	mxArray* grad_g;
 	mxArray* n;
+	mxArray* analysis_type;
 
 };
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 376)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 377)
@@ -22,5 +22,5 @@
 				double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel, 
 				int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], 
-				int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot){
+				int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_reconditioning_number){
 	
 	int i;
@@ -59,4 +59,5 @@
 	thermal_steadystate = penta_thermal_steadystate;
 	viscosity_overshoot = penta_viscosity_overshoot;
+	reconditioning_number = penta_reconditioning_number;
 
 	return;
@@ -99,4 +100,5 @@
 	printf("   thermal_steadystate: %i\n",thermal_steadystate);
 	printf("   viscosity_overshoot: %i\n",viscosity_overshoot);
+	printf("   reconditioning_number: %i\n",reconditioning_number);
 	return;
 }
@@ -146,4 +148,5 @@
 	memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
 	memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
+	memcpy(marshalled_dataset,&reconditioning_number,sizeof(reconditioning_number));marshalled_dataset+=sizeof(reconditioning_number);
 	
 	*pmarshalled_dataset=marshalled_dataset;
@@ -182,4 +185,5 @@
 		sizeof(thermal_steadystate) +
 		sizeof(viscosity_overshoot) +
+		sizeof(reconditioning_number) +
 		sizeof(int); //sizeof(int) for enum type
 }
@@ -229,4 +233,5 @@
 	memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
 	memcpy(&viscosity_overshoot,marshalled_dataset,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
+	memcpy(&reconditioning_number,marshalled_dataset,sizeof(reconditioning_number));marshalled_dataset+=sizeof(reconditioning_number);
 
 	/*nodes and materials are not pointing to correct objects anymore:*/
Index: /issm/trunk/src/c/objects/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Penta.h	(revision 376)
+++ /issm/trunk/src/c/objects/Penta.h	(revision 377)
@@ -53,4 +53,5 @@
 		int    thermal_steadystate;
 		double viscosity_overshoot;
+		double reconditioning_number;
 	
 	public:
@@ -60,5 +61,5 @@
 				double p, double q, int    shelf, int    onbed, int    onsurface, double meanvel,double epsvel, 
 				int    collapse, double melting[6], double accumulation[6], double geothermalflux[6], 
-				int    artdiff, int    thermal_steadystate,double viscosity_overshoot);
+				int    artdiff, int    thermal_steadystate,double viscosity_overshoot,double reconditioning_number);
 		~Penta();
 
Index: /issm/trunk/src/c/shared/Numerics/OptFunc.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/OptFunc.cpp	(revision 376)
+++ /issm/trunk/src/c/shared/Numerics/OptFunc.cpp	(revision 377)
@@ -23,5 +23,5 @@
 	double J;
 
-	mxArray*   inputs[6];
+	mxArray*   inputs[8];
 	mxArray*   psearch_scalar=NULL;
 	mxArray*   mxJ=NULL;
@@ -32,10 +32,12 @@
 	
 	inputs[1]=optargs->m;
-	inputs[2]=optargs->p_g;
-	inputs[3]=optargs->u_g_obs;
-	inputs[4]=optargs->grad_g;
-	inputs[5]=optargs->n;
+	inputs[2]=optargs->inputs;
+	inputs[3]=optargs->p_g;
+	inputs[4]=optargs->u_g_obs;
+	inputs[5]=optargs->grad_g;
+	inputs[6]=optargs->n;
+	inputs[7]=optargs->analysis_type;
 
-	mexCallMATLAB( 1, &mxJ, 6, (mxArray**)inputs, optargs->function_name);
+	mexCallMATLAB( 1, &mxJ, 8, (mxArray**)inputs, optargs->function_name);
 
 	/*extract misfit from mxArray*/
Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 376)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 377)
@@ -59,4 +59,5 @@
 	md.gridonpattyn=NaN;
 	md.gridonstokes=NaN;
+	md.borderstokes=NaN;
 
 	%Stokes mesh
Index: /issm/trunk/src/m/classes/public/marshall.m
===================================================================
--- /issm/trunk/src/m/classes/public/marshall.m	(revision 376)
+++ /issm/trunk/src/m/classes/public/marshall.m	(revision 377)
@@ -50,4 +50,6 @@
 WriteData(fid,md.gridonbed,'Mat','gridonbed');
 WriteData(fid,md.gridonsurface,'Mat','gridonsurface');
+WriteData(fid,md.gridonstokes,'Mat','gridonstokes');
+WriteData(fid,md.borderstokes,'Mat','borderstokes');
 
 
@@ -126,4 +128,5 @@
 WriteData(fid,md.solverstring,'String','solverstring');
 WriteData(fid,md.viscosity_overshoot,'Scalar','viscosity_overshoot');
+WriteData(fid,md.stokesreconditioning,'Scalar','reconditioning_number');
 WriteData(fid,md.waitonlock,'Integer','waitonlock');
 
Index: /issm/trunk/src/m/classes/public/setelementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/setelementstype.m	(revision 376)
+++ /issm/trunk/src/m/classes/public/setelementstype.m	(revision 377)
@@ -120,4 +120,10 @@
 md.deadgrids=deadgrids;
 
+%figure out the border stokes grids
+stokes_elements=find(md.elements_type(:,2)==stokesenum()); %find the elements on the stokes domain
+borderflags=zeros(md.numberofgrids,1); 
+borderflags(md.elements(stokes_elements,:))=1; %find all the grids of the elements on stokes domain, ie stokes grids and borderstokes
+md.borderstokes=borderflags-md.gridonstokes; %remove stokes grids from this list
+
 %figure out solution types
 md.ishutter=double(any(md.elements_type(:,1)==hutterenum()));
Index: /issm/trunk/src/m/solutions/cielo/GradJCompute.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/GradJCompute.m	(revision 376)
+++ /issm/trunk/src/m/solutions/cielo/GradJCompute.m	(revision 377)
@@ -1,3 +1,3 @@
-function grad_g=GradJCompute(m,inputs, u_g_obs);
+function grad_g=GradJCompute(m,inputs, u_g_obs,analysis_type);
 
 %Recover solution for this stiffness and right hand side: 
@@ -5,5 +5,5 @@
 	disp('         computing velocities...')
 end
-[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,'diagnostic_horiz');
+[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type);
 
 %Buid Du, difference between observed velocity and model velocity.
Index: /issm/trunk/src/m/solutions/cielo/control.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/control.m	(revision 376)
+++ /issm/trunk/src/m/solutions/cielo/control.m	(revision 377)
@@ -19,4 +19,8 @@
 	options=ControlOptions(m.parameters);
 
+	%initialize inputs, ie m.nparameters on which we invert.
+	inputs=inputlist;
+	inputs=add(inputs,'velocity',m.parameters.u_g,'doublevec',3,m.parameters.numberofnodes);
+
 	for n=1:m.parameters.nsteps,
 		
@@ -26,8 +30,6 @@
 		disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(m.parameters.nsteps)]));
 
-		%initialize inputs, ie m.nparameters on which we invert.
-		inputs=inputlist;
+		%update inputs with new parameter and fit
 		inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes);
-		inputs=add(inputs,'velocity',m.parameters.u_g,'doublevec',3,m.parameters.numberofnodes);
 		inputs=add(inputs,'fit',m.parameters.fit(n),'double');
 
@@ -36,5 +38,5 @@
 
 		disp('      computing gradJ...');
-		c(n).grad_g=GradJCompute(m,inputs,u_g_obs);
+		c(n).grad_g=GradJCompute(m,inputs,u_g_obs,md.analysis_type);
 		disp('      done.');
 
@@ -53,6 +55,6 @@
 		
 		disp('      optimizing along gradient direction...'); 
-		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,p_g,u_g_obs,c(n).grad_g,n);
-		%[search_scalar c(n).J]=fminbnd('objectivefunctionC',0,1,options,m,p_g,u_g_obs,c(n).grad_g,n);
+		md.analysis_type
+		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,inputs,p_g,u_g_obs,c(n).grad_g,n,md.analysis_type);
 		disp('      done.');
 
@@ -72,5 +74,5 @@
 		%some temporary saving 
 		if(mod(n,5)==0),
-			solution=controlfinalsol(c,m,p_g);
+			solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type);
 			save temporary_control_results solution
 		end
@@ -81,5 +83,5 @@
 	%Create final solution
 	disp('      preparing final velocity solution...');
-	solution=controlfinalsol(c,m,p_g);
+	solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type);
 	disp('      done.');
 	
Index: /issm/trunk/src/m/solutions/cielo/controlfinalsol.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/controlfinalsol.m	(revision 376)
+++ /issm/trunk/src/m/solutions/cielo/controlfinalsol.m	(revision 377)
@@ -1,9 +1,7 @@
-function solution=controlfinalsol(c,m,p_g);
+function solution=controlfinalsol(c,m,p_g,inputs,analysis_type);
 
 %From parameters, build inputs for icediagnostic_core, using the final parameters
-inputs=inputlist;
-inputs=add(inputs,'velocity',m.parameters.u_g,'doublevec',3,m.parameters.numberofnodes);
 inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes);
-u_g=diagnostic_core_nonlinear(m,inputs,'diagnostic_horiz');
+u_g=diagnostic_core_nonlinear(m,inputs,analysis_type);
 
 %Build partitioning vectors to recover solution
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 376)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 377)
@@ -72,5 +72,5 @@
 
 		%"recondition" pressure 
-		p_g=p_g/m_ds.parameters.stokesreconditioning;
+		p_g=p_g/m_ds.parameters.reconditioning_number;
 
 		displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
@@ -95,5 +95,5 @@
 	
 		%"decondition" pressure
-		p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning;
+		p_g=u_g(4:4:end)*m_dh.parameters.reconditioning_number;
 	end
 end
Index: /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m	(revision 376)
+++ /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m	(revision 377)
@@ -1,3 +1,3 @@
-function J =objectivefunctionC(search_scalar,m,p_g,u_g_obs,grad_g,n);
+function J =objectivefunctionC(search_scalar,m,inputs,p_g,u_g_obs,grad_g,n,analysis_type);
         
 %recover some parameters
@@ -10,12 +10,9 @@
 
 %Plug parameter into inputs
-inputs=inputlist;
-inputs=add(inputs,'velocity',m.parameters.u_g,'doublevec',3,m.parameters.numberofnodes);
 inputs=add(inputs,m.parameters.control_type,parameter,'doublevec',2,m.parameters.numberofnodes);
 
 %Run diagnostic with updated parameters.
-u_g=diagnostic_core_nonlinear(m,inputs,'diagnostic_horiz');
+u_g=diagnostic_core_nonlinear(m,inputs,analysis_type);
 
 %Compute misfit for this velocity field. 
-inputs=add(inputs,'fit',m.parameters.fit(n),'double');
 J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, u_g_obs,inputs);
Index: /issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp
===================================================================
--- /issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp	(revision 376)
+++ /issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp	(revision 377)
@@ -45,8 +45,10 @@
 	optargs.function_name=function_name;
 	optargs.m=MODEL;
-	optargs.p_g=PG;
+	optargs.p_g=PG; 
+	optargs.inputs=INPUTS;
 	optargs.u_g_obs=VELOBS;
 	optargs.grad_g=GRADIENT;
 	optargs.n=STEP;
+	optargs.analysis_type=ANALYSIS;
 
 	optpars.xmin=xmin;
@@ -71,5 +73,5 @@
 {
 	_printf_("\n");
-	_printf_("   usage: [search_scalar J] = %s(function_name,xmin,xmax,options,m,p_g,u_g_obs,grad_g,step)\n");
+	_printf_("   usage: [search_scalar J] = %s(function_name,xmin,xmax,options,m,inputs,p_g,u_g_obs,grad_g,step,analysis_type)\n",__FUNCT__);
 	_printf_("\n");
 }
Index: /issm/trunk/src/mex/ControlOptimization/ControlOptimization.h
===================================================================
--- /issm/trunk/src/mex/ControlOptimization/ControlOptimization.h	(revision 376)
+++ /issm/trunk/src/mex/ControlOptimization/ControlOptimization.h	(revision 377)
@@ -22,8 +22,10 @@
 #define OPTIONS (mxArray*)prhs[3]
 #define MODEL (mxArray*)prhs[4]
-#define PG (mxArray*)prhs[5]
-#define VELOBS (mxArray*)prhs[6]
-#define GRADIENT (mxArray*)prhs[7]
-#define STEP (mxArray*)prhs[8]
+#define INPUTS (mxArray*)prhs[5]
+#define PG (mxArray*)prhs[6]
+#define VELOBS (mxArray*)prhs[7]
+#define GRADIENT (mxArray*)prhs[8]
+#define STEP (mxArray*)prhs[9]
+#define ANALYSIS (mxArray*)prhs[10]
 
 /* serial output macros: */
@@ -35,5 +37,5 @@
 #define NLHS  2
 #undef NRHS
-#define NRHS  9
+#define NRHS  11