Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16843)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16844)
@@ -1114,5 +1114,4 @@
 
 	/*Intermediaries*/
-	IssmDouble  rho_ice,rho_water,gravity;
 	IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
 	IssmDouble  surface_under_water,base_under_water,pressure;
@@ -1129,9 +1128,9 @@
 
 	/*Retrieve all inputs and parameters*/
-	Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input* bed_input      =element->GetInput(BedEnum);       _assert_(bed_input);
-	rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
-	rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
-	gravity   = element->GetMaterialParameter(ConstantsGEnum);
+	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* bed_input       = element->GetInput(BedEnum);       _assert_(bed_input);
+	IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
 	element->GetVerticesCoordinates(&xyz_list);
 	element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
@@ -1574,6 +1573,5 @@
 
 	/*Intermediaries*/
-	IssmDouble  rho_ice,rho_water,gravity;
-	IssmDouble  Jdet,surface,z_g,water_pressure,ice_pressure;
+	IssmDouble  Jdet,surface,z,water_pressure,ice_pressure;
 	IssmDouble  surface_under_water,base_under_water,pressure;
 	IssmDouble* xyz_list       = NULL;
@@ -1583,5 +1581,6 @@
 
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
+	int numnodes    = element->GetNumberOfNodes();
+	int numvertices = element->GetNumberOfVertices();
 
 	/*Initialize Element vector and other vectors*/
@@ -1590,8 +1589,8 @@
 
 	/*Retrieve all inputs and parameters*/
-	Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
-	rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
-	rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
-	gravity   = element->GetMaterialParameter(ConstantsGEnum);
+	Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
+	IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
 	element->GetVerticesCoordinates(&xyz_list);
 	element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
@@ -1599,6 +1598,6 @@
 
 	/*Initialize gauss points*/
-	IssmDouble zmax=xyz_list[0*3+2]; for(int i=1;i<6;i++) if(xyz_list[i*3+2]>zmax) zmax=xyz_list[i*3+2];
-	IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i<6;i++) if(xyz_list[i*3+2]<zmin) zmin=xyz_list[i*3+2];
+	IssmDouble zmax=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]>zmax) zmax=xyz_list[i*3+2];
+	IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]<zmin) zmin=xyz_list[i*3+2];
 	if(zmax>0 && zmin<0) gauss=element->NewGauss(xyz_list,xyz_list_front,3,10);//refined in vertical because of the sea level discontinuity
 	else                 gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
@@ -1609,10 +1608,10 @@
 		gauss->GaussPoint(ig);
 		surface_input->GetInputValue(&surface,gauss);
-		z_g=element->GetZcoord(gauss);
+		z=element->GetZcoord(gauss);
 		element->NodalFunctions(basis,gauss);
 		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
 
-		water_pressure = rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
-		ice_pressure   = rho_ice*gravity*(surface-z_g);
+		water_pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
+		ice_pressure   = rho_ice*gravity*(surface-z);
 		pressure       = ice_pressure + water_pressure;
 
@@ -1998,5 +1997,5 @@
 		element->NodalFunctionsVelocity(vbasis,gauss);
 
-		element->NormalSection(&normal[0],xyz_list_base);
+		element->NormalBase(&normal[0],xyz_list_base);
 		_assert_(normal[dim-1]<0.);
 		bed_input->GetInputValue(&bed, gauss);
@@ -2024,5 +2023,82 @@
 ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
 
-	return NULL;
+	/*If no front, return NULL*/
+	if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
+
+	/*Intermediaries*/
+	int         i,meshtype,dim;
+	IssmDouble  Jdet,pressure,surface,z;
+	IssmDouble	normal[3];
+	IssmDouble *xyz_list       = NULL;
+	IssmDouble *xyz_list_front = NULL;
+	Gauss      *gauss          = NULL;
+
+	/*Get basal element*/
+	if(!element->IsOnBed() || !element->IsFloating()) return NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; break;
+		case Mesh3DEnum:         dim = 3; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes   = element->GetNumberOfNodesVelocity();
+	int pnumnodes   = element->GetNumberOfNodesPressure();
+	int numvertices = element->GetNumberOfVertices();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+	else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize vectors*/
+	ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
+	element->NormalSection(&normal[0],xyz_list_front);
+	Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
+	IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble  gravity   = element->GetMaterialParameter(ConstantsGEnum);
+
+	/*Initialize gauss points*/
+	IssmDouble zmax=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]>zmax) zmax=xyz_list[i*3+2];
+	IssmDouble zmin=xyz_list[0*3+2]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+2]<zmin) zmin=xyz_list[i*3+2];
+	if(zmax>0 && zmin<0) gauss=element->NewGauss(xyz_list,xyz_list_front,3,30);//refined in vertical because of the sea level discontinuity
+	else                 gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
+
+	/* Start  looping on the number of gaussian points: */
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
+		element->NodalFunctionsVelocity(vbasis,gauss);
+		surface_input->GetInputValue(&surface,gauss);
+		z=element->GetZcoord(gauss);
+		pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
+
+		for (int i=0;i<vnumnodes;i++){
+			pe->values[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
+			pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
+			if(dim==3) pe->values[dim*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(vbasis);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_front);
+	return pe;
 }/*}}}*/
 void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16843)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16844)
@@ -65,4 +65,5 @@
 		virtual void   NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0;
 		virtual void   NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
+		virtual void   NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0;
 		virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0;
 		virtual void   TransformLoadVectorCoord(ElementVector* pe,int transformenum)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16843)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16844)
@@ -127,4 +127,5 @@
 		void        NormalSection(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
 		void        NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
+		void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
 		int         NumberofNodesVelocity(void){_error_("not implemented yet");};
 		int         NumberofNodesPressure(void){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16843)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16844)
@@ -261,4 +261,5 @@
 		void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
 		void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
+		void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
 		IssmDouble     GetMaterialParameter(int enum_in);
 		Input*         GetInput(int inputenum);