Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 19743)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 19744)
@@ -69,4 +69,5 @@
 					./classes/Loads/Penpair.cpp\
 					./classes/Loads/Pengrid.cpp\
+					./classes/Loads/Moulin.cpp\
 					./classes/Loads/Numericalflux.cpp\
 					./classes/matrix/ElementMatrix.cpp\
Index: /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp	(revision 19743)
+++ /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp	(revision 19744)
@@ -22,5 +22,29 @@
 }/*}}}*/
 void HydrologySommersAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-	/*No loads*/
+
+	/*Fetch parameters: */
+	int  hydrology_model;
+	iomodel->Constant(&hydrology_model,HydrologyModelEnum);
+
+	/*Now, do we really want Sommers?*/
+	if(hydrology_model!=HydrologysommersEnum) return;
+
+	//create penalties for nodes: no node can have water above the max
+	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(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum));
+			}
+		}
+		else if(reCast<int>(iomodel->Data(MeshVertexonbaseEnum)[i])){
+			if(iomodel->my_vertices[i]){
+				loads->AddObject(new Moulin(iomodel->loadcounter+i+1,i,iomodel,HydrologySommersAnalysisEnum));
+			}	
+		}
+	}
+	iomodel->DeleteData(1,MeshVertexonbaseEnum);
+
 }/*}}}*/
 void HydrologySommersAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
@@ -72,4 +96,5 @@
 	iomodel->FetchDataToInput(elements,HydrologyGapHeightEnum);
 	iomodel->FetchDataToInput(elements,HydrologyEnglacialInputEnum);
