Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23962)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23963)
@@ -14,6 +14,15 @@
 #include "../classes.h"
 /*}}}*/	
+
+/*Macros*/
 #define NUMNODES    2
 #define NUMVERTICES 2
+
+#define C_W         4.22e3   /*specific heat capacity of water (J/kg/K)*/
+#define ALPHA_C     5./4.
+#define BETA_C      3./2.
+/*Make sure these are the same as in HydrologyGlaDSAnalysis::CreateKMatrix*/
+#define ALPHA_S     5./4.
+#define BETA_S      3./2.
 
 /*Channel constructors and destructor*/
@@ -338,7 +347,7 @@
 
 	/*Intermediaries */
-	IssmDouble  Jdet,v1;
-	IssmDouble  A,B,n,phi_old,phi,phi_0;
-	IssmDouble  H,h,b;
+	IssmDouble  Jdet,v1,qc,fFactor;
+	IssmDouble  A,B,n,phi_old,phi,phi_0,dPw;
+	IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_tria[3][3];
@@ -354,4 +363,5 @@
 	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);
@@ -359,4 +369,9 @@
 	IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
 	IssmDouble g         = element->FindParam(ConstantsGEnum);
+	IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
+	IssmDouble ks        = element->FindParam(HydrologySheetConductivityEnum);
+	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);
@@ -381,15 +396,48 @@
 		tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
 		tria->GetSegmentNodalFunctionsDerivatives(&dbasisdx[0],&xyz_list_tria[0][0],gauss,index1,index2,tria->FiniteElement());
+		dbasisds[0] = dbasisdx[0*2+0]*tx + dbasisdx[0*2+1]*ty;
+		dbasisds[1] = dbasisdx[1*2+0]*tx + dbasisdx[1*2+1]*ty;
 
 		/*Get input values at gauss points*/
+		phi_input->GetInputValue(&phi,gauss);
+		phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
 		h_input->GetInputValue(&h,gauss);
 		B_input->GetInputValue(&B,gauss);
 		n_input->GetInputValue(&n,gauss);
-		phi_input->GetInputValue(&phi,gauss);
 		b_input->GetInputValue(&b,gauss);
+		b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
 		H_input->GetInputValue(&H,gauss);
 
+		/*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);
+
+		/*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;
+		}
+
+		/*Diffusive term*/
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+							+Kc*dbasisds[i]*dbasisds[j] /*GlaDSCoupledSolver.F90 line 1659*/
+							);
+			}
+		}
+
 		/*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*S*(pow(fabs(phi_0 - phi),n-1.)*( - n));