Index: /issm/trunk-jpl/test/SandBox/Makefile
===================================================================
--- /issm/trunk-jpl/test/SandBox/Makefile	(revision 24757)
+++ /issm/trunk-jpl/test/SandBox/Makefile	(revision 24757)
@@ -0,0 +1,2 @@
+clean: 
+	rm -rf *.tar.gz *.errlog *.outlog *.bin *.qmu.in *.qmu.out *qmu.err *.queue *.toolkits
Index: /issm/trunk-jpl/test/SandBox/testbayes.m
===================================================================
--- /issm/trunk-jpl/test/SandBox/testbayes.m	(revision 24757)
+++ /issm/trunk-jpl/test/SandBox/testbayes.m	(revision 24757)
@@ -0,0 +1,74 @@
+%Test Name: Test Queos
+md=triangle(model(),'../Exp/Square.exp',150000.);
+md=setmask(md,'../Exp/SquareShelf.exp','');
+md=parameterize(md,'../Par/SquareSheetShelf.par');
+md=setflowequation(md,'SSA','all');
+md.cluster=generic('name',oshostname(),'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
+md.qmu.isdakota=1;
+
+%partitioning
+md.qmu.numberofpartitions=md.mesh.numberofvertices;
+md=partitioner(md,'package','linear','npart',md.qmu.numberofpartitions);
+md.qmu.vpartition=md.qmu.vpartition-1;
+
+
+%input/output:
+md.qmu.variables.drag_coefficient=normal_uncertain('scaled_FrictionCoefficient',1,0.01);
+md.qmu.responses.misfit=calibration_function('MaxVel');
+%md.qmu.responses.MaxVel=response_function('MaxVel',[],[0.0001 0.001 0.01 0.25 0.5 0.75 0.99 0.999 0.9999]);
+
+%bayesian calibration:
+md.qmu.method     =dakota_method('bayes_calibration');
+md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),...
+'seed',1234,...
+'samples',20,...
+'queso',true,...
+'metropolis_hastings',true,...
+'proposal_covariance',true,...
+'diagonal',true,...
+'values',[num2str(numel(fields(md.qmu.variables))*md.qmu.numberofpartitions) ' * 0.00025']);
+
+%direct solve: 
+md.qmu.params.direct=true;
+md.qmu.params.interval_type='forward';
+md.qmu.params.tabular_graphics_data=true;
+md.qmu.params.analysis_driver='matlab';
+md.qmu.params.evaluation_scheduling='master';
+md.qmu.params.processors_per_evaluation=2;
+
+md.qmu.isdakota=1;
+
+md=solve(md,'Stressbalance','overwrite','y');
+
+%Fields and tolerances to track changes
+md.qmu.results=md.results.dakota;
+
+error;
+
+%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=1:8,
+	md.results.dakota.montecarlo=[md.results.dakota.montecarlo md.results.dakota.dresp_out(i).mean];
+end
+for i=1:8,
+	md.results.dakota.montecarlo=[md.results.dakota.montecarlo md.results.dakota.dresp_out(i).stddev];
+end
+field_names     ={'montecarlo'};
+field_tolerances={1e-11};
+field_values={...
+         md.results.dakota.montecarlo,...
+	};