Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 7390)
+++ /issm/trunk/src/c/Makefile.am	(revision 7391)
@@ -1139,4 +1139,5 @@
 					./modules/Solverx/Solverx.cpp\
 					./modules/Solverx/Solverx.h\
+					./modules/Solverx/DofTypesToIndexSet.cpp\
 					./modules/Scotchx/Scotchx.cpp\
 					./modules/Scotchx/Scotchx.h\
Index: /issm/trunk/src/c/modules/Solverx/DofTypesToIndexSet.cpp
===================================================================
--- /issm/trunk/src/c/modules/Solverx/DofTypesToIndexSet.cpp	(revision 7391)
+++ /issm/trunk/src/c/modules/Solverx/DofTypesToIndexSet.cpp	(revision 7391)
@@ -0,0 +1,78 @@
+/*!\file Solverx
+ * \brief solver
+ */
+
+#include "./Solverx.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+        
+void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){
+
+	/*output: */
+	IS          isv=NULL;
+	IS          isp=NULL;
+
+	int         start,end;
+	double*     df_local=NULL;
+	int         df_local_size;
+	int         i;
+
+	int*     pressure_indices=NULL;
+	int*     velocity_indices=NULL;
+	int      pressure_num=0;
+	int      velocity_num=0;
+	int      pressure_count=0;
+	int      velocity_count=0;
+
+	if(typeenum==StokesSolverEnum){
+
+		/*Ok, recover doftypes vector values and indices: */
+		VecGetOwnershipRange(df,&start,&end);
+		VecGetLocalSize(df,&df_local_size);
+		VecGetArray(df,&df_local);
+
+		pressure_num=0;
+		velocity_num=0;
+		for(i=0;i<df_local_size;i++){
+			if (df_local[i]==PressureEnum)pressure_num++;
+			else velocity_num++;
+		}
+
+		/*Allocate indices: */
+		if(pressure_num)pressure_indices=(int*)xmalloc(pressure_num*sizeof(int));
+		if(velocity_num)velocity_indices=(int*)xmalloc(velocity_num*sizeof(int));
+
+		pressure_count=0;
+		velocity_count=0;
+		for(i=0;i<df_local_size;i++){
+			if (df_local[i]==PressureEnum){
+				pressure_indices[pressure_count]=start+i;
+				pressure_count++;
+			}
+			if (df_local[i]==VelocityEnum){
+				velocity_indices[velocity_count]=start+i;
+				velocity_count++;
+			}
+		}
+		VecRestoreArray(df,&df_local);
+		
+		/*Create indices sets: */
+		ISCreateGeneral(PETSC_COMM_WORLD,pressure_num,pressure_indices,&isp);
+		ISCreateGeneral(PETSC_COMM_WORLD,velocity_num,velocity_indices,&isv);
+
+	}
+	/*Free ressources:*/
+	xfree((void**)&pressure_indices);
+	xfree((void**)&velocity_indices);
+
+
+	/*Assign output pointers:*/
+	*pisv=isv;
+	*pisp=isp;
+}
Index: /issm/trunk/src/c/modules/Solverx/Solverx.cpp
===================================================================
--- /issm/trunk/src/c/modules/Solverx/Solverx.cpp	(revision 7390)
+++ /issm/trunk/src/c/modules/Solverx/Solverx.cpp	(revision 7391)
@@ -13,5 +13,5 @@
 #endif
 
-void	Solverx(Vec* puf, Mat Kff, Vec pf, Vec uf0,Parameters* parameters){
+void	Solverx(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){
 
 	/*output: */
@@ -29,4 +29,7 @@
 	int        solver_type;
 	bool       fromlocalsize    = true;
+	/*stokes: */
+	IS         isv=NULL;
+	IS         isp=NULL;
 
 	#if _PETSC_VERSION_ == 3 
@@ -35,9 +38,10 @@
 	#endif
 
+
 	/*Display message*/
 	_printf_(VerboseModule(),"   Solving\n");
 	if(VerboseSolver())PetscOptionsPrint(stdout);
 
-	/*First, check that f-set is not NULL, ie model is fully constrained: */
+	/*First, check that f-set is not NULL, ie model is fully constrained: {{{*/
 	_assert_(Kff);
 	MatGetSize(Kff,&global_m,&global_n); _assert_(global_m==global_m);
@@ -45,5 +49,6 @@
 		*puf=NULL; return;
 	}
-
+	/*}}}*/
+	/*Initial guess logic here: {{{1*/
 	/*Now, check that we are not getting an initial guess to the solver, if we are running a direct solver: */
 	#if _PETSC_VERSION_ == 3 
@@ -60,6 +65,6 @@
 		MatGetLocalSize(Kff,&local_m,&local_n);uf=NewVec(local_n,fromlocalsize);
 	}
-
-	/*Process petsc options to see if we are not using special types of external solvers: */
+	/*}}}*/
+	/*Process petsc options to see if we are not using special types of external solvers: {{{1*/
 	PetscOptionsDetermineSolverType(&solver_type);
 
