Index: /issm/trunk-jpl/src/m/classes/maskpsl.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/maskpsl.py	(revision 20225)
+++ /issm/trunk-jpl/src/m/classes/maskpsl.py	(revision 20225)
@@ -0,0 +1,102 @@
+from MatlabFuncs import *
+from model import *
+from EnumDefinitions import *
+from numpy import *
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
+
+class maskpsl(object):
+#MASKPSL class definition
+#
+#   Usage:
+#      maskpsl=maskpsl();
+	def __init__(self,*args): # {{{
+		self.groundedice_levelset = NaN
+		self.ice_levelset         = NaN
+		self.ocean_levelset = NaN
+		self.land_levelset = NaN
+		
+		if not len(args):
+			self.setdefaultparameters()
+		else:
+			raise RuntimeError('constructor not supported')
+	# }}}
+	def __repr__(self): # {{{
+		string='   masks:'
+		
+		string="%s\n%s"%(string,fielddisplay(self,'groundedice_levelset','is ice grounded ? grounded ice if > 0, grounding line position if = 0, floating ice if < 0'))
+		string="%s\n%s"%(string,fielddisplay(self,'ice_levelset','presence of ice if < 0, icefront position if = 0, no ice if > 0'))
+		string="%s\n%s"%(string,fielddisplay(self,'ocean_levelset','is the vertex on the ocean ? yes if = 1, no if = 0'))
+		string="%s\n%s"%(string,fielddisplay(self,'land_levelset','is the vertex on the land ? yes if = 1, no if = 0'))
+	
+		return string
+	# }}}	
+	def loadobj(self): # {{{
+		# This def is directly called by matlab when a model object is
+		# loaded. Update old properties here
+
+		#2014 February 5th
+		if numel(self.ice_levelset)>1 and all(self.ice_levelset>=0):
+			print('WARNING: md.mask.ice_levelset>=0, you probably need to change the sign of this levelset')
+		return self
+	# }}}
+	def setdefaultparameters(self): # {{{
+		return self
+	# }}}
+	def checkconsistency(self,md,solution,analyses): # {{{
+
+		md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices, 1])
+		md = checkfield(md,'fieldname','mask.ice_levelset'        ,'size',[md.mesh.numberofvertices, 1])
+		md = checkfield(md,'fieldname','mask.ocean_levelset','size',[md.mesh.numberofvertices, 1])
+		md = checkfield(md,'fieldname','mask.land_levelset','size',[md.mesh.numberofvertices, 1])
+		isice=(md.mask.ice_levelset<=0)
+		if sum(isice)==0:
+			print('no ice present in the domain')
+
+		if max(md.mask.ice_levelset)<0:
+			print('no ice front provided')
+
+		icefront=sum(md.mask.ice_levelset[md.mesh.elements.astype(int)-1]==0,axis=1)
+		if (amax(icefront)==3 & strcmp(md.mesh.elementtype(),'Tria')) or (amax(icefront==6) & strcmp(md.mesh.elementtype(),'Penta')):
+			raise RuntimeError('At least one element has all nodes on ice front, change md.mask.ice_levelset to fix it')
+		
+		return md
+	# }}}
+	def extrude(self,md): # {{{
+		self.groundedice_levelset=project3d(md,'vector',self.groundedice_levelset,'type','node')
+		self.ice_levelset=project3d(md,'vector',self.ice_levelset,'type','node')
+		self.ocean_levelset=project3d(md,'vector',self.ocean_levelset,'type','node')
+		self.land_levelset=project3d(md,'vector',self.land_levelset,'type','node')
+		return self
+	# }}}
+	def mask(*args): # {{{
+		if not len(args):
+			self.setdefaultparameters()
+		else:
+			raise RuntimeError('constructor not supported')
+		return self
+	# }}}
+	
+	def marshall(self,md,fid): # {{{
+		WriteData(fid,'object',self,'class','mask','fieldname','groundedice_levelset','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'class','mask','fieldname','ice_levelset','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'class','mask','fieldname','ocean_levelset','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'class','mask','fieldname','land_levelset','format','DoubleMat','mattype',1)
+
+		# get mask of vertices of elements with ice
+		isice=md.mask.ice_levelset<=0.
+		vlist = zeros(md.mesh.numberofvertices,1)
+		pos=find(sum(isice(md.mesh.elements),2)>0)
+		vlist[md.mesh.elements[pos,:]]=1
+		WriteData(fid,'data',vlist,'enum',IceMaskNodeActivationEnum(),'format','DoubleMat','mattype',1)
+
+	# }}}
+	def savemodeljs(self,fid,modelname): # {{{
+
+		writejs1Darray(fid,[modelname, '.mask.groundedice_levelset'],self.groundedice_levelset)
+		writejs1Darray(fid,[modelname, '.mask.ice_levelset'],self.ice_levelset)
+		writejs1Darray(fid,[modelname, '.mask.ocean_levelset'],self.ocean_levelset)
+		writejs1Darray(fid,[modelname, '.mask.land_levelset'],self.land_levelset)
+
+	# }}}
Index: /issm/trunk-jpl/src/m/classes/mesh3dsurface.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh3dsurface.py	(revision 20224)
+++ /issm/trunk-jpl/src/m/classes/mesh3dsurface.py	(revision 20225)
@@ -3,4 +3,7 @@
 from EnumDefinitions import *
 from numpy import *
+from fielddisplay import fielddisplay
+from checkfield import checkfield
+from WriteData import WriteData
 
 class mesh3dsurface(object):
@@ -33,5 +36,5 @@
 		self.extractedelements           = NaN
 		
-		if len(args)==0:
+		if not len(args):
 			self.setdefaultparameters()
 		elif len(args)==1:
@@ -46,9 +49,8 @@
 			raise RuntimeError('constructor not supported')	
 
-		#set defaults
-		self.setdefaultparameters()
-	#}}}
+	# }}}
 	def __repr__(self): # {{{
 		string='   2D tria Mesh (horizontal):'
+		
 		string+='\n      Elements and vertices:'
 		string="%s\n%s"%(string,fielddisplay(self,'numberofelements','number of elements'))
@@ -76,5 +78,7 @@
 		string="%s\n%s"%(string,fielddisplay(self,'extractedvertices','vertices extracted from the model()'))
 		string="%s\n%s"%(string,fielddisplay(self,'extractedelements','elements extracted from the model()')) 
-		# }}}
+		
+		return string
+	# }}}
 	def loadobj(self): # {{{
 		# This def is directly called by matlab when a model() selfect is
@@ -92,5 +96,6 @@
 				else:
 					self.epsg=3031
-		# }}}
+		return self
+	# }}}
 	def setdefaultparameters(self): # {{{
 
@@ -102,14 +107,14 @@
 		self.average_vertex_connectivity=25
 		return self
-		# }}}
+	# }}}
 	def checkconsistency(self,md,solution,analyses): # {{{
 
-		md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
-		md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
-		md = checkfield(md,'fieldname','mesh.z','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
-		md = checkfield(md,'fieldname','mesh.lat','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
-		md = checkfield(md,'fieldname','mesh.long','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
-		md = checkfield(md,'fieldname','mesh.r','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices,1])
-		md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',[md.mesh.numberofvertices])
+		md = checkfield(md,'fieldname','mesh.x','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'fieldname','mesh.y','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'fieldname','mesh.z','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'fieldname','mesh.lat','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'fieldname','mesh.long','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'fieldname','mesh.r','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'fieldname','mesh.elements','NaN',1,'Inf',1,'>',0,'values',arange(1,md.mesh.numberofvertices+1))
 		md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,3])
 		if any(not ismember([md.mesh.numberofvertices],sort(unique(md.mesh.elements)))):
@@ -124,16 +129,16 @@
 			
 		return md
-		# }}}
+	# }}}
 	def marshall(self,md,fid): # {{{
-		WriteData(fid,'enum',DomainTypeEnum(),'data',StringToEnum(['Domain', domaintype(obj)]),'format','Integer')
-		WriteData(fid,'enum',DomainDimensionEnum(),'data',dimension(obj),'format','Integer')
-		WriteData(fid,'enum',MeshElementtypeEnum(),'data',StringToEnum(elementtype(obj)),'format','Integer')
+		WriteData(fid,'enum',DomainTypeEnum(),'data',StringToEnum(['Domain', self.domaintype()]),'format','Integer')
+		WriteData(fid,'enum',DomainDimensionEnum(),'data',self.dimension(),'format','Integer')
+		WriteData(fid,'enum',MeshElementtypeEnum(),'data',StringToEnum(self.elementtype()),'format','Integer')
 		WriteData(fid,'object',obj,'class','mesh','fieldname','x','format','DoubleMat','mattype',1)
 		WriteData(fid,'object',obj,'class','mesh','fieldname','y','format','DoubleMat','mattype',1)
 		WriteData(fid,'object',obj,'class','mesh','fieldname','z','format','DoubleMat','mattype',1)
-		WriteData(fid,'enum',MeshLatEnum(),'data',obj.lat,'format','DoubleMat','mattype',1)
-		WriteData(fid,'enum',MeshLongEnum(),'data',obj.long,'format','DoubleMat','mattype',1)
+		WriteData(fid,'enum',MeshLatEnum(),'data',md.mesh.lat,'format','DoubleMat','mattype',1)
+		WriteData(fid,'enum',MeshLongEnum(),'data',md.mesh.long,'format','DoubleMat','mattype',1)
 		WriteData(fid,'object',obj,'class','mesh','fieldname','r','format','DoubleMat','mattype',1)
-		WriteData(fid,'enum',MeshZEnum(),'data',zeros(obj.numberofvertices,1),'format','DoubleMat','mattype',1)
+		WriteData(fid,'enum',MeshZEnum(),'data',zeros(md.mesh.numberofvertices,1),'format','DoubleMat','mattype',1)
 		WriteData(fid,'object',obj,'class','mesh','fieldname','elements','format','DoubleMat','mattype',2)
 		WriteData(fid,'object',obj,'class','mesh','fieldname','numberofelements','format','Integer')
@@ -141,14 +146,14 @@
 		WriteData(fid,'object',obj,'class','mesh','fieldname','average_vertex_connectivity','format','Integer')
 		WriteData(fid,'object',obj,'class','mesh','fieldname','vertexonboundary','format','DoubleMat','mattype',1)
-		# }}}
+	# }}}
 	def domaintype(self): # {{{
 		return '3Dsurface'
-		# }}}
+	# }}}
 	def dimension(self): # {{{
 		return 2
-		# }}}
+	# }}}
 	def elementtype(self): # {{{
 		return 'Tria'
-		# }}}
+	# }}}
 	def processmesh(self,options): # {{{
 	
@@ -161,5 +166,5 @@
 		z        = self.z
 		return [x, y, z, elements, is2d, isplanet]
-		# }}}
+	# }}}
 	def savemodeljs(self,fid,modelname): # {{{
 	
@@ -185,4 +190,4 @@
 		writejs1Darray(fid,[modelname, '.mesh.extractedelements'],self.extractedelements)
 
-		# }}}
+	# }}}
 	
Index: /issm/trunk-jpl/src/m/coordsystems/gmtmask.py
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/gmtmask.py	(revision 20225)
+++ /issm/trunk-jpl/src/m/coordsystems/gmtmask.py	(revision 20225)
@@ -0,0 +1,81 @@
+from MatlabFuncs import *
+from model import *
+from EnumDefinitions import *
+from numpy import *
+from os import getenv, putenv
+import subprocess
+
+def gmtmask(lat,long,*varargin):
+#GMTMASK - figure out which lat,long points are on the ocean
+#
+#   Usage:
+#      mask.ocean = gmtmask(md.mesh.lat,md.mesh.long);
+#
+	lenlat=len(lat)
+	mask=empty(lenlat)
+	
+	#are we doing a recursive call? 
+	if len(varargin)==3:
+		recursive=1
+	else:
+		recursive=0
+
+	if recursive:
+		string='             recursing: num vertices #i'+str(lenlat)
+	else:
+		string='gmtmask: num vertices #i'+str(lenlat)
+	
+	#Check lat and long size is not more than 50,000 If so, recursively call gmtmask: 
+
+	if lenlat>50000:
+		for i in range(ceil(lenlat/50000)):
+			j=(i+1)*50000-1
+			if j>lenlat:
+				j=lenlat
+			mask[i:j]=gmtmask(lat[i:j],long[i:j],1)
+		return mask
+	
+	
+	#First, write our lat,long file for gmt:
+	nv=lenlat
+	savetxt('./all_vertices.txt',[long, lat, arange(1,nv+1)],delimiter='\t')
+
+	#Avoid bypassing of the ld library path by Matlab (:()
+	try:
+		issmdir
+	except:
+		issmdir=getenv('ISSM_DIR')
+	try:
+		ismac
+	except:
+		ismac=False	
+
+	if ismac:
+		dyld_library_path_old=getenv('DYLD_LIBRARY_PATH')
+		putenv('DYLD_LIBRARY_PATH',issmdir+'/externalpackages/curl/install/lib:'+issmdir+'/externalpackages/hdf5/install/lib:'+issmdir+'/externalpackages/netcdf/install/lib')
+		
+	#figure out which vertices are on the ocean, which one on the continent:
+	subprocess.call(issmdir+'/externalpackages/gmt/install/bin/gmt gmtselect ./all_vertices.txt -h0 -Df -R0/360/-90/90  -A0 -JQ180/200 -Nk/s/s/k/s > ./oce_vertices.txt',shell=True)
+
+	#reset DYLD_LIBRARY_PATH to what it was: 
+	if ismac:
+		putenv('DYLD_LIBRARY_PATH',dyld_library_path_old)
+	
+	#read the con_vertices.txt file and flag our mesh vertices on the continent
+	fid=open('./oce_vertices.txt','r')
+	line=fid.readline()
+	line=fid.readline()
+	oce_vertices=[]
+	while line:
+		ind=str2num(line.split()[3])
+		oce_vertices.append(ind)
+		line=fid.readline()
+	fid.close()
+
+	mask=zeros([nv,1])
+	mask[oce_vertices]=1
+	
+	subprocess.call('rm -rf ./all_vertices.txt ./oce_vertices.txt ./gmt.history',shell=True)
+	if not recursive:
+		string='gmtmask: done'
+	return mask