Index: /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp	(revision 23092)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp	(revision 23093)
@@ -32,6 +32,9 @@
 
 	/*Add input to elements*/
+	iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
 	iomodel->FetchDataToInput(elements,"md.hydrology.drainage_rate",HydrologyDrainageRateEnum);
-	iomodel->FetchDataToInput(elements,"md.hydrology.watercolumn",WatercolumnEnum,0.);
+	iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",WatercolumnEnum,0.);
 }/*}}}*/
 void HydrologyPismAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
@@ -46,4 +49,10 @@
 	if(hydrology_model!=HydrologypismEnum) return;
 	parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
+
+	/*Requested outputs*/
+	iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
+	parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
+	if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
+	iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
 
 	/*Nothing else to add for now*/
@@ -98,4 +107,6 @@
 	/*Retrieve all inputs and parameters*/
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	IssmDouble  rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
 
 	/*Get water column and drainage rate*/
@@ -109,5 +120,5 @@
 
 	/*Add water*/
-	for(int i=0;i<numvertices;i++) watercolumn[i] += (meltingrate[i]-drainagerate[i])*dt;
+	for(int i=0;i<numvertices;i++) watercolumn[i] += (meltingrate[i]/rho_ice*rho_water-drainagerate[i])*dt;
 
 	/* Divide by connectivity, add degree of channelization as an input */
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23092)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23093)
@@ -771,4 +771,5 @@
 				analyses_temp[numanalyses++]=HydrologyShreveAnalysisEnum;
 				analyses_temp[numanalyses++]=HydrologyShaktiAnalysisEnum;
+				analyses_temp[numanalyses++]=HydrologyPismAnalysisEnum;
 				analyses_temp[numanalyses++]=HydrologyDCInefficientAnalysisEnum;
 				analyses_temp[numanalyses++]=HydrologyDCEfficientAnalysisEnum;
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 23092)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 23093)
@@ -150,4 +150,7 @@
 			}
 			else if(hydrology_model==HydrologyshaktiEnum){
+				/*Nothing to add*/
+			}
+			else if(hydrology_model==HydrologypismEnum){
 				/*Nothing to add*/
 			}
Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 23092)
+++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 23093)
@@ -13,11 +13,11 @@
 
 	/*intermediary*/
-	int							hydrology_model;
-	int							solution_type;
-	int							numoutputs				=	0;
-	bool						save_results;
-	bool						modify_loads			=	true;
-	char					**requested_outputs = NULL;
-	IssmDouble			ThawedNodes;
+	int          hydrology_model;
+	int          solution_type;
+	int          numoutputs        = 0;
+	bool         save_results;
+	bool         modify_loads      = true;
+	char       **requested_outputs = NULL;
+	IssmDouble   ThawedNodes;
 
 	/*first recover parameters common to all solutions*/
@@ -46,9 +46,9 @@
 	else if (hydrology_model==HydrologydcEnum){
 		/*intermediary: */
-		bool				isefficientlayer;
-		bool				isthermal;
-		int					step,hydroslices;
-		IssmDouble	time,init_time,hydrotime,yts;
-		IssmDouble	dt,hydrodt;
+		bool       isefficientlayer;
+		bool       isthermal;
+		int        step,hydroslices;
+		IssmDouble time,init_time,hydrotime,yts;
+		IssmDouble dt,hydrodt;
 
 		femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp	(revision 23092)
+++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp	(revision 23093)
@@ -204,4 +204,5 @@
 		case 2: return HydrologyshreveEnum;
 		case 3: return HydrologyshaktiEnum;
+		case 4: return HydrologypismEnum;
 		default: _error_("Marshalled hydrology code \""<<enum_in<<"\" not supported yet"); 
 	}
Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp	(revision 23092)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscSolver.cpp	(revision 23093)
@@ -115,5 +115,5 @@
 	KSPGetPC(ksp,&pc);
 	if (solver_type==MUMPSPACKAGE_LU){
-		#if defined(_HAVE_PETSCDEV_)
+		#if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=9)
 		PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
 		#else
Index: /issm/trunk-jpl/src/m/classes/hydrologypism.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/hydrologypism.m	(revision 23092)
+++ /issm/trunk-jpl/src/m/classes/hydrologypism.m	(revision 23093)
@@ -47,4 +47,5 @@
 			WriteData(fid,prefix,'name','md.hydrology.model','data',4,'format','Integer');
 			WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','drainage_rate','format','DoubleMat','mattype',1,'scale',1./(1000.*yts)); %from mm/yr to m/s
+			WriteData(fid,prefix,'data',{'Watercolumn'},'name','md.hydrology.requested_outputs','format','StringArray');
 		end % }}}
 	end
Index: /issm/trunk-jpl/src/m/classes/hydrologyshreve.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/hydrologyshreve.m	(revision 23092)
+++ /issm/trunk-jpl/src/m/classes/hydrologyshreve.m	(revision 23093)
@@ -8,5 +8,5 @@
 		spcwatercolumn     = NaN;
 		stabilization      = 0;
-    requested_outputs  = {};
+		requested_outputs  = {};
 	end
 	methods
Index: /issm/trunk-jpl/test/NightlyRun/test610.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test610.m	(revision 23093)
+++ /issm/trunk-jpl/test/NightlyRun/test610.m	(revision 23093)
@@ -0,0 +1,33 @@
+%Test Name: 79NorthPISMhydro
+md=triangle(model(),'../Exp/79North.exp',10000.);
+md=setmask(md,'../Exp/79NorthShelf.exp','');
+md=parameterize(md,'../Par/79North.par');
+md=setflowequation(md,'SSA','all');
+
+%Hydrology
+md.hydrology = hydrologypism();
+md.hydrology.drainage_rate      = 1000*ones(md.mesh.numberofvertices,1);
+md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
+md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1);
+md.basalforcings.groundedice_melting_rate = [1:md.mesh.numberofvertices]';
+md.transient.ishydrology        = 1;
+md.transient.issmb              = 0;
+md.transient.ismasstransport    = 0;
+md.transient.isstressbalance    = 0;
+md.transient.isthermal          = 0;
+md.transient.isgroundingline    = 0;
+
+md.cluster=generic('name',oshostname(),'np',3);
+md=solve(md,'Transient');
+
+%Plot to check result
+%plotmodel(md,'data',md.results.TransientSolution(3).Watercolumn - 3*(md.materials.rho_freshwater/md.materials.rho_ice*[1:md.mesh.numberofvertices]' -1),'caxis',[-1 1])
+
+%Fields and tolerances to track changes
+field_names     ={'WaterColumn1','WaterColumn2','WaterColumn3'};
+field_tolerances={1e-12,1e-12,1e-12};
+field_values={...
+	(md.results.TransientSolution(1).Watercolumn),...
+	(md.results.TransientSolution(2).Watercolumn),...
+	(md.results.TransientSolution(3).Watercolumn),...
+	};