Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 27292)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 27293)
@@ -444,4 +444,81 @@
 void           HydrologyGlaDSAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 	element->InputUpdateFromSolutionOneDof(solution,HydraulicPotentialEnum);
+
+	/*Compute Hydrology Vx and Vy for time stepping purposes (These inputs do not affect GlaDS)*/
+
+	/*Intermediaries*/
+   IssmDouble  dphi[3],h,k,phi;
+	IssmDouble  oceanLS,iceLS;
+	IssmDouble* xyz_list = NULL;
+
+	/*Hard coded coefficients*/
+	const IssmPDouble alpha = 5./4.;
+	const IssmPDouble beta  = 3./2.;
+
+	/*Fetch number vertices for this element*/
+	int numvertices = element->GetNumberOfVertices();
+
+	/*Initialize new sheet thickness*/
+	IssmDouble* vx = xNew<IssmDouble>(numvertices);
+	IssmDouble* vy = xNew<IssmDouble>(numvertices);
+
+	/*Set to 0 if inactive element*/
+	if(element->IsAllFloating() || !element->IsIceInElement()){
+		for(int iv=0;iv<numvertices;iv++) vx[iv] = 0.;
+		for(int iv=0;iv<numvertices;iv++) vy[iv] = 0.;
+		element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
+		element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
+		xDelete<IssmDouble>(vx);
+		xDelete<IssmDouble>(vy);
+		return;
+	}
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input *k_input       = element->GetInput(HydrologySheetConductivityEnum); _assert_(k_input);
+	Input *phi_input     = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
+	Input *h_input       = element->GetInput(HydrologySheetThicknessEnum);    _assert_(h_input);
+	Input *oceanLS_input = element->GetInput(MaskOceanLevelsetEnum);          _assert_(oceanLS_input);
+	Input *iceLS_input   = element->GetInput(MaskIceLevelsetEnum);            _assert_(iceLS_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*/
+      phi_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+      phi_input->GetInputValue(&phi,gauss);
+      h_input->GetInputValue(&h,gauss);
+      k_input->GetInputValue(&k,gauss);
+		oceanLS_input->GetInputValue(&oceanLS,gauss);
+		iceLS_input->GetInputValue(&iceLS,gauss);
+
+		/*Set sheet thickness to zero if floating or no ice*/
+		if(oceanLS<0. || iceLS>0.){
+			vx[iv] = 0.;
+         vy[iv] = 0.;
+		}
+		else{
+
+         /*Get norm of gradient of hydraulic potential and make sure it is >0*/
+         IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
+         if(normgradphi < AEPS) normgradphi = AEPS;
+
+         IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h); // divide by h to get speed instead of discharge
+
+			vx[iv] = -coeff*dphi[0];
+			vy[iv] = -coeff*dphi[1];
+		}
+	}
+
+	element->AddInput(HydrologyWaterVxEnum,vx,P1DGEnum);
+	element->AddInput(HydrologyWaterVyEnum,vy,P1DGEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(vy);
+	delete gauss;
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/