Index: issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 20155)
+++ issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 20156)
@@ -673,10 +673,11 @@
 	IssmDouble* newthickness   = xNew<IssmDouble>(numnodes);
 	IssmDouble* deltathickness = xNew<IssmDouble>(numnodes);
-	IssmDouble* newbed         = xNew<IssmDouble>(numnodes);
+	IssmDouble* newbase        = xNew<IssmDouble>(numnodes);
 	IssmDouble* newsurface     = xNew<IssmDouble>(numnodes);
 	IssmDouble* oldthickness   = xNew<IssmDouble>(numnodes);
-	IssmDouble* oldbed         = xNew<IssmDouble>(numnodes);
+	IssmDouble* oldbase        = xNew<IssmDouble>(numnodes);
 	IssmDouble* oldsurface     = xNew<IssmDouble>(numnodes);
 	IssmDouble* phi            = xNew<IssmDouble>(numnodes);
+	IssmDouble* sealevel       = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
@@ -689,12 +690,12 @@
 	}
 
-	/*Get previous bed, thickness and surface*/
-	basalelement->GetInputListOnNodes(&oldbed[0],BaseEnum);
+	/*Get previous base, thickness, surfac and current sealevel:*/
+	basalelement->GetInputListOnNodes(&oldbase[0],BaseEnum);
 	basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
 	basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
 	basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
-
-	
-	/*What is the delta thickness, useful for Sea-level rise:*/
+	basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum);
+
+	/*What is the delta thickness forcing the sea-level rise core: */
 	for(i=0;i<numnodes;i++) deltathickness[i]=newthickness[i]-oldthickness[i];
 
@@ -705,17 +706,16 @@
 
 	for(i=0;i<numnodes;i++) {
-		/*If shelf: hydrostatic equilibrium*/
-		if (phi[i]>0.){
-			newsurface[i] = oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
-			newbed[i]     = oldbed[i];                 //same bed: do nothing
-		}
-		else{ //this is an ice shelf
+		if (phi[i]>0.){ //this is an ice sheet: just add thickness to base.
+			newsurface[i] = oldbase[i]+newthickness[i]; //surface = oldbase + newthickness
+			newbase[i]     = oldbase[i];                 //same base: do nothing
+		}
+		else{ //this is an ice shelf: hydrostatic equilibrium*/
 			if(hydroadjustment==AbsoluteEnum){
-				newsurface[i] = newthickness[i]*(1-rho_ice/rho_water);
-				newbed[i]     = newthickness[i]*(-rho_ice/rho_water);
+				newsurface[i] = newthickness[i]*(1-rho_ice/rho_water)+sealevel[i];
+				newbase[i]     = newthickness[i]*(-rho_ice/rho_water)+sealevel[i];
 			}
 			else if(hydroadjustment==IncrementalEnum){
-				newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
-				newbed[i]     = oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed               = oldbed + di * dH
+				newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i])+sealevel[i]; //surface = oldsurface + (1-di) * dH
+				newbase[i]     = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base               = oldbed + di * dH
 			}
 			else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
@@ -727,13 +727,13 @@
 	element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum);
 	element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
-	element->AddBasalInput(BaseEnum,newbed,P1Enum);
+	element->AddBasalInput(BaseEnum,newbase,P1Enum);
 
 	/*Free ressources:*/
 	xDelete<IssmDouble>(newthickness);
-	xDelete<IssmDouble>(newbed);
+	xDelete<IssmDouble>(newbase);
 	xDelete<IssmDouble>(newsurface);
 	xDelete<IssmDouble>(oldthickness);
 	xDelete<IssmDouble>(deltathickness);
-	xDelete<IssmDouble>(oldbed);
+	xDelete<IssmDouble>(oldbase);
 	xDelete<IssmDouble>(oldsurface);
 	xDelete<IssmDouble>(phi);
Index: issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 20155)
+++ issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 20156)
@@ -1610,5 +1610,5 @@
 	/*Intermediaries*/
 	int         dim,domaintype;
