source: issm/trunk-jpl/src/m/materials/DepthAvgTempCond.m@ 17961

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;