Index: /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.m
===================================================================
--- /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.m	(revision 17961)
+++ /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.m	(revision 17961)
@@ -0,0 +1,25 @@
+function Tbar=DepthAvgTempCond(md)
+%DEPTHAVGTEMPCOND- 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;
+
+Tbar=zeros(md.mesh.numberofvertices,1);
+
+%find temperature average when we are below melting point: 
+pos=find( (Ts+alpha) <Tpmp);
+Tbar(pos)=Ts(pos)+alpha(pos)/2;
+
+pos=find( (Ts+alpha) >=Tpmp);
+Tbar(pos)=Tpmp+(Tpmp^2-Ts(pos).^2)/2./alpha(pos)+ Tpmp*(Ts(pos)-Tpmp)./alpha(pos);
+
+%on ice shelf, easier: 
+pos=find(md.mask.groundedice_levelset<=0);
+Tbar(pos)=(Ts(pos)+Tpmp)/2;
Index: /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.py
===================================================================
--- /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.py	(revision 17961)
+++ /issm/trunk-jpl/src/m/materials/DepthAvgTempCond.py	(revision 17961)
@@ -0,0 +1,28 @@
+def Tbar=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
+
+	Tbar=npy.zeros(md.mesh.numberofvertices,)
+
+	#find temperature average when we are below melting point: 
+	pos=numpy.nonzero( (Ts+alpha) < Tpmp)
+	Tbar[pos]=Ts[pos]+alpha[pos]/2
+
+	pos=numpy.nonzero( (Ts+alpha) >=Tpmp)
+	Tbar[pos]=Tpmp+(Tpmp**2-Ts[pos]**2)/2./alpha[pos]+ Tpmp*(Ts[pos]-Tpmp)./alpha[pos]
+
+	#on ice shelf, easier: 
+	pos=numpy.nonzero(md.mask.groundedice_levelset<=0)
+	Tbar[pos]=(Ts[pos]+Tpmp)/2;
+
+	return Tbar