Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23969)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 23970)
@@ -60,4 +60,18 @@
 	}
 	iomodel->DeleteData(1,"md.mesh.vertexonbase");
+
+	/*Deal with Neumann BC*/
+	int M,N;
+	int *segments = NULL;
+	iomodel->FetchData(&segments,&M,&N,"md.mesh.segments");
+
+	/*Check that the size seem right*/
+	_assert_(N==3); _assert_(M>=3);
+	for(int i=0;i<M;i++){
+		if(iomodel->my_elements[segments[i*3+2]-1]){
+			loads->AddObject(new Neumannflux(i+1,i,iomodel,segments));
+		}
+	}
+
 }/*}}}*/
 void HydrologyGlaDSAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
@@ -108,4 +122,5 @@
 	iomodel->FetchDataToInput(elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
 	iomodel->FetchDataToInput(elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum);
+	iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum,0.);
 	iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum);
 	iomodel->FetchDataToInput(elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum);
Index: /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp	(revision 23969)
+++ /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp	(revision 23970)
@@ -187,4 +187,7 @@
 			/*Nothing!*/
 			break;
+		case HydrologyGlaDSAnalysisEnum:
+			/*Nothing!*/
+			break;
 		default:
 			_error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
@@ -209,4 +212,7 @@
 		case HydrologyShaktiAnalysisEnum:
 			pe=CreatePVectorHydrologyShakti();
+			break;
+		case HydrologyGlaDSAnalysisEnum:
+			pe=CreatePVectorHydrologyGlaDS();
 			break;
 		default:
@@ -379,2 +385,46 @@
 }
 /*}}}*/
+ElementVector* Neumannflux::CreatePVectorHydrologyGlaDS(void){/*{{{*/
+
+	/* constants*/
+	const int numdof=2;
+
+	/* Intermediaries*/
+	IssmDouble Jdet,flux;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble basis[numdof];
+
+	/*Initialize Load Vector and return if necessary*/
+	Tria*  tria=(Tria*)element;
+	_assert_(tria->FiniteElement()==P1Enum); 
+	if(!tria->IsIceInElement() || tria->IsFloating()) return NULL;
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	Input* flux_input = tria->inputs->GetInput(HydrologyNeumannfluxEnum);  _assert_(flux_input); 
+
+	/*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
+	int index1=tria->GetNodeIndex(nodes[0]);
+	int index2=tria->GetNodeIndex(nodes[1]);
+
+	/* Start  looping on the number of gaussian points: */
+	GaussTria* gauss=new GaussTria(index1,index2,2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
+		tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
+		flux_input->GetInputValue(&flux,gauss);
+
+		for(int i=0;i<numdof;i++) pe->values[i] += gauss->weight*Jdet*flux*basis[i];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	return pe;
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.h	(revision 23969)
+++ /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.h	(revision 23970)
@@ -72,4 +72,5 @@
 		/*Neumannflux management:{{{*/
 		ElementVector* CreatePVectorHydrologyShakti(void);
+		ElementVector* CreatePVectorHydrologyGlaDS(void);
 		/*}}}*/