@@ -88,14 +93,16 @@
 		MatConvert(Kff,MATSUPERLU_DIST,MAT_REUSE_MATRIX,&Kff);
 	}
+	if (solver_type==StokesSolverEnum){
+		_error_("Petsc 2 does not support multi-physics solvers");
+	}
 	#endif
 	#endif
-
-	/*Prepare solver*/
+	/*}}}*/
+	/*Prepare solver:{{{1*/
 	KSPCreate(MPI_COMM_WORLD,&ksp);
 	KSPSetOperators(ksp,Kff,Kff,DIFFERENT_NONZERO_PATTERN);
 	KSPSetFromOptions(ksp);
 
-	#ifdef _SERIAL_
-	#if _PETSC_VERSION_ == 3 
+	#if defined(_SERIAL_) && _PETSC_VERSION_==3
 	/*specific solver?: */
 	KSPGetPC(ksp,&pc);
@@ -104,8 +111,22 @@
 	}
 	#endif
+
+	#if defined(_PARALLEL_) && _PETSC_VERSION_==3
+	/*Stokes: */
+	if (solver_type==StokesSolverEnum){
+		/*Make indices out of doftypes: */
+		if(!df)_error_("need doftypes for Stokes solver!\n");
+		DofTypesToIndexSet(&isv,&isp,df,StokesSolverEnum);
+
+		/*Set field splits: */
+		KSPGetPC(ksp,&pc);
+		PCFieldSplitSetIS(pc,isv);
+		PCFieldSplitSetIS(pc,isp);
+
+	}
 	#endif
 
-	/*If initial guess for solution, use it, except if we are using the MUMPS direct solver, where any initial 
-	 * guess will crash Petsc: */
+	/*}}}*/
+	/*If initial guess for solution, use it, except if we are using the MUMPS direct solver, where any initial guess will crash Petsc: {{{1*/
 	if (uf0){
 		if( (solver_type!=MUMPSPACKAGE_LU) && (solver_type!=MUMPSPACKAGE_CHOL) && (solver_type!=SPOOLESPACKAGE_LU)&& (solver_type!=SPOOLESPACKAGE_CHOL)&& (solver_type!=SUPERLUDISTPACKAGE)){
@@ -113,4 +134,5 @@
 		}
 	}
+	/*}}}*/
 	
 	if(VerboseSolver())KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);
@@ -125,5 +147,5 @@
 	/*Free ressources:*/
 	KSPFree(&ksp);
-	
+		
 	/*Assign output pointers:*/
 	*puf=uf;
Index: /issm/trunk/src/c/modules/Solverx/Solverx.h
===================================================================
--- /issm/trunk/src/c/modules/Solverx/Solverx.h	(revision 7390)
+++ /issm/trunk/src/c/modules/Solverx/Solverx.h	(revision 7391)
@@ -9,6 +9,6 @@
 
 /* local prototypes: */
-void	Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0,Parameters* parameters);
-
+void	Solverx( Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df,Parameters* parameters);
+void    DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
 #endif  /* _SOLVERX_H */
 
Index: /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 7390)
+++ /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 7391)
@@ -9,5 +9,5 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){
+void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, Vec* pdg, Vec* pdf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){
 	
 	/*intermediary: */
@@ -26,4 +26,6 @@
 	Vec    pg   = NULL;
 	Vec    pf   = NULL;
+	Vec    dg=NULL;
+	Vec    df=NULL;
 	double kmax = 0;
 
@@ -52,8 +54,12 @@
 	if(kflag){
 
-		if(!buildkff)Kgg=NewMat(gsize,gsize,connectivity,numberofdofspernode);
+		if(!buildkff){
+			Kgg=NewMat(gsize,gsize,connectivity,numberofdofspernode);
+			dg=NewVec(gsize);
+		}
 		else{
 			Kff=NewMat(fsize,fsize,connectivity,numberofdofspernode);
 			Kfs=NewMat(fsize,ssize,connectivity,numberofdofspernode);
+			df=NewVec(fsize);
 		}
 
@@ -61,5 +67,5 @@
 		for (i=0;i<elements->Size();i++){
 			element=(Element*)elements->GetObjectByOffset(i);
-			element->CreateKMatrix(Kgg,Kff,Kfs);
+			element->CreateKMatrix(Kgg,Kff,Kfs,dg,df);
 		}
 
@@ -70,5 +76,5 @@
 		}
 
-		/*Assemble matrix and compress matrix to save memory: */
+		/*Assemble matrix and doftypes and compress matrix to save memory: */
 		if(!buildkff){
 			MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY);
@@ -77,4 +83,6 @@
 			MatCompress(Kgg);
 			#endif
+			VecAssemblyBegin(dg);
+			VecAssemblyEnd(dg);
 		}
 		else{
@@ -90,7 +98,11 @@
 			MatCompress(Kfs);
 			#endif
+			VecAssemblyBegin(df);
+			VecAssemblyEnd(df);
 		}
