Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18346)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18347)
@@ -167,39 +167,94 @@
 	IssmDouble* udot_serial = NULL;
 	IssmDouble* u_serial    = NULL;
-	VecToMPISerial(&udot_serial,udot,IssmComm::GetComm());
-	VecToMPISerial(&u_serial   ,u   ,IssmComm::GetComm());
-
+	IssmDouble* ml_serial    = NULL;
+	VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
+	VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
+	VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
 
 	/*Anti diffusive fluxes*/
-	Mat F    = NULL;
+	Vec Ri_plus  = NULL;
+	Vec Ri_minus = NULL;
+	double uiLmax = 3.;
+	double uiLmin = 2.;
+	VecDuplicate(u,&Ri_plus);
+	VecDuplicate(u,&Ri_minus);
+	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*/
+	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);
-	MatDuplicate(D_petsc,MAT_SHARE_NONZERO_PATTERN,&F);
 	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);
-		mi = 0.;
-		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]) mi += d;
-			MatSetValues(F,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
-
-		}
-		VecSetValue(Fbar,row,(const double)mi,INSERT_VALUES);
+		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);
 	}
-	MatAssemblyBegin(F,MAT_FINAL_ASSEMBLY);
-	MatAssemblyEnd(  F,MAT_FINAL_ASSEMBLY);
 	VecAssemblyBegin(Fbar);
 	VecAssemblyEnd(  Fbar);
 
 	MatFree(&D_petsc);
-	MatFree(&F);
 	delete Mc;
 	xDelete<IssmDouble>(udot_serial);
 	xDelete<IssmDouble>(u_serial);
+	xDelete<IssmDouble>(ml_serial);
+	xDelete<IssmDouble>(Ri_plus_serial);
+	xDelete<IssmDouble>(Ri_minus_serial);
 
 	/*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/