Index: issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py
===================================================================
--- issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py	(revision 12852)
+++ issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py	(revision 12854)
@@ -1,44 +1,47 @@
-from numpy import *
+import os
+import numpy
+
 def SetIceShelfBC(md,icefrontfile=''):
-	#SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
-	#
-	#   Neumann BC are used on the ice front (an ANRGUS contour around the ice front
-	#   must be given in input)
-	#   Dirichlet BC are used elsewhere for diagnostic
-	#
-	#   Usage:
-	#      md=SetIceShelfBC(md,varargin)
-	#
-	#   Example:
-	#      md=SetIceShelfBC(md);
-	#      md=SetIceShelfBC(md,'Front.exp');
-	#
-	#   See also: SETICESHEETBC, SETMARINEICESHEETBC
+	"""
+	SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
+	 
+	    Neumann BC are used on the ice front (an ANRGUS contour around the ice front
+	    must be given in input)
+	    Dirichlet BC are used elsewhere for diagnostic
+	 
+	    Usage:
+	       md=SetIceShelfBC(md,varargin)
+	 
+	    Example:
+	       md=SetIceShelfBC(md);
+	       md=SetIceShelfBC(md,'Front.exp');
+	 
+	    See also: SETICESHEETBC, SETMARINEICESHEETBC
+	"""
 
 	#node on Dirichlet (boundary and ~icefront)
-	if not icefrontfile:
-		nodeonicefront=zeros(md.mesh.numberofvertices)
+	if icefrontfile:
+		if not os.path.exists(icefrontfile):
+			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=double(md.mesh.vertexonboundary and nodeinsideicefront)
 	else:
-		if not os.path.isfile(icefrontfile):
-			print 'SetIceShelfBC error message: ice front file '+icefrontfile+ ' not found'
-			return []
-		nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2)
-		nodeonicefront=double(md.mesh.vertexonboundary.astype(bool) and nodeinsideicefront.astype(bool))
-	
-	pos=argwhere(logical_and(md.mesh.vertexonboundary.astype(bool),~nodeonicefront.astype(bool)))
-	md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices);
-	md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices);
-	md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices);
-	md.diagnostic.spcvx[pos]=0;
-	md.diagnostic.spcvy[pos]=0;
-	md.diagnostic.spcvz[pos]=0;
-	md.diagnostic.referential=float('NaN')*ones((md.mesh.numberofvertices,6),float);
+		nodeonicefront=numpy.zeros(md.mesh.numberofvertices)
+
+#	pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
+	pos=[i for i,(vob,noif) in enumerate(zip(md.mesh.vertexonboundary,nodeonicefront)) if vob and not noif]
+	md.diagnostic.spcvx=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+	md.diagnostic.spcvy=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+	md.diagnostic.spcvz=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+	md.diagnostic.spcvx[pos]=0
+	md.diagnostic.spcvy[pos]=0
+	md.diagnostic.spcvz[pos]=0
+	md.diagnostic.referential=float('NaN')*numpy.ones((md.mesh.numberofvertices,6))
 
 	#Dirichlet Values
-	if ~isnan(md.inversion.vx_obs) and ~isnan(md.inversion.vy_obs):
-		if (len(md.inversion.vx_obs)==md.mesh.numberofvertices) and (len(md.inversion.vy_obs)==md.mesh.numberofvertices):
-			print '      boundary conditions for diagnostic model: spc set as observed velocities'
-			md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos]
-			md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos]
+	if numpy.size(md.inversion.vx_obs)==md.mesh.numberofvertices and numpy.size(md.inversion.vy_obs)==md.mesh.numberofvertices:
+		print '      boundary conditions for diagnostic model: spc set as observed velocities'
+		md.diagnostic.spcvx[pos]=md.inversion.vx_obs[pos]
+		md.diagnostic.spcvy[pos]=md.inversion.vy_obs[pos]
 	else:
 		print '      boundary conditions for diagnostic model: spc set as zero'
@@ -46,46 +49,50 @@
 	#segment on Ice Front
 	#segment on Neumann (Ice Front)
