Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 17421)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 17422)
@@ -86,4 +86,5 @@
 	iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
 	iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
+	iomodel->FetchDataToInput(elements,HydrologydcBasalMoulinInputEnum);
 	iomodel->FetchDataToInput(elements,SedimentHeadEnum);
 	iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum);
@@ -278,5 +279,6 @@
 	/*Intermediaries */
 	bool       active_element,isefficientlayer;
-	IssmDouble dt,scalar,water_head;
+	IssmDouble dt,scalar;
+	IssmDouble moulin_load,water_head;
 	IssmDouble water_load,transfer;
 	IssmDouble Jdet;
@@ -304,6 +306,7 @@
 	Input* thickness_input   = basalelement->GetInput(ThicknessEnum);
 	Input* bed_input         = basalelement->GetInput(BedEnum);
-	Input* water_input       = basalelement->GetInput(BasalforcingsMeltingRateEnum); _assert_(water_input);
-	if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);          _assert_(old_wh_input);}
+	Input* water_input       = basalelement->GetInput(BasalforcingsMeltingRateEnum);    _assert_(water_input);
+	Input* moulin_input      = basalelement->GetInput(HydrologydcBasalMoulinInputEnum); _assert_(moulin_input);
+	if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);             _assert_(old_wh_input);}
 
 	IssmDouble sediment_storing    = SedimentStoring(basalelement);
@@ -324,6 +327,8 @@
 		/*Loading term*/
 		water_input->GetInputValue(&water_load,gauss);
+		moulin_input->GetInputValue(&moulin_load,gauss);
 	
 		scalar = Jdet*gauss->weight*(water_load);
+		scalar = scalar + Jdet* moulin_load;
 		if(dt!=0.) scalar = scalar*dt;
 		for(int i=0;i<numnodes;i++){
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 17421)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 17422)
@@ -126,4 +126,5 @@
 	HydrologydcPenaltyFactorEnum,
 	HydrologydcPenaltyLockEnum,
+	HydrologydcBasalMoulinInputEnum,
 	HydrologyLayerEnum,
 	HydrologySedimentEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 17421)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 17422)
@@ -134,4 +134,5 @@
 		case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor";
 		case HydrologydcPenaltyLockEnum : return "HydrologydcPenaltyLock";
+		case HydrologydcBasalMoulinInputEnum : return "HydrologydcBasalMoulinInput";
 		case HydrologyLayerEnum : return "HydrologyLayer";
 		case HydrologySedimentEnum : return "HydrologySediment";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 17421)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 17422)
@@ -134,11 +134,12 @@
 	      else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
 	      else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum;
+	      else if (strcmp(name,"HydrologydcBasalMoulinInput")==0) return HydrologydcBasalMoulinInputEnum;
 	      else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum;
 	      else if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum;
-	      else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
          else stage=2;
    }
    if(stage==2){
-	      if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
+	      if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
+	      else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
 	      else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
 	      else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
@@ -259,9 +260,9 @@
 	      else if (strcmp(name,"MaxIterationConvergenceFlag")==0) return MaxIterationConvergenceFlagEnum;
 	      else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
-	      else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
+	      if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
+	      else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
 	      else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
 	      else if (strcmp(name,"Surface")==0) return SurfaceEnum;
@@ -382,9 +383,9 @@
 	      else if (strcmp(name,"Nodes")==0) return NodesEnum;
 	      else if (strcmp(name,"Contours")==0) return ContoursEnum;
-	      else if (strcmp(name,"Parameters")==0) return ParametersEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"Vertices")==0) return VerticesEnum;
+	      if (strcmp(name,"Parameters")==0) return ParametersEnum;
+	      else if (strcmp(name,"Vertices")==0) return VerticesEnum;
 	      else if (strcmp(name,"Results")==0) return ResultsEnum;
 	      else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
@@ -505,9 +506,9 @@
 	      else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
 	      else if (strcmp(name,"VzHO")==0) return VzHOEnum;
-	      else if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"VzFS")==0) return VzFSEnum;
+	      if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
+	      else if (strcmp(name,"VzFS")==0) return VzFSEnum;
 	      else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
 	      else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
