Index: issm/trunk/test/NightlyRun/runme.py
===================================================================
--- issm/trunk/test/NightlyRun/runme.py	(revision 23189)
+++ issm/trunk/test/NightlyRun/runme.py	(revision 23394)
@@ -11,5 +11,5 @@
 	"""
 	RUNME - test deck for ISSM nightly runs
- 
+
 	    In a test deck directory (tests/Vertification/NightlyRun for example)
 	    The following command will launch all the existing tests:
@@ -18,5 +18,5 @@
 	    >> runme(id=[101,102])
 	    etc...
- 
+
 	    Available options:
 	       'id'            followed by the list of ids or (parts of) test names requested
@@ -37,8 +37,8 @@
 	       'procedure'     'check' : run the test (default)
 	                       'update': update the archive
- 
+
 	    Usage:
 	       runme(varargin)
- 
+
 	    Examples:
 	       runme()
@@ -49,5 +49,4 @@
 	       runme(id=[[101,102],['Dakota','Slr']])
 	"""
-
 	from parallelrange import parallelrange
 	from IdToName import IdToName
@@ -82,17 +81,14 @@
 	flist=glob('test*.py')    #File name must start with 'test' and must end by '.py' and must be different than 'test.py'
 	list_ids=[int(file[4:-3]) for file in flist if not file == 'test.py']    #Keep test id only (skip 'test' and '.py')
-	#print 'list_ids =',list_ids
 
 	i1,i2=parallelrange(rank,numprocs,len(list_ids))    #Get tests for this cpu only
 	list_ids=list_ids[i1:i2+1]
-	#print 'list_ids after parallelrange =',list_ids
-	
-	if id:
+
+	if np.size(id) > 0 and not id==None:
 		test_ids = set(GetIds(id)).intersection(set(list_ids))
 	else:
 		# if no tests are specifically provided, do them all
 		test_ids = set(list_ids)
-		
-	#print 'test_ids after list =',test_ids
+
 	# }}}
 	#GET exclude {{{
@@ -100,6 +96,4 @@
 
 	test_ids=test_ids.difference(exclude_ids)
-	#print 'test_ids after exclude =',sorted(test_ids)
-	#return
 	# }}}
 	#Process Ids according to benchmarks {{{
@@ -124,8 +118,6 @@
 	elif benchmark=='adolc':
 		test_ids=test_ids.intersection(set(range(3001,3200)))
-	#print 'test_ids after benchmark =',test_ids
 	test_ids=list(test_ids)
 	test_ids.sort()
-	#print 'test_ids after sort =',test_ids
 	# }}}
 
@@ -190,5 +182,5 @@
 							if np.shape(field) != np.shape(archive):
 								raise RuntimeError("Field '"+archive_name+"' from test is malformed; shape is "+str(np.shape(field.T))+", should be "+str(np.shape(archive))+" (or "+str(np.shape(archive.T))+").")
-						
+
 						error_diff=np.amax(np.abs(archive-field),axis=0)/(np.amax(np.abs(archive),axis=0)+float_info.epsilon)
 
Index: issm/trunk/test/NightlyRun/test2084.m
===================================================================
--- issm/trunk/test/NightlyRun/test2084.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test2084.m	(revision 23394)
@@ -1,4 +1,3 @@
-
-%Test Name: GiaCaron. 
+%Test Name: GiaCaron
 %Forward Love number solution for a viscoelastic earth, model M3-L70-V01 from
 %Spada, G., Barletta, V. R., Klemann, V., Riva, R. E. M., Martinec, Z.,
@@ -39,5 +38,5 @@
 %loading love numbers
 field_names     ={'LoveH_loading_elastic','LoveK_loading_elastic','LoveL_loading_elastic'};
