Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/mechanicalproperties.py	(revision 19465)
+++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.py	(revision 19466)
@@ -4,5 +4,5 @@
 from averaging import averaging
 
-def mechanicalproperties(md,vx,vy):
+def mechanicalproperties(md,vx,vy,**kwargs):
 	"""
 	MECHANICALPROPERTIES - compute stress and strain rate for a goven velocity
@@ -29,4 +29,13 @@
 	if npy.any(md.flowequation.element_equation!=2):
 		print 'Warning: the model has some non SSA elements. These will be treated like SSA elements'
+
+        #unpack kwargs
+	if 'damage' in kwargs: 
+	    damage=kwargs.pop('damage')
+            if len(damage)!=md.mesh.numberofvertices:
+		raise ValueError('if damage is supplied it should be of size ' + md.mesh.numberofvertices)
+            if npy.ndim(damage)==2:
+                damage=damage.reshape(-1,)
+        else: damage=None
 
 	if npy.ndim(vx)==2:
@@ -75,8 +84,7 @@
 		location=npy.nonzero(npy.logical_and(second_inv==0,power!=0))
 		nu[location]=10^18
-	elif 'matdamageice' in md.materials.__module__:
+	elif 'matdamageice' in md.materials.__module__ and damage is not None:
 		print 'computing damage-dependent properties!'
-                # FIXME might be using damage from a solution, not the initial md.damage.D
-		Zinv=npy.dot(1-md.damage.D[index-1],summation/3.).reshape(-1,)
+		Zinv=npy.dot(1-damage[index-1],summation/3.).reshape(-1,)
 		location=npy.nonzero(second_inv)
 		nu[location]=Zinv[location]*B_bar[location]/npy.power(second_inv[location],power[location])