-	IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
+	IssmDouble  Jdet,thickness,base,sealevel,water_pressure,ice_pressure;
 	IssmDouble  surface_under_water,base_under_water,pressure;
 	IssmDouble *xyz_list = NULL;
@@ -1635,4 +1635,5 @@
 	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
 	Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
+	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
 	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
@@ -1648,10 +1649,11 @@
 		gauss->GaussPoint(ig);
 		thickness_input->GetInputValue(&thickness,gauss);
-		base_input->GetInputValue(&bed,gauss);
+		sealevel_input->GetInputValue(&sealevel,gauss);
+		base_input->GetInputValue(&base,gauss);
 		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
 		element->NodalFunctions(basis,gauss);
 
-		surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
-		base_under_water    = min(0.,bed);           // 0 if the bottom of the glacier is above water level
+		surface_under_water = min(0.,thickness+base-sealevel); // 0 if the top of the glacier is above water level
+		base_under_water    = min(0.,base-sealevel);           // 0 if the bottom of the glacier is above water level
 		water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
 		ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
@@ -2070,5 +2072,5 @@
 
 	/*Intermediaries*/
-	IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
+	IssmDouble  Jdet,thickness,bed,sealevel,water_pressure,ice_pressure;
 	IssmDouble  surface_under_water,base_under_water,pressure;
 	IssmDouble *xyz_list = NULL;
@@ -2086,4 +2088,5 @@
 	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
 	Input* base_input       = element->GetInput(BaseEnum);       _assert_(base_input);
+	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
 	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
@@ -2100,9 +2103,10 @@
 		thickness_input->GetInputValue(&thickness,gauss);
 		base_input->GetInputValue(&bed,gauss);
+		sealevel_input->GetInputValue(&sealevel,gauss);
 		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
 		element->NodalFunctions(basis,gauss);
 
-		surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
-		base_under_water    = min(0.,bed);           // 0 if the bottom of the glacier is above water level
+		surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level
+		base_under_water    = min(0.,bed-sealevel);           // 0 if the bottom of the glacier is above water level
 		water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
 		ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
@@ -2584,5 +2588,5 @@
 	/*Intermediaries*/
 	int         dim;
-	IssmDouble  Jdet,surface,z,water_pressure,ice_pressure;
+	IssmDouble  Jdet,surface,sealevel,z,water_pressure,ice_pressure;
 	IssmDouble  surface_under_water,base_under_water,pressure;
 	IssmDouble* xyz_list       = NULL;
@@ -2604,4 +2608,5 @@
 	/*Retrieve all inputs and parameters*/
 	Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
+	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
 	IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
 	IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
@@ -2622,4 +2627,5 @@
 		gauss->GaussPoint(ig);
 		surface_input->GetInputValue(&surface,gauss);
+		sealevel_input->GetInputValue(&sealevel,gauss);
 		if(dim==3) z=element->GetZcoord(xyz_list,gauss);
 		else       z=element->GetYcoord(xyz_list,gauss);
@@ -2627,5 +2633,5 @@
 		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
 
-		water_pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
+		water_pressure = rho_water*gravity*min(0.,z-sealevel);//0 if the gaussian point is above water level
 		ice_pressure   = rho_ice*gravity*(surface-z);
 		pressure       = ice_pressure + water_pressure;
@@ -3811,5 +3817,5 @@
 	/*Intermediaries*/
 	int         i,dim;
-	IssmDouble  Jdet,pressure,surface,z;
+	IssmDouble  Jdet,pressure,surface,sealevel,z;
 	IssmDouble	normal[3];
 	IssmDouble *xyz_list       = NULL;
@@ -3843,4 +3849,5 @@
 	element->NormalSection(&normal[0],xyz_list_front);
 	Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
+	Input* sealevel_input       = element->GetInput(SealevelEnum);       _assert_(sealevel_input);
 	IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
 	IssmDouble  gravity   = element->GetMaterialParameter(ConstantsGEnum);
@@ -3859,7 +3866,8 @@
 		element->NodalFunctionsVelocity(vbasis,gauss);
 		surface_input->GetInputValue(&surface,gauss);