+	iomodel->FetchDataToInput(elements,HydrologyMoulinInputEnum);
 	iomodel->FetchDataToInput(elements,HydrologyBumpSpacingEnum);
 	iomodel->FetchDataToInput(elements,HydrologyBumpHeightEnum);
Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp	(revision 19744)
+++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp	(revision 19744)
@@ -0,0 +1,350 @@
+/*!\file Moulin.c
+ * \brief: implementation of the Moulin object
+ */
+
+/*Headers*/
+/*{{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../classes.h"
+#include "shared/shared.h"
+#include "../../analyses/analyses.h"
+/*}}}*/
+
+/*Element macros*/
+#define NUMVERTICES   1
+
+/*Moulin constructors and destructor*/
+Moulin::Moulin(){/*{{{*/
+	this->parameters=NULL;
+	this->hnode=NULL;
+	this->node=NULL;
+	this->helement=NULL;
+	this->element=NULL;
+	this->hmatpar=NULL;
+	this->matpar=NULL;
+}
+/*}}}*/
+Moulin::Moulin(int id, int index, IoModel* iomodel, int in_analysis_type){ //i is the element index/*{{{*/
+
+	int pengrid_node_id;
+	int pengrid_matpar_id;
+	int pengrid_element_id;
+
+	/*Some checks if debugging activated*/
+	_assert_(iomodel->singlenodetoelementconnectivity);
+	_assert_(index>=0 && index<iomodel->numberofvertices);
+	_assert_(id);
+
+	/*id: */
+	this->id=id;
+	this->analysis_type=in_analysis_type;
+
+	/*hooks: */
+	pengrid_node_id=iomodel->nodecounter+index+1;
+	pengrid_element_id=iomodel->singlenodetoelementconnectivity[index];
+	_assert_(pengrid_element_id);
+	pengrid_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
+
+	this->hnode=new Hook(&pengrid_node_id,1);
+	this->helement=new Hook(&pengrid_element_id,1);
+	this->hmatpar=new Hook(&pengrid_matpar_id,1);
+
+	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
+	this->parameters=NULL;
+	this->node=NULL;
+	this->element=NULL;
+	this->matpar=NULL;
+}
+/*}}}*/
+Moulin::~Moulin(){/*{{{*/
+	delete hnode;
+	delete helement;
+	delete hmatpar;
+	return;
+}
+/*}}}*/
+
+/*Object virtual functions definitions:*/
+Object* Moulin::copy() {/*{{{*/
+
+	Moulin* pengrid=NULL;
+
+	pengrid=new Moulin();
+
+	/*copy fields: */
+	pengrid->id=this->id;
+	pengrid->analysis_type=this->analysis_type;
+
+	/*point parameters: */
+	pengrid->parameters=this->parameters;
+
+	/*now deal with hooks and objects: */
+	pengrid->hnode=(Hook*)this->hnode->copy();
+	pengrid->hmatpar=(Hook*)this->hmatpar->copy();
+	pengrid->helement=(Hook*)this->helement->copy();
+
+	/*corresponding fields*/
+	pengrid->node  =(Node*)pengrid->hnode->delivers();
+	pengrid->matpar =(Matpar*)pengrid->hmatpar->delivers();
+	pengrid->element=(Element*)pengrid->helement->delivers();
+
+	return pengrid;
+}
+/*}}}*/
+void    Moulin::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
+
+	_assert_(this);
+
+	/*ok, marshall operations: */
+	MARSHALLING_ENUM(MoulinEnum);
+	MARSHALLING(id);
+	MARSHALLING(analysis_type);
+
+	if(marshall_direction==MARSHALLING_BACKWARD){
+		this->hnode      = new Hook();
+		this->helement   = new Hook();
+		this->hmatpar    = new Hook();
+	}
+
+	this->hnode->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
+	this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
+	this->hmatpar->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
+
+	/*corresponding fields*/
+	node   =(Node*)this->hnode->delivers();
+	matpar =(Matpar*)this->hmatpar->delivers();
+	element=(Element*)this->helement->delivers();
+}
+/*}}}*/
+void    Moulin::DeepEcho(void){/*{{{*/
+
+	_printf_("Moulin:\n");
+	_printf_("   id: " << id << "\n");
+	_printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
+	hnode->DeepEcho();
+	helement->DeepEcho();
+	hmatpar->DeepEcho();
+	_printf_("   parameters\n");
+	parameters->DeepEcho();
+}
+/*}}}*/
+void    Moulin::Echo(void){/*{{{*/
+	this->DeepEcho();
+}
+/*}}}*/
+int     Moulin::Id(void){ return id; }/*{{{*/
+/*}}}*/
+int     Moulin::ObjectEnum(void){/*{{{*/
+
+	return MoulinEnum;
+}
+/*}}}*/
+
+/*Load virtual functions definitions:*/
+void  Moulin::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
+	/*Take care of hooking up all objects for this load, ie links the objects in the hooks to their respective 
+	 * datasets, using internal ids and offsets hidden in hooks: */
+	hnode->configure(nodesin);
+	helement->configure(elementsin);
+	hmatpar->configure(materialsin);
+
+	/*Get corresponding fields*/
+	node=(Node*)hnode->delivers();
+	element=(Element*)helement->delivers();
+	matpar=(Matpar*)hmatpar->delivers();
+
+	/*point parameters to real dataset: */
+	this->parameters=parametersin;
+}
+/*}}}*/
+void  Moulin::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
+
+	/*No loads applied, do nothing: */
+	return;
+
+}
+/*}}}*/
+void  Moulin::CreatePVector(Vector<IssmDouble>* pf){/*{{{*/
+
+	ElementVector* pe=NULL;
+	int analysis_type;
+	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	switch(analysis_type){
+		
+	case HydrologySommersAnalysisEnum:
+		pe = this->CreatePVectorMoulin();
+		break;
+	default:
+		_error_("Don't know why we should be here");
+		/*No loads applied, do nothing: */
+		return;
+	}
+	if(pe){
+		pe->AddToGlobal(pf);
+		delete pe;
+	}
+
+}
+/*}}}*/
+void  Moulin::GetNodesLidList(int* lidlist){/*{{{*/
+
+	_assert_(lidlist);
+	_assert_(node);
+
+	lidlist[0]=node->Lid();
+}
+/*}}}*/
+void  Moulin::GetNodesSidList(int* sidlist){/*{{{*/
+
+	_assert_(sidlist);
+	_assert_(node);
+
+	sidlist[0]=node->Sid();
+}
+/*}}}*/
+int   Moulin::GetNumberOfNodes(void){/*{{{*/
+
+	return NUMVERTICES;
+}
+/*}}}*/
+bool  Moulin::InAnalysis(int in_analysis_type){/*{{{*/
+	if (in_analysis_type==this->analysis_type)return true;
+	else return false;
+}
+/*}}}*/
+bool  Moulin::IsPenalty(void){/*{{{*/
+	return true;
+}
+/*}}}*/
+void  Moulin::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
+
+	/*Don't do anything for now*/
+
+}
+/*}}}*/
+void  Moulin::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
+
+	/*Don't do anything for now*/
+}
+/*}}}*/
+void  Moulin::ResetHooks(){/*{{{*/
+
+	this->node=NULL;
+	this->element=NULL;
+	this->matpar=NULL;
+	this->parameters=NULL;
+
+	/*Get Element type*/
+	this->hnode->reset();
+	this->helement->reset();
+	this->hmatpar->reset();
+
+}
+/*}}}*/
+void  Moulin::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
+
+}
+/*}}}*/
+void  Moulin::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
+
+	/*Output */
+	int d_nz = 0;
+	int o_nz = 0;
+
+	if(!flags[this->node->Lid()]){
+
+		/*flag current node so that no other element processes it*/
+		flags[this->node->Lid()]=true;
+
+		int counter=0;
+		while(flagsindices[counter]>=0) counter++;
+		flagsindices[counter]=this->node->Lid();
+
+		/*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
+		switch(set2_enum){
+			case FsetEnum:
+				if(node->indexing.fsize){
+					if(this->node->IsClone())
+					 o_nz += 1;
+					else
+					 d_nz += 1;
+				}
+				break;
+			case GsetEnum:
+				if(node->indexing.gsize){
+					if(this->node->IsClone())
+					 o_nz += 1;
+					else
+					 d_nz += 1;
+				}
+				break;
+			case SsetEnum:
+				if(node->indexing.ssize){
+					if(this->node->IsClone())
+					 o_nz += 1;
+					else
+					 d_nz += 1;
+				}
+				break;
+			default: _error_("not supported");
+		}
+	}
+
+	/*Assign output pointers: */
+	*pd_nz=d_nz;
+	*po_nz=o_nz;
+}
+/*}}}*/
+
+/*Update virtual functions definitions:*/
+void  Moulin::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/
+	/*Nothing*/
+}
+/*}}}*/
+void  Moulin::InputUpdateFromConstant(int constant, int name){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Moulin::InputUpdateFromConstant(bool constant, int name){/*{{{*/
+
+	/*Don't do anything for now*/
+}
+/*}}}*/
+void  Moulin::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Moulin::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+void  Moulin::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
+	/*Nothing updated yet*/
+}
+/*}}}*/
+
+ElementVector* Moulin::CreatePVectorMoulin(void){/*{{{*/
+
+	/*If this node is not the master node (belongs to another partition of the
+	 * mesh), don't add the moulin input a second time*/
+	if(node->IsClone()) return NULL;
+
+	IssmDouble moulin_load;
+
+	/*Initialize Element matrix*/
+	ElementVector* pe=new ElementVector(&node,1,this->parameters);
+
+	this->element->GetInputValue(&moulin_load,node,HydrologyMoulinInputEnum);
+	pe->values[0]=moulin_load;
+
+	/*Clean up and return*/
+	return pe;
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Moulin.h	(revision 19744)
+++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.h	(revision 19744)
@@ -0,0 +1,85 @@
+/*!\file Moulin.h
+ * \brief: header file for pengrid object */
+
+#ifndef _MOULIN_H_
+#define _MOULIN_H_
+
+/*Headers:*/
+/*{{{*/
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include "./Load.h"
+class Hook;
+class Inputs;
+class Parameters;
+class IoModel;
+/*}}}*/
+
+class Moulin: public Load{
+
+	private: 
+
+		int id;
+		int analysis_type;
+
+		/*Hooks*/
+		Hook* hnode;  //hook to 1 node
+		Hook* helement;  //hook to 1 element
+		Hook* hmatpar; //hook to 1 matpar
+
+		/*Corresponding fields*/
+		Node    *node;
+		Element *element;
+		Matpar  *matpar;
+
+		Parameters* parameters; //pointer to solution parameters
+
+	public:
+
+		/*Moulin constructors, destructors {{{*/
+		Moulin();
+		Moulin(int index, int id, IoModel* iomodel,int analysis_type);
+		~Moulin();
+		/*}}}*/
+		/*Object virtual functions definitions:{{{ */
+		void  Echo();
+		void  DeepEcho();
+		int   Id(); 
+		int   ObjectEnum();
+		Object* copy();
+		void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
+		/*}}}*/
+		/*Update virtual functions resolution: {{{*/
+		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
+		void  InputUpdateFromMatrixDakota(IssmDouble* matrix ,int nrows, int ncols, int name, int type);
+		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
+		void  InputUpdateFromConstant(IssmDouble constant, int name);
+		void  InputUpdateFromConstant(int constant, int name);
+		void  InputUpdateFromConstant(bool constant, int name);
+		void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
+		/*}}}*/
+		/*Load virtual functions definitions: {{{*/
+		void  Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
+		void  SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
+		void  ResetHooks();
+		void  CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
+		void  CreatePVector(Vector<IssmDouble>* pf);
+		void  CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("Not implemented yet");};
+		void  GetNodesSidList(int* sidlist);
+		void  GetNodesLidList(int* lidlist);
+		int   GetNumberOfNodes(void);
+		bool  IsPenalty(void);
+		void  PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){_error_("Not implemented yet");};
+		void  PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* kfs, IssmDouble kmax);
+		void  PenaltyCreatePVector(Vector<IssmDouble>* pf, IssmDouble kmax);
+		void  SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
+		bool  InAnalysis(int analysis_type);
+		/*}}}*/
+
+		ElementVector* CreatePVectorMoulin(void);
+};
+
+#endif  /* _MOULIN_H_ */
Index: /issm/trunk-jpl/src/c/classes/classes.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/classes.h	(revision 19743)
+++ /issm/trunk-jpl/src/c/classes/classes.h	(revision 19744)
@@ -37,4 +37,5 @@
 #include "./Loads/Penpair.h"
 #include "./Loads/Pengrid.h"
