Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23976)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23977)
@@ -381,5 +381,4 @@
 	Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
 	Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
-	Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
 	Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
@@ -402,6 +401,7 @@
 
 		/*Get input values at gauss points*/
+		phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
+		b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
 		phi_input->GetInputValue(&phi,gauss);
-		phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
 		h_input->GetInputValue(&h,gauss);
 		ks_input->GetInputValue(&ks,gauss);
@@ -409,5 +409,4 @@
 		n_input->GetInputValue(&n,gauss);
 		b_input->GetInputValue(&b,gauss);
-		b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
 		H_input->GetInputValue(&H,gauss);
 
@@ -484,6 +483,6 @@
 	/*Intermediaries */
 	IssmDouble  Jdet,v2,Afactor,Bfactor,fFactor;
-	IssmDouble  A,B,n,phi_old,phi,phi_0,dphimds;
-	IssmDouble  H,h,b,db[2],dphids,qc,dPw,ks;
+	IssmDouble  A,B,n,phi_old,phi,phi_0,dphimds,dphi[2];
+	IssmDouble  H,h,b,db[2],dphids,qc,dPw,ks,Ngrad;
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_tria[3][3];
@@ -512,5 +511,4 @@
 	Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
 	Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
-	Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
 	Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
@@ -531,4 +529,5 @@
 		/*Get input values at gauss points*/
 		b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
+		phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
 		h_input->GetInputValue(&h,gauss);
 		ks_input->GetInputValue(&ks,gauss);
@@ -540,9 +539,13 @@
 
 		/*Get values for a few potentials*/
+		phi_0   = rho_water*g*b + rho_ice*g*H;
+		dphids  = dphi[0]*tx + dphi[1]*ty;
 		dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
+		Ngrad   = fabs(dphids);
+		if(Ngrad<1.e-12) Ngrad = 1.e-12;
 
 		/*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
-		IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(fabs(dphids),BETA_C-2.);
-		IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(fabs(dphids),BETA_S-2.);
+		IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
+		IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
 
 		/*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
@@ -564,5 +567,4 @@
 		/*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);
 		v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
@@ -616,5 +618,4 @@
 	Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
 	Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
-	Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
 	Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
@@ -676,5 +677,5 @@
 	/*Compute new S based on Forward Euler (explicit)*/
 	this->S = this->S + dt*( (Xi - Pi)/(rho_ice*L) - vc);
-	if(this->S<0); this->S = 0.;
+	if(this->S<0.); this->S = 0.;
 
 	/*Clean up and return*/