Index: /issm/trunk-jpl/src/m/solve/WriteData.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/WriteData.py	(revision 22271)
+++ /issm/trunk-jpl/src/m/solve/WriteData.py	(revision 22272)
@@ -45,5 +45,5 @@
 		scale = options.getfieldvalue('scale')
 		if np.size(data) > 1 and np.ndim(data) > 1 and np.size(data,0)==timeserieslength:
-			data[0:-1,:] = scale*data[0:-1,:]
+			data[0:-2,:] = scale*data[0:-2,:]
 		else:
 			data  = scale*data
Index: /issm/trunk-jpl/test/NightlyRun/test1108.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1108.py	(revision 22271)
+++ /issm/trunk-jpl/test/NightlyRun/test1108.py	(revision 22272)
@@ -8,4 +8,5 @@
 from setflowequation import *
 from solve import *
+from squaremesh import *
 
 from PythonFuncs import *
@@ -16,4 +17,5 @@
 """
 
+#L_list=[5000.,10000.,20000.,40000.,80000.,160000.]
 L_list=[80000.]
 results=[]
@@ -35,5 +37,11 @@
 	md.stressbalance.spcvz=np.nan*np.ones((md.mesh.numberofvertices))
 	
-	pos=np.nonzero(logical_and.reduce_n(md.mesh.vertexonbase,np.logical_or.reduce(md.mesh.x==0.,md.mesh.x==np.max(md.mesh.x)),np.logical_or.reduce(md.mesh.y==0.,md.mesh.y==np.max(md.mesh.y))))
+	maxx = max(md.mesh.x)
+	maxy = max(md.mesh.y)
+	posA = np.where(md.mesh.vertexonbase)
+	posB = np.unique(np.concatenate((np.where(md.mesh.x==0.), np.where(md.mesh.x==maxx))))
+	posC = np.unique(np.concatenate((np.where(md.mesh.y==0.), np.where(md.mesh.y==maxy))))
+	pos = np.intersect1d(np.intersect1d(posA,posB),posC)
+
 	md.stressbalance.spcvx[pos]=0.
 	md.stressbalance.spcvy[pos]=0.
@@ -44,6 +52,6 @@
 	posx2=np.nonzero(md.mesh.x==np.max(md.mesh.x))[0]
 
-	posy=np.nonzero(logical_and.reduce_n(md.mesh.y==0.,md.mesh.x!=0.,md.mesh.x!=np.max(md.mesh.x)))[0]    #Don't take the same nodes two times
-	posy2=np.nonzero(logical_and.reduce_n(md.mesh.y==np.max(md.mesh.y),md.mesh.x!=0.,md.mesh.x!=np.max(md.mesh.x)))[0]
+	posy=np.intersect1d(np.intersect1d(np.where(md.mesh.y==0.),np.where(md.mesh.x!=0.)),np.where(md.mesh.x!=np.max(md.mesh.x)))[0]    #Don't take the same nodes two times
+	posy2=np.intersect1d(np.intersect1d(np.where(md.mesh.y==np.max(md.mesh.y)),np.where(md.mesh.x!=0.)),np.where(md.mesh.x!=np.max(md.mesh.x)))[0]
 
 	md.stressbalance.vertex_pairing=np.vstack((np.hstack((posx.reshape(-1,1)+1,posx2.reshape(-1,1)+1)),np.hstack((posy.reshape(-1,1)+1,posy2.reshape(-1,1)+1))))
Index: /issm/trunk-jpl/test/NightlyRun/test1109.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1109.py	(revision 22272)
+++ /issm/trunk-jpl/test/NightlyRun/test1109.py	(revision 22272)
@@ -0,0 +1,86 @@
+#Test Name: ISMIPE
+import numpy as np
+import sys
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from squaremesh import *
+
+#This test is a test from the ISMP-HOM Intercomparison project.
+#TestE 
+#Four tests to run: - Pattyn frozen
+#                   - Stokes frozen
+#                   - Pattyn with some sliding
+#                   - Stokes with some sliding
+printingflag = False
+results = []
+
+for i in range(4):
+	Lx=10. #in m
+	Ly=5000. #in m
+	nx=3 #number of nodes in x direction
+	ny=51
+	md = model()
+	md = squaremesh(md,Lx,Ly,nx,ny)
+	md = setmask(md,'','') #ice sheet test
+	md = parameterize(md,'../Par/ISMIPE.py')
+	md = md.extrude(10,1.)
+
+	if i==0 or i==2:
+		md = setflowequation(md,'HO','all')
+	elif i==1 or i==3:
+		md = setflowequation(md,'FS','all')
+
+	#Create MPCs to have periodic boundary conditions
+	posx = np.where(md.mesh.x == 0.)[0]
+	posx2 = np.where(md.mesh.x == max(md.mesh.x))[0]
+	md.stressbalance.vertex_pairing = np.column_stack((posx,posx2))
+
+	#Create spcs on the bed 
+	pos = np.where(md.mesh.vertexonbase)[0]
+	md.stressbalance.spcvx = float('NaN') * np.ones((md.mesh.numberofvertices,))
+	md.stressbalance.spcvy = float('NaN') * np.ones((md.mesh.numberofvertices,))
+	md.stressbalance.spcvz = float('NaN') * np.ones((md.mesh.numberofvertices,))
+	md.stressbalance.spcvx[pos] = 0.
+	md.stressbalance.spcvy[pos] = 0.
+	md.stressbalance.spcvz[pos] = 0.
+
+	#Remove the spc where there is some sliding (case 3 and 4):
+	if i==2 or i==3:
+		pos = np.intersect1d(np.where((md.mesh.y / max(md.mesh.y)) >= 0.44), np.where((md.mesh.y / max(md.mesh.y)) <= 0.5))[0]
+		md.stressbalance.spcvx[pos] = float('NaN')
+		md.stressbalance.spcvy[pos] = float('NaN')
+		md.stressbalance.spcvz[pos] = float('NaN')
+
+	#Compute the stressbalance
+	md.cluster = generic('name',gethostname(),'np',8)
+	md = solve(md,'Stressbalance')
+
+	vx = md.results.StressbalanceSolution.Vx
+	vy = md.results.StressbalanceSolution.Vy
+	vz = md.results.StressbalanceSolution.Vz
+	results[i] = md.results.StressbalanceSolution
+
+
+#Fields and tolerances to track changes
+field_names = [
+	'VyPattynSliding','VzPattynSliding',
+	'VxStokesSliding','VyStokesSliding','VzStokesSliding',
+	'VyPattynFrozen','VzPattynFrozen',
+	'VxStokesFrozen','VyStokesFrozen','VzStokesFrozen'
+	]
+field_tolerances = [
+	1e-05,1e-05,
+	1e-05,1e-06,1e-06,
+	1e-05,1e-04,
+	1e-05,1e-05,1e-06
+	]
+
+field_values = []
+for i in range(4):
+	result = results[i]
+	field_values += [result.Vx,result.Vy,result.Vz]
Index: /issm/trunk-jpl/test/NightlyRun/test1110.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test1110.py	(revision 22272)
+++ /issm/trunk-jpl/test/NightlyRun/test1110.py	(revision 22272)
@@ -0,0 +1,141 @@
+#Test Name: ISMIPF
+import numpy as np
+from model import *
+from socket import gethostname
+from bamg import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from squaremesh import *
+
+#This test is a test from the ISMP-HOM Intercomparison project.
+#TestF 
+printingflag = False
+results = []
+
+for i in range(4):
+	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 = setmask(md,'','') #ice sheet test
+	md = parameterize(md,'../Par/ISMIPF.py')
+	md = md.extrude(4,1.)
+
+	if (i == 0 or i == 1):
+		md = setflowequation(md,'HO','all')
+	else:
+		md = setflowequation(md,'FS','all')
+
+	md.stressbalance.spcvx = float('NaN') * np.ones((md.mesh.numberofvertices,))
+	md.stressbalance.spcvy = float('NaN') * np.ones((md.mesh.numberofvertices,))
+	md.stressbalance.spcvz = float('NaN') * np.ones((md.mesh.numberofvertices,))
+	if (i == 0 or i == 2):
+		#Create dirichlet on the bed if no slip
+		pos = np.where(md.mesh.vertexonbase)
+		md.stressbalance.spcvx[pos] = 0.
+		md.stressbalance.spcvy[pos] = 0.
+		md.stressbalance.spcvz[pos] = 0.
+	else:
+		posA = np.where(md.mesh.vertexonbase)
+		posB = np.unique(np.concatenate(np.where(md.mesh.x == 0.),np.where(md.mesh.x == max(md.mesh.x))))
+		posC = np.unique(np.concatenate(np.where(md.mesh.y == 0.),np.where(md.mesh.y == max(md.mesh.y))))
+		pos = np.intersect1d(np.intersect1d(posA,posB),posC)
+		md.stressbalance.spcvx[pos] = 100. #because we need a dirichlet somewhere
+		md.stressbalance.spcvy[pos] = 0.
+		md.stressbalance.spcvz[pos] = 0.
+	
+	pos = np.where(np.logical_not(md.mesh.vertexonbase))
+	md.thermal.spctemperature[pos] = 255.
+
+	#Create MPCs to have periodic boundary conditions
+	posx = np.where(md.mesh.x == 0.)
+	posx2 = np.where(md.mesh.x == max(md.mesh.x))
+
+	posy = np.where(md.mesh.y == 0.)
+	posy2 = np.where(md.mesh.y == max(md.mesh.y))
+
+	md.stressbalance.vertex_pairing = np.column_stack((posx,posx2,posy,posy2))
+	md.masstransport.vertex_pairing = np.column_stack((posx,posx2,posy,posy2))
+
+	md.timestepping.time_step = 3.
+	md.timestepping.final_time = 300.
+	md.settings.output_frequency = 50
+	md.masstransport.stabilization = 1
+	md.stressbalance.maxiter = 1
+	
+	#Compute the stressbalance
+	md.cluster = generic('name',gethostname(),'np',8)
+	md.verbose = verbose('convergence',True,'solution',True)
+	md = solve(md,'Transient')
+
+	#save the results
+	results[i] = md.results.TransientSolution()
+	
+	#Now plot vx and delta surface 
+	if (i == 0 or i == 2):
+		plotmodel(md,'data',(md.results.TransientSolution().Vx),'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIP100000.exp','title','','xlabel','','ylabel','Velocity (m/yr)','linewidth',3,'grid','on','unit','km','ylim',[91,100])
+	elif (i == 1 or i == 3):
+		plotmodel(md,'data',(md.results.TransientSolution().Vx),'layer',md.mesh.numberoflayers,'sectionvalue','../Exp/ISMIP100000.exp','title','','xlabel','','ylabel','Velocity (m/yr)','linewidth',3,'grid','on','unit','km','ylim',[185,200])
+	
+	if printingflag:
+		#set(gcf,'Color','w')
+		if i == 0:
+			printmodel('ismipfHOvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
+			#system(['mv ismipfHOvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF'])
+		elif i == 1:
+			printmodel('ismipfHOvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
+			#system(['mv ismipfHOvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF'])
+		elif i == 2:
+			printmodel('ismipfFSvxfrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
+			#system(['mv ismipfFSvxfrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF'])
+		elif i == 3:
+			printmodel('ismipfFSvxsliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
+			#system(['mv ismipfFSvxsliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF'])
+		
+	
+
+	plotmodel(md,'data',(md.results.TransientSolution().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:
+		#set(gcf,'Color','w')
+		if i == 0:
+			printmodel('ismipfHOdeltasurfacefrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
+			#system(['mv ismipfHOdeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF'])
+		elif i == 1:
+			printmodel('ismipfHOdeltasurfacesliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
+			#system(['mv ismipfHOdeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF'])
+		elif i == 2:
+			printmodel('ismipfFSdeltasurfacefrozen','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
+			#system(['mv ismipfFSdeltasurfacefrozen.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF'])
+		elif i == 3:
+			printmodel('ismipfFSdeltasurfacesliding','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
+			#system(['mv ismipfFSdeltasurfacesliding.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestF'])
+		
+	
+
+
+#Fields and tolerances to track changes
+field_names = [
+	'VxPattynFrozen','VyPattynFrozen','VzPattynFrozen','SurfacePattynFrozen',
+	'VxPattynSliding','VyPattynSliding','VzPattynSliding','SurfacePattynSliding',
+	'VxStokesFrozen','VyStokesFrozen','VzStokesFrozen','SurfaceStokesFrozen',
+	'VxStokesSliding','VyStokesSliding','VzStokesSliding','SurfaceStokesSliding'
+]
+field_tolerances = [
+	1e-10,1e-09,1e-09,1e-10,
+	1e-10,1e-09,1e-09,1e-10,
+	1e-08,1e-09,1e-08,1e-09,
+	1e-08,2e-09,1e-08,1e-09
+]
+field_values = []
+for i in range(4):
+	result = results[i]
+	field_values += ([
+		result.Vx,
+		result.Vy,
+		result.Vz,
+		result.Surface] - md.geometry.surface)
+