Index: issm/trunk-jpl/src/c/CMakeLists.txt
===================================================================
--- issm/trunk-jpl/src/c/CMakeLists.txt	(revision 16517)
+++ issm/trunk-jpl/src/c/CMakeLists.txt	(revision 16518)
@@ -72,5 +72,5 @@
                           $ENV{ISSM_DIR}/src/c/modules/VecMergex 
                           $ENV{ISSM_DIR}/src/c/modules/Mergesolutionfromftogx 
-                          $ENV{ISSM_DIR}/src/c/analyses 
+                          $ENV{ISSM_DIR}/src/c/cores 
                           $ENV{ISSM_DIR}/src/c/solutionsequences 
                           $ENV{ISSM_DIR}/src/c/solutionsequences)
@@ -214,12 +214,12 @@
 					./modules/VecMergex/VecMergex.cpp
 					./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp
-					./analyses/ProcessArguments.cpp
-					./analyses/ResetBoundaryConditions.cpp
-					./analyses/AnalysisConfiguration.cpp
-					./analyses/WrapperCorePointerFromSolutionEnum.cpp
-					./analyses/CorePointerFromSolutionEnum.cpp
-					./analyses/EnvironmentInit.cpp
-					./analyses/EnvironmentFinalize.cpp
-					./analyses/ad_core.cpp
+					./cores/ProcessArguments.cpp
+					./cores/ResetBoundaryConditions.cpp
+					./cores/AnalysisConfiguration.cpp
+					./cores/WrapperCorePointerFromSolutionEnum.cpp
+					./cores/CorePointerFromSolutionEnum.cpp
+					./cores/EnvironmentInit.cpp
+					./cores/EnvironmentFinalize.cpp
+					./cores/ad_core.cpp
 					./solutionsequences/solutionsequence_linear.cpp
 					./solutionsequences/solutionsequence_nonlinear.cpp
@@ -238,13 +238,13 @@
 					  ./modules/AverageOntoPartitionx/AverageOntoPartitionx.cpp
 					  				  ./modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp
-					  ./analyses/dakota_core.cpp
-					  				  ./analyses/DakotaSpawnCore.cpp)
+					  ./cores/dakota_core.cpp
+					  				  ./cores/DakotaSpawnCore.cpp)
 #}}}
 #Transient sources  {{{
 set(transient_sources ./modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp 
-							./analyses/transient_core.cpp)
+							./cores/transient_core.cpp)
 #}}}
 #Steadystate sources  {{{
-set(steadystate_sources ./analyses/steadystate_core.cpp)
+set(steadystate_sources ./cores/steadystate_core.cpp)
 #}}}
 #Prognostic sources  {{{
@@ -253,5 +253,5 @@
 					      ./modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp
 					      ./modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp
-							./analyses/prognostic_core.cpp)
+							./cores/prognostic_core.cpp)
 #}}}
 #Thermal sources  {{{
@@ -271,5 +271,5 @@
 					   ./modules/ConstraintsStatex/ThermalIsPresent.cpp
 					   ./modules/ResetConstraintsx/ThermalConstraintsReset.cpp 
-					   ./analyses/thermal_core.cpp
+					   ./cores/thermal_core.cpp
 					   ./solutionsequences/solutionsequence_thermal_nonlinear.cpp)
 #}}}
@@ -298,11 +298,11 @@
 					  ./shared/Numerics/BrentSearch.cpp
 					  ./shared/Numerics/OptimalSearch.cpp 
-					  ./analyses/control_core.cpp
-					  ./analyses/controltao_core.cpp
-					  ./analyses/objectivefunction.cpp
-					  ./analyses/gradient_core.cpp
-					  ./analyses/adjointdiagnostic_core.cpp
-					  ./analyses/adjointbalancethickness_core.cpp
-					  ./analyses/AdjointCorePointerFromSolutionEnum.cpp
+					  ./cores/control_core.cpp
+					  ./cores/controltao_core.cpp
+					  ./cores/objectivefunction.cpp
+					  ./cores/gradient_core.cpp
+					  ./cores/adjointdiagnostic_core.cpp
+					  ./cores/adjointbalancethickness_core.cpp
+					  ./cores/AdjointCorePointerFromSolutionEnum.cpp
 					  ./solutionsequences/solutionsequence_adjoint_linear.cpp)
 #}}}
@@ -323,5 +323,5 @@
 							./modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp 
 							./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp 
-							./analyses/hydrology_core.cpp
+							./cores/hydrology_core.cpp
 							./solutionsequences/solutionsequence_hydro_nonlinear.cpp)
 #}}}
@@ -339,5 +339,5 @@
 					      ./modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp 
 							./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp 
-							./analyses/diagnostic_core.cpp
+							./cores/diagnostic_core.cpp
 							./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp)
 #}}}
@@ -347,6 +347,6 @@
 					    ./modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp
 						 ./modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp
-						 ./analyses/balancethickness_core.cpp 
-						 ./analyses/dummy_core.cpp)
+						 ./cores/balancethickness_core.cpp 
+						 ./cores/dummy_core.cpp)
 #}}}
 #Slope sources  {{{
@@ -359,9 +359,9 @@
 					  ./modules/ModelProcessorx/SurfaceSlope/CreateConstraintsSurfaceSlope.cpp
 					  ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp
-					  ./analyses/surfaceslope_core.cpp
-					  ./analyses/bedslope_core.cpp)
+					  ./cores/surfaceslope_core.cpp
+					  ./cores/bedslope_core.cpp)
 #}}}
 #Gia sources  {{{
-set(gia_sources  ./analyses/gia_core.cpp
+set(gia_sources  ./cores/gia_core.cpp
 					./modules/ModelProcessorx/Gia/UpdateElementsGia.cpp
 					./modules/ModelProcessorx/Gia/CreateNodesGia.cpp 
Index: issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- issm/trunk-jpl/src/c/Makefile.am	(revision 16517)
+++ issm/trunk-jpl/src/c/Makefile.am	(revision 16518)
@@ -330,12 +330,12 @@
 					./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
 					./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.h\
-					./analyses/ProcessArguments.cpp\
-					./analyses/ResetBoundaryConditions.cpp\
-					./analyses/AnalysisConfiguration.cpp\
-					./analyses/WrapperCorePointerFromSolutionEnum.cpp\
-					./analyses/CorePointerFromSolutionEnum.cpp\
-					./analyses/EnvironmentInit.cpp\
-					./analyses/EnvironmentFinalize.cpp\
-					./analyses/ad_core.cpp\
+					./cores/ProcessArguments.cpp\
+					./cores/ResetBoundaryConditions.cpp\
+					./cores/AnalysisConfiguration.cpp\
+					./cores/WrapperCorePointerFromSolutionEnum.cpp\
+					./cores/CorePointerFromSolutionEnum.cpp\
+					./cores/ad_core.cpp\
+					./main/EnvironmentInit.cpp\
+					./main/EnvironmentFinalize.cpp\
 					./solutionsequences/solutionsequence_linear.cpp\
 					./solutionsequences/solutionsequence_nonlinear.cpp\
@@ -362,14 +362,14 @@
 					  ./modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp\
 					  ./modules/ModelProcessorx/Dakota/UpdateElementsAndMaterialsDakota.cpp\
-					  ./analyses/dakota_core.cpp\
-					  ./analyses/DakotaSpawnCore.h\
-					  ./analyses/DakotaSpawnCore.cpp
+					  ./cores/dakota_core.cpp\
+					  ./cores/DakotaSpawnCore.h\
+					  ./cores/DakotaSpawnCore.cpp
 #}}}
 #Transient sources  {{{
 transient_sources  = ./modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp \
-							./analyses/transient_core.cpp
+							./cores/transient_core.cpp
 #}}}
 #Steadystate sources  {{{
-steadystate_sources = ./analyses/steadystate_core.cpp
+steadystate_sources = ./cores/steadystate_core.cpp
 #}}}
 #Masstransport sources  {{{
@@ -394,7 +394,7 @@
 								./modules/ModelProcessorx/ExtrudeFromTop/CreateConstraintsExtrudeFromTop.cpp\
 								./modules/ModelProcessorx/ExtrudeFromTop/CreateLoadsExtrudeFromTop.cpp\
-								./analyses/masstransport_core.cpp\
-								./analyses/extrudefrombase_core.cpp\
-								./analyses/extrudefromtop_core.cpp
+								./cores/masstransport_core.cpp\
+								./cores/extrudefrombase_core.cpp\
+								./cores/extrudefromtop_core.cpp
 
 #}}}
@@ -414,5 +414,5 @@
 					   ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\
 					   ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\
-					   ./analyses/thermal_core.cpp\
+					   ./cores/thermal_core.cpp\
 					   ./solutionsequences/solutionsequence_thermal_nonlinear.cpp
 #}}}
@@ -460,11 +460,11 @@
 					  ./shared/Numerics/BrentSearch.cpp\
 					  ./shared/Numerics/OptimalSearch.cpp \
-					  ./analyses/control_core.cpp\
-					  ./analyses/controltao_core.cpp\
-					  ./analyses/objectivefunction.cpp\
-					  ./analyses/gradient_core.cpp\
-					  ./analyses/adjointstressbalance_core.cpp\
-					  ./analyses/adjointbalancethickness_core.cpp\
-					  ./analyses/AdjointCorePointerFromSolutionEnum.cpp\
+					  ./cores/control_core.cpp\
+					  ./cores/controltao_core.cpp\
+					  ./cores/objectivefunction.cpp\
+					  ./cores/gradient_core.cpp\
+					  ./cores/adjointstressbalance_core.cpp\
+					  ./cores/adjointbalancethickness_core.cpp\
+					  ./cores/AdjointCorePointerFromSolutionEnum.cpp\
 					  ./solutionsequences/solutionsequence_adjoint_linear.cpp
 
@@ -486,5 +486,5 @@
 							./modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp \
 							./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp \
-							./analyses/hydrology_core.cpp\
+							./cores/hydrology_core.cpp\
 							./solutionsequences/solutionsequence_hydro_nonlinear.cpp
 #}}}
@@ -503,5 +503,5 @@
 					      ./modules/ModelProcessorx/StressbalanceSIA/CreateConstraintsStressbalanceSIA.cpp \
 							./modules/ModelProcessorx/StressbalanceSIA/CreateLoadsStressbalanceSIA.cpp \
-							./analyses/stressbalance_core.cpp\
+							./cores/stressbalance_core.cpp\
 							./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp
 #}}}
@@ -515,7 +515,7 @@
 						 ./modules/ModelProcessorx/Balancevelocity/CreateConstraintsBalancevelocity.cpp\
 						 ./modules/ModelProcessorx/Balancevelocity/CreateLoadsBalancevelocity.cpp\
-						 ./analyses/balancethickness_core.cpp \
-						 ./analyses/balancevelocity_core.cpp \
-						 ./analyses/dummy_core.cpp
+						 ./cores/balancethickness_core.cpp \
+						 ./cores/balancevelocity_core.cpp \
+						 ./cores/dummy_core.cpp
 #}}}
 #Slope sources  {{{
@@ -528,6 +528,6 @@
 					  ./modules/ModelProcessorx/SurfaceSlope/CreateConstraintsSurfaceSlope.cpp\
 					  ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp\
-					  ./analyses/surfaceslope_core.cpp\
-					  ./analyses/bedslope_core.cpp
+					  ./cores/surfaceslope_core.cpp\
+					  ./cores/bedslope_core.cpp
 #}}}
 #MeshDeformation sources  {{{
@@ -536,8 +536,8 @@
 					  ./modules/ModelProcessorx/MeshDeformation/CreateConstraintsMeshDeformation.cpp\
 					  ./modules/ModelProcessorx/MeshDeformation/CreateLoadsMeshDeformation.cpp\
-					  ./analyses/meshdeformation_core.cpp
+					  ./cores/meshdeformation_core.cpp
 #}}}
 #Gia sources  {{{
-gia_sources =  ./analyses/gia_core.cpp\
+gia_sources =  ./cores/gia_core.cpp\
 					./modules/ModelProcessorx/Gia/UpdateElementsGia.cpp\
 					./modules/ModelProcessorx/Gia/CreateNodesGia.cpp \
@@ -557,5 +557,5 @@
 #}}}
 #Damage sources  {{{
-damage_sources =  ./analyses/damage_core.cpp\
+damage_sources =  ./cores/damage_core.cpp\
 					./modules/ModelProcessorx/Damage/UpdateElementsDamage.cpp\
 					./modules/ModelProcessorx/Damage/CreateNodesDamage.cpp \
Index: issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 16517)
+++ issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 16518)
@@ -10,5 +10,5 @@
 
 #include <stdio.h>
-#include "../analyses/analyses.h"
+#include "../cores/cores.h"
 #include "../shared/io/io.h"
 #include "./classes.h"
