Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 20009)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 20010)
@@ -2216,5 +2216,5 @@
 	IssmDouble  eartharea=0;
 	IssmDouble  eartharea_cpu=0;
-	int         ns;
+	int         ns,nsmax;
 	
 	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
@@ -2233,12 +2233,19 @@
 	ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
+	/*Figure out max of ns: */
+	ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
 	/*Call the sea level rise core: */
-	for(int i=0;i<ns;i++){
+	for(int i=0;i<nsmax;i++){
+		if(i<ns){
 		
-		if(IssmComm::GetRank()==0)if(VerboseConvergence())if(i%1000)_printf_("\r" << "      convolution progress: " << (float)i/(float)ns*100 << "\%");
-	
-		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-		element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
-		eustatic_cpu+=eustatic_cpu_e;
+			if(IssmComm::GetRank()==0)if(VerboseConvergence())if(i%1000)_printf_("\r" << "      convolution progress: " << (float)i/(float)ns*100 << "\%");
+		
+			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+			element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
+			eustatic_cpu+=eustatic_cpu_e;
+		}
+		pSgi->Assemble();
 	}
 	if(IssmComm::GetRank()==0)if(VerboseConvergence())_printf_("\n");
@@ -2263,5 +2270,5 @@
 	IssmDouble  eartharea_cpu=0;
 
-	int         ns;
+	int         ns,nsmax;
 	
 	/*Serialize vectors from previous iteration:*/
@@ -2270,5 +2277,5 @@
 	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
 	ns = elements->Size();
-
+	
 	/*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
 	for(int i=0;i<ns;i++){
@@ -2283,14 +2290,20 @@
 	ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
+	/*Figure out max of ns: */
+	ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
+	ISSM_MPI_Bcast(&nsmax,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
 	/*Call the sea level rise core: */
-	for(int i=0;i<ns;i++){
-
-		if(verboseconvolution)if(IssmComm::GetRank()==0)if(VerboseConvergence())if(i%1000)_printf_("\r" << "      convolution progress: " << (float)i/(float)ns*100 << "\%");
-		
-		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-		element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,oceanarea,eartharea);
+	for(int i=0;i<nsmax;i++){
+		if(i<ns){
+
+			if(verboseconvolution)if(IssmComm::GetRank()==0)if(VerboseConvergence())if(i%1000)_printf_("\r" << "      convolution progress: " << (float)i/(float)ns*100 << "\%");
+
+			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+			element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,oceanarea,eartharea);
+		}
+		pSgo->Assemble();
 	}
 	if(verboseconvolution)if(IssmComm::GetRank()==0)if(VerboseConvergence())_printf_("\n");
-
 
 	/*Free ressources:*/