Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16803)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16804)
@@ -812,7 +812,67 @@
 		case SSAApproximationEnum: 
 			return CreatePVectorSSA(element);
+		case HOApproximationEnum: 
+			return CreatePVectorHO(element);
 		default:
 			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
 	}
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorHODrivingStress(element);
+	ElementVector* pe2=CreatePVectorHOFront(element);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  Jdet,slope[3];
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe=element->NewElementVector(HOApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input*     surface_input = element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+
+	/* 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);
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	delete gauss;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/
+
+	return NULL;
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
@@ -849,13 +909,9 @@
 
 	/*Intermediaries */
-	int         i;
-	IssmDouble  driving_stress_baseline,thickness;
-	IssmDouble  Jdet;
-	IssmDouble  slope[2];
-	IssmDouble* xyz_list     = NULL;
+	IssmDouble  thickness,Jdet,slope[2];
+	IssmDouble* xyz_list = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
-	IssmDouble* icefrontlevel= xNew<IssmDouble>(numnodes);
 
 	/*Initialize Element vector and vectors*/
@@ -865,8 +921,7 @@
 	/*Retrieve all inputs and parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
-	Input* thickness_input=element->GetInput(ThicknessEnum);          _assert_(thickness_input); 
-	Input* surface_input  =element->GetInput(SurfaceEnum);            _assert_(surface_input);
-	Input* drag_input     =element->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
-	element->GetInputListOnVertices(icefrontlevel,BedEnum);
+	Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 
+	Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
 	/* Start  looping on the number of gaussian points: */
@@ -881,10 +936,8 @@
 		thickness_input->GetInputValue(&thickness,gauss);
 		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
-		driving_stress_baseline=element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum)*thickness;
-
-		/*Build load vector: */
-		for(i=0;i<numnodes;i++){
-			pe->values[i*2+0]+=-driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1]+=-driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i];
+
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
 		}
 	}
@@ -894,4 +947,5 @@
 
 	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(basis);
 	delete gauss;
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16803)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16804)
@@ -23,4 +23,7 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		ElementVector* CreatePVectorHO(Element* element);
+		ElementVector* CreatePVectorHODrivingStress(Element* element);
+		ElementVector* CreatePVectorHOFront(Element* element);
 		ElementVector* CreatePVectorSSA(Element* element);
 		ElementVector* CreatePVectorSSADrivingStress(Element* element);
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 16803)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 16804)
@@ -115,6 +115,5 @@
 	/*Intermediaries*/
 	int         approximation;
-	IssmDouble  Jdet;
-	IssmDouble  dudx,dvdy,dwdz;
+	IssmDouble  Jdet,dudx,dvdy,dwdz;
 	IssmDouble  du[3],dv[3],dw[3];
 	IssmDouble* xyz_list = NULL;
@@ -167,10 +166,9 @@
 
 	/*Intermediaries */
-	int         i,j,approximation;
+	int         approximation;
 	IssmDouble *xyz_list      = NULL;
 	IssmDouble *xyz_list_base = NULL;
-	IssmDouble  Jdet;
-	IssmDouble  vx,vy,vz,dbdx,dbdy,basalmeltingvalue;
-	IssmDouble  slope[3];
+	IssmDouble  Jdet,slope[3];
+	IssmDouble  vx,vy,vz=0.,dbdx,dbdy,basalmeltingvalue;
 
 	if(!element->IsOnBed()) return NULL;
@@ -208,6 +206,4 @@
 			vzFS_input->GetInputValue(&vz,gauss);
 		}
-		else vz=0.;
-
 		dbdx=slope[0];
 		dbdy=slope[1];
@@ -216,5 +212,5 @@
 		element->NodalFunctions(basis,gauss);
 
-		for(i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
 	}
 
@@ -225,5 +221,4 @@
 	xDelete<IssmDouble>(xyz_list_base);
 	return pe;
-
 }/*}}}*/
 void StressbalanceVerticalAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/