source: issm/trunk/src/m/mech/robintemperature.m@ 17806

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,z)

%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 temperature for check (Cuffey2010 p412):
% tprofile-surftemp