Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 18343)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 18344)
@@ -354,4 +354,5 @@
 					./solutionsequences/solutionsequence_nonlinear.cpp\
 					./solutionsequences/solutionsequence_newton.cpp\
+					./solutionsequences/solutionsequence_fct.cpp\
 					./solutionsequences/convergence.cpp\
 					./classes/Options/Options.h\
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 18343)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 18344)
@@ -70,4 +70,7 @@
 	if(stabilization==3){
 		iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for DG, we need the spc in the element
+	}
+	if(stabilization==4){
+		iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for FCT, we need the spc in the element (penlaties)
 	}
 
@@ -148,8 +151,9 @@
 	iomodel->Constant(&stabilization,MasstransportStabilizationEnum);
 
-	/*Do not add constraints in DG, they are weakly imposed*/
+	/*Do not add constraints in DG,  they are weakly imposed*/
 	if(stabilization!=3){
 		IoModelToConstraintsx(constraints,iomodel,MasstransportSpcthicknessEnum,MasstransportAnalysisEnum,P1Enum);
 	}
+	/*Do not add constraints in FCT, they are imposed using penalties*/
 }/*}}}*/
 void MasstransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
@@ -675,5 +679,5 @@
 }/*}}}*/
 void MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-	   _error_("not implemented yet");
+	element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
 }/*}}}*/
 void MasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
@@ -773,2 +777,147 @@
 	return;
 }/*}}}*/
+
+/*Flux Correction Transport*/
+void           MasstransportAnalysis::LumpedMassMatrix(Vector<IssmDouble>** pMlff,FemModel* femmodel){/*{{{*/
+
+	/*Intermediaries*/
+	int  configuration_type;
+
+	/*Initialize Lumped mass matrix (actually we just save its diagonal)*/
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	int fsize      = femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
+	int flocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
+	Vector<IssmDouble>* Mlff = new Vector<IssmDouble>(flocalsize,fsize);
+
+	/*Create and assemble matrix*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element*       element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ElementMatrix* MLe     = this->CreateMassMatrix(element);
+		if(MLe){
+			MLe->Lump();
+			MLe->AddDiagonalToGlobal(Mlff);
+		}
+		delete MLe;
+	}
+	Mlff->Assemble();
+
+	/*Assign output pointer*/
+	*pMlff=Mlff;
+}/*}}}*/
+ElementMatrix* MasstransportAnalysis::CreateMassMatrix(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries*/
+	IssmDouble  D,Jdet;
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Me     = element->NewElementMatrix();
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		D=gauss->weight*Jdet;
+		TripleMultiply(basis,1,numnodes,1,
+					&D,1,1,0,
+					basis,1,numnodes,0,
+					&Me->values[0],1);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return Me;
+}/*}}}*/
+void           MasstransportAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
+
+	/*Output*/
+	Matrix<IssmDouble>* Kff = NULL;
+	Matrix<IssmDouble>* Kfs = NULL;
+
+	/*Initialize Jacobian Matrix*/
+	AllocateSystemMatricesx(&Kff,&Kfs,NULL,NULL,femmodel);
+
+	/*Create and assemble matrix*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element*       element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ElementMatrix* Ke     = this->CreateFctKMatrix(element);
+		if(Ke) Ke->AddToGlobal(Kff,Kfs);
+		delete Ke;
+	}
+	Kff->Assemble();
+	Kfs->Assemble();
+
+	/*Assign output pointer*/
+	*pKff=Kff;
+	*pKfs=Kfs;
+}/*}}}*/
+ElementMatrix* MasstransportAnalysis::CreateFctKMatrix(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries */
+	IssmDouble Jdet;
+	IssmDouble vx,vy;
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int dim      = 2;
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+	IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
+	Input* vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		GetB(B,element,dim,xyz_list,gauss);
+		GetBprime(Bprime,element,dim,xyz_list,gauss);
+		vxaverage_input->GetInputValue(&vx,gauss);
+		vyaverage_input->GetInputValue(&vy,gauss);
+
+		D[0*dim+0] = -gauss->weight*vx*Jdet;
+		D[1*dim+1] = -gauss->weight*vy*Jdet;
+
+		TripleMultiply(B,dim,numnodes,1,
+					D,dim,dim,0,
+					Bprime,dim,numnodes,0,
+					&Ke->values[0],1);
+
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(D);
+	delete gauss;
+	return Ke;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 18343)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 18344)
@@ -36,4 +36,10 @@
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void UpdateConstraints(FemModel* femmodel);
+
+		/*FCT*/
+		void           LumpedMassMatrix(Vector<IssmDouble>** pMLff,FemModel* femmodel);
+		void           FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel);
+		ElementMatrix* CreateMassMatrix(Element* element);
+		ElementMatrix* CreateFctKMatrix(Element* element);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Constraints/Constraint.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Constraints/Constraint.h	(revision 18343)
+++ /issm/trunk-jpl/src/c/classes/Constraints/Constraint.h	(revision 18344)
@@ -21,4 +21,5 @@
 		virtual      ~Constraint(){};
 		virtual void ConstrainNode(Nodes* nodes,Parameters* parameters)=0;