+#include "./Loads/Moulin.h"
 
 /*Elements: */
Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 19743)
+++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 19744)
@@ -60,5 +60,4 @@
 		
 		/*Proceed now to heads computations*/
-		if(VerboseSolution()) _printf0_("   computing water head\n");
 		solutionsequence_hydro_nonlinear(femmodel);
 
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 19743)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 19744)
@@ -162,4 +162,5 @@
 	HydrologyBumpHeightEnum,
 	HydrologyEnglacialInputEnum,
+	HydrologyMoulinInputEnum,
 	HydrologyReynoldsEnum,
 	HydrologySpcheadEnum,
@@ -579,4 +580,5 @@
 	NumericalfluxTypeEnum,
 	ParamEnum,
+	MoulinEnum,
 	PengridEnum,
 	PenpairEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 19743)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 19744)
@@ -168,4 +168,5 @@
 		case HydrologyBumpHeightEnum : return "HydrologyBumpHeight";
 		case HydrologyEnglacialInputEnum : return "HydrologyEnglacialInput";
+		case HydrologyMoulinInputEnum : return "HydrologyMoulinInput";
 		case HydrologyReynoldsEnum : return "HydrologyReynolds";
 		case HydrologySpcheadEnum : return "HydrologySpchead";
@@ -571,4 +572,5 @@
 		case NumericalfluxTypeEnum : return "NumericalfluxType";
 		case ParamEnum : return "Param";
