Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m	(revision 13634)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m	(revision 13635)
@@ -58,5 +58,5 @@
 	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=50*10^-3*ones(md.mesh.numberofvertices,1); %50 mW/m^2
+		md.basalforcings.geothermalflux=50.*10^-3*ones(md.mesh.numberofvertices,1); %50 mW/m^2
 	end
 else
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py	(revision 13635)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py	(revision 13635)
@@ -0,0 +1,67 @@
+import os
+import numpy
+from ContourToMesh import *
+
+def SetIceSheetBC(md):
+	"""
+	SETICESHEETBC - Create the boundary conditions for diagnostic and thermal models for an IceSheet with no Ice Front
+
+	   Usage:
+	      md=SetIceSheetBC(md)
+
+	   See also: SETICESHELFBC, SETMARINEICESHEETBC
+	"""
+
+	#node on Dirichlet
+	pos=numpy.nonzero(md.mesh.vertexonboundary)
+	md.diagnostic.spcvx=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.diagnostic.spcvy=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.diagnostic.spcvz=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	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 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:
+		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 neumann (Ice Front) -> none
+	if md.mesh.dimension==2:
+		md.diagnostic.icefront=numpy.zeros((0,4))
+	else:
+		md.diagnostic.icefront=numpy.zeros((0,6))
+
+	#Create zeros basal melting rate and surface mass balance if not specified
+	if numpy.all(numpy.isnan(md.surfaceforcings.precipitation)) and (md.surfaceforcings.ispdd==1):
+		md.surfaceforcings.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
+		print "      no surfaceforcings.precipitation specified: values set as zero"
+	if numpy.all(numpy.isnan(md.surfaceforcings.mass_balance)) and (md.surfaceforcings.ispdd==0):
+		md.surfaceforcings.mass_balance=numpy.zeros((md.mesh.numberofvertices,1))
+		print "      no surfaceforcings.mass_balance specified: values set as zero"
+	if numpy.all(numpy.isnan(md.basalforcings.melting_rate)):
+		md.basalforcings.melting_rate=numpy.zeros((md.mesh.numberofvertices,1))
+		print "      no basalforcings.melting_rate specified: values set as zero"
+	if numpy.all(numpy.isnan(md.balancethickness.thickening_rate)):
+		md.balancethickness.thickening_rate=numpy.zeros((md.mesh.numberofvertices,1))
+		print "      no balancethickness.thickening_rate specified: values set as zero"
+
+	md.prognostic.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+	md.balancethickness.spcthickness=float('nan')*numpy.ones((md.mesh.numberofvertices,1))
+
+	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 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
+	else:
+		print "      no thermal boundary conditions created: no observed temperature found"
+
+	return md
+