Index: /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt	(revision 14134)
@@ -24,2 +24,5 @@
 propagate into the sixth digit of the thickness, which is enough to fail NR.
 
+test1109,test1110
+Ugly crashes in solver, but same behavior as Matlab.
+
Index: /issm/trunk-jpl/test/NightlyRun/test1101.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1101.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1101.m	(revision 14134)
@@ -1,7 +1,7 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %Pattyn and Payne 2006
 printingflag=false;
 
-L_list={5000,10000,20000,40000,80000,160000};
+L_list={5000.,10000.,20000.,40000.,80000.,160000.};
 results={};
 minvx=[];
@@ -26,13 +26,13 @@
 
 	pos=find(md.mesh.vertexonbed);
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
+	md.diagnostic.spcvx(pos)=0.;
+	md.diagnostic.spcvy(pos)=0.;
 
 	%Create MPCs to have periodic boundary conditions
-	posx=find(md.mesh.x==0);
+	posx=find(md.mesh.x==0.);
 	posx2=find(md.mesh.x==max(md.mesh.x));
 
-	posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
-	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
+	posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
+	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
 
 	md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
@@ -55,5 +55,5 @@
 		set(gcf,'Color','w')
 		printmodel(['ismipapattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipapattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+		system(['mv ismipapattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 	end
 	plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
@@ -61,5 +61,5 @@
 		set(gcf,'Color','w')
 		printmodel(['ismipapattynvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipapattynvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+		system(['mv ismipapattynvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 	end
 	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
@@ -67,23 +67,23 @@
 		set(gcf,'Color','w')
 		printmodel(['ismipapattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipapattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+		system(['mv ismipapattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 	end
 
-	if(L==5000),
+	if(L==5000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','')
-	elseif(L==10000),
+	elseif(L==10000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','')
-	elseif(L==20000),
+	elseif(L==20000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','')
-	elseif(L==40000),
+	elseif(L==40000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','')
-	elseif(L==80000),
+	elseif(L==80000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','')
-	elseif(L==160000),
+	elseif(L==160000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
@@ -92,5 +92,5 @@
 		set(gcf,'Color','w')
 		printmodel(['ismipapattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipapattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+		system(['mv ismipapattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 	end
 end
@@ -101,5 +101,5 @@
 	set(gcf,'Color','w')
 	printmodel('ismipapattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipapattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+	system(['mv ismipapattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 end
 plot([5 10 20 40 80 160],maxvx);ylim([0 120]);xlim([0 160])
@@ -107,8 +107,9 @@
 	set(gcf,'Color','w')
 	printmodel('ismipapattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipapattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+	system(['mv ismipapattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 end
+
 %Fields and tolerances to track changes
-field_names     ={ ...
+field_names     ={...
 	'Vx5km','Vy5km','Vz5km',...
 	'Vx10km','Vy10km','Vz10km',...
Index: /issm/trunk-jpl/test/NightlyRun/test1101.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1101.py	(revision 14134)
+++ /issm/trunk-jpl/test/NightlyRun/test1101.py	(revision 14134)
@@ -0,0 +1,152 @@
+import numpy
+import shutil
+from model import *
+from squaremesh import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This test is a test from the ISMP-HOM Intercomparison project.
+Pattyn and Payne 2006
+"""
+
+printingflag=False
+
+L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
+results=[]
+minvx=[]
+maxvx=[]
+
+for L in L_list:
+	nx=20    #numberof nodes in x direction
+	ny=20
+	md=model()
+	md=squaremesh(md,L,L,nx,ny)
+	md=setmask(md,'','')    #ice sheet test
+	md=parameterize(md,'../Par/ISMIPA.py')
+	md.extrude(9,1.)
+
+	md=setflowequation(md,'pattyn','all')
+
+	#Create dirichlet on the bed only
+	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))
+
+	pos=numpy.nonzero(md.mesh.vertexonbed)
+	md.diagnostic.spcvx[pos]=0.
+	md.diagnostic.spcvy[pos]=0.
+
+	#Create MPCs to have periodic boundary conditions
+	posx=numpy.nonzero(md.mesh.x==0.)[0]
+	posx2=numpy.nonzero(md.mesh.x==numpy.max(md.mesh.x))[0]
+
+	posy=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==0.,md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]    #Don't take the same nodes two times
+	posy2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==numpy.max(md.mesh.y),md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]
+
+	md.diagnostic.vertex_pairing=numpy.vstack((numpy.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),numpy.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))
+
+	#Compute the diagnostic
+	md.cluster=generic('name',oshostname(),'np',8)
+	md=solve(md,DiagnosticSolutionEnum())
+
+	#Plot the results and save them
+	vx=md.results.DiagnosticSolution.Vx
+	vy=md.results.DiagnosticSolution.Vy
+	vz=md.results.DiagnosticSolution.Vz
+	results.append(md.results.DiagnosticSolution)
+	minvx.append(numpy.min(vx[-md.mesh.numberofvertices2d:]))
+	maxvx.append(numpy.max(vx[-md.mesh.numberofvertices2d:]))
+
+	#Now plot vx, vy, vz and vx on a cross section
+#	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipapattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipapattynvx%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+#	plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipapattynvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipapattynvy%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+#	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipapattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipapattynvz%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+
+	if   (L==5000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','')
+	elif (L==10000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','')
+	elif (L==20000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','')
+	elif (L==40000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','')
+	elif (L==80000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','')
+	elif (L==160000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipapattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipapattynvxsec%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+
+#Now plot the min and max values of vx for each size of the square
+#plot([5 10 20 40 80 160],minvx);ylim([0 18]);xlim([0 160])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipapattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipapattynminvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+#plot([5 10 20 40 80 160],maxvx);ylim([0 120]);xlim([0 160])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipapattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipapattynmaxvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+
+#Fields and tolerances to track changes
+field_names     =[\
+	'Vx5km','Vy5km','Vz5km',\
+	'Vx10km','Vy10km','Vz10km',\
+	'Vx20km','Vy20km','Vz20km',\
+	'Vx40km','Vy40km','Vz40km',\
+	'Vx80km','Vy80km','Vz80km',\
+	'Vx160km','Vy160km','Vz160km'
+]
+field_tolerances=[\
+	1e-09,1e-09,1e-09,\
+	1e-10,1e-10,1e-09,\
+	1e-09,1e-09,1e-09,\
+	1e-09,1e-08,1e-09,\
+	1e-08,1e-08,1e-08,\
+	1e-08,1e-07,1e-08,\
+]
+field_values=[]
+for result in results:
+	field_values=field_values+[\
+		result.Vx,\
+		result.Vy,\
+		result.Vz,\
+		]
Index: /issm/trunk-jpl/test/NightlyRun/test1102.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1102.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1102.m	(revision 14134)
@@ -1,7 +1,7 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %Pattyn and Payne 2006
 printingflag=false;
 
-L_list={5000,10000,20000,40000,80000,160000};
+L_list={5000.,10000.,20000.,40000.,80000.,160000.};
 results={};
 minvx=[];
@@ -17,5 +17,5 @@
 
 %	%Find elements at the corner and extract model
-%	posnodes=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
+%	posnodes=find((md.mesh.x==0. | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0. | md.mesh.y==max(md.mesh.y)));
 %	[a,b]=find(ismember(md.mesh.elements,posnodes));
 %	elements=ones(md.mesh.numberofelements,1);
@@ -29,23 +29,23 @@
 	%Create dirichlet on the bed only
 	pos=find(md.mesh.vertexonbed);
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
-	md.diagnostic.spcvz(pos)=0;
+	md.diagnostic.spcvx(pos)=0.;
+	md.diagnostic.spcvy(pos)=0.;
+	md.diagnostic.spcvz(pos)=0.;
 
-	%Create MPCs to have periodic boundary conditions
-	%posx=find(md.mesh.x==0);
-	%posx2=find(md.mesh.x==max(md.mesh.x));
-	%posx=find(md.mesh.x==0 & md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
-	%posx2=find(md.mesh.x==max(md.mesh.x) &  md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
+%	%Create MPCs to have periodic boundary conditions
+%	posx=find(md.mesh.x==0.);
+%	posx2=find(md.mesh.x==max(md.mesh.x));
+%	posx=find(md.mesh.x==0. & md.mesh.y~=0. & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
+%	posx2=find(md.mesh.x==max(md.mesh.x) &  md.mesh.y~=0. & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
 
-	%posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed); %Don't take the same nodes two times
-	%posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed);
+%	posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed); %Don't take the same nodes two times
+%	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed);
 
-	%md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
+%	md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
 
 	%Compute the diagnostic
 	md.diagnostic.abstol=NaN;
 	md.diagnostic.reltol=NaN;
-	md.diagnostic.restol=1;
+	md.diagnostic.restol=1.;
 	md.cluster=generic('name',oshostname(),'np',8);
 	md=solve(md,DiagnosticSolutionEnum());
@@ -62,45 +62,45 @@
 	%Now plot vx, vy, vz and vx on a cross section
 	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipastokesvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipastokesvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+		system(['mv ismipastokesvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 	end
 	plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipastokesvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipastokesvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+		system(['mv ismipastokesvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 	end
 	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4)
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipastokesvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipastokesvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+		system(['mv ismipastokesvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 	end
 
-	if(L==5000),
+	if(L==5000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','')
-	elseif(L==10000),
+	elseif(L==10000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','')
-	elseif(L==20000),
+	elseif(L==20000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','')
-	elseif(L==40000),
+	elseif(L==40000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','')
-	elseif(L==80000),
+	elseif(L==80000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','')
-	elseif(L==160000),
+	elseif(L==160000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
 	end
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipastokesvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipastokesvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+		system(['mv ismipastokesvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 	end
 end
@@ -108,15 +108,16 @@
 %Now plot the min and max values of vx for each size of the square
 plot([5 10 20 40 80 160],minvx);ylim([0 18])
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('ismipastokesminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipastokesminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+	system(['mv ismipastokesminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 end
 plot([5 10 20 40 80 160],maxvx);ylim([0 120])
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('ismipastokesmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipastokesmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
+	system(['mv ismipastokesmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA']);
 end
+
 %Fields and tolerances to track changes
 field_names     ={...
@@ -127,5 +128,5 @@
 	'Vx80km','Vy80km','Vz80km',...
 	'Vx160km','Vy160km','Vz160km'
-}
+};
 field_tolerances={...
 	1e-12,1e-12,1e-12,...
Index: /issm/trunk-jpl/test/NightlyRun/test1102.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1102.py	(revision 14134)
+++ /issm/trunk-jpl/test/NightlyRun/test1102.py	(revision 14134)
@@ -0,0 +1,162 @@
+import numpy
+import shutil
+from model import *
+from squaremesh import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This test is a test from the ISMP-HOM Intercomparison project.
+Pattyn and Payne 2006
+"""
+
+printingflag=False
+
+L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
+results=[]
+minvx=[]
+maxvx=[]
+
+for L in L_list:
+	nx=20    #numberof nodes in x direction
+	ny=20
+	md=model()
+	md=squaremesh(md,L,L,nx,ny)
+	md=setmask(md,'','')    #ice sheet test
+
+#	#Find elements at the corner and extract model
+#	posnodes=numpy.nonzero(numpy.logical_and(numpy.logical_or(md.mesh.x==0.,md.mesh.x==numpy.max(md.mesh.x)),numpy.logical_or(md.mesh.y==0.,md.mesh.y==numpy.max(md.mesh.y))))
+#	a=numpy.nonzero(ismember(md.mesh.elements,posnodes))[0]
+#	elements=numpy.ones((md.mesh.numberofelements),int)
+#	elements[a]=0
+#	md.modelextract(elements)
+
+	md=parameterize(md,'../Par/ISMIPA.py')
+	md.extrude(10,1.)
+	md=setflowequation(md,'stokes','all')
+
+	#Create dirichlet on the bed only
+	pos=numpy.nonzero(md.mesh.vertexonbed)
+	md.diagnostic.spcvx[pos]=0.
+	md.diagnostic.spcvy[pos]=0.
+	md.diagnostic.spcvz[pos]=0.
+
+#	#Create MPCs to have periodic boundary conditions
+#	posx=numpy.nonzero(md.mesh.x==0.)[0]
+#	posx2=numpy.nonzero(md.mesh.x==numpy.max(md.mesh.x))[0]
+#	posx=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x==0.,md.mesh.y!=0.),numpy.logical_and(md.mesh.y!=numpy.max(md.mesh.y),numpy.logical_not(md.mesh.vertexonbed))))[0]
+#	posx2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x==numpy.max(md.mesh.x),md.mesh.y!=0.),numpy.logical_and(md.mesh.y!=numpy.max(md.mesh.y),numpy.logical_not(md.mesh.vertexonbed))))[0]
+
+#	posy=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==0.,md.mesh.x!=0.),numpy.logical_and(md.mesh.x!=numpy.max(md.mesh.x),numpy.logical_not(md.mesh.vertexonbed))))[0]    #Don't take the same nodes two times
+#	posy2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==numpy.max(md.mesh.y),md.mesh.x!=0.),numpy.logical_and(md.mesh.x!=numpy.max(md.mesh.x),numpy.logical_not(md.mesh.vertexonbed))))[0]
+
+#	md.diagnostic.vertex_pairing=numpy.vstack((numpy.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),numpy.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))
+
+	#Compute the diagnostic
+	md.diagnostic.abstol=float('NaN')
+	md.diagnostic.reltol=float('NaN')
+	md.diagnostic.restol=1.
+	md.cluster=generic('name',oshostname(),'np',8)
+	md=solve(md,DiagnosticSolutionEnum())
+
+	#Plot the results and save them
+	vx=md.results.DiagnosticSolution.Vx
+	vy=md.results.DiagnosticSolution.Vy
+	vz=md.results.DiagnosticSolution.Vz
+	pressure=md.results.DiagnosticSolution.Pressure
+	results.append(md.results.DiagnosticSolution)
+	minvx.append(numpy.min(vx[-md.mesh.numberofvertices2d:]))
+	maxvx.append(numpy.max(vx[-md.mesh.numberofvertices2d:]))
+
+	#Now plot vx, vy, vz and vx on a cross section
+#	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipastokesvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipastokesvx%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+#	plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipastokesvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipastokesvy%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+#	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipastokesvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipastokesvz%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+
+	if   (L==5000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','')
+	elif (L==10000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','')
+	elif (L==20000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','')
+	elif (L==40000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','')
+	elif (L==80000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','')
+	elif (L==160000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipastokesvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipastokesvxsec.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+
+#Now plot the min and max values of vx for each size of the square
+#plot([5 10 20 40 80 160],minvx);ylim([0 18])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipastokesminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipastokesminvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+#plot([5 10 20 40 80 160],maxvx);ylim([0 120])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipastokesmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipastokesmaxvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
+
+#Fields and tolerances to track changes
+field_names     =[\
+	'Vx5km','Vy5km','Vz5km',\
+	'Vx10km','Vy10km','Vz10km',\
+	'Vx20km','Vy20km','Vz20km',\
+	'Vx40km','Vy40km','Vz40km',\
+	'Vx80km','Vy80km','Vz80km',\
+	'Vx160km','Vy160km','Vz160km'
+]
+field_tolerances=[\
+	1e-12,1e-12,1e-12,\
+	1e-12,1e-12,1e-12,\
+	1e-12,1e-11,1e-12,\
+	1e-12,1e-11,1e-12,\
+	1e-12,1e-11,1e-12,\
+	1e-12,1e-11,1e-12,\
+]
+field_values=[]
+for result in results:
+	field_values=field_values+[\
+		result.Vx,\
+		result.Vy,\
+		result.Vz,\
+		]
Index: /issm/trunk-jpl/test/NightlyRun/test1103.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1103.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1103.m	(revision 14134)
@@ -1,8 +1,10 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %Pattyn and Payne 2006
 printingflag=false;
 
-L_list={5000,10000,20000,40000,80000,160000};
+L_list={5000.,10000.,20000.,40000.,80000.,160000.};
 results={};
+minvx=[];
+maxvx=[];
 
 for i=1:length(L_list),
@@ -23,13 +25,13 @@
 	md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
 	pos=find(md.mesh.vertexonbed);
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
+	md.diagnostic.spcvx(pos)=0.;
+	md.diagnostic.spcvy(pos)=0.;
 
 	%Create MPCs to have periodic boundary conditions
-	posx=find(md.mesh.x==0);
+	posx=find(md.mesh.x==0.);
 	posx2=find(md.mesh.x==max(md.mesh.x));
 
-	posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
-	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
+	posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
+	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
 
 	md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
@@ -49,39 +51,39 @@
 	%Now plot vx, vy, vz and vx on a cross section
 	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipbpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipbpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
+		system(['mv ismipbpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
 	end
 	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipbpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipbpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
+		system(['mv ismipbpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
 	end
 
-	if(L==5000),
+	if(L==5000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[6 16],'xlim',[0 5000],'title','','xlabel','')
-	elseif(L==10000),
+	elseif(L==10000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 40],'xlim',[0 10000],'title','','xlabel','')
-	elseif(L==20000),
+	elseif(L==20000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 60],'xlim',[0 20000],'title','','xlabel','')
-	elseif(L==40000),
+	elseif(L==40000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 100],'xlim',[0 40000],'title','','xlabel','')
-	elseif(L==80000),
+	elseif(L==80000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 120],'xlim',[0 80000],'title','','xlabel','')
-	elseif(L==160000),
+	elseif(L==160000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
 	end
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipbpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipbpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
+		system(['mv ismipbpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
 	end
 end
@@ -89,15 +91,16 @@
 %Now plot the min and max values of vx for each size of the square
 plot([5 10 20 40 80 160],minvx);ylim([0 14]);xlim([0 160])
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('ismipbpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipbpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
+	system(['mv ismipbpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
 end
 plot([5 10 20 40 80 160],maxvx);ylim([0 120]);xlim([0 160])
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('ismipbpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipbpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB ']);
+	system(['mv ismipbpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestB']);
 end
+
 %Fields and tolerances to track changes
 field_names     ={...
Index: /issm/trunk-jpl/test/NightlyRun/test1103.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1103.py	(revision 14134)
+++ /issm/trunk-jpl/test/NightlyRun/test1103.py	(revision 14134)
@@ -0,0 +1,145 @@
+import numpy
+import shutil
+from model import *
+from squaremesh import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This test is a test from the ISMP-HOM Intercomparison project.
+Pattyn and Payne 2006
+"""
+
+printingflag=False
+
+L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
+results=[]
+minvx=[]
+maxvx=[]
+
+for L in L_list:
+	nx=20    #numberof nodes in x direction
+	ny=20
+	md=model()
+	md=squaremesh(md,L,L,nx,ny)
+	md=setmask(md,'','')    #ice sheet test
+	md=parameterize(md,'../Par/ISMIPB.py')
+	md.extrude(10,1.)
+
+	md=setflowequation(md,'pattyn','all')
+
+	#Create dirichlet on the bed only
+	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))
+	pos=numpy.nonzero(md.mesh.vertexonbed)
+	md.diagnostic.spcvx[pos]=0.
+	md.diagnostic.spcvy[pos]=0.
+
+	#Create MPCs to have periodic boundary conditions
+	posx=numpy.nonzero(md.mesh.x==0.)[0]
+	posx2=numpy.nonzero(md.mesh.x==numpy.max(md.mesh.x))[0]
+
+	posy=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==0.,md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]    #Don't take the same nodes two times
+	posy2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==numpy.max(md.mesh.y),md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]
+
+	md.diagnostic.vertex_pairing=numpy.vstack((numpy.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),numpy.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))
+
+	#Compute the diagnostic
+	md.cluster=generic('name',oshostname(),'np',8)
+	md=solve(md,DiagnosticSolutionEnum())
+
+	#Plot the results and save them
+	vx=md.results.DiagnosticSolution.Vx
+	vy=md.results.DiagnosticSolution.Vy
+	vz=md.results.DiagnosticSolution.Vz
+	results.append(md.results.DiagnosticSolution)
+	minvx.append(numpy.min(vx[md.mesh.numberofvertices2d:]))
+	maxvx.append(numpy.max(vx[md.mesh.numberofvertices2d:]))
+
+	#Now plot vx, vy, vz and vx on a cross section
+#	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipbpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipbpattynvx%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestB')
+#	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipbpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipbpattynvz%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestB')
+
+	if   (L==5000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[6 16],'xlim',[0 5000],'title','','xlabel','')
+	elif (L==10000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 40],'xlim',[0 10000],'title','','xlabel','')
+	elif (L==20000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 60],'xlim',[0 20000],'title','','xlabel','')
+	elif (L==40000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 100],'xlim',[0 40000],'title','','xlabel','')
+	elif (L==80000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 120],'xlim',[0 80000],'title','','xlabel','')
+	elif (L==160000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipbpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipbpattynvxsec%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestB')
+
+#Now plot the min and max values of vx for each size of the square
+#plot([5 10 20 40 80 160],minvx);ylim([0 14]);xlim([0 160])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipbpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipbpattynminvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestB')
+#plot([5 10 20 40 80 160],maxvx);ylim([0 120]);xlim([0 160])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipbpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipbpattynmaxvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestB')
+
+#Fields and tolerances to track changes
+field_names     =[\
+	'Vx5km','Vy5km','Vz5km',\
+	'Vx10km','Vy10km','Vz10km',\
+	'Vx20km','Vy20km','Vz20km',\
+	'Vx40km','Vy40km','Vz40km',\
+	'Vx80km','Vy80km','Vz80km',\
+	'Vx160km','Vy160km','Vz160km'
+]
+field_tolerances=[\
+	1e-09,1e-09,1e-09,\
+	1e-09,1e-09,1e-09,\
+	1e-09,1e-09,1e-09,\
+	1e-08,1e-08,1e-08,\
+	1e-08,1e-07,1e-07,\
+	1e-07,1e-06,1e-07,\
+]
+field_values=[]
+for result in results:
+	field_values=field_values+[\
+		result.Vx,\
+		result.Vy,\
+		result.Vz,\
+		]
Index: /issm/trunk-jpl/test/NightlyRun/test1104.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1104.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1104.m	(revision 14134)
@@ -1,6 +1,6 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %Pattyn and Payne 2006
 
-L_list={5000,10000,20000,40000,80000,160000};
+L_list={5000.,10000.,20000.,40000.,80000.,160000.};
 results={};
 
@@ -22,13 +22,13 @@
 
 	pos=find(md.mesh.vertexonbed);
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
+	md.diagnostic.spcvx(pos)=0.;
+	md.diagnostic.spcvy(pos)=0.;
 
 	%Create MPCs to have periodic boundary conditions
-	posx=find(md.mesh.x==0);
+	posx=find(md.mesh.x==0.);
 	posx2=find(md.mesh.x==max(md.mesh.x));
 
-	posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
-	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
+	posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
+	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
 
 	md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
@@ -38,5 +38,5 @@
 	md.cluster=generic('name',oshostname(),'np',8);
 	md=solve(md,DiagnosticSolutionEnum());
-	pos=find(md.mesh.x==0 | md.mesh.y==0 | md.mesh.x==max(md.mesh.x) | md.mesh.y==max(md.mesh.y));
+	pos=find(md.mesh.x==0. | md.mesh.y==0. | md.mesh.x==max(md.mesh.x) | md.mesh.y==max(md.mesh.y));
 	md.diagnostic.spcvx(pos)=md.results.DiagnosticSolution.Vx(pos);
 	md.diagnostic.spcvy(pos)=md.results.DiagnosticSolution.Vy(pos);
Index: /issm/trunk-jpl/test/NightlyRun/test1104.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1104.py	(revision 14134)
+++ /issm/trunk-jpl/test/NightlyRun/test1104.py	(revision 14134)
@@ -0,0 +1,89 @@
+import numpy
+from model import *
+from squaremesh import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This test is a test from the ISMP-HOM Intercomparison project.
+Pattyn and Payne 2006
+"""
+
+L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
+results=[]
+
+for L in L_list:
+	nx=20    #numberof nodes in x direction
+	ny=20
+	md=model()
+	md=squaremesh(md,L,L,nx,ny)
+	md=setmask(md,'','')    #ice sheet test
+	md=parameterize(md,'../Par/ISMIPB.py')
+	md.extrude(10,1.)
+	md=setflowequation(md,'pattyn','all')
+
+	#Create dirichlet on the bed only
+	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))
+
+	pos=numpy.nonzero(md.mesh.vertexonbed)
+	md.diagnostic.spcvx[pos]=0.
+	md.diagnostic.spcvy[pos]=0.
+
+	#Create MPCs to have periodic boundary conditions
+	posx=numpy.nonzero(md.mesh.x==0.)[0]
+	posx2=numpy.nonzero(md.mesh.x==numpy.max(md.mesh.x))[0]
+
+	posy=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==0.,md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]    #Don't take the same nodes two times
+	posy2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==numpy.max(md.mesh.y),md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]
+
+	md.diagnostic.vertex_pairing=numpy.vstack((numpy.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),numpy.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))
+
+	#Compute the diagnostic
+	md.diagnostic.abstol=float('NaN')
+	md.cluster=generic('name',oshostname(),'np',8)
+	md=solve(md,DiagnosticSolutionEnum())
+	pos=numpy.nonzero(numpy.logical_or(numpy.logical_or(md.mesh.x==0.,md.mesh.y==0.),numpy.logical_or(md.mesh.x==numpy.max(md.mesh.x),md.mesh.y==numpy.max(md.mesh.y))))
+	md.diagnostic.spcvx[pos]=md.results.DiagnosticSolution.Vx[pos]
+	md.diagnostic.spcvy[pos]=md.results.DiagnosticSolution.Vy[pos]
+	md.diagnostic.vertex_pairing=numpy.empty((0,2),int)
+	md=setflowequation(md,'stokes','all')
+	md=solve(md,DiagnosticSolutionEnum())
+
+	#Plot the results and save them
+	vx=md.results.DiagnosticSolution.Vx
+	vy=md.results.DiagnosticSolution.Vy
+	vz=md.results.DiagnosticSolution.Vz
+	results.append(md.results.DiagnosticSolution)
+
+#	plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.mesh.numberoflayers)
+
+#Fields and tolerances to track changes
+field_names     =[\
+	'Vx5km','Vy5km','Vz5km',\
+	'Vx10km','Vy10km','Vz10km',\
+	'Vx20km','Vy20km','Vz20km',\
+	'Vx40km','Vy40km','Vz40km',\
+	'Vx80km','Vy80km','Vz80km',\
+	'Vx160km','Vy160km','Vz160km'
+]
+field_tolerances=[\
+	1e-08,1e-08,1e-08,\
+	1e-08,1e-08,1e-08,\
+	1e-08,1e-08,1e-08,\
+	1e-08,1e-08,1e-08,\
+	1e-08,1e-07,1e-08,\
+	1e-07,1e-07,1e-07,\
+]
+field_values=[]
+for result in results:
+	field_values=field_values+[\
+		result.Vx,\
+		result.Vy,\
+		result.Vz,\
+		]
Index: /issm/trunk-jpl/test/NightlyRun/test1105.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1105.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1105.m	(revision 14134)
@@ -1,7 +1,7 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %Pattyn and Payne 2006
 printingflag=false;
 
-L_list={5000,10000,20000,40000,80000,160000};
+L_list={5000.,10000.,20000.,40000.,80000.,160000.};
 results={};
 minvx=[];
@@ -9,5 +9,5 @@
 
 for i=1:length(L_list),
-	L=L_list{i};  %in m (3 times the desired lenght for BC problems)  
+	L=L_list{i};  %in m (3 times the desired length for BC problems)  
 	nx=30; %number of nodes in x direction
 	ny=30;
@@ -18,5 +18,5 @@
 	md=extrude(md,10,1.);
 
-	md=setflowequation(md,'pattyn','all'); 
+	md=setflowequation(md,'pattyn','all');
 
 	%Create MPCs to have periodic boundary conditions
@@ -25,32 +25,32 @@
 	md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
 
-	posx=find(md.mesh.x==0 & md.mesh.y~=0 & md.mesh.y~=L);
-	posx2=find(md.mesh.x==L & md.mesh.y~=0 & md.mesh.y~=L);
+	posx=find(md.mesh.x==0. & md.mesh.y~=0. & md.mesh.y~=L);
+	posx2=find(md.mesh.x==L & md.mesh.y~=0. & md.mesh.y~=L);
 
-	posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=L); %Don't take the same nodes two times
-	posy2=find(md.mesh.y==L & md.mesh.x~=0 & md.mesh.x~=L);
+	posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=L); %Don't take the same nodes two times
+	posy2=find(md.mesh.y==L & md.mesh.x~=0. & md.mesh.x~=L);
 
 	md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
 
 	%Add spc on the corners
-	pos=find((md.mesh.x==0 | md.mesh.x==L) & (md.mesh.y==0 | md.mesh.y==L) & md.mesh.vertexonbed);
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
-	if(L==5000),
+	pos=find((md.mesh.x==0. | md.mesh.x==L) & (md.mesh.y==0. | md.mesh.y==L) & md.mesh.vertexonbed);
+	md.diagnostic.spcvx(pos)=0.;
+	md.diagnostic.spcvy(pos)=0.;
+	if(L==5000.),
 		md.diagnostic.spcvx(pos)=15.66;
 		md.diagnostic.spcvy(pos)=-0.1967;
-	elseif(L==10000),
+	elseif(L==10000.),
 		md.diagnostic.spcvx(pos)=16.04;
 		md.diagnostic.spcvy(pos)=-0.1977;
-	elseif(L==20000),
+	elseif(L==20000.),
 		md.diagnostic.spcvx(pos)=16.53;
 		md.diagnostic.spcvy(pos)=-1.27;
-	elseif(L==40000),
+	elseif(L==40000.),
 		md.diagnostic.spcvx(pos)=17.23;
 		md.diagnostic.spcvy(pos)=-3.17;
-	elseif(L==80000),
+	elseif(L==80000.),
 		md.diagnostic.spcvx(pos)=16.68;
 		md.diagnostic.spcvy(pos)=-2.69;
-	elseif(L==160000),
+	elseif(L==160000.),
 		md.diagnostic.spcvx(pos)=16.03;
 		md.diagnostic.spcvy(pos)=-1.27;
@@ -59,5 +59,5 @@
 	%Spc the bed at zero for vz
 	pos=find(md.mesh.vertexonbed);
-	md.diagnostic.spcvz(pos)=0;
+	md.diagnostic.spcvz(pos)=0.;
 
 	%Compute the diagnostic
@@ -75,45 +75,45 @@
 	%Now plot vx, vy, vz and vx on a cross section
 	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipcpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipcpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
+		system(['mv ismipcpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
 	end
 	plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipcpattynvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipcpattynvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
+		system(['mv ismipcpattynvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
 	end
 	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4)
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipcpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipcpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
+		system(['mv ismipcpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
 	end
 
-	if(L==5000),
+	if(L==5000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',5)
-	elseif(L==10000),
+	elseif(L==10000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[13 18],'xlim',[0 10000],'title','','xlabel','')
-	elseif(L==20000),
+	elseif(L==20000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[14 22],'xlim',[0 20000],'title','','xlabel','')
-	elseif(L==40000),
+	elseif(L==40000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[10 40],'xlim',[0 40000],'title','','xlabel','')
-	elseif(L==80000),
+	elseif(L==80000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 80],'xlim',[0 80000],'title','','xlabel','')
-	elseif(L==160000),
+	elseif(L==160000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 200],'xlim',[0 160000],'title','','xlabel','')
 	end
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipcpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipcpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
+		system(['mv ismipcpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
 	end
 end
@@ -121,15 +121,16 @@
 %Now plot the min and max values of vx for each size of the square
 plot([5 10 20 40 80 160],minvx);ylim([4 18]);xlim([0 160])
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('ismipcpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipcpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
+	system(['mv ismipcpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
 end
 plot([5 10 20 40 80 160],maxvx);ylim([0 200]); xlim([0 160])
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('ismipcpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipcpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC ']);
+	system(['mv ismipcpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestC']);
 end
+
 %Fields and tolerances to track changes
 field_names     ={...
Index: /issm/trunk-jpl/test/NightlyRun/test1105.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1105.py	(revision 14134)
+++ /issm/trunk-jpl/test/NightlyRun/test1105.py	(revision 14134)
@@ -0,0 +1,174 @@
+import numpy
+import shutil
+from model import *
+from squaremesh import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This test is a test from the ISMP-HOM Intercomparison project.
+Pattyn and Payne 2006
+"""
+
+printingflag=False
+
+L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
+results=[]
+minvx=[]
+maxvx=[]
+
+for L in L_list:    #in m (3 times the desired length for BC problems)  
+	nx=30    #number of nodes in x direction
+	ny=30
+	md=model()
+	md=squaremesh(md,L,L,nx,ny)
+	md=setmask(md,'','')    #ice sheet test
+	md=parameterize(md,'../Par/ISMIPC.py')
+	md.extrude(10,1.)
+
+	md=setflowequation(md,'pattyn','all')
+
+	#Create MPCs to have periodic boundary conditions
+	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))
+
+	posx=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x==0.,md.mesh.y!=0.),md.mesh.y!=L))[0]
+	posx2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.x==L,md.mesh.y!=0.),md.mesh.y!=L))[0]
+
+	posy=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==0.,md.mesh.x!=0.),md.mesh.x!=L))[0]    #Don't take the same nodes two times
+	posy2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==L,md.mesh.x!=0.),md.mesh.x!=L))[0]
+
+	md.diagnostic.vertex_pairing=numpy.vstack((numpy.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),numpy.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))
+
+	#Add spc on the corners
+	pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(numpy.logical_or(md.mesh.x==0.,md.mesh.x==L),numpy.logical_or(md.mesh.y==0.,md.mesh.y==L)),md.mesh.vertexonbed))
+	md.diagnostic.spcvx[pos]=0.
+	md.diagnostic.spcvy[pos]=0.
+	if   (L==5000.):
+		md.diagnostic.spcvx[pos]=15.66
+		md.diagnostic.spcvy[pos]=-0.1967
+	elif (L==10000.):
+		md.diagnostic.spcvx[pos]=16.04
+		md.diagnostic.spcvy[pos]=-0.1977
+	elif (L==20000.):
+		md.diagnostic.spcvx[pos]=16.53
+		md.diagnostic.spcvy[pos]=-1.27
+	elif (L==40000.):
+		md.diagnostic.spcvx[pos]=17.23
+		md.diagnostic.spcvy[pos]=-3.17
+	elif (L==80000.):
+		md.diagnostic.spcvx[pos]=16.68
+		md.diagnostic.spcvy[pos]=-2.69
+	elif (L==160000.):
+		md.diagnostic.spcvx[pos]=16.03
+		md.diagnostic.spcvy[pos]=-1.27
+	
+	#Spc the bed at zero for vz
+	pos=numpy.nonzero(md.mesh.vertexonbed)
+	md.diagnostic.spcvz[pos]=0.
+
+	#Compute the diagnostic
+	md.cluster=generic('name',oshostname(),'np',8)
+	md=solve(md,DiagnosticSolutionEnum())
+
+	#Plot the results and save them
+	vx=md.results.DiagnosticSolution.Vx
+	vy=md.results.DiagnosticSolution.Vy
+	vz=md.results.DiagnosticSolution.Vz
+	results.append(md.results.DiagnosticSolution)
+	minvx.append(numpy.min(vx[-md.mesh.numberofvertices2d:]))
+	maxvx.append(numpy.max(vx[-md.mesh.numberofvertices2d:]))
+
+	#Now plot vx, vy, vz and vx on a cross section
+#	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipcpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipcpattynvx%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
+#	plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipcpattynvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipcpattynvy%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
+#	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipcpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipcpattynvz%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
+
+	if   (L==5000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',5)
+	elif (L==10000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[13 18],'xlim',[0 10000],'title','','xlabel','')
+	elif (L==20000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[14 22],'xlim',[0 20000],'title','','xlabel','')
+	elif (L==40000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[10 40],'xlim',[0 40000],'title','','xlabel','')
+	elif (L==80000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 80],'xlim',[0 80000],'title','','xlabel','')
+	elif (L==160000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 200],'xlim',[0 160000],'title','','xlabel','')
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipcpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipcpattynvxsec%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
+
+#Now plot the min and max values of vx for each size of the square
+#plot([5 10 20 40 80 160],minvx);ylim([4 18]);xlim([0 160])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipcpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipcpattynminvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
+#plot([5 10 20 40 80 160],maxvx);ylim([0 200]); xlim([0 160])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipcpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipcpattynmaxvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestC')
+
+#Fields and tolerances to track changes
+field_names     =[\
+	'Vx5km','Vy5km','Vz5km',\
+	'Vx10km','Vy10km','Vz10km',\
+	'Vx20km','Vy20km','Vz20km',\
+	'Vx40km','Vy40km','Vz40km',\
+	'Vx80km','Vy80km','Vz80km',\
+	'Vx160km','Vy160km','Vz160km'
+]
+field_tolerances=[\
+	1e-08,1e-07,1e-07,\
+	1e-09,1e-07,1e-07,\
+	1e-09,1e-09,1e-07,\
+	1e-09,1e-09,1e-08,\
+	1e-09,1e-08,1e-08,\
+	1e-09,1e-08,1e-08,\
+]
+field_values=[]
+for result in results:
+	field_values=field_values+[\
+		result.Vx,\
+		result.Vy,\
+		result.Vz,\
+		]
Index: /issm/trunk-jpl/test/NightlyRun/test1106.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1106.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1106.m	(revision 14134)
@@ -1,6 +1,6 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %Pattyn and Payne 2006
 
-L_list={5000,10000,20000,40000,80000,160000};
+L_list={5000.,10000.,20000.,40000.,80000.,160000.};
 results={};
 
@@ -10,32 +10,32 @@
 	md=setmask(md,'',''); %ice sheet test
 	md=parameterize(md,'../Par/ISMIPC.par');
-	md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/L).*sin(md.mesh.y*2*pi/L)));
+	md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/L).*sin(md.mesh.y*2.*pi/L)));
 	md=extrude(md,10,1.);
 
 	%Add spc on the borders
-	pos=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x) | md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
-	if(L==5000),
+	pos=find(md.mesh.x==0. | md.mesh.x==max(md.mesh.x) | md.mesh.y==0. | md.mesh.y==max(md.mesh.y));
+	md.diagnostic.spcvx(pos)=0.;
+	md.diagnostic.spcvy(pos)=0.;
+	if(L==5000.),
 		md.diagnostic.spcvx(pos)=15.66;
 		md.diagnostic.spcvy(pos)=-0.1967;
-	elseif(L==10000),
+	elseif(L==10000.),
 		md.diagnostic.spcvx(pos)=16.04;
 		md.diagnostic.spcvy(pos)=-0.1977;
-	elseif(L==20000),
+	elseif(L==20000.),
 		md.diagnostic.spcvx(pos)=16.53;
 		md.diagnostic.spcvy(pos)=-1.27;
-	elseif(L==40000),
+	elseif(L==40000.),
 		md.diagnostic.spcvx(pos)=17.23;
 		md.diagnostic.spcvy(pos)=-3.17;
-	elseif(L==80000),
+	elseif(L==80000.),
 		md.diagnostic.spcvx(pos)=16.68;
 		md.diagnostic.spcvy(pos)=-2.69;
-	elseif(L==160000),
+	elseif(L==160000.),
 		md.diagnostic.spcvx(pos)=16.03;
 		md.diagnostic.spcvy(pos)=-1.27;
 	end
 
-	md=setflowequation(md,'stokes','all'); 
+	md=setflowequation(md,'stokes','all');
 
 	%Compute the diagnostic
Index: /issm/trunk-jpl/test/NightlyRun/test1106.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1106.py	(revision 14134)
+++ /issm/trunk-jpl/test/NightlyRun/test1106.py	(revision 14134)
@@ -0,0 +1,86 @@
+import numpy
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This test is a test from the ISMP-HOM Intercomparison project.
+Pattyn and Payne 2006
+"""
+
+L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
+results=[]
+
+for L in L_list:
+	md=triangle(model(),"../Exp/Square_%d.exp" % L,L/10.)    #size 3*L 
+	md=setmask(md,'','')    #ice sheet test
+	md=parameterize(md,'../Par/ISMIPC.py')
+	md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/L)*numpy.sin(md.mesh.y.reshape(-1,1)*2.*numpy.pi/L)))
+	md.extrude(10,1.)
+
+	#Add spc on the borders
+	pos=numpy.nonzero(numpy.logical_or(numpy.logical_or(md.mesh.x==0.,md.mesh.x==numpy.max(md.mesh.x)),numpy.logical_or(md.mesh.y==0.,md.mesh.y==numpy.max(md.mesh.y))))
+	md.diagnostic.spcvx[pos]=0.
+	md.diagnostic.spcvy[pos]=0.
+	if   (L==5000.):
+		md.diagnostic.spcvx[pos]=15.66
+		md.diagnostic.spcvy[pos]=-0.1967
+	elif (L==10000.):
+		md.diagnostic.spcvx[pos]=16.04
+		md.diagnostic.spcvy[pos]=-0.1977
+	elif (L==20000.):
+		md.diagnostic.spcvx[pos]=16.53
+		md.diagnostic.spcvy[pos]=-1.27
+	elif (L==40000.):
+		md.diagnostic.spcvx[pos]=17.23
+		md.diagnostic.spcvy[pos]=-3.17
+	elif (L==80000.):
+		md.diagnostic.spcvx[pos]=16.68
+		md.diagnostic.spcvy[pos]=-2.69
+	elif (L==160000.):
+		md.diagnostic.spcvx[pos]=16.03
+		md.diagnostic.spcvy[pos]=-1.27
+
+	md=setflowequation(md,'stokes','all')
+
+	#Compute the diagnostic
+	md.cluster=generic('name',oshostname(),'np',8)
+	md=solve(md,DiagnosticSolutionEnum())
+
+	#Plot the results and save them
+	vx=md.results.DiagnosticSolution.Vx
+	vy=md.results.DiagnosticSolution.Vy
+	vz=md.results.DiagnosticSolution.Vz
+	results.append(md.results.DiagnosticSolution)
+
+#	plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.mesh.numberoflayers)
+
+#Fields and tolerances to track changes
+field_names     =[\
+	'Vx5km','Vy5km','Vz5km',\
+	'Vx10km','Vy10km','Vz10km',\
+	'Vx20km','Vy20km','Vz20km',\
+	'Vx40km','Vy40km','Vz40km',\
+	'Vx80km','Vy80km','Vz80km',\
+	'Vx160km','Vy160km','Vz160km'
+]
+field_tolerances=[\
+	1e-12,1e-12,1e-11,\
+	1e-12,1e-12,1e-12,\
+	1e-12,1e-12,1e-12,\
+	1e-12,1e-12,1e-12,\
+	1e-12,1e-12,1e-12,\
+	1e-12,1e-11,1e-12,\
+]
+field_values=[]
+for result in results:
+	field_values=field_values+[\
+		result.Vx,\
+		result.Vy,\
+		result.Vz,\
+		]
Index: /issm/trunk-jpl/test/NightlyRun/test1107.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1107.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1107.m	(revision 14134)
@@ -1,7 +1,7 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %Pattyn and Payne 2006
 printingflag=false;
 
-L_list={5000,10000,20000,40000,80000,160000};
+L_list={5000.,10000.,20000.,40000.,80000.,160000.};
 results={};
 minvx=[];
@@ -26,27 +26,27 @@
 
 	%Create MPCs to have periodic boundary conditions
-	posx=find(md.mesh.x==0 & ~(md.mesh.y==0 & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
-	posx2=find(md.mesh.x==max(md.mesh.x) & ~(md.mesh.y==0 & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
+	posx=find(md.mesh.x==0. & ~(md.mesh.y==0. & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
+	posx2=find(md.mesh.x==max(md.mesh.x) & ~(md.mesh.y==0. & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
 
-	posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
-	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
+	posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
+	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
 
 	md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
 
 	%Add spc on the corners
-	pos=find((md.mesh.x==0 | md.mesh.x==L) & (md.mesh.y==0 | md.mesh.y==L) & md.mesh.vertexonbed);
-	md.diagnostic.spcvy(:)=0;
-	md.diagnostic.spcvx(pos)=0;
-	if(L==5000),
+	pos=find((md.mesh.x==0. | md.mesh.x==L) & (md.mesh.y==0. | md.mesh.y==L) & md.mesh.vertexonbed);
+	md.diagnostic.spcvy(:)=0.;
+	md.diagnostic.spcvx(pos)=0.;
+	if(L==5000.),
 		md.diagnostic.spcvx(pos)=16.0912;
-	elseif(L==10000),
+	elseif(L==10000.),
 		md.diagnostic.spcvx(pos)=16.52;
-	elseif(L==20000),
+	elseif(L==20000.),
 		md.diagnostic.spcvx(pos)=17.77;
-	elseif(L==40000),
+	elseif(L==40000.),
 		md.diagnostic.spcvx(pos)=19.88;
-	elseif(L==80000),
+	elseif(L==80000.),
 		md.diagnostic.spcvx(pos)=18.65;
-	elseif(L==160000),
+	elseif(L==160000.),
 		md.diagnostic.spcvx(pos)=16.91;
 	end
@@ -54,5 +54,5 @@
 	%Spc the bed at zero for vz
 	pos=find(md.mesh.vertexonbed);
-	md.diagnostic.spcvz(pos)=0;
+	md.diagnostic.spcvz(pos)=0.;
 
 	%Compute the diagnostic
@@ -70,39 +70,39 @@
 	%Now plot vx, vy, vz and vx on a cross section
 	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipdpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipdpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
+		system(['mv ismipdpattynvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
 	end
 	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipdpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipdpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
+		system(['mv ismipdpattynvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
 	end
 
-	if(L==5000),
+	if(L==5000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',4)
-	elseif(L==10000),
+	elseif(L==10000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 20],'xlim',[0 10000],'title','','xlabel','','figure',4)
-	elseif(L==20000),
+	elseif(L==20000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 30],'xlim',[0 20000],'title','','xlabel','','figure',4)
-	elseif(L==40000),
+	elseif(L==40000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[10 60],'xlim',[0 40000],'title','','xlabel','','figure',4)
-	elseif(L==80000),
+	elseif(L==80000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 200],'xlim',[0 80000],'title','','xlabel','','figure',4)
-	elseif(L==160000),
+	elseif(L==160000.),
 		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
 			'resolution',[10 10],'ylim',[0 400],'xlim',[0 160000],'title','','xlabel','','figure',4)
 	end
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		printmodel(['ismipdpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-		system(['mv ismipdpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
+		system(['mv ismipdpattynvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
 	end
 end
@@ -110,15 +110,16 @@
 %Now plot the min and max values of vx for each size of the square
 plot([5 10 20 40 80 160],minvx);ylim([2 18]);xlim([0 160])
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('ismipdpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipdpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
+	system(['mv ismipdpattynminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
 end
 plot([5 10 20 40 80 160],maxvx);ylim([0 300]);xlim([0 160])
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('ismipdpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-	system(['mv ismipdpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD ']);
+	system(['mv ismipdpattynmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestD']);
 end
+
 %Fields and tolerances to track changes
 field_names     ={...
Index: /issm/trunk-jpl/test/NightlyRun/test1107.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1107.py	(revision 14134)
+++ /issm/trunk-jpl/test/NightlyRun/test1107.py	(revision 14134)
@@ -0,0 +1,165 @@
+import numpy
+import shutil
+from model import *
+from squaremesh import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This test is a test from the ISMP-HOM Intercomparison project.
+Pattyn and Payne 2006
+"""
+
+printingflag=False
+
+L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
+results=[]
+minvx=[]
+maxvx=[]
+
+for L in L_list:
+	nx=30    #numberof nodes in x direction
+	ny=30
+	md=model()
+	md=squaremesh(md,L,L,nx,ny)
+	md=setmask(md,'','')    #ice sheet test
+	md=parameterize(md,'../Par/ISMIPD.py')
+	md.extrude(10,1.)
+
+	md=setflowequation(md,'pattyn','all')
+
+	#We need one grd on dirichlet: the 4 corners are set to zero
+	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))
+
+	#Create MPCs to have periodic boundary conditions
+#	posx=find(md.mesh.x==0. & ~(md.mesh.y==0. & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
+	posx=numpy.nonzero(numpy.logical_and(md.mesh.x==0.,numpy.logical_and(numpy.logical_not(numpy.logical_and(md.mesh.y==0.,md.mesh.vertexonbed)),numpy.logical_not(numpy.logical_and(md.mesh.y==L,md.mesh.vertexonbed)))))[0]
+#	posx2=find(md.mesh.x==max(md.mesh.x) & ~(md.mesh.y==0. & md.mesh.vertexonbed) & ~(md.mesh.y==L & md.mesh.vertexonbed));
+	posx2=numpy.nonzero(numpy.logical_and(md.mesh.x==numpy.max(md.mesh.x),numpy.logical_and(numpy.logical_not(numpy.logical_and(md.mesh.y==0.,md.mesh.vertexonbed)),numpy.logical_not(numpy.logical_and(md.mesh.y==L,md.mesh.vertexonbed)))))[0]
+
+	posy=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==0.,md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]    #Don't take the same nodes two times
+	posy2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==numpy.max(md.mesh.y),md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]
+
+	md.diagnostic.vertex_pairing=numpy.vstack((numpy.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),numpy.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))
+
+	#Add spc on the corners
+	pos=numpy.nonzero(numpy.logical_and(numpy.logical_and(numpy.logical_or(md.mesh.x==0.,md.mesh.x==L),numpy.logical_or(md.mesh.y==0.,md.mesh.y==L)),md.mesh.vertexonbed))
+	md.diagnostic.spcvy[:]=0.
+	md.diagnostic.spcvx[pos]=0.
+	if   (L==5000.):
+		md.diagnostic.spcvx[pos]=16.0912
+	elif (L==10000.):
+		md.diagnostic.spcvx[pos]=16.52
+	elif (L==20000.):
+		md.diagnostic.spcvx[pos]=17.77
+	elif (L==40000.):
+		md.diagnostic.spcvx[pos]=19.88
+	elif (L==80000.):
+		md.diagnostic.spcvx[pos]=18.65
+	elif (L==160000.):
+		md.diagnostic.spcvx[pos]=16.91
+	
+	#Spc the bed at zero for vz
+	pos=numpy.nonzero(md.mesh.vertexonbed)
+	md.diagnostic.spcvz[pos]=0.
+
+	#Compute the diagnostic
+	md.cluster=generic('name',oshostname(),'np',8)
+	md=solve(md,DiagnosticSolutionEnum())
+
+	#Plot the results and save them
+	vx=md.results.DiagnosticSolution.Vx
+	vy=md.results.DiagnosticSolution.Vy
+	vz=md.results.DiagnosticSolution.Vz
+	results.append(md.results.DiagnosticSolution)
+	minvx.append(numpy.min(vx[-md.mesh.numberofvertices2d:]))
+	maxvx.append(numpy.max(vx[-md.mesh.numberofvertices2d:]))
+
+	#Now plot vx, vy, vz and vx on a cross section
+#	plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipdpattynvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipdpattynvx%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD')
+#	plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipdpattynvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipdpattynvz%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD')
+
+	if   (L==5000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 20],'xlim',[0 5000],'title','','xlabel','','figure',4)
+	elif (L==10000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 20],'xlim',[0 10000],'title','','xlabel','','figure',4)
+	elif (L==20000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 30],'xlim',[0 20000],'title','','xlabel','','figure',4)
+	elif (L==40000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[10 60],'xlim',[0 40000],'title','','xlabel','','figure',4)
+	elif (L==80000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 200],'xlim',[0 80000],'title','','xlabel','','figure',4)
+	elif (L==160000.):
+		pass
+#		plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
+#			'resolution',[10 10],'ylim',[0 400],'xlim',[0 160000],'title','','xlabel','','figure',4)
+	if printingflag:
+		pass
+#		set(gcf,'Color','w')
+#		printmodel(['ismipdpattynvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#		shutil.move("ismipdpattynvxsec%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD')
+
+#Now plot the min and max values of vx for each size of the square
+#plot([5 10 20 40 80 160],minvx);ylim([2 18]);xlim([0 160])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipdpattynminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipdpattynminvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD')
+#plot([5 10 20 40 80 160],maxvx);ylim([0 300]);xlim([0 160])
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('ismipdpattynmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
+#	shutil.move('ismipdpattynmaxvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestD')
+
+#Fields and tolerances to track changes
+field_names     =[\
+	'Vx5km','Vy5km','Vz5km',\
+	'Vx10km','Vy10km','Vz10km',\
+	'Vx20km','Vy20km','Vz20km',\
+	'Vx40km','Vy40km','Vz40km',\
+	'Vx80km','Vy80km','Vz80km',\
+	'Vx160km','Vy160km','Vz160km'
+]
+field_tolerances=[\
+	1e-07,1e-08,1e-06,\
+	1e-08,1e-08,1e-06,\
+	1e-08,1e-08,1e-07,\
+	1e-08,1e-08,1e-07,\
+	1e-08,1e-08,1e-07,\
+	1e-07,1e-08,1e-06,\
+]
+field_values=[]
+for result in results:
+	field_values=field_values+[\
+		result.Vx,\
+		result.Vy,\
+		result.Vz,\
+		]
Index: /issm/trunk-jpl/test/NightlyRun/test1108.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1108.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1108.m	(revision 14134)
@@ -1,6 +1,6 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %Pattyn and Payne 2006
 
-L_list={5000,10000,20000,40000,80000,160000};
+L_list={5000.,10000.,20000.,40000.,80000.,160000.};
 results={};
 
@@ -17,20 +17,20 @@
 	md=setflowequation(md,'pattyn','all');
 
-	%We need one grd on dirichlet: the 4 corners are set to zero
+	%We need one grid on dirichlet: the 4 corners are set to zero
 	md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
 	md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
 	md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
 	
-	pos=find(md.mesh.vertexonbed & (md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
-	md.diagnostic.spcvz(pos)=0;
+	pos=find(md.mesh.vertexonbed & (md.mesh.x==0. | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0. | md.mesh.y==max(md.mesh.y)));
+	md.diagnostic.spcvx(pos)=0.;
+	md.diagnostic.spcvy(pos)=0.;
+	md.diagnostic.spcvz(pos)=0.;
 
 	%Create MPCs to have periodic boundary conditions
-	posx=find(md.mesh.x==0);
+	posx=find(md.mesh.x==0.);
 	posx2=find(md.mesh.x==max(md.mesh.x));
 
-	posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
-	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
+	posy=find(md.mesh.y==0. & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x)); %Don't take the same nodes two times
+	posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0. & md.mesh.x~=max(md.mesh.x));
 
 	md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
@@ -47,5 +47,5 @@
 	md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
 	md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
-	pos=find(md.mesh.y==0 | md.mesh.x==0 | md.mesh.x==max(md.mesh.x) | md.mesh.y==max(md.mesh.y)); %Don't take the same nodes two times
+	pos=find(md.mesh.y==0. | md.mesh.x==0. | md.mesh.x==max(md.mesh.x) | md.mesh.y==max(md.mesh.y)); %Don't take the same nodes two times
 	md.diagnostic.spcvx(pos)=md.results.DiagnosticSolution.Vx(pos);
 	md.diagnostic.spcvy(pos)=md.results.DiagnosticSolution.Vy(pos);
Index: /issm/trunk-jpl/test/NightlyRun/test1108.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1108.py	(revision 14134)
+++ /issm/trunk-jpl/test/NightlyRun/test1108.py	(revision 14134)
@@ -0,0 +1,97 @@
+import numpy
+from model import *
+from bamg import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This test is a test from the ISMP-HOM Intercomparison project.
+Pattyn and Payne 2006
+"""
+
+L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
+results=[]
+
+for L in L_list:
+	nx=30    #numberof nodes in x direction
+	ny=30
+	md=model()
+	md=squaremesh(md,L,L,nx,ny)
+	md=setmask(md,'','')    #ice sheet test
+	md=parameterize(md,'../Par/ISMIPD.py')
+	md.extrude(10,1.)
+
+	md=setflowequation(md,'pattyn','all')
+
+	#We need one grd on dirichlet: the 4 corners are set to zero
+	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))
+	
+	pos=numpy.nonzero(numpy.logical_and(md.mesh.vertexonbed,numpy.logical_and(numpy.logical_or(md.mesh.x==0.,md.mesh.x==numpy.max(md.mesh.x)),numpy.logical_or(md.mesh.y==0.,md.mesh.y==numpy.max(md.mesh.y)))))
+	md.diagnostic.spcvx[pos]=0.
+	md.diagnostic.spcvy[pos]=0.
+	md.diagnostic.spcvz[pos]=0.
+
+	#Create MPCs to have periodic boundary conditions
+	posx=numpy.nonzero(md.mesh.x==0.)[0]
+	posx2=numpy.nonzero(md.mesh.x==numpy.max(md.mesh.x))[0]
+
+	posy=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==0.,md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]    #Don't take the same nodes two times
+	posy2=numpy.nonzero(numpy.logical_and(numpy.logical_and(md.mesh.y==numpy.max(md.mesh.y),md.mesh.x!=0.),md.mesh.x!=numpy.max(md.mesh.x)))[0]
+
+	md.diagnostic.vertex_pairing=numpy.vstack((numpy.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),numpy.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))
+
+	#Compute the diagnostic
+	md.cluster=generic('name',oshostname(),'np',8)
+	md.verbose=verbose('convergence',True)
+	md=solve(md,DiagnosticSolutionEnum())
+	md.diagnostic.reltol=float('NaN')
+	md.diagnostic.abstol=float('NaN')
+	md.diagnostic.vertex_pairing=numpy.empty((0,2))
+	#We need one grid on dirichlet: the 4 corners are set to zero
+	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))
+	pos=numpy.nonzero(numpy.logical_or(numpy.logical_or(md.mesh.y==0.,md.mesh.x==0.),numpy.logical_or(md.mesh.x==numpy.max(md.mesh.x),md.mesh.y==numpy.max(md.mesh.y))))    #Don't take the same nodes two times
+	md.diagnostic.spcvx[pos]=md.results.DiagnosticSolution.Vx[pos]
+	md.diagnostic.spcvy[pos]=md.results.DiagnosticSolution.Vy[pos]
+	md=setflowequation(md,'stokes','all')
+	md=solve(md,DiagnosticSolutionEnum())
+
+	#Plot the results and save them
+	vx=md.results.DiagnosticSolution.Vx
+	vy=md.results.DiagnosticSolution.Vy
+	vz=md.results.DiagnosticSolution.Vz
+	results.append(md.results.DiagnosticSolution)
+
+#	plotmodel(md,'data',vx,'data',vy,'data',vz,'layer#all',md.mesh.numberoflayers)
+
+#Fields and tolerances to track changes
+field_names     =[\
+	'Vx5km','Vy5km','Vz5km',\
+	'Vx10km','Vy10km','Vz10km',\
+	'Vx20km','Vy20km','Vz20km',\
+	'Vx40km','Vy40km','Vz40km',\
+	'Vx80km','Vy80km','Vz80km',\
+	'Vx160km','Vy160km','Vz160km'
+]
+field_tolerances=[\
+	1e-07,1e-07,1e-07,\
+	1e-08,1e-08,1e-08,\
+	1e-08,1e-07,1e-07,\
+	1e-08,1e-08,1e-08,\
+	1e-08,1e-07,1e-07,\
+	1e-07,1e-06,1e-07,\
+]
+field_values=[]
+for result in results:
+	field_values=field_values+[\
+		result.Vx,\
+		result.Vy,\
+		result.Vz,\
+		]
Index: /issm/trunk-jpl/test/NightlyRun/test1109.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1109.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1109.m	(revision 14134)
@@ -1,3 +1,3 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %TestE 
 %Four tests to run: - Pattyn frozen
@@ -9,6 +9,6 @@
 
 for i=1:4,
-	Lx=10; %in m
-	Ly=5000; %in m
+	Lx=10.; %in m
+	Ly=5000.; %in m
 	nx=3; %number of nodes in x direction
 	ny=51;
@@ -26,5 +26,5 @@
 
 	%Create MPCs to have periodic boundary conditions
-	posx=find(md.mesh.x==0);
+	posx=find(md.mesh.x==0.);
 	posx2=find(md.mesh.x==max(md.mesh.x));
 	md.diagnostic.vertex_pairing=[posx,posx2];
@@ -35,7 +35,7 @@
 	md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
 	md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
-	md.diagnostic.spcvx(pos)=0;
-	md.diagnostic.spcvy(pos)=0;
-	md.diagnostic.spcvz(pos)=0;
+	md.diagnostic.spcvx(pos)=0.;
+	md.diagnostic.spcvy(pos)=0.;
+	md.diagnostic.spcvz(pos)=0.;
 
 	%Remove the spc where there is some sliding (case 3 and 4):
@@ -58,34 +58,35 @@
 	if i==1,
 		plotmodel(md,'data',vy,'ylim',[-10 80],'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIPE.exp','resolution',[10 10],'title','','xlabel','')
-		if printingflag, 
+		if printingflag,
 			set(gcf,'Color','w')
 			printmodel('ismipepattynvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipepattynvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE ']);
+			system(['mv ismipepattynvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE']);
 		end
 	elseif i==2,
 		plotmodel(md,'data',vy,'ylim',[-10 80],'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIPE.exp','resolution',[10 10],'title','','xlabel','')
-		if printingflag, 
+		if printingflag,
 			set(gcf,'Color','w')
 			printmodel('ismipestokesvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipestokesvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE ']);
+			system(['mv ismipestokesvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE']);
 		end
 	elseif i==3,
 		plotmodel(md,'data',vy,'ylim',[-50 200],'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIPE.exp','resolution',[10 10],'title','','xlabel','')
-		if printingflag, 
+		if printingflag,
 			set(gcf,'Color','w')
 			printmodel('ismipepattynvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipepattynvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE ']);
+			system(['mv ismipepattynvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE']);
 		end
 	elseif i==4,
 		plotmodel(md,'data',vy,'ylim',[-50 200],'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIPE.exp','resolution',[10 10],'title','','xlabel','')
-		if printingflag, 
+		if printingflag,
 			set(gcf,'Color','w')
 			printmodel('ismipestokesvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipestokesvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE ']);
+			system(['mv ismipestokesvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestE']);
 		end
 	end
 end
+
 %Fields and tolerances to track changes
-field_names     ={ ...
+field_names     ={...
 	'VyPattynSliding','VzPattynSliding',...
 	'VxStokesSliding','VyStokesSliding','VzStokesSliding',...
Index: /issm/trunk-jpl/test/NightlyRun/test1110.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1110.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1110.m	(revision 14134)
@@ -1,13 +1,14 @@
-%This test is a test from the ISMP-HOM Intercomparison project
+%This test is a test from the ISMP-HOM Intercomparison project.
 %TestF 
 printingflag=false;
+results={};
 
 for i=1:4,
-	L=100000; %in m
+	L=100000.; %in m
 	nx=30; %numberof nodes in x direction
 	ny=30;
 	md=model();
 	md=squaremesh(md,L,L,nx,ny);
-	%md=triangle(md,'../Exp/SquareISMIP.exp',5500.);
+%	md=triangle(md,'../Exp/SquareISMIP.exp',5500.);
 	md=setmask(md,'',''); %ice sheet test
 	md=parameterize(md,'../Par/ISMIPF.par');
@@ -26,21 +27,21 @@
 		%Create dirichlet on the bed if no slip
 		pos=find(md.mesh.vertexonbed);
-		md.diagnostic.spcvx(pos)=0;
-		md.diagnostic.spcvy(pos)=0;
-		md.diagnostic.spcvz(pos)=0;
+		md.diagnostic.spcvx(pos)=0.;
+		md.diagnostic.spcvy(pos)=0.;
+		md.diagnostic.spcvz(pos)=0.;
 	else
-		pos=find(md.mesh.vertexonbed & (md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
-		md.diagnostic.spcvx(pos)=100; %because we need a dirichlet somewhere
-		md.diagnostic.spcvy(pos)=0;
-		md.diagnostic.spcvz(pos)=0;
+		pos=find(md.mesh.vertexonbed & (md.mesh.x==0. | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0. | md.mesh.y==max(md.mesh.y)));
+		md.diagnostic.spcvx(pos)=100.; %because we need a dirichlet somewhere
+		md.diagnostic.spcvy(pos)=0.;
+		md.diagnostic.spcvz(pos)=0.;
 	end
 	pos=find(~md.mesh.vertexonbed);
-	md.thermal.spctemperature(pos)=255;
+	md.thermal.spctemperature(pos)=255.;
 
 	%Create MPCs to have periodic boundary conditions
-	posx=find(md.mesh.x==0);
+	posx=find(md.mesh.x==0.);
 	posx2=find(md.mesh.x==max(md.mesh.x));
 
-	posy=find(md.mesh.y==0);
+	posy=find(md.mesh.y==0.);
 	posy2=find(md.mesh.y==max(md.mesh.y));
 
@@ -48,6 +49,6 @@
 	md.prognostic.vertex_pairing=[posx,posx2;posy,posy2];
 
-	md.timestepping.time_step=3;
-	md.timestepping.final_time=300;
+	md.timestepping.time_step=3.;
+	md.timestepping.final_time=300.;
 	md.settings.output_frequency=50;
 	md.prognostic.stabilization=1;
@@ -68,36 +69,36 @@
 		plotmodel(md,'data',(md.results.TransientSolution(end).Vx),'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIP100000.exp','title','','xlabel','','ylabel','Velocity (m/yr)','linewidth',3,'grid','on','unit','km','ylim',[185 200])
 	end
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		if i==1,
 			printmodel('ismipfpattynvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipfpattynvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
+			system(['mv ismipfpattynvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
 		elseif i==2,
 			printmodel('ismipfpattynvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipfpattynvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
+			system(['mv ismipfpattynvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
 		elseif i==3,
 			printmodel('ismipfstokesvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipfstokesvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
+			system(['mv ismipfstokesvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
 		elseif i==4,
 			printmodel('ismipfstokesvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipfstokesvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
+			system(['mv ismipfstokesvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
 		end
 	end
 
 	plotmodel(md,'data',(md.results.TransientSolution(end).Surface)-md.geometry.surface,'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIP100000.exp','title','','xlabel','','ylabel','Surface (m)','linewidth',3,'grid','on','unit','km','ylim',[-30 50])
-	if printingflag, 
+	if printingflag,
 		set(gcf,'Color','w')
 		if i==1,
 			printmodel('ismipfpattyndeltasurfacefrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipfpattyndeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
+			system(['mv ismipfpattyndeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
 		elseif i==2,
 			printmodel('ismipfpattyndeltasurfacesliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipfpattyndeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
+			system(['mv ismipfpattyndeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
 		elseif i==3,
 			printmodel('ismipfstokesdeltasurfacefrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipfstokesdeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
+			system(['mv ismipfstokesdeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
 		elseif i==4,
 			printmodel('ismipfstokesdeltasurfacesliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
-			system(['mv ismipfstokesdeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF ']);
+			system(['mv ismipfstokesdeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF']);
 		end
 	end
@@ -105,5 +106,5 @@
 
 %Fields and tolerances to track changes
-field_names     ={ ...
+field_names     ={...
 	'VxPattynFrozen','VyPattynFrozen','VzPattynFrozen','SurfacePattynFrozen',...
 	'VxPattynSliding','VyPattynSliding','VzPattynSliding','SurfacePattynSliding',...
Index: /issm/trunk-jpl/test/NightlyRun/test1201.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1201.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1201.m	(revision 14134)
@@ -1,3 +1,3 @@
-%This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996
+%This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996.
 printingflag=false;
 
Index: /issm/trunk-jpl/test/NightlyRun/test1201.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1201.py	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1201.py	(revision 14134)
@@ -10,5 +10,5 @@
 
 """
-This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996
+This test is a test from the EISMINT for Ice shelves Vincent Rommelaere 1996.
 """
 
Index: /issm/trunk-jpl/test/NightlyRun/test1202.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1202.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1202.m	(revision 14134)
@@ -1,3 +1,3 @@
-%Test on the diagnostic model and the prognostic in 2d
+%Test on the diagnostic model and the prognostic in 2d.
 printingflag=false;
 
Index: /issm/trunk-jpl/test/NightlyRun/test1205.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1205.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1205.m	(revision 14134)
@@ -1,11 +1,11 @@
-%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
+%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling.
 printingflag=false;
 
 numlayers=10;
-resolution=30000;
+resolution=30000.;
 
 %To begin with the numerical model
 md=model();
-md=roundmesh(md,750000,resolution);
+md=roundmesh(md,750000.,resolution);
 md=setmask(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
 md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
@@ -13,7 +13,7 @@
 %Calculation of the analytical 2d velocity field
 constant=0.3;
-vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1;
-vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1;
-vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2));
+vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1;
+vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1;
+vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2);
 
 %We extrude the model to have a 3d model
@@ -23,7 +23,7 @@
 %Spc the nodes on the bed
 pos=find(md.mesh.vertexonbed);
-md.diagnostic.spcvx(pos)=0;
-md.diagnostic.spcvy(pos)=0;
-md.diagnostic.spcvz(pos)=0;
+md.diagnostic.spcvx(pos)=0.;
+md.diagnostic.spcvy(pos)=0.;
+md.diagnostic.spcvz(pos)=0.;
 
 %Now we can solve the problem 
@@ -36,13 +36,12 @@
 vel=zeros(md.mesh.numberofvertices2d,1);
 
-node_vel=0;
 for i=1:md.mesh.numberofvertices2d
-	for j=1:(md.mesh.numberoflayers-1)
-		node_vel=node_vel+1/(2*(md.mesh.numberoflayers-1))*(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+...
-			vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
+	node_vel=0.;
+	for j=1:md.mesh.numberoflayers-1
+		node_vel=node_vel+1./(2.*(md.mesh.numberoflayers-1))*...
+			(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
 			sqrt(vx(i+(j-1)*md.mesh.numberofvertices2d,1).^2+vy(i+(j-1)*md.mesh.numberofvertices2d,1).^2));
 	end
 	vel(i,1)=node_vel;
-	node_vel=0;
 end
 
@@ -82,5 +81,5 @@
 caxis([0 100]);
 
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('hutterstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
@@ -92,8 +91,8 @@
 	'Vx','Vy','Vel', ...
 };
-field_tolerances={...
+field_tolerances={ ...
 	1e-13,1e-13,1e-13, ...
 };
-field_values={
+field_values={ ...
 	vx,vy,vel, ...
 };
Index: /issm/trunk-jpl/test/NightlyRun/test1206.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1206.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1206.m	(revision 14134)
@@ -1,11 +1,11 @@
-%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
+%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling.
 printingflag=false;
 
 numlayers=10;
-resolution=30000;
+resolution=30000.;
 
 %To begin with the numerical model
 md=model();
-md=roundmesh(md,750000,resolution);
+md=roundmesh(md,750000.,resolution);
 md=setmask(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
 md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
@@ -13,7 +13,7 @@
 %Calculation of the analytical 2d velocity field
 constant=0.3;
-vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1;
-vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1;
-vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2));
+vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1;
+vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1;
+vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2);
 
 %We extrude the model to have a 3d model
@@ -23,7 +23,7 @@
 %Spc the nodes on the bed
 pos=find(md.mesh.vertexonbed);
-md.diagnostic.spcvx(pos)=0;
-md.diagnostic.spcvy(pos)=0;
-md.diagnostic.spcvz(pos)=0;
+md.diagnostic.spcvx(pos)=0.;
+md.diagnostic.spcvy(pos)=0.;
+md.diagnostic.spcvz(pos)=0.;
 
 %Now we can solve the problem 
@@ -36,13 +36,12 @@
 vel=zeros(md.mesh.numberofvertices2d,1);
 
-node_vel=0;
 for i=1:md.mesh.numberofvertices2d
-	for j=1:(md.mesh.numberoflayers-1)
-		node_vel=node_vel+1/(2*(md.mesh.numberoflayers-1))*(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+...
-			vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
+	node_vel=0.;
+	for j=1:md.mesh.numberoflayers-1
+		node_vel=node_vel+1./(2.*(md.mesh.numberoflayers-1))*...
+			(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
 			sqrt(vx(i+(j-1)*md.mesh.numberofvertices2d,1).^2+vy(i+(j-1)*md.mesh.numberofvertices2d,1).^2));
 	end
 	vel(i,1)=node_vel;
-	node_vel=0;
 end
 
@@ -81,5 +80,5 @@
 caxis([0 100]);
 
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('pattynstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
@@ -91,8 +90,8 @@
 	'Vx','Vy','Vel', ...
 };
-field_tolerances={...
+field_tolerances={ ...
 	1e-12,1e-12,1e-12, ...
 };
-field_values={
+field_values={ ...
 	vx,vy,vel, ...
 };
Index: /issm/trunk-jpl/test/NightlyRun/test1207.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1207.m	(revision 14133)
+++ /issm/trunk-jpl/test/NightlyRun/test1207.m	(revision 14134)
@@ -1,11 +1,11 @@
-%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
+%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling.
 printingflag=false;
 
 numlayers=10;
-resolution=30000;
+resolution=30000.;
 
 %To begin with the numerical model
 md=model();
-md=roundmesh(md,750000,resolution);
+md=roundmesh(md,750000.,resolution);
 md=setmask(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
 md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
@@ -13,7 +13,7 @@
 %Calculation of the analytical 2d velocity field
 constant=0.3;
-vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1;
-vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1;
-vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2));
+vx_obs=constant/2.*md.mesh.x.*(md.geometry.thickness).^-1;
+vy_obs=constant/2.*md.mesh.y.*(md.geometry.thickness).^-1;
+vel_obs=sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2);
 
 %We extrude the model to have a 3d model
@@ -23,7 +23,7 @@
 %Spc the nodes on the bed
 pos=find(md.mesh.vertexonbed);
-md.diagnostic.spcvx(pos)=0;
-md.diagnostic.spcvy(pos)=0;
-md.diagnostic.spcvz(pos)=0;
+md.diagnostic.spcvx(pos)=0.;
+md.diagnostic.spcvy(pos)=0.;
+md.diagnostic.spcvz(pos)=0.;
 
 %Now we can solve the problem 
@@ -36,13 +36,12 @@
 vel=zeros(md.mesh.numberofvertices2d,1);
 
-node_vel=0;
 for i=1:md.mesh.numberofvertices2d
-	for j=1:(md.mesh.numberoflayers-1)
-		node_vel=node_vel+1/(2*(md.mesh.numberoflayers-1))*(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+...
-			vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
+	node_vel=0.;
+	for j=1:md.mesh.numberoflayers-1
+		node_vel=node_vel+1./(2.*(md.mesh.numberoflayers-1))*...
+			(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
 			sqrt(vx(i+(j-1)*md.mesh.numberofvertices2d,1).^2+vy(i+(j-1)*md.mesh.numberofvertices2d,1).^2));
 	end
 	vel(i,1)=node_vel;
-	node_vel=0;
 end
 
@@ -81,5 +80,5 @@
 caxis([0 100]);
 
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('stokesstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
@@ -91,8 +90,8 @@
 	'Vx','Vy','Vel', ...
 };
-field_tolerances={...
+field_tolerances={ ...
 	1e-12,1e-12,1e-12, ...
 };
-field_values={
+field_values={ ...
 	vx,vy,vel, ...
 };
Index: /issm/trunk-jpl/test/Par/ISMIPA.par
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPA.par	(revision 14133)
+++ /issm/trunk-jpl/test/Par/ISMIPA.par	(revision 14134)
@@ -2,21 +2,21 @@
 
 disp('      creating thickness');
-md.geometry.surface=-md.mesh.x*tan(0.5*pi/180);
-md.geometry.bed=md.geometry.surface-1000+500*sin(md.mesh.x*2*pi/max(md.mesh.x)).*sin(md.mesh.y*2*pi/max(md.mesh.x));
+md.geometry.surface=-md.mesh.x*tan(0.5*pi/180.);
+md.geometry.bed=md.geometry.surface-1000.+500.*sin(md.mesh.x*2.*pi/max(md.mesh.x)).*sin(md.mesh.y*2.*pi/max(md.mesh.x));
 md.geometry.thickness=md.geometry.surface-md.geometry.bed;
 
 disp('      creating drag');
-md.friction.coefficient=200*ones(md.mesh.numberofvertices,1); %q=1.
+md.friction.coefficient=200.*ones(md.mesh.numberofvertices,1); %q=1.
 %Take care of iceshelves: no basal drag
 pos=find(md.mask.elementonfloatingice);
-md.friction.coefficient(md.mesh.elements(pos,:))=0;
+md.friction.coefficient(md.mesh.elements(pos,:))=0.;
 md.friction.p=ones(md.mesh.numberofelements,1);
 md.friction.q=ones(md.mesh.numberofelements,1);
 
-disp('      creating flow law paramter');
+disp('      creating flow law parameter');
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
 
 disp('      boundary conditions for diagnostic model');
-%Create node on boundary fist (because we cannot use mesh)
+%Create node on boundary first (because we cannot use mesh)
 md=SetIceSheetBC(md);
Index: /issm/trunk-jpl/test/Par/ISMIPA.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPA.py	(revision 14134)
+++ /issm/trunk-jpl/test/Par/ISMIPA.py	(revision 14134)
@@ -0,0 +1,25 @@
+import numpy
+from SetIceSheetBC import *
+
+#Ok, start defining model parameters here
+
+print "      creating thickness"
+md.geometry.surface=-md.mesh.x.reshape(-1,1)*numpy.tan(0.5*numpy.pi/180.)
+md.geometry.bed=md.geometry.surface-1000.+500.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))*numpy.sin(md.mesh.y.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))
+md.geometry.thickness=md.geometry.surface-md.geometry.bed
+
+print "      creating drag"
+md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1))    #q=1.
+#Take care of iceshelves: no basal drag
+pos=numpy.nonzero(md.mask.elementonfloatingice)[0]
+md.friction.coefficient[md.mesh.elements[pos,:]-1]=0.
+md.friction.p=numpy.ones((md.mesh.numberofelements,1))
+md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+
+print "      creating flow law parameter"
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+
+print "      boundary conditions for diagnostic model"
+#Create node on boundary first (because we cannot use mesh)
+md=SetIceSheetBC(md)
Index: /issm/trunk-jpl/test/Par/ISMIPB.par
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPB.par	(revision 14133)
+++ /issm/trunk-jpl/test/Par/ISMIPB.par	(revision 14134)
@@ -2,21 +2,21 @@
 
 disp('      creating thickness');
-md.geometry.surface=-md.mesh.x*tan(0.5*pi/180);
-md.geometry.bed=md.geometry.surface-1000+500*sin(md.mesh.x*2*pi/max(md.mesh.x));
+md.geometry.surface=-md.mesh.x*tan(0.5*pi/180.);
+md.geometry.bed=md.geometry.surface-1000.+500.*sin(md.mesh.x*2.*pi/max(md.mesh.x));
 md.geometry.thickness=md.geometry.surface-md.geometry.bed;
 
 disp('      creating drag');
-md.friction.coefficient=200*ones(md.mesh.numberofvertices,1); %q=1.
+md.friction.coefficient=200.*ones(md.mesh.numberofvertices,1); %q=1.
 %Take care of iceshelves: no basal drag
 pos=find(md.mask.elementonfloatingice);
-md.friction.coefficient(md.mesh.elements(pos,:))=0;
+md.friction.coefficient(md.mesh.elements(pos,:))=0.;
 md.friction.p=ones(md.mesh.numberofelements,1);
 md.friction.q=ones(md.mesh.numberofelements,1);
 
-disp('      creating flow law paramter');
+disp('      creating flow law parameter');
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
 
 disp('      boundary conditions for diagnostic model');
-%Create node on boundary fist (because we cannot use mesh)
+%Create node on boundary first (because we cannot use mesh)
 md=SetIceSheetBC(md);
Index: /issm/trunk-jpl/test/Par/ISMIPB.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPB.py	(revision 14134)
+++ /issm/trunk-jpl/test/Par/ISMIPB.py	(revision 14134)
@@ -0,0 +1,25 @@
+import numpy
+from SetIceSheetBC import *
+
+#Ok, start defining model parameters here
+
+print "      creating thickness"
+md.geometry.surface=-md.mesh.x.reshape(-1,1)*numpy.tan(0.5*numpy.pi/180.)
+md.geometry.bed=md.geometry.surface-1000.+500.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))
+md.geometry.thickness=md.geometry.surface-md.geometry.bed
+
+print "      creating drag"
+md.friction.coefficient=200.*numpy.ones((md.mesh.numberofvertices,1))    #q=1.
+#Take care of iceshelves: no basal drag
+pos=numpy.nonzero(md.mask.elementonfloatingice)[0]
+md.friction.coefficient[md.mesh.elements[pos,:]-1]=0.
+md.friction.p=numpy.ones((md.mesh.numberofelements,1))
+md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+
+print "      creating flow law parameter"
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+
+print "      boundary conditions for diagnostic model"
+#Create node on boundary first (because we cannot use mesh)
+md=SetIceSheetBC(md)
Index: /issm/trunk-jpl/test/Par/ISMIPC.par
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPC.par	(revision 14133)
+++ /issm/trunk-jpl/test/Par/ISMIPC.par	(revision 14134)
@@ -2,22 +2,22 @@
 
 disp('      creating thickness');
-md.geometry.surface=2000-md.mesh.x*tan(0.1*pi/180); %to have z>0
-md.geometry.bed=md.geometry.surface-1000;
+md.geometry.surface=2000.-md.mesh.x*tan(0.1*pi/180.); %to have z>0
+md.geometry.bed=md.geometry.surface-1000.;
 md.geometry.thickness=md.geometry.surface-md.geometry.bed;
 
 disp('      creating drag');
-%md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x/2)).*sin(md.mesh.y*2*pi/max(md.mesh.x/2)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)));
-md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x)).*sin(md.mesh.y*2*pi/max(md.mesh.x))));
+%md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x/2.)).*sin(md.mesh.y*2.*pi/max(md.mesh.x/2.)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)));
+md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x)).*sin(md.mesh.y*2.*pi/max(md.mesh.x))));
 %Take care of iceshelves: no basal drag
 pos=find(md.mask.elementonfloatingice);
-md.friction.coefficient(md.mesh.elements(pos,:))=0;
+md.friction.coefficient(md.mesh.elements(pos,:))=0.;
 md.friction.p=ones(md.mesh.numberofelements,1);
 md.friction.q=zeros(md.mesh.numberofelements,1);
 
-disp('      creating flow law paramter');
+disp('      creating flow law parameter');
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
 
-disp('      boundary conditions for diagnostic model: ');
-%Create node on boundary fist (because wi can not use mesh)
+disp('      boundary conditions for diagnostic model:');
+%Create node on boundary first (because we can not use mesh)
 md=SetIceSheetBC(md);
Index: /issm/trunk-jpl/test/Par/ISMIPC.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPC.py	(revision 14134)
+++ /issm/trunk-jpl/test/Par/ISMIPC.py	(revision 14134)
@@ -0,0 +1,26 @@
+import numpy
+from SetIceSheetBC import *
+
+#Ok, start defining model parameters here
+
+print "      creating thickness"
+md.geometry.surface=2000.-md.mesh.x.reshape(-1,1)*numpy.tan(0.1*numpy.pi/180.)    #to have z>0
+md.geometry.bed=md.geometry.surface-1000.
+md.geometry.thickness=md.geometry.surface-md.geometry.bed
+
+print "      creating drag"
+#md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x/2.)).*sin(md.mesh.y*2.*pi/max(md.mesh.x/2.)))./(md.constants.g*(md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.bed)));
+md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))*numpy.sin(md.mesh.y.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))))
+#Take care of iceshelves: no basal drag
+pos=numpy.nonzero(md.mask.elementonfloatingice)[0]
+md.friction.coefficient[md.mesh.elements[pos,:]-1]=0.
+md.friction.p=numpy.ones((md.mesh.numberofelements,1))
+md.friction.q=numpy.zeros((md.mesh.numberofelements,1))
+
+print "      creating flow law parameter"
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+
+print "      boundary conditions for diagnostic model:"
+#Create node on boundary first (because we can not use mesh)
+md=SetIceSheetBC(md)
Index: /issm/trunk-jpl/test/Par/ISMIPD.par
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPD.par	(revision 14133)
+++ /issm/trunk-jpl/test/Par/ISMIPD.par	(revision 14134)
@@ -2,21 +2,21 @@
 
 disp('      creating thickness');
-md.geometry.surface=2000-md.mesh.x*tan(0.1*pi/180); %to have z>0
-md.geometry.bed=md.geometry.surface-1000;
+md.geometry.surface=2000.-md.mesh.x*tan(0.1*pi/180.); %to have z>0
+md.geometry.bed=md.geometry.surface-1000.;
 md.geometry.thickness=md.geometry.surface-md.geometry.bed;
 
 disp('      creating drag');
-md.friction.coefficient=sqrt(md.constants.yts.*(1000+1000*sin(md.mesh.x*2*pi/max(md.mesh.x))));
+md.friction.coefficient=sqrt(md.constants.yts.*(1000.+1000.*sin(md.mesh.x*2.*pi/max(md.mesh.x))));
 %Take care of iceshelves: no basal drag
 pos=find(md.mask.elementonfloatingice);
-md.friction.coefficient(md.mesh.elements(pos,:))=0;
+md.friction.coefficient(md.mesh.elements(pos,:))=0.;
 md.friction.p=ones(md.mesh.numberofelements,1);
 md.friction.q=zeros(md.mesh.numberofelements,1);
 
-disp('      creating flow law paramter');
+disp('      creating flow law parameter');
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
 
-disp('      boundary conditions for diagnostic model: ');
-%Create node on boundary fist (because wi can not use mesh)
+disp('      boundary conditions for diagnostic model:');
+%Create node on boundary first (because we can not use mesh)
 md=SetIceSheetBC(md);
Index: /issm/trunk-jpl/test/Par/ISMIPD.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPD.py	(revision 14134)
+++ /issm/trunk-jpl/test/Par/ISMIPD.py	(revision 14134)
@@ -0,0 +1,25 @@
+import numpy
+from SetIceSheetBC import *
+
+#Ok, start defining model parameters here
+
+print "      creating thickness"
+md.geometry.surface=2000.-md.mesh.x.reshape(-1,1)*numpy.tan(0.1*numpy.pi/180.)    #to have z>0
+md.geometry.bed=md.geometry.surface-1000.
+md.geometry.thickness=md.geometry.surface-md.geometry.bed
+
+print "      creating drag"
+md.friction.coefficient=numpy.sqrt(md.constants.yts*(1000.+1000.*numpy.sin(md.mesh.x.reshape(-1,1)*2.*numpy.pi/numpy.max(md.mesh.x))))
+#Take care of iceshelves: no basal drag
+pos=numpy.nonzero(md.mask.elementonfloatingice)[0]
+md.friction.coefficient[md.mesh.elements[pos,:]-1]=0.
+md.friction.p=numpy.ones((md.mesh.numberofelements,1))
+md.friction.q=numpy.zeros((md.mesh.numberofelements,1))
+
+print "      creating flow law parameter"
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+
+print "      boundary conditions for diagnostic model:"
+#Create node on boundary first (because we can not use mesh)
+md=SetIceSheetBC(md)
Index: /issm/trunk-jpl/test/Par/ISMIPE.par
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPE.par	(revision 14133)
+++ /issm/trunk-jpl/test/Par/ISMIPE.par	(revision 14134)
@@ -7,9 +7,9 @@
 for i=1:md.mesh.numberofvertices
 	y=md.mesh.y(i);
-	point1=floor(y/100)+1;
-	point2=min(point1+1,51);
-	coeff=(y-(point1-1)*100)/100;
-	md.geometry.bed(i)=(1-coeff)*data(point1,2)+coeff*data(point2,2);
-	md.geometry.surface(i)=(1-coeff)*data(point1,3)+coeff*data(point2,3);
+	point1=floor(y/100.)+1.;
+	point2=min(point1+1,51.);
+	coeff=(y-(point1-1.)*100.)/100.;
+	md.geometry.bed(i)=(1.-coeff)*data(point1,2)+coeff*data(point2,2);
+	md.geometry.surface(i)=(1.-coeff)*data(point1,3)+coeff*data(point2,3);
 end
 md.geometry.thickness=md.geometry.surface-md.geometry.bed;
@@ -22,9 +22,9 @@
 md.friction.q=ones(md.mesh.numberofelements,1);
 
-disp('      creating flow law paramter');
+disp('      creating flow law parameter');
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
 
-disp('      boundary conditions for diagnostic model: ');
-%Create node on boundary fist (because wi can not use mesh)
+disp('      boundary conditions for diagnostic model:');
+%Create node on boundary first (because we can not use mesh)
 md=SetIceSheetBC(md);
Index: /issm/trunk-jpl/test/Par/ISMIPE.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPE.py	(revision 14134)
+++ /issm/trunk-jpl/test/Par/ISMIPE.py	(revision 14134)
@@ -0,0 +1,35 @@
+import numpy
+import netCDF4
+from SetIceSheetBC import *
+
+#Ok, start defining model parameters here
+
+print "      creating thickness"
+f = netCDF4.Dataset('../Data/ISMIPE.nc','r')
+data = f.variables['data'][:]
+f.close()
+md.geometry.surface=numpy.zeros((md.mesh.numberofvertices,1))
+md.geometry.bed=numpy.zeros((md.mesh.numberofvertices,1))
+for i in xrange(0,md.mesh.numberofvertices):
+	y=md.mesh.y[i]
+	point1=numpy.floor(y/100.)
+	point2=numpy.minimum(point1+1,50)
+	coeff=(y-(point1-1.)*100.)/100.
+	md.geometry.bed[i]=(1.-coeff)*data[point1,1]+coeff*data[point2,1]
+	md.geometry.surface[i]=(1.-coeff)*data[point1,2]+coeff*data[point2,2]
+md.geometry.thickness=md.geometry.surface-md.geometry.bed
+md.geometry.thickness[numpy.nonzero(numpy.logical_not(md.geometry.thickness))]=0.01
+md.geometry.bed=md.geometry.surface-md.geometry.thickness
+
+print "      creating drag"
+md.friction.coefficient=numpy.zeros((md.mesh.numberofvertices,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements,1))
+md.friction.q=numpy.ones((md.mesh.numberofelements,1))
+
+print "      creating flow law parameter"
+md.materials.rheology_B=6.8067*10**7*numpy.ones((md.mesh.numberofvertices,1))
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+
+print "      boundary conditions for diagnostic model:"
+#Create node on boundary first (because we can not use mesh)
+md=SetIceSheetBC(md)
Index: /issm/trunk-jpl/test/Par/ISMIPF.par
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPF.par	(revision 14133)
+++ /issm/trunk-jpl/test/Par/ISMIPF.par	(revision 14134)
@@ -3,23 +3,23 @@
 
 disp('      creating thickness');
-md.geometry.surface=-md.mesh.x*tan(3*pi/180);
-%md.geometry.bed=md.geometry.surface-1000;
-md.geometry.bed=md.geometry.surface-1000+100*exp(-((md.mesh.x-max(md.mesh.x)/2).^2+(md.mesh.y-max(md.mesh.y)/2).^2)/(10000^2));
+md.geometry.surface=-md.mesh.x*tan(3.*pi/180.);
+%md.geometry.bed=md.geometry.surface-1000.;
+md.geometry.bed=md.geometry.surface-1000.+100.*exp(-((md.mesh.x-max(md.mesh.x)/2.).^2+(md.mesh.y-max(md.mesh.y)/2.).^2)/(10000.^2));
 md.geometry.thickness=md.geometry.surface-md.geometry.bed;
 
 disp('      creating drag');
-md.friction.coefficient=sqrt(md.constants.yts/(2.140373*10^-7*1000))*ones(md.mesh.numberofvertices,1);
+md.friction.coefficient=sqrt(md.constants.yts/(2.140373*10^-7*1000.))*ones(md.mesh.numberofvertices,1);
 md.friction.p=ones(md.mesh.numberofelements,1);
 md.friction.q=zeros(md.mesh.numberofelements,1);
 
-disp('      creating flow law paramter');
+disp('      creating flow law parameter');
 md.materials.rheology_B=1.4734*10^14*ones(md.mesh.numberofvertices,1);
-md.materials.rheology_n=1*ones(md.mesh.numberofelements,1);
+md.materials.rheology_n=1.*ones(md.mesh.numberofelements,1);
 md.materials.rheology_law='None';
 
 disp('      boundary conditions for diagnostic model');
-%Create node on boundary fist (because we cannot use mesh)
+%Create node on boundary first (because we cannot use mesh)
 md=SetIceSheetBC(md);
-md.diagnostic.spcvx=100*ones(md.mesh.numberofvertices,1);
+md.diagnostic.spcvx=100.*ones(md.mesh.numberofvertices,1);
 md.initialization.vx=zeros(md.mesh.numberofvertices,1);
 md.initialization.vy=zeros(md.mesh.numberofvertices,1);
@@ -27,5 +27,5 @@
 md.initialization.vel=zeros(md.mesh.numberofvertices,1);
 md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
-md.initialization.temperature=255*ones(md.mesh.numberofvertices,1);
+md.initialization.temperature=255.*ones(md.mesh.numberofvertices,1);
 pos=find(md.mesh.x==min(md.mesh.x) | md.mesh.x==max(md.mesh.x) | md.mesh.y==min(md.mesh.y) | md.mesh.y==max(md.mesh.y));
 md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
@@ -33,5 +33,5 @@
 md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
 md.prognostic.spcthickness(pos)=md.geometry.thickness(pos);
-md.thermal.spctemperature=255*ones(md.mesh.numberofvertices,1);
+md.thermal.spctemperature=255.*ones(md.mesh.numberofvertices,1);
 md.basalforcings.geothermalflux=0.4*ones(md.mesh.numberofvertices,1);
 
@@ -40,6 +40,6 @@
 
 %Transient options
-md.timestepping.time_step=1;
-md.timestepping.final_time=10;
+md.timestepping.time_step=1.;
+md.timestepping.final_time=10.;
 md.prognostic.stabilization=1;
 md.thermal.stabilization=1;
Index: /issm/trunk-jpl/test/Par/ISMIPF.py
===================================================================
--- /issm/trunk-jpl/test/Par/ISMIPF.py	(revision 14134)
+++ /issm/trunk-jpl/test/Par/ISMIPF.py	(revision 14134)
@@ -0,0 +1,49 @@
+import numpy
+from SetIceSheetBC import *
+
+#Ok, start defining model parameters here
+md.verbose=2
+
+print "      creating thickness"
+md.geometry.surface=-md.mesh.x.reshape(-1,1)*numpy.tan(3.*numpy.pi/180.)
+#md.geometry.bed=md.geometry.surface-1000.
+md.geometry.bed=md.geometry.surface-1000.+100.*numpy.exp(-((md.mesh.x.reshape(-1,1)-numpy.max(md.mesh.x)/2.)**2+(md.mesh.y.reshape(-1,1)-numpy.max(md.mesh.y)/2.)**2)/(10000.**2))
+md.geometry.thickness=md.geometry.surface-md.geometry.bed
+
+print "      creating drag"
+md.friction.coefficient=numpy.sqrt(md.constants.yts/(2.140373*10**-7*1000.))*numpy.ones((md.mesh.numberofvertices,1))
+md.friction.p=numpy.ones((md.mesh.numberofelements,1))
+md.friction.q=numpy.zeros((md.mesh.numberofelements,1))
+
+print "      creating flow law parameter"
+md.materials.rheology_B=1.4734*10**14*numpy.ones((md.mesh.numberofvertices,1))
+md.materials.rheology_n=1.*numpy.ones((md.mesh.numberofelements,1))
+md.materials.rheology_law='None'
+
+print "      boundary conditions for diagnostic model"
+#Create node on boundary first (because we cannot use mesh)
+md=SetIceSheetBC(md)
+md.diagnostic.spcvx=100.*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.vx=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vy=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vz=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.vel=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1))
+md.initialization.temperature=255.*numpy.ones((md.mesh.numberofvertices,1))
+pos=numpy.nonzero(numpy.logical_or(numpy.logical_or(md.mesh.x==numpy.min(md.mesh.x),md.mesh.x==numpy.max(md.mesh.x)),numpy.logical_or(md.mesh.y==numpy.min(md.mesh.y),md.mesh.y==numpy.max(md.mesh.y))))
+md.balancethickness.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+md.balancethickness.spcthickness[pos]=md.geometry.thickness[pos]
+md.prognostic.spcthickness=float('NaN')*numpy.ones((md.mesh.numberofvertices,1))
+md.prognostic.spcthickness[pos]=md.geometry.thickness[pos]
+md.thermal.spctemperature=255.*numpy.ones((md.mesh.numberofvertices,1))
+md.basalforcings.geothermalflux=0.4*numpy.ones((md.mesh.numberofvertices,1))
+
+#Parallel options
+md.mesh.average_vertex_connectivity=200
+
+#Transient options
+md.timestepping.time_step=1.
+md.timestepping.final_time=10.
+md.prognostic.stabilization=1
+md.thermal.stabilization=1
+md.thermal.penalty_threshold=10**5