+		sealevel_input->GetInputValue(&sealevel,gauss);
 		if(dim==3) z=element->GetZcoord(xyz_list,gauss);
 		else       z=element->GetYcoord(xyz_list,gauss);
-		pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
+		pressure = rho_water*gravity*min(0.,z-sealevel);//0 if the gaussian point is above water level
 
 		for (int i=0;i<vnumnodes;i++){
Index: issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 20155)
+++ issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 20156)
@@ -1653,4 +1653,5 @@
 	IssmDouble* b       = xNew<IssmDouble>(numvertices);
 	IssmDouble* r       = xNew<IssmDouble>(numvertices);
+	IssmDouble* sl      = xNew<IssmDouble>(numvertices);
 
 	/*Recover info at the vertices: */
@@ -1661,4 +1662,5 @@
 	GetInputListOnVertices(&b[0],BaseEnum);
 	GetInputListOnVertices(&r[0],BedEnum);
+	GetInputListOnVertices(&sl[0],SealevelEnum);
 	GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
 	rho_water   = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
@@ -1684,14 +1686,14 @@
 		/*Change only if AggressiveMigration or if the ice sheet is in contact with the ocean*/
 		else{ // phi>0
-			bed_hydro=-density*h[i];
+			bed_hydro=-density*h[i]+sl[i];
 			if (bed_hydro>r[i]){
 				/*Unground only if the element is connected to the ice shelf*/
 				if(migration_style==AggressiveMigrationEnum || migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
-					s[i]        = (1-density)*h[i];
-					b[i]        = -density*h[i];
+					s[i]        = (1-density)*h[i]+sl[i];
+					b[i]        = -density*h[i]+sl[i];
 				}
 				else if(migration_style==SoftMigrationEnum && phi_ungrounding[vertices[i]->Pid()]<0.){
-					s[i]        = (1-density)*h[i];
-					b[i]        = -density*h[i];
+					s[i]        = (1-density)*h[i]+sl[i];
+					b[i]        = -density*h[i]+sl[i];
 				}
 				else{
@@ -1705,10 +1707,10 @@
 	for(i=0;i<numvertices;i++){
 		if(migration_style==SoftMigrationEnum){
-			bed_hydro=-density*h[i];
+			bed_hydro=-density*h[i]+sl[i];
 			if(phi[i]<0. || bed_hydro<=r[i] || phi_ungrounding[vertices[i]->Pid()]<0.){
-				phi[i]=h[i]+r[i]/density;
+				phi[i]=h[i]+(r[i]-sl[i])/density;
 			}
 		}
-		else if(migration_style!=ContactEnum) phi[i]=h[i]+r[i]/density;
+		else if(migration_style!=ContactEnum) phi[i]=h[i]+(r[i]-sl[i])/density;
 		else{
 			/*do nothing*/
@@ -1727,4 +1729,5 @@
 	xDelete<IssmDouble>(b);
 	xDelete<IssmDouble>(s);
+	xDelete<IssmDouble>(sl);
 	xDelete<IssmDouble>(h);
 
Index: issm/trunk-jpl/src/c/classes/Loads/Friction.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 20155)
+++ issm/trunk-jpl/src/c/classes/Loads/Friction.cpp	(revision 20156)
@@ -67,5 +67,5 @@
 	IssmDouble  Neff;
 	IssmDouble  drag_coefficient;
-	IssmDouble  bed,thickness;
+	IssmDouble  bed,thickness,sealevel;
 	IssmDouble  alpha_complement;
 
@@ -75,4 +75,5 @@
 	element->GetInputValue(&thickness, gauss,ThicknessEnum);
 	element->GetInputValue(&bed, gauss,BaseEnum);
+	element->GetInputValue(&sealevel, gauss,SealevelEnum);
 	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
 	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
@@ -85,5 +86,5 @@
 
 	//From bed and thickness, compute effective pressure when drag is viscous:
-	Neff=gravity*(rho_ice*thickness+rho_water*bed);
+	Neff=gravity*(rho_ice*thickness+rho_water*(bed-sealevel));
 	if(Neff<0)Neff=0;
 
@@ -231,5 +232,5 @@
 	IssmDouble  drag_p, drag_q;
 	IssmDouble  Neff;
-	IssmDouble  thickness,base,bed,floatation_thickness;
+	IssmDouble  thickness,base,bed,floatation_thickness,sealevel;
 	IssmDouble  vx,vy,vz,vmag;
 	IssmDouble  drag_coefficient,drag_coefficient_coulomb;
@@ -241,4 +242,5 @@
 	element->GetInputValue(&thickness, gauss,ThicknessEnum);
 	element->GetInputValue(&base, gauss,BaseEnum);
+	element->GetInputValue(&sealevel, gauss,SealevelEnum);
 	element->GetInputValue(&bed, gauss,BedEnum);
 	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
@@ -253,5 +255,5 @@
 
 	//From base and thickness, compute effective pressure when drag is viscous:
-	Neff=gravity*(rho_ice*thickness+rho_water*base);
+	Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
 	if(Neff<0)Neff=0;
 
@@ -407,5 +409,5 @@
 	IssmDouble  drag_p, drag_q;
 	IssmDouble  Neff;
-	IssmDouble  thickness,bed;
+	IssmDouble  thickness,base,sealevel;
 	IssmDouble  vx,vy,vz,vmag;
 	IssmDouble  drag_coefficient;
@@ -416,5 +418,6 @@
 	element->GetInputValue(&drag_q,FrictionQEnum);
 	element->GetInputValue(&thickness, gauss,ThicknessEnum);
-	element->GetInputValue(&bed, gauss,BaseEnum);
+	element->GetInputValue(&base, gauss,BaseEnum);
+	element->GetInputValue(&sealevel, gauss,SealevelEnum);
 	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
 	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
@@ -426,6 +429,6 @@
 	s=1./drag_p;
 
-	//From bed and thickness, compute effective pressure when drag is viscous:
-	Neff=gravity*(rho_ice*thickness+rho_water*bed);
+	//From base and thickness, compute effective pressure when drag is viscous:
+	Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
 	if(Neff<0)Neff=0;
 
@@ -467,5 +470,5 @@
 	IssmDouble  drag_p, drag_q;
 	IssmDouble  Neff,F;
-	IssmDouble  thickness,bed;
+	IssmDouble  thickness,bed,sealevel;
 	IssmDouble  vx,vy,vz,vmag;
 	IssmDouble  drag_coefficient,water_layer;
@@ -478,4 +481,5 @@
 	element->GetInputValue(&thickness, gauss,ThicknessEnum);
 	element->GetInputValue(&bed, gauss,BaseEnum);
+	element->GetInputValue(&sealevel, gauss,SealevelEnum);
 	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
 	element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum);
@@ -490,5 +494,5 @@
 	//From bed and thickness, compute effective pressure when drag is viscous:
 	if(bed>0) bed=0;
-	if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*bed;
+	if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(bed-sealevel);
 	else if(water_layer>0) Neff=gravity*rho_ice*thickness*F;
 	else _error_("negative water layer thickness");
@@ -599,5 +603,5 @@
 	IssmDouble  Neff;
 	IssmDouble  drag_coefficient;
-	IssmDouble  bed,thickness,head;
+	IssmDouble  bed,thickness,head,sealevel;
 	IssmDouble  alpha2;
 
@@ -606,4 +610,5 @@
 	element->GetInputValue(&bed, gauss,BaseEnum);
 	element->GetInputValue(&head, gauss,HydrologyHeadEnum);
+	element->GetInputValue(&sealevel, gauss,SealevelEnum);
 	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
 	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
@@ -613,5 +618,5 @@
 	//From bed and thickness, compute effective pressure when drag is viscous:
 	pressure_ice   = rho_ice*gravity*thickness;
-	pressure_water = rho_water*gravity*(head-bed);
+	pressure_water = rho_water*gravity*(head-bed+sealevel);
 	Neff=pressure_ice-pressure_water;
 	if(Neff<0.) Neff=0.;
