Index: /issm/trunk-jpl/src/m/materials/TMeltingPoint.m
===================================================================
--- /issm/trunk-jpl/src/m/materials/TMeltingPoint.m	(revision 17196)
+++ /issm/trunk-jpl/src/m/materials/TMeltingPoint.m	(revision 17196)
@@ -0,0 +1,15 @@
+function Tm=TMeltingPoint(reftemp, pressure)
+%TMELTINGPOINT- calculate pressure melting point of ice
+%
+%   reftemp is the melting temperature at atmospheric pressure (initialized in md.materials.meltingpoint)   
+%
+%   pressure is in Pa   
+%
+%   Usage:
+%   Tm=TMeltingPoint(md.materials.meltingpoint,pressure)
+
+%variables
+beta=7.9e-8; % K Pa^-1
+
+Tm=reftemp-beta*pressure;
+
Index: /issm/trunk-jpl/src/m/materials/TMeltingPoint.py
===================================================================
--- /issm/trunk-jpl/src/m/materials/TMeltingPoint.py	(revision 17196)
+++ /issm/trunk-jpl/src/m/materials/TMeltingPoint.py	(revision 17196)
@@ -0,0 +1,21 @@
+import numpy as npy
+
+def TMeltingPoint(reftemp,pressure):
+	'''
+	Calculate the pressure melting point of ice at a given pressure
+
+	reftemp is the melting temperature in K at atmospheric pressure (initialized in md.materials.meltingpoint)
+
+	pressure is in Pa
+
+	Usage:
+		Tm=TMeltingPoint(md.materials.meltingpoint,pressure)
+	'''
+
+	#variables
+	beta=7.9e-8
+
+	#ensure ref is same dimension as pressure
+	ref=reftemp*npy.ones_like(pressure)
+
+	return reftemp-beta*pressure
Index: /issm/trunk-jpl/src/m/materials/arrhenius.m
===================================================================
--- /issm/trunk-jpl/src/m/materials/arrhenius.m	(revision 17195)
+++ /issm/trunk-jpl/src/m/materials/arrhenius.m	(revision 17196)
@@ -45,12 +45,4 @@
 end
 
-    function Tm=TMeltingPoint(pressure)
-        Tm=GetThom(T0,pressure);
-    end
-    
-    function Th=GetThom(T,pressure)
-        Th=T-beta*pressure;
-    end
-
 %   values for Activation energy Q and pre-exponential constants from
 %   Grewe/Blatter 2009, p54
@@ -82,5 +74,5 @@
     
     function B=GetRigidity(T,w,pressure)
-        Thom=GetThom(T, pressure);
+        Thom=TMeltingPoint(T, pressure);
         A=GetA(Thom, w);
         B=1./(A.^(1/n));