Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23941)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23942)
@@ -324,2 +324,72 @@
 	return;
 }/*}}}*/
+
+/*GlaDS specifics*/
+void HydrologyGlaDSAnalysis::UpdateSheetThickness(FemModel* femmodel){/*{{{*/
+
+	for(int j=0;j<femmodel->elements->Size();j++){
+		Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
+		UpdateSheetThickness(element);
+	}
+
+}/*}}}*/
+void HydrologyGlaDSAnalysis::UpdateSheetThickness(Element* element){/*{{{*/
+
+	/*Skip if water or ice shelf element*/
+	if(element->IsFloating()) return;
+
+	/*Intermediaries */
+	IssmDouble  Jdet,w,v,vx,vy,ub,h_old,N;
+	IssmDouble  A,B,n;
+
+	/*Fetch number vertices for this element*/
+	int numvertices = element->GetNumberOfVertices();
+
+	/*Initialize new sheet thickness*/
+	IssmDouble* h_new = xNew<IssmDouble>(numvertices);
+
+	/*Retrieve all inputs and parameters*/
+	IssmDouble  h_r = element->FindParam(HydrologyBumpHeightEnum);
+	IssmDouble  l_r = element->FindParam(HydrologyBumpSpacingEnum);
+	IssmDouble  dt  = element->FindParam(TimesteppingTimeStepEnum);
+	Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
+	Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
+	Input*  N_input = element->GetInput(EffectivePressureEnum); _assert_(N_input);
+	Input*  h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
+	Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+	Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss();
+	for(int iv=0;iv<numvertices;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Get input values at gauss points*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		h_input->GetInputValue(&h_old,gauss);
+		B_input->GetInputValue(&B,gauss);
+		n_input->GetInputValue(&n,gauss);
+		N_input->GetInputValue(&N,gauss);
+
+		/*Get basal velocity*/
+		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*/
+		A=pow(B,-n);
+		v = 2./pow(n,n)*A*h_old*pow(fabs(N),n-1.)*N;
+
+		/*Get new sheet thickness*/
+		h_new[iv] += h_old + dt*(w-v);
+	}
+
+	element->AddInput(HydrologySheetThicknessEnum,h_new,P1Enum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(h_new);
+	delete gauss;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h	(revision 23941)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.h	(revision 23942)
@@ -30,4 +30,8 @@
 		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void           UpdateConstraints(FemModel* femmodel);
+
+		/*Specific to GlaDS*/
+		void UpdateSheetThickness(FemModel* femmodel);
+		void UpdateSheetThickness(Element*  element);
 };
 #endif