Index: /issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py
===================================================================
--- /issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py	(revision 12190)
+++ /issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py	(revision 12191)
@@ -1,3 +1,4 @@
-def SetIceShelfBC(md,icefrontfile='')
+from numpy import *
+def SetIceShelfBC(md,icefrontfile=''):
 	#SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
 	#
@@ -22,40 +23,41 @@
 			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 & nodeinsideicefront);
+		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=find(md.mesh.vertexonboundary & ~nodeonicefront);
-	md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
-	md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
-	md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
-	md.diagnostic.spcvz(pos)=0;
-	md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);
+	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);
 
 	#Dirichlet Values
-	if (length(md.inversion.vx_obs)==md.mesh.numberofvertices & length(md.inversion.vy_obs)==md.mesh.numberofvertices)
-		disp('      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
-		disp('      boundary conditions for diagnostic model: spc set as zero');
-	end
+	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]
+	else:
+		print '      boundary conditions for diagnostic model: spc set as zero'
 
 	#segment on Ice Front
 	#segment on Neumann (Ice Front)
-	pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
-	if (md.mesh.dimension==2)
-		pressureload=md.mesh.segments(pos,:);
-	elseif md.mesh.dimension==3
-		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)];
+	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,:];
+	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=1:md.mesh.numberoflayers-1,
-			pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
-		end
-	end
+		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 ];
 
 	#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])];
 
 	#plug onto model
@@ -63,31 +65,27 @@
 
 	#Create zeros basalforcings and surfaceforcings
-	if isnan(md.surfaceforcings.precipitation),
-		md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
-		disp('      no surfaceforcings.precipitation specified: values set as zero');
-	end
-	if isnan(md.surfaceforcings.mass_balance),
-		md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1);
-		disp('      no surfaceforcings.mass_balance specified: values set as zero');
-	end
-	if isnan(md.basalforcings.melting_rate),
-		md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
-		disp('      no basalforcings.melting_rate specified: values set as zero');
-	end
-	if isnan(md.balancethickness.thickening_rate),
-		md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1);
-		disp('      no balancethickness.thickening_rate specified: values set as zero');
-	end
+	if isnan(md.surfaceforcings.precipitation):
+		md.surfaceforcings.precipitation=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);
+		print '      no surfaceforcings.mass_balance specified: values set as zero'
+	if isnan(md.basalforcings.melting_rate):
+		md.basalforcings.melting_rate=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);
+		print '      no balancethickness.thickening_rate specified: values set as zero'
 
-	md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
-	md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
+	md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices);
+	md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices);
 
-	if (length(md.initialization.temperature)==md.mesh.numberofvertices),
-		md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
-		pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); #impose observed temperature on surface
-		if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
-			md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1);
-		end
-	else
-		disp('      no thermal boundary conditions created: no observed temperature found');
-	end
+	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);
+	else:
+		print '      no thermal boundary conditions created: no observed temperature found'
+
+	return md