Changeset 17963

Timestamp:

05/08/14 11:38:48 (11 years ago)

Author:

Eric.Larour

Message:

CHG: issues with shape of vectors.

File:

• issm/trunk-jpl/src/m/materials/DepthAvgTempCond.py (modified) (1 diff)

issm/trunk-jpl/src/m/materials/DepthAvgTempCond.py

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