Index: issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1301.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1301.py	(revision 14101)
+++ 	(revision )
@@ -1,87 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test1301.m
-Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-from numpy import *
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-#  This file can be run to check that the melting in simple conduction is correctly modeled.
-
-#  There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
-
-#  at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower
-
-#  surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed...
-
-printingflag=false
-
-
-md=model()
-md=triangle(md,'../Exp/Square.exp',100000)
-md=setmask(md,'','')
-md=parameterize(md,'../Par/SquareThermal.py')
-md.extrude(3,2.)
-md=setflowequation(md,'Pattyn','all')
-
-
-# Some conditions specific to melting test
-
-md.initialization.pressure=zeros(md.mesh.numberofvertices,1)
-md.initialization.temperature=273.15*ones(md.mesh.numberofvertices,1)
-pos=numpy.nonzero(md.mesh.vertexonsurface)
-md.thermal.spctemperature(pos)=md.initialization.temperature(pos)
-md.materials.rheology_B=paterson(md.initialization.temperature)
-
-
-# analytical results
-
-# melting heat = geothermal flux
-
-# Mb*L*rho=G   => Mb=G/L*rho
-
-melting=md.basalforcings.geothermalflux/(md.materials.rho_ice*md.materials.latentheat)*md.constants.yts
-
-
-# modeled  results
-
-md.cluster=generic('name',oshostname(),'np',2)
-md=solve(md,ThermalSolutionEnum())
-
-
-# plot results
-
-comp_melting=md.results.ThermalSolution.BasalforcingsMeltingRate
-relative=abs((comp_melting-melting)./melting)*100
-relative(find(comp_melting==melting))=0
-plotmodel()(md,'data',comp_melting,'title','Modeled melting','data',melting,'title','Analytical melting',\
-	'data',comp_melting-melting,'title','Absolute error','data',relative,'title','Relative error [%]',\
-	'layer#all',1,'caxis#2',[1.02964 1.02966]*10^-4,'FontSize#all',20,'figposition','mathieu')
-if printingflag, 
-	set(gcf,'Color','w')
-	printmodel()('thermalmelting','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off')
-	system(['mv thermalmelting.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal '])
-end
-
-
-
-
-# Fields and tolerances to track changes
-
-field_names     =['BasalMelting']
-field_tolerances=[1e-08]
-field_values    =[comp_melting]
Index: issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1302.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1302.py	(revision 14101)
+++ 	(revision )
@@ -1,86 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test1302.m
-Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-from numpy import *
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-#  This file can be run to check that the advection-diffusion  is correctly modeled.
-
-#  There is u=v=0 and w=cst everywhere the only thermal boundary conditions are an imposed temperature
-
-#  at upper surface and an impose flux at its base.
-
-printingflag=false
-
-
-md=model()
-md=triangle(md,'../Exp/Square.exp',100000)
-md=setmask(md,'','')
-md=parameterize(md,'../Par/SquareThermal.py')
-md.extrude(30,1.)   %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%.)
-md=setflowequation(md,'Pattyn','all') 
-
-
-# Thermal boundary conditions
-
-pos1=numpy.nonzero(md.mesh.elementonbed)     md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10
-pos2=numpy.nonzero(md.mesh.elementonsurface) md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0
-md.initialization.vz=0.1*ones(md.mesh.numberofvertices,1)
-md.initialization.vel=sqrt( md.initialization.vx.^2+ md.initialization.vy.^2+ md.initialization.vz.^2)
-md.initialization.pressure=zeros(md.mesh.numberofvertices,1)
-
-
-md.thermal.stabilization=2
-# analytical results
-
-# d2T/dz2-w*rho_ice*c/k*dT/dz=0   T(surface)=0  T(bed)=10   => T=A exp(alpha z)+B
-
-alpha=0.1/md.constants.yts*md.materials.rho_ice*md.materials.heatcapacity/md.materials.thermalconductivity   %alpha=w rho_ice c /k  and w=0.1m/an
-A=10/(exp(alpha*(-1000))-1)    %A=T(bed)/(exp(alpha*bed)-1)  with bed=-1000 T(bed)=10
-B=-A
-md.initialization.temperature=A*exp(alpha*md.mesh.z)+B
-
-
-# modeled  results
-
-md.cluster=generic('name',oshostname(),'np',2)
-md=solve(md,ThermalSolutionEnum())
-
-
-# plot results
-
-comp_temp=md.results.ThermalSolution.Temperature
-relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100
-relative(find(comp_temp==md.initialization.temperature))=0
-plotmodel()(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,\
-	'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,\
-	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,\
-	'figposition','mathieu','FontSize#all',20)
-if printingflag, 
-	set(gcf,'Color','w')
-	printmodel()('thermaladvection','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off')
-	system(['mv thermaladvection.png ' ISSM_DIR '/website/doc_pdf/validation/Images/EISMINT '])
-end
-
-
-# Fields and tolerances to track changes
-
-field_names     =['AdvectionTemperature']
-field_tolerances=[1e-13]
-field_values    =[comp_temp]
Index: issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1303.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1303.py	(revision 14101)
+++ 	(revision )
@@ -1,80 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test1303.m
-Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-from numpy import *
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-#  This file can be run to check that the conduction is correctly modeled.
-
-#  There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
-
-#  at the lower and upper surface. The result must be a linear temperature from the upper to the lower
-
-#  surface. if it is not the case, something is thermal modeling has been changed...
-
-printingflag=false
-
-
-md=model()
-md=triangle(md,'../Exp/Square.exp',100000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareThermal.py')
-md.extrude(11,2.)
-md=setflowequation(md,'Pattyn','all')
-pos1=numpy.nonzero(md.mesh.elementonbed)     md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10
-pos2=numpy.nonzero(md.mesh.elementonsurface) md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0
-md.initialization.pressure=zeros(md.mesh.numberofvertices,1)
-
-
-# analytical results
-
-# d2T/dz2=0 T(bed)=10 T(surface)=0  => T=0*(z-bed)/thickness+10*(surface-z)/thickness
-
-# each layer of the 3d mesh must have a constant value
-
-md.initialization.temperature=10*(md.geometry.surface-md.mesh.z)./md.geometry.thickness
-
-
-# modeled  results
-
-md.cluster=generic('name',oshostname(),'np',2)
-md=solve(md,ThermalSolutionEnum())
-
-
-# plot results
-
-comp_temp=md.results.ThermalSolution.Temperature
-relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100
-relative(find(comp_temp==md.initialization.temperature))=0
-plotmodel()(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,\
-	'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,\
-	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,\
-	'figposition','mathieu','FontSize#all',20)
-if printingflag, 
-	set(gcf,'Color','w')
-	printmodel()('thermalconduction','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off')
-	system(['mv thermalconduction.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal '])
-end
-
-
-# Fields and tolerances to track changes
-
-field_names     =['ConductionTemperature']
-field_tolerances=[1e-13]
-field_values    =[comp_temp]
Index: issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1304.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/InNeedOfDebugging/test1304.py	(revision 14101)
+++ 	(revision )
@@ -1,82 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test1304.m
-Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-from numpy import *
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-#  This file can be run to check that the geothermal flux in simple conduction is correctly modeled.
-
-#  There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
-
-#  at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower
-
-#  surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed...
-
-printingflag=false
-
-
-md=model()
-md=triangle(md,'../Exp/Square.exp',100000)
-md=setmask(md,'','')
-md=parameterize(md,'../Par/SquareThermal.py')
-md.extrude(11,1.)
-md=setflowequation(md,'Pattyn','all')
-
-
-pos2=numpy.nonzero(md.mesh.elementonsurface) md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0
-md.initialization.pressure=zeros(md.mesh.numberofvertices,1)
-md.basalforcings.geothermalflux(:)=0.1 %100mW/m^2
-
-
-# analytical results
-
-# the result is linear with depth and is equal to 0 on the upper surface (See BC)
-
-# d2T/dz2=0  -k*dT/dz(bed)=G  T(surface)=0  => T=-G/k*(z-surface)
-
-md.initialization.temperature=-0.1/md.materials.thermalconductivity*(md.mesh.z-md.geometry.surface) %G=0.1 W/m2
-
-
-# modeled  results
-
-md.cluster=generic('name',oshostname(),'np',2)
-md=solve(md,ThermalSolutionEnum())
-
-
-# plot results
-
-comp_temp=md.results.ThermalSolution.Temperature
-relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100
-relative(find(comp_temp==md.initialization.temperature))=0
-plotmodel()(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,\
-	'title','Analytical temperature','view',3,'data',comp_temp-md.initialization.temperature,\
-	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,\
-	'figposition','mathieu','FontSize#all',20)
-if printingflag, 
-	set(gcf,'Color','w')
-	printmodel()('thermalgeothermalflux','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off')
-	system(['mv thermalgeothermalflux.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal '])
-end
-
-
-# Fields and tolerances to track changes
-
-field_names     =['GeothermalFluxTemperature']
-field_tolerances=[1e-13]
-field_values    =[comp_temp]
Index: issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt
===================================================================
--- issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt	(revision 14101)
+++ issm/trunk-jpl/test/NightlyRun/python_skipped_tests.txt	(revision 14102)
@@ -11,2 +11,4 @@
 test418    needs Dakota
 test420    needs Dakota
+test1401    roundoff error in metric causes different meshes from matlab
+test1402    roundoff error in metric causes different meshes from matlab
Index: issm/trunk-jpl/test/NightlyRun/test1301.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1301.m	(revision 14101)
+++ issm/trunk-jpl/test/NightlyRun/test1301.m	(revision 14102)
@@ -10,5 +10,5 @@
 md=parameterize(md,'../Par/SquareThermal.par');
 md=extrude(md,3,2.);
-md=setflowequation(md,'Pattyn','all');
+md=setflowequation(md,'pattyn','all');
 
 %Some conditions specific to melting test
@@ -24,5 +24,5 @@
 melting=md.basalforcings.geothermalflux/(md.materials.rho_ice*md.materials.latentheat)*md.constants.yts;
 
-%modeled  results
+%modeled results
 md.cluster=generic('name',oshostname(),'np',2);
 md=solve(md,ThermalSolutionEnum());
@@ -30,10 +30,10 @@
 %plot results
 comp_melting=md.results.ThermalSolution.BasalforcingsMeltingRate;
-relative=abs((comp_melting-melting)./melting)*100;
-relative(find(comp_melting==melting))=0;
+relative=abs((comp_melting-melting)./melting)*100.;
+relative(find(comp_melting==melting))=0.;
 plotmodel(md,'data',comp_melting,'title','Modeled melting','data',melting,'title','Analytical melting',...
 	'data',comp_melting-melting,'title','Absolute error','data',relative,'title','Relative error [%]',...
 	'layer#all',1,'caxis#2',[1.02964 1.02966]*10^-4,'FontSize#all',20,'figposition','mathieu')
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('thermalmelting','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
@@ -41,5 +41,4 @@
 end
 
-
 %Fields and tolerances to track changes
 field_names     ={'BasalMelting'};
Index: issm/trunk-jpl/test/NightlyRun/test1301.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1301.py	(revision 14102)
+++ issm/trunk-jpl/test/NightlyRun/test1301.py	(revision 14102)
@@ -0,0 +1,60 @@
+import numpy
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from paterson import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This file can be run to check that the melting in simple conduction is correctly modeled.
+There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
+at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower
+surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed...
+"""
+
+printingflag=False
+
+md=model()
+md=triangle(md,'../Exp/Square.exp',100000.)
+md=setmask(md,'','')
+md=parameterize(md,'../Par/SquareThermal.py')
+md.extrude(3,2.)
+md=setflowequation(md,'pattyn','all')
+
+#Some conditions specific to melting test
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1),int)
+md.initialization.temperature=273.15*numpy.ones((md.mesh.numberofvertices,1))
+pos=numpy.nonzero(md.mesh.vertexonsurface)[0]
+md.thermal.spctemperature[pos]=md.initialization.temperature[pos]
+md.materials.rheology_B=paterson(md.initialization.temperature)
+
+#analytical results
+#melting heat = geothermal flux
+#Mb*L*rho=G   => Mb=G/L*rho
+melting=md.basalforcings.geothermalflux/(md.materials.rho_ice*md.materials.latentheat)*md.constants.yts
+
+#modeled results
+md.cluster=generic('name',oshostname(),'np',2)
+md=solve(md,ThermalSolutionEnum())
+
+#plot results
+comp_melting=md.results.ThermalSolution.BasalforcingsMeltingRate
+relative=numpy.abs((comp_melting-melting)/melting)*100.
+relative[numpy.nonzero(comp_melting==melting)[0]]=0.
+#plotmodel(md,'data',comp_melting,'title','Modeled melting','data',melting,'title','Analytical melting',...
+#	'data',comp_melting-melting,'title','Absolute error','data',relative,'title','Relative error [%]',...
+#	'layer#all',1,'caxis#2',[1.02964 1.02966]*10^-4,'FontSize#all',20,'figposition','mathieu')
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('thermalmelting','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
+#	system(['mv thermalmelting.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']);
+
+#Fields and tolerances to track changes
+field_names     =['BasalMelting']
+field_tolerances=[1e-08]
+field_values    =[comp_melting]
Index: issm/trunk-jpl/test/NightlyRun/test1302.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1302.m	(revision 14101)
+++ issm/trunk-jpl/test/NightlyRun/test1302.m	(revision 14102)
@@ -9,9 +9,9 @@
 md=parameterize(md,'../Par/SquareThermal.par');
 md=extrude(md,30,1.);   %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%)
-md=setflowequation(md,'Pattyn','all'); 
+md=setflowequation(md,'pattyn','all');
 
 %Thermal boundary conditions
-pos1=find(md.mesh.elementonbed);     md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10;
-pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0;
+pos1=find(md.mesh.elementonbed);     md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10.;
+pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0.;
 md.initialization.vz=0.1*ones(md.mesh.numberofvertices,1);
 md.initialization.vel=sqrt( md.initialization.vx.^2+ md.initialization.vy.^2+ md.initialization.vz.^2);
@@ -22,9 +22,9 @@
 %d2T/dz2-w*rho_ice*c/k*dT/dz=0   T(surface)=0  T(bed)=10   => T=A exp(alpha z)+B
 alpha=0.1/md.constants.yts*md.materials.rho_ice*md.materials.heatcapacity/md.materials.thermalconductivity;   %alpha=w rho_ice c /k  and w=0.1m/an
-A=10/(exp(alpha*(-1000))-1);    %A=T(bed)/(exp(alpha*bed)-1)  with bed=-1000 T(bed)=10
+A=10./(exp(alpha*(-1000.))-1.);    %A=T(bed)/(exp(alpha*bed)-1)  with bed=-1000 T(bed)=10
 B=-A;
 md.initialization.temperature=A*exp(alpha*md.mesh.z)+B;
 
-%modeled  results
+%modeled results
 md.cluster=generic('name',oshostname(),'np',2);
 md=solve(md,ThermalSolutionEnum());
@@ -32,11 +32,11 @@
 %plot results
 comp_temp=md.results.ThermalSolution.Temperature;
-relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100;
-relative(find(comp_temp==md.initialization.temperature))=0;
+relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100.;
+relative(find(comp_temp==md.initialization.temperature))=0.;
 plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,...
 	'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,...
 	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,...
 	'figposition','mathieu','FontSize#all',20)
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('thermaladvection','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
Index: issm/trunk-jpl/test/NightlyRun/test1302.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1302.py	(revision 14102)
+++ issm/trunk-jpl/test/NightlyRun/test1302.py	(revision 14102)
@@ -0,0 +1,65 @@
+import numpy
+import sys
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This file can be run to check that the advection-diffusion  is correctly modeled.
+There is u=v=0 and w=cst everywhere the only thermal boundary conditions are an imposed temperature
+at upper surface and an impose flux at its base.
+"""
+
+printingflag=False
+
+md=model()
+md=triangle(md,'../Exp/Square.exp',100000.)
+md=setmask(md,'','')
+md=parameterize(md,'../Par/SquareThermal.py')
+md.extrude(30,1.)    #NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%)
+md=setflowequation(md,'pattyn','all')
+
+#Thermal boundary conditions
+pos1=numpy.nonzero(md.mesh.elementonbed)[0]
+md.thermal.spctemperature[md.mesh.elements[pos1,0:3]-1]=10.
+pos2=numpy.nonzero(md.mesh.elementonsurface)[0]
+md.thermal.spctemperature[md.mesh.elements[pos2,3:6]-1]=0.
+md.initialization.vz=0.1*numpy.ones((md.mesh.numberofvertices,1))
+md.initialization.vel=numpy.sqrt( md.initialization.vx**2+ md.initialization.vy**2+ md.initialization.vz**2)
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1),int)
+
+md.thermal.stabilization=2
+#analytical results
+#d2T/dz2-w*rho_ice*c/k*dT/dz=0   T(surface)=0  T(bed)=10   => T=A exp(alpha z)+B
+alpha=0.1/md.constants.yts*md.materials.rho_ice*md.materials.heatcapacity/md.materials.thermalconductivity    #alpha=w rho_ice c /k  and w=0.1m/an
+A=10./(numpy.exp(alpha*(-1000.))-1.)    #A=T(bed)/(exp(alpha*bed)-1)  with bed=-1000 T(bed)=10
+B=-A
+md.initialization.temperature=A*numpy.exp(alpha*md.mesh.z)+B
+
+#modeled results
+md.cluster=generic('name',oshostname(),'np',2)
+md=solve(md,ThermalSolutionEnum())
+
+#plot results
+comp_temp=md.results.ThermalSolution.Temperature
+relative=numpy.abs((comp_temp-md.initialization.temperature)/md.initialization.temperature)*100.
+relative[numpy.nonzero(comp_temp==md.initialization.temperature)[0]]=0.
+#plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,...
+#	'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,...
+#	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,...
+#	'figposition','mathieu','FontSize#all',20)
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('thermaladvection','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
+#	system(['mv thermaladvection.png ' ISSM_DIR '/website/doc_pdf/validation/Images/EISMINT ']);
+
+#Fields and tolerances to track changes
+field_names     =['AdvectionTemperature']
+field_tolerances=[1e-13]
+field_values    =[comp_temp]
Index: issm/trunk-jpl/test/NightlyRun/test1303.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1303.m	(revision 14101)
+++ issm/trunk-jpl/test/NightlyRun/test1303.m	(revision 14102)
@@ -10,7 +10,8 @@
 md=parameterize(md,'../Par/SquareThermal.par');
 md=extrude(md,11,2.);
-md=setflowequation(md,'Pattyn','all');
-pos1=find(md.mesh.elementonbed);     md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10;
-pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0;
+md=setflowequation(md,'pattyn','all');
+
+pos1=find(md.mesh.elementonbed);     md.thermal.spctemperature(md.mesh.elements(pos1,1:3))=10.;
+pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0.;
 md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
 
@@ -18,7 +19,7 @@
 %d2T/dz2=0 T(bed)=10 T(surface)=0  => T=0*(z-bed)/thickness+10*(surface-z)/thickness
 %each layer of the 3d mesh must have a constant value
-md.initialization.temperature=10*(md.geometry.surface-md.mesh.z)./md.geometry.thickness;
+md.initialization.temperature=10.*(md.geometry.surface-md.mesh.z)./md.geometry.thickness;
 
-%modeled  results
+%modeled results
 md.cluster=generic('name',oshostname(),'np',2);
 md=solve(md,ThermalSolutionEnum());
@@ -26,11 +27,11 @@
 %plot results
 comp_temp=md.results.ThermalSolution.Temperature;
-relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100;
-relative(find(comp_temp==md.initialization.temperature))=0;
+relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100.;
+relative(find(comp_temp==md.initialization.temperature))=0.;
 plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,...
 	'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,...
 	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,...
 	'figposition','mathieu','FontSize#all',20)
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('thermalconduction','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
Index: issm/trunk-jpl/test/NightlyRun/test1303.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1303.py	(revision 14102)
+++ issm/trunk-jpl/test/NightlyRun/test1303.py	(revision 14102)
@@ -0,0 +1,60 @@
+import numpy
+import sys
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This file can be run to check that the conduction is correctly modeled.
+There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
+at the lower and upper surface. The result must be a linear temperature from the upper to the lower
+surface. if it is not the case, something is thermal modeling has been changed...
+"""
+
+printingflag=False
+
+md=model()
+md=triangle(md,'../Exp/Square.exp',100000.)
+md=setmask(md,'all','')
+md=parameterize(md,'../Par/SquareThermal.py')
+md.extrude(11,2.)
+md=setflowequation(md,'pattyn','all')
+
+pos1=numpy.nonzero(md.mesh.elementonbed)[0]
+md.thermal.spctemperature[md.mesh.elements[pos1,0:3]-1]=10.
+pos2=numpy.nonzero(md.mesh.elementonsurface)[0]
+md.thermal.spctemperature[md.mesh.elements[pos2,3:6]-1]=0.
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1),int)
+
+#analytical results
+#d2T/dz2=0 T(bed)=10 T(surface)=0  => T=0*(z-bed)/thickness+10*(surface-z)/thickness
+#each layer of the 3d mesh must have a constant value
+md.initialization.temperature=10.*(md.geometry.surface-md.mesh.z)/md.geometry.thickness
+
+#modeled results
+md.cluster=generic('name',oshostname(),'np',2)
+md=solve(md,ThermalSolutionEnum())
+
+#plot results
+comp_temp=md.results.ThermalSolution.Temperature
+relative=numpy.abs((comp_temp-md.initialization.temperature)/md.initialization.temperature)*100.
+relative[numpy.nonzero(comp_temp==md.initialization.temperature)[0]]=0.
+#plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,...
+#	'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,...
+#	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,...
+#	'figposition','mathieu','FontSize#all',20)
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('thermalconduction','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
+#	system(['mv thermalconduction.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']);
+
+#Fields and tolerances to track changes
+field_names     =['ConductionTemperature']
+field_tolerances=[1e-13]
+field_values    =[comp_temp]
Index: issm/trunk-jpl/test/NightlyRun/test1304.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1304.m	(revision 14101)
+++ issm/trunk-jpl/test/NightlyRun/test1304.m	(revision 14102)
@@ -10,7 +10,7 @@
 md=parameterize(md,'../Par/SquareThermal.par');
 md=extrude(md,11,1.);
-md=setflowequation(md,'Pattyn','all');
+md=setflowequation(md,'pattyn','all');
 
-pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0;
+pos2=find(md.mesh.elementonsurface); md.thermal.spctemperature(md.mesh.elements(pos2,4:6))=0.;
 md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
 md.basalforcings.geothermalflux(:)=0.1; %100mW/m^2
@@ -21,5 +21,5 @@
 md.initialization.temperature=-0.1/md.materials.thermalconductivity*(md.mesh.z-md.geometry.surface); %G=0.1 W/m2
 
-%modeled  results
+%modeled results
 md.cluster=generic('name',oshostname(),'np',2);
 md=solve(md,ThermalSolutionEnum());
@@ -27,11 +27,11 @@
 %plot results
 comp_temp=md.results.ThermalSolution.Temperature;
-relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100;
-relative(find(comp_temp==md.initialization.temperature))=0;
+relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100.;
+relative(find(comp_temp==md.initialization.temperature))=0.;
 plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,...
 	'title','Analytical temperature','view',3,'data',comp_temp-md.initialization.temperature,...
 	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,...
 	'figposition','mathieu','FontSize#all',20)
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('thermalgeothermalflux','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
Index: issm/trunk-jpl/test/NightlyRun/test1304.py
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1304.py	(revision 14102)
+++ issm/trunk-jpl/test/NightlyRun/test1304.py	(revision 14102)
@@ -0,0 +1,59 @@
+import numpy
+import sys
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+"""
+This file can be run to check that the geothermal flux in simple conduction is correctly modeled.
+There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
+at upper surface and an impose flux at its base. The result must be a linear temperature from the upper to the lower
+surface with an imposed slope (Geothermal flux). if it is not the case, something is thermal modeling has been changed...
+"""
+
+printingflag=False
+
+md=model()
+md=triangle(md,'../Exp/Square.exp',100000.)
+md=setmask(md,'','')
+md=parameterize(md,'../Par/SquareThermal.py')
+md.extrude(11,1.)
+md=setflowequation(md,'pattyn','all')
+
+pos2=numpy.nonzero(md.mesh.elementonsurface)[0]
+md.thermal.spctemperature[md.mesh.elements[pos2,3:6]-1]=0.
+md.initialization.pressure=numpy.zeros((md.mesh.numberofvertices,1),int)
+md.basalforcings.geothermalflux[:]=0.1    #100mW/m^2
+
+#analytical results
+#the result is linear with depth and is equal to 0 on the upper surface (See BC)
+#d2T/dz2=0  -k*dT/dz(bed)=G  T(surface)=0  => T=-G/k*(z-surface)
+md.initialization.temperature=-0.1/md.materials.thermalconductivity*(md.mesh.z-md.geometry.surface)    #G=0.1 W/m2
+
+#modeled results
+md.cluster=generic('name',oshostname(),'np',2)
+md=solve(md,ThermalSolutionEnum())
+
+#plot results
+comp_temp=md.results.ThermalSolution.Temperature
+relative=numpy.abs((comp_temp-md.initialization.temperature)/md.initialization.temperature)*100.
+relative[numpy.nonzero(comp_temp==md.initialization.temperature)[0]]=0.
+#plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data',md.initialization.temperature,'view',3,...
+#	'title','Analytical temperature','view',3,'data',comp_temp-md.initialization.temperature,...
+#	'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,...
+#	'figposition','mathieu','FontSize#all',20)
+if printingflag:
+	pass
+#	set(gcf,'Color','w')
+#	printmodel('thermalgeothermalflux','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
+#	system(['mv thermalgeothermalflux.png ' ISSM_DIR '/website/doc_pdf/validation/Images/Thermal ']);
+
+#Fields and tolerances to track changes
+field_names     =['GeothermalFluxTemperature']
+field_tolerances=[1e-13]
+field_values    =[comp_temp]
Index: issm/trunk-jpl/test/NightlyRun/test1401.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1401.m	(revision 14101)
+++ issm/trunk-jpl/test/NightlyRun/test1401.m	(revision 14102)
@@ -1,8 +1,8 @@
 %test the anisotropic mesh adaptation
-%function to capture = exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2;
+%function to capture = exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.;
 printingflag=false;
 
 %create square mesh
-L=1; %in m
+L=1.; %in m
 nx=70; %numberof nodes in x direction
 ny=70;
@@ -11,7 +11,7 @@
 %mesh adaptation loop YAMS
 md=squaremesh(md,L,L,nx,ny);
-md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2;
+md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.;
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh1_yams1','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -19,8 +19,8 @@
 end
 
-md=YamsCall(md,md.inversion.vel_obs,0.001,0.3,1.3,10^-4);
-md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2;
+md=YamsCall(md,md.inversion.vel_obs,0.001,0.3,1.3,10.^-4);
+md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.;
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh1_yams2','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -29,7 +29,7 @@
 
 md=YamsCall(md,md.inversion.vel_obs,0.001,0.3,2.5,0.008);
-md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2;
+md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.;
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh1_yams3','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -41,7 +41,7 @@
 %mesh adaptation loop BAMG
 md=squaremesh(md,L,L,nx,ny);
-md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2;
+md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.;
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh1_bamg1','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -50,8 +50,8 @@
 
 md.private.bamg=NaN;
-md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',1.3,'err',10^-4);
-md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2;
+md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',1.3,'err',10.^-4);
+md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.;
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh1_bamg2','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -61,7 +61,7 @@
 md.private.bamg=NaN;
 md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',2.5,'err',0.008);
-md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2;
+md.inversion.vel_obs=exp(-(sqrt((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)-0.75).^2*10.^6)+((md.mesh.x+0.1).^2+(md.mesh.y+0.1).^2)/2.;
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh1_bamg3','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
Index: issm/trunk-jpl/test/NightlyRun/test1402.m
===================================================================
--- issm/trunk-jpl/test/NightlyRun/test1402.m	(revision 14101)
+++ issm/trunk-jpl/test/NightlyRun/test1402.m	(revision 14102)
@@ -3,5 +3,5 @@
 
 %create square mesh
-L=1; %in m
+L=1.; %in m
 nx=30; %numberof nodes in x direction
 ny=30;
@@ -10,10 +10,10 @@
 %mesh adaptation loop YAMS
 md=squaremesh(md,L,L,nx,ny);
-u=4*md.mesh.x-2; v=4*md.mesh.y-2;
-md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ...
-	+tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ...
-	+tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ;
+u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.;
+md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ...
+	+tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ...
+	+tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25));
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh2_yams1','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -21,11 +21,11 @@
 end
 
-md=YamsCall(md,md.inversion.vel_obs,0.005,0.3,2.3,10^-2);
-u=4*md.mesh.x-2; v=4*md.mesh.y-2;
-md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ...
-	+tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ...
-	+tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ;
+md=YamsCall(md,md.inversion.vel_obs,0.005,0.3,2.3,10.^2);
+u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.;
+md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ...
+	+tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ...
+	+tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25));
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh2_yams2','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -34,10 +34,10 @@
 
 md=YamsCall(md,md.inversion.vel_obs,0.005,0.3,3,0.005);
-u=4*md.mesh.x-2; v=4*md.mesh.y-2;
-md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ...
-	+tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ...
-	+tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ;
+u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.;
+md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ...
+	+tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ...
+	+tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25));
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh2_yams3','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -49,10 +49,10 @@
 %mesh adaptation loop BAMG
 md=squaremesh(md,L,L,nx,ny);
-u=4*md.mesh.x-2; v=4*md.mesh.y-2;
-md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ...
-	+tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ...
-	+tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ;
+u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.;
+md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ...
+	+tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ...
+	+tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25));
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh2_bamg1','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -61,11 +61,11 @@
 
 md.private.bamg=NaN;
-md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',2.3,'err',10^-2);
-u=4*md.mesh.x-2; v=4*md.mesh.y-2;
-md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ...
-	+tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ...
-	+tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ;
+md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',2.3,'err',10.^2);
+u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.;
+md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ...
+	+tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ...
+	+tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25));
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh2_bamg2','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -75,10 +75,10 @@
 md.private.bamg=NaN;
 md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',3,'err',0.005);
-u=4*md.mesh.x-2; v=4*md.mesh.y-2;
-md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ...
-	+tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ...
-	+tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ;
+u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.;
+md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ...
+	+tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ...
+	+tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25));
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh2_bamg3','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
@@ -88,10 +88,10 @@
 md.private.bamg=NaN;
 md=bamg(md,'field',md.inversion.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',1.5,'err',0.003,'anisomax',1);
-u=4*md.mesh.x-2; v=4*md.mesh.y-2;
-md.inversion.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ...
-	+tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ...
-	+tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ;
+u=4.*md.mesh.x-2.; v=4.*md.mesh.y-2.;
+md.inversion.vel_obs=tanh(30.*(u.^2+v.^2-0.25)) ...
+	+tanh(30.*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u-0.75).^2+(v+0.75).^2-0.25)) ...
+	+tanh(30.*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30.*((u+0.75).^2+(v+0.75).^2-0.25));
 plotmodel(md,'data',md.inversion.vel_obs,'data',md.inversion.vel_obs,'nlines',1,'ncols',2,'title','','figposition',[500 500 1000 500],'axis#all','equal','xlim#all',[0 1],'ylim#all',[0 1],'edgecolor#1','w'); pause(0.5);
-if printingflag, 
+if printingflag,
 	set(gcf,'Color','w')
 	printmodel('mesh2_bamgiso','png','margin','on','marginsize',25,'frame','off','resolution',1,'hardcopy','off');
Index: issm/trunk-jpl/test/Par/Pig.py
===================================================================
--- issm/trunk-jpl/test/Par/Pig.py	(revision 14101)
+++ issm/trunk-jpl/test/Par/Pig.py	(revision 14102)
@@ -48,5 +48,5 @@
 md.prognostic.stabilization=1.
 md.verbose=verbose(0)
-md.settings.waitonlock=30.
+md.settings.waitonlock=30
 md.timestepping.time_step=1.
 md.timestepping.final_time=2.
Index: issm/trunk-jpl/test/Par/SquareSheetConstrained.py
===================================================================
--- issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 14101)
+++ issm/trunk-jpl/test/Par/SquareSheetConstrained.py	(revision 14102)
@@ -51,5 +51,5 @@
 md.thermal.stabilization=1.
 md.verbose=verbose(0)
-md.settings.waitonlock=30.
+md.settings.waitonlock=30
 md.diagnostic.restol=0.05
 md.steadystate.reltol=0.05
Index: issm/trunk-jpl/test/Par/SquareSheetShelf.py
===================================================================
--- issm/trunk-jpl/test/Par/SquareSheetShelf.py	(revision 14101)
+++ issm/trunk-jpl/test/Par/SquareSheetShelf.py	(revision 14102)
@@ -57,5 +57,5 @@
 md.thermal.stabilization=1
 md.verbose=verbose(0)
-md.settings.waitonlock=30.
+md.settings.waitonlock=30
 md.diagnostic.restol=0.05
 md.steadystate.reltol=0.05
Index: issm/trunk-jpl/test/Par/SquareShelf.py
===================================================================
--- issm/trunk-jpl/test/Par/SquareShelf.py	(revision 14101)
+++ issm/trunk-jpl/test/Par/SquareShelf.py	(revision 14102)
@@ -70,5 +70,5 @@
 md.prognostic.stabilization = 1.
 md.thermal.stabilization = 1.
-md.settings.waitonlock = 30.
+md.settings.waitonlock = 30
 md.verbose=verbose()
 md.diagnostic.restol = 0.10
Index: issm/trunk-jpl/test/Par/SquareShelfConstrained.py
===================================================================
--- issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 14101)
+++ issm/trunk-jpl/test/Par/SquareShelfConstrained.py	(revision 14102)
@@ -55,5 +55,5 @@
 md.thermal.stabilization=1.
 md.verbose = verbose(0)
-md.settings.waitonlock=30.
+md.settings.waitonlock=30
 md.diagnostic.restol=0.05
 md.diagnostic.reltol=0.05
Index: issm/trunk-jpl/test/Par/SquareThermal.par
===================================================================
--- issm/trunk-jpl/test/Par/SquareThermal.par	(revision 14101)
+++ issm/trunk-jpl/test/Par/SquareThermal.par	(revision 14102)
@@ -4,7 +4,7 @@
 
 disp('      creating thickness');
-h=1000;
+h=1000.;
 md.geometry.thickness=h*ones(md.mesh.numberofvertices,1);
-md.geometry.bed=-1000*ones(md.mesh.numberofvertices,1);
+md.geometry.bed=-1000.*ones(md.mesh.numberofvertices,1);
 md.geometry.surface=md.geometry.bed+md.geometry.thickness;
 
@@ -15,21 +15,21 @@
 
 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 temperatures');
-md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
+md.initialization.temperature=(273.-20.)*ones(md.mesh.numberofvertices,1);
 
-disp('      creating flow law paramter');
+disp('      creating flow law parameter');
 md.materials.rheology_B=paterson(md.initialization.temperature);
-md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_n=3.*ones(md.mesh.numberofelements,1);
 
 disp('      creating surface mass balance');
 md.surfaceforcings.mass_balance=ones(md.mesh.numberofvertices,1)/md.constants.yts; %1m/a
-md.basalforcings.melting_rate=0*ones(md.mesh.numberofvertices,1)/md.constants.yts; %1m/a
+md.basalforcings.melting_rate=0.*ones(md.mesh.numberofvertices,1)/md.constants.yts; %1m/a
 
 %Deal with boundary conditions:
@@ -41,3 +41,3 @@
 md.thermal.spctemperature(:)=md.initialization.temperature;
 md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1); 
-pos=find(md.mask.elementongroundedice);md.basalforcings.geothermalflux(md.mesh.elements(pos,:))=1*10^-3; %1 mW/m^2
+pos=find(md.mask.elementongroundedice);md.basalforcings.geothermalflux(md.mesh.elements(pos,:))=1.*10^-3; %1 mW/m^2
Index: issm/trunk-jpl/test/Par/SquareThermal.py
===================================================================
--- issm/trunk-jpl/test/Par/SquareThermal.py	(revision 14102)
+++ issm/trunk-jpl/test/Par/SquareThermal.py	(revision 14102)
@@ -0,0 +1,48 @@
+import numpy
+from paterson import *
+from SetMarineIceSheetBC import *
+
+#Ok, start defining model parameters here
+
+md.timestepping.time_step=0
+
+print "      creating thickness"
+h=1000.
+md.geometry.thickness=h*numpy.ones((md.mesh.numberofvertices,1))
+md.geometry.bed=-1000.*numpy.ones((md.mesh.numberofvertices,1))
+md.geometry.surface=md.geometry.bed+md.geometry.thickness;
+
+print "      creating velocities"
+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))
+
+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 temperatures"
+md.initialization.temperature=(273.-20.)*numpy.ones((md.mesh.numberofvertices,1))
+
+print "      creating flow law parameter"
+md.materials.rheology_B=paterson(md.initialization.temperature)
+md.materials.rheology_n=3.*numpy.ones((md.mesh.numberofelements,1))
+
+print "      creating surface mass balance"
+md.surfaceforcings.mass_balance=numpy.ones((md.mesh.numberofvertices,1))/md.constants.yts    #1m/a
+md.basalforcings.melting_rate=0.*numpy.ones((md.mesh.numberofvertices,1))/md.constants.yts    #1m/a
+
+#Deal with boundary conditions:
+
+print "      boundary conditions for diagnostic model"
+md=SetMarineIceSheetBC(md,'../Exp/SquareFront.exp')
+
+print "      boundary conditions for thermal model"
+md.thermal.spctemperature[:]=md.initialization.temperature
+md.basalforcings.geothermalflux=numpy.zeros((md.mesh.numberofvertices,1)) 
+pos=numpy.nonzero(md.mask.elementongroundedice)[0]
+md.basalforcings.geothermalflux[md.mesh.elements[pos,:]-1]=1.*10**-3    #1 mW/m^2
