Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 19266)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 19267)
@@ -112,4 +112,5 @@
 	iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
 	iomodel->FetchDataToInput(elements,BasalforcingsFloatingiceMeltingRateEnum);
+	iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
 	iomodel->FetchDataToInput(elements,VxEnum,0.);
 	iomodel->FetchDataToInput(elements,VyEnum,0.);
@@ -137,7 +138,11 @@
 ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
+	bool hack = false;
+
 	/*compute all stiffness matrices for this element*/
 	ElementMatrix* Ke1=CreateKMatrixVolume(element);
-	ElementMatrix* Ke2=CreateKMatrixSurface(element);
+	ElementMatrix* Ke2=NULL;
+	if(hack) Ke2=CreateKMatrixBase(element);
+	else Ke2=CreateKMatrixSurface(element);
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
 
@@ -147,4 +152,45 @@
 	return Ke;
 
+}/*}}}*/
+ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixBase(Element* element){/*{{{*/
+
+
+	if(!element->IsOnBase()) return NULL;
+
+	/*Intermediaries*/
+	IssmDouble  D,Jdet,normal[3];
+	IssmDouble *xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke    = element->NewElementMatrix(NoneApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss = element->NewGaussBase(2);
+	element->NormalBase(&normal[0],xyz_list);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantBase(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		D = -gauss->weight*Jdet*normal[2];
+
+		TripleMultiply( basis,1,numnodes,1,
+					&D,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	return Ke;
 }/*}}}*/
 ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
@@ -232,7 +278,11 @@
 ElementVector* StressbalanceVerticalAnalysis::CreatePVector(Element* element){/*{{{*/
 
+	bool hack = false;
+
 	/*compute all load vectors for this element*/
 	ElementVector* pe1=CreatePVectorVolume(element);
-	ElementVector* pe2=CreatePVectorBase(element);
+	ElementVector* pe2=NULL;
+	if(hack) pe2=CreatePVectorSurface(element);
+	else     pe2=CreatePVectorBase(element);
 	ElementVector* pe =new ElementVector(pe1,pe2);
 
@@ -308,4 +358,64 @@
 	return pe;
 }/*}}}*/
+ElementVector* StressbalanceVerticalAnalysis::CreatePVectorSurface(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	int         approximation;
+	IssmDouble *xyz_list      = NULL;
+	IssmDouble *xyz_list_surface= NULL;
+	IssmDouble  Jdet,slope[3];
+	IssmDouble  vx,vy,vz=0.,dsdx,dsdy;
+	IssmDouble  smb,smbvalue;
+
+	if(!element->IsOnSurface()) return NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetVerticesCoordinatesTop(&xyz_list_surface);
+	element->GetInputValue(&approximation,ApproximationEnum);
+	Input* surface_input    =element->GetInput(SurfaceEnum);               _assert_(surface_input);
+	Input* smb_input=element->GetInput(SurfaceforcingsMassBalanceEnum);    _assert_(smb_input);
+	Input* vx_input=element->GetInput(VxEnum);                             _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum);                             _assert_(vy_input);
+	Input* vzFS_input=NULL;
+	if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+		vzFS_input=element->GetInput(VzFSEnum);       _assert_(vzFS_input);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussTop(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		smb_input->GetInputValue(&smb,gauss);
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+			vzFS_input->GetInputValue(&vz,gauss);
+		}
+		dsdx=slope[0];
+		dsdy=slope[1];
+
+		element->JacobianDeterminantTop(&Jdet,xyz_list_surface,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dsdx+vy*dsdy-vz+smb)*basis[i];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_surface);
+	return pe;
+}/*}}}*/
 ElementVector* StressbalanceVerticalAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h	(revision 19266)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h	(revision 19267)
@@ -25,8 +25,10 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
+		ElementMatrix* CreateKMatrixBase(Element* element);
 		ElementMatrix* CreateKMatrixSurface(Element* element);
 		ElementMatrix* CreateKMatrixVolume(Element* element);
 		ElementVector* CreatePVector(Element* element);
 		ElementVector* CreatePVectorBase(Element* element);
+		ElementVector* CreatePVectorSurface(Element* element);
 		ElementVector* CreatePVectorVolume(Element* element);
 		void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);