+		
 	}
-if(pflag){
+	
+	if(pflag){
 
 		if(!buildkff)pg=NewVec(gsize);
@@ -154,5 +166,7 @@
 		}
 	}
-if(penalty_pflag){
+
+	
+	if(penalty_pflag){
 
 		/*Fill right hand side vector, from loads: */
@@ -178,4 +192,6 @@
 	if(ppg)  *ppg=pg;
 	if(ppf)  *ppf=pf;
+	if(pdg)  *pdg=dg;
+	if(pdf)  *pdf=df;
 	if(pkmax) *pkmax=kmax;
 }
Index: /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.h
===================================================================
--- /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.h	(revision 7390)
+++ /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.h	(revision 7391)
@@ -10,5 +10,5 @@
 
 /* local prototypes: */
-void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,
+void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, Vec* pdg, Vec* pdf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,
 			bool kflag=true,bool pflag=true,bool penalty_kflag=true,bool penalty_pflag=true);
 
Index: /issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Element.h	(revision 7390)
+++ /issm/trunk/src/c/objects/Elements/Element.h	(revision 7391)
@@ -28,5 +28,5 @@
 		virtual void   Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0;
 		virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0;
-		virtual void   CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs)=0;
+		virtual void   CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,Vec dg, Vec df)=0;
 		virtual void   CreatePVector(Vec pg, Vec pf)=0;
 		virtual void   GetSolutionFromInputs(Vec solution)=0;
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 7390)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 7391)
@@ -529,8 +529,9 @@
 }/*}}}*/
 /*FUNCTION Penta::CreateKMatrix {{{1*/
-void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
+void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg, Vec df){
 
 	/*retrieve parameters: */
 	ElementMatrix* Ke=NULL;
+	ElementVector* De=NULL;
 	int analysis_type;
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
@@ -546,5 +547,5 @@
 	switch(analysis_type){
 		case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
-			Ke=CreateKMatrixDiagnosticHoriz();
+			Ke=CreateKMatrixDiagnosticHoriz(); De=CreateDVectorDiagnosticHoriz();
 			break;
 		case DiagnosticHutterAnalysisEnum:
@@ -581,4 +582,9 @@
 		delete Ke;
 	}
+	/*Add to global Vector*/
+	if(De){
+		De->InsertIntoGlobal(dg,df);
+		delete De;
+	}
 }
 /*}}}*/
@@ -1006,4 +1012,18 @@
 }
 /*}}}*/
+/*FUNCTION Penta::CreateDVectorDiagnosticHoriz {{{1*/
+ElementVector* Penta::CreateDVectorDiagnosticHoriz(void){
+
+	int approximation;
+	inputs->GetParameterValue(&approximation,ApproximationEnum);
+
+	switch(approximation){
+		case StokesApproximationEnum:
+			return CreateDVectorDiagnosticStokes();
+		default:
+			return NULL; //no need for doftypes outside of stokes approximation
+	}
+}
+/*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/
 ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){
@@ -1365,4 +1385,29 @@
 	delete Ke2;
 	return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateDVectorDiagnosticStokes{{{1*/
+ElementVector* Penta::CreateDVectorDiagnosticStokes(void){
+
+	/*output: */
+	ElementVector* De=NULL;
+	/*intermediary: */
+	int approximation;
+	int i;
+
+	/*Initialize Element vector and return if necessary*/
+	inputs->GetParameterValue(&approximation,ApproximationEnum);
+	if(approximation!=StokesApproximationEnum) return NULL;
+
+	De=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+
+	for (i=0;i<NUMVERTICES;i++){
+		De->values[i*4+0]=VelocityEnum;
+		De->values[i*4+1]=VelocityEnum;
+		De->values[i*4+2]=VelocityEnum;
+		De->values[i*4+3]=PressureEnum;
+	}
+
+	return De;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 7390)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 7391)
@@ -80,5 +80,5 @@
 		void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
 		double RegularizeInversion(void);
-		void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
+		void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df);
 		void   CreatePVector(Vec pg, Vec pf);
 		void   DeleteResults(void);
@@ -147,4 +147,6 @@
 		ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
 		ElementMatrix* CreateKMatrixDiagnosticHoriz(void);
+		ElementVector* CreateDVectorDiagnosticHoriz(void);
+		ElementVector* CreateDVectorDiagnosticStokes(void);
 		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 7390)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 7391)
@@ -377,5 +377,5 @@
 /*}}}*/
 /*FUNCTION Tria::CreateKMatrix {{{1*/
-void  Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
+void  Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df){
 
 	/*retreive parameters: */
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 7390)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 7391)
@@ -77,5 +77,5 @@
 		void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
 		double RegularizeInversion(void);
-		void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs);
+		void   CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs,Vec dg,Vec df);
 		void   CreatePVector(Vec pg, Vec pf);
 		int    GetNodeIndex(Node* node);
