Index: issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 16811)
+++ issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 16812)
@@ -188,5 +188,155 @@
 }/*}}}*/
 ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorVolume(element);
+	ElementVector* pe2=CreatePVectorSheet(element);
+	ElementVector* pe3=CreatePVectorShelf(element);
+	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	delete pe3;
+	return pe;
+}/*}}}*/
+ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int         stabilization;
+	IssmDouble  Jdet,phi,dt;
+	IssmDouble  enthalpy;
+	IssmDouble  kappa,tau_parameter,diameter;
+	IssmDouble  u,v,w;
+	IssmDouble  scalar_def,scalar_transient;
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes    = element->GetNumberOfNodes();
+	int numvertices = element->GetNumberOfVertices();
+
+	/*Initialize Element vector*/
+	ElementVector* pe             = element->NewElementVector();
+	IssmDouble*    basis          = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis         = xNew<IssmDouble>(3*numnodes);
+	IssmDouble*    pressure       = xNew<IssmDouble>(numvertices);
+	IssmDouble*    enthalpypicard = xNew<IssmDouble>(numvertices);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	element->FindParam(&stabilization,ThermalStabilizationEnum);
+	Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
+	Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
+	Input* enthalpy_input = NULL;
+	if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
+	if(stabilization==2){
+		element->GetInputListOnVertices(enthalpypicard,EnthalpyPicardEnum);
+		element->GetInputListOnVertices(pressure,PressureEnum);
+		diameter=element->MinEdgeLength(xyz_list);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
+
+		scalar_def=phi/rho_ice*Jdet*gauss->weight;
+		if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
+
+		/*TODO: add -beta*laplace T_m(p)*/
+		for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
+
+		/* Build transient now */
+		if(reCast<bool,IssmDouble>(dt)){
+			enthalpy_input->GetInputValue(&enthalpy, gauss);
+			scalar_transient=enthalpy*Jdet*gauss->weight;
+			for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
+		}
+
+		if(stabilization==2){
+			element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+			vx_input->GetInputValue(&u,gauss);
+			vy_input->GetInputValue(&v,gauss);
+			vz_input->GetInputValue(&w,gauss);
+			kappa          = element->EnthalpyDiffusionParameterVolume(numvertices,enthalpypicard,pressure) / rho_ice;
+			tau_parameter  = element->StabilizationParameter(u,v,w,diameter,kappa);
+
+			for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
+			if(reCast<bool,IssmDouble>(dt)){
+				for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]);
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(enthalpypicard);
+	xDelete<IssmDouble>(xyz_list);
+	delete gauss;
+	return pe;
+
+}/*}}}*/
+ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
+	return NULL;
+}/*}}}*/
+ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
+
+	IssmDouble  t_pmp,dt,Jdet,scalar_ocean,pressure;
+	IssmDouble *xyz_list_base = NULL;
+
+	/*Get basal element*/
+	if(!element->IsOnBed() || !element->IsFloating()) return NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input*      pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
+	IssmDouble  gravity             = element->GetMaterialParameter(ConstantsGEnum);
+	IssmDouble  rho_water           = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+	IssmDouble  mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
+	IssmDouble  thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussBase(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		pressure_input->GetInputValue(&pressure,gauss);
+		t_pmp=element->TMeltingPoint(pressure);
+
+		scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*(t_pmp)/(heatcapacity*rho_ice);
+		if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_ocean*basis[i];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list_base);
+	return pe;
 }/*}}}*/
 void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h
===================================================================
--- issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 16811)
+++ issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 16812)
@@ -23,4 +23,7 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		ElementVector* CreatePVectorVolume(Element* element);
+		ElementVector* CreatePVectorSheet(Element* element);
+		ElementVector* CreatePVectorShelf(Element* element);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
Index: issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16811)
+++ issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16812)
@@ -49,4 +49,6 @@
 		virtual void   DeleteMaterials(void)=0;
 		virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
+		virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0;
+		virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0;
 		virtual void   FindParam(int* pvalue,int paramenum)=0;
 		virtual void   FindParam(IssmDouble* pvalue,int paramenum)=0;
Index: issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16811)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16812)
@@ -873,4 +873,14 @@
 void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){
 	matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
+}
+/*}}}*/
+/*FUNCTION Penta::EnthalpyDiffusionParameter{{{*/
+IssmDouble Penta::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
+	return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
+}
+/*}}}*/
+/*FUNCTION Penta::EnthalpyDiffusionParameterVolume{{{*/
+IssmDouble Penta::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){
+	return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
 }
 /*}}}*/
