Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m	(revision 16317)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m	(revision 16318)
@@ -53,5 +53,8 @@
 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 isprop(md.mesh,'vertexonsurface')
+		pos=find(md.mesh.vertexonsurface);
+		md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
+	end
 	if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
 		md.basalforcings.geothermalflux=50.*10^-3*ones(md.mesh.numberofvertices,1); %50 mW/m^2
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py	(revision 16317)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py	(revision 16318)
@@ -54,7 +54,7 @@
 	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))
-#		pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
-		pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
-		md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
+		if hasattr(md.mesh,'vertexonsurface'):
+			pos=numpy.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
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m	(revision 16317)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m	(revision 16318)
@@ -38,8 +38,10 @@
 
 %First find segments that are not completely on the front
-if(md.mesh.dimension==2),
+if strcmp(meshtype(md.mesh),'3D'),
+	numbernodesfront=4;
+elseif strcmp(meshtype(md.mesh),'2Dhorizontal'),
 	numbernodesfront=2;
-else 
-	numbernodesfront=4;
+else
+	error('mesh type not supported yet');
 end
 segmentsfront=md.mask.ice_levelset(md.mesh.segments(:,1:numbernodesfront))==0;
@@ -84,5 +86,8 @@
 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 isprop(md.mesh,'vertexonsurface'),
+		pos=find(md.mesh.vertexonsurface);
+		md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
+	end
 	if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
 		md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1);
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py	(revision 16317)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py	(revision 16318)
@@ -2,4 +2,5 @@
 import numpy
 from ContourToMesh import *
+from MatlabFuncs import *
 
 def SetIceShelfBC(md,icefrontfile=''):
@@ -43,8 +44,10 @@
 
 	#First find segments that are not completely on the front
-	if md.mesh.dimension==2:
-		numbernodesfront=2
+	if strcmp(md.mesh.meshtype(),'3D'):
+		numbernodesfront=4;
+	elif strcmp(md.mesh.meshtype(),'2Dhorizontal'):
+		numbernodesfront=2;
 	else:
-		numbernodesfront=4
+		raise	error('mesh type not supported yet')
 	if any(md.mask.ice_levelset<=0):
 		values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1]
@@ -93,7 +96,7 @@
 	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))
-#		pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
-		pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
-		md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
+		if hasattr(md.mesh,'vertexonsurface'):
+			pos=numpy.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))
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m	(revision 16317)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m	(revision 16318)
@@ -102,5 +102,8 @@
 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 isprop(md.mesh,'vertexonsurface'),
+		pos=find(md.mesh.vertexonsurface);
+		md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
+	end
 	if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
 		md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1);
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 16317)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 16318)
@@ -101,7 +101,7 @@
 	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))
-#		pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
-		pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
-		md.thermal.spctemperature[pos]=md.initialization.temperature[pos]    #impose observed temperature on surface
+		if hasattr(md.mesh,'vertexonsurface'):
+			pos=numpy.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))