Index: /issm/trunk/src/c/objects/Numerics/ElementVector.cpp
===================================================================
--- /issm/trunk/src/c/objects/Numerics/ElementVector.cpp	(revision 7390)
+++ /issm/trunk/src/c/objects/Numerics/ElementVector.cpp	(revision 7391)
@@ -197,4 +197,30 @@
 }
 /*}}}*/
+/*FUNCTION ElementVector::InsertIntoGlobal(Vec pg, Vec pf){{{1*/
+void ElementVector::InsertIntoGlobal(Vec pg, Vec pf){
+
+	int i;
+	double* localvalues=NULL;
+
+	if(!pf){
+		/*add local values into global  vector, using the fglobaldoflist: */
+		VecSetValues(pg,this->nrows,this->gglobaldoflist,(const double*)values,INSERT_VALUES);
+	}
+	else{
+		if(this->fsize){
+			/*first, retrieve values that are in the f-set from the g-set values vector: */
+			localvalues=(double*)xmalloc(this->fsize*sizeof(double));
+			for(i=0;i<this->fsize;i++){
+				localvalues[i]=this->values[this->flocaldoflist[i]];
+			}
+			/*add local values into global  vector, using the fglobaldoflist: */
+			VecSetValues(pf,this->fsize,this->fglobaldoflist,(const double*)localvalues,INSERT_VALUES);
+
+			/*Free ressources:*/
+			xfree((void**)&localvalues);
+		}
+	}
+}
+/*}}}*/
 /*FUNCTION ElementVector::Echo{{{1*/
 void ElementVector::Echo(void){
Index: /issm/trunk/src/c/objects/Numerics/ElementVector.h
===================================================================
--- /issm/trunk/src/c/objects/Numerics/ElementVector.h	(revision 7390)
+++ /issm/trunk/src/c/objects/Numerics/ElementVector.h	(revision 7391)
@@ -42,4 +42,5 @@
 		/*ElementVector specific routines {{{1*/
 		void AddToGlobal(Vec pg, Vec pf);
+		void InsertIntoGlobal(Vec pg, Vec pf);
 		void Echo(void);
 		void Init(ElementVector* pe);
Index: /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp	(revision 7390)
+++ /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp	(revision 7391)
@@ -16,4 +16,5 @@
 	Vec  ug  = NULL, uf  = NULL;
 	Vec  pg  = NULL, pf  = NULL;
+	Vec  dg  = NULL, df  = NULL;
 	bool kffpartitioning = false;
 
@@ -22,14 +23,15 @@
 
 	if(kffpartitioning){
-		SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 		Reduceloadx(pf, Kfs, femmodel->ys,true); MatFree(&Kfs); //true means spc = 0
 	}
 	else{
-		SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, &dg, NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets); MatFree(&Kgg);
 		Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,true);VecFree(&pg); MatFree(&Kfs);//true means spc = 0
+		Reducevectorgtofx(&df, dg, femmodel->nodesets,femmodel->parameters); VecFree(&dg);
 	}
 
-	Solverx(&uf, Kff, pf, NULL, femmodel->parameters); MatFree(&Kff); VecFree(&pf);
+	Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); MatFree(&Kff); VecFree(&pf); VecFree(&df);
 	Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,true); VecFree(&uf);//true means spc0
 	InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); VecFree(&ug); VecFree(&uf);
Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 7390)
+++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 7391)
@@ -16,4 +16,6 @@
 	Vec ug  = NULL, uf  = NULL, old_ug= NULL, old_uf = NULL;
 	Vec pg  = NULL, pf  = NULL;
+	Vec dg  = NULL, df  = NULL;
+	
 	Loads* loads=NULL;
 	int converged;
@@ -53,14 +55,15 @@
 
 		if(kffpartitioning){
-			SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
+			SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
 			Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
 		}
 		else{
-			SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
+			SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, &dg, NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
 			Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets); MatFree(&Kgg);
 			Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets); VecFree(&pg); MatFree(&Kfs);
+			Reducevectorgtofx(&df, dg, femmodel->nodesets,femmodel->parameters); VecFree(&dg);
 		}
 		
-		Solverx(&uf, Kff, pf, old_uf, femmodel->parameters);
+		Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
 		Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets);
 		InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
@@ -69,5 +72,5 @@
 		_printf_(VerboseConvergence(),"   number of unstable constraints: %i\n",num_unstable_constraints);
 
-		convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf);
+		convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); VecFree(&df);
 		InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
 
Index: /issm/trunk/src/c/solvers/solver_linear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_linear.cpp	(revision 7390)
+++ /issm/trunk/src/c/solvers/solver_linear.cpp	(revision 7391)
@@ -14,4 +14,5 @@
 	Vec  ug  = NULL, uf  = NULL;
 	Vec  pg  = NULL, pf  = NULL;
+	Vec dg  = NULL, df  = NULL;
 	bool kffpartitioning;
 
