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++){
