Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18364)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp	(revision 18365)
@@ -94,4 +94,10 @@
 	}
 
+	/*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;
@@ -209,4 +215,12 @@
 	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*/
@@ -234,6 +248,8 @@
 		//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.);
@@ -372,5 +388,5 @@
 	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);