Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 24056)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 24057)
@@ -425,6 +425,7 @@
 
 	/*Intermediaries */
-	IssmDouble  Jdet,w,v,vx,vy,ub,h_old,N,h_r;
+	IssmDouble  Jdet,vx,vy,ub,h_old,N,h_r;
 	IssmDouble  A,B,n;
+	IssmDouble  alpha,beta;
 
 	/*Fetch number vertices for this element*/
@@ -462,14 +463,21 @@
 		ub = sqrt(vx*vx + vy*vy);
 
-		/*Compute cavity opening w*/
-		w  = 0.;
-		if(h_old<h_r) w = ub*(h_r-h_old)/l_r;
-
-		/*Compute closing rate*/
+		/*Get A from B and n*/
 		A=pow(B,-n);
-		v = 2./pow(n,n)*A*h_old*pow(fabs(N),n-1.)*N;
+
+		/*Define alpha and beta*/
+		if(h_old<h_r){
+			alpha = -ub/l_r - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+			beta  = ub*h_r/l_r;
+		}
+		else{
+			alpha = - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+			beta  = 0.;
+		}
 
 		/*Get new sheet thickness*/
-		h_new[iv] = h_old + dt*(w-v);
+		h_new[iv] = ODE1(alpha,beta,h_old,dt,0);
+
+		/*Make sure it is positive*/
 		if(h_new[iv]<1.e-12) h_new[iv] = 1.e-12;
 	}
Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 24056)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 24057)
@@ -648,13 +648,9 @@
 	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.);
-
-	/*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
-	IssmDouble qc = - Ks * dphids;
-
 	/*d(phi - phi_m)/ds*/
 	IssmDouble dPw = dphids - dphimds;
+
+	/*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
+	IssmDouble qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids;
 
 	/*Compute f factor*/
@@ -664,18 +660,19 @@
 	}
 
-	/*Compute total discharge*/
-	IssmDouble Q = -Kc*dphids;
-
-	/*Compute Pi and Xi*/
-	IssmDouble Pi = -C_W*c_t*rho_water*(Q+fFactor)*dPw;
-	IssmDouble Xi = fabs(Q*dphids) + fabs(lc * qc * dphids);
-
-	/*Compute closing rate*/
 	A=pow(B,-n);
-	IssmDouble vc = 2./pow(n,n)*A*this->S*pow(fabs(phi_0 - phi),n-1.)*(phi_0 - phi);
-
-	/*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.;
+
+	IssmDouble C = C_W*c_t*rho_water;
+	IssmDouble Qprime = -kc * pow(Ngrad,BETA_C-2.)*dphids;
+	IssmDouble N = phi_0 - phi;
+
+	IssmDouble alpha = 1./(rho_ice*L)*(
+				fabs(Qprime*pow(this->S,ALPHA_C-1.)*dphids)
+				+ C*Qprime*pow(this->S,ALPHA_C-1.)*dPw
+				) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+
+	IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
+
+	/*Solve ODE*/
+	this->S = ODE1(alpha,beta,this->S,dt,0);
 
 	/*Clean up and return*/