Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23967)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23968)
@@ -42,7 +42,22 @@
 			loads->AddObject(new Channel(i+1,i,i,iomodel));
 		}
-
-		/*Free data: */
-	}
+	}
+
+	/*Create discrete loads for Moulins*/
+	CreateSingleNodeToElementConnectivity(iomodel);
+	for(int i=0;i<iomodel->numberofvertices;i++){
+		if (iomodel->domaintype!=Domain3DEnum){
+			/*keep only this partition's nodes:*/
+			if(iomodel->my_vertices[i]){
+				loads->AddObject(new Moulin(i+1,i,iomodel));
+			}
+		}
+		else if(reCast<int>(iomodel->Data("md.mesh.vertexonbase")[i])){
+			if(iomodel->my_vertices[i]){
+				loads->AddObject(new Moulin(i+1,i,iomodel));
+			}	
+		}
+	}
+	iomodel->DeleteData(1,"md.mesh.vertexonbase");
 }/*}}}*/
 void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23967)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 23968)
@@ -585,16 +585,8 @@
 
 	/*Intermediaries */
-	IssmDouble  Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
-	IssmDouble  A,B,n,phi_old,phi,phi_0,dPw;
+	IssmDouble  A,B,n,phi,phi_0;
 	IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_tria[3][3];
-	const int   numnodes = NUMNODES;
-
-	/*Initialize Element vector and other vectors*/
-	ElementMatrix* Ke=new ElementMatrix(this->nodes,NUMNODES,this->parameters);
-	IssmDouble     basis[NUMNODES];
-	IssmDouble     dbasisdx[2*NUMNODES];
-	IssmDouble     dbasisds[NUMNODES];
 
 	/*Retrieve all inputs and parameters*/
@@ -610,4 +602,5 @@
 	IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
 	IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
+	IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
 
 	Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
@@ -628,10 +621,4 @@
 	GaussTria* gauss=new GaussTria();
 	gauss->GaussEdgeCenter(index1,index2);
-
-	tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-	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*/
@@ -655,39 +642,29 @@
 
 	/*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
-	qc = - Ks * dphids;
+	IssmDouble qc = - Ks * dphids;
 
 	/*d(phi - phi_m)/ds*/
-	dPw = dphids - dphimds;
+	IssmDouble dPw = dphids - dphimds;
 
 	/*Compute f factor*/
-	fFactor = 0.;
+	IssmDouble 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);
-	if(dphids>0){
-		Xifactor = + Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
-	}
-	else{
-		Xifactor = - Bfactor * (fabs(-Kc*dphids) + fabs(lc*qc));
-	}
-
-	_error_("STOP");
-
-	/*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]                               /*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*/
-						);
-		}
-	}
+	/*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.;
 
 	/*Clean up and return*/