+		virtual void PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters)=0;
 		virtual bool InAnalysis(int analysis_type)=0;
 
Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.h	(revision 18343)
+++ /issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.h	(revision 18344)
@@ -38,4 +38,5 @@
 		void   ConstrainNode(Nodes* nodes,Parameters* parameters);
 		bool   InAnalysis(int analysis_type);
+		void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters){_error_("not implemented yet");};
 		/*}}}*/
 		/*SpcDynamic management:{{{ */
Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.h	(revision 18343)
+++ /issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.h	(revision 18344)
@@ -37,4 +37,5 @@
 		void   ConstrainNode(Nodes* nodes,Parameters* parameters);
 		bool   InAnalysis(int analysis_type);
+		void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters){_error_("not implemented yet");};
 		/*}}}*/
 		/*SpcStatic management:{{{ */
Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp	(revision 18343)
+++ /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp	(revision 18344)
@@ -15,11 +15,12 @@
 /*SpcTransient constructors and destructor*/
 SpcTransient::SpcTransient(){/*{{{*/
-	sid=-1;
-	nodeid=-1;
-	dof=-1;
-	values=NULL;
-	times=NULL;
-	nsteps=-1;
-	analysis_type=-1;
+	penalty       = false;
+	sid           = -1;
+	nodeid        = -1;
+	dof           = -1;
+	values        = NULL;
+	times         = NULL;
+	nsteps        = -1;
+	analysis_type = -1;
 	return;
 }