-field_tolerances={7.0e-10,7.0e-10,7.0e-10};
+field_tolerances={4.3e-9,4.3e-9,4.3e-9};
 field_values={...
 	(md.results.LoveSolution.LoveHr(:,1)),...
@@ -57,5 +56,5 @@
 %Fields and tolerances to track changes
 field_names     ={field_names{:},'LoveH_loading_realpart','LoveK_loading_realpart','LoveL_loading_realpart','LoveH_loading_imagpart','LoveK_loading_imagpart','LoveL_loading_imagpart'};
-field_tolerances={field_tolerances{:},7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10};
+field_tolerances={field_tolerances{:},5e-7,5e-7,5e-7,5e-7,5e-7,5e-7};
 field_values={field_values{:},...
 	(md.results.LoveSolution.LoveHr(:,:)),...
@@ -75,5 +74,5 @@
 %tidal love numbers
 field_names     ={field_names{:},'LoveH_tidal_elastic','LoveK_tidal_elastic','LoveL_tidal_elastic','LoveH_tidal_realpart','LoveK_tidal_realpart','LoveL_tidal_realpart','LoveH_tidal_imagpart','LoveK_tidal_imagpart','LoveL_tidal_imagpart'};
-field_tolerances={field_tolerances{:},7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10};
+field_tolerances={field_tolerances{:},8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6};
 field_values={field_values{:},...
 	(md.results.LoveSolution.LoveHr(:,1)),...
@@ -115,5 +114,5 @@
 %
 %field_names     ={field_names{:},'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart'};
-%field_tolerances={field_tolerances{:},7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10};
+%field_tolerances={field_tolerances{:},2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6};
 %field_values={field_values{:},...
 %	(md.results.LoveSolution.LoveHr(:,1)),...
@@ -147,5 +146,5 @@
 
 field_names     ={field_names{:},'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart'};
-field_tolerances={field_tolerances{:},7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10,7.0e-10};
+field_tolerances={field_tolerances{:},2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6};
 field_values={field_values{:},...
 	(md.results.LoveSolution.LoveHr(:,1)),...
Index: issm/trunk/test/NightlyRun/test2084.py
===================================================================
--- issm/trunk/test/NightlyRun/test2084.py	(revision 23189)
+++ issm/trunk/test/NightlyRun/test2084.py	(revision 23394)
@@ -1,8 +1,9 @@
-#Test Name: GiaCaron. Forward Love number solution for a viscoelastic earth,
-#	model M3-L70-V01 from Spada, G., Barletta, V. R., Klemann, V., Riva, R. E. M.,
-#	Martinec, Z., Gasperini, P., Lund, B., Wolf, D., Vermeersen, L. L. A.
-#	and King, M. A. (2011), A benchmark study for glacial isostatic
-#	adjustment codes. Geophysical Journal International,
-#	185: 106-132. doi:10.1111/j.1365-246X.2011.04952.x
+#Test Name: GiaCaron
+#Forward Love number solution for a viscoelastic earth,
+#model M3-L70-V01 from Spada, G., Barletta, V. R., Klemann, V., Riva, R. E. M.,
+#Martinec, Z., Gasperini, P., Lund, B., Wolf, D., Vermeersen, L. L. A.
+#and King, M. A. (2011), A benchmark study for glacial isostatic
+#adjustment codes. Geophysical Journal International,
+#185: 106-132. doi:10.1111/j.1365-246X.2011.04952.x
 
 from model import *
@@ -48,6 +49,6 @@
 #Fields and tolerances to track changes
 #loading love numbers
-field_names=['LoveHe','LoveKe','LoveLe']
-field_tolerances=[3.7e-9,3.7e-9,3.7e-9]
+field_names=['LoveH_loading_elastic','LoveK_loading_elastic','LoveL_loading_elastic']
+field_tolerances=[4.3e-9,4.3e-9,4.3e-9]
 field_values=[
 np.array(md.results.LoveSolution.LoveHr)[:,0],
@@ -67,6 +68,6 @@
 #Fields and tolerances to track changes
 #loading love numbers
-field_names+=['LoveHlr','LoveKlr','LoveLlr','LoveHli','LoveKli','LoveLli']
-field_tolerances+=[3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9]
+field_names+=['LoveH_loading_realpart','LoveK_loading_realpart','LoveL_loading_realpart','LoveH_loading_imagpart','LoveK_loading_imagpart','LoveL_loading_imagpart']
+field_tolerances+=[5e-7,5e-7,5e-7,5e-7,5e-7,5e-7]
 field_values+=[
 np.array(md.results.LoveSolution.LoveHr),
@@ -87,6 +88,6 @@
 
 #tidal love numbers, check
-field_names+=['LoveHtr','LoveKtr','LoveLtr','LoveHti','LoveKti','LoveLti']
-field_tolerances+=[3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9]
+field_names+=['LoveH_tidal_elastic','LoveK_tidal_elastic','LoveL_tidal_elastic','LoveH_tidal_realpart','LoveK_tidal_realpart','LoveL_tidal_realpart','LoveH_tidal_imagpart','LoveK_tidal_imagpart','LoveL_tidal_imagpart']
+field_tolerances+=[8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6]
 field_values+=[
 np.array(md.results.LoveSolution.LoveHr)[:,0],
@@ -127,6 +128,6 @@
 #md=solve(md,'lv')
 #
-#field_names=[field_names,'LoveHmr','LoveKmr','LoveLmr','LoveHmi','LoveKmi','LoveLmi']
-#field_tolerances=[field_tolerances,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
+#field_names=[field_names,'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart']
+#field_tolerances=[field_tolerances,4.3e-9,4.3e-9,4.3e-9,4.3e-9,4.3e-9,4.3e-9,4.3e-9,4.3e-9,4.3e-9]
 #field_values=[field_values,\
 #	(md.results.LoveSolution.LoveHr[:][0]),\
@@ -186,6 +187,6 @@
 md.love.nfreq=len(md.love.frequencies)
 
-field_names+=['LoveHmr','LoveKmr','LoveLmr','LoveHmi','LoveKmi','LoveLmi']
-field_tolerances+=[3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9,3.7e-9]
+field_names+=['LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart']
+field_tolerances+=[2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6]
 field_values+=[
 np.array(md.results.LoveSolution.LoveHr)[:,0],
Index: issm/trunk/test/NightlyRun/test211.js
===================================================================
--- issm/trunk/test/NightlyRun/test211.js	(revision 23189)
+++ issm/trunk/test/NightlyRun/test211.js	(revision 23394)
@@ -12,10 +12,10 @@
 //Fields and tolerances to track changes
 field_names=[
-	'Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsGroundediceMeltingRate1', 
-	'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Temperature2','BasalforcingsGroundediceMeltingRate2', 
+	'Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsGroundediceMeltingRate1',
+	'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Temperature2','BasalforcingsGroundediceMeltingRate2',
 	'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Temperature3','BasalforcingsGroundediceMeltingRate3'];
 field_tolerances=[
 	2e-08,2e-08,5e-05,2e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,
-	5e-06,5e-06,8e-05,5e-06,1e-07,5e-07,5e-07,5e-07,3e-06,5e-05,
+	5e-06,5e-06,8e-05,5e-06,5e-07,5e-07,5e-07,5e-07,3e-06,5e-05,
 	8e-06,8e-06,8e-05,8e-06,5e-07,8e-07,8e-07,8e-07,5e-06,8e-05];
 field_values=[
Index: issm/trunk/test/NightlyRun/test211.m
===================================================================
--- issm/trunk/test/NightlyRun/test211.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test211.m	(revision 23394)
@@ -16,5 +16,5 @@
 field_tolerances={...
 	2e-08,2e-08,5e-05,2e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,...
-	5e-06,5e-06,8e-05,5e-06,1e-07,5e-07,5e-07,5e-07,3e-06,5e-05,...
+	5e-06,5e-06,8e-05,5e-06,5e-07,5e-07,5e-07,5e-07,3e-06,5e-05,...
 	8e-06,8e-06,8e-05,8e-06,5e-07,8e-07,8e-07,8e-07,5e-06,8e-05};
 field_values={...
Index: issm/trunk/test/NightlyRun/test211.py
===================================================================
--- issm/trunk/test/NightlyRun/test211.py	(revision 23189)
+++ issm/trunk/test/NightlyRun/test211.py	(revision 23394)
@@ -26,6 +26,6 @@
 						 'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Temperature3','BasalforcingsGroundediceMeltingRate3']
 field_tolerances=[2e-08,2e-08,5e-05,2e-08,1e-08,1e-08,1e-08,1e-08,1e-08,1e-08,
-									5e-06,5e-06,8e-05,5e-06,1e-07,5e-07,5e-07,5e-07,3e-06,5e-05,
-									8e-06,8e-06,8e-05,8e-06,5e-07,8e-07,8e-07,8e-07,5e-06,8e-05]
+        5e-06,5e-06,8e-05,5e-06,5e-07,5e-07,5e-07,5e-07,3e-06,5e-05,
+        8e-06,8e-06,8e-05,8e-06,5e-07,8e-07,8e-07,8e-07,5e-06,8e-05]
 field_values=[md.results.TransientSolution[0].Vx,
 							md.results.TransientSolution[0].Vy,
Index: issm/trunk/test/NightlyRun/test215.m
===================================================================
--- issm/trunk/test/NightlyRun/test215.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test215.m	(revision 23394)
@@ -24,5 +24,5 @@
 %Fields and tolerances to track changes
 field_names     ={'Gradient','Misfits','MaterialsRheologyBbar','Pressure','Vel','Vx','Vy'};
-field_tolerances={4.6e-08,1e-08,2e-08,1e-09,2e-09,2e-08,2e-08};
+field_tolerances={4.6e-08,1e-08,2e-08,2e-09,3e-09,2e-08,2e-08};
 field_values={...
 	(md.results.StressbalanceSolution.Gradient1),...
Index: issm/trunk/test/NightlyRun/test215.py
===================================================================
--- issm/trunk/test/NightlyRun/test215.py	(revision 23189)
+++ issm/trunk/test/NightlyRun/test215.py	(revision 23394)
@@ -37,5 +37,5 @@
 # Fields and tolerances to track changes
 field_names     =['Gradient','Misfits','MaterialsRheologyBbar','Pressure','Vel','Vx','Vy']
-field_tolerances=[4.6e-08,1e-08,2e-09,1e-08,2e-09,2e-08,2e-08]
+field_tolerances=[4.6e-08,1e-08,2e-08,2e-09,3e-09,2e-08,2e-08]
 field_values=[md.results.StressbalanceSolution.Gradient1,
         md.results.StressbalanceSolution.J,
Index: issm/trunk/test/NightlyRun/test234.py
===================================================================
--- issm/trunk/test/NightlyRun/test234.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test234.py	(revision 23394)
@@ -0,0 +1,101 @@
+#Test Name: SquareShelfTranForceNeg2dDakotaSamp
+import numpy as np
+import scipy.io as spio
+from os import getcwd
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from SMBgemb import *
+from partitioner import *
+from dmeth_params_set import *
+
+md = triangle(model(),'../Exp/Square.exp',180000.)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+
+md.timestepping.time_step = 1
+md.settings.output_frequency = 1
+md.timestepping.final_time = 4
+
+smb = np.ones((md.mesh.numberofvertices,))*3.6
+smb = np.array([smb,smb*-1]).T
+
+md.smb.mass_balance = smb
+md.smb.mass_balance = np.concatenate((md.smb.mass_balance,[[1.5,3]]))
+md.transient.isthermal = 0
+
+#Dakota options
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#partitioning
+md.qmu.numberofpartitions = 20
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on')
+md.qmu.partition = md.qmu.partition - 1
+
+#variables
+md.qmu.variables.surface_mass_balance = normal_uncertain.normal_uncertain('scaled_SmbMassBalance',1,0.1)
+
+#responses
+md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.IceVolume = response_function.response_function('IceVolume',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux1 = response_function.response_function('indexed_MassFlux_1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux2 = response_function.response_function('indexed_MassFlux_2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux3 = response_function.response_function('indexed_MassFlux_3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux4 = response_function.response_function('indexed_MassFlux_4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux5 = response_function.response_function('indexed_MassFlux_5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.massFlux6 = response_function.response_function('indexed_MassFlux_6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+
+#mass flux profiles
+md.qmu.mass_flux_profiles = ['../Exp/MassFlux1.exp','../Exp/MassFlux2.exp','../Exp/MassFlux3.exp','../Exp/MassFlux4.exp','../Exp/MassFlux5.exp','../Exp/MassFlux6.exp']
+md.qmu.mass_flux_profile_directory = getcwd()
+
+##  nond_sampling study
+md.qmu.method = dakota_method.dakota_method('nond_samp')
+md.qmu.method = dmeth_params_set(md.qmu.method,'seed',1234,'samples',20,'sample_type','lhs')
+dver = str(version)
+if ((int(dver[0]) == 4 and int(dver[2])>2) or int(dver[0])>4):
+	md.qmu.method = dmeth_params_set(md.qmu.method,'rng','rnum2')
+
+#parameters
+md.qmu.params.direct = True
+md.qmu.params.analysis_components = ''
+md.qmu.params.interval_type = 'forward'
+md.qmu.params.tabular_graphics_data = True
+md.qmu.isdakota = 1
+
+if version >= 6:
+	md.qmu.params.analysis_driver = 'matlab'
+	md.qmu.params.evaluation_scheduling = 'master'
+	md.qmu.params.processors_per_evaluation = 2
+else:
+	md.qmu.params.analysis_driver = 'stressbalance'
+	md.qmu.params.evaluation_concurrency = 1
+
+md.stressbalance.reltol = 10**-5 #tighten for qmu analyses
+md.transient.requested_outputs = ['IceVolume']
+
+#solve
+md.verbose = verbose('000000000')	# this line is recommended
+md = solve(md,'Transient','overwrite','y')
+md.qmu.results = md.results.dakota
+
+#Fields and tolerances to track changes
+md.results.dakota.moments = []
+for i in range(8):
+	md.results.dakota.moments.append(md.results.dakota.dresp_out[i].mean)
+
+for i in range(8):
+	md.results.dakota.moments.append(md.results.dakota.dresp_out[i].stddev)
+
+field_names      = ['moments']
+field_tolerances = [1e-11]
+field_values = [md.results.dakota.moments]
Index: issm/trunk/test/NightlyRun/test235.py
===================================================================
--- issm/trunk/test/NightlyRun/test235.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test235.py	(revision 23394)
@@ -0,0 +1,97 @@
+#Test Name: SquareShelfTranForceNeg2dDakotaLocal
+import numpy as np
+import scipy.io as spio
+from os import getcwd
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from SMBgemb import *
+from partitioner import *
+from dmeth_params_set import *
+
+md = triangle(model(),'../Exp/Square.exp',180000.)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+
+md.timestepping.time_step = 1
+md.settings.output_frequency = 1
+md.timestepping.final_time = 4
+
+smb = np.ones((md.mesh.numberofvertices,))*3.6
+smb = np.array([smb,smb*-1]).T
+
+md.smb.mass_balance = smb
+md.smb.mass_balance = np.concatenate((md.smb.mass_balance,[[1.5,3]]))
+md.transient.isthermal = 0
+
+#Dakota options
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#partitioning
+md.qmu.numberofpartitions = 20
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on')
+md.qmu.partition = md.qmu.partition - 1
+
+#variables
+md.qmu.variables.surface_mass_balance = normal_uncertain.normal_uncertain('scaled_SmbMassBalance',1,100)
+
+#responses
+md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.IceVolume = response_function.response_function('IceVolume',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux1 = response_function.response_function('indexed_MassFlux_1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux2 = response_function.response_function('indexed_MassFlux_2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux3 = response_function.response_function('indexed_MassFlux_3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux4 = response_function.response_function('indexed_MassFlux_4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux5 = response_function.response_function('indexed_MassFlux_5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.massFlux6 = response_function.response_function('indexed_MassFlux_6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+
+#mass flux profiles
+md.qmu.mass_flux_profiles = ['../Exp/MassFlux1.exp','../Exp/MassFlux2.exp','../Exp/MassFlux3.exp','../Exp/MassFlux4.exp','../Exp/MassFlux5.exp','../Exp/MassFlux6.exp']
+md.qmu.mass_flux_profile_directory = getcwd()
+
+#method
+md.qmu.method = dakota_method.dakota_method('nond_l')
+
+#parameters
+md.qmu.params.direct = True
+md.qmu.params.analysis_components = ''
+md.qmu.params.interval_type = 'forward'
+md.qmu.params.fd_gradient_step_size = '0.1'
+md.qmu.isdakota = 1
+
+if version >= 6:
+	md.qmu.params.analysis_driver = 'matlab'
+	md.qmu.params.evaluation_scheduling = 'master'
+	md.qmu.params.processors_per_evaluation = 2
+else:
+	md.qmu.params.analysis_driver = 'stressbalance'
+	md.qmu.params.evaluation_concurrency = 1
+
+md.stressbalance.reltol = 10**-5 #tighten for qmu analyses
+md.transient.requested_outputs = ['IceVolume']
+
+#solve
+md.verbose = verbose('000000000')	# this line is recommended
+md = solve(md,'Transient','overwrite','y')
+md.qmu.results = md.results.dakota
+
+#Fields and tolerances to track changes
+md.results.dakota.moments = []
+for i in range(8):
+	md.results.dakota.moments.append(md.results.dakota.dresp_out[i].mean)
+
+for i in range(8):
+	md.results.dakota.moments.append(md.results.dakota.dresp_out[i].stddev)
+
+field_names      = ['moments']
+field_tolerances = [1e-11]
+field_values = [md.results.dakota.moments]
Index: issm/trunk/test/NightlyRun/test245.m
===================================================================
--- issm/trunk/test/NightlyRun/test245.m	(revision 23394)
+++ issm/trunk/test/NightlyRun/test245.m	(revision 23394)
@@ -0,0 +1,47 @@
+%Test Name: SquareShelfTranIspddSicopolisSSA2d
+md=triangle(model(),'../Exp/Square.exp',150000.);
+md=setmask(md,'all','');
+md=parameterize(md,'../Par/SquareShelf.par');
+
+% Use of SMBpddSicopolis
+md.smb = SMBpddSicopolis();
+% initalize pdd fields
+md.smb=initialize(md.smb,md);
+md.smb.s0p=md.geometry.surface;
+md.smb.s0t=md.geometry.surface;
+
+% 
+md.smb.monthlytemperatures=[]; md.smb.precipitation=[]; md.smb.precipitation=[];
+temp_ma_present=-10*ones(md.mesh.numberofvertices,1)-md.smb.rlaps*md.geometry.surface/1000;
+temp_mj_present=10*ones(md.mesh.numberofvertices,1)-md.smb.rlaps*md.geometry.surface/1000;
+precipitation=5*ones(md.mesh.numberofvertices,1);
+for imonth=0:11
+    md.smb.monthlytemperatures(1:md.mesh.numberofvertices,imonth+1)=md.materials.meltingpoint+temp_ma_present+(temp_mj_present-temp_ma_present)*sin(double(imonth+1-4)*pi/6.0);
+    md.smb.monthlytemperatures(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
+    md.smb.precipitation(1:md.mesh.numberofvertices,imonth+1)=precipitation;
+    md.smb.precipitation(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
+end
+
+% time steps and resolution
+md.timestepping.time_step=1;
+md.settings.output_frequency=1;
+md.timestepping.final_time=2;
+
+md.transient.issmb=1;
+md.transient.ismasstransport=1;
+md.transient.isstressbalance=0;
+md.transient.isthermal=0;
+
+md.transient.requested_outputs={'default','TemperaturePDD'};
+md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster
+md=solve(md,'Transient');
+
+%Fields and tolerances to track changes
+field_names     ={'TemperaturePDD1','SmbMassBalance1','TemperaturePDD2','SmbMassBalance2'};
+field_tolerances={1e-13,1e-13,1e-13,1e-13};
+field_values={...
+	(md.results.TransientSolution(1).TemperaturePDD),...
+	(md.results.TransientSolution(1).SmbMassBalance),...
+	(md.results.TransientSolution(2).TemperaturePDD),...
+	(md.results.TransientSolution(2).SmbMassBalance),...
+	};
Index: issm/trunk/test/NightlyRun/test245.py
===================================================================
--- issm/trunk/test/NightlyRun/test245.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test245.py	(revision 23394)
@@ -0,0 +1,59 @@
+#Test Name: SquareShelfTranIspddSicopolisSSA2d
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+
+from SMBpddSicopolis import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelf.py')
+
+# Use of SMBpddSicopolis
+md.smb  =  SMBpddSicopolis()
+# initalize pdd fields
+md.smb.initialize(md)
+md.smb.s0p = md.geometry.surface
+md.smb.s0t = md.geometry.surface
+
+# 
+md.smb.monthlytemperatures = np.empty((md.mesh.numberofvertices+1,12))
+md.smb.precipitation = np.empty((md.mesh.numberofvertices+1,12))
+temp_ma_present = -10. * np.ones((md.mesh.numberofvertices,)) - md.smb.rlaps * md.geometry.surface / 1000.
+temp_mj_present = 10. * np.ones((md.mesh.numberofvertices,)) - md.smb.rlaps * md.geometry.surface / 1000.
+precipitation = 5. * np.ones((md.mesh.numberofvertices,))
+for imonth in range(12): 
+	md.smb.monthlytemperatures[0:md.mesh.numberofvertices,imonth] = md.materials.meltingpoint + temp_ma_present + (temp_mj_present - temp_ma_present) * np.sin((imonth + 1. - 4.) * np.pi / 6.0)
+	md.smb.monthlytemperatures[md.mesh.numberofvertices,imonth] = ((imonth+1)/12.)
+	md.smb.precipitation[0:md.mesh.numberofvertices,imonth] = precipitation
+	md.smb.precipitation[md.mesh.numberofvertices,imonth] = ((imonth+1)/12.)
+
+# time steps and resolution
+md.timestepping.time_step = 1
+md.settings.output_frequency = 1
+md.timestepping.final_time = 2
+
+md.transient.issmb = 1
+md.transient.ismasstransport = 1
+md.transient.isstressbalance = 0
+md.transient.isthermal = 0
+
+md.transient.requested_outputs = ['default','TemperaturePDD']
+md.cluster = generic('name',gethostname(),'np',1) # 3 for the cluster
+md = solve(md,'Transient')
+
+#Fields and tolerances to track changes
+field_names      = ['TemperaturePDD1','SmbMassBalance1','TemperaturePDD2','SmbMassBalance2']
+field_tolerances = [1e-13,1e-13,1e-13,1e-13]
+field_values = [
+	md.results.TransientSolution[0].TemperaturePDD,
+	md.results.TransientSolution[0].SmbMassBalance,
+	md.results.TransientSolution[1].TemperaturePDD,
+	md.results.TransientSolution[1].SmbMassBalance,
+	]
+
Index: issm/trunk/test/NightlyRun/test290.m
===================================================================
--- issm/trunk/test/NightlyRun/test290.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test290.m	(revision 23394)
@@ -11,5 +11,5 @@
 %Fields and tolerances to track changes
 field_names     ={'Vx','Vy','Vz','Vel','Pressure'};
-field_tolerances={5e-5,5e-5,8e-5,5e-5,1e-7};
+field_tolerances={5e-5,5e-5,8e-4,5e-4,2e-7};
 field_values={...
 	(md.results.StressbalanceSolution.Vx),...
Index: issm/trunk/test/NightlyRun/test290.py
===================================================================
--- issm/trunk/test/NightlyRun/test290.py	(revision 23189)
+++ issm/trunk/test/NightlyRun/test290.py	(revision 23394)
@@ -21,5 +21,5 @@
 #Fields and tolerances to track changes
 field_names     =['Vx', 'Vy', 'Vz', 'Vel','Pressure']
-field_tolerances=[5e-5,5e-5,8e-5,5e-5,1e-7]
+field_tolerances=[5e-5,5e-5,8e-4,5e-4,2e-7]
 field_values=[\
 	md.results.StressbalanceSolution.Vx,\
Index: issm/trunk/test/NightlyRun/test3015.m
===================================================================
--- issm/trunk/test/NightlyRun/test3015.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test3015.m	(revision 23394)
@@ -15,12 +15,23 @@
 %setup autodiff parameters
 index=1; %this is the scalar component we are checking against
-md.autodiff.independents={...
-	independent('name','md.geometry.thickness','type','vertex','nods',md.mesh.numberofvertices,'fos_forward_index',index)
+
+if IssmConfig('_HAVE_CODIPACK_')
+	md.autodiff.independents={...
+		independent('name','md.geometry.thickness','type','vertex','nods',md.mesh.numberofvertices)
+	};
+	md.autodiff.dependents={...
+		dependent('name','IceVolume','type','scalar','fos_reverse_index',index)...
+		};
+	md.autodiff.driver='fos_reverse';
+else
+	md.autodiff.independents={...
+		independent('name','md.geometry.thickness','type','vertex','nods',md.mesh.numberofvertices,'fos_forward_index',index)
 	};
 
-md.autodiff.dependents={...
-	dependent('name','IceVolume','type','scalar')...
-	};
-md.autodiff.driver='fos_forward';
+	md.autodiff.dependents={...
+		dependent('name','IceVolume','type','scalar')...
+		};
+	md.autodiff.driver='fos_forward';
+end
 
 %parameters for the step-wise derivative
@@ -70,7 +81,11 @@
 md=solve(md,'Masstransport');
 %retrieve directly
-dVdh_ad=md.results.MasstransportSolution.AutodiffJacobian;
+if IssmConfig('_HAVE_CODIPACK_')
+	dVdh_ad=md.results.MasstransportSolution.AutodiffJacobian(index);
+else
+	dVdh_ad=md.results.MasstransportSolution.AutodiffJacobian;
+end
 
-disp(sprintf('dV/dh: analytical:  %16.16g\n       using adolc:  %16.16g\n',dVdh_an,dVdh_ad));
+disp(sprintf('dV/dh: analytical:  %16.16g\n       using ad:    %16.16g\n',dVdh_an,dVdh_ad));
 
 %Fields and tolerances to track changes
Index: issm/trunk/test/NightlyRun/test3480.m
===================================================================
--- issm/trunk/test/NightlyRun/test3480.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test3480.m	(revision 23394)
@@ -18,4 +18,5 @@
 md.inversion.iscontrol=1;
 md.autodiff.isautodiff=1;
+md.autodiff.driver='fos_reverse';
 
 md.friction.coefficient(1:md.mesh.numberofvertices,1)=50;
Index: issm/trunk/test/NightlyRun/test3481.m
===================================================================
--- issm/trunk/test/NightlyRun/test3481.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test3481.m	(revision 23394)
@@ -18,4 +18,5 @@
 md.inversion.iscontrol=1;
 md.autodiff.isautodiff=1;
+md.autodiff.driver='fos_reverse';
 
 md.friction.coefficient(1:md.mesh.numberofvertices,1)=50;
Index: issm/trunk/test/NightlyRun/test413.py
===================================================================
--- issm/trunk/test/NightlyRun/test413.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test413.py	(revision 23394)
@@ -0,0 +1,66 @@
+#Test Name: SquareSheetShelfDiadSSA3dDakotaPart
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from partitioner import *
+from importancefactors import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+
+#Dakota options
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#partitioning
+md.qmu.numberofpartitions = 20
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on')
+md.qmu.partition = md.qmu.partition-1
+
+#variables
+md.qmu.variables.rho_ice = normal_uncertain.normal_uncertain('MaterialsRhoIce',md.materials.rho_ice,0.01)
+md.qmu.variables.drag_coefficient = normal_uncertain.normal_uncertain('scaled_FrictionCoefficient',1,0.01)
+
+#responses
+md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+
+#method
+md.qmu.method = dakota_method.dakota_method('nond_l')
+
+#parameters
+md.qmu.params.direct = True
+md.qmu.params.interval_type = 'forward'
+
+if version >= 6:
+	md.qmu.params.analysis_driver = 'matlab'
+	md.qmu.params.evaluation_scheduling = 'master'
+	md.qmu.params.processors_per_evaluation = 2
+else:
+	md.qmu.params.analysis_driver = 'stressbalance'
+	md.qmu.params.evaluation_concurrency = 1
+
+
+#imperative!
+md.stressbalance.reltol = 10**-5 #tighten for qmu analyses
+md.qmu.isdakota = 1
+
+#solve
+md.verbose = verbose('000000000')	# this line is recommended
+md = solve(md,'Stressbalance','overwrite','y')
+
+#Fields and tolerances to track changes
+md.qmu.results = md.results.dakota
+md.results.dakota.importancefactors = importancefactors(md,'scaled_FrictionCoefficient','MaxVel').T
+field_names = ['importancefactors']
+field_tolerances = [1e-10]
+field_values = [md.results.dakota.importancefactors]
Index: issm/trunk/test/NightlyRun/test414.py
===================================================================
--- issm/trunk/test/NightlyRun/test414.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test414.py	(revision 23394)
@@ -0,0 +1,92 @@
+#Test Name: SquareSheetShelfDiadSSA3dDakotaMassFlux
+import numpy as np
+from os import getcwd
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from partitioner import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.materials.rho_ice = 10**7 #involved in the mass flux, make it easy
+md.geometry.thickness[:] = 1 #make it easy
+md.geometry.surface = md.geometry.base+md.geometry.thickness
+
+#constrain all velocities to 1 m/yr, in the y-direction
+md.stressbalance.spcvx = np.zeros((md.mesh.numberofvertices,))
+md.stressbalance.spcvy = np.ones((md.mesh.numberofvertices,))
+md.stressbalance.spcvz = np.zeros((md.mesh.numberofvertices,))
+
+#Dakota options
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#partitioning
+md.qmu.numberofpartitions = 20
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on')
+md.qmu.partition = md.qmu.partition-1
+
+#variables
+md.qmu.variables.drag_coefficient = normal_uncertain.normal_uncertain('scaled_FrictionCoefficient',1,0.01)
+
+#responses
+md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux1 = response_function.response_function('indexed_MassFlux_1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux2 = response_function.response_function('indexed_MassFlux_2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux3 = response_function.response_function('indexed_MassFlux_3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux4 = response_function.response_function('indexed_MassFlux_4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux5 = response_function.response_function('indexed_MassFlux_5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux6 = response_function.response_function('indexed_MassFlux_6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux7 = response_function.response_function('indexed_MassFlux_7',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+
+#mass flux profiles
+md.qmu.mass_flux_profiles = ['../Exp/MassFlux1.exp','../Exp/MassFlux2.exp','../Exp/MassFlux3.exp','../Exp/MassFlux4.exp','../Exp/MassFlux5.exp','../Exp/MassFlux6.exp','../Exp/Square.exp']
+md.qmu.mass_flux_profile_directory = getcwd()
+
+#method
+md.qmu.method = dakota_method.dakota_method('nond_l')
+
+#parameters
+md.qmu.params.direct = True
+md.qmu.params.interval_type = 'forward'
+md.qmu.isdakota = 1
+md.stressbalance.reltol = 10**-5 #tighten for qmu analyses
+
+if version >= 6:
+	md.qmu.params.analysis_driver = 'matlab'
+	md.qmu.params.evaluation_scheduling = 'master'
+	md.qmu.params.processors_per_evaluation = 2
+else:
+	md.qmu.params.analysis_driver = 'stressbalance'
+	md.qmu.params.evaluation_concurrency = 1
+
+#solve
+md.verbose = verbose('000000000')	# this line is recommended
+md = solve(md,'Stressbalance','overwrite','y')
+md.qmu.results = md.results.dakota
+
+#Fields and tolerances to track changes
+#ok, mass flux of 3 profiles should be -3 Gt/yr -3 Gt/yr and the sum, which is -6 Gt/yr
+#we recover those mass fluxes through the mean of the response.
+#also, we recover the max velo, which should be 1m/yr. 
+#we put all that data in the moments, which we will use to test for success.
+#also, check that the stddev are 0.
+md.results.dakota.moments = []
+for i in range(8):
+	md.results.dakota.moments.append(md.results.dakota.dresp_out[i].mean)
+
+for i in range(8):
+	md.results.dakota.moments.append(md.results.dakota.dresp_out[i].stddev)
+
+field_names      = ['moments']
+field_tolerances = [1e-11]
+field_values = [md.results.dakota.moments]
Index: issm/trunk/test/NightlyRun/test417.py
===================================================================
--- issm/trunk/test/NightlyRun/test417.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test417.py	(revision 23394)
@@ -0,0 +1,106 @@
+#Test Name: SquareSheetShelfDiadSSA3dDakotaSamp
+import numpy as np
+from os import getcwd
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from partitioner import *
+from dmeth_params_set import *
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+md.materials.rho_ice = 10**7 #involved in the mass flux, make it easy
+md.geometry.thickness[:] = 1 #make it easy
+md.geometry.surface = md.geometry.base+md.geometry.thickness
+
+#constrain all velocities to 1 m/yr, in the y-direction
+md.stressbalance.spcvx[:] = 0
+md.stressbalance.spcvy[:] = 1
+md.stressbalance.spcvz[:] = 0
+
+#Dakota options
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#partitioning
+md.qmu.numberofpartitions = 20
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on')
+md.qmu.partition = md.qmu.partition-1
+
+#variables
+md.qmu.variables.drag_coefficient = normal_uncertain.normal_uncertain('scaled_FrictionCoefficient',1,0.01)
+
+#responses
+md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux1 = response_function.response_function('indexed_MassFlux_1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux2 = response_function.response_function('indexed_MassFlux_2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux3 = response_function.response_function('indexed_MassFlux_3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux4 = response_function.response_function('indexed_MassFlux_4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux5 = response_function.response_function('indexed_MassFlux_5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux6 = response_function.response_function('indexed_MassFlux_6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.MassFlux7 = response_function.response_function('indexed_MassFlux_7',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+
+#mass flux profiles
+md.qmu.mass_flux_profiles = ['../Exp/MassFlux1.exp','../Exp/MassFlux2.exp','../Exp/MassFlux3.exp','../Exp/MassFlux4.exp','../Exp/MassFlux5.exp','../Exp/MassFlux6.exp','../Exp/Square.exp']
+md.qmu.mass_flux_profile_directory = getcwd()
+
+#method
+md.qmu.method = dakota_method.dakota_method('nond_samp')
+md.qmu.method = dmeth_params_set(md.qmu.method,'seed',1234,'samples',20,'sample_type','lhs')
+
+#parameters
+md.qmu.params.direct = True
+md.qmu.params.analysis_components = ''
+md.qmu.params.tabular_graphics_data = True;
+
+if version >= 6:
+	md.qmu.params.analysis_driver = 'matlab'
+	md.qmu.params.evaluation_scheduling = 'master'
+	md.qmu.params.processors_per_evaluation = 2
+else:
+	md.qmu.params.analysis_driver = 'stressbalance'
+	md.qmu.params.evaluation_concurrency = 1
+
+#partitioning
+md.qmu.numberofpartitions = 20
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on')
+md.qmu.partition = md.qmu.partition-1
+md.qmu.isdakota = 1
+
+md.stressbalance.reltol = 10**-5	#tighten for qmu analyses
+
+#solve
+md.verbose = verbose('000000000')	# this line is recommended
+
+# There may be a pair of numpy warnings in the function true_divide,
+#	this is normal and will not affect the results
+#	See src/m/qmu/dakota_out_parse.py, function "dak_tab_out" for details
+md = solve(md,'Stressbalance','overwrite','y')
+
+#Fields and tolerances to track changes
+md.qmu.results = md.results.dakota
+
+#ok, mass flux of 3 profiles should be -3 Gt/yr -3 Gt/yr and the sum, which is -6 Gt/yr
+#we recover those mass fluxes through the mean of the response.
+#also, we recover the max velo, which should be 1m/yr. 
+#we put all that data in the montecarlo field, which we will use to test for success.
+#also, check that the stddev are 0.
+md.results.dakota.montecarlo = []
+for i in range(8):
+	md.results.dakota.montecarlo.append(md.results.dakota.dresp_out[i].mean)
+
+for i in range(8):
+	md.results.dakota.montecarlo.append(md.results.dakota.dresp_out[i].stddev)
+
+field_names      = ['moments']
+field_tolerances = [1e-11]
+field_values     = [md.results.dakota.montecarlo]
Index: issm/trunk/test/NightlyRun/test418.py
===================================================================
--- issm/trunk/test/NightlyRun/test418.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test418.py	(revision 23394)
@@ -0,0 +1,43 @@
+#Test Name: SquareSheetShelfDiadSSA3dDakotaAreaAverage
+
+# this test may crash
+
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from partitioner import *
+from AreaAverageOntoPartition import *
+
+#test partitioning, and partition averaging
+# python cannot handle resolutions greater than 30010
+md = triangle(model(),'../Exp/Square.exp',30000.)
+#print md.mesh.numberofvertices
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+
+#partitioning
+md.qmu.numberofpartitions = 100
+
+#corrupted size vs. prev_size
+#Aborted (core dumped)
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions)
+
+
+md.qmu.partition=md.qmu.partition-1
+
+vector = np.arange(1,1+md.mesh.numberofvertices,1).reshape(-1,1)
+# double check this before committing:
+#print 'before AreaAverageOntoPartition'
+vector_on_partition = AreaAverageOntoPartition(md,vector)
+vector_on_nodes = vector_on_partition[md.qmu.partition+1]
+
+field_names      = ['vector_on_nodes']
+field_tolerances = [1e-11]
+field_values     = [vector_on_nodes]
Index: issm/trunk/test/NightlyRun/test420.py
===================================================================
--- issm/trunk/test/NightlyRun/test420.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test420.py	(revision 23394)
@@ -0,0 +1,71 @@
+#Test Name: SquareSheetShelfDakotaScaledResponse
+
+import numpy as np
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from partitioner import *
+
+md = triangle(model(),'../Exp/Square.exp',200000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+
+#partitioning
+md.qmu.numberofpartitions = 10
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions)
+md.qmu.partition = md.qmu.partition-1
+md.qmu.isdakota = 1
+
+#Dakota options
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#variables
+md.qmu.variables.rho_ice = normal_uncertain.normal_uncertain('MaterialsRhoIce',md.materials.rho_ice,0.01)
+
+#responses
+md.qmu.responses.MaxVel = response_function.response_function('scaled_Thickness',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+
+#method
+md.qmu.method = dakota_method.dakota_method('nond_l')
+
+#parameters
+md.qmu.params.direct = True
+md.qmu.params.interval_type = 'forward'
+
+if version >= 6:
+	md.qmu.params.analysis_driver = 'matlab'
+	md.qmu.params.evaluation_scheduling = 'master'
+	md.qmu.params.processors_per_evaluation = 2
+else:
+	md.qmu.params.analysis_driver = 'stressbalance'
+	md.qmu.params.evaluation_concurrency = 1
+
+#imperative! 
+md.stressbalance.reltol = 10**-5 #tighten for qmu analysese
+
+#solve
+md.verbose = verbose('000000000')	# this line is recommended
+md = solve(md,'Stressbalance','overwrite','y')
+md.qmu.results = md.results.dakota
+
+#test on thickness
+h = np.zeros((md.qmu.numberofpartitions,))
+for i in range(md.qmu.numberofpartitions):
+	h[i] = md.qmu.results.dresp_out[i].mean
+
+#project onto grid
+thickness = h[(md.qmu.partition).flatten()]
+
+#Fields and tolerances to track changes
+field_names      = ['Thickness']
+field_tolerances = [1e-10]
+field_values     = [thickness]
Index: issm/trunk/test/NightlyRun/test421.m
===================================================================
--- issm/trunk/test/NightlyRun/test421.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test421.m	(revision 23394)
@@ -10,5 +10,5 @@
 %Fields and tolerances to track changes
 field_names     ={'Vx','Vy','Vz','Vel','Pressure'};
-field_tolerances={1e-06,9e-07,2e-05,2e-06,2e-06};
+field_tolerances={2e-06,9e-07,2e-05,2e-06,2e-06};
 field_values={...
 	(md.results.StressbalanceSolution.Vx),...
Index: issm/trunk/test/NightlyRun/test421.py
===================================================================
--- issm/trunk/test/NightlyRun/test421.py	(revision 23189)
+++ issm/trunk/test/NightlyRun/test421.py	(revision 23394)
@@ -20,5 +20,5 @@
 #Fields and tolerances to track changes
 field_names     =['Vx','Vy','Vz','Vel','Pressure']
-field_tolerances=[1e-06,9e-07,2e-05,2e-06,2e-06]
+field_tolerances=[2e-06,9e-07,2e-05,2e-06,2e-06]
 field_values=[\
 	md.results.StressbalanceSolution.Vx,\
Index: issm/trunk/test/NightlyRun/test444.py
===================================================================
--- issm/trunk/test/NightlyRun/test444.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test444.py	(revision 23394)
@@ -0,0 +1,127 @@
+#Test Name: SquareShelfTranForceNeg2dDakotaLocal
+import numpy as np
+from os import getcwd
+from model import *
+from socket import gethostname
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from solve import *
+from SMBgemb import *
+from partitioner import *
+from dmeth_params_set import *
+from ContourToMesh import *
+from regionaloutput import *
+
+#model not consistent:  equality thickness=surface-base violated
+
+md = triangle(model(),'../Exp/Square.exp',150000.)
+md = setmask(md,'../Exp/SquareShelf.exp','')
+md = parameterize(md,'../Par/SquareSheetShelf.py')
+md.geometry.bed = md.geometry.base.copy()
+pos = np.where(md.mask.groundedice_levelset < 0)
+md.geometry.bed[pos] = md.geometry.base[pos]-10
+md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices,))
+md.friction.p = np.ones((md.mesh.numberofelements,))
+md.friction.q = np.ones((md.mesh.numberofelements,))
+md.transient.isthermal = 0
+md.transient.isgroundingline = 1
+md.groundingline.migration = 'AggressiveMigration'
+
+md.settings.output_frequency = 3
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',oshostname(),'np',3)
+
+regionalmask = np.zeros((md.mesh.numberofvertices,))
+c_in = ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,'../Exp/SquareHalfRight.exp','node',1)
+regionalmask[np.where(c_in)] = 1
+md.transient.requested_outputs = ['default','GroundedArea','FloatingArea','IceVolumeAboveFloatation','GroundedArea1','FloatingArea1','TotalFloatingBmb1','TotalGroundedBmb1','TotalSmb1',
+	'IceMass1','IceVolume1','IceVolumeAboveFloatation1','IceVolumeAboveFloatation']
+md.outputdefinition.definitions.append(regionaloutput('name','GroundedArea1','outputnamestring','GroundedArea','mask',regionalmask,
+	'definitionstring','Outputdefinition1'))
+md.outputdefinition.definitions.append(regionaloutput('name','FloatingArea1','outputnamestring','FloatingArea','mask',regionalmask,
+	'definitionstring','Outputdefinition2'))
+md.outputdefinition.definitions.append(regionaloutput('name','TotalFloatingBmb1','outputnamestring','TotalFloatingBmb','mask',regionalmask,
+	'definitionstring','Outputdefinition3'))
+md.outputdefinition.definitions.append(regionaloutput('name','TotalGroundedBmb1','outputnamestring','TotalGroundedBmb','mask',regionalmask,
+	'definitionstring','Outputdefinition4'))
+md.outputdefinition.definitions.append(regionaloutput('name','IceMass1','outputnamestring','IceMass','mask',regionalmask,
+	'definitionstring','Outputdefinition5'))
+md.outputdefinition.definitions.append(regionaloutput('name','IceVolume1','outputnamestring','IceVolume','mask',regionalmask,
+	'definitionstring','Outputdefinition6'))
+md.outputdefinition.definitions.append(regionaloutput('name','IceVolumeAboveFloatation1','outputnamestring','IceVolumeAboveFloatation','mask',regionalmask,
+	'definitionstring','Outputdefinition7'))
+md.outputdefinition.definitions.append(regionaloutput('name','TotalSmb1','outputnamestring','TotalSmb','mask',regionalmask,
+	'definitionstring','Outputdefinition8'))
+md.outputdefinition.definitions.append(regionaloutput('name','TotalSmb2','outputnamestring','TotalSmb','mask',regionalmask,
+	 'definitionstring','Outputdefinition9'))
+
+md.extrude(3,1.)
+md.collapse()
+
+#Dakota options
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#variables
+md.qmu.variables.drag_coefficient = normal_uncertain.normal_uncertain('scaled_BasalforcingsFloatingiceMeltingRate',1,0.1)
+
+#responses
+md.qmu.responses.IceMass1 = response_function.response_function('Outputdefinition5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.IceVolume1 = response_function.response_function('Outputdefinition6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.IceVolumeAboveFloatation1 = response_function.response_function('Outputdefinition7',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.IceVolumeAboveFloatation = response_function.response_function('IceVolumeAboveFloatation',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.GroundedArea1 = response_function.response_function('Outputdefinition1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.FloatingArea1 = response_function.response_function('Outputdefinition2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.TotalFloatingBmb1 = response_function.response_function('Outputdefinition3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.TotalGroundedBmb1 = response_function.response_function('Outputdefinition4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.TotalSmb1 = response_function.response_function('Outputdefinition8',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.TotalSmb2 = response_function.response_function('Outputdefinition9',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.FloatingArea = response_function.response_function('FloatingArea',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+
+#method
+md.qmu.method = dakota_method.dakota_method('nond_samp')
+md.qmu.method = dmeth_params_set(md.qmu.method,'seed',1234,'samples',20,'sample_type','random')
+
+#parameters
+md.qmu.params.direct = True
+md.qmu.params.analysis_components = ''
+md.qmu.params.tabular_graphics_data = True
+
+if version >= 6:
+	md.qmu.params.analysis_driver = 'matlab'
+	md.qmu.params.evaluation_scheduling = 'master'
+	md.qmu.params.processors_per_evaluation = 2
+else:
+	md.qmu.params.analysis_driver = 'stressbalance'
+	md.qmu.params.evaluation_concurrency = 1
+
+#partitioning
+md.qmu.numberofpartitions = 10
+md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on')
+md.qmu.partition = md.qmu.partition-1
+md.qmu.isdakota = 1
+
+md.stressbalance.reltol = 10**-5 #tighten for qmu analyses
+
+#solve
+md.verbose = verbose('000000000')	# this line is recommended
+md = solve(md,'Transient','overwrite','y')
+
+#Fields and tolerances to track changes
+md.qmu.results = md.results.dakota
+
+#we put all the mean and stdev data in the montecarlo field, which we will use to test for success.
+md.results.dakota.montecarlo = []
+for i in range(11):
+	md.results.dakota.montecarlo.append(md.results.dakota.dresp_out[i].mean)
+
+for i in range(11):
+	md.results.dakota.montecarlo.append(md.results.dakota.dresp_out[i].stddev)
+
+field_names      = ['montecarlo']
+field_tolerances = [1e-11]
+field_values = [md.results.dakota.montecarlo]
Index: issm/trunk/test/NightlyRun/test455.m
===================================================================
--- issm/trunk/test/NightlyRun/test455.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test455.m	(revision 23394)
@@ -14,5 +14,5 @@
 	md=solve(md,'Stressbalance');
 	field_names     ={field_names{:},['Vx' i{1}],['Vy' i{1}],['Vz' i{1}],['Vel' i{1}],['Pressure' i{1}]};
-	field_tolerances={field_tolerances{:},6e-08,6e-08,6e-08,6e-08,3e-13};
+	field_tolerances={field_tolerances{:},7e-08,6e-08,6e-08,6e-08,3e-13};
 	field_values={field_values{:},...
 	(md.results.StressbalanceSolution.Vx),...
Index: issm/trunk/test/NightlyRun/test455.py
===================================================================
--- issm/trunk/test/NightlyRun/test455.py	(revision 23189)
+++ issm/trunk/test/NightlyRun/test455.py	(revision 23394)
@@ -24,5 +24,5 @@
 	md=solve(md,'Stressbalance')
 	field_names     =field_names+['Vx'+i,'Vy'+i,'Vz'+i,'Vel'+i,'Pressure'+i]
-	field_tolerances=field_tolerances+[6e-08,6e-08,6e-08,6e-08,3e-13]
+	field_tolerances=field_tolerances+[7e-08,6e-08,6e-08,6e-08,3e-13]
 	field_values=field_values+[\
 			md.results.StressbalanceSolution.Vx,\
Index: issm/trunk/test/NightlyRun/test701.m
===================================================================
--- issm/trunk/test/NightlyRun/test701.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test701.m	(revision 23394)
@@ -35,4 +35,5 @@
 md.stressbalance.spcvx(find(vertexflags(md.mesh,4)))=0;
 md.stressbalance.spcvy(find(vertexflags(md.mesh,4)))=0;
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
 
 %Misc
Index: issm/trunk/test/NightlyRun/test701.py
===================================================================
--- issm/trunk/test/NightlyRun/test701.py	(revision 23189)
+++ issm/trunk/test/NightlyRun/test701.py	(revision 23394)
@@ -42,4 +42,5 @@
 md.stressbalance.spcvx[np.where(md.mesh.vertexflags(4))] = 0.
 md.stressbalance.spcvy[np.where(md.mesh.vertexflags(4))] = 0.
+md.basalforcings.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices,))
 
 #Misc
Index: issm/trunk/test/NightlyRun/test702.m
===================================================================
--- issm/trunk/test/NightlyRun/test702.m	(revision 23189)
+++ issm/trunk/test/NightlyRun/test702.m	(revision 23394)
@@ -34,4 +34,5 @@
 md.stressbalance.spcvx(find(vertexflags(md.mesh,4)))=800;
 md.stressbalance.spcvy(find(vertexflags(md.mesh,4)))=0;
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
 
 %Misc
Index: issm/trunk/test/NightlyRun/test702.py
===================================================================
--- issm/trunk/test/NightlyRun/test702.py	(revision 23394)
+++ issm/trunk/test/NightlyRun/test702.py	(revision 23394)
@@ -0,0 +1,85 @@
+#Test Name: FlowbandFSsheetshelf
+import numpy as np
+#from scipy.interpolate import interp1d
+from model import  * 
+from setflowequation import  * 
+from solve import  *
+from NowickiProfile import *
+from bamgflowband import *
+from paterson import *
+
+#mesh parameters
+x = np.arange(-5,5.5,.5).T
+[b,h,sea] = NowickiProfile(x)
+x = x * 10**3
+h = h * 10**3
+b = (b-sea) * 10**3
+
+#mesh domain
+md = bamgflowband(model(),x,b+h,b,'hmax',150.)
+
+#parameterize
+md.geometry.surface = np.interp(md.mesh.x,x,b+h)
+md.geometry.base = np.interp(md.mesh.x,x,b)
+md.geometry.thickness = md.geometry.surface - md.geometry.base
+
+#mask
+md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices,))
+md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 0
+md.mask.groundedice_levelset = -0.5*np.ones((md.mesh.numberofvertices))
+md.mask.groundedice_levelset[np.where(md.mesh.x<0)]=0.5
+
+#materials
+md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices,))
+md.materials.rheology_B = paterson(md.initialization.temperature)
+md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements,))
+
+#damage
+md.damage.D = np.zeros((md.mesh.numberofvertices,))
+md.damage.spcdamage = float('NaN') * np.ones((md.mesh.numberofvertices,))
+
+#friciton
+md.friction.coefficient = np.zeros((md.mesh.numberofvertices,))
+md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 20
+md.friction.p = np.ones((md.mesh.numberofelements,))
+md.friction.q = np.ones((md.mesh.numberofelements,))
+
+#boundary conditions
+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.referential = float('NaN') * np.ones((md.mesh.numberofvertices,6))
+md.stressbalance.loadingforce = np.zeros((md.mesh.numberofvertices,3))
+md.stressbalance.spcvx[np.where(md.mesh.vertexflags(4))] = 800.
+md.stressbalance.spcvy[np.where(md.mesh.vertexflags(4))] = 0.
+md.basalforcings.floatingice_melting_rate=np.zeros((md.mesh.numberofvertices,))
+
+#misc
+md = setflowequation(md,'FS','all')
+md.stressbalance.abstol = float('NaN')
+md.stressbalance.FSreconditioning = 1
+md.stressbalance.maxiter = 20
+md.flowequation.augmented_lagrangian_r = 10000.
+md.flowequation.augmented_lagrangian_rhop = 10000.
+md.initialization.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.y)
+md.miscellaneous.name = 'test702'
+md.groundingline.migration = 'None'
+md.cluster = generic('np',2)
+
+#Fields and tolerances to track changes
+field_names = []
+field_tolerances = []
+field_values = []
+for i in ['MINI','MINIcondensed','TaylorHood','XTaylorHood','LATaylorHood']:
+	print ' '
+	print '======Testing ' +i+ ' Full-Stokes Finite element====='
+	md.flowequation.fe_FS = i
+	md = solve(md,'Stressbalance')
+	field_names.extend(['Vx'+i,'Vy'+i,'Vel'+i,'Pressure'+i])
+	field_tolerances.extend([8e-5,8e-5,8e-5,1e-08])
+	field_values.extend([
+		md.results.StressbalanceSolution.Vx,
+		md.results.StressbalanceSolution.Vy,
+		md.results.StressbalanceSolution.Vel,
+		md.results.StressbalanceSolution.Pressure
+		])
