Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26150)
+++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26151)
@@ -411,4 +411,5 @@
 	Vector<IssmDouble>*    oceanareas=NULL;
 	IssmDouble             oceanarea;
+	IssmDouble             oceanaverage;
 	bool                   scaleoceanarea=false;
 	IssmDouble             rho_water;
@@ -496,7 +497,8 @@
 	}
 	sealevelloads->Assemble();
-
+	oceanareas->Assemble();
+	
 	//Get ocean area: 
-	oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
+	oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
 	if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
 
@@ -523,5 +525,6 @@
 	
 		//substract ocean average and barystatic contribution
-		sealevelloads->Shift(barycontrib->Total()/oceanarea- SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
+		oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
+		sealevelloads->Shift(barycontrib->Total()/oceanarea- oceanaverage);
 		
 		//broadcast sea level loads 
@@ -554,5 +557,5 @@
 	/*Update bedrock motion and geoid:*/
 	if(computesealevel){
-		femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/oceanarea- SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
+		femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water);
 
 		/*cumulate barystatic contributions and save to results: */