Index: /issm/trunk-jpl/jenkins/ross-debian_linux-solid_earth
===================================================================
--- /issm/trunk-jpl/jenkins/ross-debian_linux-solid_earth	(revision 25514)
+++ /issm/trunk-jpl/jenkins/ross-debian_linux-solid_earth	(revision 25515)
@@ -80,3 +80,3 @@
 #
 MATLAB_NROPTIONS="'benchmark','slr','exclude',[2021]"
-PYTHON_NROPTIONS="--benchmark slr --exclude 2004 2021"
+PYTHON_NROPTIONS="--benchmark slr --exclude 2021"
Index: /issm/trunk-jpl/src/m/geometry/delaunay.py
===================================================================
--- /issm/trunk-jpl/src/m/geometry/delaunay.py	(revision 25514)
+++ /issm/trunk-jpl/src/m/geometry/delaunay.py	(revision 25515)
@@ -9,11 +9,20 @@
 
     Usage:
-        # Creates a 2-D or 3-D Delaunay triangulation from the points in a matrix P. The output DT is a three-column (for two dimensions) or four-column (for three dimensions) matrix where each row contains the row indices of the input points that make up a triangle or tetrahedron in the triangulation.
+        # Creates a 2-D or 3-D Delaunay triangulation from the points in a 
+        matrix P. The output DT is a three-column (for two dimensions) or 
+        four-column (for three dimensions) matrix where each row contains the 
+        row indices of the input points that make up a triangle or tetrahedron 
+        in the triangulation.
+
         DT = delaunay(P)
 
-        # Creates a 2-D Delaunay triangulation from the points in vectors x and y.
+        # Creates a 2-D Delaunay triangulation from the points in vectors x and 
+        y.
+
         DT = delaunay(x, y)
 
-        # Creates a 3-D Delaunay triangulation from the points in vectors x, y, and z.
+        # Creates a 3-D Delaunay triangulation from the points in vectors x, y, 
+        and z.
+
         DT = delaunay(x, y, z)
 
Index: /issm/trunk-jpl/src/m/geometry/delaunayn.py
===================================================================
--- /issm/trunk-jpl/src/m/geometry/delaunayn.py	(revision 25515)
+++ /issm/trunk-jpl/src/m/geometry/delaunayn.py	(revision 25515)
@@ -0,0 +1,52 @@
+import numpy as np
+from scipy.spatial import Delaunay # See also: https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html?highlight=delaunay#scipy.spatial.Delaunay
+
+
+def delaunayn(X, options=[]):
+    """DELAUNAYN - N-D Delaunay triangulation
+
+    Replicates behaviour of MATLAB's 'delaunayn'.
+
+    Usage:
+        # Computes a set of simplices such that no data points of X are 
+        contained in any circumspheres of the simplices. The set of simplices 
+        forms the Delaunay triangulation. X is an m-by-n array representing m 
+        points in n-dimensional space. T is a numt-by-(n+1) array where each 
+        row contains the indices into X of the vertices of the corresponding 
+        simplex.
+
+        T = delaunayn(X)
+
+        # Specifies a cell array of options. The default options are:
+        - {'Qt','Qbb','Qc'} for 2- and 3-dimensional input
+        - {'Qt','Qbb','Qc','Qx'} for 4 and higher-dimensional input
+        If options is [], the default options used. If options is {''}, no options are used, not even the default.
+
+        T = delaunayn(X, options)
+
+    Sources:
+    - https://www.mathworks.com/help/matlab/ref/delaunayn.html
+    - https://stackoverflow.com/questions/36604172/difference-between-matlab-delaunayn-and-scipy-delaunay
+    """
+
+    # Get dimensions of points
+    ndims = X.shape[1]
+
+    # Set up default options
+    if options == []:
+        options = 'Qt Qbb Qc' if ndims <= 3 else 'Qt Qbb Qc Qx'
+
+    # Triangulate
+    triangles = Delaunay(X, qhull_options=options).simplices
+
+    # Strip the zero area/volume simplices that may have been created in the 
+    # presence of degeneracy.
+    keep = np.ones((len(triangles), ), dtype=bool)
+
+    for i, t in enumerate(triangles):
+        if abs(np.linalg.det(np.hstack((X[t], np.ones([1, ndims + 1]).T)))) < 1e-15:
+            keep[i] = False
+
+    triangles = triangles[keep].astype(int) # NOTE: Need to covert to array of int so that it can be fetched properly by core
+
+    return triangles
Index: /issm/trunk-jpl/src/m/solve/WriteData.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/WriteData.py	(revision 25514)
+++ /issm/trunk-jpl/src/m/solve/WriteData.py	(revision 25515)
@@ -8,9 +8,8 @@
 
 def WriteData(fid, prefix, *args):
