Index: /issm/trunk-jpl/src/m/mech/robintemperature.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/robintemperature.py	(revision 17715)
+++ /issm/trunk-jpl/src/m/mech/robintemperature.py	(revision 17715)
@@ -0,0 +1,41 @@
+import numpy as npy
+from scipy.special import erf
+
+def robintemperature(heatflux,accumrate,thickness,surftemp,z):
+	'''
+	Compute vertical temperature profile of an ice sheet (Robin, 1955)
+
+	This routine computes the vertical temperature profile of an ice sheet
+	according to the solution of Robin (1955), neglecting friction and
+	horizontal advection.  The solution is thus most appropriate at an ice
+	divide.
+
+	The coordinate system for the solution runs from z=0 at the base 
+	to z=H at the surface of the ice.
+
+	Parameters (SI units):
+		-heatflux	Geothermal heat flux (W m^-2)
+		-accumrate	Surface accumulation rate (m s^-1 ice equivalent)
+		-thickness	Ice thickness (m)
+		-surftemp	Surface temperature (K)
+		-z				Vertical position at which to calculate temperature
+						(z can be a scalar or a vector)
+
+	Returns a vector the same length as z containing the temperature in K
+
+	Usage:
+		tprofile=robintemperature(heatflux,accumrate,thickness,surftemp,z)
+	'''
+
+	# some constants (from Holland and Jenkins, 1999)
+	alphaT=1.14e-6 # thermal diffusivity (m^2 s^-1)
+	c=2009. # specific heat capacity (J kg^-1 K^-1)
+	rho=917.  # ice density (kg m^-3)
+	
+	#create vertical coordinate variable
+	zstar=npy.sqrt(2.*alphaT*thickness/accumrate)
+	
+	tprofile=surftemp+npy.sqrt(2.*thickness*npy.pi/accumrate/alphaT)*(-heatflux)/2./rho/c*(erf(z/zstar)-erf(thickness/zstar))
+	
+	# difference between surface and base temperature for check (Cuffey2010 p412):
+	# print tprofile-surftemp