Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23960)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23961)
@@ -169,4 +169,5 @@
 	/*Intermediaries */
 	IssmDouble  Jdet,dphi[3],h,k;
+	IssmDouble  A,B,n,phi_old,phi,phi_0,H,b,v1;
 	IssmDouble* xyz_list = NULL;
 
@@ -189,4 +190,5 @@
 	IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
 	IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
 	IssmDouble g         = element->FindParam(ConstantsGEnum);
 	IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
@@ -194,4 +196,8 @@
 	Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_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* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+	Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -205,6 +211,11 @@
 
 		phi_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+		phi_input->GetInputValue(&phi,gauss);
 		h_input->GetInputValue(&h,gauss);
 		k_input->GetInputValue(&k,gauss);
+		B_input->GetInputValue(&B,gauss);
+		n_input->GetInputValue(&n,gauss);
+		b_input->GetInputValue(&b,gauss);
+		H_input->GetInputValue(&H,gauss);
 
 		/*Get norm of gradient of hydraulic potential and make sure it is >0*/
@@ -220,4 +231,14 @@
 							coeff*dbasis[0*numnodes+i]*dbasis[0*numnodes+j]
 							+ coeff*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
+			}
+		}
+
+		/*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
+		phi_0 = rho_water*g*b + rho_ice*g*H;
+		A=pow(B,-n);
+		v1 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*( - n));
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j];
 			}
 		}
@@ -248,7 +269,8 @@
 
 	/*Intermediaries */
-	IssmDouble  Jdet,w,v,vx,vy,ub,h,N,h_r;
+	IssmDouble  Jdet,w,v2,vx,vy,ub,h,h_r;
 	IssmDouble  G,m,frictionheat,alpha2;
-	IssmDouble  A,B,n,phi_old;
+	IssmDouble  A,B,n,phi_old,phi,phi_0;
+	IssmDouble  H,b;
 	IssmDouble* xyz_list = NULL;
 
@@ -272,10 +294,12 @@
 	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* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
 	Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
 	Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
 	Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
+	Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
 	/*Build friction element, needed later: */
@@ -297,6 +321,8 @@
 		B_input->GetInputValue(&B,gauss);
 		n_input->GetInputValue(&n,gauss);
-		N_input->GetInputValue(&N,gauss);
 		hr_input->GetInputValue(&h_r,gauss);
+		phi_input->GetInputValue(&phi,gauss);
+		b_input->GetInputValue(&b,gauss);
+		H_input->GetInputValue(&H,gauss);
 
 		/*Get basal velocity*/
@@ -315,8 +341,10 @@
 
 		/*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);
-		v = 2./pow(n,n)*A*h*pow(fabs(N),n-1.)*N;
-
-		for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v-m)*basis[i];
+		v2 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v2-m)*basis[i];
 
 		/*Transient term if dt>0*/
@@ -413,4 +441,5 @@
 		/*Get new sheet thickness*/
 		h_new[iv] = h_old + dt*(w-v);
+		if(h_new[iv]<1.e-12) h_new[iv] = 1.e-12;
 	}