@@ -628,9 +629,9 @@
 	      else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum;
 	      else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
-	      else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
+	      if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
+	      else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
 	      else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
 	      else if (strcmp(name,"Regular")==0) return RegularEnum;
Index: /issm/trunk-jpl/src/m/classes/hydrologydc.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/hydrologydc.m	(revision 17421)
+++ /issm/trunk-jpl/src/m/classes/hydrologydc.m	(revision 17422)
@@ -16,4 +16,5 @@
 		transfer_flag            = 0;
 		leakage_factor           = 0;
+		basal_moulin_input       = NaN;
 
 		spcsediment_head         = NaN;
@@ -42,4 +43,11 @@
 		end 
 		% }}}
+		function self = initialize(self,md) % {{{
+			if isnan(self.basal_moulin_input),
+				self.basal_moulin_input=zeros(md.mesh.numberofvertices,1);
+				disp('      no hydrology.basal_moulin_input specified: values set as zero');
+			end
+
+		end % }}}
 		% {{{ function obj = setdefaultparameters(obj) 
 		function obj = setdefaultparameters(obj) 
@@ -88,4 +96,5 @@
 				md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1);
 	    end
+			md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'forcing',1);
 
 			md = checkfield(md,'fieldname','hydrology.spcsediment_head','forcing',1);
@@ -129,4 +138,5 @@
 				fielddisplay(obj,'leakage_factor','user defined leakage factor [m]');
 	    end
+			fielddisplay(obj,'basal_moulin_input','Figure out what it is');
 			disp(sprintf('   - for the sediment layer'));
 			fielddisplay(obj,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]');
@@ -165,4 +175,5 @@
 				WriteData(fid,'object',obj,'fieldname','leakage_factor','format','Double');
 	    end
+			WriteData(fid,'object',obj,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
 
 			WriteData(fid,'object',obj,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 17421)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 17422)
@@ -727,4 +727,6 @@
 			md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node');
 			md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node');
+
+			% Hydrologydc variables
 			if isa(md.hydrology,'hydrologydc');
 				md.hydrology.spcsediment_head=project3d(md,'vector',md.hydrology.spcsediment_head,'type','node','layer',1);
@@ -732,4 +734,5 @@
 				md.hydrology.mask_eplactive_node=project3d(md,'vector',md.hydrology.mask_eplactive_node,'type','node','layer',1);
 				md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1);
+				md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1);
 	    end
 
Index: /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.m	(revision 17421)
+++ /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.m	(revision 17422)
@@ -60,11 +60,11 @@
 	timestep=step;
 
-	FID = fopen(strcat(path,filesep,name,filesep,name,'.vtk',int2str(timestep),'.vtk'),'w+');
-	fprintf(FID,'# vtk DataFile Version 2.0 \n');
-	fprintf(FID,'Data for run %s \n',model.miscellaneous.name);
-	fprintf(FID,'ASCII \n');
-	fprintf(FID,'DATASET UNSTRUCTURED_GRID \n');
+	fid = fopen(strcat(path,filesep,name,filesep,name,'.vtk',int2str(timestep),'.vtk'),'w+');
+	fprintf(fid,'# vtk DataFile Version 2.0 \n');
+	fprintf(fid,'Data for run %s \n',model.miscellaneous.name);
+	fprintf(fid,'ASCII \n');
+	fprintf(fid,'DATASET UNSTRUCTURED_GRID \n');
 	
-	fprintf(FID,'POINTS %d float\n',num_of_points);
+	fprintf(fid,'POINTS %d float\n',num_of_points);
 	if(dim==3);
 		s='%f %f %f \n';
@@ -73,7 +73,7 @@
   end
 	P=[points zeros(num_of_points,3-dim)];
