Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16363)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16364)
@@ -3456,5 +3456,66 @@
 	if(!isfront) return NULL;
 
-	_error_("STOP");
+	/*Intermediaries*/
+	IssmDouble  rho_ice,rho_water,gravity,y_g;
+	IssmDouble  Jdet,water_pressure,air_pressure,pressure;
+	IssmDouble  xyz_list_front[2][3];
+	IssmDouble  area_coordinates[2][3];
+	IssmDouble  normal[2];
+	GaussTria*  gauss = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize Element matrix and vectors*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	rho_water = matpar->GetRhoWater();
+	gravity   = matpar->GetG();
+	GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
+	GetSegmentNormal(&normal[0],xyz_list_front);
+
+	/*Start looping on Gaussian points*/
+	IssmDouble ymax=max(xyz_list_front[0][1],xyz_list_front[1][1]);
+	IssmDouble ymin=min(xyz_list_front[0][1],xyz_list_front[1][1]);
+	if(ymax>0. && ymin<0.) gauss=new GaussTria(area_coordinates,30); //refined in vertical because of the sea level discontinuity
+	else                   gauss=new GaussTria(area_coordinates,3);
+
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		y_g=GetYcoord(gauss);
+		GetNodalFunctionsVelocity(vbasis,gauss);
+		GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
+
+		water_pressure = rho_water*gravity*min(0.,y_g); //0 if the gaussian point is above water level
+		air_pressure   = 0.;
+		pressure       = water_pressure + air_pressure;
+
+		for(i=0;i<vnumnodes;i++){
+			pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
+			pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
+		}
+	}
+
+
+	/*Transform coordinate system*/
+	TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+
+	/*Clean up and return*/
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(vbasis);
+	delete gauss;
+	return pe;
+
 }
 /*}}}*/
@@ -3579,6 +3640,8 @@
 	}
 
+	pe->CheckConsistency();
 	/*Transform coordinate system*/
 	TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+	pe->CheckConsistency();
 
 	/*Clean up and return*/
@@ -3985,4 +4048,18 @@
 	delete gauss;
 	xDelete<int>(doflist);
+}
+/*}}}*/
+/*FUNCTION Tria::GetYcoord {{{*/
+IssmDouble Tria::GetYcoord(GaussTria* gauss){
+
+	IssmDouble y;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble y_list[NUMVERTICES];
+
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	for(int i=0;i<NUMVERTICES;i++) y_list[i]=xyz_list[i][1];
+	TriaRef::GetInputValue(&y,&y_list[0],gauss,P1Enum);
+
+	return y;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16363)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16364)
@@ -248,10 +248,11 @@
 		ElementVector* CreatePVectorStressbalanceFSShelf(void);
 		ElementMatrix* CreateJacobianStressbalanceSSA(void);
-		void	  GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution);
-		void	  GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solution);
-		void	  GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);
-		void	  InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
-		void	  InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
-		void	  InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
+		void	         GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution);
+		void	         GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solution);
+		void	         GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);
+		IssmDouble     GetYcoord(GaussTria* gauss);
+		void	         InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
+		void	         InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
+		void	         InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 16363)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 16364)
@@ -533,11 +533,20 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctions{{{*/
+/*FUNCTION TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){{{*/
 void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
 	_assert_(basis);
-
-	switch(this->element_type){
+	GetNodalFunctions(basis,gauss,this->element_type);
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation){{{*/
+void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	_assert_(basis);
+
+	switch(interpolation){
 		case P1Enum: case P1DGEnum:
 			basis[0]=gauss->coord1;
@@ -867,6 +876,12 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetInputValue{{{*/
+/*FUNCTION TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss){{{*/
 void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss){
+
+	GetInputValue(p,plist,gauss,this->element_type);
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss,int interpolation){{{*/
+void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss,int interpolation){
 
 	/*Output*/
@@ -874,9 +889,9 @@
 
 	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
+	int numnodes = this->NumberofNodes(interpolation);
 
 	/*Get nodal functions*/
 	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis, gauss);
+	GetNodalFunctions(basis, gauss,interpolation);
 
 	/*Calculate parameter for this Gauss point*/
@@ -888,8 +903,14 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::NumberofNodes{{{*/
+/*FUNCTION TriaRef::NumberofNodes(){{{*/
 int TriaRef::NumberofNodes(void){
 
-	switch(this->element_type){
+	return this->NumberofNodes(this->element_type);
+}
+/*}}}*/
+/*FUNCTION TriaRef::NumberofNodes(int interpolation){{{*/
+int TriaRef::NumberofNodes(int interpolation){
+
+	switch(interpolation){
 		case P1Enum:                return NUMNODESP1;
 		case P1DGEnum:              return NUMNODESP1;
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 16363)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 16364)
@@ -38,4 +38,5 @@
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
 		void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
+		void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int interpolation);
 		void GetNodalFunctionsVelocity(IssmDouble* basis, GaussTria* gauss);
 		void GetNodalFunctionsPressure(IssmDouble* basis, GaussTria* gauss);
@@ -48,7 +49,9 @@
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
+		void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss,int interpolation);
 		void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
 
 		int  NumberofNodes(void);
+		int  NumberofNodes(int interpolation);
 		int  NumberofNodesVelocity(void);
 		int  NumberofNodesPressure(void);