+		case MoulinEnum : return "Moulin";
 		case PengridEnum : return "Pengrid";
 		case PenpairEnum : return "Penpair";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 19743)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 19744)
@@ -171,4 +171,5 @@
 	      else if (strcmp(name,"HydrologyBumpHeight")==0) return HydrologyBumpHeightEnum;
 	      else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
+	      else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
 	      else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
 	      else if (strcmp(name,"HydrologySpchead")==0) return HydrologySpcheadEnum;
@@ -259,9 +260,9 @@
 	      else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
 	      else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum;
-	      else if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
+	      if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum;
+	      else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
 	      else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
 	      else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
@@ -382,9 +383,9 @@
 	      else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
 	      else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
-	      else if (strcmp(name,"SmbP")==0) return SmbPEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SmbSwf")==0) return SmbSwfEnum;
+	      if (strcmp(name,"SmbP")==0) return SmbPEnum;
+	      else if (strcmp(name,"SmbSwf")==0) return SmbSwfEnum;
 	      else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
 	      else if (strcmp(name,"SmbPAir")==0) return SmbPAirEnum;
@@ -505,9 +506,9 @@
 	      else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
 	      else if (strcmp(name,"SurfaceNormalVelocity")==0) return SurfaceNormalVelocityEnum;
-	      else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
+	      if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
+	      else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
 	      else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
 	      else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
