Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23963)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23964)
@@ -347,5 +347,5 @@
 
 	/*Intermediaries */
-	IssmDouble  Jdet,v1,qc,fFactor;
+	IssmDouble  Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
 	IssmDouble  A,B,n,phi_old,phi,phi_0,dPw;
 	IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
@@ -367,5 +367,4 @@
 	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
 	IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
-	IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
 	IssmDouble g         = element->FindParam(ConstantsGEnum);
 	IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
@@ -430,9 +429,23 @@
 		}
 
+		/*Compute Afactor and Bfactor*/
+		Afactor = C_W*c_t*rho_water;
+		Bfactor = 1./L * (1./rho_ice - 1./rho_water);
+		if(dphids>0){
+			Xifactor = + Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
+		}
+		else{
+			Xifactor = - Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
+		}
+
 		/*Diffusive term*/
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
+				/*GlaDSCoupledSolver.F90 line 1659*/
 				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
-							+Kc*dbasisds[i]*dbasisds[j] /*GlaDSCoupledSolver.F90 line 1659*/
+							+Kc*dbasisds[i]*dbasisds[j]                               /*Diffusion term*/
+							- Afactor * Bfactor* Kc * dPw * basis[i] * dbasisds[j]    /*First part of Pi*/
+							+ Afactor * fFactor * Bfactor * basis[i] * dbasisds[j]    /*Second part of Pi*/
+							+ Xifactor* basis[i] * dbasisds[j]                        /*Xi term*/
 							);
 			}
@@ -447,8 +460,4 @@
 			}
 		}
-
-		/*Transient term if dt>0*/
-		if(dt>0.){
-		}
 	}
 
@@ -468,8 +477,9 @@
 
 	/*Intermediaries */
-	IssmDouble  Jdet,v2;
-	IssmDouble  A,B,n,phi_old,phi,phi_0;
-	IssmDouble  H,h,b;
+	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;
 	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  xyz_list_tria[3][3];
 	const int   numnodes = NUMNODES;
 
@@ -480,9 +490,15 @@
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+	GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
+
 	IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
 	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
 	IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
-	IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
+	IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
+	IssmDouble ks        = element->FindParam(HydrologySheetConductivityEnum);
 	IssmDouble g         = element->FindParam(ConstantsGEnum);
+	IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
+	IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
+
 	Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
 	Input* H_input      = element->GetInput(ThicknessEnum); _assert_(H_input);
@@ -493,4 +509,10 @@
 	Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
+	/*Get tangent vector*/
+	IssmDouble tx = xyz_list_tria[index2][0] - xyz_list_tria[index1][0];
+	IssmDouble ty = xyz_list_tria[index2][1] - xyz_list_tria[index1][1];
+	tx = tx/sqrt(tx*tx+ty*ty);
+	ty = ty/sqrt(tx*tx+ty*ty);
+
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=new GaussTria(index1,index2,2);
@@ -502,4 +524,5 @@
 
 		/*Get input values at gauss points*/
+		b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
 		h_input->GetInputValue(&h,gauss);
 		B_input->GetInputValue(&B,gauss);
@@ -509,4 +532,27 @@
 		H_input->GetInputValue(&H,gauss);
 
+		/*Get values for a few potentials*/
+		dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
+
+		/*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.);
+
+		/*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
+		qc = - Ks * dphids;
+
+		/*d(phi - phi_m)/ds*/
+		dPw = dphids - dphimds;
+
+		/*Compute f factor*/
+		fFactor = 0.;
+		if(this->S>0. || qc*dPw>0.){
+			fFactor = lc * qc;
+		}
+
+		/*Compute Afactor and Bfactor*/
+		Afactor = C_W*c_t*rho_water;
+		Bfactor = 1./L * (1./rho_ice - 1./rho_water);
+
 		/*Compute closing rate*/
 		/*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/
@@ -515,10 +561,7 @@
 		v2 = 2./pow(n,n)*A*this->S*(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*(-v2)*basis[i];
-
-		/*Transient term if dt>0*/
-		if(dt>0.){
-			//phiold_input->GetInputValue(&phi_old,gauss);
-			//for(int i=0;i<numnodes;i++) pe->values[i] += gauss->weight*Jdet*e_v/(rho_water*g*dt)*phi_old*basis[i];
+		for(int i=0;i<numnodes;i++){
+			pe->values[i]+= - Jdet*gauss->weight*(-v2)*basis[i];
+			pe->values[i]+= + Jdet*gauss->weight*Afactor*Bfactor*fFactor*dphimds*basis[i];
 		}
 	}