Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 24062)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 24063)
@@ -33,4 +33,5 @@
 		CreateFaces(iomodel);
 		for(int i=0;i<iomodel->numberoffaces;i++){
+
 			/*Get left and right elements*/
 			int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
@@ -303,4 +304,5 @@
 
 	/*Intermediaries */
+	bool        meltflag;
 	IssmDouble  Jdet,w,v2,vx,vy,ub,h,h_r;
 	IssmDouble  G,m,frictionheat,alpha2;
@@ -317,7 +319,6 @@
 
 	/*Retrieve all inputs and parameters*/
-	bool meltflag;
+	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&meltflag,HydrologyMeltFlagEnum);
-	element->GetVerticesCoordinates(&xyz_list);
 	IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
 	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
@@ -362,5 +363,4 @@
 		b_input->GetInputValue(&b,gauss);
 		H_input->GetInputValue(&H,gauss);
-		m_input->GetInputValue(&m,gauss);
 
 		/*Get basal velocity*/
@@ -375,11 +375,13 @@
 		frictionheat=alpha2*ub*ub;
 
-		/*Compute melt*/
+		/*Compute melt (if necessary)*/
 		if(!meltflag){
 			m = (G + frictionheat)/(rho_ice*L);
 		}
+		else{
+			m_input->GetInputValue(&m,gauss);
+		}
 
 		/*Compute closing rate*/
-		/*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/
 		phi_0 = rho_water*g*b + rho_ice*g*H;
 		A=pow(B,-n);
@@ -433,6 +435,6 @@
 
 	/*Intermediaries */
-	IssmDouble  Jdet,vx,vy,ub,h_old,N,h_r;
-	IssmDouble  A,B,n;
+	IssmDouble  Jdet,vx,vy,ub,h_old,N,h_r,H,b;
+	IssmDouble  A,B,n,phi,phi_0;
 	IssmDouble  alpha,beta;
 
@@ -444,13 +446,18 @@
 
 	/*Retrieve all inputs and parameters*/
-	IssmDouble  dt  = element->FindParam(TimesteppingTimeStepEnum);
-	IssmDouble  l_r = element->FindParam(HydrologyCavitySpacingEnum);
+	IssmDouble  dt       = element->FindParam(TimesteppingTimeStepEnum);
+	IssmDouble  l_r      = element->FindParam(HydrologyCavitySpacingEnum);
+	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
+	IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+	IssmDouble g         = element->FindParam(ConstantsGEnum);
 	Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
 	Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
 	Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
-	Input*  N_input = element->GetInput(EffectivePressureEnum); _assert_(N_input);
-	Input*  h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
+	Input* H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
+	Input* b_input  = element->GetInput(BedEnum); _assert_(b_input);
+	Input* h_input  = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
 	Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
 	Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+	Input* phi_input = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -460,4 +467,5 @@
 
 		/*Get input values at gauss points*/
+		phi_input->GetInputValue(&phi,gauss);
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
@@ -465,6 +473,11 @@
 		B_input->GetInputValue(&B,gauss);
 		n_input->GetInputValue(&n,gauss);
-		N_input->GetInputValue(&N,gauss);
 		hr_input->GetInputValue(&h_r,gauss);
+		b_input->GetInputValue(&b,gauss);
+		H_input->GetInputValue(&H,gauss);
+
+		/*Get values for a few potentials*/
+		phi_0 = rho_water*g*b + rho_ice*g*H;
+		N = phi_0 - phi;
 
 		/*Get basal velocity*/
@@ -485,5 +498,5 @@
 
 		/*Get new sheet thickness*/
-		h_new[iv] = ODE1(alpha,beta,h_old,dt,0);
+		h_new[iv] = ODE1(alpha,beta,h_old,dt,1);
 
 		/*Make sure it is positive*/