@@ -27,11 +28,12 @@
 SpcTransient::SpcTransient(int spc_sid,int spc_nodeid, int spc_dof,int spc_nsteps, IssmDouble* spc_times, IssmDouble* spc_values,int spc_analysis_type){/*{{{*/
 
-	sid=spc_sid;
-	nodeid=spc_nodeid;
-	dof=spc_dof;
-	nsteps=spc_nsteps;
+	penalty= false;
+	sid    = spc_sid;
+	nodeid = spc_nodeid;
+	dof    = spc_dof;
+	nsteps = spc_nsteps;
 	if(spc_nsteps){
-		values=xNew<IssmDouble>(spc_nsteps);
-		times=xNew<IssmDouble>(spc_nsteps);
+		values = xNew<IssmDouble>(spc_nsteps);
+		times  = xNew<IssmDouble>(spc_nsteps);
 		xMemCpy<IssmDouble>(values,spc_values,nsteps);
 		xMemCpy<IssmDouble>(times,spc_times,nsteps);
@@ -69,5 +71,7 @@
 }		
 /*}}}*/
-int    SpcTransient::Id(void){ return sid; }/*{{{*/
+int  SpcTransient::Id(void){/*{{{*/
+	return sid;
+}
 /*}}}*/
 int SpcTransient::ObjectEnum(void){/*{{{*/
@@ -101,5 +105,5 @@
 	node=(Node*)nodes->GetObjectById(NULL,nodeid);
 
-	if(node){ //in case the spc is dealing with a node on another cpu
+	if(!this->penalty && node){ //in case the spc is dealing with a node on another cpu
 
 		/*Retrieve time in parameters: */
@@ -137,4 +141,60 @@
 }
 /*}}}*/
+void SpcTransient::PenaltyDofAndValue(int* pdof,IssmDouble* pvalue,Nodes* nodes,Parameters* parameters){/*{{{*/
+
+	if(!this->penalty) _error_("cannot return dof and value for non penalty constraint");
+
+	Node       *node  = NULL;
+	IssmDouble  time  = 0.;
+	int         i,gdof;
+	IssmDouble  alpha = -1.;
+	IssmDouble  value;
+	bool        found = false;
+
+	/*Chase through nodes and find the node to which this SpcTransient applys: */
+	node=(Node*)nodes->GetObjectById(NULL,nodeid);
+
+	if(node){ //in case the spc is dealing with a node on another cpu
+
+		/*Retrieve time in parameters: */
+		parameters->FindParam(&time,TimeEnum);
+
+		/*Now, go fetch value for this time: */
+		if (time<=times[0]){
+			value=values[0];
+			found=true;
+		}
+		else if (time>=times[nsteps-1]){
+			value=values[nsteps-1];
+			found=true;
+		}
+		else{
+			for(i=0;i<nsteps-1;i++){
+				if (times[i]<=time && time<times[i+1]){
+					alpha=(time-times[i])/(times[i+1]-times[i]);
+					value=(1-alpha)*values[i]+alpha*values[i+1];
+					found=true;
+					break;
+				}
+			}
+		}
+		if(!found)_error_("could not find time segment for constraint");
+
+		/*Get gdof */
+		gdof = node->GetDof(dof,GsetEnum);
+		if(xIsNan<IssmDouble>(value)){
+			gdof = -1;
+		}
+	}
+	else{
+		value = NAN;
+		gdof = -1;
+	}
+
+	/*Assign output pointers*/
+	*pdof   = gdof;
+	*pvalue = value;
+}
+/*}}}*/
 
 /*SpcTransient functions*/
Index: /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.h	(revision 18343)
+++ /issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.h	(revision 18344)
@@ -14,11 +14,12 @@
 
 	private: 
-		int	sid; /*! id, to track it*/
-		int	nodeid; /*!node id*/
-		int dof; /*!component*/
-		IssmDouble* values; /*different values in time*/
-		IssmDouble* times; /*different time steps*/
-		int nsteps; /*number of time steps*/
-		int analysis_type;
+		bool        penalty;         /*Is this a penalty constraint */
+		int         sid;             /* id, to track it             */
+		int         nodeid;          /*node id                      */
+		int         dof;             /*component                    */
+		IssmDouble *values;          /*different values in time     */
+		IssmDouble *times;           /*different time steps         */
+		int         nsteps;          /*number of time steps         */
+		int         analysis_type;
 
 	public:
@@ -39,4 +40,5 @@
 		void   ConstrainNode(Nodes* nodes,Parameters* parameters);
 		bool   InAnalysis(int analysis_type);
+		void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters);
 		/*}}}*/
 		/*SpcTransient management:{{{ */
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp	(revision 18343)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp	(revision 18344)
@@ -322,4 +322,37 @@
 		}
 
+	}
+	else{
+		_error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
+	}
+
+}
+/*}}}*/
+void ElementMatrix::AddDiagonalToGlobal(Vector<IssmDouble>* pf){/*{{{*/
+
+	IssmDouble* localvalues=NULL;
+
+	/*Check that pf is not NULL*/
+	_assert_(pf); 
+
+	/*In debugging mode, check consistency (no NaN, and values not too big)*/
+	this->CheckConsistency();
+
+	if(this->dofsymmetrical){
+		/*only use row dofs to add values into global matrices: */
+
+		if(this->row_fsize){
+			/*first, retrieve values that are in the f-set from the g-set values matrix: */
+			localvalues=xNew<IssmDouble>(this->row_fsize);
+			for(int i=0;i<this->row_fsize;i++){
+				localvalues[i] = this->values[this->ncols*this->row_flocaldoflist[i]+ this->row_flocaldoflist[i]];
+			}
+
+			/*add local values into global  matrix, using the fglobaldoflist: */
+			pf->SetValues(this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
+
+			/*Free ressources:*/
+			xDelete<IssmDouble>(localvalues);
+		}
 	}
 	else{
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h	(revision 18343)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h	(revision 18344)
@@ -60,4 +60,5 @@
 		void AddToGlobal(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
 		void AddToGlobal(Matrix<IssmDouble>* Jff);
+		void AddDiagonalToGlobal(Vector<IssmDouble>* pf);
 		void Echo(void);
 		void CheckConsistency(void);
Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp	(revision 18343)
+++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp	(revision 18344)
@@ -16,5 +16,5 @@
 	bool   save_results;
 	bool   isFS,isfreesurface,dakota_analysis;
-	int    solution_type;
+	int    solution_type,stabilization;
 	char** requested_outputs = NULL;
 
@@ -29,4 +29,5 @@
 	femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
 	femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
+	femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
 	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
 
@@ -49,5 +50,10 @@
 	else{
 		if(VerboseSolution()) _printf0_("   call computational core\n");
-		solutionsequence_linear(femmodel);
+		if(stabilization==4){
+			solutionsequence_fct(femmodel);
+		}
+		else{
+			solutionsequence_linear(femmodel);
+		}
 		femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
 		extrudefrombase_core(femmodel);
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18344)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18344)
@@ -0,0 +1,161 @@
+/*!\file: solutionsequence_linear.cpp
+ * \brief: numerical core of linear solutions
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+
+void solutionsequence_fct(FemModel* femmodel){
+
+	/*intermediary: */
+	Vector<IssmDouble>*  Mlf = NULL;
+	Matrix<IssmDouble>*  Kff = NULL;
+	Matrix<IssmDouble>*  Kfs = NULL;
+	Vector<IssmDouble>*  ug  = NULL;
+	Vector<IssmDouble>*  uf  = NULL;
+	Vector<IssmDouble>*  ys  = NULL;
+
+	IssmDouble theta,deltat;
+	int        dof,ncols,ncols2,ncols3,rstart,rend;
+	int        configuration_type,analysis_type;
+	double     d,diagD,mi;
+	double     dmax = 0.;
+	int*       cols  = NULL;
+	int*       cols2 = NULL;
+	int*       cols3 = NULL;
+	double*    vals  = NULL;
+	double*    vals2 = NULL;
+	double*    vals3 = NULL;
+
+	/*Create analysis*/
+	MasstransportAnalysis* analysis = new MasstransportAnalysis();
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	femmodel->UpdateConstraintsx();
+	theta = 0.5;
+
+	/*Create lumped mass matrix*/
+	analysis->LumpedMassMatrix(&Mlf,femmodel);
+	analysis->FctKMatrix(&Kff,&Kfs,femmodel);
+
+	/*Create Dff Matrix*/
+	#ifdef _HAVE_PETSC_
+	Mat Kff_transp = NULL;
+	Mat Dff_petsc  = NULL;
+	Mat LHS        = NULL;
+	Mat Kff_petsc  = Kff->pmatrix->matrix;
+	Vec Mlf_petsc  = Mlf->pvector->vector;
+	MatTranspose(Kff_petsc,MAT_INITIAL_MATRIX,&Kff_transp);
+	MatDuplicate(Kff_petsc,MAT_SHARE_NONZERO_PATTERN,&Dff_petsc);
+	MatGetOwnershipRange(Kff_transp,&rstart,&rend);
+	for(int row=rstart; row<rend; row++){
+		diagD = 0.;
+		MatGetRow(Kff_petsc ,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatGetRow(Kff_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+		_assert_(ncols==ncols2);
+		for(int j=0; j<ncols; j++) {
+			_assert_(cols[j]==cols2[j]);
+			d = max(max(-vals[j],-vals2[j]),0.);
+			MatSetValues(Dff_petsc,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
+			if(cols[j]!=row) diagD -= d;
+		}
+		MatSetValues(Dff_petsc,1,&row,1,&row,(const double*)&diagD,INSERT_VALUES);
+		MatRestoreRow(Kff_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatRestoreRow(Kff_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+	}
+	MatAssemblyBegin(Dff_petsc,MAT_FINAL_ASSEMBLY);
+	MatAssemblyEnd(  Dff_petsc,MAT_FINAL_ASSEMBLY);
+	MatFree(&Kff_transp);
+
+	/*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
+	MatDuplicate(Kff_petsc,MAT_SHARE_NONZERO_PATTERN,&LHS);
+	for(int row=rstart; row<rend; row++){
+		MatGetRow(Kff_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatGetRow(Dff_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+		_assert_(ncols==ncols2);
+		for(int j=0; j<ncols; j++) {
+			_assert_(cols[j]==cols2[j]);
+			d = -theta*deltat*(vals[j] + vals2[j]);
+			if(cols[j]==row){
+				VecGetValues(Mlf_petsc,1,(const int*)&cols[j],&mi);
+				d += mi;
+			}
+			if(fabs(d)>dmax) dmax = fabs(d);
+			MatSetValues(LHS,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
+		}
+		MatRestoreRow(Kff_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatRestoreRow(Dff_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+	}
+
+	/*Penalize Dirichlet boundary*/
+	dmax = dmax * 1.e+3;
+	for(int i=0;i<femmodel->constraints->Size();i++){
+		Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+		if(constraint->InAnalysis(analysis_type)){
+			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+			if(dof!=-1){
+				MatSetValues(LHS,1,&dof,1,&dof,(const double*)&dmax,INSERT_VALUES);
+			}
+		}
+	}
+	MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
+	MatAssemblyEnd(  LHS,MAT_FINAL_ASSEMBLY);
+
+	/*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
+	GetSolutionFromInputsx(&ug,femmodel);
+	Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
+	delete ug;
+	Vec u  = uf->pvector->vector;
+	Vec Ku = NULL;
+	Vec Du = NULL;
+	Vec RHS = NULL;
+	VecDuplicate(u,&Ku);
+	VecDuplicate(u,&Du);
+	VecDuplicate(u,&RHS);
+	MatMult(Kff_petsc,u,Ku);
+	MatMult(Dff_petsc,u,Du);
+	VecPointwiseMult(RHS,Mlf_petsc,u);
+	VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);// RHS = M*u + (1-theta)*deltat*K*u + (1-theta)*deltat*D*u
+	VecFree(&Ku);
+	VecFree(&Du);
+	MatFree(&Dff_petsc);
+	delete uf;
+
+	/*Penalize Dirichlet boundary*/
+	for(int i=0;i<femmodel->constraints->Size();i++){
+		Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
+		if(constraint->InAnalysis(analysis_type)){
+			constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
+			d = d*dmax;
+			if(dof!=-1){
+				VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
+			}
+		}
+	}
+
+	/*Go solve!*/
+	SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
+	MatFree(&LHS);
+	VecFree(&RHS);
+	uf =new Vector<IssmDouble>(u);
+	VecFree(&u);
+
+	Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
+	InputUpdateFromSolutionx(femmodel,ug); 
+	delete ug;  
+
+	#else
+	_error_("PETSc needs to be installed");
+	#endif
+
+	delete Mlf;
+	delete Kff;
+	delete Kfs;
+	delete analysis;
+
+}
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h	(revision 18343)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h	(revision 18344)
@@ -16,4 +16,5 @@
 void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solutionsequence_newton(FemModel* femmodel);
+void solutionsequence_fct(FemModel* femmodel);
 void solutionsequence_FScoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solutionsequence_linear(FemModel* femmodel);