Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 21078)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 21079)
@@ -3631,5 +3631,5 @@
 
 	/*Initialize eustatic component: do not skip this step :):*/
-	IssmDouble eustatic = 0;
+	IssmDouble eustatic = 0.;
 
 	/*Computational flags:*/
@@ -3710,4 +3710,5 @@
 
 	/*Compute eustatic compoent:*/
+	_assert_(oceanarea>0.);
 	eustatic += rho_ice*area*I/(oceanarea*rho_water); 
 
@@ -3749,4 +3750,6 @@
 	
 	/*Assign output pointer:*/
+	_assert_(!xIsNan<IssmDouble>(eustatic));
+	_assert_(!xIsInf<IssmDouble>(eustatic));
 	*peustatic=eustatic;
 	return;
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21078)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21079)
@@ -2296,11 +2296,11 @@
 
 	/*serialized vectors:*/
-	IssmDouble  eustatic=0;
-	IssmDouble  eustatic_cpu=0;
-	IssmDouble  eustatic_cpu_e=0;
-	IssmDouble  oceanarea=0;
-	IssmDouble  oceanarea_cpu=0;
-	IssmDouble  eartharea=0;
-	IssmDouble  eartharea_cpu=0;
+	IssmDouble  eustatic       = 0.;
+	IssmDouble  eustatic_cpu   = 0.;
+	IssmDouble  eustatic_cpu_e = 0.;
+	IssmDouble  oceanarea      = 0.;
+	IssmDouble  oceanarea_cpu  = 0.;
+	IssmDouble  eartharea      = 0.;
+	IssmDouble  eartharea_cpu  = 0.;
 	int         ns,nsmax;
 	
@@ -2312,8 +2312,9 @@
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 		oceanarea_cpu += element->OceanArea();
-		eartharea_cpu+= element->GetAreaSpherical();
+		eartharea_cpu += element->GetAreaSpherical();
 	}
 	ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
 	ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	_assert_(oceanarea>0.);
 
 	ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
@@ -2344,4 +2345,5 @@
 	ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
 	ISSM_MPI_Bcast(&eustatic,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	_assert_(!xIsNan<IssmDouble>(eustatic));
 
 	/*Assign output pointers:*/
@@ -2372,5 +2374,5 @@
 		Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
 		oceanarea_cpu += element->OceanArea();
-		eartharea_cpu+= element->GetAreaSpherical();
+		eartharea_cpu += element->GetAreaSpherical();
 	}
 	ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp	(revision 21078)
+++ /issm/trunk-jpl/src/c/cores/sealevelrise_core_eustatic.cpp	(revision 21079)
@@ -19,10 +19,10 @@
 	int  gsize;
 	bool spherical=true;
-	IssmDouble          *latitude    = NULL;
-	IssmDouble          *longitude    = NULL;
-	IssmDouble          *radius    = NULL;
+	IssmDouble *latitude  = NULL;
+	IssmDouble *longitude = NULL;
+	IssmDouble *radius    = NULL;
 
 	/*outputs:*/
-	IssmDouble eustatic=0;
+	IssmDouble eustatic;
 
 	/*recover parameters:*/
@@ -33,5 +33,5 @@
 
 	/*Figure out size of g-set deflection vector and allocate solution vector: */
-	gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
+	gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
 	
 	/*Initialize:*/
@@ -51,9 +51,8 @@
 	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelEustaticEnum,-eustatic));
 
+	/*clean up and return:*/
 	xDelete<IssmDouble>(latitude);
 	xDelete<IssmDouble>(longitude);
 	xDelete<IssmDouble>(radius);
-
-	/*return:*/
 	return Sgi;
 }
Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp	(revision 21078)
+++ /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp	(revision 21079)
@@ -46,5 +46,5 @@
 
 	/*Figure out size of g-set deflection vector and allocate solution vector: */
-	gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
+	gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
 	
 	/*Initialize:*/
@@ -55,7 +55,4 @@
 	Sg_old = new Vector<IssmDouble>(gsize);
 	Sg_old->Assemble();
-
-	Sgo = new Vector<IssmDouble>(gsize);
-	Sgo->Assemble();
 
 	count=1;
@@ -69,9 +66,9 @@
 
 		/*Initialize solution vector: */
-		Sg = new Vector<IssmDouble>(gsize); Sg->Assemble();
+		Sg  = new Vector<IssmDouble>(gsize); Sg->Assemble();
 		Sgo = new Vector<IssmDouble>(gsize); Sgo->Assemble();
 
 		/*call the non eustatic module: */
-        femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution);
+		femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution);
 	
 		/*assemble solution vector: */