Index: ../trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 24063)
+++ ../trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 24064)
@@ -54,6 +54,13 @@
 	int e1 = iomodel->faces[4*index+2];
 	int e2 = iomodel->faces[4*index+3];
 
+	if(e2==-1){
+		this->boundary = true;
+	}
+	else{
+		this->boundary = false;
+	}
+
 	/*Set Element hook (4th column may be -1 for boundary edges)*/
 	this->helement  = new Hook(&e1,1);
 
@@ -544,9 +551,8 @@
 		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(Ngrad,BETA_C-2.);
-		IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
+		/*Compute the effective conductivity Ks = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
+		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*/
 		qc = - Ks * dphids;
@@ -584,7 +590,7 @@
 
 	/*Initialize Element matrix and return if necessary*/
 	Tria*  tria=(Tria*)element;
-	if(!tria->IsIceInElement()){
+	if(!tria->IsIceInElement() || this->boundary){
 		this->S = 0.;
 		return;
 	}
@@ -673,8 +679,11 @@
 	IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
 
 	/*Solve ODE*/
-	this->S = ODE1(alpha,beta,this->S,dt,0);
+	this->S = ODE1(alpha,beta,this->S,dt,2);
 
+	/*Make sure Area > 0*/
+	if(this->S<0.) this->S = 0.;
+
 	/*Clean up and return*/
 	delete gauss;
 }
Index: ../trunk-jpl/src/c/classes/Loads/Channel.h
===================================================================
--- ../trunk-jpl/src/c/classes/Loads/Channel.h	(revision 24063)
+++ ../trunk-jpl/src/c/classes/Loads/Channel.h	(revision 24064)
@@ -19,6 +19,7 @@
 
 	private: 
 		IssmDouble S;
+		bool       boundary;
 
 	public: 
 		int sid;