Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17948)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17949)
@@ -212,4 +212,5 @@
 	iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
 	iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
+	iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
 	iomodel->FetchDataToInput(elements,VxEnum,0.);
 	if(dakota_analysis)elements->InputDuplicate(VxEnum,QmuVxEnum);
@@ -3597,4 +3598,25 @@
 	element->TransformLoadVectorCoord(pe,cs_list);
 
+	/* shelf dampening*/
+	int shelf_dampening;
+	element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
+	if(shelf_dampening) {
+		Input*      mb_input=element->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input);
+		IssmDouble dt,mb,normal_b;
+		element->FindParam(&dt,TimesteppingTimeStepEnum);
+		for(int ig=gauss->begin();ig<gauss->end();ig++){
+			gauss->GaussPoint(ig);
+			element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+			element->NodalFunctionsVelocity(vbasis,gauss);
+			element->NormalBase(&normal[0],xyz_list_base);
+			if (dim==2) normal_b=normal[1];
+			else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
+			mb_input->GetInputValue(&mb, gauss);
+			for(i=0;i<vnumnodes;i++){
+				pe->values[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b;
+			}
+		}
+	}
+
 	/*Clean up and return*/
 	delete gauss;