-	segs1=md.mesh.segments[:,0].astype(int)-1
-	segs2=md.mesh.segments[:,1].astype(int)-1
-
-	pos=argwhere(logical_and(nodeonicefront[segs1].astype(bool),nodeonicefront[segs2].astype(bool)))
-	if (md.mesh.dimension==2):
-		pressureload=md.mesh.segments[pos,:];
+#	pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
+	pos=[i for i,(noif1,noif2) in enumerate(zip(nodeonicefront[md.mesh.segments[:,0].astype('int')-1],nodeonicefront[md.mesh.segments[:,1].astype('int')-1])) if noif1 or noif2]
+	if   md.mesh.dimension==2:
+		pressureload=md.mesh.segments[pos,:]
 	elif md.mesh.dimension==3:
-		pressureload_layer=[md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d, md.mesh.segments[pos,2]];
-		pressureload=[];
-		for i in range(1,md.mesh.numberoflayers-1):
-			pressureload=[pressureload ,pressureload_layer1[:,1:4]+(i-1)*md.mesh.numberofvertices2d, pressureload_layer1[:,5]+(i-1)*md.mesh.numberofelements2d ];
+#		pressureload_layer1=[md.mesh.segments(pos,1:2)  md.mesh.segments(pos,2)+md.mesh.numberofvertices2d  md.mesh.segments(pos,1)+md.mesh.numberofvertices2d  md.mesh.segments(pos,3)];
+		pressureload_layer1=numpy.concatenate((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]),axis=1)
+		pressureload=numpy.zeros((0,5))
+		for i in xrange(1,md.mesh.numberoflayers):
+#			pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
+			pressureload=numpy.concatenate((pressureload,numpy.concatenate((pressureload_layer1[:,0:3]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d),axis=1)),axis=0)
 
 	#Add water or air enum depending on the element
-	pressureload=[pressureload, 1*md.mask.elementonfloatingice(pressureload[:,end])];
+#	pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
+	pressureload=numpy.concatenate((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape((-1,1))),axis=1)
 
 	#plug onto model
-	md.diagnostic.icefront=pressureload;
+	md.diagnostic.icefront=pressureload
 
 	#Create zeros basalforcings and surfaceforcings
-	if isnan(md.surfaceforcings.precipitation):
-		md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices);
+	if numpy.isnan(md.surfaceforcings.precipitation).all():
+		md.surfaceforcings.precipitation=numpy.zeros(md.mesh.numberofvertices)
 		print '      no surfaceforcings.precipitation specified: values set as zero'
-	if isnan(md.surfaceforcings.mass_balance):
-		md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices);
+	if numpy.isnan(md.surfaceforcings.mass_balance).all():
+		md.surfaceforcings.mass_balance=numpy.zeros(md.mesh.numberofvertices)
 		print '      no surfaceforcings.mass_balance specified: values set as zero'
-	if isnan(md.basalforcings.melting_rate):
-		md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices)
+	if numpy.isnan(md.basalforcings.melting_rate).all():
+		md.basalforcings.melting_rate=numpy.zeros(md.mesh.numberofvertices)
 		print '      no basalforcings.melting_rate specified: values set as zero'
-	if isnan(md.balancethickness.thickening_rate):
-		md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices);
+	if numpy.isnan(md.balancethickness.thickening_rate).all():
+		md.balancethickness.thickening_rate=numpy.zeros(md.mesh.numberofvertices)
 		print '      no balancethickness.thickening_rate specified: values set as zero'
 
-	md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices);
-	md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices);
+	md.prognostic.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+	md.balancethickness.spcthickness=float('NaN')*numpy.ones(md.mesh.numberofvertices)
 
-	if (len(md.initialization.temperature)==md.mesh.numberofvertices):
-		md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices);
-		pos=argwhere(md.mesh.vertexonsurface); md.thermal.spctemperature[pos]=md.initialization.temperature[pos]; #impose observed temperature on surface
-		if (len(md.basalforcings.geothermalflux) !=md.mesh.numberofvertices):
-			md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices);
+	if numpy.size(md.initialization.temperature)==md.mesh.numberofvertices:
+		md.thermal.spctemperature=float('NaN')*numpy.ones(md.mesh.numberofvertices)
+#		pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
+		pos=[i for i,vos in enumerate(md.mesh.vertexonsurface) if vos]
+		md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    # impose observed temperature on surface
+		if not numpy.size(md.basalforcings.geothermalflux)==md.mesh.numberofvertices:
+			md.basalforcings.geothermalflux=numpy.zeros(md.mesh.numberofvertices)
 	else:
 		print '      no thermal boundary conditions created: no observed temperature found'
 
 	return md
+
