Index: /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.py
===================================================================
--- /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.py	(revision 17962)
+++ /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.py	(revision 17963)
@@ -1,28 +1,34 @@
-def Tbar=DepthAvgTempCond(md)
+import numpy as npy
+from TMeltingPoint  import TMeltingPoint
+
+def DepthAvgTempCond(md):
    ''' compute conduction dependent temperature profile for an ice sheet. 
-
    Usage:
    Tbar=DepthAvgTempCond(md)
    '''
 
-	Tpmp=TMeltingPoint(md.materials.meltingpoint,0); #pressure melting point at 0 pressure.
-	k=md.materials.thermalconductivity
-	G=md.basalforcings.geothermalflux 
-	H=md.geometry.thickness
-	Ts=md.initialization.temperature
-	alpha=G*H/k
+   Tpmp=TMeltingPoint(md.materials.meltingpoint,0) #pressure melting point at 0 pressure.
 
-	Tbar=npy.zeros(md.mesh.numberofvertices,)
+   k=md.materials.thermalconductivity
+   G=md.basalforcings.geothermalflux
+   H=md.geometry.thickness[:,0]
+   Ts=md.initialization.temperature
+   alpha=G*H/k
 
-	#find temperature average when we are below melting point: 
-	pos=numpy.nonzero( (Ts+alpha) < Tpmp)
-	Tbar[pos]=Ts[pos]+alpha[pos]/2
+   Tbar=npy.zeros(md.mesh.numberofvertices,)
 
-	pos=numpy.nonzero( (Ts+alpha) >=Tpmp)
-	Tbar[pos]=Tpmp+(Tpmp**2-Ts[pos]**2)/2./alpha[pos]+ Tpmp*(Ts[pos]-Tpmp)./alpha[pos]
+   #find temperature average when we are below melting point: 
+   pos=npy.nonzero( Ts+alpha < Tpmp)
+   if pos:
+	   Tbar[pos]=Ts[pos]+alpha[pos]/2 
 
-	#on ice shelf, easier: 
-	pos=numpy.nonzero(md.mask.groundedice_levelset<=0)
-	Tbar[pos]=(Ts[pos]+Tpmp)/2;
+   pos=npy.nonzero( Ts+alpha>= Tpmp)
+   if pos:
+	   Tbar[pos]=Tpmp+(Tpmp**2-Ts[pos]**2)/2/alpha[pos]+ Tpmp*(Ts[pos]-Tpmp)/alpha[pos]
+   
+   #on ice shelf, easier: 
+   pos=npy.nonzero(md.mask.groundedice_levelset[0]<=0)
+   if pos:
+	   Tbar[pos]=(Ts[pos]+Tpmp)/2
 
-	return Tbar
+   return Tbar