@@ -20,14 +21,15 @@
 
 	if(kffpartitioning){
-		SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 		Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
 	}
 	else{
-		SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, &dg, NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets); MatFree(&Kgg);
 		Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets);VecFree(&pg); MatFree(&Kfs);
+		Reducevectorgtofx(&df, dg, femmodel->nodesets,femmodel->parameters); VecFree(&dg);
 	}
 
-	Solverx(&uf, Kff, pf, NULL, femmodel->parameters); MatFree(&Kff); VecFree(&pf);
+	Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters); MatFree(&Kff); VecFree(&pf); VecFree(&df);
 	Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets);VecFree(&uf);
 	InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); VecFree(&ug); VecFree(&uf);
Index: /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp	(revision 7390)
+++ /issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp	(revision 7391)
@@ -16,7 +16,9 @@
 	Vec ug_horiz  = NULL, uf_horiz  = NULL, old_ug_horiz= NULL, old_uf_horiz = NULL;
 	Vec pg_horiz  = NULL, pf_horiz  = NULL;
+	Vec dg_horiz  = NULL, df_horiz  = NULL;
 	Mat Kgg_vert  = NULL, Kff_vert  = NULL, Kfs_vert    = NULL;
 	Vec ug_vert   = NULL, uf_vert   = NULL;
 	Vec pg_vert   = NULL, pf_vert   = NULL;
+	Vec dg_vert   = NULL, df_vert   = NULL;
 	int converged;
 	int constraints_converged;
@@ -56,32 +58,34 @@
 
 		if(kffpartitioning){
-			SystemMatricesx(NULL,&Kff_horiz, &Kfs_horiz, NULL,&pf_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			SystemMatricesx(NULL,&Kff_horiz, &Kfs_horiz, NULL,&pf_horiz, NULL, &df_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 			Reduceloadx(pf_horiz, Kfs_horiz, femmodel->ys); MatFree(&Kfs_horiz);
 		}
 		else{
-			SystemMatricesx(&Kgg_horiz, NULL, NULL, &pg_horiz,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			SystemMatricesx(&Kgg_horiz, NULL, NULL, &pg_horiz,NULL, &dg_horiz,NULL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 			Reducematrixfromgtofx(&Kff_horiz,&Kfs_horiz,Kgg_horiz,femmodel->nodesets); MatFree(&Kgg_horiz);
 			Reduceloadfromgtofx(&pf_horiz, pg_horiz, Kfs_horiz, femmodel->ys, femmodel->nodesets); VecFree(&pg_horiz); MatFree(&Kfs_horiz);
+			Reducevectorgtofx(&df_horiz, dg_horiz, femmodel->nodesets,femmodel->parameters); VecFree(&dg_horiz);
 		}
 		
-		Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, femmodel->parameters);
+		Solverx(&uf_horiz, Kff_horiz, pf_horiz, old_uf_horiz, df_horiz,femmodel->parameters);
 		Mergesolutionfromftogx(&ug_horiz, uf_horiz,femmodel->ys,femmodel->nodesets);
 		InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_horiz);
 
-		convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); MatFree(&Kff_horiz);VecFree(&pf_horiz);
+		convergence(&converged,Kff_horiz,pf_horiz,uf_horiz,old_uf_horiz,femmodel->parameters); MatFree(&Kff_horiz);VecFree(&pf_horiz); VecFree(&df_horiz);
 
 		/*Second compute vertical velocity: */
 		femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
 		if(kffpartitioning){
-			SystemMatricesx(NULL,&Kff_vert, &Kfs_vert, NULL,&pf_vert, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			SystemMatricesx(NULL,&Kff_vert, &Kfs_vert, NULL,&pf_vert, NULL, &df_vert,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 			Reduceloadx(pf_vert, Kfs_vert, femmodel->ys); MatFree(&Kfs_vert);
 		}
 		else{
-			SystemMatricesx(&Kgg_vert, NULL, NULL, &pg_vert,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			SystemMatricesx(&Kgg_vert, NULL, NULL, &pg_vert,NULL, &dg_horiz, NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 			Reducematrixfromgtofx(&Kff_vert,&Kfs_vert,Kgg_vert,femmodel->nodesets); MatFree(&Kgg_vert);
 			Reduceloadfromgtofx(&pf_vert, pg_vert, Kfs_vert, femmodel->ys, femmodel->nodesets);VecFree(&pg_vert); MatFree(&Kfs_vert);
+			Reducevectorgtofx(&df_vert, dg_vert, femmodel->nodesets,femmodel->parameters); VecFree(&dg_vert);
 		}
 
-		Solverx(&uf_vert, Kff_vert, pf_vert, NULL, femmodel->parameters); MatFree(&Kff_vert); VecFree(&pf_vert);
+		Solverx(&uf_vert, Kff_vert, pf_vert, NULL, df_vert,femmodel->parameters); MatFree(&Kff_vert); VecFree(&pf_vert); VecFree(&df_vert);
 		Mergesolutionfromftogx(&ug_vert, uf_vert,femmodel->ys,femmodel->nodesets);VecFree(&uf_vert);
 		InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug_vert); VecFree(&ug_vert); VecFree(&uf_vert);
Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 7390)
+++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 7391)
@@ -22,4 +22,6 @@
 	Vec pg=NULL;
 	Vec pf=NULL;
+	Vec dg=NULL;
+	Vec df=NULL;
 
 	bool converged;
@@ -50,17 +52,18 @@
 
 		if(kffpartitioning){
-			SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf,NULL, &df, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 			Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
 		}
 		else{
-			SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, &dg, NULL, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 			Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets); MatFree(&Kgg);
 			Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets);VecFree(&pg); MatFree(&Kfs);
+			Reducevectorgtofx(&df, dg, femmodel->nodesets,femmodel->parameters); VecFree(&dg);
 		}
 
 		VecFree(&tf);
