Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 23161)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 23162)
@@ -3243,4 +3243,54 @@
 	return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreatePressureMassMatrix(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int         i,dim;
+	IssmDouble  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);
+
+	/* 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);
+
+		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*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){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 23161)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 23162)
@@ -74,4 +74,5 @@
 		ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
+		ElementMatrix* CreatePressureMassMatrix(Element* element);
 		ElementVector* CreatePVectorFS(Element* element);
 		ElementVector* CreatePVectorFSFriction(Element* element);
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp	(revision 23161)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp	(revision 23162)
@@ -11,5 +11,5 @@
 
 #ifdef _HAVE_PETSC_
-void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff, Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/
+void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff,Mat Mff,Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/
 
 	/*Initialize output*/
@@ -35,4 +35,5 @@
 	Vector<IssmDouble>* df  = NULL;
 	Vector<IssmDouble>* ys  = NULL;
+	Matrix<IssmDouble>* Mff = NULL;
 
 	/*parameters:*/
@@ -72,4 +73,18 @@
 		Reduceloadx(pf, Kfs, ys); delete Kfs;
 
+		/*Create mass matrix*/
+		int fsize; Kff->GetSize(&fsize,&fsize);
+		Mff=new Matrix<IssmDouble>(fsize,fsize,100,4);
+		StressbalanceAnalysis* analysis = new StressbalanceAnalysis();
+		/*Get complete stiffness matrix without penalties*/
+		for(int i=0;i<femmodel->elements->Size();i++){
+			Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+			ElementMatrix* Me = analysis->CreatePressureMassMatrix(element);
+			if(Me) Me->AddToGlobal(Mff,NULL);
+			delete Me;
+		}
+		Mff->Assemble();
+		delete analysis;
+
 		/*Solve*/
 		femmodel->profiler->Start(SOLVER);
@@ -77,4 +92,5 @@
 		SchurCGSolver(&uf,
 					Kff->pmatrix->matrix,
+					Mff->pmatrix->matrix,
 					pf->pvector->vector,
 					old_uf->pvector->vector,
@@ -82,5 +98,5 @@
 					femmodel->parameters);
 		femmodel->profiler->Stop(SOLVER);
-		delete pf0;
+		delete pf0; delete Mff;
 
 		/*Merge solution from f set to g set*/