Index: /issm/trunk-jpl/test/NightlyRun/runme.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/runme.m	(revision 13970)
+++ /issm/trunk-jpl/test/NightlyRun/runme.m	(revision 13971)
@@ -67,5 +67,5 @@
 end
 % }}}
-%GET RANK and NUMPROCS for mutlithreaded runs  {{{
+%GET RANK and NUMPROCS for multithreaded runs  {{{
 rank=getfieldvalue(options,'rank',1);
 numprocs=getfieldvalue(options,'numprocs',1);
Index: /issm/trunk-jpl/test/NightlyRun/test1501.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1501.m	(revision 13970)
+++ /issm/trunk-jpl/test/NightlyRun/test1501.m	(revision 13971)
@@ -8,10 +8,10 @@
 md.transient.isthermal=0;
 
-md.timestepping.time_step=1;
+md.timestepping.time_step=1.;
 md.settings.output_frequency=1;
-md.timestepping.final_time=2000;
+md.timestepping.final_time=2000.;
 
 %Solve for thinning rate -> -1 * surface mass balance
-smb= 2*ones(md.mesh.numberofvertices,1);   
+smb= 2.*ones(md.mesh.numberofvertices,1);   
 md.surfaceforcings.mass_balance= smb;
 md.basalforcings.melting_rate= smb;
@@ -27,9 +27,9 @@
 smb = md.surfaceforcings.mass_balance;
 
-tooth= [ [ones(400,1)*(smb') - 10]' [ones(400,1)*(smb')]' ];
+tooth= [ [ones(400,1)*(smb') - 10.]' [ones(400,1)*(smb')]' ];
 smb=[ [ones(399,1)*(smb')]' smb  tooth tooth];
 
 md.surfaceforcings.mass_balance= smb;
-md.surfaceforcings.mass_balance(end+1,:)=[1:2000];
+md.surfaceforcings.mass_balance(end+1,:)=[1.:2000.];
 
 md=solve(md,TransientSolutionEnum());
Index: /issm/trunk-jpl/test/NightlyRun/test1502.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1502.m	(revision 13970)
+++ /issm/trunk-jpl/test/NightlyRun/test1502.m	(revision 13971)
@@ -9,10 +9,10 @@
 md.transient.isthermal=0;
 
-md.timestepping.time_step=1;
+md.timestepping.time_step=1.;
 md.settings.output_frequency=1;
-md.timestepping.final_time=2000;
+md.timestepping.final_time=2000.;
 
 %Solve for thinning rate -> -1 * surface mass balance
-smb= 2*ones(md.mesh.numberofvertices,1);   
+smb= 2.*ones(md.mesh.numberofvertices,1);   
 md.surfaceforcings.mass_balance= smb;
 md.basalforcings.melting_rate= smb;
@@ -28,9 +28,9 @@
 smb = md.surfaceforcings.mass_balance;
 
-tooth= [ [ones(400,1)*(smb') - 10]' [ones(400,1)*(smb')]' ];
+tooth= [ [ones(400,1)*(smb') - 10.]' [ones(400,1)*(smb')]' ];
 smb=[ [ones(399,1)*(smb')]' smb  tooth tooth];
 
 md.surfaceforcings.mass_balance= smb;
-md.surfaceforcings.mass_balance(end+1,:)=[1:2000];
+md.surfaceforcings.mass_balance(end+1,:)=[1.:2000.];
 
 md=solve(md,TransientSolutionEnum());
Index: /issm/trunk-jpl/test/NightlyRun/test1601.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1601.m	(revision 13970)
+++ /issm/trunk-jpl/test/NightlyRun/test1601.m	(revision 13971)
@@ -3,14 +3,14 @@
 md=parameterize(md,'../Par/SquareShelf.par');
 md=setflowequation(md,'macayeal','all');
-md.diagnostic.spcvx(find(md.mesh.y>0))=NaN;
-md.initialization.vx(:)=0;
-md.initialization.vy(:)=0;
-md.initialization.vel(:)=0;
+md.diagnostic.spcvx(find(md.mesh.y>0.))=NaN;
+md.initialization.vx(:)=0.;
+md.initialization.vy(:)=0.;
+md.initialization.vel(:)=0.;
 
-md.cluster=generic('name',oshostname,'np',2);
+md.cluster=generic('name',oshostname(),'np',2);
 md=solve(md,DiagnosticSolutionEnum());
 vel0=md.results.DiagnosticSolution.Vel;
 
-theta=30*pi/180;
+theta=30.*pi/180.;
 x=md.mesh.x;
 y=md.mesh.y;
@@ -27,5 +27,5 @@
 
 %Now, put CS back to normal except on the side where the spc are applied
-pos=find(x==0 | x==1000000);
+pos=find(x==0. | x==1000000.);
 md.diagnostic.referential(:)=NaN;
 md.diagnostic.referential(pos,1:3)=repmat([cos(theta),sin(theta),0],size(pos,1),1);
Index: /issm/trunk-jpl/test/NightlyRun/test1601.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1601.py	(revision 13971)
+++ /issm/trunk-jpl/test/NightlyRun/test1601.py	(revision 13971)
@@ -0,0 +1,56 @@
+import numpy
+import sys
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+md=triangle(model(),'../Exp/Square.exp',150000.)
+md=setmask(md,'all','')
+md=parameterize(md,'../Par/SquareShelf.py')
+md=setflowequation(md,'macayeal','all')
+md.diagnostic.spcvx[numpy.nonzero(md.mesh.y>0.)]=float('NaN')
+md.initialization.vx[:]=0.
+md.initialization.vy[:]=0.
+md.initialization.vel=numpy.zeros_like(md.initialization.vx)
+
+md.cluster=generic('name',oshostname(),'np',2)
+md=solve(md,DiagnosticSolutionEnum())
+vel0=md.results.DiagnosticSolution.Vel
+
+theta=30.*numpy.pi/180.
+x=md.mesh.x
+y=md.mesh.y
+md.mesh.x=numpy.cos(theta)*x-numpy.sin(theta)*y
+md.mesh.y=numpy.sin(theta)*x+numpy.cos(theta)*y
+
+md.diagnostic.referential[:,0:3]=numpy.tile([numpy.cos(theta),numpy.sin(theta),0],(md.mesh.numberofvertices,1))
+md.diagnostic.referential[:,3: ]=numpy.tile([0,0,1],(md.mesh.numberofvertices,1))
+md=solve(md,DiagnosticSolutionEnum())
+vel1=md.results.DiagnosticSolution.Vel
+
+#plotmodel(md,'data',vel0,'data',vel1,'data',vel1-vel0,'title','Cartesian CS','title','Rotated CS','title','difference')
+print "Error between Cartesian and rotated CS: %g" % (numpy.max(numpy.abs(vel0-vel1))/(numpy.max(numpy.abs(vel0))+sys.float_info.epsilon))
+
+#Now, put CS back to normal except on the side where the spc are applied
+pos=numpy.nonzero(numpy.logical_or(x==0.,x==1000000.))[0]
+md.diagnostic.referential[:]=float('NaN')
+md.diagnostic.referential[pos,0:3]=numpy.tile([numpy.cos(theta),numpy.sin(theta),0],(len(pos),1))
+md.diagnostic.referential[pos,3: ]=numpy.tile([0,0,1],(len(pos),1))
+md=solve(md,DiagnosticSolutionEnum())
+vel2=md.results.DiagnosticSolution.Vel
+
+#plotmodel(md,'data',vel0,'data',vel2,'data',vel2-vel0,'title','Cartesian CS','title','Rotated CS','title','difference')
+print "Error between Cartesian and rotated CS: %g" % (numpy.max(numpy.abs(vel0-vel2))/(numpy.max(numpy.abs(vel0))+sys.float_info.epsilon))
+
+#Fields and tolerances to track changes
+field_names     =['vel1','vel2']
+field_tolerances=[1e-11,1e-11]
+field_values=[\
+	vel1, \
+	vel2, \
+	]
Index: /issm/trunk-jpl/test/NightlyRun/test1602.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1602.m	(revision 13970)
+++ /issm/trunk-jpl/test/NightlyRun/test1602.m	(revision 13971)
@@ -4,14 +4,14 @@
 md=extrude(md,5,1.);
 md=setflowequation(md,'pattyn','all');
-md.diagnostic.spcvx(find(md.mesh.y>0))=NaN;
-md.initialization.vx(:)=0;
-md.initialization.vy(:)=0;
-md.initialization.vel(:)=0;
+md.diagnostic.spcvx(find(md.mesh.y>0.))=NaN;
+md.initialization.vx(:)=0.;
+md.initialization.vy(:)=0.;
+md.initialization.vel(:)=0.;
 
-md.cluster=generic('name',oshostname,'np',3);
+md.cluster=generic('name',oshostname(),'np',3);
 md=solve(md,DiagnosticSolutionEnum());
 vel0=md.results.DiagnosticSolution.Vel;
 
-theta=30*pi/180;
+theta=30.*pi/180.;
 x=md.mesh.x;
 y=md.mesh.y;
Index: /issm/trunk-jpl/test/NightlyRun/test1602.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1602.py	(revision 13971)
+++ /issm/trunk-jpl/test/NightlyRun/test1602.py	(revision 13971)
@@ -0,0 +1,45 @@
+import numpy
+import sys
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+md=triangle(model(),'../Exp/Square.exp',150000.)
+md=setmask(md,'../Exp/SquareShelf.exp','')
+md=parameterize(md,'../Par/SquareSheetShelf.py')
+md.extrude(5,1.)
+md=setflowequation(md,'pattyn','all')
+md.diagnostic.spcvx[numpy.nonzero(md.mesh.y>0.)]=float('NaN')
+md.initialization.vx[:]=0.
+md.initialization.vy[:]=0.
+md.initialization.vel=numpy.zeros_like(md.initialization.vx)
+
+md.cluster=generic('name',oshostname(),'np',3)
+md=solve(md,DiagnosticSolutionEnum())
+vel0=md.results.DiagnosticSolution.Vel
+
+theta=30.*numpy.pi/180.
+x=md.mesh.x
+y=md.mesh.y
+md.mesh.x=numpy.cos(theta)*x-numpy.sin(theta)*y
+md.mesh.y=numpy.sin(theta)*x+numpy.cos(theta)*y
+
+md.diagnostic.referential[:,0:3]=numpy.tile([numpy.cos(theta),numpy.sin(theta),0],(md.mesh.numberofvertices,1))
+md.diagnostic.referential[:,3: ]=numpy.tile([0,0,1],(md.mesh.numberofvertices,1))
+md=solve(md,DiagnosticSolutionEnum())
+vel1=md.results.DiagnosticSolution.Vel
+
+#plotmodel(md,'data',vel0,'data',vel1,'data',vel1-vel0,'title','Cartesian CS','title','Rotated CS','title','difference','view#all',2)
+print "Error between Cartesian and rotated CS: %g" % (numpy.max(numpy.abs(vel0-vel1))/(numpy.max(numpy.abs(vel0))+sys.float_info.epsilon))
+
+#Fields and tolerances to track changes
+field_names     =['vel1']
+field_tolerances=[1e-9]
+field_values=[\
+	vel1, \
+	];