@@ -583,4 +584,5 @@
 	      else if (strcmp(name,"NumericalfluxType")==0) return NumericalfluxTypeEnum;
 	      else if (strcmp(name,"Param")==0) return ParamEnum;
+	      else if (strcmp(name,"Moulin")==0) return MoulinEnum;
 	      else if (strcmp(name,"Pengrid")==0) return PengridEnum;
 	      else if (strcmp(name,"Penpair")==0) return PenpairEnum;
@@ -627,10 +629,10 @@
 	      else if (strcmp(name,"Free")==0) return FreeEnum;
 	      else if (strcmp(name,"Open")==0) return OpenEnum;
-	      else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
-	      else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
+	      if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
+	      else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
+	      else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
 	      else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
 	      else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
@@ -750,10 +752,10 @@
 	      else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
 	      else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
-	      else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
-	      else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
+	      if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
+	      else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
+	      else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
 	      else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
 	      else if (strcmp(name,"J")==0) return JEnum;
@@ -873,10 +875,10 @@
 	      else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
 	      else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
-	      else if (strcmp(name,"MisfitLocal")==0) return MisfitLocalEnum;
-	      else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
+	      if (strcmp(name,"MisfitLocal")==0) return MisfitLocalEnum;
+	      else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
+	      else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
 	      else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
 	      else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
Index: /issm/trunk-jpl/src/m/classes/hydrologysommers.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/hydrologysommers.m	(revision 19743)
+++ /issm/trunk-jpl/src/m/classes/hydrologysommers.m	(revision 19744)
@@ -11,4 +11,5 @@
 		bump_height     = NaN;
 		englacial_input = NaN;
+		moulin_input    = NaN;
 		reynolds        = NaN;
 		spchead         = NaN;
@@ -40,5 +41,6 @@
 			md = checkfield(md,'fieldname','hydrology.bump_spacing','>',0,'size',[md.mesh.numberofelements 1],'NaN',1);
 			md = checkfield(md,'fieldname','hydrology.bump_height','>=',0,'size',[md.mesh.numberofelements 1],'NaN',1);
-			md = checkfield(md,'fieldname','hydrology.englacial_input','>=',0,'size',[md.mesh.numberofvertices 1],'NaN',1);
+			md = checkfield(md,'fieldname','hydrology.englacial_input','>=',0,'NaN',1,'timeseries',1);
+			md = checkfield(md,'fieldname','hydrology.moulin_input','>=',0,'NaN',1,'timeseries',1);
 			md = checkfield(md,'fieldname','hydrology.reynolds','>',0,'size',[md.mesh.numberofelements 1],'NaN',1);
 			md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices 1]);
@@ -51,4 +53,5 @@
 			fielddisplay(self,'bump_height','characteristic bedrock bump height (m)');
 			fielddisplay(self,'englacial_input','liquid water input from englacial to subglacial system (m/yr)');
+			fielddisplay(self,'moulin_input','liquid water input from moulins (at the vertices) to subglacial system (m^3/s)');
 			fielddisplay(self,'reynolds','Reynolds'' number');
 			fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)');
