source: issm/oecreview/Archive/16554-17801/ISSM-17714-17715.diff@ 17802

../trunk-jpl/src/m/mech/robintemperature.py

@@ +1 @@
+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