Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 27288)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 27289)
@@ -518,4 +518,5 @@
 	IssmDouble  A,B,n,phi,phi_0;
 	IssmDouble  alpha,beta;
+	IssmDouble  oceanLS,iceLS;
 
 	/*Fetch number vertices for this element*/
@@ -548,4 +549,7 @@
 	Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
 	Input* phi_input = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
+	Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input);
+	Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input);
+
 
 	/* Start  looping on the number of gaussian points: */
@@ -564,4 +568,12 @@
 		b_input->GetInputValue(&b,gauss);
 		H_input->GetInputValue(&H,gauss);
+		oceanLS_input->GetInputValue(&oceanLS,gauss);
+		iceLS_input->GetInputValue(&iceLS,gauss);
+
+		/*Set sheet thickness to zero if floating or no ice*/
+		if(oceanLS<0. || iceLS>0.){
+			h_new[iv] = 0.;
+		}
+		else{
 
 		/*Get values for a few potentials*/
@@ -591,9 +603,10 @@
 		if(h_new[iv]<AEPS) h_new[iv] = AEPS;
 	}
-	
+	}
+
 	/*Force floating ice to have zero sheet thickness*/
-	if(element->IsAllFloating() || !element->IsIceInElement()){
+	/*if(!element->IsAllGrounded() || !element->IsIceInElement()){
 		for(int iv=0;iv<numvertices;iv++) h_new[iv] = 0.;
-	}
+	}*/
 
 	element->AddInput(HydrologySheetThicknessEnum,h_new,P1Enum);
@@ -616,4 +629,5 @@
 	IssmDouble phi_0, phi_m, p_i;
 	IssmDouble H,b,phi;
+	IssmDouble oceanLS,iceLS;
 
 	int numnodes = element->GetNumberOfNodes();
@@ -627,4 +641,6 @@
 	Input* b_input   = element->GetInput(BedEnum); _assert_(b_input);
 	Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input);
+	Input* oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input);
+	Input* iceLS_input = element->GetInput(MaskIceLevelsetEnum); _assert_(iceLS_input);
 
 	/*Set to 0 if inactive element*/
@@ -645,4 +661,6 @@
 		b_input->GetInputValue(&b,gauss);
 		phi_input->GetInputValue(&phi,gauss);
+		oceanLS_input->GetInputValue(&oceanLS,gauss);
+		iceLS_input->GetInputValue(&iceLS,gauss);
 
 		/*Elevation potential*/
@@ -657,4 +675,8 @@
 		/*Calculate effective pressure*/
 		N[iv] = phi_0 - phi;
+
+		/*Make sure that all floating ice and ice free areas have zero effective pressure*/
+		if(oceanLS<0.0) N[iv] = 0.0;
+		if(iceLS>0.0) N[iv] = 0.0;
 
 		if(xIsNan<IssmDouble>(N[iv])) _error_("NaN found in solution vector");