Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 23198)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 23199)
@@ -74,6 +74,7 @@
 		ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
 		ElementMatrix* CreatePressureMassMatrix(Element* element);
+		ElementMatrix* CreateSchurPrecondMatrix(Element* element);
 		ElementVector* CreatePVectorFS(Element* element);
 		ElementVector* CreatePVectorFSFriction(Element* element);
 		ElementVector* CreatePVectorFSFront(Element* element);
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 23198)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 23199)
@@ -944,9 +944,11 @@
 		#ifdef _HAVE_PETSC_
 		int solver_type;
 		PetscOptionsDetermineSolverType(&solver_type);
+
 		if(solver_type==FSSolverEnum) is_schur_cg_solver = true;
 		#endif
 
+		
 		if(is_schur_cg_solver)
 		 solutionsequence_schurcg(femmodel);
 		else if (fe_FS==XTaylorHoodEnum)
@@ -3288,6 +3290,64 @@
 	xDelete<IssmDouble>(pbasis);
 	return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateSchurPrecondMatrix(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int         i,dim;
+	IssmDouble  viscosity,FSreconditioning,Jdet;
+	IssmDouble *xyz_list = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&dim,DomainDimensionEnum);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumnodes = element->NumberofNodesPressure();
+	int numdof    = vnumnodes*dim + pnumnodes;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke  = element->NewElementMatrix(FSvelocityEnum);
+	IssmDouble* pbasis = xNew<IssmDouble>(pnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	//element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
+	Input* vz_input = NULL;
+	if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
+
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss = element->NewGauss(5);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctionsPressure(pbasis,gauss);
+		element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+
+
+		if(dim==3){
+			/*Pressure mass matrix*/
+			for(int k=0;k<pnumnodes;k++){
+				for(int j=0;j<pnumnodes;j++){
+					Ke->values[(3*vnumnodes+k)*numdof+3*vnumnodes+j] += gauss->weight*1./viscosity*Jdet*(pbasis[j]*pbasis[k]);
+				}
+			}
+		}
+		else{
+			_error_("STOP");
+		}
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(pbasis);
+	return Ke;
+}/*}}}*/
+
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLA(Element* element){/*{{{*/
 
 	/*Intermediaries*/