Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17410)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17411)
@@ -7,5 +7,5 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 4
+//#define FSANALYTICAL 10
 
 /*Model processing*/
@@ -2810,9 +2810,11 @@
 	ElementMatrix* Ke1=CreateKMatrixFSViscous(element);
 	ElementMatrix* Ke2=CreateKMatrixFSFriction(element);
-	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+	ElementMatrix* Ke3=CreateKMatrixFSShelf(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
 
 	/*clean-up and return*/
 	delete Ke1;
 	delete Ke2;
+	delete Ke3;
 	return Ke;
 }/*}}}*/
@@ -2979,4 +2981,74 @@
 	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(D);
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/
+
+	if(!element->IsFloating() || !element->IsOnBed()) return NULL;
+
+	/*If on not water or not FS, skip stiffness: */
+	int approximation,shelf_dampening;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+	element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
+	if(shelf_dampening==0) return NULL;
+
+	/*Intermediaries*/
+	bool        mainlyfloating;
+	int         j,i,meshtype,dim;
+	IssmDouble  Jdet,slope2,scalar,dt;
+	IssmDouble  slope[2];
+	IssmDouble *xyz_list_base = NULL;
+	Gauss*      gauss         = 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 numdof    = vnumnodes*dim + pnumnodes;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	if(dt==0) dt=pow(10,5);
+	IssmDouble  rho_water =element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
+	Input*      surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=element->NewGaussBase(3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list_base,gauss);
+		element->NodalFunctionsVelocity(vbasis,gauss);
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		if(dim==2) slope2=slope[0]*slope[0];
+		else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1];
+		scalar  = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; 
+		for(i=0;i<vnumnodes;i++){
+			for(j=0;j<vnumnodes;j++){
+				Ke->values[numdof*(i*dim+1)+j*dim+1] += scalar*vbasis[i]*vbasis[j];
+			}
+		}
+	}
+
+	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list_base);
+	xDelete<IssmDouble>(vbasis);
 	return Ke;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 17410)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 17411)
@@ -69,4 +69,5 @@
 		ElementMatrix* CreateKMatrixFSViscous(Element* element);
 		ElementMatrix* CreateKMatrixFSFriction(Element* element);
+		ElementMatrix* CreateKMatrixFSShelf(Element* element);
 		ElementVector* CreatePVectorFS(Element* element);
 		ElementVector* CreatePVectorFSViscous(Element* element);