Index: /issm/trunk-jpl/test/NightlyRun/runme.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/runme.py	(revision 23129)
+++ /issm/trunk-jpl/test/NightlyRun/runme.py	(revision 23130)
@@ -183,5 +183,11 @@
 						if str(archive) == 'None':
 							raise NameError("Field name '"+archive_name+'_field'+str(k+1)+"' does not exist in archive file.")
+						if np.shape(field) != np.shape(archive) and not np.shape(field) in [(1,1),(0,0),(1,0),(0,1)]:
+							field = field.T
+							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)
+
                                                 if not np.isscalar(error_diff): error_diff=error_diff[0]
 
Index: /issm/trunk-jpl/test/NightlyRun/test218.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test218.py	(revision 23130)
+++ /issm/trunk-jpl/test/NightlyRun/test218.py	(revision 23130)
@@ -0,0 +1,120 @@
+#Test Name: SquareShelfConstrainedDakotaB
+import numpy as np
+from model import *
+from socket import gethostname
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from paterson import *
+from solve import *
+from generic import generic
+from squaremesh import *
+from partitioner import *
+from IssmConfig import *
+from importancefactors import *
+
+from normal_uncertain import *
+from response_function import *
+from dakota_method import *
+
+from helpers import *
+
+from verbose import *
+
+md = squaremesh(model(),1000000,1000000,5,5)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelf2.py')
+md = setflowequation(md,'SSA','all')
+md.cluster = generic('name',gethostname(),'np',3)
+
+#redo the parameter file for this special shelf. 
+#constant thickness, constrained (vy=0) flow into an icefront, 
+#from 0 m/yr at the grounding line.
+
+#needed later
+ymin = min(md.mesh.y)
+ymax = max(md.mesh.y)
+xmin = min(md.mesh.x)
+xmax = max(md.mesh.x)
+
+di = md.materials.rho_ice / md.materials.rho_water
+
+h = 1000.
+md.geometry.thickness = h * np.ones((md.mesh.numberofvertices,))
+md.geometry.base = -md.materials.rho_ice / md.materials.rho_water * md.geometry.thickness
+md.geometry.surface = md.geometry.base + md.geometry.thickness
+
+#Initial velocity and pressure
+md.initialization.vx = np.zeros((md.mesh.numberofvertices,))
+md.initialization.vy = np.zeros((md.mesh.numberofvertices,))
+md.initialization.vz = np.zeros((md.mesh.numberofvertices,))
+md.initialization.pressure = np.zeros((md.mesh.numberofvertices,))
+
+#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,))
+
+#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,))
+
+#constrain flanks to 0 normal velocity
+pos = np.where((md.mesh.x == xmin) | (md.mesh.x == xmax))
+md.stressbalance.spcvx[pos] = 0.
+md.stressbalance.spcvz[pos] = float('Nan')
+
+#constrain grounding line to 0 velocity
+pos = np.where(md.mesh.y ==ymin)
+md.stressbalance.spcvx[pos] = 0.
+md.stressbalance.spcvy[pos] = 0.
+
+#partitioning
+md.qmu.numberofpartitions = md.mesh.numberofvertices
+md = partitioner(md,'package','linear','npart',md.qmu.numberofpartitions)
+md.qmu.partition = (md.qmu.partition - 1.).reshape(-1,1).T
+
+#Dakota options
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+# returns tuple "(u'6.2',)" -> unicode string '6.2', convert to float
+version=float(version[0])
+
+#variables
+md.qmu.variables.rheology_B = normal_uncertain.normal_uncertain('scaled_MaterialsRheologyB',1,.05)
+
+#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**-10 #tighten for qmu analysis
+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_MaterialsRheologyB','MaxVel').reshape(-1,1)
+field_names = ['importancefactors']
+field_tolerances = [1e-10]
+field_values=[md.results.dakota.importancefactors]
Index: /issm/trunk-jpl/test/NightlyRun/test244.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test244.py	(revision 23130)
+++ /issm/trunk-jpl/test/NightlyRun/test244.py	(revision 23130)
@@ -0,0 +1,137 @@
+#Test Name: SquareShelfSMBGembDakota
+import numpy as np
+import scipy.io as spio
+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 IssmConfig import *
+
+from partitioner import *
+from dakota_method import *
+from normal_uncertain import *
+from uniform_uncertain import *
+from response_function import *
+from dmeth_params_set import *
+
+md = triangle(model(),'../Exp/Square.exp',200000.)
+md = setmask(md,'all','')
+md = parameterize(md,'../Par/SquareShelf.py')
+md = setflowequation(md,'SSA','all')
+md.materials.rho_ice = 910
+md.cluster = generic('name',gethostname(),'np',3)
+md.geometry.bed = md.geometry.base
+
+# Use of Gemb method for SMB computation
+md.smb = SMBgemb()
+md.smb.setdefaultparameters(md.mesh,md.geometry)
+
+#load hourly surface forcing date from 1979 to 2009:
+inputs = spio.loadmat('../Data/gemb_input.mat',squeeze_me = True)
+
+#setup the inputs: 
+md.smb.Ta = np.append(np.tile(np.conjugate(inputs['Ta0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.V = np.append(np.tile(np.conjugate(inputs['V0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.dswrf = np.append(np.tile(np.conjugate(inputs['dsw0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.dlwrf = np.append(np.tile(np.conjugate(inputs['dlw0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.P = np.append(np.tile(np.conjugate(inputs['P0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.eAir = np.append(np.tile(np.conjugate(inputs['eAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.pAir = np.append(np.tile(np.conjugate(inputs['pAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.pAir = np.append(np.tile(np.conjugate(inputs['pAir0']),(md.mesh.numberofelements,1)),np.conjugate([inputs['dateN']]),axis=0)
+md.smb.Vz = np.tile(np.conjugate(inputs['LP']['Vz']),(md.mesh.numberofelements,1)).flatten()
+md.smb.Tz = np.tile(np.conjugate(inputs['LP']['Tz']),(md.mesh.numberofelements,1)).flatten()
+md.smb.Tmean = np.tile(np.conjugate(inputs['LP']['Tmean']),(md.mesh.numberofelements,1)).flatten()
+md.smb.C = np.tile(np.conjugate(inputs['LP']['C']),(md.mesh.numberofelements,1)).flatten()
+
+#smb settings
+md.smb.requested_outputs = ['SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance']
+
+#only run smb core: 
+md.transient.isstressbalance = 0
+md.transient.ismasstransport = 1
+md.transient.isthermal = 0
+
+#time stepping: 
+md.timestepping.start_time = 1965.
+md.timestepping.final_time = 1965.75
+md.timestepping.time_step = 1./365.0
+md.timestepping.interp_forcings = 0.
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#partitioning
+md.qmu.numberofpartitions = md.mesh.numberofelements
+md = partitioner(md,'package','linear')
+md.qmu.partition = (md.qmu.partition - 1).reshape(-1,1).T
+
+#variables
+md.qmu.variables.surface_mass_balanceC = normal_uncertain.normal_uncertain('scaled_SmbC',1,0.5)
+Tmin = 273.
+telms = np.atleast_2d(np.min(md.smb.Ta[0:-1,:],1))
+mint_on_partition = telms.flatten()
+for pa in range(np.size(mint_on_partition)):
+	vi = np.where(md.qmu.partition == pa)
+	mint = telms[vi]*1.05
+	pos = np.where(mint < Tmin)
+	mint[pos] = Tmin
+	mint_on_partition[pa] = max(mint/telms[vi])
+
+mint_on_partition[np.where(np.isnan(mint_on_partition))] = 10**-10
+md.qmu.variables.surface_mass_balanceTa = uniform_uncertain.uniform_uncertain('scaled_SmbTa',1,0.05)
+md.qmu.variables.surface_mass_balanceTa[0].lower = 0.95
+md.qmu.variables.surface_mass_balanceTa[0].upper = np.maximum(np.minimum(np.maximum(1.05,mint_on_partition),0.9999),0.0001)
+
+#responses
+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.IceMass = response_function.response_function('IceMass',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+md.qmu.responses.TotalSmb = response_function.response_function('TotalSmb',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999])
+
+#  nond_sampling study
+md.qmu.method = dakota_method.dakota_method('nond_samp')
+md.qmu.method = dmeth_params_set(md.qmu.method,'seed',1234,'samples',3,'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','TotalSmb','IceMass']
+
+#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(3):
+	md.results.dakota.moments.append(md.results.dakota.dresp_out[i].mean)
+
+for i in range(3):
+	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-jpl/test/NightlyRun/test250.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test250.py	(revision 23130)
+++ /issm/trunk-jpl/test/NightlyRun/test250.py	(revision 23130)
@@ -0,0 +1,100 @@
+#Test Name: SquareShelfTranForceNeg2dDakotaSampLinearPart
+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 = md.mesh.numberofvertices
+md = partitioner(md,'package','linear')
+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-jpl/test/NightlyRun/test251.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test251.py	(revision 23130)
+++ /issm/trunk-jpl/test/NightlyRun/test251.py	(revision 23130)
@@ -0,0 +1,97 @@
+#Test Name: SquareShelfTranForceNeg2dDakotaLocalLinearPart
+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 = md.mesh.numberofvertices
+md = partitioner(md,'package','linear')
+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-jpl/test/NightlyRun/test412.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test412.py	(revision 23130)
+++ /issm/trunk-jpl/test/NightlyRun/test412.py	(revision 23130)
@@ -0,0 +1,68 @@
+#Test Name: SquareSheetShelfDiadSSA3dDakota
+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 *
+
+from normal_uncertain import *
+from response_function import *
+
+md = triangle(model(),'../Exp/Square.exp',300000.)
+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 = md.mesh.numberofvertices
+md = partitioner(md,'package','linear','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)
+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
+
+#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-jpl/test/NightlyRun/test440.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test440.py	(revision 23130)
+++ /issm/trunk-jpl/test/NightlyRun/test440.py	(revision 23130)
@@ -0,0 +1,75 @@
+#Test Name: SquareSheetShelfSteaHO
+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 normal_uncertain import *
+from response_function import *
+
+#Test Name: SquareSheetShelfDakotaScaledResponseLinearPart
+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',oshostname(),'np',3)
+
+#partitioning
+md.qmu.numberofpartitions = md.mesh.numberofvertices
+md = partitioner(md,'package','linear')
+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]
+
+#Fields and tolerances to track changes
+field_names     = ['Thickness']
+field_tolerances= [1e-10]
+field_values = [thickness]
Index: /issm/trunk-jpl/test/NightlyRun/test511.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test511.py	(revision 23129)
+++ /issm/trunk-jpl/test/NightlyRun/test511.py	(revision 23130)
@@ -8,5 +8,4 @@
 from setflowequation import *
 from solve import *
-
 
 md=triangle(model(),'../Exp/Pig.exp',11000.)
@@ -41,11 +40,11 @@
 field_names     =['Gradient','Misfits','MaterialsRheologyB','Pressure','Vel','Vx','Vy']
 field_tolerances=[5e-11,5e-11,5e-11,1e-09,1e-11,5e-11,1e-11]
-field_values=[\
-	md.results.StressbalanceSolution.Gradient1,\
-	md.results.StressbalanceSolution.J,\
-	md.results.StressbalanceSolution.MaterialsRheologyBbar,\
-	md.results.StressbalanceSolution.Pressure,\
-	md.results.StressbalanceSolution.Vel,\
-	md.results.StressbalanceSolution.Vx,\
-	md.results.StressbalanceSolution.Vy,\
+field_values=[
+	md.results.StressbalanceSolution.Gradient1,
+	md.results.StressbalanceSolution.J,
+	md.results.StressbalanceSolution.MaterialsRheologyBbar,
+	md.results.StressbalanceSolution.Pressure,
+	md.results.StressbalanceSolution.Vel,
+	md.results.StressbalanceSolution.Vx,
+	md.results.StressbalanceSolution.Vy,
 ]
Index: /issm/trunk-jpl/test/NightlyRun/test613.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test613.py	(revision 23129)
+++ /issm/trunk-jpl/test/NightlyRun/test613.py	(revision 23130)
@@ -19,5 +19,5 @@
 
 #Ice sheet only
-md=md.extract(md.mask.groundedice_levelset>0.)
+md=model.extract(md,md.mask.groundedice_levelset>0.)
 pos=np.nonzero(md.mesh.vertexonboundary)
 md.balancethickness.spcthickness[pos]=md.geometry.thickness[pos]
Index: /issm/trunk-jpl/test/NightlyRun/test703.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test703.py	(revision 23129)
+++ /issm/trunk-jpl/test/NightlyRun/test703.py	(revision 23130)
@@ -59,5 +59,5 @@
 md.stressbalance.spcvz = np.nan * np.ones((md.mesh.numberofvertices))
 md.stressbalance.referential = np.nan * np.ones((md.mesh.numberofvertices,6))
-md.stressbalance.loadingforce = 0 * np.ones((md.mesh.numberofvertices,3))
+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.