-		Solverx(&tf, Kff, pf,tf_old, femmodel->parameters);
+		Solverx(&tf, Kff, pf,tf_old, df, femmodel->parameters);
 		VecFree(&tf_old); VecDuplicatePatch(&tf_old,tf);
-		MatFree(&Kff);VecFree(&pf);VecFree(&tg);
+		MatFree(&Kff);VecFree(&pf);VecFree(&tg); VecFree(&df);
 
 		Mergesolutionfromftogx(&tg, tf,femmodel->ys,femmodel->nodesets);
Index: /issm/trunk/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp	(revision 7390)
+++ /issm/trunk/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp	(revision 7391)
@@ -45,4 +45,9 @@
 		solver_type=SUPERLUDISTPACKAGE;
 	}
+	
+	PetscOptionsGetString(PETSC_NULL,"-issm_option_solver",&option[0],100,&flag);
+	if (strcmp(option,"stokes")==0){
+		solver_type=StokesSolverEnum;
+	}
 
 	*psolver_type=solver_type;
Index: /issm/trunk/src/m/model/solvers/stokesoptions.m
===================================================================
--- /issm/trunk/src/m/model/solvers/stokesoptions.m	(revision 7391)
+++ /issm/trunk/src/m/model/solvers/stokesoptions.m	(revision 7391)
@@ -0,0 +1,39 @@
+function options=stokesoptions(varargin)
+%STOKESOPTIONS - return STOKES multi-physics solver petsc options
+%
+%   Usage:
+%      options=stokesoptions;
+
+%retrieve options provided in varargin
+arguments=pairoptions(varargin{:});
+
+%default stokes options
+PETSC_VERSION=petscversion();
+if PETSC_VERSION==2,
+	error('stokesoptions error message: multi-physics options not supported in Petsc 2');
+end
+if PETSC_VERSION==3,
+	options={{'mat_type','mpiaij'},{'ksp_type','gmres'},{'pc_type','fieldsplit'},{'pc_field_split_type','multiplicative'},...
+	{'fieldsplit_0_pc_type','asm'},{'fieldsplit_0_ksp_type','gmres'},{'fieldsplit_0_sub_pc_type','jacobi'},{'fieldsplit_0_pc_asm_overlap',2},...
+	{'fieldsplit_1_pc_type','jacobi'},{'fieldsplit_1_ksp_type','preonly'},{'issm_option_solver','stokes'}};
+end
+
+%now, go through our arguments, and write over default options.
+for i=1:size(arguments.list,1),
+	arg1=arguments.list{i,1};
+	arg2=arguments.list{i,2};
+	found=0;
+	for j=1:size(options,2),
+		joption=options{j};
+		if strcmpi(joption{1},arg1),
+			joption{2}=arg2;
+			options{j}=joption;
+			found=1;
+			break;
+		end
+	end
+	if ~found,
+		%this option did not exist, add it: 
+		options{end+1}={arg1,arg2};
+	end
+end
Index: /issm/trunk/src/m/solvers/solver_adjoint_linear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_adjoint_linear.m	(revision 7390)
+++ /issm/trunk/src/m/solvers/solver_adjoint_linear.m	(revision 7391)
@@ -9,14 +9,16 @@
 
 	if kffpartitioning,
-		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+		[K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 		p_f = Reduceload( p_f, K_fs, femmodel.ys,true);
+
 	else
-		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+		[K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets); 
 		p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,true);
+		df=Reducevectorgtof( dg, femmodel.nodesets,femmodel.parameters);
 	end
 
 	issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
-	u_f=Solver(K_ff,p_f,[],femmodel.parameters);
+	u_f=Solver(K_ff,p_f,[],df,femmodel.parameters);
 	u_g= Mergesolutionfromftog( u_f, femmodel.ys, femmodel.nodesets,true); 
 
Index: /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m	(revision 7390)
+++ /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m	(revision 7391)
@@ -29,13 +29,14 @@
 
 		if kffpartitioning,
