Index: /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp	(revision 16859)
+++ /issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp	(revision 16860)
@@ -63,5 +63,95 @@
 }/*}}}*/
 ElementVector* BalancevelocityAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Intermediaries */
+	IssmDouble dhdt,mb,ms,Jdet;
+	IssmDouble gamma,thickness;
+	IssmDouble hnx,hny,dhnx[2],dhny[2];
+	IssmDouble *xyz_list  = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementVector* pe    = basalelement->NewElementVector();
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(numnodes*2);
+	IssmDouble*    H      = xNew<IssmDouble>(numnodes);
+	IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
+	IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	Input* ms_input   = basalelement->GetInput(SurfaceforcingsMassBalanceEnum);     _assert_(ms_input);
+	Input* mb_input   = basalelement->GetInput(BasalforcingsMeltingRateEnum);       _assert_(mb_input);
+	Input* dhdt_input = basalelement->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
+	Input* H_input    = basalelement->GetInput(ThicknessEnum);                      _assert_(H_input);
+	IssmDouble h = basalelement->CharacteristicLength();
+
+	/*Get vector N for all nodes*/
+	basalelement->GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
+	basalelement->GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
+	basalelement->GetInputListOnNodes(H,ThicknessEnum);
+	for(int i=0;i<numnodes;i++){
+		IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
+		Nx[i] = -H[i]*Nx[i]/norm;
+		Ny[i] = -H[i]*Ny[i]/norm;
+	}
+
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=basalelement->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
+		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		element->ValueP1DerivativesOnGauss(&dhnx[0],Nx,xyz_list,gauss);
+		element->ValueP1DerivativesOnGauss(&dhny[0],Ny,xyz_list,gauss);
+		element->ValueP1OnGauss(&hnx,Nx,gauss);
+		element->ValueP1OnGauss(&hny,Ny,gauss);
+
+		ms_input->GetInputValue(&ms,gauss);
+		mb_input->GetInputValue(&mb,gauss);
+		dhdt_input->GetInputValue(&dhdt,gauss);
+		H_input->GetInputValue(&thickness,gauss);
+		if(thickness<50.) thickness=50.;
+
+		gamma=h/(2.*thickness+1.e-10);
+
+		for(int i=0;i<numnodes;i++){
+			pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*( basis[i] + gamma*(basis[i]*(dhnx[0]+dhny[1])+hnx*dbasis[0*numnodes+i] + hny*dbasis[1*numnodes+i]));
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(H);
+	xDelete<IssmDouble>(Nx);
+	xDelete<IssmDouble>(Ny);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	return pe;
 }/*}}}*/
 void BalancevelocityAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp	(revision 16859)
+++ /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp	(revision 16860)
@@ -49,5 +49,72 @@
 }/*}}}*/
 ElementVector* SmoothedSurfaceSlopeXAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Intermediaries */
+	int         input_enum;
+	IssmDouble  Jdet,thickness,slope[2];
+	IssmDouble  taud_x,norms,normv,vx,vy;
+	IssmDouble *xyz_list  = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementVector* pe    = basalelement->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
+	Input* H_input       = basalelement->GetInput(ThicknessEnum); _assert_(H_input);
+	Input* surface_input = basalelement->GetInput(SurfaceEnum);   _assert_(surface_input);
+	Input* vx_input      = basalelement->GetInput(VxEnum);
+	Input* vy_input      = basalelement->GetInput(VyEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=basalelement->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
+
+		H_input->GetInputValue(&thickness,gauss);
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+		if(vx_input && vy_input){
+			vx_input->GetInputValue(&vx,gauss);
+			vy_input->GetInputValue(&vy,gauss);
+			norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10);
+			normv = sqrt(vx*vx + vy*vy);
+			if(normv>15./(365.*24.*3600.)) slope[0] = -vx/normv*norms;
+		}
+		taud_x = rho_ice*gravity*thickness*slope[0];
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*taud_x*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	return pe;
 }/*}}}*/
 void SmoothedSurfaceSlopeXAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp	(revision 16859)
+++ /issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp	(revision 16860)
@@ -49,5 +49,72 @@
 }/*}}}*/
 ElementVector* SmoothedSurfaceSlopeYAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Intermediaries */
+	int         input_enum;
+	IssmDouble  Jdet,thickness,slope[2];
+	IssmDouble  taud_y,norms,normv,vx,vy;
+	IssmDouble *xyz_list  = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementVector* pe    = basalelement->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
+	Input* H_input       = basalelement->GetInput(ThicknessEnum); _assert_(H_input);
+	Input* surface_input = basalelement->GetInput(SurfaceEnum);   _assert_(surface_input);
+	Input* vx_input      = basalelement->GetInput(VxEnum);
+	Input* vy_input      = basalelement->GetInput(VyEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=basalelement->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
+
+		H_input->GetInputValue(&thickness,gauss);
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+		if(vx_input && vy_input){
+			vx_input->GetInputValue(&vx,gauss);
+			vy_input->GetInputValue(&vy,gauss);
+			norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10);
+			normv = sqrt(vx*vx + vy*vy);
+			if(normv>15./(365.*24.*3600.)) slope[1] = -vy/normv*norms;
+		}
+		taud_y = rho_ice*gravity*thickness*slope[1];
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*taud_y*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	return pe;
 }/*}}}*/
 void SmoothedSurfaceSlopeYAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16859)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16860)
@@ -158,4 +158,6 @@
 		virtual void   ResetCoordinateSystem()=0;
 		virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0;
+		virtual void   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss)=0;
+		virtual void   ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss)=0;
 		virtual void   ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
 		virtual void   ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16859)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16860)
@@ -3473,4 +3473,14 @@
 			break;
 	}
+}
+/*}}}*/
+/*FUNCTION Penta::ValueP1OnGauss{{{*/
+void Penta::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){
+	PentaRef::GetInputValue(pvalue,values,gauss);
+}
+/*}}}*/
+/*FUNCTION Penta::ValueP1DerivativesOnGauss{{{*/
+void Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){
+	PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16859)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16860)
@@ -132,4 +132,6 @@
 		int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
 		IssmDouble TimeAdapt();
+		void   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
+		void   ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
 		int    VertexConnectivity(int vertexindex);
 		void   VerticalSegmentIndices(int** pindices,int* pnumseg);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16859)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16860)
@@ -138,4 +138,6 @@
 		IssmDouble  PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};
 		int         PressureInterpolation(void){_error_("not implemented yet");};
+		void        ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");};
+		void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		int         VelocityInterpolation(void){_error_("not implemented yet");};
 		Input*      GetInput(int inputenum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16859)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16860)
@@ -2946,4 +2946,14 @@
 	delete gauss;
 
+}
+/*}}}*/
+/*FUNCTION Tria::ValueP1OnGauss{{{*/
+void Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){
+	TriaRef::GetInputValue(pvalue,values,gauss);
+}
+/*}}}*/
+/*FUNCTION Tria::ValueP1DerivativesOnGauss{{{*/
+void Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){
+	TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16859)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16860)
@@ -139,4 +139,6 @@
 		void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
 		IssmDouble  TimeAdapt();
+		void   ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
+		void   ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
 		int         VertexConnectivity(int vertexindex);
 		void   VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};