Index: /issm/trunk-jpl/test/NightlyRun/test445.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test445.m	(revision 23440)
+++ /issm/trunk-jpl/test/NightlyRun/test445.m	(revision 23441)
@@ -20,4 +20,5 @@
 %variables
 md.qmu.variables.neff=normal_uncertain('scaled_FrictionEffectivePressure',1,.05);
+md.qmu.variables.geoflux=normal_uncertain('scaled_BasalforcingsGeothermalflux',1,.05);
 
 %responses
Index: /issm/trunk-jpl/test/NightlyRun/test445.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test445.py	(revision 23441)
+++ /issm/trunk-jpl/test/NightlyRun/test445.py	(revision 23441)
@@ -0,0 +1,97 @@
+#Test Name: SquareSheetShelfSteaEnthalpyHO3dDakotaSampNeff 
+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 *
+from ContourToMesh 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.extrude(3,2.)
+md=setflowequation(md,'HO','all')
+md.cluster=generic('name',gethostname(),'np',3)
+md.timestepping.time_step=0.
+md.thermal.isenthalpy=1
+md.thermal.isdynamicbasalspc=1
+md.initialization.waterfraction=np.zeros((md.mesh.numberofvertices))
+md.initialization.watercolumn=np.zeros((md.mesh.numberofvertices))
+
+md.friction.coupling=3
+md.friction.effective_pressure=md.materials.rho_ice*md.constants.g*md.geometry.thickness+md.materials.rho_water*md.constants.g*md.geometry.base
+
+#dakota version
+version = IssmConfig('_DAKOTA_VERSION_')
+version = float(version[0])
+
+#variables
+md.qmu.variables.neff=normal_uncertain.normal_uncertain('scaled_FrictionEffectivePressure',1,.05)
+md.qmu.variables.geoflux=normal_uncertain.normal_uncertain('scaled_BasalforcingsGeothermalflux',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])
+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','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 = 'steadystate'
+	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,'Steadystate','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(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      = ['montecarlo']
+field_tolerances = [1e-11]
+field_values = [md.results.dakota.montecarlo]