-    """
-    WRITEDATA - write model field in binary file
-
-       Usage:
-          WriteData(fid, varargin)
+    """WRITEDATA - write model field in binary file
+
+    Usage:
+        WriteData(fid, varargin)
     """
 
Index: /issm/trunk-jpl/test/NightlyRun/runme.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/runme.m	(revision 25514)
+++ /issm/trunk-jpl/test/NightlyRun/runme.m	(revision 25515)
@@ -2,5 +2,5 @@
 %RUNME - test deck for ISSM nightly runs
 %
-%   In a test deck directory (tests/Vertification/NightlyRun for example)
+%   In a test deck directory (test/Vertification/NightlyRun for example)
 %   The following command will launch all the existing tests:
 %   >> runme
Index: /issm/trunk-jpl/test/NightlyRun/runme.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/runme.py	(revision 25514)
+++ /issm/trunk-jpl/test/NightlyRun/runme.py	(revision 25515)
@@ -20,22 +20,23 @@
 
 def runme(id=None, exclude=None, benchmark='nightly', procedure='check', output='none', rank=1, numprocs=1):
-    """
-    RUNME - test deck for ISSM nightly runs
-
-        In a test deck directory (tests / Vertification / NightlyRun for example)
-        The following command will launch all the existing tests:
-    >>     runme()
-        To run the tests 101 and 102:
-    >>     runme(id = [101, 102])
-        etc...
-
-        Available options:
-            'id'            followed by the list of ids or (parts of) test names requested
-                            Note: runs all tests by default
-
-            'exclude'       ids or (parts of) test names to be excluded (same format as id)
-                            Note: exclude does nothing if 'id' is specified with different values
-
-            'benchmark'     'all' (all of the tests)
+    """RUNME - test deck for ISSM nightly runs
+
+    In a test deck directory (tests/Vertification / NightlyRun for example) the 
+    following command will launch all the existing tests:
+
+        >> ./runme.py
+
+    To run the tests 101 and 102:
+
+        >> ./runme.py -i [101, 102]
+
+    Available options:
+        -i/--id             followed by the list of ids or (parts of) test names requested
+                            NOTE: runs all tests by default
+
+        -e/--exclude        ids or (parts of) test names to be excluded (same format as id)
+                            NOTE: exclude does nothing if 'id' is specified with different values
+
+        -b/--benchmark      'all' (all of the tests)
                             'nightly' (nightly run / daily run)
                             'ismip'  : validation of ismip - hom tests
@@ -46,23 +47,24 @@
                             'slr'    : validation of slr tests
 
-            'procedure'     'check' : run the test (default)
+        -p/--procedure      'check' : run the test (default)
                             'update': update the archive
 
-        Usage:
-            runme(varargin)
-
-        Examples:
-            runme()
-            runme(101)
-            runme('SquareShelf')
-            runme(exclude = 2001)
-            runme(exclude = 'Dakota', benchmark = 'all')
-            runme(id = [[101, 102], ['Dakota', 'Slr']])
-
-        TODO:
-        - At '#disp test result', make sure precision of output matches that of 
+    Usage:
+        ./runme.py [option [args]]
+
+    Examples:
+        ./runme.py
+        ./runme.py -i 101
+        ./runme.py -i 'SquareShelf'
+        ./runme.py -e 2001
+        ./runme.py -e 'Dakota' --benchmark 'all'
+        ./runme.py -i [[101, 102], ['Dakota', 'Slr']]
+
+    TODO:
+    - At '#disp test result', make sure precision of output matches that of 
         MATLAB.
-        - Check for failures that do not raise exceptions (for example, 'Standard exception'; see also jenkins/jenkins.sh). These should be counted as failures.
-        """
+    - Check for failures that do not raise exceptions (for example, 'Standard exception'; see also jenkins/jenkins.sh). These should be counted as failures.
+    """
+
     #Get ISSM_DIR variable
     ISSM_DIR = os.environ['ISSM_DIR']
@@ -75,20 +77,20 @@
         print(("runme warning: benchmark '{}' not supported, defaulting to test 'nightly'.".format(benchmark)))
         benchmark = 'nightly'
-        # }}}
-        #GET procedure {{{
+    # }}}
+    #GET procedure {{{
     if procedure not in ['check', 'update']:
         print(("runme warning: procedure '{}' not supported, defaulting to test 'check'.".format(procedure)))
         procedure = 'check'
-        # }}}
-        #GET output {{{
+    # }}}
+    #GET output {{{
     if output not in ['nightly', 'none']:
         print(("runme warning: output '{}' not supported, defaulting to test 'none'.".format(output)))
         output = 'none'
-        # }}}
-        #GET RANK and NUMPROCS for multithreaded runs {{{
+    # }}}
+    #GET RANK and NUMPROCS for multithreaded runs {{{
     if (numprocs < rank):
         numprocs = 1
-        # }}}
-        #GET ids  {{{
+    # }}}
+    #GET ids  {{{
     flist = glob('test*.py')  #File name must start with 'test' and must end by '.py' and must be different than 'test.py'
     list_ids = [int(file[4:-3]) for file in flist if not file == 'test.py']  #Keep test id only (skip 'test' and '.py')
Index: /issm/trunk-jpl/test/NightlyRun/test2004.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2004.m	(revision 25514)
+++ /issm/trunk-jpl/test/NightlyRun/test2004.m	(revision 25515)
@@ -154,5 +154,6 @@
 			latAIS=delH{:,2};
 			delHAIS=delH{:,3};
-			index=delaunay(longAIS,latAIS);
+			points=[longAIS,latAIS];
+			index=delaunayn(points);
 			lat=md.mesh.lat; 
 			long=md.mesh.long+360;
@@ -190,5 +191,5 @@
 
 	%recover lat,long: 
-	[md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); 
+	[md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326');
 
 	%mask:  %{{{
@@ -262,5 +263,6 @@
 		latGIS=delH(:,2);
 		delHGIS=delH(:,3);
-		index=delaunay(longGIS,latGIS);
+		points=[longGIS,latGIS];
+		index=delaunayn(points);
 		lat=md.mesh.lat;
 		long=md.mesh.long+360;
@@ -274,5 +276,6 @@
 		latGLA=delH(:,2);
 		delHGLA=sum(delH(:,3:end),2);
-		index=delaunay(longGLA,latGLA);
+		points=[longGLA,latGLA];
+		index=delaunayn(points);
 		lat=md.mesh.lat; 
 		long=md.mesh.long+360; 
@@ -409,5 +412,4 @@
 md=solve(md,'Sealevelrise');
 Srotation=md.results.SealevelriseSolution.Sealevel;
-
 %}}}