-			[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
+			[K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
 			p_f = Reduceload( p_f, K_fs, femmodel.ys);
 		else
-			[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
+			[K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
 			[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets); 
 			p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets);
+			df=Reducevectorgtof( dg, femmodel.nodesets,femmodel.parameters);
 		end
 
-		uf=Solver(K_ff,p_f,old_uf,femmodel.parameters);
+		uf=Solver(K_ff,p_f,old_uf,df,femmodel.parameters);
 		ug= Mergesolutionfromftog( uf, femmodel.ys, femmodel.nodesets); 
 
Index: /issm/trunk/src/m/solvers/solver_linear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_linear.m	(revision 7390)
+++ /issm/trunk/src/m/solvers/solver_linear.m	(revision 7391)
@@ -9,14 +9,15 @@
 
 	if kffpartitioning,
-		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+		[K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 		p_f = Reduceload( p_f, K_fs, femmodel.ys);
 	else
-		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+		[K_gg,K_ff,K_fs,p_g,p_f,dg,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets); 
 		p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets);
+		df=Reducevectorgtof( dg, femmodel.nodesets,femmodel.parameters);
 	end
 
 	issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
-	u_f=Solver(K_ff,p_f,[],femmodel.parameters);
+	u_f=Solver(K_ff,p_f,[],df,femmodel.parameters);
 	u_g= Mergesolutionfromftog( u_f, femmodel.ys, femmodel.nodesets); 
 
Index: /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m	(revision 7390)
+++ /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m	(revision 7391)
@@ -1,3 +1,3 @@
-function femmodel=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
+fpunction femmodel=solver_couplingstokes_nonlinear(femmodel,conserve_loads)
 %SOLVER_COUPLINGSTOKES_NONLINEAR - core solver of coupling run
 %
@@ -29,13 +29,15 @@
 
 		if kffpartitioning,
-			[K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax_horiz]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			[K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,d_g_horiz,d_f_horiz,kmax_horiz]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 			p_f_horiz = Reduceload( p_f_horiz, K_fs_horiz, femmodel.ys);
 		else
-			[K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			[K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,d_g_horiz,d_f_horiz,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 			[K_ff_horiz, K_fs_horiz] = Reducematrixfromgtof( K_gg_horiz, femmodel.nodesets); 
 			p_f_horiz = Reduceloadfromgtof( p_g_horiz, K_fs_horiz, femmodel.ys, femmodel.nodesets);
+			d_f_horiz=Reducevectorgtof( d_g_horiz, femmodel.nodesets,femmodel.parameters);
+
 		end
 
-		uf_horiz=Solver(K_ff_horiz,p_f_horiz,old_uf_horiz,femmodel.parameters);
+		uf_horiz=Solver(K_ff_horiz,p_f_horiz,old_uf_horiz,d_f_horiz,femmodel.parameters);
 		ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets); 
 
@@ -49,13 +51,14 @@
 
 		if kffpartitioning,
-			[K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			[K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,d_g_vert,d_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 			p_f_vert = Reduceload( p_f_vert, K_fs_vert, femmodel.ys);
 		else
-			[K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			[K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,d_g_vert,d_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 			[K_ff_vert, K_fs_vert] = Reducematrixfromgtof( K_gg_vert, femmodel.nodesets); 
 			p_f_vert = Reduceloadfromgtof( p_g_vert, K_fs_vert, femmodel.ys, femmodel.nodesets);
+			d_f_vert=Reducevectorgtof( d_g_vert, femmodel.nodesets,femmodel.parameters);
 		end
 
-		uf_vert=Solver(K_ff_vert,p_f_vert,[],femmodel.parameters);
+		uf_vert=Solver(K_ff_vert,p_f_vert,[],d_f_vert,femmodel.parameters);
 		ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets); 
 
Index: /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 7390)
+++ /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 7391)
@@ -24,14 +24,16 @@
 
 		if kffpartitioning,
-			[K_gg,K_ff,K_fs,p_g,p_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			[K_gg,K_ff,K_fs,p_g,p_f,d_g,d_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 			p_f = Reduceload( p_f, K_fs, femmodel.ys);
 		else
-			[K_gg,K_ff,K_fs,p_g,p_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			[K_gg,K_ff,K_fs,p_g,p_f,d_g,d_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 			[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets); 
 			p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets);
+			d_f=Reducevectorgtof( d_g, femmodel.nodesets,femmodel.parameters);
+
 		end
 
 		issmprintf(VerboseSolver(),'%s%g','   condition number of stiffness matrix: ',condest(K_ff));
-		t_f=Solver(K_ff,p_f,[],femmodel.parameters);
+		t_f=Solver(K_ff,p_f,[],d_f,femmodel.parameters);
 		t_g= Mergesolutionfromftog( t_f, femmodel.ys, femmodel.nodesets); 
 
Index: /issm/trunk/src/mex/Solver/Solver.cpp
===================================================================
--- /issm/trunk/src/mex/Solver/Solver.cpp	(revision 7390)
+++ /issm/trunk/src/mex/Solver/Solver.cpp	(revision 7391)
@@ -12,4 +12,5 @@
 	Vec         uf0           = NULL;
 	Vec         uf            = NULL;
+	Vec         df            = NULL;
 	Parameters *parameters    = NULL;
 	int         verbose;
@@ -47,7 +48,8 @@
 		FetchData(&pf,PF);
 		FetchData(&uf0,UF0);
+		FetchData(&df,DF);
 
 		/*Core module: */
-		Solverx(&uf, Kff, pf, uf0, parameters);
+		Solverx(&uf, Kff, pf, uf0, df,parameters);
 
 		/*Write output*/
@@ -67,4 +69,5 @@
 	VecFree(&uf0);
 	VecFree(&uf);
+	VecFree(&df);
 	delete parameters;
 
@@ -76,5 +79,5 @@
 {
 	_printf_(true,"\n");
-	_printf_(true,"   usage: [uf] = %s(Kff,pf,uf0,parameters);\n",__FUNCT__);
+	_printf_(true,"   usage: [uf] = %s(Kff,pf,uf0,df,parameters);\n",__FUNCT__);
 	_printf_(true,"\n");
 }
Index: /issm/trunk/src/mex/Solver/Solver.h
===================================================================
--- /issm/trunk/src/mex/Solver/Solver.h	(revision 7390)
+++ /issm/trunk/src/mex/Solver/Solver.h	(revision 7391)
@@ -23,5 +23,6 @@
 #define PF (mxArray*)prhs[1]
 #define UF0 (mxArray*)prhs[2]
-#define PARAMETERS (mxArray*)prhs[3]
+#define DF (mxArray*)prhs[3]
+#define PARAMETERS (mxArray*)prhs[4]
 
 /* serial output macros: */
@@ -32,5 +33,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  4
+#define NRHS  5
 
 
Index: /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp
===================================================================
--- /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp	(revision 7390)
+++ /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp	(revision 7391)
@@ -22,4 +22,7 @@
 	Vec    pg   = NULL;
 	Vec    pf   = NULL;
+	Vec    dg   = NULL;
+	Vec    df   = NULL;
+
 	double kmax;
 
@@ -53,8 +56,8 @@
 		FetchData(&penalty_kflag,PENALTYKFLAG);
 		FetchData(&penalty_pflag,PENALTYPFLAG);
-		SystemMatricesx(&Kgg,&Kff,&Kfs,&pg,&pf,&kmax,elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);
+		SystemMatricesx(&Kgg,&Kff,&Kfs,&pg,&pf,&dg,&df,&kmax,elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);
 	}
 	else
-	 SystemMatricesx(&Kgg,&Kff,&Kfs,&pg,&pf,&kmax,elements,nodes,vertices,loads,materials,parameters);
+	 SystemMatricesx(&Kgg,&Kff,&Kfs,&pg,&pf,&dg,&df,&kmax,elements,nodes,vertices,loads,materials,parameters);
 
 	/*write output datasets: */
@@ -64,4 +67,6 @@
 	WriteData(PG,pg);
 	WriteData(PF,pf);
+	WriteData(DG,dg);
+	WriteData(DF,df);
 	WriteData(KMAX,kmax);
 	
@@ -78,4 +83,6 @@
 	VecFree(&pg);
 	VecFree(&pf);
+	VecFree(&dg);
+	VecFree(&df);
 
 	/*end module: */
@@ -86,6 +93,6 @@
 {
 	_printf_(true,"\n");
-	_printf_(true,"   usage: [Kgg,Kff,Kfs,pg,pf,kmax] = %s(elements,nodes,vertices,loads,materials,parameters);\n",__FUNCT__);
-	_printf_(true,"   usage: [Kgg,Kff,Kfs,pg,pf,kmax] = %s(elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);\n",__FUNCT__);
+	_printf_(true,"   usage: [Kgg,Kff,Kfs,pg,pf,dg,df,kmax] = %s(elements,nodes,vertices,loads,materials,parameters);\n",__FUNCT__);
+	_printf_(true,"   usage: [Kgg,Kff,Kfs,pg,pf,dg,df,kmax] = %s(elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);\n",__FUNCT__);
 	_printf_(true,"\n");
 }
Index: /issm/trunk/src/mex/SystemMatrices/SystemMatrices.h
===================================================================
--- /issm/trunk/src/mex/SystemMatrices/SystemMatrices.h	(revision 7390)
+++ /issm/trunk/src/mex/SystemMatrices/SystemMatrices.h	(revision 7391)
@@ -35,9 +35,11 @@
 #define PG   (mxArray**)&plhs[3]
 #define PF   (mxArray**)&plhs[4]
-#define KMAX (mxArray**)&plhs[5]
+#define DG   (mxArray**)&plhs[5]
+#define DF   (mxArray**)&plhs[6]
+#define KMAX (mxArray**)&plhs[7]
 
 /* serial arg counts: */
 #undef NLHS
-#define NLHS  6
+#define NLHS  8
 #undef NRHS
 #define NRHS  10