Changeset 21759 for issm/branches/trunk-larour-NatGeoScience2016/src/py3/mech/steadystateiceshelftemp.py

Timestamp:

06/07/17 10:50:54 (8 years ago)

Author:

Eric.Larour

Message:

CHG: merged branch back to trunk-jpl 21754.

File:

• issm/branches/trunk-larour-NatGeoScience2016/src/py3/mech/steadystateiceshelftemp.py (modified) (2 diffs)

issm/branches/trunk-larour-NatGeoScience2016/src/py3/mech/steadystateiceshelftemp.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 
 def steadystateiceshelftemp(md,surfacetemp,basaltemp):
@@ -42 +42 @@
         temperature=(Ts+Tb)/2  # where wi~=0
         
-        pos=npy.nonzero(abs(wi)>=1e-4) # to avoid division by zero
+        pos=np.nonzero(abs(wi)>=1e-4) # to avoid division by zero
 
-        npy.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation
+        np.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation
         #calculate depth-averaged temperature (in Celsius)
         try:
-                temperature[pos]=-( (Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos] - (Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])*npy.exp(Hi[pos]*wi[pos]/ki) )/( Hi[pos]*(npy.exp(Hi[pos]*wi[pos]/ki)-1))
+                temperature[pos]=-( (Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos] - (Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])*np.exp(Hi[pos]*wi[pos]/ki) )/( Hi[pos]*(np.exp(Hi[pos]*wi[pos]/ki)-1))
         except FloatingPointError:
                 print('steadystateiceshelf warning: overflow encountered in multipy/divide/exp, trying another formulation.') 
-                temperature[pos]=-( ((Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos])/npy.exp(Hi[pos]*wi[pos]/ki) - Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])/( Hi[pos]*(1-npy.exp(-Hi[pos]*wi[pos]/ki)))
+                temperature[pos]=-( ((Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos])/np.exp(Hi[pos]*wi[pos]/ki) - Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])/( Hi[pos]*(1-np.exp(-Hi[pos]*wi[pos]/ki)))
         
         #temperature should not be less than surface temp
-        pos=npy.nonzero(temperature<Ts)
+        pos=np.nonzero(temperature<Ts)
         temperature[pos]=Ts[pos]
         
         # NaN where melt rates are too high (infinity/infinity in exponential)
-        pos=npy.nonzero(npy.isnan(temperature))
+        pos=np.nonzero(np.isnan(temperature))
         temperature[pos]=Ts[pos]