@@ -63,5 +66,6 @@
 			WriteData(fid,'object',self,'class','hydrology','fieldname','bump_spacing','format','DoubleMat','mattype',2);
 			WriteData(fid,'object',self,'class','hydrology','fieldname','bump_height','format','DoubleMat','mattype',2);
-			WriteData(fid,'object',self,'class','hydrology','fieldname','englacial_input','format','DoubleMat','mattype',1,'scale',1./yts);
+			WriteData(fid,'object',self,'class','hydrology','fieldname','englacial_input','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1);
+			WriteData(fid,'object',self,'class','hydrology','fieldname','moulin_input','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
 			WriteData(fid,'object',self,'class','hydrology','fieldname','reynolds','format','DoubleMat','mattype',2);
 			WriteData(fid,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1);
Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m
===================================================================
--- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m	(revision 19743)
+++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m	(revision 19744)
@@ -80,5 +80,5 @@
 
 	case TransientSolutionEnum(),
-		analyses=[StressbalanceAnalysisEnum();StressbalanceVerticalAnalysisEnum();StressbalanceSIAAnalysisEnum();L2ProjectionBaseAnalysisEnum();ThermalAnalysisEnum();MeltingAnalysisEnum();EnthalpyAnalysisEnum();MasstransportAnalysisEnum()];
+		analyses=[StressbalanceAnalysisEnum();StressbalanceVerticalAnalysisEnum();StressbalanceSIAAnalysisEnum();L2ProjectionBaseAnalysisEnum();ThermalAnalysisEnum();MeltingAnalysisEnum();EnthalpyAnalysisEnum();MasstransportAnalysisEnum();HydrologySommersAnalysisEnum];
 
 	case FlaimSolutionEnum(),
Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 19743)
+++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 19744)
@@ -160,4 +160,5 @@
 def HydrologyBumpHeightEnum(): return StringToEnum("HydrologyBumpHeight")[0]
 def HydrologyEnglacialInputEnum(): return StringToEnum("HydrologyEnglacialInput")[0]
+def HydrologyMoulinInputEnum(): return StringToEnum("HydrologyMoulinInput")[0]
 def HydrologyReynoldsEnum(): return StringToEnum("HydrologyReynolds")[0]
 def HydrologySpcheadEnum(): return StringToEnum("HydrologySpchead")[0]
@@ -563,4 +564,5 @@
 def NumericalfluxTypeEnum(): return StringToEnum("NumericalfluxType")[0]
 def ParamEnum(): return StringToEnum("Param")[0]
+def MoulinEnum(): return StringToEnum("Moulin")[0]
 def PengridEnum(): return StringToEnum("Pengrid")[0]
 def PenpairEnum(): return StringToEnum("Penpair")[0]
Index: /issm/trunk-jpl/src/m/enum/HydrologyMoulinInputEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/HydrologyMoulinInputEnum.m	(revision 19744)
+++ /issm/trunk-jpl/src/m/enum/HydrologyMoulinInputEnum.m	(revision 19744)
@@ -0,0 +1,11 @@
+function macro=HydrologyMoulinInputEnum()
+%HYDROLOGYMOULININPUTENUM - Enum of HydrologyMoulinInput
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=HydrologyMoulinInputEnum()
+
+macro=StringToEnum('HydrologyMoulinInput');
Index: /issm/trunk-jpl/src/m/enum/MoulinEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/MoulinEnum.m	(revision 19744)
+++ /issm/trunk-jpl/src/m/enum/MoulinEnum.m	(revision 19744)
@@ -0,0 +1,11 @@
+function macro=MoulinEnum()
+%MOULINENUM - Enum of Moulin
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MoulinEnum()
+
+macro=StringToEnum('Moulin');
Index: /issm/trunk-jpl/src/m/plot/plot_transient_movie.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_transient_movie.m	(revision 19743)
+++ /issm/trunk-jpl/src/m/plot/plot_transient_movie.m	(revision 19744)
@@ -62,5 +62,5 @@
 
 			clf;
-			titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
+			titlestring=[field ' at time ' num2str(results(i).time) ' year'];
 			plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
 			apply_options_movie(md,options,titlestring);