Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18364)
+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18365)
@@ -93,6 +93,12 @@
 		MatRestoreRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
 	}
 
+	/*Broadcast max(dmax)*/
+	IssmDouble dmax_all;
+	ISSM_MPI_Reduce(&dmax,&dmax_all,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
+	ISSM_MPI_Bcast(&dmax_all,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	dmax = dmax_all;
+
 	/*Penalize Dirichlet boundary*/
 	dmax = dmax * 1.e+3;
 	for(int i=0;i<femmodel->constraints->Size();i++){
@@ -209,6 +215,14 @@
 	VecDuplicate(uvec,&Ri_plus);
 	VecDuplicate(uvec,&Ri_minus);
 
+	/*Get global extremas*/
+	IssmDouble uLmin_global =  ulmin[0];
+	IssmDouble uLmax_global =  ulmax[0];
+	VecGetSize(uvec,&ncols);
+	for(int i=1;i<ncols;i++) if(ulmin[i]<uLmin_global) uLmin_global = ulmin[i];
+	for(int i=1;i<ncols;i++) if(ulmax[i]>uLmax_global) uLmax_global = ulmax[i];
+	//printf("%g %g\n",uLmin_global,uLmax_global);
+
 	/*Go through D and calculate Ris*/
 	MatGetOwnershipRange(D,&rstart,&rend);
 	for(int row=rstart; row<rend; row++){
@@ -233,8 +247,10 @@
 		/*Compute Qis and Ris*/
 		//double Qi_plus  = ml_serial[row]/deltat*(3. - u[row]);
 		//double Qi_minus = ml_serial[row]/deltat*(2. - u[row]);
-		double Qi_plus  = ml_serial[row]/deltat*(ulmax[row] - u[row]);
-		double Qi_minus = ml_serial[row]/deltat*(ulmin[row] - u[row]);
+		//double Qi_plus  = ml_serial[row]/deltat*(ulmax[row] - u[row]);
+		//double Qi_minus = ml_serial[row]/deltat*(ulmin[row] - u[row]);
+		double Qi_plus  = ml_serial[row]/deltat*(uLmax_global - u[row]);
+		double Qi_minus = ml_serial[row]/deltat*(uLmin_global - u[row]);
 		_assert_(Qi_plus  >= 0.);
 		_assert_(Qi_minus <= 0.);
 		d = 1.;
@@ -371,7 +387,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);
 	delete uf;
-
+	
 	/*Go solve lower order solution*/
 	SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
 	MatFree(&LHS);
