Index: /issm/trunk-jpl/src/jl/md/exp.jl
===================================================================
--- /issm/trunk-jpl/src/jl/md/exp.jl	(revision 26708)
+++ /issm/trunk-jpl/src/jl/md/exp.jl	(revision 26709)
@@ -136,2 +136,99 @@
 	return contours
 end# }}}
+#ContourToNodes{{{
+"""
+	ContourToNodes - Flag points that are in contour
+
+   More doc to come later....
+
+	Usage:
+		exp=expread(filename)
+
+# Examples:
+```julia-repl
+julia> exp=expread('domainoutline.exp')
+```
+
+# Arguments:
+- filename: the ARGUS file to read
+"""
+function ContourToNodes(x::Vector{Float64},y::Vector{Float64},filename::String,edgevalue::Float64)
+
+	#Read input file
+	contours = expread(filename)
+
+	#Initialize output
+	nbpts = length(x)
+	flags = zeros(Bool,nbpts)
+
+	#Loop over contours
+	for c in 1:length(contours)
+
+		#Get current contours
+		contour = contours[c]
+		xp      = contour.x
+		yp      = contour.y
+
+		#Check that we are within box
+		xmin = minimum(xp); xmax = maximum(xp)
+		ymin = minimum(yp); ymax = maximum(yp)
+
+		#Loop over all points provided
+		for ii in 1:nbpts
+
+			#If this node is already within one of the contours, do not change it
+			if(flags[ii]) continue end
+
+			#Are we within bounds?
+			if(x[ii]<xmin || x[ii]>xmax || y[ii]<ymin || y[ii]>ymax) continue end
+
+			#we are potentially inside... perform pnpoly test
+			flags[ii] = pnpoly(xp, yp, x[ii], y[ii], edgevalue)
+		end
+	end
+	
+	return flags
+end# }}}
+
+function pnpoly(xp::Vector{Float64},yp::Vector{Float64},x::Float64,y::Float64,edgevalue::Float64) #{{{
+
+	npol = length(xp)
+
+	#Do we need to test for colinearity?
+	if(edgevalue!=2)
+		i = 1
+		j = npol
+		while(i<=npol)
+
+			n1 = (yp[i]-yp[j])^2 + (xp[i]-xp[j])^2
+			n2 = (y-yp[j])^2 + (x-xp[j])^2
+
+			normp=sqrt(n1*n2)
+			scalar=(yp[i]-yp[j])*(y-yp[j])+(xp[i]-xp[j])*(x-xp[j])
+
+			if (scalar == normp)
+				if (n2<=n1)
+					return edgevalue
+				end
+			end
+
+			j =  i
+			i += 1
+		end
+	end
+
+	#second test : point is neither on a vertex, nor on a side, where is it ?
+	i = 1
+	j = npol
+	c = false
+	while(i<=npol)
+		if (((yp[i]<=y && y<yp[j]) || (yp[j]<=y && y<yp[i])) &&
+            (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
+         c = !c
+		end
+
+		j =  i
+		i += 1
+	end
+	return c
+end# }}}