Index: issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py
===================================================================
--- issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py	(revision 21254)
+++ issm/trunk-jpl/src/m/boundaryconditions/PattynSMB.py	(revision 21303)
@@ -1,4 +1,4 @@
 import os
-import numpy as np
+import numpy as  np
 def PattynSMB(md,Tf):
 	"""
Index: issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py
===================================================================
--- issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py	(revision 21254)
+++ issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py	(revision 21303)
@@ -1,4 +1,4 @@
 import os
-import numpy
+import numpy as np
 from ContourToMesh import ContourToMesh
 
@@ -14,16 +14,16 @@
 
 	#node on Dirichlet
-	pos=numpy.nonzero(md.mesh.vertexonboundary)
-	md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices)
-	md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices)
-	md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices)
+	pos=np.nonzero(md.mesh.vertexonboundary)
+	md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices)
+	md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices)
+	md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices)
 	md.stressbalance.spcvx[pos]=0
 	md.stressbalance.spcvy[pos]=0
 	md.stressbalance.spcvz[pos]=0
-	md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
-	md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
+	md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6))
+	md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3))
 
 	#Dirichlet Values
-	if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
+	if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
 		print "      boundary conditions for stressbalance model: spc set as observed velocities"
 		md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos]
@@ -39,18 +39,18 @@
 
 	#Deal with other boundary conditions
-	if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
+	if np.all(np.isnan(md.balancethickness.thickening_rate)):
+		md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices))
 		print "      no balancethickness.thickening_rate specified: values set as zero"
-	md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-	md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-	md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+	md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+	md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices))
 
-	if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-		md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
+		md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices))
 		if hasattr(md.mesh,'vertexonsurface'):
-			pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+			pos=np.nonzero(md.mesh.vertexonsurface)[0]
 			md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
-		if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
-			md.basalforcings.geothermalflux=50.*10**-3*numpy.ones((md.mesh.numberofvertices,1))    #50 mW/m^2
+		if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
+			md.basalforcings.geothermalflux=50.*10**-3*np.ones((md.mesh.numberofvertices))    #50 mW/m^2
 	else:
 		print "      no thermal boundary conditions created: no observed temperature found"
Index: issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py
===================================================================
--- issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py	(revision 21254)
+++ issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py	(revision 21303)
@@ -1,4 +1,4 @@
 import os
-import numpy
+import numpy as np
 from ContourToMesh import ContourToMesh
 import MatlabFuncs as m
@@ -27,18 +27,18 @@
 			raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
 		nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
-		nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
+		nodeonicefront=np.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
 	else:
-		nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool)
+		nodeonicefront=np.zeros((md.mesh.numberofvertices),bool)
 
 #	pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
-	pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(nodeonicefront)))[0]
-	md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices)
-	md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices)
-	md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices)
-	md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
-	md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
+	pos=np.nonzero(np.logical_and(md.mesh.vertexonboundary,np.logical_not(nodeonicefront)))[0]
+	md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices)
+	md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices)
+	md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices)
+	md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6))
+	md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3))
 
 	#Icefront position
-	pos=numpy.nonzero(nodeonicefront)[0]
+	pos=np.nonzero(nodeonicefront)[0]
 	md.mask.ice_levelset[pos]=0
 
@@ -53,10 +53,10 @@
 		values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1]
 		segmentsfront=1-values
-		numpy.sum(segmentsfront,axis=1)!=numbernodesfront
-		segments=numpy.nonzero(numpy.sum(segmentsfront,axis=1)!=numbernodesfront)[0]
+		np.sum(segmentsfront,axis=1)!=numbernodesfront
+		segments=np.nonzero(np.sum(segmentsfront,axis=1)!=numbernodesfront)[0]
 		#Find all nodes for these segments and spc them
 		pos=md.mesh.segments[segments,0:-1]-1
 	else:
-		pos=numpy.nonzero(md.mesh.vertexonboundary)[0]
+		pos=np.nonzero(md.mesh.vertexonboundary)[0]
 	md.stressbalance.spcvx[pos]=0
 	md.stressbalance.spcvy[pos]=0
@@ -64,10 +64,10 @@
 																													   
 	#Dirichlet Values
-	if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
+	if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
 		#reshape to rank-2 if necessary to match spc arrays
-		if numpy.ndim(md.inversion.vx_obs)==1:
-			md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1,1)
-		if numpy.ndim(md.inversion.vy_obs)==1:
-			md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1,1)
+		if np.ndim(md.inversion.vx_obs)==1:
+			md.inversion.vx_obs=md.inversion.vx_obs.reshape(-1,)
+		if np.ndim(md.inversion.vy_obs)==1:
+			md.inversion.vy_obs=md.inversion.vy_obs.reshape(-1,)
 		print "      boundary conditions for stressbalance model: spc set as observed velocities"
 		md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos]
@@ -81,18 +81,18 @@
 
 	#Deal with other boundary conditions
-	if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
+	if np.all(np.isnan(md.balancethickness.thickening_rate)):
+		md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices))
 		print "      no balancethickness.thickening_rate specified: values set as zero"
-	md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-	md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-	md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+	md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+	md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices))
 
-	if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-		md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
+		md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices))
 		if hasattr(md.mesh,'vertexonsurface'):
-			pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+			pos=np.nonzero(md.mesh.vertexonsurface)[0]
 			md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
-		if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
-			md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
+		if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
+			md.basalforcings.geothermalflux=np.zeros((md.mesh.numberofvertices))
 	else:
 		print "      no thermal boundary conditions created: no observed temperature found"
Index: issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py
===================================================================
--- issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 21254)
+++ issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 21303)
@@ -1,4 +1,4 @@
 import os
-import numpy
+import numpy as np
 from ContourToMesh import ContourToMesh
 import MatlabFuncs as m
@@ -29,25 +29,25 @@
 			raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile)
 		incontour=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
-		vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1))
+		vertexonicefront=np.logical_and(md.mesh.vertexonboundary,incontour.reshape(-1))
 	else:
 		#Guess where the ice front is
-		vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices,1))
-		pos=numpy.nonzero(numpy.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0]
+		vertexonfloatingice=np.zeros((md.mesh.numberofvertices))
+		pos=np.nonzero(np.sum(md.mask.groundedice_levelset[md.mesh.elements-1]<0.,axis=1) >0.)[0]
 		vertexonfloatingice[md.mesh.elements[pos].astype(int)-1]=1.
-		vertexonicefront=numpy.logical_and(numpy.reshape(md.mesh.vertexonboundary,(-1,1)),vertexonfloatingice>0.)
+		vertexonicefront=np.logical_and(np.reshape(md.mesh.vertexonboundary,(-1,)),vertexonfloatingice>0.)
 
 #	pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
-	pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonboundary,numpy.logical_not(vertexonicefront)))[0]
-	if not numpy.size(pos):
+	pos=np.nonzero(np.logical_and(md.mesh.vertexonboundary,np.logical_not(vertexonicefront)))[0]
+	if not np.size(pos):
 		print "SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually."
 
-	md.stressbalance.spcvx=float('nan')*numpy.ones(md.mesh.numberofvertices)
-	md.stressbalance.spcvy=float('nan')*numpy.ones(md.mesh.numberofvertices)
-	md.stressbalance.spcvz=float('nan')*numpy.ones(md.mesh.numberofvertices)
-	md.stressbalance.referential=float('nan')*numpy.ones((md.mesh.numberofvertices,6))
-	md.stressbalance.loadingforce=0*numpy.ones((md.mesh.numberofvertices,3))
+	md.stressbalance.spcvx=float('nan')*np.ones(md.mesh.numberofvertices)
+	md.stressbalance.spcvy=float('nan')*np.ones(md.mesh.numberofvertices)
+	md.stressbalance.spcvz=float('nan')*np.ones(md.mesh.numberofvertices)
+	md.stressbalance.referential=float('nan')*np.ones((md.mesh.numberofvertices,6))
+	md.stressbalance.loadingforce=0*np.ones((md.mesh.numberofvertices,3))
 
 	#Position of ice front
-	pos=numpy.nonzero(vertexonicefront)[0]
+	pos=np.nonzero(vertexonicefront)[0]
 	md.mask.ice_levelset[pos]=0
 
@@ -62,10 +62,10 @@
 		values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1]
 		segmentsfront=1-values
-		numpy.sum(segmentsfront,axis=1)!=numbernodesfront
-		segments=numpy.nonzero(numpy.sum(segmentsfront,axis=1)!=numbernodesfront)[0]
+		np.sum(segmentsfront,axis=1)!=numbernodesfront
+		segments=np.nonzero(np.sum(segmentsfront,axis=1)!=numbernodesfront)[0]
 		#Find all nodes for these segments and spc them
 		pos=md.mesh.segments[segments,0:-1]-1
 	else:
-		pos=numpy.nonzero(md.mesh.vertexonboundary)[0]
+		pos=np.nonzero(md.mesh.vertexonboundary)[0]
 	md.stressbalance.spcvx[pos]=0
 	md.stressbalance.spcvy[pos]=0
@@ -73,5 +73,5 @@
 
 	#Dirichlet Values
-	if isinstance(md.inversion.vx_obs,numpy.ndarray) and numpy.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,numpy.ndarray) and numpy.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
+	if isinstance(md.inversion.vx_obs,np.ndarray) and np.size(md.inversion.vx_obs,axis=0)==md.mesh.numberofvertices and isinstance(md.inversion.vy_obs,np.ndarray) and np.size(md.inversion.vy_obs,axis=0)==md.mesh.numberofvertices:
 		print "      boundary conditions for stressbalance model: spc set as observed velocities"
 		md.stressbalance.spcvx[pos]=md.inversion.vx_obs[pos]
@@ -80,6 +80,6 @@
 		print "      boundary conditions for stressbalance model: spc set as zero"
 
-	md.hydrology.spcwatercolumn=numpy.zeros((md.mesh.numberofvertices,2))
-	pos=numpy.nonzero(md.mesh.vertexonboundary)[0]
+	md.hydrology.spcwatercolumn=np.zeros((md.mesh.numberofvertices,2))
+	pos=np.nonzero(md.mesh.vertexonboundary)[0]
 	md.hydrology.spcwatercolumn[pos,0]=1
 
@@ -89,20 +89,20 @@
 
 	#Deal with other boundary conditions
-	if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
-		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
+	if np.all(np.isnan(md.balancethickness.thickening_rate)):
+		md.balancethickness.thickening_rate=np.zeros((md.mesh.numberofvertices))
 		print "      no balancethickness.thickening_rate specified: values set as zero"
 
-	md.masstransport.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-	md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
-	md.damage.spcdamage=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.masstransport.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+	md.balancethickness.spcthickness=float('nan')*np.ones((md.mesh.numberofvertices))
+	md.damage.spcdamage=float('nan')*np.ones((md.mesh.numberofvertices))
 
-	if isinstance(md.initialization.temperature,numpy.ndarray) and numpy.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
-		md.thermal.spctemperature=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	if isinstance(md.initialization.temperature,np.ndarray) and np.size(md.initialization.temperature,axis=0)==md.mesh.numberofvertices:
+		md.thermal.spctemperature=float('nan')*np.ones((md.mesh.numberofvertices))
 		if hasattr(md.mesh,'vertexonsurface'):
-			pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+			pos=np.nonzero(md.mesh.vertexonsurface)[0]
 			md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
-		if not isinstance(md.basalforcings.geothermalflux,numpy.ndarray) or not numpy.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
-			md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1))
-			md.basalforcings.geothermalflux[numpy.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3    #50mW/m2
+		if not isinstance(md.basalforcings.geothermalflux,np.ndarray) or not np.size(md.basalforcings.geothermalflux,axis=0)==md.mesh.numberofvertices:
+			md.basalforcings.geothermalflux=np.zeros((md.mesh.numberofvertices))
+			md.basalforcings.geothermalflux[np.nonzero(md.mask.groundedice_levelset>0.)]=50.*10.**-3    #50mW/m2
 	else:
 		print "      no thermal boundary conditions created: no observed temperature found"
Index: issm/trunk-jpl/src/m/boundaryconditions/love_numbers.py
===================================================================
--- issm/trunk-jpl/src/m/boundaryconditions/love_numbers.py	(revision 21254)
+++ issm/trunk-jpl/src/m/boundaryconditions/love_numbers.py	(revision 21303)
@@ -1,5 +1,5 @@
 from MatlabFuncs import *
 from model import *
-from numpy import *
+from np.import *
 
 def love_numbers(value,*varargin):
