Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 24730)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 24731)
@@ -144,5 +144,5 @@
 	*pLHS  = LHS;
 }/*}}}*/
-void CreateRHS(Vec* pRHS,Mat K,Mat D,Vec Ml,Vec u,IssmDouble theta,IssmDouble deltat,IssmDouble dmax,FemModel* femmodel,int configuration_type){/*{{{*/
+void CreateRHS(Vec* pRHS,Mat K,Mat D,Vec Ml,Vec p,Vec u,IssmDouble theta,IssmDouble deltat,IssmDouble dmax,FemModel* femmodel,int configuration_type){/*{{{*/
 	/*Create Right Hand side of Lower order solution
 	 *
@@ -170,4 +170,5 @@
 	VecPointwiseMult(RHS,Ml,u);
 	VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);
+	VecAXPBY(RHS,deltat,1,p);
 	VecFree(&Ku);
 	VecFree(&Du);
@@ -392,4 +393,5 @@
 	Matrix<IssmDouble>*  K  = NULL;
 	Matrix<IssmDouble>*  Mc = NULL;
+	Vector<IssmDouble>*  p  = NULL;
 	Vector<IssmDouble>*  ug = NULL;
 	Vector<IssmDouble>*  uf = NULL;
@@ -401,5 +403,5 @@
 	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
 	femmodel->UpdateConstraintsx();
-	theta = 0.5;
+	theta = 0.5; //0.5==Crank-Nicolson, 1==Backward Euler
 
 	/*Create analysis*/
@@ -410,4 +412,5 @@
 			manalysis->MassMatrix(&Mc,femmodel);
 			manalysis->FctKMatrix(&K,NULL,femmodel);
+			manalysis->FctPVector(&p,femmodel);
 			break;
 		case DamageEvolutionAnalysisEnum:
@@ -431,4 +434,5 @@
 	Mat K_petsc  = K->pmatrix->matrix;
 	Mat Mc_petsc = Mc->pmatrix->matrix;
+	Vec p_petsc	 = p->pvector->vector;
 
 	/*Get previous solution u^n*/
@@ -449,6 +453,7 @@
 
 	/*Create RHS: [ML + (1 - theta) deltaT L^n] u^n */
-	CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
+	CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,p_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
 	delete uf;
+	delete p;
 
 	/*Go solve lower order solution*/