Changeset 21254

Timestamp:

10/11/16 01:24:07 (9 years ago)

Author:

bdef

Message:

CHG: replacing npy by np for numpy abreviation

Location:

issm/trunk-jpl/src/m

Files:

• boundaryconditions/PattynSMB.py (modified) (4 diffs)
• classes/outputdefinition.py (modified) (2 diffs)
• classes/qmu.py (modified) (6 diffs)
• classes/settings.py (modified) (1 diff)
• coordsystems/ll2xy.py (modified) (2 diffs)
• coordsystems/xy2ll.py (modified) (3 diffs)
• exp/expcoarsen.py (modified) (4 diffs)
• exp/expdisp.py (modified) (1 diff)
• extrusion/DepthAverage.py (modified) (4 diffs)
• extrusion/project2d.py (modified) (2 diffs)
• geometry/slope.py (modified) (2 diffs)
• interp/SectionValues.py (modified) (9 diffs)
• interp/averaging.py (modified) (4 diffs)
• interp/holefiller.py (modified) (3 diffs)
• interp/interp.py (modified) (7 diffs)
• inversions/parametercontroldrag.py (modified) (7 diffs)
• materials/TMeltingPoint.py (modified) (2 diffs)
• mech/analyticaldamage.py (modified) (4 diffs)
• mech/backstressfrominversion.py (modified) (2 diffs)
• mech/calcbackstress.py (modified) (2 diffs)
• mech/damagefrominversion.py (modified) (2 diffs)
• mech/mechanicalproperties.py (modified) (10 diffs)
• mech/robintemperature.py (modified) (2 diffs)
• mech/steadystateiceshelftemp.py (modified) (2 diffs)
• mech/thomasparams.py (modified) (6 diffs)
• plot/applyoptions.py (modified) (2 diffs)
• plot/colormaps/cmaptools.py (modified) (2 diffs)
• plot/plot_contour.py (modified) (1 diff)
• plot/plot_overlay.py (modified) (5 diffs)
• plot/plot_streamlines.py (modified) (3 diffs)
• plot/plotmodel.py (modified) (2 diffs)
• plot/processmesh.py (modified) (2 diffs)
• solve/parseresultsfromdisk.py (modified) (5 diffs)

issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py

@@ -1 +1 @@
 import os
