Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26151)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26152)
@@ -398,4 +398,5 @@
 		virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector,  SealevelMasks* masks)=0;
 		virtual void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0;
+		virtual void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks)=0;
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26151)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26152)
@@ -235,4 +235,5 @@
 		void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
 		void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 26151)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 26152)
@@ -190,4 +190,5 @@
 		void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
 		void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
 
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 26151)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 26152)
@@ -196,4 +196,5 @@
 		void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
 		void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
 
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26151)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26152)
@@ -6965,4 +6965,11 @@
 	return;
 }/*}}}*/
+void       Tria::SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){ /*{{{*/
+	
+	if(masks->isoceanin[this->lid]){
+		loads->SetValue(this->sid,offset,ADD_VAL);
+	}
+
+} /*}}}*/
 #endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 26151)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 26152)
@@ -180,4 +180,5 @@
 		void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);
 		void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks);
+		void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks);
 		#endif
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/SealevelMasks.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/SealevelMasks.h	(revision 26151)
+++ /issm/trunk-jpl/src/c/classes/SealevelMasks.h	(revision 26152)
@@ -18,10 +18,10 @@
 		
 		/*SealevelMasks constructors, destructors :*/
-		SealevelMasks(int nel){
+		SealevelMasks(int localnel){
 			/*allocate fields:*/
-			this->isiceonly=xNew<bool>(nel);
-			this->isfullyfloating=xNew<bool>(nel);
-			this->notfullygrounded=xNew<bool>(nel);
-			this->isoceanin=xNew<bool>(nel);
+			this->isiceonly=xNew<bool>(localnel);
+			this->isfullyfloating=xNew<bool>(localnel);
+			this->notfullygrounded=xNew<bool>(localnel);
+			this->isoceanin=xNew<bool>(localnel);
 		};
 		~SealevelMasks(){
Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26151)
+++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26152)
@@ -12,4 +12,5 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+
 /*support routines local definitions:{{{*/
 void TransferForcing(FemModel* femmodel,int forcingenum);
@@ -19,4 +20,5 @@
 SealevelMasks* sealevel_masks(FemModel* femmodel);
 void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads);
+void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, IssmDouble offset, SealevelMasks* masks);
 /*}}}*/
 
@@ -503,9 +505,12 @@
 	if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
 
-	//substract ocean average and barystatic contributionfrom sea level loads: 
-	sealevelloads->Shift(barycontrib->Total()/oceanarea - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
-	
+	//conserve ocean mass: 
+	oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
+	ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks);
+
 	//broadcast sea level loads 
 	allsealevelloads=sealevelloads->ToMPISerial();
+
+	femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,DummyEnum,allsealevelloads,nel,1,1,1));
 
 	//compute rotation axis motion:
@@ -523,9 +528,9 @@
 		}
 		sealevelloads->Assemble();
-	
-		//substract ocean average and barystatic contribution
+			
+		//Conserve ocean mass: 
 		oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
-		sealevelloads->Shift(barycontrib->Total()/oceanarea- oceanaverage);
-		
+		ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks);
+
 		//broadcast sea level loads 
 		allsealevelloads=sealevelloads->ToMPISerial();
@@ -1318,2 +1323,12 @@
 	m[2]=m3;
 } /*}}}*/
+void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, IssmDouble offset, SealevelMasks* masks){ /*{{{*/
+
+	/*Shift sealevel loads by ocean average, only on ocean! :*/
+	for(Object* & object : femmodel->elements->objects){
+		Element* element = xDynamicCast<Element*>(object);
+		element->SealevelchangeShift(sealevelloads,offset,masks);
+	}
+	sealevelloads->Assemble();
+
+} /*}}}*/