Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 23160)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 23161)
@@ -251,4 +251,5 @@
 					./solutionsequences/solutionsequence_newton.cpp\
 					./solutionsequences/solutionsequence_fct.cpp\
+					./solutionsequences/solutionsequence_schurcg.cpp\
 					./solutionsequences/convergence.cpp\
 					./classes/Options/Options.cpp\
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 23160)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 23161)
@@ -944,5 +944,15 @@
 
 		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
-		if (fe_FS==XTaylorHoodEnum)
+
+		bool is_schur_cg_solver = false;
+		#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)
 		 solutionsequence_la_theta(femmodel);
 		else if (fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum)
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 23161)
@@ -0,0 +1,116 @@
+/*!\file: solutionsequence_schurcg.cpp
+ * \brief: numerical core of 
+ */ 
+
+#include "./solutionsequences.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../analyses/analyses.h"
+
+#ifdef _HAVE_PETSC_
+void SchurCGSolver(Vector<IssmDouble>** puf,Mat Kff, Vec pf, Vec uf0,Vec df,Parameters* parameters){/*{{{*/
+
+	/*Initialize output*/
+	Vec uf = NULL;
+
+	_error_("not implemented yet");
+
+	/*return output pointer*/
+	Vector<IssmDouble>* out_uf=new Vector<IssmDouble>(uf);
+	VecFree(&uf);
+	*puf=out_uf;
+}/*}}}*/
+void solutionsequence_schurcg(FemModel* femmodel){/*{{{*/
+
+	/*intermediary: */
+	Matrix<IssmDouble>* Kff = NULL;
+	Matrix<IssmDouble>* Kfs = NULL;
+	Vector<IssmDouble>* ug  = NULL;
+	Vector<IssmDouble>* uf  = NULL;
+	Vector<IssmDouble>* old_uf = NULL;
+	Vector<IssmDouble>* pf  = NULL;
+	Vector<IssmDouble>* pf0 = NULL;
+	Vector<IssmDouble>* df  = NULL;
+	Vector<IssmDouble>* ys  = NULL;
+
+	/*parameters:*/
+	int max_nonlinear_iterations;
+	int configuration_type;
+	IssmDouble eps_res,eps_rel,eps_abs;
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&max_nonlinear_iterations,StressbalanceMaxiterEnum);
+	femmodel->parameters->FindParam(&eps_res,StressbalanceRestolEnum);
+	femmodel->parameters->FindParam(&eps_rel,StressbalanceReltolEnum);
+	femmodel->parameters->FindParam(&eps_abs,StressbalanceAbstolEnum);
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->UpdateConstraintsx();
+
+	int  count=0;
+	bool converged=false;
+
+	/*Start non-linear iteration using input velocity: */
+	GetSolutionFromInputsx(&ug,femmodel);
+	Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
+
+	/*Update once again the solution to make sure that vx and vxold are similar*/
+	InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+	InputUpdateFromSolutionx(femmodel,ug);
+
+	for(;;){
+
+		/*save pointer to old velocity*/
+		delete old_uf; old_uf=uf;
+		delete ug;
+
+		/*Get stiffness matrix and Load vector*/
+		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
+		pf0=pf->Duplicate(); pf->Copy(pf0);
+		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
+		Reduceloadx(pf, Kfs, ys); delete Kfs;
+
+		/*Solve*/
+		femmodel->profiler->Start(SOLVER);
+		_assert_(Kff->type==PetscMatType); 
+		SchurCGSolver(&uf,
+					Kff->pmatrix->matrix,
+					pf->pvector->vector,
+					old_uf->pvector->vector,
+					df->pvector->vector,
+					femmodel->parameters);
+		femmodel->profiler->Stop(SOLVER);
+		delete pf0;
+
+		/*Merge solution from f set to g set*/
+		Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
+
+		/*Check for convergence and update inputs accordingly*/
+		convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df;
+		count++;
+		if(count>=max_nonlinear_iterations){
+			_printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
+			converged=true;
+		}
+		InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
+		InputUpdateFromSolutionx(femmodel,ug);
+
+		/*Increase count: */
+		if(converged==true){
+			femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,StressbalanceConvergenceNumStepsEnum,count));
+			break;
+		}
+	}
+
+	if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count << "\n");
+
+	/*clean-up*/
+	delete uf;
+	delete ug;
+	delete old_uf;
+
+}/*}}}*/
+#else
+void solutionsequence_schurcg(FemModel* femmodel){_error_("PETSc needs to be installed");}
+#endif
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h	(revision 23160)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h	(revision 23161)
@@ -23,4 +23,5 @@
 void solutionsequence_la_theta(FemModel* femmodel);
 void solutionsequence_adjoint_linear(FemModel* femmodel);
+void solutionsequence_schurcg(FemModel* femmodel);
 
 /*convergence*/
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp	(revision 23160)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp	(revision 23161)
@@ -72,5 +72,5 @@
 	PetscOptionsGetString(PETSC_NULL,"-issm_option_solver",&option[0],100,&flag);
 	#endif
-	if (strcmp(option,"FS")==0){
+	if(strcmp(option,"FS")==0 || strcmp(option,"stokes")==0){
 		solver_type=FSSolverEnum;
 	}