Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 23161)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 23162)
@@ -73,6 +73,7 @@
 		ElementMatrix* CreateKMatrixFSViscous(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
+		ElementMatrix* CreatePressureMassMatrix(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 23161)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 23162)
@@ -3242,6 +3242,56 @@
 	xDelete<int>(cs_list);
 	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){/*{{{*/
 
 	/*Intermediaries*/
Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp	(revision 23161)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_schurcg.cpp	(revision 23162)
@@ -10,7 +10,7 @@
 #include "../analyses/analyses.h"
 
 #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*/
 	Vec uf = NULL;
@@ -34,6 +34,7 @@
 	Vector<IssmDouble>* pf0 = NULL;
 	Vector<IssmDouble>* df  = NULL;
 	Vector<IssmDouble>* ys  = NULL;
+	Matrix<IssmDouble>* Mff = NULL;
 
 	/*parameters:*/
 	int max_nonlinear_iterations;
@@ -71,17 +72,32 @@
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		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);
 		_assert_(Kff->type==PetscMatType); 
 		SchurCGSolver(&uf,
 					Kff->pmatrix->matrix,
+					Mff->pmatrix->matrix,
 					pf->pvector->vector,
 					old_uf->pvector->vector,
 					df->pvector->vector,
 					femmodel->parameters);
 		femmodel->profiler->Stop(SOLVER);
-		delete pf0;
+		delete pf0; delete Mff;
 
 		/*Merge solution from f set to g set*/
 		Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;