Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 5245)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 5246)
@@ -346,5 +346,5 @@
 		error('model not consistent: control_type should be a string');
 	end
-	if ~(strcmpi(md.control_type,'rheology_B') | strcmpi(md.control_type,'drag_coefficient')),
+	if ~(strcmpi(md.control_type,'rheology_B') | strcmpi(md.control_type,'drag_coefficient') | strcmpi(md.control_type,'dhdt')),
 		error('model not consistent: control_type should be rheology_B or drag_coefficient');
 	end
Index: /issm/trunk/src/m/solutions/CorePointerFromSolutionEnum.m
===================================================================
--- /issm/trunk/src/m/solutions/CorePointerFromSolutionEnum.m	(revision 5246)
+++ /issm/trunk/src/m/solutions/CorePointerFromSolutionEnum.m	(revision 5246)
@@ -0,0 +1,21 @@
+function [solutioncore]=CorePointerFromSolutionEnum(solutiontype),
+%COREPOINTERFROMSOLUTIONENUM - returns solution_core function
+%
+%   Usage:
+%      [solutioncore]=CorePointerFromSolutionEnum(solutiontype);
+
+switch solutiontype,
+
+	case DiagnosticSolutionEnum,         solutioncore='diagnostic_core';
+	case SteadystateSolutionEnum,        solutioncore='steadystate_core';
+	case ThermalSolutionEnum,            solutioncore='thermal_core';
+	case PrognosticSolutionEnum,         solutioncore='prognostic_core';
+	case BalancedthicknessSolutionEnum,  solutioncore='balancedthickness_core';
+	case BalancedvelocitiesSolutionEnum, solutioncore='balancedvelocities_core';
+	case SurfaceSlopeSolutionEnum,       solutioncore='surfaceslope_core';
+	case BedSlopeSolutionEnum,           solutioncore='bedslope_core';
+	case Transient2DSolutionEnum,        solutioncore='transient2d_core';
+	case Transient3DSolutionEnum,        solutioncore='transient3d_core';
+	otherwise error('%s%s%s',' solution type: ',EnumToString(solutiontype),' not supported yet!');
+
+end
Index: /issm/trunk/src/m/solutions/SolutionConfiguration.m
===================================================================
--- /issm/trunk/src/m/solutions/SolutionConfiguration.m	(revision 5245)
+++ /issm/trunk/src/m/solutions/SolutionConfiguration.m	(revision 5246)
@@ -5,4 +5,7 @@
 %      [analyses, numanalyses, solutioncore]=SolutionConfiguration(solutiontype);
 
+%get solution_core string
+solutioncore=CorePointerFromSolutionEnum(solutiontype);
+
 switch solutiontype,
 
@@ -10,50 +13,40 @@
 		numanalyses=6;
 		analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticStokesAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum];
-		solutioncore='diagnostic_core';
 
 	case SteadystateSolutionEnum,
 		numanalyses=8; 
 		analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticStokesAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum];
-		solutioncore='steadystate_core';
 
 	case ThermalSolutionEnum,
 		numanalyses=2; 
 		analyses=[ThermalAnalysisEnum;MeltingAnalysisEnum];
-		solutioncore='thermal_core';
 
 	case PrognosticSolutionEnum,
 		numanalyses=1; 
 		analyses=[PrognosticAnalysisEnum];
-		solutioncore='prognostic_core';
 
 	case BalancedthicknessSolutionEnum,
 		numanalyses=1; 
 		analyses=[BalancedthicknessAnalysisEnum];
-		solutioncore='balancedthickness_core';
 
 	case BalancedvelocitiesSolutionEnum,
 		numanalyses=1; 
 		analyses=[BalancedvelocitiesAnalysisEnum];
-		solutioncore='balancedvelocities_core';
 
 	case SurfaceSlopeSolutionEnum,
 		numanalyses=1; 
 		analyses=[SurfaceSlopeAnalysisEnum];
-		solutioncore='surfaceslope_core';
 
 	case BedSlopeSolutionEnum,
 		numanalyses=1; 
 		analyses=[BedSlopeAnalysisEnum];
-		solutioncore='bedslope_core';
 
 	case Transient2DSolutionEnum,
 		numanalyses=7; 
 		analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticStokesAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;PrognosticAnalysisEnum];
-		solutioncore='transient2d_core';
 
 	case Transient3DSolutionEnum,
 		numanalyses=9; 
 		analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticStokesAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;PrognosticAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum];
-		solutioncore='transient3d_core';
 
 	otherwise
Index: /issm/trunk/src/m/solutions/balancedvelocities_core.m
===================================================================
--- /issm/trunk/src/m/solutions/balancedvelocities_core.m	(revision 5245)
+++ /issm/trunk/src/m/solutions/balancedvelocities_core.m	(revision 5246)
@@ -18,6 +18,5 @@
 	if solution_type==BalancedvelocitiesSolutionEnum,
 		displaystring(verbose,'\n%s',['saving results...']);
-		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum);
-		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum);
+		femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VelEnum);
 	end
 
Index: /issm/trunk/src/m/solutions/control_core.m
===================================================================
--- /issm/trunk/src/m/solutions/control_core.m	(revision 5245)
+++ /issm/trunk/src/m/solutions/control_core.m	(revision 5246)
@@ -6,11 +6,9 @@
 %
 
-	%Preprocess models
-	femmodel=stokescontrolinit(femmodel);
-
 	%recover parameters common to all solutions
 	verbose=femmodel.parameters.Verbose;
 	nsteps=femmodel.parameters.NSteps;
 	control_type=femmodel.parameters.ControlType;
+	solution_type=femmodel.parameters.SolutionType;
 	fit=femmodel.parameters.Fit;
 	optscal=femmodel.parameters.OptScal;
@@ -20,5 +18,4 @@
 	tolx=femmodel.parameters.TolX;
 	cm_gradient=femmodel.parameters.CmGradient;
-	control_steady=femmodel.parameters.ControlSteady;
 	dim=femmodel.parameters.Dim;
 
@@ -30,4 +27,12 @@
 	J=zeros(nsteps,1);
 
+	%Get core from solution type
+	solutioncore=CorePointerFromSolutionEnum(solution_type);
+
+	%Preprocess models
+	if(solution_type==SteadystateSolutionEnum || solution_type==DiagnosticSolutionEnum)
+		femmodel=stokescontrolinit(femmodel);
+	end
+
 	for n=1:nsteps,
 
@@ -36,5 +41,5 @@
 
 		%In case we are running a steady state control method, compute new temperature field using new parameter distribution: 
-		if control_steady;
+		if (solution_type==SteadystateSolutionEnum)
 			femmodel=steadystate_core(femmodel);
 		end
@@ -67,9 +72,5 @@
 	%generate output
 	displaystring(verbose,'\n%s',['      preparing final velocity solution...']);
-	if control_steady,
-		femmodel=steadystate_core(femmodel);
-	else
-		femmodel=diagnostic_core(femmodel);
-	end
+	eval(['femmodel=' solutioncore '(femmodel);']);
 
 	%Some results not computed by diagnostic or steadystate