Index: /issm/trunk-jpl/src/m/mech/robintemperature.m
===================================================================
--- /issm/trunk-jpl/src/m/mech/robintemperature.m	(revision 16316)
+++ /issm/trunk-jpl/src/m/mech/robintemperature.m	(revision 16316)
@@ -0,0 +1,42 @@
+function tprofile=robintemperature(heatflux,accumrate,thickness,surftemp,z)
+%ROBINTEMPERATURE - 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)
+
+%checks
+if nargin~=5
+	help robintemperature
+	error('bad usage - wrong number of arguments.')
+end
+
+%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=sqrt(2*alphaT*thickness/accumrate);
+
+tprofile=surftemp+sqrt(2*thickness*pi/accumrate/alphaT)*heatflux/2/rho/c*(erf(z/zstar)-erf(thickness/zstar));
+
+% difference between surface and base temperautre for check (Cuffey2010 p412):
+% tprofile-surftemp