@@ -3781,5 +3791,5 @@
 		GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum);
 		GetInputListOnVertices(&pressure[0],PressureEnum);
-		kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.);
+		kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpy,pressure); _assert_(kappa>0.);
 
 		D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
@@ -4250,5 +4260,5 @@
 			GetInputListOnVertices(&enthalpypicard[0],EnthalpyPicardEnum);
 			GetInputListOnVertices(&pressure_p[0],PressureEnum);
-			kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpypicard,pressure_p);
+			kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpypicard,pressure_p);
 			kappa/=rho_ice;
 			tau_parameter=StabilizationParameter(u,v,w,diameter,kappa); 
@@ -4808,5 +4818,5 @@
 			else{
 				enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
-				kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy, pressure); _assert_(kappa>0.);
+				kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,&enthalpy[0],&pressure[0]); _assert_(kappa>0.);
 				for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
 			}
Index: issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16811)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16812)
@@ -204,4 +204,6 @@
 		ElementVector* CreatePVectorFreeSurfaceBase(void);
 		ElementVector* CreatePVectorL2ProjectionBase(void);
+		IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
+		IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
 		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints);
 
Index: issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16811)
+++ issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16812)
@@ -84,4 +84,6 @@
 		void        Delta18oParameterization(void){_error_("not implemented yet");};
 		void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
+		IssmDouble  EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
+		IssmDouble  EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
 		void        FindParam(int* pvalue,int paramenum);
 		void        FindParam(IssmDouble* pvalue,int paramenum);
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16811)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16812)
@@ -239,4 +239,6 @@
 		ElementVector* CreatePVectorL2Projection(void);
 		ElementVector* CreatePVectorL2ProjectionBase(void);
+		IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
+		IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
 		IssmDouble     GetArea(void);
 		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints);
Index: issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16811)
+++ issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp	(revision 16812)
@@ -417,13 +417,14 @@
 /*}}}*/
 /*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/
-IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]){
-
-	IssmDouble lambda; // fraction of cold ice
-	IssmDouble kappa,kappa_c,kappa_t;  //enthalpy conductivities
-	IssmDouble Hc,Ht;
-	IssmDouble PIE[6],dHpmp[6];
-	int iv;
-
-	for (iv=0; iv<6; iv++){
+IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){
+
+	int         iv;
+	IssmDouble  lambda;                 // fraction of cold ice
+	IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
+	IssmDouble  Hc,Ht;
+	IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
+	IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);
+
+	for(iv=0; iv<numvertices; iv++){
 		PIE[iv]=PureIceEnthalpy(pressure[iv]);
 		dHpmp[iv]=enthalpy[iv]-PIE[iv];
@@ -432,9 +433,7 @@
 	bool allequalsign=true;
 	if(dHpmp[0]<0)
-		for(iv=1; iv<6;iv++)
-			allequalsign=(allequalsign && (dHpmp[iv]<0));
+		for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0));
 	else
-		for(iv=1; iv<6;iv++)
-			allequalsign=(allequalsign && (dHpmp[iv]>=0));
+		for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0));
 
 	if(allequalsign){
@@ -447,14 +446,18 @@
 		kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
 		Hc=0.; Ht=0.;
-		for(iv=0; iv<6;iv++){
+		for(iv=0; iv<numvertices;iv++){
 			if(enthalpy[iv]<PIE[iv])
-				Hc+=(PIE[iv]-enthalpy[iv]);
+			 Hc+=(PIE[iv]-enthalpy[iv]);
 			else
-				Ht+=(enthalpy[iv]-PIE[iv]);
+			 Ht+=(enthalpy[iv]-PIE[iv]);
 		}
 		_assert_((Hc+Ht)>0.);
-		lambda=Hc/(Hc+Ht);
-		kappa=1./(lambda/kappa_c + (1.-lambda)/kappa_t);
-	}
+		lambda = Hc/(Hc+Ht);
+		kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(PIE);
+	xDelete<IssmDouble>(dHpmp);
 	return kappa;
 }
Index: issm/trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 16811)
+++ issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 16812)
@@ -129,5 +129,5 @@
 		IssmDouble PureIceEnthalpy(IssmDouble pressure);
 		IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
-		IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]);
+		IssmDouble GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
 		IssmDouble GetLithosphereShearModulus();
 		IssmDouble GetLithosphereDensity();