-	fprintf(FID,s,P');
+	fprintf(fid,s,P');
 	
-	fprintf(FID,'CELLS %d %d\n',num_of_elt,num_of_elt*(point_per_elt+1));
+	fprintf(fid,'CELLS %d %d\n',num_of_elt,num_of_elt*(point_per_elt+1));
 	s='%d';
 	for j=1:point_per_elt
@@ -81,10 +81,10 @@
   end
 	s=cell2mat(horzcat(s,{'\n'}));
-	fprintf(FID,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]');
+	fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]');
 	
-	fprintf(FID,'CELL_TYPES %d\n',num_of_elt);
+	fprintf(fid,'CELL_TYPES %d\n',num_of_elt);
 	s='%d\n';
-	fprintf(FID,s,celltype*ones(num_of_elt,1));
-	fprintf(FID,'POINT_DATA %s \n',num2str(num_of_points));
+	fprintf(fid,s,celltype*ones(num_of_elt,1));
+	fprintf(fid,'POINT_DATA %s \n',num2str(num_of_points));
 
 	%loop over the different solution structures
@@ -112,8 +112,8 @@
 					smallval=(abs(sol_struct{j}(timestep).(fieldnames{k}))<1.0e-20);
 					sol_struct{j}(timestep).(fieldnames{k})(smallval)=0.0;
-					fprintf(FID,'SCALARS %s float 1 \n',fieldnames{k});
-					fprintf(FID,'LOOKUP_TABLE default\n');
+					fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k});
+					fprintf(fid,'LOOKUP_TABLE default\n');
 					s='%e\n';
-					fprintf(FID,s,sol_struct{j}(timestep).(fieldnames{k}));
+					fprintf(fid,s,sol_struct{j}(timestep).(fieldnames{k}));
 		    end		
 	    end 
@@ -134,11 +134,11 @@
 				smallval=(abs(res_struct.(fieldnames{k}))<1.0e-20);
 				res_struct.(fieldnames{k})(smallval)=0.0;
-				fprintf(FID,'SCALARS %s float 1 \n',fieldnames{k});
-				fprintf(FID,'LOOKUP_TABLE default\n');
+				fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k});
+				fprintf(fid,'LOOKUP_TABLE default\n');
 				s='%e\n';
-				fprintf(FID,s,res_struct.(fieldnames{k}));
+				fprintf(fid,s,res_struct.(fieldnames{k}));
 	    end		
 		end 
 	end
-	fclose(FID);
+	fclose(fid);
 end
Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 17421)
+++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 17422)
@@ -126,4 +126,5 @@
 def HydrologydcPenaltyFactorEnum(): return StringToEnum("HydrologydcPenaltyFactor")[0]
 def HydrologydcPenaltyLockEnum(): return StringToEnum("HydrologydcPenaltyLock")[0]
+def HydrologydcBasalMoulinInputEnum(): return StringToEnum("HydrologydcBasalMoulinInput")[0]
 def HydrologyLayerEnum(): return StringToEnum("HydrologyLayer")[0]
 def HydrologySedimentEnum(): return StringToEnum("HydrologySediment")[0]
Index: /issm/trunk-jpl/test/NightlyRun/test3300.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test3300.m	(revision 17421)
+++ /issm/trunk-jpl/test/NightlyRun/test3300.m	(revision 17422)
@@ -5,4 +5,5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=1;
 md.hydrology.sedimentlimit_flag=1;
Index: /issm/trunk-jpl/test/NightlyRun/test332.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test332.m	(revision 17421)
+++ /issm/trunk-jpl/test/NightlyRun/test332.m	(revision 17422)
@@ -5,4 +5,5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=0;
 md.hydrology.sedimentlimit_flag=1;
Index: /issm/trunk-jpl/test/NightlyRun/test333.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test333.m	(revision 17421)
+++ /issm/trunk-jpl/test/NightlyRun/test333.m	(revision 17422)
@@ -5,4 +5,5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=1;
 md.hydrology.sedimentlimit_flag=1;
Index: /issm/trunk-jpl/test/NightlyRun/test334.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test334.m	(revision 17421)
+++ /issm/trunk-jpl/test/NightlyRun/test334.m	(revision 17422)
@@ -5,4 +5,5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=0;
 md.hydrology.sedimentlimit_flag=1;
Index: /issm/trunk-jpl/test/NightlyRun/test335.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test335.m	(revision 17421)
+++ /issm/trunk-jpl/test/NightlyRun/test335.m	(revision 17422)
@@ -5,4 +5,5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=1;
 md.hydrology.sedimentlimit_flag=1;