-import numpy as npy
+import numpy as np
 def PattynSMB(md,Tf):
         """
@@ -25 +25 @@
 
         #Double check lat and long exist:
-        if npy.any(npy.isnan(md.mesh.lat)): 
+        if np.any(np.isnan(md.mesh.lat)): 
                 raise IOError('PattynSMB error message: md.mesh.lat field required')
 
@@ -36 +36 @@
         # Calculate summer temperature, Eqn (12)
         # No melting at PIG in mean conditions - need about 6 degress Tf to start having a negative yearly SMB
-        Tms = 16.81 - 0.00692*md.geometry.surface - 0.27937*npy.abs(md.mesh.lat) + Tf
+        Tms = 16.81 - 0.00692*md.geometry.surface - 0.27937*np.abs(md.mesh.lat) + Tf
         Tms= Tms[0]
 
@@ -43 +43 @@
 
         # Calculate Ablation, Eqn (10) (use for both Antarctica & Greenland), max melt is 10m/a
-        ABLdot=0.*npy.ones(md.mesh.numberofvertices)
-        pos=npy.nonzero(Tms>=0)
-        ABLdot[pos]=npy.minimum(1.4*Tms[pos],10)
+        ABLdot=0.*np.ones(md.mesh.numberofvertices)
+        pos=np.nonzero(Tms>=0)
+        ABLdot[pos]=np.minimum(1.4*Tms[pos],10)
 
         smb=ACCdot-ABLdot

issm/trunk-jpl/src/m/classes/outputdefinition.py

@@ -2 +2 @@
 from checkfield import checkfield
 from WriteData import WriteData
-import numpy as npy
+import numpy as np
 
 class outputdefinition(object):
@@ -40 +40 @@
                         data.append(classdefinition)
 
-                data=npy.unique(data);
+                data=np.unique(data);
                 WriteData(fid,prefix,'data',data,'name','md.outputdefinition.list','format','StringArray');
         # }}}

issm/trunk-jpl/src/m/classes/qmu.py

@@ -5 +5 @@
 from checkfield import checkfield
 from WriteData import WriteData
-import MatlabFuncs as m
 
 class qmu(object):
@@ -41 +40 @@
 
                 s+="%s\n" % fielddisplay(self,'isdakota','is qmu analysis activated?')
-                for i,variable in enumerate(self.variables.iteritems()):
-                        s+="         variables%s:  (arrays of each variable class)\n" % \
-                                        string_dim(self.variables,i)
+                for variable in self.variables:
+                        s+="         variables%s:  (arrays of each variable class)\n" % len(variable)
                         fnames=vars(variable)
                         maxlen=0
@@ -53 +51 @@
                                                 (maxlen+1,fname,size(getattr(variable,fname)),type(getattr(variable,fname)))
 
-                for i,response in enumerate(self.responses.iteritems()):
-                        s+="         responses%s:  (arrays of each response class)\n" % \
-                                        string_dim(self.responses,i)
+                for response in self.responses:
+                        s+="         responses%s:  (arrays of each response class)\n" % len(responses)
                         fnames=vars(response)
                         maxlen=0
@@ -67 +64 @@
                 s+="%s\n" % fielddisplay(self,'numberofresponses','number of responses') 
 
-                for i,method in enumerate(self.method.iteritems()):
+                for method in self.method:
                         if isinstance(method,'dakota_method'):
                                 s+="            method%s :    '%s'\n" % \
-                                                (string_dim(method,i),method.method)
+                                                (len(method),method.method)
 
-                for i,param in enumerate(self.params.iteritems()):
-                        s+="         params%s:  (array of method-independent parameters)\n" % \
-                                        string_dim(self.params,i)
+                for param in self.params:
+                        s+="         params%s:  (array of method-independent parameters)\n" % len(param)
                         fnames=vars(param)
                         maxlen=0
@@ -84 +80 @@
                                                 (maxlen+1,fname,any2str(getattr(param,fname)))
 
-                for i,result in enumerate(self.results.iteritems()):
-                        s+="         results%s:  (information from dakota files)\n" % \
-                                        string_dim(self.results,i)
+                for result in self.results:
+                        s+="         results%s:  (information from dakota files)\n" % len(self.result)
                         fnames=vars(result)
                         maxlen=0
@@ -132 +127 @@
                                 md.checkmessage("for qmu analysis, partitioning vector cannot go over npart, number of partition areas")
 
-                if not m.strcmpi(md.cluster.name,'none'):
+                if md.cluster.name!='none':
                         if not md.settings.waitonlock:
                                 md.checkmessage("waitonlock should be activated when running qmu in parallel mode!")

issm/trunk-jpl/src/m/classes/settings.py

@@ -16 +16 @@
                 self.lowmem              = 0
                 self.output_frequency    = 0
-                self.recording_frequency    = 0
+                self.recording_frequency = 0
                 self.waitonlock          = 0
 

issm/trunk-jpl/src/m/coordsystems/ll2xy.py

@@ -1 @@
-import numpy as npy 
+import numpy as np 
 
 def ll2xy(lat,lon,sgn=-1,central_meridian=0,standard_parallel=71):
@@ -36 +36 @@
         ex2 = .006693883
         # Eccentricity of the Hughes ellipsoid
-        ex = npy.sqrt(ex2)
+        ex = np.sqrt(ex2)
         
-        latitude = npy.abs(lat) * npy.pi/180.
-        longitude = (lon + delta) * npy.pi/180.
+        latitude = np.abs(lat) * np.pi/180.
+        longitude = (lon + delta) * np.pi/180.
         
         # compute X and Y in grid coordinates.
-        T = npy.tan(npy.pi/4-latitude/2) / ((1-ex*npy.sin(latitude))/(1+ex*npy.sin(latitude)))**(ex/2)
+        T = np.tan(np.pi/4-latitude/2) / ((1-ex*np.sin(latitude))/(1+ex*np.sin(latitude)))**(ex/2)
         
         if (90 - slat) <  1.e-5:
-                rho = 2.*re*T/npy.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex))
+                rho = 2.*re*T/np.sqrt((1.+ex)**(1.+ex)*(1.-ex)**(1.-ex))
         else:
-                sl  = slat*npy.pi/180.
-                tc  = npy.tan(npy.pi/4.-sl/2.)/((1.-ex*npy.sin(sl))/(1.+ex*npy.sin(sl)))**(ex/2.)
-                mc  = npy.cos(sl)/npy.sqrt(1.0-ex2*(npy.sin(sl)**2))
+                sl  = slat*np.pi/180.
+                tc  = np.tan(np.pi/4.-sl/2.)/((1.-ex*np.sin(sl))/(1.+ex*np.sin(sl)))**(ex/2.)
+                mc  = np.cos(sl)/np.sqrt(1.0-ex2*(np.sin(sl)**2))
                 rho = re*mc*T/tc
         
-        y = -rho * sgn * npy.cos(sgn*longitude)
-        x =  rho * sgn * npy.sin(sgn*longitude)
+        y = -rho * sgn * np.cos(sgn*longitude)
+        x =  rho * sgn * np.sin(sgn*longitude)
 
-        cnt1=npy.nonzero(latitude>= npy.pi/2.)[0]
+        cnt1=np.nonzero(latitude>= np.pi/2.)[0]
         
         if cnt1:

issm/trunk-jpl/src/m/coordsystems/xy2ll.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from math import pi
 
@@ -39 +39 @@
         # if x,y passed as lists, convert to numpy arrays
         if type(x) != "numpy.ndarray":
-                x=npy.array(x)
+                x=np.array(x)
         if type(y) != "numpy.ndarray":
-                y=npy.array(y)
+                y=np.array(y)
 
         ## Conversion constant from degrees to radians
@@ -50 +50 @@
         ex2 = .006693883
         ## Eccentricity of the Hughes ellipsoid
-        ex = npy.sqrt(ex2)
+        ex = np.sqrt(ex2)
         
         sl = slat*pi/180.
-        rho = npy.sqrt(x**2 + y**2)
-        cm = npy.cos(sl) / npy.sqrt(1.0 - ex2 * (npy.sin(sl)**2))
-        T = npy.tan((pi/4.0) - (sl/2.0)) / ((1.0 - ex*npy.sin(sl)) / (1.0 + ex*npy.sin(sl)))**(ex / 2.0)
+        rho = np.sqrt(x**2 + y**2)
+        cm = np.cos(sl) / np.sqrt(1.0 - ex2 * (np.sin(sl)**2))
+        T = np.tan((pi/4.0) - (sl/2.0)) / ((1.0 - ex*np.sin(sl)) / (1.0 + ex*np.sin(sl)))**(ex / 2.0)
         
         if abs(slat-90.) < 1.e-5:
-                T = rho*npy.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re
+                T = rho*np.sqrt((1. + ex)**(1. + ex) * (1. - ex)**(1. - ex)) / 2. / re
         else:
                 T = rho * T / (re * cm)
         
-        chi = (pi / 2.0) - 2.0 * npy.arctan(T)
+        chi = (pi / 2.0) - 2.0 * np.arctan(T)
         lat = chi + ((ex2 / 2.0) + (5.0 * ex2**2.0 / 24.0) + (ex2**3.0 / 12.0)) * \
-                npy.sin(2 * chi) + ((7.0 * ex2**2.0 / 48.0) + (29.0 * ex2**3 / 240.0)) * \
-                npy.sin(4.0 * chi) + (7.0 * ex2**3.0 / 120.0) * npy.sin(6.0 * chi) 
+                np.sin(2 * chi) + ((7.0 * ex2**2.0 / 48.0) + (29.0 * ex2**3 / 240.0)) * \
+                np.sin(4.0 * chi) + (7.0 * ex2**3.0 / 120.0) * np.sin(6.0 * chi) 
         
         lat = sgn * lat
-        lon = npy.arctan2(sgn * x,-sgn * y)
+        lon = np.arctan2(sgn * x,-sgn * y)
         lon = sgn * lon
         
-        res1 = npy.nonzero(rho <= 0.1)[0]
+        res1 = np.nonzero(rho <= 0.1)[0]
         if len(res1) > 0:
                 lat[res1] = 90. * sgn

issm/trunk-jpl/src/m/exp/expcoarsen.py

@@ -1 +1 @@
 import os.path
-import numpy as npy
+import numpy as np
 from collections import OrderedDict
 from expread import expread
@@ -34 +34 @@
         for contour in  contours:
                 
-                numpoints=npy.size(contour['x'])
+                numpoints=np.size(contour['x'])
 
                 j=0
@@ -44 +44 @@
 
                         #see whether we keep this point or not
-                        distance=npy.sqrt((x[j]-x[j+1])**2+(y[j]-y[j+1])**2)
+                        distance=np.sqrt((x[j]-x[j+1])**2+(y[j]-y[j+1])**2)
                         if distance<resolution and j<numpoints-2:   #do not remove last point
-                                x=npy.delete(x,j+1,0)
-                                y=npy.delete(y,j+1,0)
+                                x=np.delete(x,j+1,0)
+                                y=np.delete(y,j+1,0)
                                 numpoints=numpoints-1
                         else:
-                                division=int(npy.floor(distance/resolution)+1)
+                                division=int(np.floor(distance/resolution)+1)
                                 if division>=2:
-                                        xi=npy.linspace(x[j],x[j+1],division)
-                                        yi=npy.linspace(y[j],y[j+1],division)
+                                        xi=np.linspace(x[j],x[j+1],division)
+                                        yi=np.linspace(y[j],y[j+1],division)
                                         
-                                        x=npy.hstack((x[0:j+1],xi[1:-1],x[j+1:]))
-                                        y=npy.hstack((y[0:j+1],yi[1:-1],y[j+1:]))
+                                        x=np.hstack((x[0:j+1],xi[1:-1],x[j+1:]))
+                                        y=np.hstack((y[0:j+1],yi[1:-1],y[j+1:]))
 
                                         #update current point
@@ -65 +65 @@
                                         j=j+1
                 
-                if npy.size(x)>1:
+                if np.size(x)>1:
                         #keep the (x,y) contour arond
                         newcontour=OrderedDict()
-                        newcontour['nods']=npy.size(x)
+                        newcontour['nods']=np.size(x)
                         newcontour['density']=contour['density']
                         newcontour['x']=x

issm/trunk-jpl/src/m/exp/expdisp.py

@@ -1 +1 @@
 from expread import expread
-import numpy as npy
+import numpy as np
 
 def expdisp(domainoutline,ax,linestyle='--k',linewidth=1,unitmultiplier=1.):

issm/trunk-jpl/src/m/extrusion/DepthAverage.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from project2d import project2d
 
@@ -19 +19 @@
 
         # coerce to array in case float is passed
-        if type(vector)!=npy.ndarray:
+        if type(vector)!=np.ndarray:
                 print 'coercing array'
-                vector=npy.array(value)
+                vector=np.array(value)
 
         vec2d=False
@@ -30 +30 @@
         #nods data
         if vector.shape[0]==md.mesh.numberofvertices:
-                vector_average=npy.zeros(md.mesh.numberofvertices2d)
+                vector_average=np.zeros(md.mesh.numberofvertices2d)
                 for i in xrange(1,md.mesh.numberoflayers):
                         vector_average=vector_average+(project2d(md,vector,i)+project2d(md,vector,i+1))/2.*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i))
@@ -37 +37 @@
         #element data
         elif vector.shape[0]==md.mesh.numberofelements:
-                vector_average=npy.zeros(md.mesh.numberofelements2d)
+                vector_average=np.zeros(md.mesh.numberofelements2d)
                 for i in xrange(1,md.mesh.numberoflayers):
                         vector_average=vector_average+project2d(md,vector,i)*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i))

issm/trunk-jpl/src/m/extrusion/project2d.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 
 def project2d(md3d,value,layer):
@@ -25 +25 @@
         
         # coerce to array in case float is passed
-        if type(value)!=npy.ndarray:
+        if type(value)!=np.ndarray:
                 print 'coercing array'
-                value=npy.array(value)
+                value=np.array(value)
 
         vec2d=False

issm/trunk-jpl/src/m/geometry/slope.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from GetNodalFunctionsCoeff import  GetNodalFunctionsCoeff
 
@@ -33 +33 @@
         alpha,beta=GetNodalFunctionsCoeff(index,x,y)[0:2]
 
-        summation=npy.array([[1],[1],[1]])
-        sx=npy.dot(surf[index-1]*alpha,summation).reshape(-1,)
-        sy=npy.dot(surf[index-1]*beta,summation).reshape(-1,)
+        summation=np.array([[1],[1],[1]])
+        sx=np.dot(surf[index-1]*alpha,summation).reshape(-1,)
+        sy=np.dot(surf[index-1]*beta,summation).reshape(-1,)
 
-        s=npy.sqrt(sx**2+sy**2)
+        s=np.sqrt(sx**2+sy**2)
 
         if md.mesh.dimension()==3:
                 sx=project3d(md,'vector',sx,'type','element')
                 sy=project3d(md,'vector',sy,'type','element')
-                s=npy.sqrt(sx**2+sy**2)
+                s=np.sqrt(sx**2+sy**2)
 
         return (sx,sy,s)

issm/trunk-jpl/src/m/interp/SectionValues.py

@@ -1 +1 @@
 import os
 from expread import expread
-import numpy as npy
+import numpy as np
 from project2d import project2d
 #from InterpFromMesh2d import InterpFromMesh2d
@@ -41 +41 @@
 
         #initialization
-        X=npy.array([]) #X-coordinate
-        Y=npy.array([]) #Y-coordinate
-        S=npy.array([0.])  #curvilinear coordinate
+        X=np.array([]) #X-coordinate
+        Y=np.array([]) #Y-coordinate
+        S=np.array([0.])  #curvilinear coordinate
         
         for i in xrange(nods-1):
@@ -53 +53 @@
                 s_start=S[-1]
         
-                length_segment=npy.sqrt((x_end-x_start)**2+(y_end-y_start)**2)
-                portion=npy.ceil(length_segment/res_h)
+                length_segment=np.sqrt((x_end-x_start)**2+(y_end-y_start)**2)
+                portion=np.ceil(length_segment/res_h)
         
-                x_segment=npy.zeros(portion)
-                y_segment=npy.zeros(portion)
-                s_segment=npy.zeros(portion)
+                x_segment=np.zeros(portion)
+                y_segment=np.zeros(portion)
+                s_segment=np.zeros(portion)
 
                 for j in xrange(int(portion)):
@@ -66 +66 @@
         
                 #plug into X and Y
-                X=npy.append(X,x_segment)
-                Y=npy.append(Y,y_segment)
-                S=npy.append(S,s_segment)
+                X=np.append(X,x_segment)
+                Y=np.append(Y,y_segment)
+                S=np.append(S,s_segment)
 
-        X=npy.append(X,x[nods-1])
-        Y=npy.append(Y,y[nods-1])
+        X=np.append(X,x[nods-1])
+        Y=np.append(Y,y[nods-1])
         
         #Number of nodes:
@@ -77 +77 @@
         
         #Compute Z
-        Z=npy.zeros(numberofnodes)
+        Z=np.zeros(numberofnodes)
         
         #New mesh and Data interpolation
@@ -88 +88 @@
         
                 #Compute index
-                index=npy.array([range(1,numberofnodes),range(2,numberofnodes+1)]).T
+                index=np.array([range(1,numberofnodes),range(2,numberofnodes+1)]).T
         
         else:
@@ -102 +102 @@
         
                 #Some useful parameters
-                layers=int(npy.ceil(npy.mean(md.geometry.thickness)/res_v))
+                layers=int(np.ceil(np.mean(md.geometry.thickness)/res_v))
                 nodesperlayer=int(numberofnodes)
                 nodestot=int(nodesperlayer*layers)
@@ -109 +109 @@
         
                 #initialization
-                X3=npy.zeros(nodesperlayer*layers) 
-                Y3=npy.zeros(nodesperlayer*layers) 
-                Z3=npy.zeros(nodesperlayer*layers) 
-                S3=npy.zeros(nodesperlayer*layers) 
-                index3=npy.zeros((elementstot,4))
+                X3=np.zeros(nodesperlayer*layers) 
+                Y3=np.zeros(nodesperlayer*layers) 
+                Z3=np.zeros(nodesperlayer*layers) 
+                S3=np.zeros(nodesperlayer*layers) 
+                index3=np.zeros((elementstot,4))
         
                 #Get new coordinates in 3d
@@ -123 +123 @@
         
                         if i<layers-1:  #Build index3 with quads
-                                ids=npy.vstack((npy.arange(i,nodestot-layers,layers),npy.arange(i+1,nodestot-layers,layers),npy.arange(i+layers+1,nodestot,layers),npy.arange(i+layers,nodestot,layers))).T
+                                ids=np.vstack((np.arange(i,nodestot-layers,layers),np.arange(i+1,nodestot-layers,layers),np.arange(i+layers+1,nodestot,layers),np.arange(i+layers,nodestot,layers))).T
                                 index3[(i-1)*elementsperlayer:i*elementsperlayer,:]=ids
 
                 #Interpolation of data on specified points
-                data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,npy.nan)[0]
+                data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,np.nan)[0]
         
                 #build outputs

issm/trunk-jpl/src/m/interp/averaging.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from GetAreas import GetAreas
 import MatlabFuncs as m
@@ -39 +39 @@
         #initialization
         if layer==0:
-                weights=npy.zeros(md.mesh.numberofvertices,)
+                weights=np.zeros(md.mesh.numberofvertices,)
                 data=data.flatten(1)
         else:
-                weights=npy.zeros(md.mesh.numberofvertices2d,)
+                weights=np.zeros(md.mesh.numberofvertices2d,)
                 data=data[(layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:]
         
@@ -69 +69 @@
         index=index-1 # since python indexes starting from zero
         line=index.flatten(1)
-        areas=npy.vstack(areas).reshape(-1,)
-        summation=1./rep*npy.ones(rep,)
+        areas=np.vstack(areas).reshape(-1,)
+        summation=1./rep*np.ones(rep,)
         linesize=rep*numberofelements
         
         #update weights that holds the volume of all the element holding the node i
-        weights=csc_matrix( (npy.tile(areas,(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
+        weights=csc_matrix( (np.tile(areas,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
         
         #initialization
         if len(data)==numberofelements:
-                average_node=csc_matrix( (npy.tile(areas*data,(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
+                average_node=csc_matrix( (np.tile(areas*data,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
                 average_node=average_node/weights
                 average_node = csc_matrix(average_node)
@@ -85 +85 @@
 
         #loop over iteration
-        for i in npy.arange(1,iterations+1):
-                average_el=npy.asarray(npy.dot(average_node.todense()[index].reshape(numberofelements,rep),npy.vstack(summation))).reshape(-1,)
-                average_node=csc_matrix( (npy.tile(areas*average_el.reshape(-1),(rep,1)).reshape(-1,),(line,npy.zeros(linesize,))), shape=(numberofnodes,1))
+        for i in np.arange(1,iterations+1):
+                average_el=np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements,rep),np.vstack(summation))).reshape(-1,)
+                average_node=csc_matrix( (np.tile(areas*average_el.reshape(-1),(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1))
                 average_node=average_node/weights
                 average_node=csc_matrix(average_node)
         
         #return output as a full matrix (C code do not like sparse matrices)
-        average=npy.asarray(average_node.todense()).reshape(-1,)
+        average=np.asarray(average_node.todense()).reshape(-1,)
 
         return average

issm/trunk-jpl/src/m/interp/holefiller.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from scipy.spatial import cKDTree
 
@@ -30 +30 @@
         filled=data
         
-        XYGood=npy.dstack([x[goodids],y[goodids]])[0]
-        XYBad=npy.dstack([x[badids],y[badids]])[0]
+        XYGood=np.dstack([x[goodids],y[goodids]])[0]
+        XYBad=np.dstack([x[badids],y[badids]])[0]
         tree=cKDTree(XYGood)
         nearest=tree.query(XYBad,k=knn)[1]
@@ -42 +42 @@
                         for j in range(knn):
                                 neardat.append(filled[goodids][nearest[i][j]])
-                                filled[badids[i]]=npy.mean(neardat)
+                                filled[badids[i]]=np.mean(neardat)
                                 
         return filled

issm/trunk-jpl/src/m/interp/interp.py

@@ -1 +1 @@
 # module for inperpolating/smoothing data
-import numpy as npy
+import numpy as np
 from scipy.interpolate import CloughTocher2DInterpolator, Rbf
 from scipy.spatial import cKDTree
@@ -51 +51 @@
         xlim=[min(xi)-dx,max(xi)+dx]
         ylim=[min(yi)-dy,max(yi)+dy]
-        xflag=npy.logical_and(x>xlim[0],x<xlim[1])
-        yflag=npy.logical_and(y>ylim[0],y<ylim[1])
-        bothind=npy.nonzero(npy.logical_and(xflag,yflag))
+        xflag=np.logical_and(x>xlim[0],x<xlim[1])
+        yflag=np.logical_and(y>ylim[0],y<ylim[1])
+        bothind=np.nonzero(np.logical_and(xflag,yflag))
         subdata=data[bothind]
         subx=x[bothind]
         suby=y[bothind]
-        points=npy.array([subx,suby]).T
+        points=np.array([subx,suby]).T
 
         # mask out any nan's in the data and corresponding coordinate points
-        mask=npy.isnan(subdata)
-        ind=npy.nonzero(mask)[0]
+        mask=np.isnan(subdata)
+        ind=np.nonzero(mask)[0]
         if len(ind) and fill_nans:
                 print "         WARNING: filling nans using spline fit through good data points, which may or may not be appropriate. Check results carefully."
-        subdata=npy.delete(subdata,ind)
-        points=npy.delete(points,ind,axis=0)
+        subdata=np.delete(subdata,ind)
+        points=np.delete(points,ind,axis=0)
 
         if maxiter:
@@ -76 +76 @@
         if not fill_nans:
                 # identify nan's in xi,yi using nearest neighbors
-                xyinterp=npy.dstack([xi,yi])[0]
-                xg,yg=npy.meshgrid(subx,suby)
-                xydata=npy.dstack([subx,suby])[0]
+                xyinterp=np.dstack([xi,yi])[0]
+                xg,yg=np.meshgrid(subx,suby)
+                xydata=np.dstack([subx,suby])[0]
                 tree=cKDTree(xydata)
                 nearest=tree.query(xyinterp)[1]
-                pos=npy.nonzero(npy.isnan(subdata[nearest]))
+                pos=np.nonzero(np.isnan(subdata[nearest]))
                 interpdata[pos]=subdata[nearest][pos]
 
         return interpdata
 #}}}
-def GridSplineToMesh2d(x,y,data,xi,yi,default_value=npy.nan,plotonly=False,fill_nans=False):#{{{
+def GridSplineToMesh2d(x,y,data,xi,yi,default_value=np.nan,plotonly=False,fill_nans=False):#{{{
         '''
         python analog to InterpFromGridToMesh.  This routine uses
@@ -108 +108 @@
 
         Usage:
-                interpdata=GridToMesh(x,y,data,xi,yi,default_value=npy.nan,plotonly=False,fill_nans=False)
+                interpdata=GridToMesh(x,y,data,xi,yi,default_value=np.nan,plotonly=False,fill_nans=False)
 
         Examples:
@@ -114 +114 @@
         '''
 
-        if npy.ndim(x)==2:
+        if np.ndim(x)==2:
                 x=x.reshape(-1,)
-        if npy.ndim(y)==2:
+        if np.ndim(y)==2:
                 y=y.reshape(-1,)
         if len(x) != data.shape[1]+1 and len(x) != data.shape[1]:
@@ -134 +134 @@
         if len(x)==data.shape[1] and len(y)==data.shape[0]:
                 print '         x,y taken to define the center of data grid cells'
-                xind=npy.nonzero(npy.logical_and(x>xlim[0],x<xlim[1]))[0]
-                yind=npy.nonzero(npy.logical_and(y>ylim[0],y<ylim[1]))[0]
-                xg,yg=npy.meshgrid(x[xind],y[yind])
+                xind=np.nonzero(np.logical_and(x>xlim[0],x<xlim[1]))[0]
+                yind=np.nonzero(np.logical_and(y>ylim[0],y<ylim[1]))[0]
+                xg,yg=np.meshgrid(x[xind],y[yind])
                 subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
         elif len(x)==data.shape[1]+1 and len(y)==data.shape[0]+1:
                 print '         x,y taken to define the corners of data grid cells'
-                xcenter=npy.fromiter(((x[i]+x[i+1])/2 for i in range(len(x)-1)),npy.float)
-                ycenter=npy.fromiter(((y[i]+y[i+1])/2 for i in range(len(y)-1)),npy.float)
-                xind=npy.nonzero(npy.logical_and(xcenter>xlim[0],xcenter<xlim[1]))[0]
-                yind=npy.nonzero(npy.logical_and(ycenter>ylim[0],ycenter<ylim[1]))[0]
-                xg,yg=npy.meshgrid(xcenter[xind],ycenter[yind])
+                xcenter=np.fromiter(((x[i]+x[i+1])/2 for i in range(len(x)-1)),np.float)
+                ycenter=np.fromiter(((y[i]+y[i+1])/2 for i in range(len(y)-1)),np.float)
+                xind=np.nonzero(np.logical_and(xcenter>xlim[0],xcenter<xlim[1]))[0]
+                yind=np.nonzero(np.logical_and(ycenter>ylim[0],ycenter<ylim[1]))[0]
+                xg,yg=np.meshgrid(xcenter[xind],ycenter[yind])
                 subdata=data[yind[0]:yind[-1]+1,xind[0]:xind[-1]+1]
         else:
                 raise ValueError('x and y have inconsistent sizes: both should have length ncols(data)/nrows(data) or ncols(data)+1/nrows(data)+1')
 
-        points=npy.array([xg.ravel(),yg.ravel()]).T
+        points=np.array([xg.ravel(),yg.ravel()]).T
         flatsubdata=subdata.ravel()
 
         if plotonly:
-                plt.imshow(npy.flipud(subdata),origin='upper')
+                plt.imshow(np.flipud(subdata),origin='upper')
                 plt.show()
                 return
 
         # mask out any nan's in the data and corresponding coordinate points
-        mask=npy.isnan(flatsubdata)
-        ind=npy.nonzero(mask)[0]
+        mask=np.isnan(flatsubdata)
+        ind=np.nonzero(mask)[0]
         if len(ind) and fill_nans:
                 print "         WARNING: filling nans using spline fit through good data points, which may or may not be appropriate. Check results carefully."
-        goodsubdata=npy.delete(flatsubdata,ind)
-        goodpoints=npy.delete(points,ind,axis=0)
+        goodsubdata=np.delete(flatsubdata,ind)
+        goodpoints=np.delete(points,ind,axis=0)
 
         # create spline and index spline at mesh points
@@ -171 +171 @@
         if not fill_nans:
                 # identify nan's in xi,yi using nearest neighbors
-                xyinterp=npy.dstack([xi,yi])[0]
-                xydata=npy.dstack([xg.ravel(),yg.ravel()])[0]
+                xyinterp=np.dstack([xi,yi])[0]
+                xydata=np.dstack([xg.ravel(),yg.ravel()])[0]
                 tree=cKDTree(xydata)
                 nearest=tree.query(xyinterp)[1]
-                pos=npy.nonzero(npy.isnan(flatsubdata[nearest]))
+                pos=np.nonzero(np.isnan(flatsubdata[nearest]))
                 interpdata[pos]=flatsubdata[nearest][pos]
 

issm/trunk-jpl/src/m/inversions/parametercontroldrag.py

@@ -26 +26 @@
 
         #weights
-        weights=options.getfieldvalue('weights',npy.ones(md.mesh.numberofvertices))
-        if npy.size(weights)!=md.mesh.numberofvertices:
+        weights=options.getfieldvalue('weights',np.ones(md.mesh.numberofvertices))
+        if np.size(weights)!=md.mesh.numberofvertices:
                 md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices)
         else:
@@ -34 +34 @@
         #nsteps
         nsteps=options.getfieldvalue('nsteps',100);
-        if (npy.size(nsteps)!=1) | (nsteps<=0) | (floor(nsteps)!=nsteps):
+        if (np.size(nsteps)!=1) | (nsteps<=0) | (floor(nsteps)!=nsteps):
                 md.inversion.nsteps=100
         else:
@@ -41 +41 @@
         #cm_min
         cm_min=options.getfieldvalue('cm_min',ones(md.mesh.numberofvertices))
-        if (npy.size(cm_min)==1):
+        if (np.size(cm_min)==1):
                 md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices)
-        elif (npy.size(cm_min)==md.mesh.numberofvertices):
+        elif (np.size(cm_min)==md.mesh.numberofvertices):
                 md.inversion.min_parameters=cm_min
         else:
@@ -50 +50 @@
         #cm_max
         cm_max=options.getfieldvalue('cm_max',250*ones(md.mesh.numberofvertices))
-        if (npy.size(cm_max)==1):
+        if (np.size(cm_max)==1):
                 md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices)
-        elif (npy.size(cm_max)==md.mesh.numberofvertices):
+        elif (np.size(cm_max)==md.mesh.numberofvertices):
                 md.inversion.max_parameters=cm_max
         else:
@@ -59 +59 @@
         #eps_cm
         eps_cm=optoins.getfieldvalue('eps_cm',float('nan'))
-        if (npy.size(eps_cm)~=1 | eps_cm<0 ):
+        if (np.size(eps_cm)~=1 | eps_cm<0 ):
                 md.inversion.cost_function_threshold=float('nan')
         else:
@@ -66 +66 @@
         #maxiter
         maxiter=options.getfieldvalue('maxiter',10*ones(md.inversion.nsteps))
-        if (npy.any(maxiter<0) | npy.any(floor(maxiter)~=maxiter)):
+        if (np.any(maxiter<0) | np.any(floor(maxiter)~=maxiter)):
                 md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps)
         else:
@@ -75 +75 @@
         #cm_jump
         cm_jump=options.getfieldvalue('cm_jump',0.8*ones(md.inversion.nsteps))
-        if !npy.isreal(cm_jump):
+        if !np.isreal(cm_jump):
                 md.inversion.step_threshold=0.8*ones(md.inversion.nsteps)
         else:

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

@@ -1 @@
-import numpy as npy
+import numpy as np
 
 def TMeltingPoint(reftemp,pressure):
@@ -17 +17 @@
 
         #ensure ref is same dimension as pressure
-        ref=reftemp*npy.ones_like(pressure)
+        ref=reftemp*np.ones_like(pressure)
 
         return reftemp-beta*pressure

issm/trunk-jpl/src/m/mech/analyticaldamage.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from averaging import averaging
 #from plotmodel import plotmodel
@@ -54 +54 @@
 
         if isinstance(sigmab,(int,float)):
-                sigmab=sigmab*npy.ones((md.mesh.numberofvertices,))
+                sigmab=sigmab*np.ones((md.mesh.numberofvertices,))
 
         # check inputs
@@ -61 +61 @@
         if not '2d' in md.mesh.__doc__:
                 raise StandardError('only 2d (planview) model supported currently')
-        if npy.any(md.flowequation.element_equation!=2):
+        if np.any(md.flowequation.element_equation!=2):
                 print 'Warning: the model has some non SSA elements. These will be treated like SSA elements'
 
@@ -76 +76 @@
         n=averaging(md,md.materials.rheology_n,0)
         
-        D=1.-(1.+a+a**2+b**2)**((n-1.)/(2.*n))/npy.abs(ex)**(1./n)*(T-sigmab)/B/(2.+a)/npy.sign(ex)
+        D=1.-(1.+a+a**2+b**2)**((n-1.)/(2.*n))/np.abs(ex)**(1./n)*(T-sigmab)/B/(2.+a)/np.sign(ex)
         
         # D>1 where (2+a).*sign(ex)<0, compressive regions where high backstress needed
-        pos=npy.nonzero(D>1)
+        pos=np.nonzero(D>1)
         D[pos]=0
         
-        backstress=npy.zeros((md.mesh.numberofvertices,))
+        backstress=np.zeros((md.mesh.numberofvertices,))
 
         # backstress to bring D down to one 
-        backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*npy.sign(ex[pos])*(2.+a[pos])*npy.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
+        backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*np.sign(ex[pos])*(2.+a[pos])*np.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
         
-        pos=npy.nonzero(D<0)
+        pos=np.nonzero(D<0)
         #mask=ismember(1:md.mesh.numberofvertices,pos);
         D[pos]=0
         
         # backstress to bring negative damage to zero
-        backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*npy.sign(ex[pos])*(2.+a[pos])*npy.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
+        backstress[pos]=T[pos]-(1.-D[pos])*B[pos]*np.sign(ex[pos])*(2.+a[pos])*np.abs(ex[pos])**(1./n[pos])/(1.+a[pos]+a[pos]**2)**((n[pos]-1.)/2./n[pos])
         
-        pos=npy.nonzero(backstress<0)
+        pos=np.nonzero(backstress<0)
         backstress[pos]=0
         
         # rigidity from Thomas relation for D=0 and backstress=0
-        B=npy.sign(ex)/(2.+a)*(1.+a+a**2)**((n-1.)/2./n)*T/(npy.abs(ex)**(1./n))
-        pos=npy.nonzero(B<0)
+        B=np.sign(ex)/(2.+a)*(1.+a+a**2)**((n-1.)/2./n)*T/(np.abs(ex)**(1./n))
+        pos=np.nonzero(B<0)
         B[pos]=md.materials.rheology_B[pos]
         

issm/trunk-jpl/src/m/mech/backstressfrominversion.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from averaging import averaging
 from thomasparams import thomasparams
@@ -65 +65 @@
         if tempmask:
                 Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar
-                pos=npy.nonzero(Bi>md.materials.rheology_B)
+                pos=np.nonzero(Bi>md.materials.rheology_B)
                 Bi[pos]=md.materials.rheology_B[pos]
         
         # analytical backstress solution
-        backstress=T-Bi*npy.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
-        backstress[npy.nonzero(backstress<0)]=0
+        backstress=T-Bi*np.sign(ex0)*(2+a0)*np.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
+        backstress[np.nonzero(backstress<0)]=0
 
         return backstress

issm/trunk-jpl/src/m/mech/calcbackstress.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from averaging import averaging
 from thomasparams import thomasparams
@@ -61 +61 @@
         
         # analytical backstress solution
-        backstress=T-(1.-D)*B*npy.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
-        backstress[npy.nonzero(backstress<0)]=0
+        backstress=T-(1.-D)*B*np.sign(ex0)*(2+a0)*np.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
+        backstress[np.nonzero(backstress<0)]=0
 
         return backstress

issm/trunk-jpl/src/m/mech/damagefrominversion.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 
 def damagefrominversion(md):
@@ -24 +24 @@
         if any(md.flowequation.element_equation!=2):
                 raise StandardError('Warning: the model has some non-SSA elements.  These will be treated like SSA elements')
-        if npy.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2:
+        if np.ndim(md.results.StressbalanceSolution.MaterialsRheologyBbar)==2:
                 Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar.reshape(-1,)
         else:
                 Bi=md.results.StressbalanceSolution.MaterialsRheologyBbar
-        if npy.ndim(md.materials.rheology_B)==2:
+        if np.ndim(md.materials.rheology_B)==2:
                 BT=md.materials.rheology_B.reshape(-1,)
         else:
                 BT=md.materials.rheology_B
 
-        damage=npy.zeros_like(Bi)
+        damage=np.zeros_like(Bi)
 
         # Damage where Bi softer than B(T)
-        pos=npy.nonzero(Bi<BT)[0]
+        pos=np.nonzero(Bi<BT)[0]
         damage[pos]=1.-Bi[pos]/BT[pos]
         
-        pos=npy.nonzero(damage<0)
+        pos=np.nonzero(damage<0)
         damage[pos]=0
 

issm/trunk-jpl/src/m/mech/mechanicalproperties.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from GetNodalFunctionsCoeff import GetNodalFunctionsCoeff
 from results import results
@@ -27 +27 @@
         #       raise StandardError('only 2D model supported currently')
 
-        if npy.any(md.flowequation.element_equation!=2):
+        if np.any(md.flowequation.element_equation!=2):
                 print 'Warning: the model has some non SSA elements. These will be treated like SSA elements'
 
@@ -35 +35 @@
             if len(damage)!=md.mesh.numberofvertices:
                 raise ValueError('if damage is supplied it should be of size ' + md.mesh.numberofvertices)
-            if npy.ndim(damage)==2:
+            if np.ndim(damage)==2:
                 damage=damage.reshape(-1,)
         else: damage=None
 
-        if npy.ndim(vx)==2:
+        if np.ndim(vx)==2:
                 vx=vx.reshape(-1,)
-        if npy.ndim(vy)==2:
+        if np.ndim(vy)==2:
                 vy=vy.reshape(-1,)
         
@@ -48 +48 @@
         numberofvertices=md.mesh.numberofvertices
         index=md.mesh.elements
-        summation=npy.array([[1],[1],[1]])
-        directionsstress=npy.zeros((numberofelements,4))
-        directionsstrain=npy.zeros((numberofelements,4))
-        valuesstress=npy.zeros((numberofelements,2))
-        valuesstrain=npy.zeros((numberofelements,2))
+        summation=np.array([[1],[1],[1]])
+        directionsstress=np.zeros((numberofelements,4))
+        directionsstrain=np.zeros((numberofelements,4))
+        valuesstress=np.zeros((numberofelements,2))
+        valuesstrain=np.zeros((numberofelements,2))
         
         #compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
@@ -60 +60 @@
         vxlist=vx[index-1]/md.constants.yts
         vylist=vy[index-1]/md.constants.yts
-        ux=npy.dot((vxlist*alpha),summation).reshape(-1,)
-        uy=npy.dot((vxlist*beta),summation).reshape(-1,)
-        vx=npy.dot((vylist*alpha),summation).reshape(-1,)
-        vy=npy.dot((vylist*beta),summation).reshape(-1,)
+        ux=np.dot((vxlist*alpha),summation).reshape(-1,)
+        uy=np.dot((vxlist*beta),summation).reshape(-1,)
+        vx=np.dot((vylist*alpha),summation).reshape(-1,)
+        vy=np.dot((vylist*beta),summation).reshape(-1,)
         uyvx=(vx+uy)/2.
         #clear vxlist vylist
         
         #compute viscosity
-        nu=npy.zeros((numberofelements,))
-        B_bar=npy.dot(md.materials.rheology_B[index-1],summation/3.).reshape(-1,)
+        nu=np.zeros((numberofelements,))
+        B_bar=np.dot(md.materials.rheology_B[index-1],summation/3.).reshape(-1,)
         power=((md.materials.rheology_n-1.)/(2.*md.materials.rheology_n)).reshape(-1,)
         second_inv=(ux**2.+vy**2.+((uy+vx)**2.)/4.+ux*vy).reshape(-1,)
         
         #some corrections
-        location=npy.nonzero(npy.logical_and(second_inv==0,power!=0))
+        location=np.nonzero(np.logical_and(second_inv==0,power!=0))
         nu[location]=10^18      #arbitrary maximum viscosity to apply where there is no effective shear
         
         if 'matice' in md.materials.__module__:
-                location=npy.nonzero(second_inv)
+                location=np.nonzero(second_inv)
                 nu[location]=B_bar[location]/(second_inv[location]**power[location])
-                location=npy.nonzero(npy.logical_and(second_inv==0,power==0))
+                location=np.nonzero(np.logical_and(second_inv==0,power==0))
                 nu[location]=B_bar[location]
-                location=npy.nonzero(npy.logical_and(second_inv==0,power!=0))
+                location=np.nonzero(np.logical_and(second_inv==0,power!=0))
                 nu[location]=10^18
         elif 'matdamageice' in md.materials.__module__ and damage is not None:
                 print 'computing damage-dependent properties!'
-                Zinv=npy.dot(1-damage[index-1],summation/3.).reshape(-1,)
-                location=npy.nonzero(second_inv)
-                nu[location]=Zinv[location]*B_bar[location]/npy.power(second_inv[location],power[location])
-                location=npy.nonzero(npy.logical_and(second_inv==0,power==0))
+                Zinv=np.dot(1-damage[index-1],summation/3.).reshape(-1,)
+                location=np.nonzero(second_inv)
+                nu[location]=Zinv[location]*B_bar[location]/np.power(second_inv[location],power[location])
+                location=np.nonzero(np.logical_and(second_inv==0,power==0))
                 nu[location]=Zinv[location]*B_bar[location]
                 #clear Zinv
@@ -101 +101 @@
         
         #compute principal properties of stress
-        for i in npy.arange(numberofelements):
+        for i in np.arange(numberofelements):
         
                 #compute stress and strainrate matrices
-                stress=npy.array([ [tau_xx[i], tau_xy[i]], [tau_xy[i], tau_yy[i]] ])
-                strain=npy.array([ [ux[i], uyvx[i]], [uyvx[i], vy[i]] ])
+                stress=np.array([ [tau_xx[i], tau_xy[i]], [tau_xy[i], tau_yy[i]] ])
+                strain=np.array([ [ux[i], uyvx[i]], [uyvx[i], vy[i]] ])
         
                 #eigenvalues and vectors for stress
-                value,directions=npy.linalg.eig(stress);
+                value,directions=np.linalg.eig(stress);
                 idx=abs(value).argsort()[::-1] # sort in descending order
                 value=value[idx]
@@ -116 +116 @@
 
                 #eigenvalues and vectors for strain
-                value,directions=npy.linalg.eig(strain);
+                value,directions=np.linalg.eig(strain);
                 idx=abs(value).argsort()[::-1] # sort in descending order
                 value=value[idx]
@@ -133 +133 @@
         stress.principalvalue2=valuesstress[:,1]
         stress.principalaxis2=directionsstress[:,2:4]
-        stress.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
+        stress.effectivevalue=1./np.sqrt(2.)*np.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
         md.results.stress=stress
         
@@ -144 +144 @@
         strainrate.principalvalue2=valuesstrain[:,1]*md.constants.yts 
         strainrate.principalaxis2=directionsstrain[:,2:4]
-        strainrate.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(strainrate.xx**2+strainrate.yy**2+2.*strainrate.xy**2)
+        strainrate.effectivevalue=1./np.sqrt(2.)*np.sqrt(strainrate.xx**2+strainrate.yy**2+2.*strainrate.xy**2)
         md.results.strainrate=strainrate
         
@@ -155 +155 @@
         deviatoricstress.principalvalue2=valuesstress[:,1]
         deviatoricstress.principalaxis2=directionsstress[:,2:4]
-        deviatoricstress.effectivevalue=1./npy.sqrt(2.)*npy.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
+        deviatoricstress.effectivevalue=1./np.sqrt(2.)*np.sqrt(stress.xx**2+stress.yy**2+2.*stress.xy**2)
         md.results.deviatoricstress=deviatoricstress
 

issm/trunk-jpl/src/m/mech/robintemperature.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from scipy.special import erf
 
@@ -34 +34 @@
         
         #create vertical coordinate variable
-        zstar=npy.sqrt(2.*alphaT*thickness/accumrate)
+        zstar=np.sqrt(2.*alphaT*thickness/accumrate)
         
-        tprofile=surftemp+npy.sqrt(2.*thickness*npy.pi/accumrate/alphaT)*(-heatflux)/2./rho/c*(erf(z/zstar)-erf(thickness/zstar))
+        tprofile=surftemp+np.sqrt(2.*thickness*np.pi/accumrate/alphaT)*(-heatflux)/2./rho/c*(erf(z/zstar)-erf(thickness/zstar))
 
         return tprofile 

issm/trunk-jpl/src/m/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]
         

issm/trunk-jpl/src/m/mech/thomasparams.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from averaging import averaging
 
@@ -75 +75 @@
         
         # checks: any of e1 or e2 equal to zero?
-        pos=npy.nonzero(e1==0)
-        if npy.any(pos==1):
+        pos=np.nonzero(e1==0)
+        if np.any(pos==1):
                 print 'WARNING: first principal strain rate equal to zero.  Value set to 1e-13 s^-1'
                 e1[pos]=1.e-13
-        pos=npy.nonzero(e2==0)
-        if npy.any(pos==1):
+        pos=np.nonzero(e2==0)
+        if np.any(pos==1):
                 print 'WARNING: second principal strain rate equal to zero.  Value set to 1e-13 s^-1'
                 e2[pos]=1.e-13
@@ -89 +89 @@
         
         if coordsys=='principal':
-                b=npy.zeros((md.mesh.numberofvertices,))
+                b=np.zeros((md.mesh.numberofvertices,))
                 ex=e1
                 a=e2/e1
-                pos=npy.nonzero(npy.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
+                pos=np.nonzero(np.logical_and(e1<0,e2>0)) # longitudinal compression and lateral tension
                 a[pos]=e1[pos]/e2[pos]
                 ex[pos]=e2[pos]
-                pos2=npy.nonzero(e1<0 & e2<0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal compression
+                pos2=np.nonzero(e1<0 & e2<0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal compression
                 a[pos2]=e1[pos2]/e2[pos2]
                 ex[pos2]=e2[pos2]
-                pos3=npy.nonzero(e1>0 & e2>0 & npy.abs(e1)<npy.abs(e2)) # lateral and longitudinal tension
+                pos3=np.nonzero(e1>0 & e2>0 & np.abs(e1)<np.abs(e2)) # lateral and longitudinal tension
                 a[pos3]=e1[pos3]/e2[pos3]
                 ex[pos3]=e2[pos3]
-                ind=npy.nonzero(e1<0 & e2<0)
+                ind=np.nonzero(e1<0 & e2<0)
                 a[ind]=-a[ind] # where both strain rates are compressive, enforce negative alpha
-                sigxx=(npy.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
+                sigxx=(np.abs(ex)/((1.+a+a**2)**((n-1.)/2.)))**(1./n)*B
         elif coordsys=='xy':
                 ex=exx
@@ -110 +110 @@
         elif coordsys=='longitudinal':
                 # using longitudinal strain rates defined by observed velocity vector
-                velangle=npy.arctan(md.initialization.vy/md.initialization.vx)
-                pos=npy.nonzero(md.initialization.vx==0)
-                velangle[pos]=npy.pi/2
-                ex=0.5*(exx+eyy)+0.5*(exx-eyy)*npy.cos(2.*velangle)+exy*npy.sin(2.*velangle)
+                velangle=np.arctan(md.initialization.vy/md.initialization.vx)
+                pos=np.nonzero(md.initialization.vx==0)
+                velangle[pos]=np.pi/2
+                ex=0.5*(exx+eyy)+0.5*(exx-eyy)*np.cos(2.*velangle)+exy*np.sin(2.*velangle)
                 ey=exx+eyy-ex # trace of strain rate tensor is invariant
-                exy=-0.5*(exx-eyy)*npy.sin(2.*velangle)+exy*npy.cos(2.*velangle)
+                exy=-0.5*(exx-eyy)*np.sin(2.*velangle)+exy*np.cos(2.*velangle)
                 a=ey/ex
                 b=exy/ex
@@ -124 +124 @@
         # a < -1 in areas of strong lateral compression or longitudinal compression and
         # theta flips sign at a = -2
-        pos=npy.nonzero(npy.abs((npy.abs(a)-2.))<1.e-3)
+        pos=np.nonzero(np.abs((np.abs(a)-2.))<1.e-3)
         if len(pos)>0:
                 print 'Warning: ', len(pos), ' vertices have alpha within 1e-3 of -2'
@@ -131 +131 @@
         if eq=='Weertman1D':
                 theta=1./8
-                a=npy.zeros((md.mesh.numberofvertices,))
+                a=np.zeros((md.mesh.numberofvertices,))
         elif eq=='Weertman2D':
                 theta=1./9
-                a=npy.ones((md.mesh.numberofvertices,))
+                a=np.ones((md.mesh.numberofvertices,))
         elif eq=='Thomas':
-                theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(npy.abs(2.+a)**n)
+                theta=((1.+a+a**2+b**2)**((n-1.)/2.))/(np.abs(2.+a)**n)
         else:
                 raise ValueError('argument passed to "eq" not valid')

issm/trunk-jpl/src/m/plot/applyoptions.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from cmaptools import truncate_colormap
 from plot_contour import plot_contour
@@ -223 +223 @@
                         cb.locator=MaxNLocator(nbins=5) # default 5 ticks
                 if options.exist('colorbartickspacing'):
-                        locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing'))
+                        locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing'))
                         cb.set_ticks(locs)
                 if options.exist('colorbarlines'):
-                        locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines'))
-                        cb.add_lines(locs,['k' for i in range(len(locs))],npy.ones_like(locs))
+                        locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines'))
+                        cb.add_lines(locs,['k' for i in range(len(locs))],np.ones_like(locs))
                 if options.exist('colorbarlineatvalue'):
                         locs=options.getfieldvalue('colorbarlineatvalue')
                         colors=options.getfieldvalue('colorbarlineatvaluecolor',['k' for i in range (len(locs))])
-                        widths=options.getfieldvalue('colorbarlineatvaluewidth',npy.ones_like(locs))
+                        widths=options.getfieldvalue('colorbarlineatvaluewidth',np.ones_like(locs))
                         cb.add_lines(locs,colors,widths)
                 if options.exist('colorbartitle'):

issm/trunk-jpl/src/m/plot/colormaps/cmaptools.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 
 try:
@@ -21 +21 @@
 
         new_cmap = mpl.colors.LinearSegmentedColormap.from_list('trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name,
-                a=minval, b=maxval), cmap(npy.linspace(minval, maxval, n)))
+                a=minval, b=maxval), cmap(np.linspace(minval, maxval, n)))
         
         return new_cmap

issm/trunk-jpl/src/m/plot/plot_contour.py

@@ -26 +26 @@
                 pass
         elif datatype==3: # quiver (vector) data
-                data=npy.sqrt(datain**2)
+                data=np.sqrt(datain**2)
         else:
                 raise ValueError('datatype not supported in call to plot_contour')

issm/trunk-jpl/src/m/plot/plot_overlay.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from processmesh import processmesh
 from processdata import processdata
@@ -28 +28 @@
         if data=='none' or data==None:
                 imageonly=1
-                data=npy.float('nan')*npy.ones((md.mesh.numberofvertices,))
+                data=np.float('nan')*np.ones((md.mesh.numberofvertices,))
                 datatype=1
         else:
@@ -70 +70 @@
 
         # normalize array
-        arr=arr/npy.float(npy.max(arr.ravel()))
+        arr=arr/np.float(np.max(arr.ravel()))
         arr=1.-arr # somehow the values got flipped
 
@@ -96 +96 @@
         dy=trans[5]     
         
-        xarr=npy.arange(xmin,xmax,dx)
-        yarr=npy.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is)
-        xg,yg=npy.meshgrid(xarr,yarr)
+        xarr=np.arange(xmin,xmax,dx)
+        yarr=np.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is)
+        xg,yg=np.meshgrid(xarr,yarr)
         overlaylims=options.getfieldvalue('overlaylims',[min(arr.ravel()),max(arr.ravel())])
         norm=mpl.colors.Normalize(vmin=overlaylims[0],vmax=overlaylims[1])
@@ -118 +118 @@
                 #test
                 #m.ax=ax
-                meridians=npy.arange(-180.,181.,1.)
-                parallels=npy.arange(-80.,80.,1.)
+                meridians=np.arange(-180.,181.,1.)
+                parallels=np.arange(-80.,80.,1.)
                 m.drawparallels(parallels,labels=[0,0,1,1]) # labels=[left,right,top,bottom]
                 m.drawmeridians(meridians,labels=[1,1,0,0])
                 m.drawcoastlines()
-                pc=m.pcolormesh(xg, yg, npy.flipud(arr), cmap=mpl.cm.Greys, norm=norm, ax=ax)
+                pc=m.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm, ax=ax)
 
         else:
-                pc=ax.pcolormesh(xg, yg, npy.flipud(arr), cmap=mpl.cm.Greys, norm=norm)
+                pc=ax.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm)
         
         #rasterization? 

issm/trunk-jpl/src/m/plot/plot_streamlines.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from processmesh import processmesh
 from processdata import processdata
@@ -45 +45 @@
 
     # format data for matplotlib streamplot function
-    yg,xg=npy.mgrid[min(md.mesh.y):max(md.mesh.y):100j,min(md.mesh.x):max(md.mesh.x):100j]
+    yg,xg=np.mgrid[min(md.mesh.y):max(md.mesh.y):100j,min(md.mesh.x):max(md.mesh.x):100j]
     ug=griddata((x,y),u,(xg,yg),method='linear')
     vg=griddata((x,y),v,(xg,yg),method='linear')
@@ -60 +60 @@
     if linewidth=='vel':
         scale=options.getfieldvalue('streamlineswidthscale',3)
-        vel=npy.sqrt(ug**2+vg**2)
-        linewidth=scale*vel/npy.amax(vel)
+        vel=np.sqrt(ug**2+vg**2)
+        linewidth=scale*vel/np.amax(vel)
         linewidth[linewidth<0.5]=0.5
 

issm/trunk-jpl/src/m/plot/plotmodel.py

@@ -1 @@
-import numpy as npy
+import numpy as np
 from plotoptions import plotoptions
 
@@ -34 +34 @@
                 nr=True
         else:
-                nrows=npy.ceil(numberofplots/subplotwidth)
+                nrows=np.ceil(numberofplots/subplotwidth)
                 nr=False
         

issm/trunk-jpl/src/m/plot/processmesh.py

@@ -1 +1 @@
 from math import isnan
-import numpy as npy
+import numpy as np
 
 def processmesh(md,data,options):
@@ -35 +35 @@
                         z=md.mesh.z
                 except AttributeError:
-                        z=npy.zeros_like(md.mesh.x)
+                        z=np.zeros_like(md.mesh.x)
                 
                 try:

issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py

@@ -3 +3 @@
 from collections import OrderedDict
 import results as resultsclass
-import MatlabFuncs as m
 
 def parseresultsfromdisk(md,filename,iosplit):
         if iosplit:
-                results=parseresultsfromdiskiosplit(md,filename)
+                saveres=parseresultsfromdiskiosplit(md,filename)
         else:
-                results=parseresultsfromdiskioserial(md,filename)
-
-        return results
+                saveres=parseresultsfromdiskioserial(md,filename)
+
+        return saveres
 
 def parseresultsfromdiskioserial(md,filename):    # {{{
-
         #Open file
         try:
@@ -22 +20 @@
 
         #initialize results: 
-        results=[]
-        results.append(None)
+        saveres=[]
 
         #Read fields until the end of the file.
-        result=ReadData(fid,md)
+        loadres=ReadData(fid,md)
 
         counter=0
         check_nomoresteps=0
-        step=result['step']
-
-        while result:
-
+        step=loadres['step']
+        print ('step is now {}'.format(step))
+
+        while loadres:
+                #check that the new result does not add a step, which would be an error: 
                 if check_nomoresteps:
-                        #check that the new result does not add a step, which would be an error: 
-                        if result['step']>=1:
+                        if loadres['step']>=1:
                                 raise TypeError("parsing results for a steady-state core, which incorporates transient results!")
 
                 #Check step, increase counter if this is a new step
-                if(step!=result['step'] and result['step']>1):
+                if(step!=loadres['step'] and loadres['step']>1):
                         counter = counter + 1
-                        step    = result['step']
-
-                #Add result
-                if result['step']==0:
+                        step    = loadres['step']
+
+                #Add result
+                if loadres['step']==0:
                         #if we have a step = 0, this is a steady state solution, don't expect more steps. 
                         index = 0;
                         check_nomoresteps=1
-        
-                elif result['step']==1:
+                elif loadres['step']==1:
                         index = 0
                 else:
                         index = counter;
-        
-                if index > len(results)-1:
-                        for i in xrange(len(results)-1,index-1):
-                                results.append(None)
-                        results.append(resultsclass.results())
                 
-                elif results[index] is None:
-                        results[index]=resultsclass.results()
-
+                if index > len(saveres)-1:
+                        for i in xrange(len(saveres)-1,index-1):
+                                saveres.append(None)
+                        saveres.append(resultsclass.results())
+                elif saveres[index] is None:
+                        saveres[index]=resultsclass.results()
                         
                 #Get time and step
-                if result['step'] != -9999.:
-                        setattr(results[index],'step',result['step'])
-                if result['time'] != -9999.:
-                        setattr(results[index],'time',result['time']) 
-        
-                #Add result
-                if hasattr(results[index],result['fieldname']) and not m.strcmp(result['fieldname'],'SolutionType'):
-                        setattr(results[index],result['fieldname'],numpy.vstack((getattr(results[index],result['fieldname']),result['field'])))
-                else:
-                        setattr(results[index],result['fieldname'],result['field'])
-
-                #read next result
-                result=ReadData(fid,md)
+                if loadres['step'] != -9999.:
+                        saveres[index].__dict__['step']=loadres['step']
+                if loadres['time'] != -9999.:
+                        saveres[index].__dict__['time']=loadres['time']
+
+                #Add result
+                saveres[index].__dict__[loadres['fieldname']]=loadres['field']
+
+                #read next result
+                loadres=ReadData(fid,md)
 
         fid.close()
 
-        return results
+        return saveres
         # }}}
 def parseresultsfromdiskiosplit(md,filename):    # {{{
@@ -91 +82 @@
                 raise IOError("loadresultsfromdisk error message: could not open '%s' for binary reading." % filename)
 
-        results=[]
+        saveres=[]
 
         #if we have done split I/O, ie, we have results that are fragmented across patches, 
         #do a first pass, and figure out the structure of results
-        result=ReadDataDimensions(fid)
-        while result:
+        loadres=ReadDataDimensions(fid)
+        while loadres:
 
                 #Get time and step
-                if result['step'] > len(results):
-                        for i in xrange(len(results),result['step']-1):
-                                results.append(None)
-                        results.append(resultsclass.results())
-                setattr(results[result['step']-1],'step',result['step'])
-                setattr(results[result['step']-1],'time',result['time']) 
-
-                #Add result
-                setattr(results[result['step']-1],result['fieldname'],float('NaN'))
-
-                #read next result
-                result=ReadDataDimensions(fid)
+                if loadres['step'] > len(saveres):
+                        for i in xrange(len(saveres),loadres['step']-1):
+                                saveres.append(None)
+                        saveres.append(resultsclass.results())
+                setattr(saveres[loadres['step']-1],'step',loadres['step'])
+                setattr(saveres[loadres['step']-1],'time',loadres['time']) 
+
+                #Add result
+                setattr(saveres[loadres['step']-1],loadres['fieldname'],float('NaN'))
+
+                #read next result
+                loadres=ReadDataDimensions(fid)
 
         #do a second pass, and figure out the size of the patches
         fid.seek(0)    #rewind
-        result=ReadDataDimensions(fid)
-        while result:
-
-                #read next result
-                result=ReadDataDimensions(fid)
+        loadres=ReadDataDimensions(fid)
+        while loadres:
+
+                #read next result
+                loadres=ReadDataDimensions(fid)
 
         #third pass, this time to read the real information
         fid.seek(0)    #rewind
-        result=ReadData(fid,md)
-        while result:
+        loadres=ReadData(fid,md)
+        while loadres:
 
                 #Get time and step
-                if result['step']> len(results):
-                        for i in xrange(len(results),result['step']-1):
-                                results.append(None)
-                        results.append(resultsclass.results())
-                setattr(results[result['step']-1],'step',result['step'])
-                setattr(results[result['step']-1],'time',result['time']) 
-
-                #Add result
-                setattr(results[result['step']-1],result['fieldname'],result['field'])
-
-                #read next result
-                result=ReadData(fid,md)
+                if loadres['step']> len(saveres):
+                        for i in xrange(len(saveres),loadres['step']-1):
+                                saveres.append(None)
+                        saveres.append(saveresclass.saveres())
+                setattr(saveres[loadres['step']-1],'step',loadres['step'])
+                setattr(saveres[loadres['step']-1],'time',loadres['time']) 
+
+                #Add result
+                setattr(saveres[loadres['step']-1],loadres['fieldname'],loadres['field'])
+
+                #read next result
+                loadres=ReadData(fid,md)
 
         #close file
         fid.close()
 
-        return results
+        return saveres
         # }}}
 def ReadData(fid,md):    # {{{
@@ -177 +168 @@
                 #Process units here FIXME: this should not be done here!
                 yts=md.constants.yts
-                if m.strcmp(fieldname,'BalancethicknessThickeningRate'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'Time'):
+                if fieldname=='BalancethicknessThickeningRate':
+                        field = field*yts
+                elif fieldname=='Time':
                         field = field/yts
-                elif m.strcmp(fieldname,'HydrologyWaterVx'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'HydrologyWaterVy'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'Vx'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'Vy'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'Vz'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'Vel'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'TotalFloatingBmb'):
+                elif fieldname=='HydrologyWaterVx':
+                        field = field*yts
+                elif fieldname=='HydrologyWaterVy':
+                        field = field*yts
+                elif fieldname=='Vx':
+                        field = field*yts
+                elif fieldname=='Vy':
+                        field = field*yts
+                elif fieldname=='Vz':
+                        field = field*yts
+                elif fieldname=='Vel':
+                        field = field*yts
+                elif fieldname=='BasalforcingsGroundediceMeltingRate':
+                        field = field*yts
+                elif fieldname=='TotalFloatingBmb':
                         field = field/10.**12.*yts #(GigaTon/year)
-                elif m.strcmp(fieldname,'TotalGroundedBmb'):
+                elif fieldname=='TotalGroundedBmb':
                         field = field/10.**12.*yts #(GigaTon/year)
-                elif m.strcmp(fieldname,'TotalSmb'):
+                elif fieldname=='TotalSmb':
                         field = field/10.**12.*yts #(GigaTon/year)
-                elif m.strcmp(fieldname,'SmbMassBalance'):
-                        field = field*yts
-                elif m.strcmp(fieldname,'CalvingCalvingrate'):
-                        field = field*yts
-
-                result=OrderedDict()
-                result['fieldname']=fieldname
-                result['time']=time
-                result['step']=step
-                result['field']=field
+                elif fieldname=='SmbMassBalance':
+                        field = field*yts
+                elif fieldname=='CalvingCalvingrate':
+                        field = field*yts
+
+                saveres=OrderedDict()
+                saveres['fieldname']=fieldname
+                saveres['time']=time
+                saveres['step']=step
+                saveres['field']=field
 
         except struct.error as e:
-                result=None
-
-        return result
+                saveres=None
+
+        return saveres
         # }}}
 def ReadDataDimensions(fid):    # {{{
@@ -246 +237 @@
                         raise TypeError("cannot read data of type %d" % type)
 
-                result=OrderedDict()
-                result['fieldname']=fieldname
-                result['time']=time
-                result['step']=step
-                result['M']=M
-                result['N']=N
+                saveres=OrderedDict()
+                saveres['fieldname']=fieldname
+                saveres['time']=time
+                saveres['step']=step
+                saveres['M']=M
+                saveres['N']=N
 
         except struct.error as e:
-                result=None
-
-        return result
-        # }}}
+                saveres=None
+
+        return saveres
+        # }}}