Index: issm/trunk-jpl/src/c/cores/AdjointCorePointerFromSolutionEnum.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/AdjointCorePointerFromSolutionEnum.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/AdjointCorePointerFromSolutionEnum.cpp	(revision 16518)
@@ -0,0 +1,46 @@
+/*!\file:  AdjointCorePointerFromSolutionEnum.cpp
+ * \brief: return type of analyses, number of analyses and core solution function.
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void AdjointCorePointerFromSolutionEnum(void (**padjointcore)(FemModel*),int solutiontype){
+
+	/*output: */
+	void (*adjointcore)(FemModel*)=NULL;
+
+	switch(solutiontype){
+
+		case StressbalanceSolutionEnum:
+			adjointcore=&adjointstressbalance_core;
+			break;
+		case SteadystateSolutionEnum:
+			adjointcore=&adjointstressbalance_core;
+			break;
+		case BalancethicknessSolutionEnum:
+			adjointcore=&adjointbalancethickness_core;
+			break;
+		case BalancethicknessSoftSolutionEnum:
+			adjointcore=&dummy_core;
+			break;
+		default:
+			_error_("No adjoint has been implemented for solution " << EnumToStringx(solutiontype) << " yet");
+			break;
+	}
+
+	/*Assign output pointer:*/
+	_assert_(padjointcore);
+	*padjointcore=adjointcore;
+
+}
Index: issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp	(revision 16518)
@@ -0,0 +1,144 @@
+/*!\file:  AnalysisConfiguration.cpp
+ * \brief: return type of analyses, number of analyses 
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void AnalysisConfiguration(int** panalyses,int* pnumanalyses, int solutiontype){
+
+	/*output: */
+	int  numanalyses;
+	int* analyses=NULL;
+
+	/*Analyses lists*/
+	switch(solutiontype){
+
+		case StressbalanceSolutionEnum:
+			numanalyses=5;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=StressbalanceAnalysisEnum;
+			analyses[1]=StressbalanceVerticalAnalysisEnum;
+			analyses[2]=StressbalanceSIAAnalysisEnum;
+			analyses[3]=SurfaceSlopeAnalysisEnum;
+			analyses[4]=BedSlopeAnalysisEnum;
+			break;
+
+		case SteadystateSolutionEnum:
+			numanalyses=8;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=StressbalanceAnalysisEnum;
+			analyses[1]=StressbalanceVerticalAnalysisEnum;
+			analyses[2]=StressbalanceSIAAnalysisEnum;
+			analyses[3]=SurfaceSlopeAnalysisEnum;
+			analyses[4]=BedSlopeAnalysisEnum;
+			analyses[5]=EnthalpyAnalysisEnum;
+			analyses[6]=ThermalAnalysisEnum;
+			analyses[7]=MeltingAnalysisEnum;
+			break;
+
+		case ThermalSolutionEnum:
+			numanalyses=3;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=ThermalAnalysisEnum;
+			analyses[1]=MeltingAnalysisEnum;
+			analyses[2]=EnthalpyAnalysisEnum;
+			break;
+
+		case HydrologySolutionEnum:
+			numanalyses=5;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=HydrologyShreveAnalysisEnum;
+			analyses[1]=HydrologyDCInefficientAnalysisEnum;
+			analyses[2]=HydrologyDCEfficientAnalysisEnum;
+			analyses[3]=SurfaceSlopeAnalysisEnum;
+			analyses[4]=BedSlopeAnalysisEnum;
+			break;
+
+		case MasstransportSolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=MasstransportAnalysisEnum;
+			break;
+
+		case BalancethicknessSolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=BalancethicknessAnalysisEnum;
+			break;
+
+		case BalancethicknessSoftSolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=BalancethicknessAnalysisEnum;
+			break;
+
+		case BalancevelocitySolutionEnum:
+			numanalyses=2;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=BalancevelocityAnalysisEnum;
+			analyses[1]=SurfaceSlopeAnalysisEnum;
+			break;
+
+		case SurfaceSlopeSolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=SurfaceSlopeAnalysisEnum;
+			break;
+
+		case BedSlopeSolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=BedSlopeAnalysisEnum;
+			break;
+
+		case GiaSolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=GiaAnalysisEnum;
+			break;
+		
+		case DamageEvolutionSolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=DamageEvolutionAnalysisEnum;
+			break;
+
+		case TransientSolutionEnum:
+			numanalyses=13;
+			analyses=xNew<int>(numanalyses);
+			analyses[ 0]=StressbalanceAnalysisEnum;
+			analyses[ 1]=StressbalanceVerticalAnalysisEnum;
+			analyses[ 2]=StressbalanceSIAAnalysisEnum;
+			analyses[ 3]=SurfaceSlopeAnalysisEnum;
+			analyses[ 4]=BedSlopeAnalysisEnum;
+			analyses[ 5]=ThermalAnalysisEnum;
+			analyses[ 6]=MeltingAnalysisEnum;
+			analyses[ 7]=EnthalpyAnalysisEnum;
+			analyses[ 8]=MasstransportAnalysisEnum;
+			analyses[ 9]=FreeSurfaceBaseAnalysisEnum;
+			analyses[10]=FreeSurfaceTopAnalysisEnum;
+			analyses[11]=ExtrudeFromBaseAnalysisEnum;
+			analyses[12]=ExtrudeFromTopAnalysisEnum;
+			break;
+
+		default:
+			_error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");
+			break;
+	}
+
+	/*Assign output pointers:*/
+	if(pnumanalyses) *pnumanalyses=numanalyses;
+	if(panalyses)    *panalyses=analyses;
+	else              xDelete<int>(analyses);
+}
Index: issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp	(revision 16518)
@@ -0,0 +1,127 @@
+/*!\file:  CorePointerFromSolutionEnum.cpp
+ * \brief: return type of analyses, number of analyses and core solution function.
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype){
+
+	/*output: */
+	void (*solutioncore)(FemModel*)=NULL;
+
+	switch(solutiontype){
+
+		case StressbalanceSolutionEnum:
+			#ifdef _HAVE_STRESSBALANCE_
+			solutioncore=&stressbalance_core;
+			#else
+			_error_("ISSM was not compiled with stressbalance capabilities. Exiting");
+			#endif
+			break;
+		case SteadystateSolutionEnum:
+			#ifdef _HAVE_STEADYSTATE_
+			solutioncore=&steadystate_core;
+			#else
+			_error_("ISSM was not compiled with steady state capabilities. Exiting");
+			#endif
+			break;
+		case ThermalSolutionEnum:
+			#ifdef _HAVE_THERMAL_
+			solutioncore=&thermal_core;
+			#else
+			_error_("ISSM was not compiled with thermal capabilities. Exiting");
+			#endif
+			break;
+		case BalancethicknessSolutionEnum:
+			#ifdef _HAVE_BALANCED_
+			solutioncore=&balancethickness_core;
+			#else
+			_error_("ISSM was not compiled with balanced capabilities. Exiting");
+			#endif
+			break;
+		case BalancethicknessSoftSolutionEnum:
+			#ifdef _HAVE_BALANCED_
+			solutioncore=&dummy_core;
+			#else
+			_error_("ISSM was not compiled with balanced capabilities. Exiting");
+			#endif
+			break;
+		case BalancevelocitySolutionEnum:
+			#ifdef _HAVE_BALANCED_
+			solutioncore=&balancevelocity_core;
+			#else
+			_error_("ISSM was not compiled with balanced capabilities. Exiting");
+			#endif
+			break;
+		case HydrologySolutionEnum:
+			#ifdef _HAVE_HYDROLOGY_
+			solutioncore=&hydrology_core;
+			#else
+			_error_("ISSM was not compiled with hydrology capabilities. Exiting");
+			#endif
+			break;
+		case SurfaceSlopeSolutionEnum:
+			#ifdef _HAVE_SLOPE_
+			solutioncore=&surfaceslope_core;
+			#else
+			_error_("ISSM was not compiled with slope capabilities. Exiting");
+			#endif
+			break;
+		case BedSlopeSolutionEnum:
+			#ifdef _HAVE_SLOPE_
+			solutioncore=&bedslope_core;
+			#else
+			_error_("ISSM was not compiled with slope capabilities. Exiting");
+			#endif
+			break;
+		case TransientSolutionEnum:
+			#ifdef _HAVE_TRANSIENT_
+			solutioncore=&transient_core;
+			#else
+			_error_("ISSM was not compiled with transient capabilities. Exiting");
+			#endif
+			break;
+		case MasstransportSolutionEnum:
+			#ifdef _HAVE_MASSTRANSPORT_
+			solutioncore=&masstransport_core;
+			#else
+			_error_("ISSM was not compiled with masstransport capabilities. Exiting");
+			#endif
+			break;
+
+		case GiaSolutionEnum:
+			#ifdef _HAVE_GIA_
+			solutioncore=&gia_core;
+			#else
+			_error_("ISSM was not compiled with gia capabilities. Exiting");
+			#endif
+			break;
+		case DamageEvolutionSolutionEnum:
+			#ifdef _HAVE_DAMAGE_
+			solutioncore=&damage_core;
+			#else
+			_error_("ISSM was not compiled with damage evolution capabilities. Exiting");
+			#endif
+			break;
+
+		default:
+			_error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");
+			break;
+	}
+
+	/*Assign output pointer:*/
+	_assert_(psolutioncore);
+	*psolutioncore=solutioncore;
+
+}
Index: issm/trunk-jpl/src/c/cores/DakotaSpawnCore.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/DakotaSpawnCore.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/DakotaSpawnCore.cpp	(revision 16518)
@@ -0,0 +1,187 @@
+/*!\file:  DakotaSpawnCore.cpp
+
+ * \brief: run core ISSM solution using Dakota inputs coming from CPU 0.
+ * \sa qmu.cpp DakotaPlugin.cpp
+ *
+ * This routine needs to be understood simultaneously with qmu.cpp and DakotaPlugin. 
+ * DakotaSpawnCoreParallel is called by all CPUS, with CPU 0 holding Dakota variable values, along 
+ * with variable descriptors. 
+ *
+ * DakotaSpawnCoreParallel takes care of broadcasting the variables and their descriptors across the MPI 
+ * ring. Once this is done, we use the variables to modify the inputs for the solution core. 
+ * For ex, if "rho_ice" is provided, for ex 920, we include "rho_ice" in the inputs, then 
+ * call the core with the modified inputs. This is the way we get Dakota to explore the parameter 
+ * spce of the core. 
+ *
+ * Once the core is called, we process the results of the core, and using the processed results, 
+ * we compute response functions. The responses are computed on all CPUS, but they are targeted 
+ * for CPU 0, which will get these values back to the Dakota engine. 
+ *
+ */ 
+
+/*Includes and prototypes: {{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./cores.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+
+void DakotaMPI_Bcast(double** pvariables, char*** pvariables_descriptors,int* pnumvariables, int* pnumresponses);
+void DakotaFree(double** pvariables,char*** pvariables_descriptors,char*** presponses_descriptors,int numvariables,int numresponses);
+/*}}}*/
+
+/*Notice the d_, which prefixes anything that is being provided to us by the Dakota pluggin. Careful. some things are ours, some are dakotas!: */
+int DakotaSpawnCore(double* d_responses, int d_numresponses, double* d_variables, char** d_variables_descriptors,int d_numvariables, void* void_femmodel,int counter){
+
+	char     **responses_descriptors    = NULL;      //these are our! there are only numresponsedescriptors of them, not d_numresponses!!!
+	int        numresponsedescriptors;
+	int        solution_type;
+	bool       control_analysis         = false;
+	void     (*solutioncore)(FemModel*) = NULL;
+	FemModel  *femmodel                 = NULL;
+	bool       nodakotacore             = true;
+
+	/*If counter==-1 on cpu0, it means that the dakota runs are done. In which case, bail out and return 0: */
+	ISSM_MPI_Bcast(&counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
+	if(counter==-1)return 0;
+
+	/*cast void_femmodel to FemModel: */
+	femmodel=reinterpret_cast<FemModel*>(void_femmodel);
+
+	/*retrieve parameters: */
+	femmodel->parameters->FindParam(&responses_descriptors,&numresponsedescriptors,QmuResponsedescriptorsEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
+
+	if(VerboseQmu()) _printf0_("qmu iteration: " << counter << "\n");
+
+	/* only cpu 0, running dakota is providing us with variables and variables_descriptors and numresponses: broadcast onto other cpus: */
+	DakotaMPI_Bcast(&d_variables,&d_variables_descriptors,&d_numvariables,&d_numresponses);
+
+	/*Modify core inputs in objects contained in femmodel, to reflect the dakota variables inputs: */
+	InputUpdateFromDakotax(femmodel,d_variables,d_variables_descriptors,d_numvariables);
+
+	/*Determine solution sequence: */
+	if(VerboseQmu()) _printf0_("Starting " << EnumToStringx(solution_type) << " core:\n");
+	WrapperCorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type,nodakotacore);
+
+	/*Run the core solution sequence: */
+	solutioncore(femmodel);
+
+	/*compute responses: */
+	if(VerboseQmu()) _printf0_("compute dakota responses:\n");
+	femmodel->DakotaResponsesx(d_responses,responses_descriptors,numresponsedescriptors,d_numresponses);
+
+	/*Free ressources:*/
+	DakotaFree(&d_variables,&d_variables_descriptors,&responses_descriptors, d_numvariables, numresponsedescriptors);
+
+	return 1; //this is critical! do not return 0, otherwise, dakota_core will stop running!
+}
+
+void DakotaMPI_Bcast(double** pvariables, char*** pvariables_descriptors,int* pnumvariables, int* pnumresponses){ /*{{{*/
+
+	/* * \brief: broadcast variables_descriptors, variables, numvariables and numresponses
+	 * from cpu 0 to all other cpus.
+	 */ 
+
+	int i;
+	int my_rank;
+
+	/*inputs and outputs: */
+	double* variables=NULL;
+	char**  variables_descriptors=NULL;
+	int     numvariables;
+	int     numresponses;
+
+	/*intermediary: */
+	char* string=NULL;
+	int   string_length;
+
+	/*recover my_rank:*/
+	my_rank=IssmComm::GetRank();
+
+	/*recover inputs from pointers: */
+	variables=*pvariables;
+	variables_descriptors=*pvariables_descriptors;
+	numvariables=*pnumvariables;
+	numresponses=*pnumresponses;
+
+	/*numvariables: */
+	ISSM_MPI_Bcast(&numvariables,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
+
+	/*variables:*/
+	if(my_rank!=0)variables=xNew<double>(numvariables);
+	ISSM_MPI_Bcast(variables,numvariables,MPI_DOUBLE,0,IssmComm::GetComm()); 
+
+	/*variables_descriptors: */
+	if(my_rank!=0){
+		variables_descriptors=xNew<char*>(numvariables);
+	}
+	for(i=0;i<numvariables;i++){
+		if(my_rank==0){
+			string=variables_descriptors[i];
+			string_length=(strlen(string)+1)*sizeof(char);
+		}
+		ISSM_MPI_Bcast(&string_length,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
+		if(my_rank!=0)string=xNew<char>(string_length);
+		ISSM_MPI_Bcast(string,string_length,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
+		if(my_rank!=0)variables_descriptors[i]=string;
+	}
+
+	/*numresponses: */
+	ISSM_MPI_Bcast(&numresponses,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
+
+	/*Assign output pointers:*/
+	*pnumvariables=numvariables;
+	*pvariables=variables;
+	*pvariables_descriptors=variables_descriptors;
+	*pnumresponses=numresponses;
+} /*}}}*/
+void DakotaFree(double** pvariables,char*** pvariables_descriptors,char*** presponses_descriptors,int numvariables,int numresponses){ /*{{{*/
+
+	/*\brief DakotaFree: free allocations on other cpus, not done by Dakota.*/
+
+	int i;
+	int my_rank;
+
+	double  *variables             = NULL;
+	char   **variables_descriptors = NULL;
+	char   **responses_descriptors = NULL;
+	char    *string                = NULL;
+
+	/*recover pointers: */
+	variables=*pvariables;
+	variables_descriptors=*pvariables_descriptors;
+	responses_descriptors=*presponses_descriptors;
+
+	/*recover my_rank:*/
+	my_rank=IssmComm::GetRank();
+
+	/*Free variables and variables_descriptors only on cpu !=0*/
+	if(my_rank!=0){
+		xDelete<double>(variables);
+		for(i=0;i<numvariables;i++){
+			string=variables_descriptors[i];
+			xDelete<char>(string);
+		}
+		xDelete<char*>(variables_descriptors);
+	}
+
+	//responses descriptors on every cpu
+	for(i=0;i<numresponses;i++){
+		string=responses_descriptors[i];
+		xDelete<char>(string);
+	}
+	//rest of dynamic allocations.
+	xDelete<char*>(responses_descriptors);
+
+	/*Assign output pointers:*/
+	*pvariables=variables;
+	*pvariables_descriptors=variables_descriptors;
+	*presponses_descriptors=responses_descriptors;
+} /*}}}*/
Index: issm/trunk-jpl/src/c/cores/DakotaSpawnCore.h
===================================================================
--- issm/trunk-jpl/src/c/cores/DakotaSpawnCore.h	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/DakotaSpawnCore.h	(revision 16518)
@@ -0,0 +1,4 @@
+#ifndef _DAKOTA_SPAWN_CORE_
+#define _DAKOTA_SPAWN_CORE_
+int  DakotaSpawnCore(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, void* femmodel,int counter);
+#endif
Index: issm/trunk-jpl/src/c/cores/ProcessArguments.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/ProcessArguments.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/ProcessArguments.cpp	(revision 16518)
@@ -0,0 +1,57 @@
+/*!\file:  ProcessArguments.cpp
+ * \brief: process arguments
+ */ 
+
+#include <stdio.h>
+#include <cstring>
+
+#include "../shared/shared.h"
+
+void ProcessArguments(int* solution_type,char** pbinfilename,char** poutbinfilename,char** ptoolkitsfilename,char** plockfilename,char** prootpath, int argc,char **argv){
+
+	char *modelname      = NULL;
+	char *binfilename    = NULL;
+	char *outbinfilename = NULL;
+	char *toolkitsfilename  = NULL;
+	char *lockfilename   = NULL;
+	char *rootpath       = NULL;
+	char *rootpatharg    = NULL;
+
+	/*Check input arguments*/
+	if(argc<2)_error_("Usage error: no solution requested");
+	if(argc<3)_error_("Usage error: missing execution directory");
+	if(argc<4)_error_("Usage error: missing model name");
+
+	/*Get requested solution*/
+	*solution_type=StringToEnumx(argv[1]);
+
+	rootpatharg=argv[2];
+	if(strcmp(strstr(rootpatharg,"/"),"/")!=0){ 
+		rootpath       = xNew<char>(strlen(rootpatharg)+2); sprintf(rootpath,"%s/",rootpatharg);
+	}
+	else{
+		rootpath       = xNew<char>(strlen(rootpatharg)+1); sprintf(rootpath,"%s",rootpatharg);
+	}
+
+	modelname=argv[3];
+	if(strstr(modelname,rootpath)==NULL){
+		binfilename    = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".bin")   +1); sprintf(binfilename,   "%s%s%s",rootpath,modelname,".bin");
+		outbinfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".outbin")+1); sprintf(outbinfilename,"%s%s%s",rootpath,modelname,".outbin");
+		toolkitsfilename  = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".toolkits") +1); sprintf(toolkitsfilename, "%s%s%s",rootpath,modelname,".toolkits");
+		lockfilename   = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".lock")  +1); sprintf(lockfilename,  "%s%s%s",rootpath,modelname,".lock");
+	}
+	else{
+		binfilename    = xNew<char>(strlen(modelname)+strlen(".bin")   +1); sprintf(binfilename,   "%s%s",modelname,".bin");
+		outbinfilename = xNew<char>(strlen(modelname)+strlen(".outbin")+1); sprintf(outbinfilename,"%s%s",modelname,".outbin");
+		toolkitsfilename  = xNew<char>(strlen(modelname)+strlen(".toolkits") +1); sprintf(toolkitsfilename, "%s%s",modelname,".toolkits");
+		lockfilename   = xNew<char>(strlen(modelname)+strlen(".lock")  +1); sprintf(lockfilename,  "%s%s",modelname,".lock");
+	}
+
+	/*Clean up and assign output pointer*/
+	*pbinfilename=binfilename;
+	*poutbinfilename=outbinfilename;
+	*ptoolkitsfilename=toolkitsfilename;
+	*plockfilename=lockfilename;
+	*prootpath=rootpath;
+
+}
Index: issm/trunk-jpl/src/c/cores/ResetBoundaryConditions.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/ResetBoundaryConditions.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/ResetBoundaryConditions.cpp	(revision 16518)
@@ -0,0 +1,27 @@
+/*!\file: ResetBoundaryConditions.cpp
+ * \brief: change boundary conditions of a model, using a solution vector from another analysis
+ */ 
+
+#include "../classes/classes.h"
+#include "../modules/modules.h"
+#include "../shared/io/io.h"
+
+void ResetBoundaryConditions(FemModel* femmodel, int analysis_type){
+
+	/*variables: */
+	Vector<IssmDouble>* yg = NULL;
+
+	if(VerboseSolution()) _printf0_("   updating boundary conditions...\n");
+
+	/*set current analysis: */
+	femmodel->SetCurrentConfiguration(analysis_type);
+
+	/*retrieve boundary conditions from element inputs :*/
+	GetSolutionFromInputsx(&yg,femmodel);
+
+	/*update spcs using this new vector of constraints: */
+	UpdateDynamicConstraintsx(femmodel->constraints,femmodel->nodes,femmodel->parameters,yg);
+
+	/*Free ressources:*/
+	delete yg;
+}
Index: issm/trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp	(revision 16518)
@@ -0,0 +1,57 @@
+/*!\file:  WrapperCorePointerFromSolutionEnum.cpp
+ * \brief: return type of analyses, number of analyses and core solution function.
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void WrapperCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype,bool nodakotacore){
+
+	/*output: */
+	void (*solutioncore)(FemModel*)=NULL;
+
+	/*parameters: */
+	bool control_analysis=false;
+	bool tao_analysis=false;
+	bool dakota_analysis=false;
+
+	/* retrieve some parameters that tell us whether wrappers are allowed, or whether we return 
+	 * a pure core. Wrappers can be dakota_core (which samples many solution_cores) or control_core (which 
+	 * carries out adjoint based inversion on a certain core: */
+	parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
+	parameters->FindParam(&control_analysis,InversionIscontrolEnum);
+	parameters->FindParam(&tao_analysis,InversionTaoEnum);
+
+	if(nodakotacore)dakota_analysis=false;
+
+	if(dakota_analysis){
+		#ifdef _HAVE_DAKOTA_
+		solutioncore=dakota_core;
+		#else
+		_error_("ISSM was not compiled with dakota support, cannot carry out dakota analysis!");
+		#endif
+	}
+	else if(control_analysis){
+		#ifdef _HAVE_CONTROL_
+		if(tao_analysis) solutioncore=controltao_core;
+		else solutioncore=control_core;
+		#else
+		_error_("ISSM was not compiled with control support, cannot carry out control analysis!");
+		#endif
+	}
+	else CorePointerFromSolutionEnum(&solutioncore,parameters,solutiontype);  /*This means we retrieve a core solution that is not a wrapper*/
+
+	/*Assign output pointer:*/
+	_assert_(psolutioncore);
+	*psolutioncore=solutioncore;
+
+}
Index: issm/trunk-jpl/src/c/cores/ad_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/ad_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/ad_core.cpp	(revision 16518)
@@ -0,0 +1,356 @@
+/*!\file ad_core
+ * \brief: compute outputs from the AD mode,  using our dependents and independents, and drivers available in Adolc.
+ */
+
+/*Includes: {{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <set>
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+/*}}}*/
+
+void ad_core(FemModel* femmodel){
+
+	/*diverse: */
+	int     i;
+	int     dummy;
+	int     num_dependents=0;
+	int     num_independents=0;
+	bool    isautodiff       = false;
+	char   *driver           = NULL;
+	size_t  tape_stats[11];
+
+	/*state variables: */
+	IssmDouble *axp = NULL;
+	double     *xp  = NULL;
+	int my_rank=IssmComm::GetRank();
+
+	/*AD mode on?: */
+	femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
+
+	if(isautodiff){
+
+		#ifdef _HAVE_ADOLC_
+			if(VerboseAutodiff())_printf0_("   start ad core\n");
+
+			/*First, stop tracing: */
+			trace_off();
+		
+			/*retrieve parameters: */
+			femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
+			femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
+	
+			/*if no dependents, no point in running a driver: */
+			if(!(num_dependents*num_independents)) return;
+
+			/*for adolc to run in parallel, we 0 out on rank~=0:*/
+			if (my_rank!=0){
+				num_dependents=0; num_independents=0;
+			}
+			
+			/*Print tape statistics so that user can kill this run if something is off already:*/
+			tapestats(my_rank,tape_stats); //reading of tape statistics
+			if(VerboseAutodiff()){
+				int commSize=IssmComm::GetSize();
+				int *sstats=new int[7];
+				sstats[0]=tape_stats[NUM_OPERATIONS];
+				sstats[1]=tape_stats[OP_FILE_ACCESS];
+				sstats[2]=tape_stats[NUM_LOCATIONS];
+				sstats[3]=tape_stats[LOC_FILE_ACCESS];
+				sstats[4]=tape_stats[NUM_VALUES];
+				sstats[5]=tape_stats[VAL_FILE_ACCESS];
+				sstats[6]=tape_stats[TAY_STACK_SIZE];
+				int *rstats=NULL;
+				if (my_rank==0) rstats=new int[commSize*7];
+				ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
+				if (my_rank==0) {
+					int offset=50;
+					int rOffset=(commSize/10)+1;
+					_printf_("   ADOLC statistics: \n");
+					_printf_("     "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
+					_printf_("     "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
+					_printf_("     "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
+					_printf_("     operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
+					_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
+					for (int r=0;r<commSize;++r)
+					_printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
+					_printf_("     locations: entry size " << sizeof(locint) << " Bytes\n");
+					_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
+					for (int r=0;r<commSize;++r)
+					_printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
+					_printf_("     constant values: entry size " << sizeof(double) << " Bytes\n");
+					_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
+					for (int r=0;r<commSize;++r)
+					_printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
+					_printf_("     Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
+					_printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
+					for (int r=0;r<commSize;++r)
+					_printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
+					delete []rstats;
+				}
+				delete [] sstats;
+			}
+
+			/*retrieve state variable: */
+			femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
+
+			/* driver argument */
+			xp=xNew<double>(num_independents);
+			for(i=0;i<num_independents;i++){
+				xp[i]=reCast<double,IssmDouble>(axp[i]);
+			}
+
+			/*get the EDF pointer:*/
+			ext_diff_fct *anEDF_for_solverx_p=dynamic_cast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
+
+			/*Branch according to AD driver: */
+			femmodel->parameters->FindParam(&driver,AutodiffDriverEnum);
+
+			/* these are always needed regardless of the interpreter */
+			anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
+			anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
+
+			if (strcmp(driver,"fos_forward")==0){ /*{{{*/
+
+				int     anIndepIndex;
+				double *tangentDir         = NULL;
+				double *jacTimesTangentDir = NULL;
+				double *theOutput          = NULL;
+
+				/*retrieve direction index: */
+				femmodel->parameters->FindParam(&anIndepIndex,AutodiffFosForwardIndexEnum);
+
+				if (anIndepIndex<0 || anIndepIndex>=num_independents) _error_("index value for AutodiffFosForwardIndexEnum should be in [0,num_independents-1]");
+
+				tangentDir=xNewZeroInit<double>(num_independents);
+				tangentDir[anIndepIndex]=1.0;
+
+				jacTimesTangentDir=xNew<double>(num_dependents);
+				theOutput=xNew<double>(num_dependents);
+
+				/*set the forward method function pointer: */
+#ifdef _HAVE_GSL_
+				anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;
+#endif
+
+				/*allocate the space for the parameters to invoke the EDF fos_forward:*/
+				anEDF_for_solverx_p->dp_X=xNew<double>(anEDF_for_solverx_p->max_n);
+				anEDF_for_solverx_p->dp_Y=xNew<double>(anEDF_for_solverx_p->max_m);
+
+				/*call driver: */
+				fos_forward(my_rank,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir );
+
+				/*add to results*/
+				femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,jacTimesTangentDir,num_dependents,1,1,0.0));
+
+				/*free resources :*/
+				xDelete(theOutput);
+				xDelete(jacTimesTangentDir);
+				xDelete(tangentDir);
+			} /*}}}*/
+			else if ((strcmp(driver,"fov_forward")==0) || (strcmp(driver,"fov_forward_all")==0)){ /*{{{*/
+
+				int      tangentDirNum;
+				int      dummy;
+				int     *indepIndices  = NULL;
+				double **jacTimesSeed  = NULL;
+				double **seed          = NULL;
+				double  *theOutput     = NULL;
+				std::set<unsigned int> anIndexSet;
+
+				/*retrieve directions:*/
+				if (strcmp(driver,"fov_forward_all")==0){
+					tangentDirNum=num_independents;
+					indepIndices=xNewZeroInit<int>(tangentDirNum);
+					for(i=0;i<num_independents;i++)indepIndices[i]=1;
+				}
+				else{
+					femmodel->parameters->FindParam(&indepIndices,&tangentDirNum,&dummy,AutodiffFovForwardIndicesEnum);
+				}
+
+				/*Some checks: */
+				if (tangentDirNum<1 || tangentDirNum>num_independents) _error_("tangentDirNum should be in [1,num_independents]");
+
+				/* full Jacobian or Jacobian projection:*/
+				jacTimesSeed=xNew<double>(num_dependents,tangentDirNum);
+
+				/*set the forward method function pointers: */
+#ifdef _HAVE_GSL_
+				anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx;
+#endif
+				// anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
+
+				/*allocate the space for the parameters to invoke EDF fov_forward:*/
+				anEDF_for_solverx_p->dpp_X=xNew<double>(anEDF_for_solverx_p->max_n, tangentDirNum);
+				anEDF_for_solverx_p->dpp_Y=xNew<double>(anEDF_for_solverx_p->max_m, tangentDirNum);
+
+				/*seed matrix: */
+				seed=xNewZeroInit<double>(num_independents,tangentDirNum);
+
+				/*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/
+				for (int i=0; i<tangentDirNum; ++i) {
+					/* make sure the index is in range*/
+					if (indepIndices[i]>num_independents) {
+						_error_("indepIndices values must be in [0,num_independents-1]");
+					}
+					if (anIndexSet.find(indepIndices[i])!=anIndexSet.end()) {
+						_error_("duplicate indepIndices values are not allowed until we implement Jacobian decompression");
+					}
+					anIndexSet.insert(indepIndices[i]);
+					/* now populate the seed matrix from the set of independent indices;
+					 * simple setup with a single 1.0 per column and at most a single 1.0 per row*/
+					seed[indepIndices[i]][i]=1.0;
+				}
+
+				/*allocate output: */
+				theOutput=xNew<double>(num_dependents);
+
+				/*call driver: */
+				fov_forward(my_rank,num_dependents,num_independents, tangentDirNum, xp, seed, theOutput, jacTimesSeed );
+				/*Free resources: */
+				xDelete(theOutput);
+				xDelete(indepIndices);
+				xDelete(seed);
+
+				/*add to results: */
+				femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,*jacTimesSeed,num_dependents*tangentDirNum,1,1,0.0));
+
+				/*Free resources: */
+				xDelete(jacTimesSeed);
+				xDelete(indepIndices);
+			} /*}}}*/
+			else if (strcmp(driver,"fos_reverse")==0) { /*{{{*/
+
+				int     aDepIndex=0;
+				double *aWeightVector=NULL;
+				double *weightVectorTimesJac=NULL;
+
+				/*retrieve direction index: */
+				femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum);
+				aWeightVector=xNewZeroInit<double>(num_dependents);
+				if (my_rank==0) {
+					if (aDepIndex<0 || aDepIndex>=num_dependents) _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]");
+					aWeightVector[aDepIndex]=1.0;
+				}
+				weightVectorTimesJac=xNew<double>(num_independents);
+
+				/*set the forward method function pointer: */
+#ifdef _HAVE_GSL_
+				anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
+#endif
+#ifdef _HAVE_MUMPS_
+				anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
+#endif
+
+				/*allocate the space for the parameters to invoke the EDF fos_reverse :*/
+				anEDF_for_solverx_p->dp_U=xNew<double>(anEDF_for_solverx_p->max_m);
+				anEDF_for_solverx_p->dp_Z=xNew<double>(anEDF_for_solverx_p->max_n);
+
+				/*call driver: */
+				fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
+				if(VerboseAutodiff())_printf0_("   done with fos_reverse\n");
+
+				/*add to results*/
+				femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,weightVectorTimesJac,num_independents,1,1,0.0));
+
+				/*free resources :*/
+				xDelete(weightVectorTimesJac);
+				xDelete(aWeightVector);
+			} /*}}}*/
+			else if ((strcmp(driver,"fov_reverse")==0) || (strcmp(driver,"fov_reverse_all")==0)){ /*{{{*/
+
+				int* depIndices=NULL;
+				int weightNum;
+				int dummy;
+				double **weightsTimesJac=NULL;
+				double **weights=NULL;
+				std::set<unsigned int> anIndexSet;
+
+				/*retrieve directions:*/
+				if (strcmp(driver,"fov_reverse_all")==0){
+					weightNum=num_dependents;
+					depIndices=xNewZeroInit<int>(weightNum);
+					for(i=0;i<num_dependents;i++)depIndices[i]=1;
+				}
+				else{
+					femmodel->parameters->FindParam(&depIndices,&weightNum,&dummy,AutodiffFovForwardIndicesEnum);
+				}
+
+				/*Some checks: */
+				if (weightNum<1 || weightNum>num_dependents) _error_("tangentDirNum should be in [1,num_dependents]");
+
+				/* full Jacobian or Jacobian projection:*/
+				weightsTimesJac=xNew<double>(weightNum,num_independents);
+
+				/*set the forward method function pointers: */
+				#ifdef _HAVE_GSL_
+				anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
+				#endif
+
+				/*allocate the space for the parameters to invoke the EDF fos_reverse :*/
+				anEDF_for_solverx_p->dpp_U=xNew<double>(weightNum,anEDF_for_solverx_p->max_m);
+				anEDF_for_solverx_p->dpp_Z=xNew<double>(weightNum,anEDF_for_solverx_p->max_n);
+
+				/*seed matrix: */
+				weights=xNewZeroInit<double>(weightNum,num_dependents);
+
+				/*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/
+				for (int i=0; i<weightNum; ++i) {
+					/* make sure the index is in range*/
+					if (depIndices[i]>num_dependents) {
+						_error_("depIndices values must be in [0,num_dependents-1]");
+					}
+					if (anIndexSet.find(depIndices[i])!=anIndexSet.end()) {
+						_error_("duplicate depIndices values are not allowed until we implement Jacobian decompression");
+					}
+					anIndexSet.insert(depIndices[i]);
+					/* now populate the seed matrix from the set of independent indices;
+					 * simple setup with a single 1.0 per column and at most a single 1.0 per row*/
+					weights[depIndices[i]][i]=1.0;
+				}
+
+				/*call driver: */
+				fov_reverse(my_rank,num_dependents,num_independents, weightNum, weights, weightsTimesJac );
+
+				/*add to results: */
+				femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,*weightsTimesJac,weightNum*num_independents,1,1,0.0));
+
+				/*Free resources: */
+				xDelete(weights);
+				xDelete(weightsTimesJac);
+				xDelete(depIndices);
+			} /*}}}*/
+			else _error_("driver: " << driver << " not yet supported!");
+
+			/* delete the allocated space for the parameters:*/
+			xDelete(anEDF_for_solverx_p->dp_x);
+			xDelete(anEDF_for_solverx_p->dp_X);
+			xDelete(anEDF_for_solverx_p->dpp_X);
+			xDelete(anEDF_for_solverx_p->dp_y);
+			xDelete(anEDF_for_solverx_p->dp_Y);
+			xDelete(anEDF_for_solverx_p->dpp_Y);
+			xDelete(anEDF_for_solverx_p->dp_U);
+			xDelete(anEDF_for_solverx_p->dpp_U);
+			xDelete(anEDF_for_solverx_p->dp_Z);
+			xDelete(anEDF_for_solverx_p->dpp_Z);
+
+			if(VerboseAutodiff())_printf0_("   end AD core\n");
+
+			/*Free resources: */
+			xDelete(xp);
+			xDelete(axp); 
+			xDelete(driver);
+		#else
+			_error_("Should not be requesting AD drivers when an AD library is not available!");
+		#endif
+	}
+}
Index: issm/trunk-jpl/src/c/cores/adjointbalancethickness_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/adjointbalancethickness_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/adjointbalancethickness_core.cpp	(revision 16518)
@@ -0,0 +1,39 @@
+/*!\file:  adjointbalancethickness_core.cpp
+ * \brief compute inverse method adjoint state
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void adjointbalancethickness_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool save_results;
+
+	/*retrieve parameters:*/
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+	/*compute thickness */
+	if(VerboseSolution()) _printf0_("   computing thickness\n");
+	femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+	solutionsequence_linear(femmodel);
+
+	/*Call SurfaceAreax, because some it might be needed by PVector*/
+	SurfaceAreax(NULL,femmodel);
+
+	/*compute adjoint*/
+	if(VerboseSolution()) _printf0_("   computing adjoint\n");
+	femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum,AdjointBalancethicknessAnalysisEnum);
+	solutionsequence_adjoint_linear(femmodel);
+
+	/*Save results*/
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		int outputs[1] = {AdjointEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
+	}
+}
Index: issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp	(revision 16518)
@@ -0,0 +1,48 @@
+/*!\file:  adjointstressbalance_core.cpp
+ * \brief compute inverse method adjoint state
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void adjointstressbalance_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool isFS;
+	bool save_results;
+	bool conserve_loads   = true;
+
+	/*retrieve parameters:*/
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+	/*Compute velocities*/
+	if(VerboseSolution()) _printf0_("   computing velocities\n");
+	femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+	solutionsequence_nonlinear(femmodel,conserve_loads); 
+
+	/*Call SurfaceAreax, because some it might be needed by PVector*/
+	SurfaceAreax(NULL,femmodel);
+
+	/*Compute adjoint*/
+	if(VerboseSolution()) _printf0_("   computing adjoint\n");
+	femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum,AdjointHorizAnalysisEnum);
+	solutionsequence_adjoint_linear(femmodel);
+
+	/*Save results*/
+	if(save_results || true){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		if (isFS){
+			int outputs[4] = {AdjointxEnum,AdjointyEnum,AdjointzEnum,AdjointpEnum};
+			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
+		}
+		else{
+			int outputs[2] = {AdjointxEnum,AdjointyEnum};
+			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+		}
+	}
+}
Index: issm/trunk-jpl/src/c/cores/balancethickness_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/balancethickness_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/balancethickness_core.cpp	(revision 16518)
@@ -0,0 +1,32 @@
+/*!\file: balancethickness_core.cpp
+ * \brief: core of the balancethickness solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void balancethickness_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool save_results;
+
+	/*activate formulation: */
+	femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+	if(VerboseSolution()) _printf0_("call computational core:\n");
+	solutionsequence_linear(femmodel);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		int outputs = ThicknessEnum;
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+	}
+
+}
Index: issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp	(revision 16518)
@@ -0,0 +1,37 @@
+/*!\file: balancevelocity_core.cpp
+ * \brief: core of the balancevelocity solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void balancevelocity_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool save_results;
+
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+	if(VerboseSolution()) _printf0_("computing smoothed slopes:\n");
+	femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum);
+	solutionsequence_linear(femmodel);
+	femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum);
+	solutionsequence_linear(femmodel);
+	//surfaceslope_core(femmodel);
+
+	if(VerboseSolution()) _printf0_("call computational core:\n");
+	femmodel->SetCurrentConfiguration(BalancevelocityAnalysisEnum);
+	solutionsequence_linear(femmodel);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		int outputs[3] = {SurfaceSlopeXEnum,SurfaceSlopeYEnum,VelEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
+	}
+
+}
Index: issm/trunk-jpl/src/c/cores/bedslope_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/bedslope_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/bedslope_core.cpp	(revision 16518)
@@ -0,0 +1,44 @@
+/*!\file: bedslope_core.cpp
+ * \brief: core of the slope solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void bedslope_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool save_results;
+	int  meshtype;
+
+	/*Recover some parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
+
+	if(VerboseSolution()) _printf0_("   computing slope\n");
+
+	/*Call on core computations: */
+	femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeXAnalysisEnum);
+	solutionsequence_linear(femmodel);
+	if(meshtype!=Mesh2DverticalEnum){
+		femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeYAnalysisEnum);
+		solutionsequence_linear(femmodel);
+	}
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		if(meshtype!=Mesh2DverticalEnum){
+			int outputs[2] = {BedSlopeXEnum,BedSlopeYEnum};
+			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+		}
+		else{
+			int outputs[1] = {BedSlopeXEnum};
+			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
+		}
+	}
+
+}
Index: issm/trunk-jpl/src/c/cores/control_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/control_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/control_core.cpp	(revision 16518)
@@ -0,0 +1,129 @@
+/*!\file: control_core.cpp
+ * \brief: core of the control solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+/*Local prototypes*/
+bool controlconvergence(IssmDouble J, IssmDouble tol_cm);
+
+void control_core(FemModel* femmodel){
+
+	int     i;
+
+	/*parameters: */
+	int        num_controls;
+	int        nsteps;
+	IssmDouble tol_cm;
+	int        solution_type;
+	bool       isFS;
+	bool       dakota_analysis = false;
+
+	int        *control_type   = NULL;
+	IssmDouble *maxiter        = NULL;
+	IssmDouble *cm_jump        = NULL;
+
+	/*intermediary: */
+	IssmDouble search_scalar = 1.;
+	OptArgs    optargs;
+	OptPars    optpars;
+
+	/*Solution and Adjoint core pointer*/
+	void (*solutioncore)(FemModel*) = NULL;
+	void (*adjointcore)(FemModel*)  = NULL;
+
+	/*output: */
+	IssmDouble* J=NULL;
+
+	/*Recover parameters used throughout the solution*/
+	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+	femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
+	femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
+	femmodel->parameters->FindParam(&maxiter,NULL,InversionMaxiterPerStepEnum);
+	femmodel->parameters->FindParam(&cm_jump,NULL,InversionStepThresholdEnum);
+	femmodel->parameters->FindParam(&tol_cm,InversionCostFunctionThresholdEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
+	femmodel->parameters->SetParam(false,SaveResultsEnum);
+
+	/*out of solution_type, figure out solution core and adjoint function pointer*/
+	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+	AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
+
+	/*Launch once a complete solution to set up all inputs*/
+	if(VerboseControl()) _printf0_("   preparing initial solution\n");
+	if(isFS) solutioncore(femmodel);
+
+	/*Initialize cost function: */
+	J=xNew<IssmDouble>(nsteps);
+
+	/*Initialize some of the BrentSearch arguments: */
+	optargs.femmodel=femmodel;
+	optpars.xmin=0; optpars.xmax=1;
+
+	/*Start looping: */
+	for(int n=0;n<nsteps;n++){
+
+		/*Display info*/
+		if(VerboseControl()) _printf0_("\n" << "   control method step " << n+1 << "/" << nsteps << "\n");
+
+
+		/*In steady state inversion, compute new temperature field now*/
+		if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
+
+		if(VerboseControl()) _printf0_("   compute adjoint state:\n");
+		adjointcore(femmodel);
+		gradient_core(femmodel,n,search_scalar==0.);
+
+		if(VerboseControl()) _printf0_("   optimizing along gradient direction\n");
+		optpars.maxiter=reCast<int,IssmDouble>(maxiter[n]); optpars.cm_jump=cm_jump[n];
+		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,&optargs);
+
+		if(VerboseControl()) _printf0_("   updating parameter using optimized search scalar\n"); //true means update save controls
+		InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
+
+		if(controlconvergence(J[n],tol_cm)) break;
+	}
+
+	if(VerboseControl()) _printf0_("   preparing final solution\n");
+	femmodel->parameters->SetParam(true,SaveResultsEnum);
+	solutioncore(femmodel);
+
+	/*some results not computed by steadystate_core or stressbalance_core: */
+	if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run!
+		femmodel->OutputControlsx(&femmodel->results);
+
+		#ifdef _HAVE_ADOLC_
+		IssmPDouble* J_passive=xNew<IssmPDouble>(nsteps);
+		for(i=0;i<nsteps;i++) J_passive[i]=reCast<IssmPDouble>(J[i]);
+		femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,1,0));
+		xDelete<IssmPDouble>(J_passive);
+		#else
+		femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J,nsteps,1,1,0));
+		#endif
+	}
+
+	/*Free ressources: */
+	xDelete<int>(control_type);
+	xDelete<IssmDouble>(maxiter);
+	xDelete<IssmDouble>(cm_jump);
+	xDelete<IssmDouble>(J);
+}
+bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
+
+	bool converged=false;
+
+	/*Has convergence been reached?*/
+	if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
+		converged=true;
+		if(VerboseConvergence()) _printf0_("      Convergence criterion reached: J = " << J << " < " << tol_cm);
+	}
+
+	return converged;
+}
Index: issm/trunk-jpl/src/c/cores/controltao_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/controltao_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/controltao_core.cpp	(revision 16518)
@@ -0,0 +1,179 @@
+/*!\file: control_core.cpp
+ * \brief: core of the control solution 
+ */ 
+#include <config.h>
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+#if defined (_HAVE_TAO_) && defined (_HAVE_PETSC_) && (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ > 1)
+#include <tao.h>
+
+/*Local prototype*/
+int FormFunctionGradient(TaoSolver,Vec,IssmDouble*,Vec,void*);
+int IssmMonitor(TaoSolver,void*);
+typedef struct {
+	FemModel* femmodel;
+	double*   J;
+} AppCtx;
+
+void controltao_core(FemModel* femmodel){
+
+	/*TAO*/
+	int                 ierr;
+	int                 num_controls,solution_type;
+	int                 nsteps,maxiter;
+	AppCtx              user;
+	TaoSolver           tao = 0;
+	int                *control_list = NULL;
+	Vector<IssmDouble> *X            = NULL;
+	Vector<IssmDouble> *G            = NULL;
+	Vector<IssmDouble> *XL           = NULL;
+	Vector<IssmDouble> *XU           = NULL;
+
+	/*Initialize TAO*/
+	int argc; char **args=NULL;
+	PetscGetArgs(&argc,&args);
+	ierr = TaoInitialize(&argc,&args,(char*)0,"");
+	if(ierr) _error_("Could not initialize Tao");
+
+	/*Recover some parameters*/
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+	femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum);
+	femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
+	femmodel->parameters->SetParam(false,SaveResultsEnum);
+	maxiter=nsteps*10;
+
+	/*Prepare Toolkit*/
+	ToolkitsOptionsFromAnalysis(femmodel->parameters,DefaultAnalysisEnum);
+
+	/*Initialize TAO*/
+	TaoCreate(IssmComm::GetComm(),&tao);
+	if(VerboseControl()) _printf0_("   Initializing the Toolkit for Advanced Optimization (TAO)\n");
+	TaoSetFromOptions(tao);
+	TaoSetType(tao,"tao_blmvm");
+	//TaoSetType(tao,"tao_cg");
+	//TaoSetType(tao,"tao_lmvm");
+
+	/*Prepare all TAO parameters*/
+	TaoSetMonitor(tao,IssmMonitor,&user,NULL);
+	TaoSetMaximumFunctionEvaluations(tao,maxiter);
+	TaoSetMaximumIterations(tao,nsteps);
+	TaoSetTolerances(tao,0.,0.,0.,0.,0.);
+
+	GetVectorFromControlInputsx(&X, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
+	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+	TaoSetInitialVector(tao,X->pvector->vector);
+	TaoSetVariableBounds(tao,XL->pvector->vector,XU->pvector->vector);
+	delete XL;
+	delete XU;
+
+	user.J=xNewZeroInit<double>(maxiter+5);
+	user.femmodel=femmodel;
+	TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user); 
+
+	/*Solver optimization problem*/
+	if(VerboseControl()) _printf0_("   Starting optimization\n");
+	TaoSolve(tao);
+	TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);
+
+	/*Save results*/
+	TaoGetSolutionVector(tao,&X->pvector->vector);
+	G=new Vector<IssmDouble>(0); VecFree(&G->pvector->vector);
+	TaoGetGradientVector(tao,&G->pvector->vector);
+	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
+	femmodel->OutputControlsx(&femmodel->results);
+	femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,1,0));
+
+	/*Finalize*/
+	if(VerboseControl()) _printf0_("   preparing final solution\n");
+	femmodel->parameters->SetParam(true,SaveResultsEnum);
+	void (*solutioncore)(FemModel*)=NULL;
+	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+	solutioncore(femmodel);
+
+	/*Clean up and return*/
+	xDelete<int>(control_list);
+	xDelete<double>(user.J);
+	delete X;
+	TaoDestroy(&tao);
+	TaoFinalize();
+}
+int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *userCtx){
+
+	/*Retreive arguments*/
+	int                  solution_type;
+	AppCtx              *user            = (AppCtx *)userCtx;
+	FemModel            *femmodel        = user->femmodel;
+	Vector<IssmDouble>  *gradient        = NULL;
+	Vector<IssmDouble>  *X               = NULL;
+
+	/*Convert input to Vec*/
+	X=new Vector<IssmDouble>(Xpetsc);
+
+	/*Set new variable*/
+	//VecView(X,PETSC_VIEWER_STDOUT_WORLD);
+	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	delete X;
+
+	/*Recover some parameters*/
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+
+	/*Compute solution and adjoint*/
+	void (*solutioncore)(FemModel*)=NULL;
+	void (*adjointcore)(FemModel*)=NULL;
+	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
+	AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
+	solutioncore(femmodel);
+	adjointcore(femmodel);
+
+	/*Compute objective function*/
+	femmodel->CostFunctionx(fcn);
+
+	/*Compute gradient*/
+	Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	VecCopy(gradient->pvector->vector,G); delete gradient;
+	VecScale(G,-1.);
+
+	/*Clean-up and return*/
+	return 0;
+}
+int IssmMonitor(TaoSolver tao, void *userCtx){
+
+	int         its,num_responses;
+	IssmDouble  f,gnorm,cnorm,xdiff;
+	AppCtx     *user      = (AppCtx *)userCtx;
+	FemModel   *femmodel  = user->femmodel;
+	int        *responses = NULL;
+
+	femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
+
+	TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, NULL);
+	if(its==0) _printf0_("Iter       Function      Residual  |  List of contributions\n");
+	if(its==0) _printf0_("___________________________________________________________\n");
+	_printf0_(setw(4)<<its<<"   "<<setw(12)<<setprecision(7)<<f<<"  "<<setw(12)<<setprecision(7)<<gnorm<<"  | ");
+	user->J[its]=f;
+
+	/*Retrieve objective functions independently*/
+	for(int i=0;i<num_responses;i++){
+		femmodel->Responsex(&f,EnumToStringx(responses[i]));
+		_printf0_(" "<<setw(12)<<setprecision(7)<<f);
+	}
+	_printf0_("\n");
+
+	/*Clean-up and return*/
+	xDelete<int>(responses);
+	return 0;
+}
+
+#else
+void controltao_core(FemModel* femmodel){
+	_error_("TAO not installed or PETSc version not supported");
+}
+#endif //_HAVE_TAO_ 
Index: issm/trunk-jpl/src/c/cores/cores.h
===================================================================
--- issm/trunk-jpl/src/c/cores/cores.h	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/cores.h	(revision 16518)
@@ -0,0 +1,60 @@
+/*
+ * cores.h: 
+ */
+
+#ifndef _ANALYSES_H_
+#define _ANALYSES_H_
+
+/*forward declarations: */
+struct OptArgs;
+class FemModel;
+class Parameters;
+template <class doubletype> class Matrix;
+template <class doubletype> class Vector;
+
+#include "../shared/io/Comm/IssmComm.h"
+#include "../shared/Numerics/types.h"
+
+/*cores: */
+void adjointstressbalance_core(FemModel* femmodel);
+void adjointbalancethickness_core(FemModel* femmodel);
+void gradient_core(FemModel* femmodel,int n=0,bool orthogonalize=false);
+void stressbalance_core(FemModel* femmodel);
+void hydrology_core(FemModel* femmodel);
+void thermal_core(FemModel* femmodel);
+void surfaceslope_core(FemModel* femmodel);
+void bedslope_core(FemModel* femmodel);
+void meshdeformation_core(FemModel* femmodel);
+void control_core(FemModel* femmodel);
+void controltao_core(FemModel* femmodel);
+void masstransport_core(FemModel* femmodel);
+void extrudefrombase_core(FemModel* femmodel);
+void extrudefromtop_core(FemModel* femmodel);
+void balancethickness_core(FemModel* femmodel);
+void balancevelocity_core(FemModel* femmodel);
+void slopecompute_core(FemModel* femmodel);
+void steadystate_core(FemModel* femmodel);
+void transient_core(FemModel* femmodel);
+void dakota_core(FemModel* femmodel);
+void ad_core(FemModel* femmodel);
+void dummy_core(FemModel* femmodel);
+void gia_core(FemModel* femmodel);
+void damage_core(FemModel* femmodel);
+IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs);
+
+//optimization
+int GradJSearch(IssmDouble* search_vector,FemModel* femmodel,int step);
+
+//diverse
+void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prootpath,int argc,char **argv);
+void WriteLockFile(char* filename);
+void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);
+void PrintBanner(void);
+
+//solution configuration
+void AnalysisConfiguration(int** pcores,int* pnumcores, int solutiontype);
+void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype);
+void WrapperCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype,bool nodakotacore=false);
+void AdjointCorePointerFromSolutionEnum(void (**padjointcore)(FemModel*),int solutiontype);
+
+#endif
Index: issm/trunk-jpl/src/c/cores/dakota_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/dakota_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/dakota_core.cpp	(revision 16518)
@@ -0,0 +1,230 @@
+/*!\file:  dakota_core.cpp
+ * \brief: wrapper to the Dakota capabilities. qmu fires up Dakota, and registers a Dakota Pluggin
+ * which will be in charge of running the solution sequences repeteadly, to garner statistics. 
+ *
+ * This routine deals with running ISSM and Dakota in library mode. In library mode, Dakota does not 
+ * run as an execuatble. Its capabilities are linked into the ISSM software. ISSM calls dakota routines 
+ * directly from the dakota library. qmu.cpp is the code that is in charge of calling those routines. 
+ *
+ * Dakota has its own way of running in parallel (for embarassingly parallel jobs). We do not want that, 
+ * as ISSM knows exactly how to run "really parallel" jobs that use all CPUS. To bypass Dakota's parallelism, 
+ * we overloaded the constructor for the parallel library (see the Dakota patch in the externalpackages/dakota
+ * directory). This overloaded constructor fires up Dakota serially on CPU 0 only! We take care of broadcasting 
+ * to the other CPUS, hence ISSM is running in parallel, and Dakota serially on CPU0. 
+ *
+ * Now, how does CPU 0 drive all other CPUS to carry out sensitivity analysese? By synchronizing its call to 
+ * our ISSM cores (stressbalance_core, thermal_core, transient_core, etc ...) on CPU 0 with all other CPUS. 
+ * This explains the structure of qmu.cpp, where cpu 0 runs Dakota, the Dakota pluggin fires up DakotaSpawnCore.cpp, 
+ * while the other CPUS are waiting for a broadcast from CPU0, once they get it, they also fire up 
+ * DakotaSpawnCore. In the end, DakotaSpawnCore is fired up on all CPUS, with CPU0 having Dakota inputs, that it will 
+ * broacast to other CPUS. 
+ *
+ * Now, how does dakota call the DakotaSpawnCore routine? The DakotaSpawnCore is embedded into the DakotaPlugin object 
+ * which is derived from the Direct Interface Dakota objct. This is the only way to run Dakota in library 
+ * mode (see their developper guide for more info). Dakota registers the DakotaPlugin object into its own 
+ * database, and calls on the embedded DakotaSpawnCore from CPU0. 
+ *
+ */ 
+
+/*include files: {{{*/
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../shared/shared.h"
+#include "../classes/classes.h"
+
+#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in.
+#if _DAKOTA_MAJOR_ < 5 || (_DAKOTA_MAJOR_ == 5 && _DAKOTA_MINOR_ < 3)
+#include <ParallelLibrary.H>
+#include <ProblemDescDB.H>
+#include <DakotaStrategy.H>
+#include <DakotaModel.H>
+#include <DakotaInterface.H>
+#else
+#include <ParallelLibrary.hpp>
+#include <ProblemDescDB.hpp>
+#include <DakotaStrategy.hpp>
+#include <DakotaModel.hpp>
+#include <DakotaInterface.hpp>
+#endif
+#include "./DakotaSpawnCore.h"
+#endif
+/*}}}*/
+/*DakotaPlugin class {{{*/
+#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in.
+#if _DAKOTA_MAJOR_ < 5 || (_DAKOTA_MAJOR_ == 5 && _DAKOTA_MINOR_ < 3)
+#include <DirectApplicInterface.H>
+#include <DakotaResponse.H>
+#include <ParamResponsePair.H>
+#include <system_defs.h>
+#include <ProblemDescDB.H>
+#include <ParallelLibrary.H>
+#else
+#include <DirectApplicInterface.hpp>
+#include <DakotaResponse.hpp>
+#include <ParamResponsePair.hpp>
+#include <ProblemDescDB.hpp>
+#include <ParallelLibrary.hpp>
+#endif
+namespace SIM {
+	class DakotaPlugin: public Dakota::DirectApplicInterface{
+		public:
+			/*these fields are used by core solutions: */
+			void *femmodel;
+			int   counter;
+			/*Constructors/Destructors*/
+			DakotaPlugin(const Dakota::ProblemDescDB& problem_db,void* in_femmodel):Dakota::DirectApplicInterface(problem_db){/*{{{*/
+				femmodel = in_femmodel;
+				counter  = 0;
+			}/*}}}*/
+			~DakotaPlugin(){/*{{{*/
+				/* Virtual destructor handles referenceCount at Interface level. */ 
+			}/*}}}*/
+		protected:
+			/*execute the input filter portion of a direct evaluation invocation*/
+			//int derived_map_if(const Dakota::String& if_name);
+			/*execute an analysis code portion of a direct evaluation invocation*/
+			int derived_map_ac(const Dakota::String& driver){/*{{{*/
+
+				int i;
+				IssmDouble* variables=NULL;
+				char** variable_descriptors=NULL;
+				char*  variable_descriptor=NULL;
+				IssmDouble* responses=NULL;
+
+				/*increae counter: */
+				counter++;
+
+				/*Before launching analysis, we need to transfer the dakota inputs into Issm 
+				 *readable variables: */
+
+				/*First, the variables: */
+				variables=xNew<IssmDouble>(numACV);
+				for(i=0;i<numACV;i++){
+					variables[i]=xC[i];
+				}
+				/*The descriptors: */
+				variable_descriptors=xNew<char*>(numACV);
+				for(i=0;i<numACV;i++){
+					std::string label=xCLabels[i];
+					variable_descriptor=xNew<char>(strlen(label.c_str())+1);
+					memcpy(variable_descriptor,label.c_str(),(strlen(label.c_str())+1)*sizeof(char));
+
+					variable_descriptors[i]=variable_descriptor;
+				}
+
+				/*Initialize responses: */
+				responses=xNewZeroInit<IssmDouble>(numFns);
+
+				/*run core solution: */
+				DakotaSpawnCore(responses,numFns, variables,variable_descriptors,numACV,femmodel,counter);
+
+				/*populate responses: */
+				for(i=0;i<numFns;i++){
+					fnVals[i]=responses[i];
+				}
+
+				/*Free ressources:*/
+				xDelete<IssmDouble>(variables);
+				for(i=0;i<numACV;i++){
+					variable_descriptor=variable_descriptors[i];
+					xDelete<char>(variable_descriptor);
+				}
+				xDelete<char*>(variable_descriptors);
+				xDelete<IssmDouble>(responses);
+
+				return 0;
+			}/*}}}*/
+			/*execute the output filter portion of a direct evaluation invocation*/
+			//int derived_map_of(const Dakota::String& of_name);
+			/*add for issm: */
+			int GetCounter(){/*{{{*/
+				return counter;
+			}/*}}}*/
+		private:
+	};
+} 
+#endif
+/*}}}*/
+
+void dakota_core(FemModel* femmodel){ 
+
+	#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in.
+
+	int                my_rank;
+	char              *dakota_input_file  = NULL;
+	char              *dakota_output_file = NULL;
+	char              *dakota_error_file  = NULL;
+	Dakota::ModelLIter ml_iter;
+
+	/*Recover dakota_input_file, dakota_output_file and dakota_error_file, in the parameters dataset in parallel */
+	femmodel->parameters->FindParam(&dakota_input_file,QmuInNameEnum);
+	femmodel->parameters->FindParam(&dakota_output_file,QmuOutNameEnum);
+	femmodel->parameters->FindParam(&dakota_error_file,QmuErrNameEnum);
+
+	/*recover my_rank:*/
+	my_rank=IssmComm::GetRank();
+
+	if(my_rank==0){
+
+		// Instantiate/initialize the parallel library and problem description
+		// database objects.
+		char* dakotamode=xNew<char>(strlen("serial")+1);
+		xMemCpy<char>(dakotamode,"serial",strlen("serial")+1);
+		Dakota::ParallelLibrary parallel_lib(dakotamode); //use our own ISSM Dakota library mode constructor, which only fires up Dakota on CPU 0. 
+		Dakota::ProblemDescDB problem_db(parallel_lib); 
+		xDelete<char>(dakotamode);
+
+		// Manage input file parsing, output redirection, and restart processing
+		// without a CommandLineHandler.  This version relies on parsing of an
+		// input file.
+		problem_db.manage_inputs(dakota_input_file);
+		// specify_outputs_restart() is only necessary if specifying non-defaults
+		parallel_lib.specify_outputs_restart(dakota_output_file,dakota_error_file,NULL,NULL);
+
+		// Instantiate the Strategy object (which instantiates all Model and
+		// Iterator objects) using the parsed information in problem_db.
+		Dakota::Strategy selected_strategy(problem_db);
+
+		// convenience function for iterating over models and performing any
+		// interface plug-ins
+		Dakota::ModelList& models = problem_db.model_list();
+
+		for (ml_iter = models.begin(); ml_iter != models.end(); ml_iter++) {
+
+			Dakota::Interface& interface = ml_iter->interface();
+
+			//set DB nodes to the existing Model specification
+			problem_db.set_db_model_nodes(ml_iter->model_id());
+
+			// Serial case: plug in derived Interface object without an analysisComm
+			interface.assign_rep(new SIM::DakotaPlugin(problem_db,(void*)femmodel), false);
+		}
+
+		// Execute the strategy
+		problem_db.lock(); // prevent run-time DB queries
+		selected_strategy.run_strategy();
+
+		//Warn other cpus that we are done running the dakota iterator, by setting the counter to -1:
+		DakotaSpawnCore(NULL,0, NULL,NULL,0,femmodel,-1);
+
+	}
+	else{
+
+		for(;;){
+			if(!DakotaSpawnCore(NULL,0, NULL,NULL,0,femmodel,0))break; //counter came in at -1 on cpu0, bail out.
+		}
+	}
+
+	/*Free ressources:*/
+	xDelete<char>(dakota_input_file);
+	xDelete<char>(dakota_error_file);
+	xDelete<char>(dakota_output_file);
+
+	#endif //#ifdef _HAVE_DAKOTA_
+}
Index: issm/trunk-jpl/src/c/cores/damage_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/damage_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/damage_core.cpp	(revision 16518)
@@ -0,0 +1,39 @@
+/* 
+ * \brief: damgage_core.cpp: core for the damage solution
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void damage_core(FemModel* femmodel){
+
+	/*intermediary*/
+	bool   save_results;
+	bool   dakota_analysis  = false;
+	int    solution_type;
+
+	if(VerboseSolution()) _printf0_("   computing damage\n");
+	
+	//first recover parameters common to all solutions
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+
+	if(dakota_analysis && solution_type!=TransientSolutionEnum){
+		femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
+		ResetConstraintsx(femmodel);
+	}
+
+	femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
+	solutionsequence_linear(femmodel);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		int outputs = DamageDEnum;
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+	}
+}
Index: issm/trunk-jpl/src/c/cores/dummy_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/dummy_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/dummy_core.cpp	(revision 16518)
@@ -0,0 +1,11 @@
+/*!\file: dummy_core.cpp
+ * \brief: dummy core (nothing done)
+ */ 
+
+class FemModel;
+
+void dummy_core(FemModel* femmodel){
+
+	/*We do not do anything*/
+
+}
Index: issm/trunk-jpl/src/c/cores/extrudefrombase_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/extrudefrombase_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/extrudefrombase_core.cpp	(revision 16518)
@@ -0,0 +1,21 @@
+/*!\file: extrudefrombase_core.cpp
+ * \brief: core of the extrusion solution
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../modules/modules.h"
+
+void extrudefrombase_core(FemModel* femmodel){
+
+	if(VerboseSolution()) _printf0_("extruding solution from base...\n");
+
+	/*Call on core computations: */
+	femmodel->SetCurrentConfiguration(ExtrudeFromBaseAnalysisEnum);
+	femmodel->UpdateConstraintsExtrudeFromBasex();
+	solutionsequence_linear(femmodel);
+
+}
Index: issm/trunk-jpl/src/c/cores/extrudefromtop_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/extrudefromtop_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/extrudefromtop_core.cpp	(revision 16518)
@@ -0,0 +1,21 @@
+/*!\file: extrudefromtop_core.cpp
+ * \brief: core of the extrusion solution
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../modules/modules.h"
+
+void extrudefromtop_core(FemModel* femmodel){
+
+	if(VerboseSolution()) _printf0_("extruding solution from top...\n");
+
+	/*Call on core computations: */
+	femmodel->SetCurrentConfiguration(ExtrudeFromTopAnalysisEnum);
+	femmodel->UpdateConstraintsExtrudeFromTopx();
+	solutionsequence_linear(femmodel);
+
+}
Index: issm/trunk-jpl/src/c/cores/gia_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/gia_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/gia_core.cpp	(revision 16518)
@@ -0,0 +1,56 @@
+/*!\file: gia_core.cpp
+ * \brief: core of the GIA solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+void gia_core(FemModel* femmodel){
+
+	Vector<IssmDouble> *wg    = NULL;
+	Vector<IssmDouble> *dwdtg = NULL;
+	IssmDouble          *x    = NULL;
+	IssmDouble          *y    = NULL;
+
+	/*parameters: */
+	bool save_results;
+	int  gsize;
+	int  configuration_type;
+
+	/*Recover some parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+
+	if(VerboseSolution()) _printf0_("   computing GIA\n");
+
+	/*Call on core computations: */
+	femmodel->SetCurrentConfiguration(GiaAnalysisEnum);
+
+	/*Figure out size of g-set deflection vector and allocate solution vector: */
+	gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
+	wg = new Vector<IssmDouble>(gsize);
+	dwdtg = new Vector<IssmDouble>(gsize);
+
+	/*first, recover x and y vectors from vertices: */
+	VertexCoordinatesx(&x,&y,NULL,femmodel->vertices); //no need for z coordinate
+
+	/*call the main module: */
+	femmodel->Deflection(wg,dwdtg,x,y);
+
+	/*assemble vector: */
+	wg->Assemble();
+	dwdtg->Assemble();
+
+	InputUpdateFromVectorx(femmodel,wg,GiaWEnum,VertexSIdEnum);
+	InputUpdateFromVectorx(femmodel,dwdtg,GiadWdtEnum,VertexSIdEnum);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		const int outputs[2] = {GiaWEnum,GiadWdtEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+	}
+
+}
Index: issm/trunk-jpl/src/c/cores/gradient_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/gradient_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/gradient_core.cpp	(revision 16518)
@@ -0,0 +1,48 @@
+/*!\file:  gradient_core.cpp
+ * \brief compute inverse method gradient direction.
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void gradient_core(FemModel* femmodel,int step,bool orthogonalize){ 
+
+	/*Intermediaries*/
+	IssmDouble  norm_inf;
+	IssmDouble *norm_list    = NULL;
+	Vector<IssmDouble>*     new_gradient = NULL;
+	Vector<IssmDouble>*     gradient     = NULL;
+	Vector<IssmDouble>*     old_gradient = NULL;
+
+	/*Compute gradient*/
+	if(VerboseControl()) _printf0_("   compute cost function gradient\n");
+	Gradjx(&gradient,&norm_list,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters);
+
+	if (orthogonalize){
+		if(VerboseControl()) _printf0_("   orthogonalization\n");
+		ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters);
+		Orthx(&new_gradient,gradient,old_gradient); delete old_gradient; delete gradient;
+	}
+	else{ 
+		new_gradient=gradient;
+	}
+
+	/*Check that gradient is clean*/
+	norm_inf=new_gradient->Norm(NORM_INF);
+	if(norm_inf<=0)    _error_("||∂J/∂α||∞ = 0    gradient norm is zero");
+	if(xIsNan<IssmDouble>(norm_inf))_error_("||∂J/∂α||∞ = NaN  gradient norm is NaN");
+
+	/*plug back into inputs: */
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,new_gradient);
+	delete new_gradient;
+
+	/*Scale Gradients*/
+	ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,norm_list,step);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(norm_list);
+}
Index: issm/trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 16518)
@@ -0,0 +1,105 @@
+/*!\file: hydrology_core.cpp
+ * \brief: core of the hydrology solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void hydrology_core(FemModel* femmodel){
+
+	int i;
+
+	/*intermediary*/
+	int        step,nsteps;
+	int        output_frequency,hydrology_model;
+	bool       save_results;
+	bool       modify_loads=true;
+	bool       isefficientlayer;
+	IssmDouble starttime,final_time;
+	IssmDouble time,dt;
+
+	/*first recover parameters common to all solutions*/
+	femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+	femmodel->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
+	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
+	femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
+
+	/*first compute slopes: */
+	if (hydrology_model==HydrologyshreveEnum){
+		surfaceslope_core(femmodel);
+		bedslope_core(femmodel);
+	}
+
+	/*Compute number of time steps: */
+	if((dt==0)|| (final_time==0)){
+		dt=0;
+		nsteps=1;
+	}
+	else nsteps=reCast<int,IssmDouble>((final_time-starttime)/dt);
+
+	/*initialize: */
+	step=0;
+	time=starttime;
+
+	/*Loop through time: */
+	for(i=0;i<nsteps;i++){
+
+		if(nsteps)if(VerboseSolution()) _printf0_("time step:" << i+1 << "/" << nsteps << "\n");
+		time+=dt;
+		step+=1;
+		femmodel->parameters->SetParam(time,TimeEnum);
+		femmodel->parameters->SetParam(step,StepEnum);
+
+		if (hydrology_model==HydrologyshreveEnum){
+			if(VerboseSolution()) _printf0_("   computing water column\n");
+			femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
+			solutionsequence_nonlinear(femmodel,modify_loads);
+
+			/*transfer water column thickness to old water column thickness: */
+
+			InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
+
+			if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
+				if(VerboseSolution()) _printf0_("   saving results \n");
+				int outputs[3] = {WatercolumnEnum,HydrologyWaterVxEnum,HydrologyWaterVyEnum};
+				femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
+
+				/*unload results*/
+				if(VerboseSolution()) _printf0_("   saving temporary results\n");
+				OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
+			}
+		}
+
+		else if (hydrology_model==HydrologydcEnum){
+			InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum);
+			femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
+			if (isefficientlayer){
+				InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
+			}
+
+			if(VerboseSolution()) _printf0_("   computing water head\n");
+			solutionsequence_hydro_nonlinear(femmodel);
+			if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
+				if(VerboseSolution()) _printf0_("   saving results \n");
+				if(isefficientlayer){
+					int outputs[4] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum};
+					femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
+				}
+				else{
+					int outputs[2] = {SedimentHeadEnum,SedimentHeadResidualEnum};
+					femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+				}
+				/*unload results*/
+				if(VerboseSolution()) _printf0_("   saving temporary results\n");
+				OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
+			}
+		}
+
+	}
+}
Index: issm/trunk-jpl/src/c/cores/masstransport_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/masstransport_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/masstransport_core.cpp	(revision 16518)
@@ -0,0 +1,87 @@
+/*!\file: masstransport_core.cpp
+ * \brief: core of the masstransport solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void masstransport_core(FemModel* femmodel){
+
+	/*parameters: */
+	int    i;
+	int    numoutputs,meshtype;
+	bool   save_results;
+	bool   issmbgradients,ispdd,isdelta18o,isFS,isfreesurface,dakota_analysis;
+	int    solution_type;
+	char** requested_outputs = NULL;
+
+	/*activate configuration*/
+	femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum);
+
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&issmbgradients,SurfaceforcingsIssmbgradientsEnum);
+	femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum);
+	femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum);
+	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
+	femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
+	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
+
+	/*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/
+	if(dakota_analysis && solution_type==MasstransportSolutionEnum){
+		InputDuplicatex(femmodel,QmuSurfaceEnum,SurfaceEnum);
+		InputDuplicatex(femmodel,QmuThicknessEnum,ThicknessEnum);
+		InputDuplicatex(femmodel,QmuBedEnum,BedEnum);
+	}
+
+	if(issmbgradients){
+	  if(VerboseSolution())_printf_("	call smb gradients module\n\n");
+	  SmbGradientsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	}
+	if(ispdd){
+		if(isdelta18o){
+			if(VerboseSolution()) _printf0_("   call Delta18oParametrization module\n");
+			Delta18oParameterizationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		} 
+		if(VerboseSolution()) _printf0_("   call positive degree day module\n");
+		PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	}
+
+	if(isFS && isfreesurface){
+		if(VerboseSolution()) _printf0_("   call free surface computational core\n");
+		femmodel->SetCurrentConfiguration(FreeSurfaceBaseAnalysisEnum);
+		solutionsequence_linear(femmodel);
+		if(meshtype==Mesh2DverticalEnum){
+			femmodel->parameters->SetParam(BedEnum,InputToExtrudeEnum);
+			extrudefrombase_core(femmodel);
+		}
+		femmodel->SetCurrentConfiguration(FreeSurfaceTopAnalysisEnum);
+		solutionsequence_linear(femmodel);
+		if(meshtype==Mesh2DverticalEnum){
+			femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
+			extrudefromtop_core(femmodel);
+		}
+	}
+	else{
+		if(VerboseSolution()) _printf0_("   call computational core\n");
+		solutionsequence_linear(femmodel);
+	}
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+	}
+
+	if(solution_type==MasstransportSolutionEnum)femmodel->RequestedDependentsx();
+
+	/*Free ressources:*/
+	if(numoutputs){for(int i=0;i<numoutputs;i++){char* string=requested_outputs[i];xDelete<char>(string);} xDelete<char*>(requested_outputs);} 
+}
Index: issm/trunk-jpl/src/c/cores/meshdeformation_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/meshdeformation_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/meshdeformation_core.cpp	(revision 16518)
@@ -0,0 +1,30 @@
+/*!\file: meshdeformation_core.cpp
+ * \brief: core of the slope solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../modules/modules.h"
+
+void meshdeformation_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool save_results;
+
+	/*Recover some parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+	if(VerboseSolution()) _printf0_("computing mesh deformation (elasticity)...\n");
+
+	/*Call on core computations: */
+	femmodel->SetCurrentConfiguration(MeshdeformationAnalysisEnum);
+	solutionsequence_linear(femmodel);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("saving results:\n");
+	}
+
+}
Index: issm/trunk-jpl/src/c/cores/objectivefunction.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/objectivefunction.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/objectivefunction.cpp	(revision 16518)
@@ -0,0 +1,77 @@
+/*!\file:  objectivefunction
+ * \brief  objective function that returns a misfit, for a certain parameter.
+ */ 
+
+/*include files: {{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../modules/modules.h"
+/*}}}*/
+
+IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs){
+
+	/*output: */
+	IssmDouble J;
+
+	/*parameters: */
+	int        solution_type,analysis_type;
+	bool       isFS       = false;
+	bool       conserve_loads = true;
+	FemModel  *femmodel       = NULL;
+
+	/*Recover finite element model: */
+	femmodel=optargs->femmodel;
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+
+	/*set analysis type to compute velocity: */
+	if (solution_type==SteadystateSolutionEnum || solution_type==StressbalanceSolutionEnum){
+		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+	}
+	else if (solution_type==BalancethicknessSolutionEnum){
+		femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+	}
+	else if (solution_type==BalancethicknessSoftSolutionEnum){
+		femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+	}
+	else{
+		_error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
+	}
+
+	/*update parameter according to scalar: */ //false means: do not save control
+	InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false);
+
+	/*Run stressbalance with updated inputs: */
+	if (solution_type==SteadystateSolutionEnum){
+		stressbalance_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
+	}
+	else if (solution_type==StressbalanceSolutionEnum){
+		solutionsequence_nonlinear(femmodel,conserve_loads); 
+	}
+	else if (solution_type==BalancethicknessSolutionEnum){
+		solutionsequence_linear(femmodel); 
+	}
+	else if (solution_type==BalancethicknessSoftSolutionEnum){
+		/*Don't do anything*/
+	}
+	else{
+		_error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
+	}
+
+	/*Compute misfit for this velocity field.*/
+	femmodel->CostFunctionx(&J);
+
+	/*Free ressources:*/
+	return J;
+}
Index: issm/trunk-jpl/src/c/cores/steadystate_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/steadystate_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/steadystate_core.cpp	(revision 16518)
@@ -0,0 +1,130 @@
+/*!\file: steadystate_core.cpp
+ * \brief: core of the steadystate solution 
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+/*Local prototypes*/
+bool steadystateconvergence(Vector<IssmDouble>* tg,Vector<IssmDouble>* tg_old,Vector<IssmDouble>* ug,Vector<IssmDouble>* ug_old,IssmDouble reltol);
+
+void steadystate_core(FemModel* femmodel){
+
+	/*intermediary: */
+	int i;
+	int step; 
+	Vector<IssmDouble>* ug     = NULL;
+	Vector<IssmDouble>* ug_old = NULL;
+	Vector<IssmDouble>* tg     = NULL;
+	Vector<IssmDouble>* tg_old = NULL;
+
+	/*parameters: */
+	bool        save_results,isenthalpy;
+	int         maxiter;
+	IssmDouble  reltol;
+	int         numoutputs        = 0;
+	char** requested_outputs = NULL;
+
+	/* recover parameters:*/
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&maxiter,SteadystateMaxiterEnum);
+	femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
+	femmodel->parameters->FindParam(&reltol,SteadystateReltolEnum);
+	femmodel->parameters->SetParam(false,SaveResultsEnum);
+	femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum);
+	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum);
+
+	/*intialize counters: */
+	step=1;
+
+	for(;;){
+
+		if(VerboseSolution()) _printf0_("   computing temperature and velocity for step: " << step << "\n");
+		#ifdef _HAVE_THERMAL_
+		thermal_core(femmodel);
+		if(!isenthalpy)femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);/*Could be MeltingAnalysis...*/
+		GetSolutionFromInputsx(&tg,femmodel);
+		#else
+		_error_("ISSM was not compiled with thermal capabilities. Exiting");
+		#endif
+
+		if(VerboseSolution()) _printf0_("   computing new velocity\n");
+		stressbalance_core(femmodel);
+		GetSolutionFromInputsx(&ug,femmodel);
+
+		if(step>1){
+			if(VerboseSolution()) _printf0_("   checking steadystate convergence\n");
+			if(steadystateconvergence(tg,tg_old,ug,ug_old,reltol)) break;
+		}
+		if(step>maxiter){
+			if(VerboseSolution()) _printf0_("   maximum number steadystate iterations " << maxiter << " reached\n");
+			break;
+		}
+
+		/*update results and increase counter*/
+		delete tg_old;tg_old=tg;
+		delete ug_old;ug_old=ug;
+		step++;
+	}
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+	}
+
+	/*Free ressources:*/
+	delete tg_old;
+	delete ug_old;
+	delete tg;
+	delete ug;	
+	if(numoutputs){ for (i=0;i<numoutputs;i++){char* string=requested_outputs[i];xDelete<char>(string);} xDelete<char*>(requested_outputs);}
+}
+bool steadystateconvergence(Vector<IssmDouble>* tg,Vector<IssmDouble>* tg_old,Vector<IssmDouble>* ug,Vector<IssmDouble>* ug_old,IssmDouble reltol){
+
+	/*Output*/
+	bool converged = true;
+
+	/*Intermediary*/
+	Vector<IssmDouble>* dug    = NULL;
+	Vector<IssmDouble>* dtg    = NULL;
+	IssmDouble          ndt,nt;
+	IssmDouble          ndu,nu;
+
+	/*compute norm(du)/norm(u)*/
+	dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0);
+	ndu=dug->Norm(NORM_TWO); nu=ug_old->Norm(NORM_TWO);
+	if (xIsNan<IssmDouble>(ndu) || xIsNan<IssmDouble>(nu)) _error_("convergence criterion is NaN!");
+	if((ndu/nu)<reltol){
+		if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<ndu/nu*100<<" < "<<reltol*100<<" %\n");
+	}
+	else{ 
+		if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<ndu/nu*100<<" > "<<reltol*100<<" %\n");
+		converged=false;
+	}
+
+	/*compute norm(dt)/norm(t)*/
+	dtg=tg_old->Duplicate(); tg_old->Copy(dtg); dtg->AYPX(tg,-1.0);
+	ndt=dtg->Norm(NORM_TWO); nt=tg_old->Norm(NORM_TWO);
+	if (xIsNan<IssmDouble>(ndt) || xIsNan<IssmDouble>(nt)) _error_("convergence criterion is NaN!");
+	if((ndt/nt)<reltol){
+		if(VerboseConvergence()) _printf0_(setw(50)<<left<<"   Temperature convergence: norm(dt)/norm(t)"<<ndt/nt*100<<" < "<<reltol*100<<" %\n");
+	}
+	else{ 
+		if(VerboseConvergence()) _printf0_(setw(50)<<left<<"   Temperature convergence: norm(dt)/norm(t)"<<ndt/nt*100<<" > "<<reltol*100<<" %\n");
+		converged=false;
+	}
+
+	/*clean up and return*/
+	delete dtg;
+	delete dug;
+	return converged;
+}
Index: issm/trunk-jpl/src/c/cores/stressbalance_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/stressbalance_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/stressbalance_core.cpp	(revision 16518)
@@ -0,0 +1,97 @@
+/*!\file: stressbalance_core.cpp
+ * \brief: core of the stressbalance solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void stressbalance_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool  dakota_analysis;
+	int   meshtype;
+	bool  isSIA,isSSA,isL1L2,isHO,isFS;
+	bool  conserve_loads    = true;
+	bool  save_results;
+	int   newton;
+	int   solution_type;
+	int   numoutputs        = 0;
+	char** requested_outputs = NULL;
+	int    i;
+
+
+	/* recover parameters:*/
+	femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
+	femmodel->parameters->FindParam(&isSIA,FlowequationIsSIAEnum);
+	femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
+	femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum);
+	femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
+	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&numoutputs,StressbalanceNumRequestedOutputsEnum);
+	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum);
+
+	/*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/
+	if(dakota_analysis && solution_type==StressbalanceSolutionEnum){
+		InputDuplicatex(femmodel,QmuVxEnum,VxEnum);
+		InputDuplicatex(femmodel,QmuVyEnum,VyEnum);
+		InputDuplicatex(femmodel,QmuVzEnum,VzEnum);
+		if(isFS) InputDuplicatex(femmodel,QmuPressureEnum,PressureEnum);
+	}
+
+	/*Compute slopes: */
+	if(isSIA) surfaceslope_core(femmodel);
+	if(isFS){
+		bedslope_core(femmodel);
+		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+		ResetCoordinateSystemx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	}
+
+	if(isSIA){
+		if(VerboseSolution()) _printf0_("   computing SIA velocities\n");
+
+		/*Take the last velocity into account so that the velocity on the SSA domain is not zero*/
+		if(isSSA || isL1L2 || isHO ) ResetBoundaryConditions(femmodel,StressbalanceSIAAnalysisEnum);
+		femmodel->SetCurrentConfiguration(StressbalanceSIAAnalysisEnum);
+		solutionsequence_linear(femmodel);
+		if(isSSA || isL1L2 || isHO) ResetBoundaryConditions(femmodel,StressbalanceAnalysisEnum);
+	}
+
+	if ((isSSA || isHO || isL1L2) ^ isFS){ // ^ = xor
+		if(VerboseSolution()) _printf0_("   computing velocities\n");
+
+		femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+		if(newton>0)
+		 solutionsequence_newton(femmodel);
+		else
+		 solutionsequence_nonlinear(femmodel,conserve_loads); 
+	}
+
+	if ((isSSA || isL1L2 || isHO) && isFS){
+		if(VerboseSolution()) _printf0_("   computing coupling betweem lower order models and full-FS\n");
+		solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads);
+	}
+
+	if (meshtype==Mesh3DEnum && (isSIA || isSSA || isL1L2 || isHO)){
+		if(VerboseSolution()) _printf0_("   computing vertical velocities\n");
+		femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum);
+		solutionsequence_linear(femmodel);
+	}
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+	}
+
+	if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx();
+
+	/*Free ressources:*/	
+	if(numoutputs){ for (i=0;i<numoutputs;i++){char* string=requested_outputs[i];xDelete<char>(string);} xDelete<char*>(requested_outputs);}
+}
Index: issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp	(revision 16518)
@@ -0,0 +1,43 @@
+/*!\file: surfaceslope_core.cpp
+ * \brief: core of the slope solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../solutionsequences/solutionsequences.h"
+#include "../modules/modules.h"
+
+void surfaceslope_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool save_results;
+	int  meshtype;
+
+	/*Recover some parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
+
+	if(VerboseSolution()) _printf0_("computing slope...\n");
+
+	/*Call on core computations: */
+	femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum);
+	solutionsequence_linear(femmodel);
+	femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum);
+	solutionsequence_linear(femmodel);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("saving results:\n");
+		if(meshtype!=Mesh2DverticalEnum){
+			int outputs[2] = {SurfaceSlopeXEnum,SurfaceSlopeYEnum};
+			femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+
+		}
+		else{
+			int outputs = SurfaceSlopeXEnum;
+			femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
+		}
+	}
+
+}
Index: issm/trunk-jpl/src/c/cores/thermal_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/thermal_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/thermal_core.cpp	(revision 16518)
@@ -0,0 +1,65 @@
+/*!\file: thermal_core.cpp
+ * \brief: core of the thermal solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void thermal_core(FemModel* femmodel){
+
+	/*intermediary*/
+	bool   save_results,isenthalpy;
+	bool   dakota_analysis;
+	int    solution_type,numoutputs;
+	char** requested_outputs = NULL;
+
+	/*first recover parameters common to all solutions*/
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
+	femmodel->parameters->FindParam(&numoutputs,ThermalNumRequestedOutputsEnum);
+	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,ThermalRequestedOutputsEnum);
+
+	if(dakota_analysis && solution_type!=TransientSolutionEnum){
+		InputDuplicatex(femmodel,QmuVxMeshEnum,VxMeshEnum);
+		InputDuplicatex(femmodel,QmuVyMeshEnum,VyMeshEnum);
+		InputDuplicatex(femmodel,QmuVzMeshEnum,VzMeshEnum);
+		InputDuplicatex(femmodel,QmuTemperatureEnum,TemperatureEnum);
+		InputDuplicatex(femmodel,QmuMeltingEnum,BasalforcingsMeltingRateEnum);
+		InputDuplicatex(femmodel,QmuMaterialsRheologyBEnum,MaterialsRheologyBEnum);
+		femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
+		ResetConstraintsx(femmodel);
+	}
+
+	if(isenthalpy){
+
+		if(VerboseSolution()) _printf0_("   computing enthalpy\n");
+		femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
+		solutionsequence_nonlinear(femmodel,true);
+
+		/*transfer enthalpy to enthalpy picard for the next step: */
+		InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum);
+
+		/*Post process*/
+		if(solution_type!=SteadystateSolutionEnum) PostprocessingEnthalpyx(femmodel);
+	}
+	else{
+		if(VerboseSolution()) _printf0_("   computing temperatures\n");
+		femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
+		solutionsequence_thermal_nonlinear(femmodel);
+
+		if(VerboseSolution()) _printf0_("   computing melting\n");
+		femmodel->SetCurrentConfiguration(MeltingAnalysisEnum);
+		solutionsequence_linear(femmodel);
+	}
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+	}
+}
Index: issm/trunk-jpl/src/c/cores/transient_core.cpp
===================================================================
--- issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/cores/transient_core.cpp	(revision 16518)
@@ -0,0 +1,182 @@
+/*!\file: transient_3d_core.cpp
+ * \brief: core of the transient_3d solution 
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <float.h>
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void transient_core(FemModel* femmodel){
+
+	/*parameters: */
+	int    i;
+	IssmDouble starttime,finaltime,dt,yts;
+	bool   isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isdelta18o,isgia;
+	bool   save_results,dakota_analysis;
+	bool   time_adapt=false;
+	int    output_frequency;
+	int    meshtype,groundingline_migration;
+	int    numoutputs         = 0;
+	char** requested_outputs = NULL;
+
+
+	/*intermediary: */
+	int    step;
+	IssmDouble time;
+
+	//first recover parameters common to all solutions
+	femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
+	femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+	femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
+	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
+	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
+	femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
+	femmodel->parameters->FindParam(&time_adapt,TimesteppingTimeAdaptEnum);
+	femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum);
+	femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
+	femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
+	femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
+	femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
+	femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
+	if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
+	femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
+
+	/*initialize: */
+	step=0;
+	time=starttime;
+
+	/*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/
+	if(dakota_analysis){
+		if(isstressbalance){
+			InputDuplicatex(femmodel,QmuVxEnum,VxEnum);
+			InputDuplicatex(femmodel,QmuVyEnum,VyEnum);
+			if(meshtype==Mesh3DEnum){
+				InputDuplicatex(femmodel,QmuVzEnum,VzEnum);
+				if(isFS)InputDuplicatex(femmodel,QmuPressureEnum,PressureEnum);
+			}
+		}
+		if(ismasstransport || isgroundingline){
+			InputDuplicatex(femmodel,QmuThicknessEnum,ThicknessEnum);
+			InputDuplicatex(femmodel,QmuSurfaceEnum,SurfaceEnum);
+			InputDuplicatex(femmodel,QmuBedEnum,BedEnum);
+			InputDuplicatex(femmodel,QmuMaskIceLevelsetEnum,MaskIceLevelsetEnum);
+		}
+		if(isgroundingline) InputDuplicatex(femmodel,QmuMaskGroundediceLevelsetEnum,MaskGroundediceLevelsetEnum);
+		if(isthermal && meshtype==Mesh3DEnum){
+			//Update Vertex Position after updating Thickness and Bed
+			femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum);
+			femmodel->UpdateVertexPositionsx();
+			InputDuplicatex(femmodel,QmuVxMeshEnum,VxMeshEnum);
+			InputDuplicatex(femmodel,QmuVyMeshEnum,VyMeshEnum);
+			InputDuplicatex(femmodel,QmuVzMeshEnum,VzMeshEnum);
+			InputDuplicatex(femmodel,QmuTemperatureEnum,TemperatureEnum);
+			InputDuplicatex(femmodel,QmuMeltingEnum,BasalforcingsMeltingRateEnum);
+			InputDuplicatex(femmodel,QmuMaterialsRheologyBEnum,MaterialsRheologyBEnum);
+			//Reset Thermal Constraints
+			femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
+			ResetConstraintsx(femmodel);
+		}
+	}
+
+	while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.
+
+		/*Increment*/
+		if(time_adapt){
+			femmodel->TimeAdaptx(&dt);
+			if(time+dt>finaltime) dt=finaltime-time;
+			femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
+		}
+		step+=1;
+		time+=dt;
+		femmodel->parameters->SetParam(time,TimeEnum);
+		femmodel->parameters->SetParam(step,StepEnum);
+
+		if(VerboseSolution()) _printf0_("iteration " << step << "/" << floor((finaltime-time)/dt)+step << "  time [yr]: " << time/yts << " (time step: " << dt/yts << ")\n");
+		if(step%output_frequency==0 || time==finaltime)
+		 save_results=true;
+		else
+		 save_results=false;
+		femmodel->parameters->SetParam(save_results,SaveResultsEnum);
+
+		if(isthermal && meshtype==Mesh3DEnum){
+			if(VerboseSolution()) _printf0_("   computing temperatures\n");
+			#ifdef _HAVE_THERMAL_
+			thermal_core(femmodel);
+			#else
+			_error_("ISSM was not compiled with thermal capabilities. Exiting");
+			#endif
+		}
+
+		if(isstressbalance){
+			if(VerboseSolution()) _printf0_("   computing new velocity\n");
+			#ifdef _HAVE_STRESSBALANCE_
+			stressbalance_core(femmodel);
+			#else
+			_error_("ISSM was not compiled with stressbalance capabilities. Exiting");
+			#endif
+		}
+
+		if(ismasstransport){
+			if(VerboseSolution()) _printf0_("   computing new thickness\n");
+			masstransport_core(femmodel);
+			if(VerboseSolution()) _printf0_("   updating vertices positions\n");
+			femmodel->UpdateVertexPositionsx();
+		}
+
+		if(isgroundingline){
+			if(VerboseSolution()) _printf0_("   computing new grounding line position\n");
+			#ifdef _HAVE_GROUNDINGLINE_
+			GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			#else
+			_error_("ISSM was not compiled with grounding line migration capabilities. Exiting");
+			#endif
+			if(save_results){
+				int outputs[3] = {SurfaceEnum,BedEnum,MaskGroundediceLevelsetEnum};
+				femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
+			}
+		}
+		if(isgia){
+			if(VerboseSolution()) _printf0_("   computing glacial isostatic adjustment\n");
+			#ifdef _HAVE_GIA_
+			gia_core(femmodel);
+			#else
+			_error_("ISSM was not compiled with gia capabilities. Exiting");
+			#endif
+
+		}
+
+		/*unload results*/
+		if(save_results){
+			if(VerboseSolution()) _printf0_("   saving transient results\n");
+			femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+			if(isdelta18o){
+				int outputs[2] = {SurfaceforcingsMonthlytemperaturesEnum,SurfaceforcingsPrecipitationEnum};
+				femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+			}
+			if(isgroundingline && (groundingline_migration==SubelementMigrationEnum || groundingline_migration==SubelementMigration2Enum)){
+				int outputs[1] = {MaskGroundediceLevelsetEnum};
+				femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
+			}
+			if(VerboseSolution()) _printf0_("   saving temporary results\n");
+			OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
+		}
+	}
+
+	femmodel->RequestedDependentsx();
+
+	/*Free ressources:*/	
+	if(numoutputs){ for (i=0;i<numoutputs;i++){char* string=requested_outputs[i];xDelete<char>(string);} xDelete<char*>(requested_outputs);}
+}
Index: issm/trunk-jpl/src/c/main/EnvironmentFinalize.cpp
===================================================================
--- issm/trunk-jpl/src/c/main/EnvironmentFinalize.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/main/EnvironmentFinalize.cpp	(revision 16518)
@@ -0,0 +1,25 @@
+/*!\file:  EnvironmentFinalize.cpp
+ * \brief: finalize Petsc, MPI, you name it
+ */ 
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include "../toolkits/toolkits.h"
+#include "../shared/shared.h"
+
+void EnvironmentFinalize(void){
+
+	int my_rank;
+
+	/*Make sure we are all here*/
+	ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD);
+
+	/*Print closing statement*/
+	ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&my_rank);
+	if(!my_rank) printf("closing MPI\n");
+
+	/*Finalize: */
+	ISSM_MPI_Finalize();
+}
Index: issm/trunk-jpl/src/c/main/EnvironmentInit.cpp
===================================================================
--- issm/trunk-jpl/src/c/main/EnvironmentInit.cpp	(revision 16518)
+++ issm/trunk-jpl/src/c/main/EnvironmentInit.cpp	(revision 16518)
@@ -0,0 +1,35 @@
+/*!\file:  EnvironmentInit.cpp
+ * \brief: initialize Petsc, MPI, you name it
+ */ 
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+#include <stdio.h>
+#include "../toolkits/toolkits.h"
+
+ISSM_MPI_Comm EnvironmentInit(int argc,char** argv){
+
+	/*Output*/
+	ISSM_MPI_Comm comm = 0;
+
+	/*Initialize MPI environment: */
+	#if defined(_HAVE_MPI_)
+	ISSM_MPI_Init(&argc,&argv);
+	comm = ISSM_MPI_COMM_WORLD;
+	#else
+	comm = 1; //bogus number for comm, which does not exist anyway.
+	#endif
+
+	/*Print Banner*/
+	int my_rank = 0;
+	ISSM_MPI_Comm_rank(comm,&my_rank);
+	if(!my_rank) printf("\n");
+	if(!my_rank) printf("%s version  %s\n",PACKAGE_NAME,PACKAGE_VERSION);
+	if(!my_rank) printf("(website: %s contact: %s)\n",PACKAGE_URL,PACKAGE_BUGREPORT);
+	if(!my_rank) printf("\n");
+
+	/*Return communicator*/
+	return comm;
+}
Index: issm/trunk-jpl/src/c/main/issm.h
===================================================================
--- issm/trunk-jpl/src/c/main/issm.h	(revision 16517)
+++ issm/trunk-jpl/src/c/main/issm.h	(revision 16518)
@@ -16,7 +16,11 @@
 #include "../classes/classes.h"
 #include "../toolkits/toolkits.h"
-#include "../analyses/analyses.h"
+#include "../cores/cores.h"
 #include "../solutionsequences/solutionsequences.h"
 #include "../modules/modules.h"
 
+/*Environment*/
+ISSM_MPI_Comm EnvironmentInit(int argc,char** argv);
+void EnvironmentFinalize(void);
+
 #endif //ifndef _ISSM_H_
