Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18352)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18353)
@@ -194,4 +194,113 @@
 
 }/*}}}*/
+void CreateRis(IssmDouble** pRi_plus_serial,IssmDouble** pRi_minus_serial,Mat Mc,Mat D,IssmDouble* ml_serial,Vec uvec,IssmDouble* u,IssmDouble* udot,IssmDouble deltat){/*{{{*/
+
+	/*Intermediaries*/
+	Vec         Ri_plus  = NULL;
+	Vec         Ri_minus = NULL;
+	IssmDouble  uiLmax   = 3.;
+	IssmDouble  uiLmin   = 2.;
+	int         ncols,ncols2,rstart,rend;
+	double      d;
+	int        *cols     = NULL;
+	int        *cols2    = NULL;
+	double     *vals     = NULL;
+	double     *vals2    = NULL;
+
+	/*Initialize vectors*/
+	VecDuplicate(uvec,&Ri_plus);
+	VecDuplicate(uvec,&Ri_minus);
+
+	/*Go through D and calculate Ris*/
+	MatGetOwnershipRange(D,&rstart,&rend);
+	for(int row=rstart; row<rend; row++){
+		double Pi_plus  = 0.;
+		double Pi_minus = 0.;
+		MatGetRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatGetRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+		_assert_(ncols==ncols2);
+		for(int j=0; j<ncols; j++) {
+			_assert_(cols[j]==cols2[j]);
+			d = vals[j]*(udot[row] - udot[cols[j]]) + vals2[j]*(u[row] - u[cols[j]]);
+			if(row!=cols[j]){
+				if(d>0.){
+					Pi_plus  += d;
+				}
+				else{
+					Pi_minus += d;
+				}
+			}
+		}
+
+		/*Compute Qis and Ris*/
+		double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u[row]);
+		double Qi_minus = ml_serial[row]/deltat*(uiLmin - u[row]);
+		d = 1.;
+		if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
+		VecSetValue(Ri_plus,row,(const double)d,INSERT_VALUES);
+		d = 1.;
+		if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
+		VecSetValue(Ri_minus,row,(const double)d,INSERT_VALUES);
+
+		MatRestoreRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatRestoreRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+	}
+	VecAssemblyBegin(Ri_plus);
+	VecAssemblyEnd(  Ri_plus);
+	VecAssemblyBegin(Ri_minus);
+	VecAssemblyEnd(  Ri_minus);
+
+	/*Serialize Ris*/
+	IssmDouble* Ri_plus_serial  = NULL;
+	IssmDouble* Ri_minus_serial = NULL;
+	VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
+	VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
+	VecFree(&Ri_plus);
+	VecFree(&Ri_minus);
+
+	/*Assign output pointer*/
+	*pRi_plus_serial  = Ri_plus_serial;
+	*pRi_minus_serial = Ri_minus_serial;
+}/*}}}*/
+void CreateFbar(Vec* pFbar,IssmDouble* Ri_plus,IssmDouble* Ri_minus,Mat Mc,Mat D,IssmDouble* udot,IssmDouble* u,Vec uvec){/*{{{*/
+
+	/*Intermediaries*/
+	Vec         Fbar = NULL;
+	int         ncols,ncols2,rstart,rend;
+	double      d,f_ij;
+	int        *cols     = NULL;
+	int        *cols2    = NULL;
+	double     *vals     = NULL;
+	double     *vals2    = NULL;
+
+	/*Build fbar*/
+	VecDuplicate(uvec,&Fbar);
+	MatGetOwnershipRange(D,&rstart,&rend);
+	for(int row=rstart; row<rend; row++){
+		MatGetRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatGetRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+		_assert_(ncols==ncols2);
+		d = 0.;
+		for(int j=0; j<ncols; j++) {
+			_assert_(cols[j]==cols2[j]);
+			if(row==cols[j]) continue;
+			f_ij = vals[j]*(udot[row] - udot[cols[j]]) + vals2[j]*(u[row] - u[cols[j]]);
+			if(f_ij>0){
+				d += min(Ri_plus[row],Ri_minus[cols[j]]) * f_ij;
+			}
+			else{
+				d += min(Ri_minus[row],Ri_plus[cols[j]]) * f_ij;
+			}
+		}
+		VecSetValue(Fbar,row,(const double)d,INSERT_VALUES);
+		MatRestoreRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
+		MatRestoreRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
+	}
+	VecAssemblyBegin(Fbar);
+	VecAssemblyEnd(  Fbar);
+
+	/*Assign output pointer*/
+	*pFbar = Fbar;
+}/*}}}*/
 #endif
 	void solutionsequence_fct(FemModel* femmodel){
@@ -253,5 +362,5 @@
 		delete uf;
 
-		/*Go solve!*/
+		/*Go solve lower order solution*/
 		SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
 		MatFree(&LHS);
@@ -264,5 +373,5 @@
 		IssmDouble* udot_serial = NULL;
 		IssmDouble* u_serial    = NULL;
-		IssmDouble* ml_serial    = NULL;
+		IssmDouble* ml_serial   = NULL;
 		VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
 		VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
@@ -270,82 +379,11 @@
 
 		/*Anti diffusive fluxes*/
-		Vec Ri_plus  = NULL;
-		Vec Ri_minus = NULL;
-		double uiLmax = 3.;
-		double uiLmin = 2.;
-		VecDuplicate(u,&Ri_plus);
-		VecDuplicate(u,&Ri_minus);
-		MatGetOwnershipRange(K_petsc,&rstart,&rend);
-		for(int row=rstart; row<rend; row++){
-			double Pi_plus  = 0.;
-			double Pi_minus = 0.;
-			MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
-			MatGetRow(D_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 = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
-				if(row!=cols[j]){
-					if(d>0.){
-						Pi_plus  += d;
-					}
-					else{
-						Pi_minus += d;
-					}
-				}
-			}
-
-			/*Compute Qis and Ris*/
-			double Qi_plus  = ml_serial[row]/deltat*(uiLmax - u_serial[row]);
-			double Qi_minus = ml_serial[row]/deltat*(uiLmin - u_serial[row]);
-			d = 1.;
-			if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
-			VecSetValue(Ri_plus,row,(const double)d,INSERT_VALUES);
-			d = 1.;
-			if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
-			VecSetValue(Ri_minus,row,(const double)d,INSERT_VALUES);
-
-			MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
-			MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-		}
-		VecAssemblyBegin(Ri_plus);
-		VecAssemblyEnd(  Ri_plus);
-		VecAssemblyBegin(Ri_minus);
-		VecAssemblyEnd(  Ri_minus);
-
-		/*Serialize Ris*/
+		Vec         Fbar            = NULL;
+		IssmDouble* Ri_minus_serial = NULL;
 		IssmDouble* Ri_plus_serial  = NULL;
-		IssmDouble* Ri_minus_serial = NULL;
-		VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
-		VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
-		VecFree(&Ri_plus);
-		VecFree(&Ri_minus);
-
-		/*Build fbar*/
-		Vec Fbar = NULL;
-		VecDuplicate(u,&Fbar);
-		for(int row=rstart; row<rend; row++){
-			MatGetRow(Mc_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
-			MatGetRow(D_petsc ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-			_assert_(ncols==ncols2);
-			d = 0.;
-			for(int j=0; j<ncols; j++) {
-				_assert_(cols[j]==cols2[j]);
-				if(row==cols[j]) continue;
-				double f_ij = vals[j]*(udot_serial[row] - udot_serial[cols[j]]) + vals2[j]*(u_serial[row] - u_serial[cols[j]]);
-				if(f_ij>0){
-					d += min(Ri_plus_serial[row],Ri_minus_serial[cols[j]]) * f_ij;
-				}
-				else{
-					d += min(Ri_minus_serial[row],Ri_plus_serial[cols[j]]) * f_ij;
-				}
-			}
-			VecSetValue(Fbar,row,(const double)d,INSERT_VALUES);
-			MatRestoreRow(Mc_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
-			MatRestoreRow(D_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
-		}
-		VecAssemblyBegin(Fbar);
-		VecAssemblyEnd(  Fbar);
-
+		CreateRis(&Ri_plus_serial,&Ri_minus_serial,Mc_petsc,D_petsc,ml_serial,u,u_serial,udot_serial,deltat);
+		CreateFbar(&Fbar,Ri_plus_serial,Ri_minus_serial,Mc_petsc,D_petsc,udot_serial,u_serial,u);
+
+		/*Clean up*/
 		MatFree(&D_petsc);
 		delete Mc;
@@ -357,5 +395,5 @@
 
 		/*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
-		Vec temp1 = NULL; VecDuplicate(u,&temp1);
+		Vec temp1 = NULL;
 		VecDuplicate(u,&temp1);
 		VecPointwiseDivide(temp1,Fbar,Ml_petsc); //temp1 = Ml^-1 temp2