Index: issm/trunk/src/c/DataSet/DataSet.cpp
===================================================================
--- issm/trunk/src/c/DataSet/DataSet.cpp	(revision 464)
+++ issm/trunk/src/c/DataSet/DataSet.cpp	(revision 465)
@@ -900,5 +900,5 @@
 #undef __FUNCT__
 #define __FUNCT__ "DataSet::CreateKMatrix"
-void  DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+void  DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	vector<Object*>::iterator object;
@@ -911,10 +911,10 @@
 
 			element=(Element*)(*object);
-			element->CreateKMatrix(Kgg,inputs,analysis_type);
+			element->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
 		}
 		if(EnumIsLoad((*object)->Enum())){
 
 			load=(Load*)(*object);
-			load->CreateKMatrix(Kgg,inputs,analysis_type);
+			load->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
 		}
 	}
@@ -924,5 +924,5 @@
 #undef __FUNCT__
 #define __FUNCT__ "DataSet::CreatePVector"
-void  DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type){
+void  DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	vector<Object*>::iterator object;
@@ -935,10 +935,10 @@
 
 			element=(Element*)(*object);
-			element->CreatePVector(pg,inputs,analysis_type);
+			element->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
 		}
 		if(EnumIsLoad((*object)->Enum())){
 
 			load=(Load*)(*object);
-			load->CreatePVector(pg,inputs,analysis_type);
+			load->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
 		}		
 	}
@@ -984,5 +984,5 @@
 #undef __FUNCT__
 #define __FUNCT__ "DataSet::PenaltyCreateKMatrix"
-void  DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
+void  DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
 
 	vector<Object*>::iterator object;
@@ -994,5 +994,5 @@
 
 			load=(Load*)(*object);
-			load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type);
+			load->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type);
 		}
 	}
@@ -1002,5 +1002,5 @@
 #undef __FUNCT__
 #define __FUNCT__ "DataSet::PenaltyCreatePVector"
-void  DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
+void  DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
 
 	vector<Object*>::iterator object;
@@ -1012,5 +1012,5 @@
 
 			load=(Load*)(*object);
-			load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type);
+			load->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type);
 		}		
 	}
@@ -1028,5 +1028,5 @@
 #undef __FUNCT__
 #define __FUNCT__ "RiftConstraints"
-void  DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){
+void  DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
 	
 	throw ErrorException(__FUNCT__," not implemented yet!");
@@ -1045,5 +1045,5 @@
 #undef __FUNCT__
 #define __FUNCT__ "MeltingConstraints"
-void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){
+void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){
 
 	throw ErrorException(__FUNCT__," not implemented yet!");
@@ -1084,5 +1084,5 @@
 
 
-void  DataSet::Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type){
+void  DataSet::Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){
 
 
@@ -1095,12 +1095,12 @@
 
 			element=(Element*)(*object);
-			element->Du(du_g,u_g,u_g_obs,inputs,analysis_type);
-		}
-	}
-
-
-}
-
-void  DataSet::Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type){
+			element->Du(du_g,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
+		}
+	}
+
+
+}
+
+void  DataSet::Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
 
 
@@ -1113,5 +1113,5 @@
 
 			element=(Element*)(*object);
-			element->Gradj(grad_g,u_g,lambda_g,inputs,analysis_type,control_type);
+			element->Gradj(grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type);
 		}
 	}
@@ -1120,5 +1120,5 @@
 }		
 		
-void  DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type){
+void  DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){
 
 	double J=0;;
@@ -1132,5 +1132,5 @@
 
 			element=(Element*)(*object);
-			J+=element->Misfit(u_g,u_g_obs,inputs,analysis_type);
+			J+=element->Misfit(u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
 
 		}
Index: issm/trunk/src/c/DataSet/DataSet.h
===================================================================
--- issm/trunk/src/c/DataSet/DataSet.h	(revision 464)
+++ issm/trunk/src/c/DataSet/DataSet.h	(revision 465)
@@ -47,5 +47,5 @@
 		void  DistributeDofs(int numberofnodes,int numdofspernode);
 		void  CreatePartitioningVector(Vec* ppartition,int numnods,int numdofspernode);
-		void  DistributeNumDofs(int** pnumdofspernode,int numberofnodes,int analysis_type);
+		void  DistributeNumDofs(int** pnumdofspernode,int numberofnodes,int analysis_type,int sub_analysis_type);
 		void  FlagClones(int numberofnodes);
 		int   NumberOfDofs();
@@ -61,17 +61,17 @@
 		void  SetSorting(int* in_sorted_ids,int* in_id_offsets);
 		void  Sort();
-		void  CreateKMatrix(Mat Kgg,void* inputs, int analysis_type);
-		void  CreatePVector(Vec pg,void* inputs, int analysis_type);
+		void  CreateKMatrix(Mat Kgg,void* inputs, int analysis_type,int sub_analysis_type);
+		void  CreatePVector(Vec pg,void* inputs, int analysis_type,int sub_analysis_type);
 		void  UpdateFromInputs(void* inputs);
-		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
-		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
+		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
+		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
 		int   RiftIsPresent();
-		void  RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type);
+		void  RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type);
 		int   MeltingIsPresent();
-		void  MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type);
+		void  MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type);
 		DataSet* Copy(void);
-		void  Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type);
-		void  Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
-		void  Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type);
+		void  Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
+		void  Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
+		void  Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
 		void  VelocityExtrude(Vec ug,double* ug_serial);
 		void  SlopeExtrude(Vec sg,double* sg_serial);
Index: issm/trunk/src/c/Dux/Dux.cpp
===================================================================
--- issm/trunk/src/c/Dux/Dux.cpp	(revision 464)
+++ issm/trunk/src/c/Dux/Dux.cpp	(revision 465)
@@ -14,5 +14,5 @@
 
 void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, 
-		double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type){
+		double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -33,5 +33,5 @@
 
 	/*Compute velocity difference: */
-	elements->Du(du_g, u_g,u_g_obs,inputs,analysis_type);
+	elements->Du(du_g, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
 
 	/*Assemble vector: */
Index: issm/trunk/src/c/Dux/Dux.h
===================================================================
--- issm/trunk/src/c/Dux/Dux.h	(revision 464)
+++ issm/trunk/src/c/Dux/Dux.h	(revision 465)
@@ -10,5 +10,5 @@
 /* local prototypes: */
 void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, 
-		double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type);
+		double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
 
 #endif  /* _DUX_H */
Index: issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp
===================================================================
--- issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp	(revision 464)
+++ issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp	(revision 465)
@@ -15,47 +15,50 @@
 	}
 
-	if (strcmp(analysis_type,"control")==0){
+	if (strcmp(analysis_type,"diagnostic")==0){
+		return DiagnosticAnalysisEnum();
+	}
+	else if (strcmp(analysis_type,"control")==0){
 		return ControlAnalysisEnum();
 	}
-	else if (strcmp(analysis_type,"diagnostic_horiz")==0){
-		return DiagnosticHorizAnalysisEnum();
-	}
-	else if (strcmp(analysis_type,"diagnostic_vert")==0){
-		return DiagnosticVertAnalysisEnum();
+	else if (strcmp(analysis_type,"thermal")==0){
+		return ThermalAnalysisEnum();
 	}
 	else if (strcmp(analysis_type,"prognostic")==0){
 		return PrognosticAnalysisEnum();
 	}
-	else if (strcmp(analysis_type,"thermalsteady")==0){
-		return ThermalSteadyAnalysisEnum();
-	}
-	else if (strcmp(analysis_type,"thermaltransient")==0){
-		return ThermalTransientAnalysisEnum();
-	}
 	else if (strcmp(analysis_type,"melting")==0){
 		return MeltingAnalysisEnum();
-	}
-	else if (strcmp(analysis_type,"diagnostic_stokes")==0){
-		return DiagnosticStokesAnalysisEnum();
-	}
-	else if (strcmp(analysis_type,"diagnostic_hutter")==0){
-		return DiagnosticHutterAnalysisEnum();
 	}
 	else if (strcmp(analysis_type,"slope_compute")==0){
 		return SlopeComputeAnalysisEnum();
 	}
-	else if (strcmp(analysis_type,"surface_slope_compute_x")==0){
-		return SurfaceSlopeComputeXAnalysisEnum();
+	else if (strcmp(analysis_type,"stokes")==0){
+		return StokesAnalysisEnum();
 	}
-	else if (strcmp(analysis_type,"surface_slope_compute_y")==0){
-		return SurfaceSlopeComputeYAnalysisEnum();
+	else if (strcmp(analysis_type,"hutter")==0){
+		return HutterAnalysisEnum();
 	}
-	else if (strcmp(analysis_type,"bed_slope_compute_x")==0){
-		return BedSlopeComputeXAnalysisEnum();
+	else if (strcmp(analysis_type,"surfacex")==0){
+		return SurfaceXAnalysisEnum();
 	}
-	else if (strcmp(analysis_type,"bed_slope_compute_y")==0){
-		return BedSlopeComputeYAnalysisEnum();
+	else if (strcmp(analysis_type,"surfacey")==0){
+		return SurfaceYAnalysisEnum();
 	}
-	else throw ErrorException(__FUNCT__," could not recognized analysis type!");
+	else if (strcmp(analysis_type,"bedx")==0){
+		return BedXAnalysisEnum();
+	}
+	else if (strcmp(analysis_type,"bedy")==0){
+		return BedYAnalysisEnum();
+	}
+	else if (strcmp(analysis_type,"horiz")==0){
+		return HorizAnalysisEnum();
+	}
+	else if (strcmp(analysis_type,"vert")==0){
+		return VertAnalysisEnum();
+	}
+	else if (strcmp(analysis_type,"")==0){
+		return NoneAnalysisEnum();
+	}
+	else throw ErrorException(__FUNCT__,exprintf("%s%s"," could not recognize analysis type: ",analysis_type));
 
 }
Index: issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp
===================================================================
--- issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 464)
+++ issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 465)
@@ -15,20 +15,19 @@
 
 /*analysis types: */
-int StructuralAnalysisEnum(void){       return          20; }
-int ThermalAnalysisEnum(void){          return          21; }
-int ControlAnalysisEnum(void){          return          22; }
-int DiagnosticHorizAnalysisEnum(void){  return          23; }
-int DiagnosticVertAnalysisEnum(void){   return          24; }
-int PrognosticAnalysisEnum(void){       return          26; }
-int ThermalSteadyAnalysisEnum(void)	{   return          27; }
-int ThermalTransientAnalysisEnum(void){ return          28; }
-int MeltingAnalysisEnum(void){          return          29; }
-int DiagnosticStokesAnalysisEnum(void){ return          30; }
-int DiagnosticHutterAnalysisEnum(void){ return          31; }
-int SlopeComputeAnalysisEnum(void){ return              32; }
-int SurfaceSlopeComputeXAnalysisEnum(void){ return      33; }
-int SurfaceSlopeComputeYAnalysisEnum(void){ return      34; }
-int BedSlopeComputeXAnalysisEnum(void){ return          35; }
-int BedSlopeComputeYAnalysisEnum(void){ return          36; }
+int DiagnosticAnalysisEnum(void){  return               20; }
+int ControlAnalysisEnum(void){          return          21; }
+int ThermalAnalysisEnum(void){          return          22; }
+int PrognosticAnalysisEnum(void){          return       23; }
+int MeltingAnalysisEnum(void){          return          24; }
+int SlopeComputeAnalysisEnum(void){ return              25; }
+int StokesAnalysisEnum(void){ return                    26; }
+int HutterAnalysisEnum(void){ return                    27; }
+int SurfaceXAnalysisEnum(void){ return                  28; }
+int SurfaceYAnalysisEnum(void){ return                  29; }
+int BedXAnalysisEnum(void){ return                      30; }
+int BedYAnalysisEnum(void){ return                      31; }
+int HorizAnalysisEnum(void){ return                     32; }
+int VertAnalysisEnum(void){ return                      33; }
+int NoneAnalysisEnum(void){ return                      34; }
 
 /*datasets: */
Index: issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 464)
+++ issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 465)
@@ -34,18 +34,19 @@
 
 /* analysis enums: */
+int DiagnosticAnalysisEnum(void);
 int ControlAnalysisEnum(void);
-int DiagnosticHorizAnalysisEnum(void);
-int DiagnosticVertAnalysisEnum(void);
+int ThermalAnalysisEnum(void);
 int PrognosticAnalysisEnum(void);
-int ThermalSteadyAnalysisEnum(void);
-int ThermalTransientAnalysisEnum(void);
 int MeltingAnalysisEnum(void);
-int DiagnosticStokesAnalysisEnum(void);
-int DiagnosticHutterAnalysisEnum(void);
 int SlopeComputeAnalysisEnum(void);
-int SurfaceSlopeComputeXAnalysisEnum(void);
-int SurfaceSlopeComputeYAnalysisEnum(void);
-int BedSlopeComputeXAnalysisEnum(void);
-int BedSlopeComputeYAnalysisEnum(void);
+int StokesAnalysisEnum(void);
+int HutterAnalysisEnum(void);
+int SurfaceXAnalysisEnum(void);
+int SurfaceYAnalysisEnum(void);
+int BedXAnalysisEnum(void);
+int BedYAnalysisEnum(void);
+int HorizAnalysisEnum(void);
+int VertAnalysisEnum(void);
+int NoneAnalysisEnum(void);
 
 
Index: issm/trunk/src/c/Gradjx/Gradjx.cpp
===================================================================
--- issm/trunk/src/c/Gradjx/Gradjx.cpp	(revision 464)
+++ issm/trunk/src/c/Gradjx/Gradjx.cpp	(revision 465)
@@ -14,5 +14,5 @@
 
 void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, 
-		double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,char* control_type){
+		double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){
 
 	int i;
@@ -33,5 +33,5 @@
 
 	/*Compute gradients: */
-	elements->Gradj(grad_g, u_g,lambda_g,inputs,analysis_type,control_type);
+	elements->Gradj(grad_g, u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type);
 
 	/*Assemble vector: */
Index: issm/trunk/src/c/Gradjx/Gradjx.h
===================================================================
--- issm/trunk/src/c/Gradjx/Gradjx.h	(revision 464)
+++ issm/trunk/src/c/Gradjx/Gradjx.h	(revision 465)
@@ -10,5 +10,5 @@
 /* local prototypes: */
 void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, 
-		double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,char* control_type);
+		double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type);
 
 #endif  /* _GRADJX_H */
Index: issm/trunk/src/c/Makefile.am
===================================================================
--- issm/trunk/src/c/Makefile.am	(revision 464)
+++ issm/trunk/src/c/Makefile.am	(revision 465)
@@ -79,4 +79,5 @@
 					./shared/Dofs/dofs.h\
 					./shared/Dofs/dofsetgen.cpp\
+					./shared/Dofs/DistributeNumDofs.cpp\
 					./shared/Numerics/numerics.h\
 					./shared/Numerics/GaussPoints.h\
@@ -86,251 +87,4 @@
 					./shared/Numerics/BrentSearch.cpp\
 					./shared/Numerics/OptFunc.cpp\
-					./shared/Numerics/extrema.cpp\
-					./shared/Numerics/DistributeNumDofs.cpp\
-					./shared/Exceptions/exceptions.h\
-					./shared/Exceptions/Exceptions.cpp\
-					./shared/Exceptions/exprintf.cpp\
-					./shared/Exp/exp.h\
-					./shared/Exp/IsInPoly.cpp\
-					./shared/Exp/IsInPolySerial.cpp\
-					./shared/Exp/DomainOutlineRead.cpp\
-					./shared/TriMesh/trimesh.h\
-					./shared/TriMesh/AssociateSegmentToElement.cpp\
-					./shared/TriMesh/GridInsideHole.cpp\
-					./shared/TriMesh/OrderSegments.cpp\
-					./shared/TriMesh/SplitMeshForRifts.cpp\
-					./shared/TriMesh/TriMeshUtils.cpp\
-					./shared/Sorting/binary_search.cpp\
-					./shared/Sorting/sorting.h\
-					./shared/Elements/elements.h\
-					./shared/Elements/ResolvePointers.cpp\
-					./shared/Elements/Paterson.cpp\
-					./shared/Elements/GetElementNodeData.cpp\
-					./toolkits/petsc\
-					./toolkits/petsc/patches\
-					./toolkits/petsc/patches/petscpatches.h\
-					./toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp\
-					./toolkits/petsc/patches/MatlabVectorToPetscVector.cpp\
-					./toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp\
-					./toolkits/petsc/patches/PetscVectorToMatlabVector.cpp\
-					./toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp\
-					./toolkits/petsc/patches/MatlabVectorToDoubleVector.cpp\
-					./toolkits/petsc/patches/PetscDetermineLocalSize.cpp\
-					./toolkits/petsc/patches/VecTranspose.cpp\
-					./toolkits/petsc/patches/VecToMPISerial.cpp\
-					./toolkits/petsc/patches/MatToSerial.cpp\
-					./toolkits/petsc/patches/VecMerge.cpp\
-					./toolkits/petsc/patches/NewVec.cpp\
-					./toolkits/petsc/patches/NewVecFromLocalSize.cpp\
-					./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp\
-					./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\
-					./toolkits/petsc/patches/NewMat.cpp\
-					./toolkits/petsc/patches/VecFree.cpp\
-					./toolkits/petsc/patches/VecDuplicatePatch.cpp\
-					./toolkits/petsc/patches/KSPFree.cpp\
-					./toolkits/petsc/patches/ISFree.cpp\
-					./toolkits/petsc/patches/MatFree.cpp\
-					./toolkits/petsc/patches/GetOwnershipBoundariesFromRange.cpp\
-					./toolkits/petsc/patches/VecPartition.cpp\
-					./toolkits/petsc/patches/MatPartition.cpp\
-					./toolkits/petsc/patches/MatInvert.cpp\
-					./toolkits/petsc/patches/MatMultPatch.cpp\
-					./toolkits/petsc/petscincludes.h\
-					./toolkits/mpi/mpiincludes.h\
-					./toolkits/mpi/patches/mpipatches.h\
-					./toolkits/mpi/patches/MPI_Upperrow.cpp\
-					./toolkits/mpi/patches/MPI_Lowerrow.cpp\
-					./toolkits/mpi/patches/MPI_Boundariesfromrange.cpp\
-					./toolkits/metis/metisincludes.h\
-					./toolkits/triangle/triangleincludes.h\
-					./toolkits.h\
-					./io/io.h\
-					./io/FetchData.cpp\
-					./io/WriteData.cpp\
-					./io/WriteDataToDisk.cpp\
-					./io/SerialFetchData.cpp\
-					./io/SerialWriteData.cpp\
-					./io/ParallelFetchData.cpp\
-					./io/ParallelFetchInteger.cpp\
-					./io/ParallelFetchMat.cpp\
-					./io/ParallelFetchScalar.cpp\
-					./io/ParallelFetchString.cpp\
-					./io/ModelFetchData.cpp\
-					./io/WriteNodeSets.cpp\
-					./io/WriteParams.cpp\
-					./io/FetchNodeSets.cpp\
-					./io/FetchRifts.cpp\
-					./io/ParameterInputsInit.cpp\
-					./EnumDefinitions/EnumDefinitions.h\
-					./EnumDefinitions/EnumDefinitions.cpp\
-					./EnumDefinitions/AnalysisTypeAsEnum.cpp\
-					./ModelProcessorx/Model.h\
-					./ModelProcessorx/Model.cpp\
-					./ModelProcessorx/CreateElementsNodesAndMaterials.cpp\
-					./ModelProcessorx/CreateConstraints.cpp \
-					./ModelProcessorx/CreateLoads.cpp\
-					./ModelProcessorx/CreateParameters.cpp\
-					./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\
-					./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
-					./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
-					./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\
-					./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
-					./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
-					./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
-					./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\
-					./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
-					./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\
-					./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\
-					./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
-					./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
-					./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\
-					./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \
-					./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\
-					./ModelProcessorx/Control/CreateParametersControl.cpp\
-					./Dofx/Dofx.h\
-					./Dofx/Dofx.cpp\
-					./Dux/Dux.h\
-					./Dux/Dux.cpp\
-					./GriddataMeshToGridx/GriddataMeshToGridx.h\
-					./GriddataMeshToGridx/GriddataMeshToGridx.cpp\
-					./ControlConstrainx/ControlConstrainx.h\
-					./ControlConstrainx/ControlConstrainx.cpp\
-					./Misfitx/Misfitx.h\
-					./Misfitx/Misfitx.cpp\
-					./Orthx/Orthx.h\
-					./Orthx/Orthx.cpp\
-					./Gradjx/Gradjx.h\
-					./Gradjx/Gradjx.cpp\
-					./UpdateFromInputsx/UpdateFromInputsx.h\
-					./UpdateFromInputsx/UpdateFromInputsx.cpp\
-					./ConfigureObjectsx/ConfigureObjectsx.h\
-					./ConfigureObjectsx/ConfigureObjectsx.cpp\
-					./ComputePressurex/ComputePressurex.h\
-					./ComputePressurex/ComputePressurex.cpp\
-					./BuildNodeSetsx/BuildNodeSetsx.h\
-					./BuildNodeSetsx/BuildNodeSetsx.cpp\
-					./BuildNodeSetsx/PartitionSets.cpp\
-					./SpcNodesx/SpcNodesx.h\
-					./SpcNodesx/SpcNodesx.cpp\
-					./MpcNodesx/MpcNodesx.h\
-					./MpcNodesx/MpcNodesx.cpp\
-					./DataInterpx/DataInterpx.cpp\
-					./DataInterpx/DataInterpx.h\
-					./HoleFillerx/HoleFillerx.cpp\
-					./HoleFillerx/HoleFillerx.h\
-					./MeshPartitionx/MeshPartitionx.cpp\
-					./MeshPartitionx/MeshPartitionx.h\
-					./ContourToMeshx/ContourToMeshx.cpp\
-					./ContourToMeshx/ContourToMeshx.h\
-					./Reducevectorgtosx/Reducevectorgtosx.cpp\
-					./Reducevectorgtosx/Reducevectorgtosx.h\
-					./Reducematrixfromgtofx/Reducematrixfromgtofx.cpp\
-					./Reducematrixfromgtofx/Reducematrixfromgton.cpp\
-					./Reducematrixfromgtofx/Reducematrixfromgtofx.h\
-					./Reduceloadfromgtofx/Reduceloadfromgtofx.h\
-					./Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\
-					./NormalizeConstraintsx/NormalizeConstraintsx.cpp\
-					./NormalizeConstraintsx/NormalizeConstraintsx.h\
-					./SystemMatricesx/SystemMatricesx.cpp\
-					./SystemMatricesx/SystemMatricesx.h\
-					./PenaltyConstraintsx/PenaltyConstraintsx.cpp\
-					./PenaltyConstraintsx/PenaltyConstraintsx.h\
-					./PenaltySystemMatricesx/PenaltySystemMatricesx.cpp\
-					./PenaltySystemMatricesx/PenaltySystemMatricesx.h\
-					./Solverx/Solverx.cpp\
-					./Solverx/Solverx.h\
-					./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
-					./Mergesolutionfromftogx/Mergesolutionfromftogx.h\
-					./ProcessParamsx/ProcessParamsx.cpp\
-					./ProcessParamsx/ProcessParamsx.h\
-					./VelocityExtrudex/VelocityExtrudex.cpp\
-					./VelocityExtrudex/VelocityExtrudex.h\
-					./SlopeExtrudex/SlopeExtrudex.cpp\
-					./SlopeExtrudex/SlopeExtrudex.h
-
-
-
-
-libISSM_a_CXXFLAGS = -fPIC -DMATLAB -D_SERIAL_ -ansi -D_GNU_SOURCE -fno-omit-frame-pointer -pthread -D_CPP_
-if LARGEARRAYS
-libISSM_a_CXXFLAGS += -D__GCC4BUILD__  
-else
-libISSM_a_CXXFLAGS += -DMX_COMPAT_32
-endif
-
-
-
-
-#Parallel compilation
-libpISSM_a_SOURCES = ./objects/objects.h\
-					./objects/Object.h\
-					./objects/Element.h\
-					./objects/Element.cpp\
-					./objects/Material.h\
-					./objects/Material.cpp\
-					./objects/Load.h\
-					./objects/Load.cpp\
-					./objects/SolverEnum.h\
-					./objects/Contour.h\
-					./objects/Contour.cpp\
-					./objects/OptArgs.h\
-					./objects/OptPars.h\
-					./objects/Friction.h\
-					./objects/Friction.cpp\
-					./objects/Node.h\
-					./objects/Node.cpp\
-					./objects/Tria.h\
-					./objects/Tria.cpp\
-					./objects/Sing.h\
-					./objects/Sing.cpp\
-					./objects/Beam.h\
-					./objects/Beam.cpp\
-					./objects/Penta.h\
-					./objects/Penta.cpp\
-					./objects/Matice.h\
-					./objects/Matice.cpp\
-					./objects/Matpar.h\
-					./objects/Matpar.cpp\
-					./objects/Input.h\
-					./objects/Input.cpp\
-					./objects/ParameterInputs.h\
-					./objects/ParameterInputs.cpp\
-					./objects/Spc.cpp\
-					./objects/Spc.h\
-					./objects/Rgb.cpp\
-					./objects/Rgb.h\
-					./objects/Penpair.cpp\
-					./objects/Penpair.h\
-					./objects/Pengrid.cpp\
-					./objects/Pengrid.h\
-					./objects/Icefront.cpp\
-					./objects/Icefront.h\
-					./objects/Param.cpp\
-					./objects/Param.h\
-					./objects/NodeSets.cpp\
-					./objects/NodeSets.h\
-					./DataSet/DataSet.cpp\
-					./DataSet/DataSet.h\
-					./shared/shared.h\
-					./shared/Alloc/alloc.h\
-					./shared/Alloc/alloc.cpp\
-					./shared/Matlab/matlabshared.h\
-					./shared/Matlab/PrintfFunction.cpp\
-					./shared/Matlab/ModuleBoot.cpp\
-					./shared/Matlab/mxGetAssignedField.cpp\
-					./shared/Matlab/mxGetField.cpp\
-					./shared/Matlab/CheckNumMatlabArguments.cpp\
-					./shared/Matrix/matrix.h\
-					./shared/Matrix/MatrixUtils.cpp\
-					./shared/Dofs/dofs.h\
-					./shared/Dofs/dofsetgen.cpp\
-					./shared/Numerics/numerics.h\
-					./shared/Numerics/GaussPoints.h\
-					./shared/Numerics/GaussPoints.cpp\
-					./shared/Numerics/cross.cpp\
-					./shared/Numerics/norm.cpp\
-					./shared/Numerics/BrentSearch.cpp\
-					./shared/Numerics/OptFunc.cpp\
-					./shared/Numerics/DistributeNumDofs.cpp\
 					./shared/Numerics/extrema.cpp\
 					./shared/Exceptions/exceptions.h\
@@ -413,7 +167,249 @@
 					./ModelProcessorx/Model.h\
 					./ModelProcessorx/Model.cpp\
-					./ModelProcessorx/CreateElementsNodesAndMaterials.cpp\
-					./ModelProcessorx/CreateConstraints.cpp \
-					./ModelProcessorx/CreateLoads.cpp\
+					./ModelProcessorx/CreateDataSets.cpp\
+					./ModelProcessorx/CreateParameters.cpp\
+					./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\
+					./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
+					./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
+					./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\
+					./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
+					./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
+					./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
+					./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\
+					./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
+					./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\
+					./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\
+					./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
+					./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
+					./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\
+					./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \
+					./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\
+					./ModelProcessorx/Control/CreateParametersControl.cpp\
+					./Dofx/Dofx.h\
+					./Dofx/Dofx.cpp\
+					./Dux/Dux.h\
+					./Dux/Dux.cpp\
+					./GriddataMeshToGridx/GriddataMeshToGridx.h\
+					./GriddataMeshToGridx/GriddataMeshToGridx.cpp\
+					./ControlConstrainx/ControlConstrainx.h\
+					./ControlConstrainx/ControlConstrainx.cpp\
+					./Misfitx/Misfitx.h\
+					./Misfitx/Misfitx.cpp\
+					./Orthx/Orthx.h\
+					./Orthx/Orthx.cpp\
+					./Gradjx/Gradjx.h\
+					./Gradjx/Gradjx.cpp\
+					./UpdateFromInputsx/UpdateFromInputsx.h\
+					./UpdateFromInputsx/UpdateFromInputsx.cpp\
+					./ConfigureObjectsx/ConfigureObjectsx.h\
+					./ConfigureObjectsx/ConfigureObjectsx.cpp\
+					./ComputePressurex/ComputePressurex.h\
+					./ComputePressurex/ComputePressurex.cpp\
+					./BuildNodeSetsx/BuildNodeSetsx.h\
+					./BuildNodeSetsx/BuildNodeSetsx.cpp\
+					./BuildNodeSetsx/PartitionSets.cpp\
+					./SpcNodesx/SpcNodesx.h\
+					./SpcNodesx/SpcNodesx.cpp\
+					./MpcNodesx/MpcNodesx.h\
+					./MpcNodesx/MpcNodesx.cpp\
+					./DataInterpx/DataInterpx.cpp\
+					./DataInterpx/DataInterpx.h\
+					./HoleFillerx/HoleFillerx.cpp\
+					./HoleFillerx/HoleFillerx.h\
+					./MeshPartitionx/MeshPartitionx.cpp\
+					./MeshPartitionx/MeshPartitionx.h\
+					./ContourToMeshx/ContourToMeshx.cpp\
+					./ContourToMeshx/ContourToMeshx.h\
+					./Reducevectorgtosx/Reducevectorgtosx.cpp\
+					./Reducevectorgtosx/Reducevectorgtosx.h\
+					./Reducematrixfromgtofx/Reducematrixfromgtofx.cpp\
+					./Reducematrixfromgtofx/Reducematrixfromgton.cpp\
+					./Reducematrixfromgtofx/Reducematrixfromgtofx.h\
+					./Reduceloadfromgtofx/Reduceloadfromgtofx.h\
+					./Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\
+					./NormalizeConstraintsx/NormalizeConstraintsx.cpp\
+					./NormalizeConstraintsx/NormalizeConstraintsx.h\
+					./SystemMatricesx/SystemMatricesx.cpp\
+					./SystemMatricesx/SystemMatricesx.h\
+					./PenaltyConstraintsx/PenaltyConstraintsx.cpp\
+					./PenaltyConstraintsx/PenaltyConstraintsx.h\
+					./PenaltySystemMatricesx/PenaltySystemMatricesx.cpp\
+					./PenaltySystemMatricesx/PenaltySystemMatricesx.h\
+					./Solverx/Solverx.cpp\
+					./Solverx/Solverx.h\
+					./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
+					./Mergesolutionfromftogx/Mergesolutionfromftogx.h\
+					./ProcessParamsx/ProcessParamsx.cpp\
+					./ProcessParamsx/ProcessParamsx.h\
+					./VelocityExtrudex/VelocityExtrudex.cpp\
+					./VelocityExtrudex/VelocityExtrudex.h\
+					./SlopeExtrudex/SlopeExtrudex.cpp\
+					./SlopeExtrudex/SlopeExtrudex.h
+
+
+
+
+libISSM_a_CXXFLAGS = -fPIC -DMATLAB -D_SERIAL_ -ansi -D_GNU_SOURCE -fno-omit-frame-pointer -pthread -D_CPP_
+if LARGEARRAYS
+libISSM_a_CXXFLAGS += -D__GCC4BUILD__  
+else
+libISSM_a_CXXFLAGS += -DMX_COMPAT_32
+endif
+
+
+
+
+#Parallel compilation
+libpISSM_a_SOURCES = ./objects/objects.h\
+					./objects/Object.h\
+					./objects/Element.h\
+					./objects/Element.cpp\
+					./objects/Material.h\
+					./objects/Material.cpp\
+					./objects/Load.h\
+					./objects/Load.cpp\
+					./objects/SolverEnum.h\
+					./objects/Contour.h\
+					./objects/Contour.cpp\
+					./objects/OptArgs.h\
+					./objects/OptPars.h\
+					./objects/Friction.h\
+					./objects/Friction.cpp\
+					./objects/Node.h\
+					./objects/Node.cpp\
+					./objects/Tria.h\
+					./objects/Tria.cpp\
+					./objects/Sing.h\
+					./objects/Sing.cpp\
+					./objects/Beam.h\
+					./objects/Beam.cpp\
+					./objects/Penta.h\
+					./objects/Penta.cpp\
+					./objects/Matice.h\
+					./objects/Matice.cpp\
+					./objects/Matpar.h\
+					./objects/Matpar.cpp\
+					./objects/Input.h\
+					./objects/Input.cpp\
+					./objects/ParameterInputs.h\
+					./objects/ParameterInputs.cpp\
+					./objects/Spc.cpp\
+					./objects/Spc.h\
+					./objects/Rgb.cpp\
+					./objects/Rgb.h\
+					./objects/Penpair.cpp\
+					./objects/Penpair.h\
+					./objects/Pengrid.cpp\
+					./objects/Pengrid.h\
+					./objects/Icefront.cpp\
+					./objects/Icefront.h\
+					./objects/Param.cpp\
+					./objects/Param.h\
+					./objects/NodeSets.cpp\
+					./objects/NodeSets.h\
+					./DataSet/DataSet.cpp\
+					./DataSet/DataSet.h\
+					./shared/shared.h\
+					./shared/Alloc/alloc.h\
+					./shared/Alloc/alloc.cpp\
+					./shared/Matlab/matlabshared.h\
+					./shared/Matlab/PrintfFunction.cpp\
+					./shared/Matlab/ModuleBoot.cpp\
+					./shared/Matlab/mxGetAssignedField.cpp\
+					./shared/Matlab/mxGetField.cpp\
+					./shared/Matlab/CheckNumMatlabArguments.cpp\
+					./shared/Matrix/matrix.h\
+					./shared/Matrix/MatrixUtils.cpp\
+					./shared/Dofs/dofs.h\
+					./shared/Dofs/dofsetgen.cpp\
+					./shared/Dofs/DistributeNumDofs.cpp\
+					./shared/Numerics/numerics.h\
+					./shared/Numerics/GaussPoints.h\
+					./shared/Numerics/GaussPoints.cpp\
+					./shared/Numerics/cross.cpp\
+					./shared/Numerics/norm.cpp\
+					./shared/Numerics/BrentSearch.cpp\
+					./shared/Numerics/OptFunc.cpp\
+					./shared/Numerics/extrema.cpp\
+					./shared/Exceptions/exceptions.h\
+					./shared/Exceptions/Exceptions.cpp\
+					./shared/Exceptions/exprintf.cpp\
+					./shared/Exp/exp.h\
+					./shared/Exp/IsInPoly.cpp\
+					./shared/Exp/IsInPolySerial.cpp\
+					./shared/Exp/DomainOutlineRead.cpp\
+					./shared/TriMesh/trimesh.h\
+					./shared/TriMesh/AssociateSegmentToElement.cpp\
+					./shared/TriMesh/GridInsideHole.cpp\
+					./shared/TriMesh/OrderSegments.cpp\
+					./shared/TriMesh/SplitMeshForRifts.cpp\
+					./shared/TriMesh/TriMeshUtils.cpp\
+					./shared/Sorting/binary_search.cpp\
+					./shared/Sorting/sorting.h\
+					./shared/Elements/elements.h\
+					./shared/Elements/ResolvePointers.cpp\
+					./shared/Elements/Paterson.cpp\
+					./shared/Elements/GetElementNodeData.cpp\
+					./toolkits/petsc\
+					./toolkits/petsc/patches\
+					./toolkits/petsc/patches/petscpatches.h\
+					./toolkits/petsc/patches/MatlabMatrixToPetscMatrix.cpp\
+					./toolkits/petsc/patches/MatlabVectorToPetscVector.cpp\
+					./toolkits/petsc/patches/PetscMatrixToMatlabMatrix.cpp\
+					./toolkits/petsc/patches/PetscVectorToMatlabVector.cpp\
+					./toolkits/petsc/patches/MatlabMatrixToDoubleMatrix.cpp\
+					./toolkits/petsc/patches/MatlabVectorToDoubleVector.cpp\
+					./toolkits/petsc/patches/PetscDetermineLocalSize.cpp\
+					./toolkits/petsc/patches/VecTranspose.cpp\
+					./toolkits/petsc/patches/VecToMPISerial.cpp\
+					./toolkits/petsc/patches/MatToSerial.cpp\
+					./toolkits/petsc/patches/VecMerge.cpp\
+					./toolkits/petsc/patches/NewVec.cpp\
+					./toolkits/petsc/patches/NewVecFromLocalSize.cpp\
+					./toolkits/petsc/patches/PetscOptionsDetermineSolverType.cpp\
+					./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\
+					./toolkits/petsc/patches/NewMat.cpp\
+					./toolkits/petsc/patches/VecFree.cpp\
+					./toolkits/petsc/patches/VecDuplicatePatch.cpp\
+					./toolkits/petsc/patches/KSPFree.cpp\
+					./toolkits/petsc/patches/ISFree.cpp\
+					./toolkits/petsc/patches/MatFree.cpp\
+					./toolkits/petsc/patches/GetOwnershipBoundariesFromRange.cpp\
+					./toolkits/petsc/patches/VecPartition.cpp\
+					./toolkits/petsc/patches/MatPartition.cpp\
+					./toolkits/petsc/patches/MatInvert.cpp\
+					./toolkits/petsc/patches/MatMultPatch.cpp\
+					./toolkits/petsc/petscincludes.h\
+					./toolkits/mpi/mpiincludes.h\
+					./toolkits/mpi/patches/mpipatches.h\
+					./toolkits/mpi/patches/MPI_Upperrow.cpp\
+					./toolkits/mpi/patches/MPI_Lowerrow.cpp\
+					./toolkits/mpi/patches/MPI_Boundariesfromrange.cpp\
+					./toolkits/metis/metisincludes.h\
+					./toolkits/triangle/triangleincludes.h\
+					./toolkits.h\
+					./io/io.h\
+					./io/FetchData.cpp\
+					./io/WriteData.cpp\
+					./io/WriteDataToDisk.cpp\
+					./io/SerialFetchData.cpp\
+					./io/SerialWriteData.cpp\
+					./io/ParallelFetchData.cpp\
+					./io/ParallelFetchInteger.cpp\
+					./io/ParallelFetchMat.cpp\
+					./io/ParallelFetchScalar.cpp\
+					./io/ParallelFetchString.cpp\
+					./io/ModelFetchData.cpp\
+					./io/WriteNodeSets.cpp\
+					./io/WriteParams.cpp\
+					./io/FetchNodeSets.cpp\
+					./io/FetchRifts.cpp\
+					./io/ParameterInputsInit.cpp\
+					./EnumDefinitions/EnumDefinitions.h\
+					./EnumDefinitions/EnumDefinitions.cpp\
+					./EnumDefinitions/AnalysisTypeAsEnum.cpp\
+					./ModelProcessorx/Model.h\
+					./ModelProcessorx/Model.cpp\
+					./ModelProcessorx/CreateDataSets.cpp\
 					./ModelProcessorx/CreateParameters.cpp\
 					./ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp\
@@ -510,6 +506,5 @@
 bin_PROGRAMS = 
 else 
-bin_PROGRAMS = diagnostic.exe control.exe 
-#thermalsteady.exe
+bin_PROGRAMS = diagnostic.exe control.exe thermal.exe
 endif
 
@@ -522,4 +517,4 @@
 control_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
 
-thermalsteady_exe_SOURCES = parallel/thermalsteady.cpp
-thermalsteady_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
+thermal_exe_SOURCES = parallel/thermal.cpp
+thermal_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
Index: issm/trunk/src/c/Misfitx/Misfitx.cpp
===================================================================
--- issm/trunk/src/c/Misfitx/Misfitx.cpp	(revision 464)
+++ issm/trunk/src/c/Misfitx/Misfitx.cpp	(revision 465)
@@ -14,5 +14,5 @@
 
 void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, 
-		double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type){
+		double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 	
 	/*output: */
@@ -24,5 +24,5 @@
 
 	/*Compute gradients: */
-	elements->Misfit(&J, u_g,u_g_obs,inputs,analysis_type);
+	elements->Misfit(&J, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
 
 	/*Sum all J from all cpus of the cluster:*/
Index: issm/trunk/src/c/Misfitx/Misfitx.h
===================================================================
--- issm/trunk/src/c/Misfitx/Misfitx.h	(revision 464)
+++ issm/trunk/src/c/Misfitx/Misfitx.h	(revision 465)
@@ -10,5 +10,5 @@
 /* local prototypes: */
 void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, 
-		double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type);
+		double* u_g, double* u_g_obs, ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
 
 #endif  /* _MISFITX_H */
Index: issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 465)
@@ -13,8 +13,9 @@
 #include "../Model.h"
 
-void CreateParametersControl(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount){
+void CreateParametersControl(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
 	
 	int i;
 	
+	DataSet* parameters=NULL;
 	Param*   param = NULL;
 	int      count;
@@ -29,6 +30,7 @@
 	double* vy_obs=NULL;
 
-	/*Get counter : */
-	count=*pcount;
+	/*Get parameters: */
+	parameters=*pparameters;
+	count=parameters->Size();
 	
 	/*control_type: */
@@ -132,5 +134,5 @@
 
 	/*Assign output pointer: */
-	*pcount=count;
+	*pparameters=parameters;
 }
 
Index: issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp	(revision 464)
+++ 	(revision )
@@ -1,48 +1,0 @@
-/*!\file: CreateConstraints.cpp
- * \brief general driver for creating constraints
- */ 
-
-#ifdef HAVE_CONFIG_H
-	#include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#undef __FUNCT__ 
-#define __FUNCT__ "CreateConstraints"
-
-#include "./Model.h"
-#include "../shared/shared.h"
-
-void CreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
-	
-	/*This is just a high level driver: */
-	if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){
-		CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"diagnostic_vert")==0){
-		CreateConstraintsDiagnosticVert(pconstraints,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"diagnostic_stokes")==0){
-		CreateConstraintsDiagnosticStokes(pconstraints,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"diagnostic_hutter")==0){
-		CreateConstraintsDiagnosticHutter(pconstraints,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"slope_compute")==0){
-		CreateConstraintsSlopeCompute(pconstraints,model,model_handle);
-	}
-	/*
-	else if (strcmp(model->analysis_type,"melting")==0){
-		CreateConstraintsMelting(pconstraints,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"prognostic")==0){
-		CreateConstraintsPrognostic(pconstraints,model,model_handle);
-	}
-	else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){
-		CreateConstraintsThermal(pconstraints,model,model_handle);
-	}*/
-	else{
-		throw ErrorException(__FUNCT__," analysis_type not supported yet!");
-	}
-}
Index: issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 465)
+++ issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 465)
@@ -0,0 +1,104 @@
+/*!\file: CreateDataSets
+ * \brief general driver for creating all datasets that make a finite element model
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateElementsNodesAndMaterials"
+
+#include "./Model.h"
+#include "../shared/shared.h"
+
+
+void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,Model* model,ConstDataHandle model_handle){
+
+
+	/*create parameters common to all solutions: */
+	CreateParameters(pparameters,model,model_handle);
+
+	/*This is just a high level driver: */
+	if (strcmp(model->analysis_type,"control")==0){
+
+		CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, model,model_handle);
+		CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle);
+		CreateLoadsDiagnosticHoriz(ploads,model,model_handle);
+		CreateParametersControl(pparameters,model,model_handle);
+					
+	}
+	else if (strcmp(model->analysis_type,"diagnostic")==0){
+
+		if (strcmp(model->sub_analysis_type,"horiz")==0){
+			
+			CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, model,model_handle);
+			CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle);
+			CreateLoadsDiagnosticHoriz(ploads,model,model_handle);
+			CreateParametersDiagnosticHoriz(pparameters,model,model_handle);
+				
+		}
+		else if (strcmp(model->sub_analysis_type,"vert")==0){
+		
+			CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, model,model_handle);
+			CreateConstraintsDiagnosticVert(pconstraints,model,model_handle);
+			CreateLoadsDiagnosticVert(ploads,model,model_handle);
+			
+		}
+		else if (strcmp(model->sub_analysis_type,"stokes")==0){
+
+			CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes,pmaterials, model,model_handle);
+			CreateConstraintsDiagnosticStokes(pconstraints,model,model_handle);
+			CreateLoadsDiagnosticStokes(ploads,model,model_handle);
+			
+		}
+		else if (strcmp(model->sub_analysis_type,"hutter")==0){
+
+			CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pmaterials, model,model_handle);
+			CreateConstraintsDiagnosticHutter(pconstraints,model,model_handle);
+			CreateLoadsDiagnosticHutter(ploads,model,model_handle);
+			
+		}
+	}
+	else if (strcmp(model->analysis_type,"slope_compute")==0){
+
+		CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes,pmaterials, model,model_handle);
+		CreateConstraintsSlopeCompute(pconstraints,model,model_handle);
+		CreateLoadsSlopeCompute(ploads,model,model_handle);
+			
+	
+	}
+	else if (strcmp(model->analysis_type,"thermal")==0){
+
+		if ((strcmp(model->sub_analysis_type,"steady")==0) || (strcmp(model->sub_analysis_type,"transient")==0)){
+		
+			//CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pmaterials, model,model_handle);
+			//CreateConstraintsThermal(pconstraints,model,model_handle);
+			//CreateLoadsThermal(ploads,model,model_handle);
+					
+		}
+		else if (strcmp(model->sub_analysis_type,"melting")==0){
+			
+			//CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pmaterials, model,model_handle);
+			//CreateConstraintsMelting(pconstraints,model,model_handle);
+			//CreateLoadsMelting(ploads,model,model_handle);
+
+		}
+	}
+	else if (strcmp(model->analysis_type,"prognostic")==0){
+			
+		//CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pmaterials, model,model_handle);
+		//CreateConstraintsPrognostic(pconstraints,model,model_handle);
+		//CreateLoadsPrognostic(ploads,model,model_handle);
+	}
+	else{
+		throw ErrorException(__FUNCT__,exprintf("%s%s%s%s"," analysis_type: ",model->analysis_type," sub_analysis_type: ",model->sub_analysis_type," not supported yet!"));
+	}
+
+
+
+
+
+}
Index: issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp	(revision 464)
+++ 	(revision )
@@ -1,63 +1,0 @@
-/*!\file: CreateElementsNodesAndMaterials.cpp
- * \brief general driver for creating elements, nodes and materials
- */ 
-
-#ifdef HAVE_CONFIG_H
-	#include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#undef __FUNCT__ 
-#define __FUNCT__ "CreateElementsNodesAndMaterials"
-
-#include "./Model.h"
-#include "../shared/shared.h"
-
-
-void CreateElementsNodesAndMaterials(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){
-
-	/*This is just a high level driver: */
-	if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){
-
-		CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, model,model_handle);
-	
-	}
-	else if ((strcmp(model->analysis_type,"diagnostic_vert")==0)){
-
-		CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, model,model_handle);
-	
-	}
-	else if ((strcmp(model->analysis_type,"diagnostic_stokes")==0)){
-
-		CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes,pmaterials, model,model_handle);
-	
-	}
-	else if ((strcmp(model->analysis_type,"diagnostic_hutter")==0)){
-
-		CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pmaterials, model,model_handle);
-	
-	}
-	else if ((strcmp(model->analysis_type,"slope_compute")==0)){
-
-		CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes,pmaterials, model,model_handle);
-	
-	}
-	/*
-	else if (strcmp(model->analysis_type,"melting")==0){
-
-		CreateElementsNodesAndMaterialsDataSetsMelting(pelements,pnodes,pmaterials, model,model->analysis_type);
-	}
-	else if (strcmp(model->analysis_type,"prognostic")==0){
-
-		CreateElementsNodesAndMaterialsDataSetsPrognostic(pelements,pnodes,pmaterials, model,model->analysis_type);
-	}
-	else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){
-
-		CreateElementsNodesAndMaterialsDataSetsThermal(pelements,pnodes,pmaterials, model,model->analysis_type);
-	}*/
-	else{
-		throw ErrorException(__FUNCT__,exprintf("%s%s%s"," analysis_type: ",model->analysis_type," not supported yet!"));
-	}
-
-}
Index: issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp	(revision 464)
+++ 	(revision )
@@ -1,57 +1,0 @@
-/*!\file:  CreateLoads.cpp
- * \brief create loads datasets. General driver.
- */ 
-
-#ifdef HAVE_CONFIG_H
-	#include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#undef __FUNCT__ 
-#define __FUNCT__ "CreateLoads"
-
-#include "./Model.h"
-#include "../shared/shared.h"
-
-void CreateLoads(DataSet** ploads, Model* model,ConstDataHandle model_handle){
-
-	/*This is just a high level driver: */
-	if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){
-
-		CreateLoadsDiagnosticHoriz(ploads,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"diagnostic_vert")==0){
-
-		CreateLoadsDiagnosticVert(ploads,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"diagnostic_stokes")==0){
-
-		CreateLoadsDiagnosticStokes(ploads,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"diagnostic_hutter")==0){
-
-		CreateLoadsDiagnosticHutter(ploads,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"slope_compute")==0){
-
-		CreateLoadsSlopeCompute(ploads,model,model_handle);
-	}
-	/*
-	else if (strcmp(model->analysis_type,"melting")==0){
-
-		CreateLoadsMelting(ploads,model,model_handle);
-	}
-	else if (strcmp(model->analysis_type,"prognostic")==0){
-
-		CreateLoadsPrognostic(ploads,model,model_handle);
-	}
-	else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){
-	
-		CreateLoadsThermal(ploads,model,model_handle);
-
-	}*/
-	else{
-		throw ErrorException(__FUNCT__," analysis_type not supported yet!");
-	}
-}
Index: issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 465)
@@ -21,4 +21,5 @@
 	int      count=1;
 	int      analysis_type;
+	int      sub_analysis_type;
 	int      numberofdofspernode;
 	int      dim;
@@ -29,8 +30,14 @@
 	//Get analysis_type:
 	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+	sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
 	
 	count++;
 	param= new Param(count,"analysis_type",INTEGER);
 	param->SetInteger(analysis_type);
+	parameters->AddObject(param);
+
+	count++;
+	param= new Param(count,"sub_analysis_type",INTEGER);
+	param->SetInteger(sub_analysis_type);
 	parameters->AddObject(param);
 
@@ -187,5 +194,5 @@
 
 	/*Deal with numberofdofspernode: */
-	DistributeNumDofs(&numberofdofspernode,analysis_type);
+	DistributeNumDofs(&numberofdofspernode,analysis_type,sub_analysis_type);
 
 	count++;
@@ -194,11 +201,4 @@
 	parameters->AddObject(param);
 
-	/*Deal with control parameters: */
-	if (analysis_type==ControlAnalysisEnum()){
-		CreateParametersControl(parameters,model,model_handle,&count);
-	}
-	if (analysis_type==DiagnosticHorizAnalysisEnum()){
-		CreateParametersDiagnosticHoriz(parameters,model,model_handle,&count);
-	}
 
 	
Index: issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 465)
@@ -38,4 +38,5 @@
 
 	int         analysis_type;
+	int         sub_analysis_type;
 	
 	/*output: */
@@ -154,4 +155,5 @@
 	/*Get analysis_type: */
 	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+	sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
 
 	
@@ -571,5 +573,5 @@
 
 	/*Get number of dofs per node: */
-	DistributeNumDofs(&node_numdofs,analysis_type);
+	DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
 
 	for (i=0;i<model->numberofnodes;i++){
Index: issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp	(revision 465)
@@ -13,7 +13,8 @@
 #include "../Model.h"
 
-void CreateParametersDiagnosticHoriz(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount){
+void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
 	
 	Param*   param = NULL;
+	DataSet* parameters=NULL;
 	int      count;
 	int i;
@@ -23,9 +24,11 @@
 	double* vz=NULL;
 
+	/*recover parameters : */
+	parameters=*pparameters;
+
+	count=parameters->Size();
+
 	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
 	if (!model->ismacayealpattyn)return;
-
-	/*Get counter : */
-	count=*pcount;
 
 	/*Get vx and vy: */
@@ -61,5 +64,5 @@
 	
 	/*Assign output pointer: */
-	*pcount=count;
+	*pparameters=parameters;
 }
 
Index: issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 465)
@@ -40,4 +40,5 @@
 
 	int         analysis_type;
+	int         sub_analysis_type;
 	
 	/*output: */
@@ -106,4 +107,5 @@
 	/*Get analysis_type: */
 	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+	sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
 
 	/*Width of elements: */
@@ -297,5 +299,5 @@
 	
 	/*Get number of dofs per node: */
-	DistributeNumDofs(&node_numdofs,analysis_type);
+	DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
 
 	for (i=0;i<model->numberofnodes;i++){
Index: issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 465)
@@ -38,4 +38,5 @@
 
 	int         analysis_type;
+	int         sub_analysis_type;
 	
 	/*output: */
@@ -134,4 +135,5 @@
 	/*Get analysis_type: */
 	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+	sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
 
 	/*Width of elements: */
@@ -428,5 +430,5 @@
 
 	/*Get number of dofs per node: */
-	DistributeNumDofs(&node_numdofs,analysis_type);
+	DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
 
 	for (i=0;i<model->numberofnodes;i++){
Index: issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 465)
@@ -37,4 +37,5 @@
 
 	int         analysis_type;
+	int         sub_analysis_type;
 	
 	/*output: */
@@ -143,4 +144,5 @@
 	/*Get analysis_type: */
 	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+	sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
 
 	/*Width of elements: */
@@ -354,5 +356,5 @@
 	
 	/*Get number of dofs per node: */
-	DistributeNumDofs(&node_numdofs,analysis_type);
+	DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
 
 	for (i=0;i<model->numberofnodes;i++){
Index: issm/trunk/src/c/ModelProcessorx/Model.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Model.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/Model.cpp	(revision 465)
@@ -37,4 +37,5 @@
 	model->meshtype=NULL;
 	model->analysis_type=NULL;
+	model->sub_analysis_type=NULL;
 	model->numberofelements=0;
 	model->numberofnodes=0;
@@ -243,4 +244,5 @@
 	xfree((void**)&model->meshtype);
 	xfree((void**)&model->analysis_type);
+	xfree((void**)&model->sub_analysis_type);
 	
 	if(model->numrifts){
@@ -289,4 +291,5 @@
 	/*In ModelInit, we get all the data that is not difficult to get, and that is small: */
 	ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL); 
+	ModelFetchData((void**)&model->sub_analysis_type,NULL,NULL,model_handle,"sub_analysis_type","String",NULL); 
 
 	ModelFetchData((void**)&model->meshtype,NULL,NULL,model_handle,"type","String",NULL);
Index: issm/trunk/src/c/ModelProcessorx/Model.h
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Model.h	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/Model.h	(revision 465)
@@ -16,4 +16,5 @@
 	char*   meshtype;
 	char*   analysis_type;
+	char*   sub_analysis_type;
 	char*   solverstring;
 
@@ -174,7 +175,5 @@
 
 	/*Creation of fem datasets: general drivers*/
-	void	CreateElementsNodesAndMaterials(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
-	void    CreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle);
-	void    CreateLoads(DataSet** ploads, Model* model, ConstDataHandle model_handle);
+	void    CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,Model* model,ConstDataHandle model_handle);
 	void    CreateParameters(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
 
@@ -186,5 +185,5 @@
 	void	CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
 	void    CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model, ConstDataHandle model_handle);
-	void    CreateParametersDiagnosticHoriz(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
+	void    CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
 
 	/*diagnostic vertical*/
@@ -192,5 +191,4 @@
 	void	CreateConstraintsDiagnosticVert(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
 	void    CreateLoadsDiagnosticVert(DataSet** ploads, Model* model, ConstDataHandle model_handle);
-	void    CreateParametersDiagnosticVert(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
 
 	/*diagnostic hutter*/
@@ -198,5 +196,4 @@
 	void	CreateConstraintsDiagnosticHutter(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
 	void    CreateLoadsDiagnosticHutter(DataSet** ploads, Model* model, ConstDataHandle model_handle);
-	void    CreateParametersDiagnosticHutter(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
 
 	/*diagnostic stokes*/
@@ -204,5 +201,4 @@
 	void	CreateConstraintsDiagnosticStokes(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
 	void    CreateLoadsDiagnosticStokes(DataSet** ploads, Model* model, ConstDataHandle model_handle);
-	void    CreateParametersDiagnosticStokes(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
 
 	/*slope compute*/
@@ -210,8 +206,7 @@
 	void	CreateConstraintsSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
 	void    CreateLoadsSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle);
-	void    CreateParametersSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
 
 	/*control:*/
-	void    CreateParametersControl(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount);
+	void    CreateParametersControl(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
 
 
Index: issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 464)
+++ issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 465)
@@ -36,4 +36,5 @@
 
 	int         analysis_type;
+	int         sub_analysis_type;
 	
 	/*output: */
@@ -131,4 +132,5 @@
 	/*Get analysis_type: */
 	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+	sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
 
 	/*Width of elements: */
@@ -397,5 +399,5 @@
 
 	/*Get number of dofs per node: */
-	DistributeNumDofs(&node_numdofs,analysis_type);
+	DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type);
 
 	for (i=0;i<model->numberofnodes;i++){
Index: issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp
===================================================================
--- issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp	(revision 464)
+++ issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp	(revision 465)
@@ -14,5 +14,5 @@
 
 void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes,
-		DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type){
+		DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -31,8 +31,8 @@
 	 * management routine, otherwise, skip : */
 	if (loads->RiftIsPresent()){
-		loads->RiftConstraints(&num_unstable_constraints,inputs,analysis_type);
+		loads->RiftConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type);
 	}
 	else if(loads->MeltingIsPresent()){
-		loads->MeltingConstraints(&num_unstable_constraints,inputs,analysis_type);
+		loads->MeltingConstraints(&num_unstable_constraints,inputs,analysis_type,sub_analysis_type);
 	}
 	else{
Index: issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h
===================================================================
--- issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h	(revision 464)
+++ issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h	(revision 465)
@@ -11,5 +11,5 @@
 /* local prototypes: */
 void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes,
-		DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type); 
+		DataSet* loads,DataSet* materials, ParameterInputs* inputs,int analysis_type,int sub_analysis_type); 
 
 #endif  /* _PENALTYCONSTRAINTSX_H */
Index: issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp
===================================================================
--- issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp	(revision 464)
+++ issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp	(revision 465)
@@ -14,5 +14,5 @@
 
 void PenaltySystemMatricesx(Mat Kgg, Vec pg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
-		int kflag,int pflag,ParameterInputs* inputs,int analysis_type){
+		int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 	
 	int i;
@@ -33,6 +33,6 @@
 
 	/*Add penalties to stiffnesses, from loads: */
-	if(kflag)loads->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type);
-	if(pflag)loads->PenaltyCreatePVector(pg,inputs,kmax,analysis_type);
+	if(kflag)loads->PenaltyCreateKMatrix(Kgg,inputs,kmax,analysis_type,sub_analysis_type);
+	if(pflag)loads->PenaltyCreatePVector(pg,inputs,kmax,analysis_type,sub_analysis_type);
 	
 	/*Assemble matrices: */
Index: issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h
===================================================================
--- issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h	(revision 464)
+++ issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h	(revision 465)
@@ -11,5 +11,5 @@
 /* local prototypes: */
 void PenaltySystemMatricesx(Mat Kgg, Vec pg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
-		int kflag,int pflag,ParameterInputs* inputs,int analysis_type); 
+		int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); 
 
 #endif  /* _PENALTYSYSTEMMATRICESX_H */
Index: issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp
===================================================================
--- issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp	(revision 464)
+++ issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp	(revision 465)
@@ -54,9 +54,5 @@
 
 	
-	if (   (analysis_type==ControlAnalysisEnum()) || 
-	       (analysis_type==DiagnosticHorizAnalysisEnum()) || 
-	       (analysis_type==DiagnosticVertAnalysisEnum()) || 
-	       (analysis_type==DiagnosticStokesAnalysisEnum()) 
-	   ){
+	if (   (analysis_type==ControlAnalysisEnum()) ||  (analysis_type==DiagnosticAnalysisEnum())){
 
 		parameters->FindParam((void*)&vx,"vx");
Index: issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp	(revision 464)
+++ issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp	(revision 465)
@@ -14,5 +14,5 @@
 
 void SystemMatricesx(Mat* pKgg, Vec* ppg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
-		int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type){
+		int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 	
 	int i;
@@ -40,10 +40,10 @@
 
 	/*Fill stiffness matrix and right hand side vector, from elements: */
-	if(kflag)elements->CreateKMatrix(Kgg,inputs,analysis_type);
-	if(pflag)elements->CreatePVector(pg,inputs,analysis_type);
+	if(kflag)elements->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
+	if(pflag)elements->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
 	
 	/*Fill stiffness matrix and right hand side vector, from loads: */
-	if(kflag)loads->CreateKMatrix(Kgg,inputs,analysis_type);
-	if(pflag)loads->CreatePVector(pg,inputs,analysis_type);
+	if(kflag)loads->CreateKMatrix(Kgg,inputs,analysis_type,sub_analysis_type);
+	if(pflag)loads->CreatePVector(pg,inputs,analysis_type,sub_analysis_type);
 
 	/*Assemble matrices: */
Index: issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h
===================================================================
--- issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h	(revision 464)
+++ issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h	(revision 465)
@@ -11,5 +11,5 @@
 /* local prototypes: */
 void SystemMatricesx(Mat* pKgg, Vec* ppg,DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
-		int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type); 
+		int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); 
 
 #endif  /* _SYSTEMMATRICESX_H */
Index: issm/trunk/src/c/io/SerialFetchData.cpp
===================================================================
--- issm/trunk/src/c/io/SerialFetchData.cpp	(revision 464)
+++ issm/trunk/src/c/io/SerialFetchData.cpp	(revision 465)
@@ -27,12 +27,12 @@
 	int      outdataset_size;
 
-	double*  outmatrix=NULL;
-	Mat      outpetscmatrix=NULL;
+	double*  outmatrix;
+	Mat      outpetscmatrix;
 	double*  outmatrix_workspace=NULL;;
 	int      outmatrix_rows,outmatrix_cols;
 	int      petsc;
 
-	double*  outvector=NULL;
-	Vec      outpetscvector=NULL;
+	double*  outvector;
+	Vec      outpetscvector;
 	double*  outvector_workspace=NULL;
 	int      outvector_rows;
Index: issm/trunk/src/c/objects/Beam.cpp
===================================================================
--- issm/trunk/src/c/objects/Beam.cpp	(revision 464)
+++ issm/trunk/src/c/objects/Beam.cpp	(revision 465)
@@ -212,11 +212,15 @@
 #define __FUNCT__ "Beam::CreateKMatrix"
 
-void  Beam::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+void  Beam::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==DiagnosticHutterAnalysisEnum()) {
-
-		CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type);
-
+	if (analysis_type==DiagnosticAnalysisEnum()) {
+	
+		if (sub_analysis_type==HutterAnalysisEnum()) {
+
+			CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type);
+		}
+		else 
+			throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis_type: ",sub_analysis_type," not supported yet"));
 	}
 	else{
@@ -230,5 +234,5 @@
 #define __FUNCT__ "Beam::CreateKMatrixDiagnosticHutter"
 
-void  Beam::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type){
+void  Beam::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 	
 	
@@ -265,9 +269,13 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::CreatePVector"
-void  Beam::CreatePVector(Vec pg,void* inputs,int analysis_type){
+void  Beam::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
 	
 	/*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==DiagnosticHutterAnalysisEnum()) {
-		CreatePVectorDiagnosticHutter( pg,inputs,analysis_type);
+	if (analysis_type==DiagnosticAnalysisEnum()) {
+		if (sub_analysis_type==HutterAnalysisEnum()) {
+			CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type);
+		}
+		else
+			throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
 	}
 	else{
@@ -281,5 +289,5 @@
 #define __FUNCT__ "Beam::CreatePVectorDiagnosticHutter"
 
-void Beam::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type){
+void Beam::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
 
 	int i,j,k;
@@ -534,25 +542,25 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::Du"
-void  Beam::Du(_p_Vec*, double*, double*, void*, int){
+void  Beam::Du(_p_Vec*, double*, double*, void*, int,int){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::Gradj"
-void  Beam::Gradj(_p_Vec*, double*, double*, void*, int, char*){
+void  Beam::Gradj(_p_Vec*, double*, double*, void*, int, int,char*){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::GradjDrag"
-void  Beam::GradjDrag(_p_Vec*, double*, double*, void*, int){
+void  Beam::GradjDrag(_p_Vec*, double*, double*, void*, int,int ){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::GradjB"
-void  Beam::GradjB(_p_Vec*, double*, double*, void*, int){
+void  Beam::GradjB(_p_Vec*, double*, double*, void*, int, int){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::Misfit"
-double Beam::Misfit(double*, double*, void*, int){
+double Beam::Misfit(double*, double*, void*, int,int ){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
Index: issm/trunk/src/c/objects/Beam.h
===================================================================
--- issm/trunk/src/c/objects/Beam.h	(revision 464)
+++ issm/trunk/src/c/objects/Beam.h	(revision 465)
@@ -56,12 +56,12 @@
 		int   MyRank();
 		void  Configure(void* loads,void* nodes,void* materials);
-		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
-		void  CreatePVector(Vec pg, void* inputs, int analysis_type);
+		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
 		void  UpdateFromInputs(void* inputs);
 		void  GetDofList(int* doflist,int* pnumberofdofs);
 		void  GetDofList1(int* doflist);
-		void  CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type);
+		void  CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
 		void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
-		void  CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type);
+		void  CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
 		void* GetMatPar();
 
@@ -78,9 +78,9 @@
 		void  GetBedList(double*);
 		void  GetThicknessList(double* thickness_list);
-		void  Du(_p_Vec*, double*, double*, void*, int);
-		void  Gradj(_p_Vec*, double*, double*, void*, int, char*);
-		void  GradjDrag(_p_Vec*, double*, double*, void*, int);
-		void  GradjB(_p_Vec*, double*, double*, void*, int);
-		double Misfit(double*, double*, void*, int);
+		void  Du(_p_Vec*, double*, double*, void*, int,int);
+		void  Gradj(_p_Vec*, double*, double*, void*, int, int,char*);
+		void  GradjDrag(_p_Vec*, double*, double*, void*, int,int );
+		void  GradjB(_p_Vec*, double*, double*, void*, int,int );
+		double Misfit(double*, double*, void*, int,int );
 		void  GetNodalFunctions(double* l1l2, double gauss_coord);
 		void  GetParameterValue(double* pvalue, double* value_list,double gauss_coord);
Index: issm/trunk/src/c/objects/Element.h
===================================================================
--- issm/trunk/src/c/objects/Element.h	(revision 464)
+++ issm/trunk/src/c/objects/Element.h	(revision 465)
@@ -24,6 +24,6 @@
 		virtual void  Demarshall(char** pmarshalled_dataset)=0;
 		virtual void  Configure(void* loads,void* nodes,void* materials)=0;
-		virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0;
-		virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type)=0;
+		virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0;
+		virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0;
 		virtual void  UpdateFromInputs(void* inputs)=0;
 		virtual void  GetNodes(void** nodes)=0;
@@ -33,9 +33,9 @@
 		virtual void  GetThicknessList(double* thickness_list)=0;
 		virtual void  GetBedList(double* bed_list)=0;
-		virtual void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0;
-		virtual void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type)=0;
-		virtual void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0;
-		virtual void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0;
-        virtual double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0;
+		virtual void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type)=0;
+		virtual void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type)=0;
+		virtual void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type)=0;
+		virtual void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type)=0;
+        virtual double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type)=0;
         virtual void  ComputePressure(Vec p_g)=0;
 
Index: issm/trunk/src/c/objects/Icefront.cpp
===================================================================
--- issm/trunk/src/c/objects/Icefront.cpp	(revision 464)
+++ issm/trunk/src/c/objects/Icefront.cpp	(revision 465)
@@ -199,5 +199,5 @@
 	return my_rank; 
 }
-void  Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
+void  Icefront::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){return;}
 		
 #undef __FUNCT__ 
@@ -240,5 +240,5 @@
 #define __FUNCT__ "Icefront::CreateKMatrix"
 
-void  Icefront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+void  Icefront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*No stiffness loads applied, do nothing: */
@@ -250,13 +250,26 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::CreatePVector"
-void  Icefront::CreatePVector(Vec pg, void* inputs, int analysis_type){
+void  Icefront::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
 
 	/*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
-	if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
-		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
-	}
-	else if ((analysis_type==DiagnosticStokesAnalysisEnum())){
-
-		CreatePVectorDiagnosticStokes( pg,inputs,analysis_type);
+	if (analysis_type==ControlAnalysisEnum()){
+		
+		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
+
+	}
+	else if (analysis_type==DiagnosticAnalysisEnum()){
+	
+		if (sub_analysis_type==HorizAnalysisEnum()){
+		
+			CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
+
+		}
+		else if (sub_analysis_type==StokesAnalysisEnum()){
+			
+			CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
+
+		}
+		else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," sub_analysis_type: ",sub_analysis_type," not supported yet"));
+
 	}
 	else{
@@ -268,12 +281,12 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::CreatePVectorDiagnosticHoriz"
-void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type){
+void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
 
 	/*Branck on the type of icefront: */
 	if (strcmp(type,"segment")==0){
-		CreatePVectorDiagnosticHorizSegment(pg,inputs,analysis_type);
+		CreatePVectorDiagnosticHorizSegment(pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else{
-		CreatePVectorDiagnosticHorizQuad(pg,inputs,analysis_type);
+		CreatePVectorDiagnosticHorizQuad(pg,inputs,analysis_type,sub_analysis_type);
 	}
 }	
@@ -282,5 +295,5 @@
 #define __FUNCT__ "Icefront::CreatePVectorDiagnosticHorizSegment"
 
-void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,void* vinputs, int analysis_type){
+void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){
 
 	int i,j;
@@ -395,5 +408,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefont::CreatePVectorDiagnosticHorizQuad"
-void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,void* vinputs, int analysis_type){
+void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){
 
 	int i,j;
@@ -557,5 +570,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefont::CreatePVectorDiagnosticStokes"
-void Icefront::CreatePVectorDiagnosticStokes( Vec pg,void* vinputs, int analysis_type){
+void Icefront::CreatePVectorDiagnosticStokes( Vec pg,void* vinputs, int analysis_type,int sub_analysis_type){
 
 	int i,j;
@@ -853,5 +866,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::PenaltyCreateKMatrix"
-void  Icefront::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
+void  Icefront::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
 	/*do nothing: */
 }
@@ -859,5 +872,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Icefront::PenaltyCreatePVector"
-void  Icefront::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
+void  Icefront::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
 	/*do nothing: */
 }
Index: issm/trunk/src/c/objects/Icefront.h
===================================================================
--- issm/trunk/src/c/objects/Icefront.h	(revision 464)
+++ issm/trunk/src/c/objects/Icefront.h	(revision 465)
@@ -55,14 +55,14 @@
 		int   GetId(); 
 		int   MyRank();
-		void  DistributeNumDofs(int* numdofspernode,int analysis_type);
+		void  DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type);
 		void  Configure(void* elements,void* nodes,void* materials);
-		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
+		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
 		void  UpdateFromInputs(void* inputs);
 		
-		void  CreatePVector(Vec pg, void* inputs, int analysis_type);
-		void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type);
-		void  CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type);
-		void  CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type);
-		void  CreatePVectorDiagnosticStokes( Vec pg,void* inputs, int analysis_type);
+		void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
+		void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
+		void  CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type,int sub_analysis_type);
+		void  CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type,int sub_analysis_type);
+		void  CreatePVectorDiagnosticStokes( Vec pg,void* inputs, int analysis_type,int sub_analysis_type);
 		void  GetDofList(int* doflist,int* pnumberofdofs);
 		void  SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill);
@@ -71,6 +71,6 @@
 		void  QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 
 		                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill);
-		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
-		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
+		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
+		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
 		Object* copy();
 };
Index: issm/trunk/src/c/objects/Load.h
===================================================================
--- issm/trunk/src/c/objects/Load.h	(revision 464)
+++ issm/trunk/src/c/objects/Load.h	(revision 465)
@@ -24,9 +24,9 @@
 		virtual void  Demarshall(char** pmarshalled_dataset)=0;
 		virtual void  Configure(void* elements,void* nodes,void* materials)=0;
-		virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0;
-		virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type)=0;
+		virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0;
+		virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0;
 		virtual void  UpdateFromInputs(void* inputs)=0;
-		virtual void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type)=0;
-		virtual void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type)=0;
+		virtual void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type)=0;
+		virtual void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type)=0;
 
 		int           Enum();
Index: issm/trunk/src/c/objects/OptArgs.h
===================================================================
--- issm/trunk/src/c/objects/OptArgs.h	(revision 464)
+++ issm/trunk/src/c/objects/OptArgs.h	(revision 465)
@@ -20,4 +20,5 @@
 	mxArray* n;
 	mxArray* analysis_type;
+	mxArray* sub_analysis_type;
 
 };
Index: issm/trunk/src/c/objects/Pengrid.cpp
===================================================================
--- issm/trunk/src/c/objects/Pengrid.cpp	(revision 464)
+++ issm/trunk/src/c/objects/Pengrid.cpp	(revision 465)
@@ -135,5 +135,5 @@
 	return my_rank; 
 }
-void  Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type){return;}
+void  Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
 
 #undef __FUNCT__ 
@@ -156,5 +156,5 @@
 #define __FUNCT__ "Pengrid::CreateKMatrix"
 
-void  Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+void  Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*No loads applied, do nothing: */
@@ -165,5 +165,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Pengrid::CreatePVector"
-void  Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type){
+void  Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
 
 	/*No loads applied, do nothing: */
@@ -179,13 +179,13 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Pengrid::PenaltyCreateKMatrix"
-void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
-
-	if ((analysis_type==DiagnosticStokesAnalysisEnum())){
-
-		PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type);
+void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+
+	if ((analysis_type==DiagnosticAnalysisEnum()) && ((sub_analysis_type==StokesAnalysisEnum()))){
+
+		PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type);
 
 	}
 	else{
-		throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
+		throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
 	}
 
@@ -194,5 +194,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes"
-void  Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type){
+void  Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
 	
 	const int numgrids=1;
Index: issm/trunk/src/c/objects/Pengrid.h
===================================================================
--- issm/trunk/src/c/objects/Pengrid.h	(revision 464)
+++ issm/trunk/src/c/objects/Pengrid.h	(revision 465)
@@ -39,12 +39,12 @@
 		int   GetId(); 
 		int   MyRank();
-		void  DistributeNumDofs(int* numdofspernode,int analysis_type);
+		void  DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type);
 		void  Configure(void* elements,void* nodes,void* materials);
-		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
-		void  CreatePVector(Vec pg, void* inputs, int analysis_type);
+		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
 		void  UpdateFromInputs(void* inputs);
-		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
-		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
-		void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* inputs,double kmax,int analysis_type);
+		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
+		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
+		void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
 		void  GetDofList(int* doflist,int* pnumberofdofspernode);
 		Object* copy();
Index: issm/trunk/src/c/objects/Penpair.cpp
===================================================================
--- issm/trunk/src/c/objects/Penpair.cpp	(revision 464)
+++ issm/trunk/src/c/objects/Penpair.cpp	(revision 465)
@@ -182,5 +182,5 @@
 	return my_rank; 
 }
-void  Penpair::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
+void  Penpair::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){return;}
 
 #undef __FUNCT__ 
@@ -208,5 +208,5 @@
 #define __FUNCT__ "Penpair::CreateKMatrix"
 
-void  Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+void  Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*No loads applied, do nothing: */
@@ -217,5 +217,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penpair::CreatePVector"
-void  Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type){
+void  Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
 
 	/*No loads applied, do nothing: */
@@ -231,5 +231,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penpair::PenaltyCreateKMatrix"
-void  Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
+void  Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
 	if(numdofs==1){
 		
@@ -245,5 +245,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penpair::PenaltyCreatePVector"
-void  Penpair::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
+void  Penpair::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
 	if(numdofs==1){
 		
Index: issm/trunk/src/c/objects/Penpair.h
===================================================================
--- issm/trunk/src/c/objects/Penpair.h	(revision 464)
+++ issm/trunk/src/c/objects/Penpair.h	(revision 465)
@@ -51,11 +51,11 @@
 		int   GetId(); 
 		int   MyRank();
-		void  DistributeNumDofs(int* numdofspernode,int analysis_type);
+		void  DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type);
 		void  Configure(void* elements,void* nodes,void* materials);
-		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
-		void  CreatePVector(Vec pg, void* inputs, int analysis_type);
+		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
 		void  UpdateFromInputs(void* inputs);
-		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
-		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
+		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
+		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
 		Object* copy();
 };
Index: issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Penta.cpp	(revision 464)
+++ issm/trunk/src/c/objects/Penta.cpp	(revision 465)
@@ -284,30 +284,32 @@
 #define __FUNCT__ "Penta::CreateKMatrix"
 
-void  Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+void  Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
-
-		CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type);
-
-	}
-	else if ((analysis_type==DiagnosticVertAnalysisEnum())){
-
-		CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type);
-
-	}
-	else if ((analysis_type==DiagnosticStokesAnalysisEnum())){
-
-		CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type);
-
-	}
-	else if (
-			(analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || 
-			(analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || 
-			(analysis_type==BedSlopeComputeXAnalysisEnum()) || 
-			(analysis_type==BedSlopeComputeYAnalysisEnum()) 
-			){
-		
-		CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type);
+	if (analysis_type==ControlAnalysisEnum()){
+
+		CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
+
+	}
+	else if (analysis_type==DiagnosticAnalysisEnum()){
+	
+		if (sub_analysis_type==HorizAnalysisEnum()){
+		
+			CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
+		}
+		else if (sub_analysis_type==VertAnalysisEnum()){
+		
+			CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type);
+		}
+		else if (sub_analysis_type==StokesAnalysisEnum()){
+		
+			CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type);
+		
+		}
+		else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
+	}
+	else if (analysis_type==SlopeComputeAnalysisEnum()){
+		
+		CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
 	}
 	else{
@@ -320,5 +322,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta:CreateKMatrixDiagnosticHoriz"
-void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type){
+void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){
 
 
@@ -421,5 +423,5 @@
 		 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-		tria->CreateKMatrix(Kgg,inputs, analysis_type);
+		tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
 		delete tria;
 		return;
@@ -552,5 +554,5 @@
 
 			tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-			tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type);
+			tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type);
 			delete tria;
 		}
@@ -573,5 +575,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta:CreateKMatrixDiagnosticVert"
-void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type){
+void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){
 
 	/* local declarations */
@@ -623,5 +625,5 @@
 	if(onsurface){
 		tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
-		tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type);
+		tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type,sub_analysis_type);
 		delete tria;
 	}
@@ -710,5 +712,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta:CreateKMatrixDiagnosticStokes"
-void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type){
+void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){
 
 	int i,j;
@@ -985,27 +987,30 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::CreatePVector"
-void  Penta::CreatePVector(Vec pg, void* inputs, int analysis_type){
+void  Penta::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
 
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
-
-		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
-	}
-	else if ((analysis_type==DiagnosticVertAnalysisEnum())){
-
-		CreatePVectorDiagnosticVert( pg,inputs,analysis_type);
-	}
-	else if ((analysis_type==DiagnosticStokesAnalysisEnum())){
-
-		CreatePVectorDiagnosticStokes( pg,inputs,analysis_type);
-	}
-	else if (
-			(analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || 
-			(analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || 
-			(analysis_type==BedSlopeComputeXAnalysisEnum()) || 
-			(analysis_type==BedSlopeComputeYAnalysisEnum()) 
-			){
-		
-		CreatePVectorSlopeCompute( pg,inputs,analysis_type);
+	if (analysis_type==ControlAnalysisEnum()){
+		
+		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
+	}
+	else if (analysis_type==DiagnosticAnalysisEnum()){
+	
+		if (sub_analysis_type==HorizAnalysisEnum()){
+
+			CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
+		}
+		else if (sub_analysis_type==VertAnalysisEnum()){
+		
+			CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type);
+		}
+		else if (sub_analysis_type==StokesAnalysisEnum()){
+		
+			CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
+		}
+		else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
+	}
+	else if (analysis_type==SlopeComputeAnalysisEnum()){
+		
+		CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else{
@@ -1091,5 +1096,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::Du"
-void  Penta::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type){
+void  Penta::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -1122,11 +1127,11 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::Gradj"
-void  Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){
+void  Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
 
 	if (strcmp(control_type,"drag")==0){
-		GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type);
+		GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (strcmp(control_type,"B")==0){
-		GradjB( grad_g, velocity, adjoint, inputs,analysis_type);
+		GradjB( grad_g, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
 	}
 	else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type));
@@ -1135,5 +1140,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GradjDrag"
-void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){
+void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){
 	
 	
@@ -1147,5 +1152,5 @@
 		
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-		tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type);
+		tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type);
 		delete tria;
 		return;
@@ -1155,5 +1160,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::GradjB"
-void  Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){
+void  Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
@@ -1161,5 +1166,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::Misfit"
-double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type){
+double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type,int sub_analysis_type){
 	
 	double J;
@@ -1641,5 +1646,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::CreatePVectorDiagnosticHoriz"
-void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type){
+void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int i,j;
@@ -1707,5 +1712,5 @@
 		 *and use its CreatePVector functionality to return an elementary load vector: */
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-		tria->CreatePVector(pg,inputs, analysis_type);
+		tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
 		delete tria;
 		return;
@@ -1997,5 +2002,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta:CreatePVectorDiagnosticVert"
-void  Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type){
+void  Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -2055,5 +2060,5 @@
 	if(onbed){
 		tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
-		tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
+		tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type,sub_analysis_type);
 		delete tria;
 	}
@@ -2176,5 +2181,5 @@
 #define __FUNCT__ "Penta::CreateKMatrixSlopeCompute"
 
-void  Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type){
+void  Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*Collapsed formulation: */
@@ -2186,5 +2191,5 @@
 	/*Spawn Tria element from the base of the Penta: */
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-	tria->CreateKMatrix(Kgg,inputs, analysis_type);
+	tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type);
 	delete tria;
 	return;
@@ -2195,5 +2200,5 @@
 #define __FUNCT__ "Penta::CreatePVectorSlopeCompute"
 
-void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type){
+void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
 	
 	/*Collapsed formulation: */
@@ -2205,5 +2210,5 @@
 	/*Spawn Tria element from the base of the Penta: */
 	tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-	tria->CreatePVector(pg,inputs, analysis_type);
+	tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type);
 	delete tria;
 	return;
@@ -2826,5 +2831,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::CreatePVectorDiagnosticStokes"
-void Penta::CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type){
+void Penta::CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){
 
 
Index: issm/trunk/src/c/objects/Penta.h
===================================================================
--- issm/trunk/src/c/objects/Penta.h	(revision 464)
+++ issm/trunk/src/c/objects/Penta.h	(revision 465)
@@ -73,8 +73,8 @@
 		int   MyRank();
 		void  Configure(void* loads,void* nodes,void* materials);
-		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
-		void  CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type);
-		void  CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type);
-		void  CreatePVector(Vec pg, void* inputs, int analysis_type);
+		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type,int sub_analysis_type);
+		void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
 		void  UpdateFromInputs(void* inputs);
 		void  GetDofList(int* doflist,int* pnumberofdofs);
@@ -84,9 +84,9 @@
 		void  GetNodes(void** nodes);
 		int   GetOnBed();
-		void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type);
-		void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
-		void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
-		void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
-        double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type);
+		void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
+		void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
+		void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type);
+        double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
 		
 		void          GetThicknessList(double* thickness_list);
@@ -105,6 +105,6 @@
 		void  GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4);
 		void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4);
-		void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type);
-		void  CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type);
+		void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type,int sub_analysis_type);
 		void  GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4);
 		void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);
@@ -113,9 +113,9 @@
 		void  SlopeExtrude(Vec sg,double* sg_serial);
 		void  ComputePressure(Vec p_g);
-		void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type);
-		void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type);
+		void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
 
-		void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type);
-		void  CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type);
+		void  CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type);
+		void  CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type);
 		void  ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
 		void  GetMatrixInvert(double*  Ke_invert, double* Ke);
Index: issm/trunk/src/c/objects/Sing.cpp
===================================================================
--- issm/trunk/src/c/objects/Sing.cpp	(revision 464)
+++ issm/trunk/src/c/objects/Sing.cpp	(revision 465)
@@ -193,14 +193,14 @@
 #define __FUNCT__ "Sing::CreateKMatrix"
 
-void  Sing::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+void  Sing::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==DiagnosticHutterAnalysisEnum()){
-
-		CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type);
+	if ((analysis_type==DiagnosticAnalysisEnum()) && (sub_analysis_type==HutterAnalysisEnum())){
+
+		CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type);
 
 	}
 	else{
-		throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
+		throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s\n","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet"));
 	}
 
@@ -211,5 +211,5 @@
 #define __FUNCT__ "Sing::CreateKMatrixDiagnosticHutter"
 
-void  Sing::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type){
+void  Sing::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 	
 	const int numgrids=1;
@@ -234,9 +234,11 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Sing::CreatePVector"
-void  Sing::CreatePVector(Vec pg,void* inputs,int analysis_type){
+void  Sing::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
 	
 	/*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==DiagnosticHutterAnalysisEnum()){
-		CreatePVectorDiagnosticHutter( pg,inputs,analysis_type);
+	if ((analysis_type==DiagnosticAnalysisEnum()) && (sub_analysis_type==HutterAnalysisEnum())){
+	
+			CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type);
+
 	}
 	else{
@@ -251,5 +253,5 @@
 #define __FUNCT__ "Sing::CreatePVectorDiagnosticHutter"
 
-void Sing::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type){
+void Sing::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
 	
 	
@@ -427,25 +429,25 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Sing::Du"
-void  Sing::Du(_p_Vec*, double*, double*, void*, int){
+void  Sing::Du(_p_Vec*, double*, double*, void*, int,int){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
 #undef __FUNCT__ 
 #define __FUNCT__ "Sing::Gradj"
-void  Sing::Gradj(_p_Vec*, double*, double*, void*, int, char*){
+void  Sing::Gradj(_p_Vec*, double*, double*, void*, int, int ,char*){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
 #undef __FUNCT__ 
 #define __FUNCT__ "Sing::GradjDrag"
-void  Sing::GradjDrag(_p_Vec*, double*, double*, void*, int){
+void  Sing::GradjDrag(_p_Vec*, double*, double*, void*, int,int){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
 #undef __FUNCT__ 
 #define __FUNCT__ "Sing::GradjB"
-void  Sing::GradjB(_p_Vec*, double*, double*, void*, int){
+void  Sing::GradjB(_p_Vec*, double*, double*, void*, int,int){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
 #undef __FUNCT__ 
 #define __FUNCT__ "Sing::Misfit"
-double Sing::Misfit(double*, double*, void*, int){
+double Sing::Misfit(double*, double*, void*, int,int){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
Index: issm/trunk/src/c/objects/Sing.h
===================================================================
--- issm/trunk/src/c/objects/Sing.h	(revision 464)
+++ issm/trunk/src/c/objects/Sing.h	(revision 465)
@@ -51,12 +51,12 @@
 		int   MyRank();
 		void  Configure(void* loads,void* nodes,void* materials);
-		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
-		void  CreatePVector(Vec pg, void* inputs, int analysis_type);
+		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
 		void  UpdateFromInputs(void* inputs);
 		void  GetDofList(int* doflist,int* pnumberofdofs);
 		void  GetDofList1(int* doflist);
-		void  CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type);
+		void  CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
 		void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
-		void  CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type);
+		void  CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
 		void* GetMatPar();
 
@@ -73,9 +73,9 @@
 		void  GetBedList(double*);
 		void  GetThicknessList(double* thickness_list);
-		void  Du(_p_Vec*, double*, double*, void*, int);
-		void  Gradj(_p_Vec*, double*, double*, void*, int, char*);
-		void  GradjDrag(_p_Vec*, double*, double*, void*, int);
-		void  GradjB(_p_Vec*, double*, double*, void*, int);
-		double Misfit(double*, double*, void*, int);
+		void  Du(_p_Vec*, double*, double*, void*, int,int);
+		void  Gradj(_p_Vec*, double*, double*, void*, int, int,char*);
+		void  GradjDrag(_p_Vec*, double*, double*, void*, int,int);
+		void  GradjB(_p_Vec*, double*, double*, void*, int,int);
+		double Misfit(double*, double*, void*, int,int);
 
 
Index: issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Tria.cpp	(revision 464)
+++ issm/trunk/src/c/objects/Tria.cpp	(revision 465)
@@ -266,20 +266,24 @@
 #define __FUNCT__ "Tria::CreateKMatrix"
 
-void  Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+void  Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
-
-		CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type);
-
-	}
-	else if (
-			(analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || 
-			(analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || 
-			(analysis_type==BedSlopeComputeXAnalysisEnum()) || 
-			(analysis_type==BedSlopeComputeYAnalysisEnum()) 
-			){
-		
-		CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type);
+	if (analysis_type==ControlAnalysisEnum()){
+		
+		CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
+	}
+	else if (analysis_type==DiagnosticAnalysisEnum()){
+	
+		if (sub_analysis_type==HorizAnalysisEnum()){
+
+			CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
+		}
+		else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
+
+	}
+	else if (analysis_type==SlopeComputeAnalysisEnum()){
+
+		CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
+
 	}
 	else{
@@ -293,5 +297,5 @@
 #define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz"
 
-void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type){
+void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 
@@ -477,5 +481,5 @@
 	/*Do not forget to include friction: */
 	if(!shelf){
-		CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type);
+		CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type);
 	}
 
@@ -508,5 +512,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction"
-void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type){
+void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 
@@ -677,5 +681,5 @@
 #define __FUNCT__ "Tria::CreateKMatrixSlopeCompute"
 
-void  Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type){
+void  Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	/* local declarations */
@@ -758,18 +762,23 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::CreatePVector"
-void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type){
+void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
 	
 	/*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
-	if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
-		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
-	}
-	else if (
-			(analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || 
-			(analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || 
-			(analysis_type==BedSlopeComputeXAnalysisEnum()) || 
-			(analysis_type==BedSlopeComputeYAnalysisEnum()) 
-			){
-		
-		CreatePVectorSlopeCompute( pg,inputs,analysis_type);
+	if (analysis_type==ControlAnalysisEnum()){
+		
+		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
+	
+	}
+	else if (analysis_type==DiagnosticAnalysisEnum()){
+		if (sub_analysis_type==HorizAnalysisEnum()){
+		
+			CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
+		
+		}
+		else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
+	}
+	else if (analysis_type==SlopeComputeAnalysisEnum()){
+		
+		CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else{
@@ -779,10 +788,8 @@
 }
 
-
-
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz"
 
-void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type){
+void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
 
 	int             i,j;
@@ -952,5 +959,5 @@
 #define __FUNCT__ "Tria::CreatePVectorSlopeCompute"
 
-void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type){
+void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
 
 	int             i,j;
@@ -992,8 +999,8 @@
 	for(i=0;i<numdofs;i++) pe_g[i]=0.0;
 
-	if ( (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || (analysis_type==SurfaceSlopeComputeYAnalysisEnum())){
+	if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==SurfaceYAnalysisEnum())){
 		for(i=0;i<numdofs;i++) param[i]=s[i];
 	}
-	if ( (analysis_type==BedSlopeComputeXAnalysisEnum()) || (analysis_type==BedSlopeComputeYAnalysisEnum())){
+	if ( (sub_analysis_type==BedXAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
 		for(i=0;i<numdofs;i++) param[i]=b[i];
 	}
@@ -1020,8 +1027,8 @@
 
 		/*Build pe_g_gaussian vector: */
-		if ( (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || (analysis_type==BedSlopeComputeXAnalysisEnum())){
+		if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==BedXAnalysisEnum())){
 			for(i=0;i<numdofs;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[0]*l1l2l3[i];
 		}
-		if ( (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || (analysis_type==BedSlopeComputeYAnalysisEnum())){
+		if ( (sub_analysis_type==SurfaceYAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
 			for(i=0;i<numdofs;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[1]*l1l2l3[i];
 		}
@@ -1543,5 +1550,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::Du"
-void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
+void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -1745,11 +1752,11 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::Gradj"
-void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){
+void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
 
 	if (strcmp(control_type,"drag")==0){
-		GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type);
+		GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (strcmp(control_type,"B")==0){
-		GradjB( grad_g, velocity, adjoint, inputs,analysis_type);
+		GradjB( grad_g, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
 	}
 	else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type));
@@ -1758,5 +1765,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::GradjDrag"
-void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){
+void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type,int sub_analysis_type){
 
 
@@ -1952,5 +1959,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::GradjB"
-void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){
+void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -2093,5 +2100,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::Misfit"
-double Tria::Misfit(double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
+double Tria::Misfit(double* velocity,double* obs_velocity,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int i;
@@ -2277,5 +2284,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
-void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type){
+void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int i,j;
@@ -2405,5 +2412,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::CreatePVectorDiagnosticBaseVert"
-void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,void* vinputs,int analysis_type){
+void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	int             i,j;
Index: issm/trunk/src/c/objects/Tria.h
===================================================================
--- issm/trunk/src/c/objects/Tria.h	(revision 464)
+++ issm/trunk/src/c/objects/Tria.h	(revision 465)
@@ -65,13 +65,13 @@
 		int   MyRank();
 		void  Configure(void* loads,void* nodes,void* materials);
-		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
-		void  CreatePVector(Vec pg, void* inputs, int analysis_type);
+		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
 		void  UpdateFromInputs(void* inputs);
 		void  GetDofList(int* doflist,int* pnumberofdofs);
 		void  GetDofList1(int* doflist);
-		void  CreateKMatrixDiagnosticHoriz(Mat Kgg,void* inputs,int analysis_type);
-		void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* inputs,int analysis_type);
-		void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* inputs,int analysis_type);
-		void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type);
+		void  CreateKMatrixDiagnosticHoriz(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
 		void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
 		void  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
@@ -87,13 +87,13 @@
 		void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3);
 		void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);
-		void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type);
-		void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
-		void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
-		void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
-        double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type);
+		void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
+		void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type);
+		void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type);
+		void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type);
+        double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
 
-		void  CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type);
-		void  CreatePVectorDiagnosticBaseVert(Vec pg,void* inputs,int analysis_type);
-		void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type);
+		void  CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorDiagnosticBaseVert(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
 		void* GetMatPar();
 		int   GetShelf();
Index: issm/trunk/src/c/objects/cielo/PenpairLoad.c
===================================================================
--- issm/trunk/src/c/objects/cielo/PenpairLoad.c	(revision 464)
+++ issm/trunk/src/c/objects/cielo/PenpairLoad.c	(revision 465)
@@ -689,5 +689,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PenpairLoadPenaltyCreateKMatrix"
-int PenpairLoadPenaltyCreateKMatrix( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){
+int PenpairLoadPenaltyCreateKMatrix( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type,int sub_analysis_type){
 
 	
@@ -702,8 +702,8 @@
 
 	if (this->penpair.numdofs==1){
-		PenpairLoadPenaltyCreateKMatrixOneDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type);
+		PenpairLoadPenaltyCreateKMatrixOneDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type,sub_analysis_type);
 	}
 	else if (this->penpair.numdofs==2){
-		PenpairLoadPenaltyCreateKMatrixTwoDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type);
+		PenpairLoadPenaltyCreateKMatrixTwoDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type,sub_analysis_type);
 	}
 	else{
@@ -719,5 +719,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PenpairLoadPenaltyCreateKMatrixOneDof"
-int PenpairLoadPenaltyCreateKMatrixOneDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type){
+int PenpairLoadPenaltyCreateKMatrixOneDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type,int sub_analysis_type){
 
 
@@ -769,5 +769,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PenpairLoadPenaltyCreateKMatrixTwoDof"
-int PenpairLoadPenaltyCreateKMatrixTwoDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){
+int PenpairLoadPenaltyCreateKMatrixTwoDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type,int sub_analysis_type){
 
 
@@ -916,5 +916,5 @@
 #define __FUNCT__ "PenpairLoadPenaltyCreatePVector"
 
-int PenpairLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){
+int PenpairLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type,int sub_analysis_type){
 
 	PenpairLoad* this   = NULL;
@@ -931,5 +931,5 @@
 	}
 	else if (this->penpair.numdofs==2){
-		PenpairLoadPenaltyCreatePVectorTwoDof( ppe_g,vpthis,inputs,kmax,analysis_type);
+		PenpairLoadPenaltyCreatePVectorTwoDof( ppe_g,vpthis,inputs,kmax,analysis_type,sub_analysis_type);
 	}
 	else{
@@ -947,5 +947,5 @@
 #define __FUNCT__ "PenpairLoadPenaltyCreatePVectorTwoDof"
 
-int PenpairLoadPenaltyCreatePVectorTwoDof( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){
+int PenpairLoadPenaltyCreatePVectorTwoDof( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type,int sub_analysis_type){
 
 	int             i,j;
@@ -1092,5 +1092,4 @@
 	return noerr;
 }
-int PotentialUnstableConstraints(DataSet* lst,ParameterInputs* inputs,int analysis_type_enum);
 
 /*--------------------------------------------------
@@ -1100,5 +1099,6 @@
 #define __FUNCT__ "PenpairLoadPotentialUnstableConstraint"
 
-int PenpairLoadPotentialUnstableConstraint(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){
+int PotentialUnstableConstraints(DataSet* lst,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
+int PenpairLoadPotentialUnstableConstraint(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	/* vpthis for polymorphic function compatibility */
@@ -1186,5 +1186,5 @@
 #define __FUNCT__ "PenpairLoadPenaltyPreConstrain"
 
-int PenpairLoadPenaltyPreConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){
+int PenpairLoadPenaltyPreConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	/* vpthis for polymorphic function compatibility */
@@ -1290,5 +1290,5 @@
 #define __FUNCT__ "PenpairLoadPenaltyConstrain"
 
-int PenpairLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){
+int PenpairLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	/* vpthis for polymorphic function compatibility */
@@ -1417,5 +1417,5 @@
 #define __FUNCT__ "PenpairLoadPenetration"
 
-int PenpairLoadPenetration(double* ppenetration, void* vpthis,ParameterInputs* inputs, int analysis_type){
+int PenpairLoadPenetration(double* ppenetration, void* vpthis,ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	/* vpthis for polymorphic function compatibility */
Index: issm/trunk/src/c/objects/cielo/PentaElement.c
===================================================================
--- issm/trunk/src/c/objects/cielo/PentaElement.c	(revision 464)
+++ issm/trunk/src/c/objects/cielo/PentaElement.c	(revision 465)
@@ -818,5 +818,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PentaElementCreateKMatrix"
-int PentaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type){
+int PentaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type,int sub_analysis_type){
 
 
@@ -2564,5 +2564,5 @@
 #define __FUNCT__ "PentaElementCreatePVector"
 
-int PentaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type){
+int PentaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	int noerr=1;
@@ -5002,5 +5002,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PentaElementCreateDuVector"
-int PentaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){
+int PentaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	int noerr=1;
@@ -5077,5 +5077,5 @@
 
 		/*Ok, now triaelement is correctly configured, call on its method CreateDuVector: */
-		noerr=TriaElementCreateDuVector( pdue_g, (void*)triaelement, velocity, obs_velocity, inputs, analysis_type);
+		noerr=TriaElementCreateDuVector( pdue_g, (void*)triaelement, velocity, obs_velocity, inputs, analysis_type,sub_analysis_type);
 
 		/*Now delete tria and triaelement: */
@@ -5092,13 +5092,13 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PentaElementCreateGradVectors"
-int PentaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,char* control_type){
+int PentaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){
 
 	int noerr=1;
 	
 	if strcasecmp_eq(control_type,"drag"){
-		noerr=PentaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type);
+		noerr=PentaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
 	}
 	else if strcasecmp_eq(control_type,"B"){
-		noerr=PentaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type);
+		noerr=PentaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
 	}
 	else{
@@ -5114,5 +5114,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PentaElementMisfit"
-int PentaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){
+int PentaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	int noerr=1;
@@ -5189,5 +5189,5 @@
 
 		/*Ok, now triaelement is correctly configured, call on its method Misfit: */
-		noerr=TriaElementMisfit(pJelem,(void*)triaelement, velocity, obs_velocity, inputs, analysis_type);
+		noerr=TriaElementMisfit(pJelem,(void*)triaelement, velocity, obs_velocity, inputs, analysis_type,sub_analysis_type);
 
 		/*Now delete tria and triaelement: */
@@ -5204,5 +5204,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PentaElementCreateGradVectorsDrag"
-int PentaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){
+int PentaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 	int noerr=1;
@@ -5279,5 +5279,5 @@
 
 		/*Ok, now triaelement is correctly configured, call on its method CreateGradjVectorsdrag: */
-		noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, (void*)triaelement,velocity, adjoint, inputs,analysis_type);
+		noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, (void*)triaelement,velocity, adjoint, inputs,analysis_type,sub_analysis_type);
 
 		/*Now delete tria and triaelement: */
@@ -5295,5 +5295,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "PentaElementCreateGradjVectorsB"
-int PentaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){
+int PentaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 	printf("Not supported yet!\n");
Index: issm/trunk/src/c/objects/cielo/TriaElement.c
===================================================================
--- issm/trunk/src/c/objects/cielo/TriaElement.c	(revision 464)
+++ issm/trunk/src/c/objects/cielo/TriaElement.c	(revision 465)
@@ -777,5 +777,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "TriaElementCreateKMatrix"
-int TriaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type){
+int TriaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type,int sub_analysis_type){
 
 
@@ -800,5 +800,5 @@
 #define __FUNCT__ "TriaElementCreatePVector"
 
-int TriaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type){
+int TriaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	int noerr=1;
@@ -2044,5 +2044,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "TriaElementCreateDuVector"
-int TriaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){
+int TriaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 					
 	int             i,j;
@@ -2263,13 +2263,13 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "TriaElementCreateGradVectors"
-int TriaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,char* control_type){
+int TriaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){
 
 	int noerr=1;
 	
 	if strcasecmp_eq(control_type,"drag"){
-		noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type);
+		noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
 	}
 	else if strcasecmp_eq(control_type,"B"){
-		noerr=TriaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type);
+		noerr=TriaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type,sub_analysis_type);
 	}
 	else{
@@ -2285,5 +2285,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "TriaElementCreateGradVectorsDrag"
-int TriaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){
+int TriaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 	int             i,j;
@@ -2566,5 +2566,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "TriaElementCreateGradjVectorsB"
-int TriaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){
+int TriaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 	int             i,j;
@@ -2770,5 +2770,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "TriaElementMisfit"
-int TriaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){
+int TriaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type,int sub_analysis_type){
 
 	int             i,j;
Index: issm/trunk/src/c/parallel/CreateFemModel.cpp
===================================================================
--- issm/trunk/src/c/parallel/CreateFemModel.cpp	(revision 464)
+++ issm/trunk/src/c/parallel/CreateFemModel.cpp	(revision 465)
@@ -10,5 +10,5 @@
 #include "../issm.h"
 
-void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type){
+void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type){
 
 	/*Model output: */
@@ -37,16 +37,8 @@
 	_printf_("   specifying analysis\n");
 	model->analysis_type=(char*)xmalloc((strlen(analysis_type)+1)*sizeof(char)); strcpy(model->analysis_type,analysis_type);
+	model->sub_analysis_type=(char*)xmalloc((strlen(sub_analysis_type)+1)*sizeof(char)); strcpy(model->sub_analysis_type,sub_analysis_type);
 
-	_printf_("   create elements, nodes and materials:\n");
-	CreateElementsNodesAndMaterials(&elements,&nodes,&materials,model,MODEL);
-
-	_printf_("   create constraints: \n");
-	CreateConstraints(&constraints,model,MODEL);
-	
-	_printf_("   create loads: \n");
-	CreateLoads(&loads,model,MODEL);
-
-	_printf_("   create parameters: \n");
-	CreateParameters(&parameters,model,MODEL);
+	_printf_("   create datasets:\n");
+	CreateDataSets(&elements,&nodes,&materials,&constraints,&loads,&parameters,model,MODEL);
 
 	_printf_("   create degrees of freedom: \n");
Index: issm/trunk/src/c/parallel/GradJCompute.cpp
===================================================================
--- issm/trunk/src/c/parallel/GradJCompute.cpp	(revision 464)
+++ issm/trunk/src/c/parallel/GradJCompute.cpp	(revision 465)
@@ -20,4 +20,5 @@
 	/*intermediary: */
 	int analysis_type;
+	int sub_analysis_type;
 	char* solverstring=NULL;
 	char* control_type=NULL;
@@ -42,13 +43,14 @@
 	/*some parameters:*/
 	femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type");
+	femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
 	femmodel->parameters->FindParam((void*)&solverstring,"solverstring");
 	femmodel->parameters->FindParam((void*)&control_type,"control_type");
 
 	_printf_("%s\n","      recover solution for this stiffness and right hand side:");
-	diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,analysis_type);
+	diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,analysis_type,sub_analysis_type);
 	VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
 
 	_printf_("%s\n","      buid Du, difference between observed velocity and model velocity:");
-	Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,u_g_double,u_g_obs, inputs,analysis_type);
+	Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,u_g_double,u_g_obs, inputs,analysis_type,sub_analysis_type);
 
 	_printf_("%s\n","      reduce adjoint load from g-set to f-set:");
@@ -68,5 +70,5 @@
 
 	Gradjx( &grad_g, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, 
-		u_g_double,lambda_g_double, inputs,analysis_type,control_type);
+		u_g_double,lambda_g_double, inputs,analysis_type,sub_analysis_type,control_type);
 	
 	/*Free ressources:*/
Index: issm/trunk/src/c/parallel/control.cpp
===================================================================
--- issm/trunk/src/c/parallel/control.cpp	(revision 464)
+++ issm/trunk/src/c/parallel/control.cpp	(revision 465)
@@ -26,4 +26,5 @@
 	char* lockname=NULL;
 	int   analysis_type;
+	int   sub_analysis_type;
 
 	/*Intermediary: */
@@ -78,5 +79,5 @@
 	
 	_printf_("read and create finite element model:\n");
-	CreateFemModel(&femmodel,fid,"control");
+	CreateFemModel(&femmodel,fid,"control","");
 
 	/*Recover parameters used throughout the solution:*/
@@ -93,4 +94,5 @@
 	femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs");
 	femmodel.parameters->FindParam((void*)&analysis_type,"analysis_type");
+	femmodel.parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
 	gsize=femmodel.nodes->NumberOfDofs();
 
@@ -156,5 +158,5 @@
 			inputs->Add(control_type,p_g,2,numberofnodes);
 			inputs->Add("fit",fit[n]);
-			diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type);
+			diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
 			OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
 			_printf_("%s\n","      done.");
@@ -175,5 +177,5 @@
 	UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs);
 	
-	diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type);
+	diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
 	
 	_printf_("%s\n","      saving final results...");
Index: issm/trunk/src/c/parallel/diagnostic.cpp
===================================================================
--- issm/trunk/src/c/parallel/diagnostic.cpp	(revision 464)
+++ issm/trunk/src/c/parallel/diagnostic.cpp	(revision 465)
@@ -59,13 +59,13 @@
 	_printf_("read and create finite element model:\n");
 	_printf_("\n   reading diagnostic horiz model data:\n");
-	CreateFemModel(&femmodels[0],fid,"diagnostic_horiz");
+	CreateFemModel(&femmodels[0],fid,"diagnostic","horiz");
 	_printf_("\n   reading diagnostic vert model data:\n");
-	CreateFemModel(&femmodels[1],fid,"diagnostic_vert");
+	CreateFemModel(&femmodels[1],fid,"diagnostic","vert");
 	_printf_("\n   reading diagnostic stokes model data:\n");
-	CreateFemModel(&femmodels[2],fid,"diagnostic_stokes");
+	CreateFemModel(&femmodels[2],fid,"diagnostic","stokes");
 	_printf_("\n   reading diagnostic hutter model data:\n");
-	CreateFemModel(&femmodels[3],fid,"diagnostic_hutter");
+	CreateFemModel(&femmodels[3],fid,"diagnostic","hutter");
 	_printf_("\n   reading surface and bed slope computation model data:\n");
-	CreateFemModel(&femmodels[4],fid,"slope_compute");
+	CreateFemModel(&femmodels[4],fid,"slope_compute","");
 
 	_printf_("initialize inputs:\n");
Index: issm/trunk/src/c/parallel/diagnostic_core.cpp
===================================================================
--- issm/trunk/src/c/parallel/diagnostic_core.cpp	(revision 464)
+++ issm/trunk/src/c/parallel/diagnostic_core.cpp	(revision 465)
@@ -74,6 +74,6 @@
 			
 		if(debug)_printf_("%s\n","computing surface slope (x and y derivatives)...");
-		diagnostic_core_linear(&slopex,fem_sl,inputs,SurfaceSlopeComputeXAnalysisEnum());
-		diagnostic_core_linear(&slopey,fem_sl,inputs,SurfaceSlopeComputeYAnalysisEnum());
+		diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),SurfaceXAnalysisEnum());
+		diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),SurfaceYAnalysisEnum());
 
 		if (dim==3){
@@ -89,5 +89,5 @@
 
 		if(debug)_printf_("%s\n"," computing hutter velocities...");
-		diagnostic_core_linear(&ug,fem_dhu,inputs,DiagnosticHutterAnalysisEnum());
+		diagnostic_core_linear(&ug,fem_dhu,inputs,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
 
 		if(debug)_printf_("%s\n"," computing pressure according to MacAyeal...");
@@ -106,5 +106,5 @@
 		
 		if(debug)_printf_("%s\n"," computing horizontal velocities...");
-		diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticHorizAnalysisEnum());
+		diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
 
 		if(debug)_printf_("%s\n"," computing pressure according to MacAyeal...");
@@ -121,5 +121,5 @@
 		if(debug)_printf_("%s\n"," computing vertical velocities...");
 		inputs->Add("velocity",ug_horiz,numberofdofspernode_dh,numberofnodes);
-		diagnostic_core_linear(&ug_vert,fem_dv,inputs,DiagnosticVertAnalysisEnum());
+		diagnostic_core_linear(&ug_vert,fem_dv,inputs,DiagnosticAnalysisEnum(),VertAnalysisEnum());
 
 		if(debug)_printf_("%s\n"," combining horizontal and vertical velocities...");
@@ -138,6 +138,6 @@
 
 			if(debug)_printf_("%s\n","computing bed slope (x and y derivatives)...");
-			diagnostic_core_linear(&slopex,fem_sl,inputs,BedSlopeComputeXAnalysisEnum());
-			diagnostic_core_linear(&slopey,fem_sl,inputs,BedSlopeComputeYAnalysisEnum());
+			diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedXAnalysisEnum());
+			diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedYAnalysisEnum());
 			SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials);
 			SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials);
@@ -160,5 +160,5 @@
 			if(debug)_printf_("%s\n"," computing stokes velocities and pressure ...");
 			VecFree(&ug);
-			diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticStokesAnalysisEnum());
+			diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
 		
 			//decondition" pressure
Index: issm/trunk/src/c/parallel/diagnostic_core_linear.cpp
===================================================================
--- issm/trunk/src/c/parallel/diagnostic_core_linear.cpp	(revision 464)
+++ issm/trunk/src/c/parallel/diagnostic_core_linear.cpp	(revision 465)
@@ -11,5 +11,5 @@
 #include "../issm.h"
 
-void diagnostic_core_linear(Vec* pug,FemModel* fem,ParameterInputs* inputs,int analysis_type){
+void diagnostic_core_linear(Vec* pug,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 	/*parameters:*/
@@ -41,5 +41,5 @@
 	//*Generate system matrices
 	if (debug) _printf_("   Generating matrices\n");
-	SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); 
+	SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
 
 	/*!Reduce matrix from g to f size:*/
Index: issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp
===================================================================
--- issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp	(revision 464)
+++ issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp	(revision 465)
@@ -11,5 +11,5 @@
 #include "../issm.h"
 
-void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type){
+void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
 
 
@@ -84,9 +84,9 @@
 		if (debug) _printf_("   Generating matrices\n");
 		//*Generate system matrices
-		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); 
+		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
 
 		if (debug) _printf_("   Generating penalty matrices\n");
 		//*Generate penalty system matrices
-		PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type); 
+		PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); 
 
 		if (debug) _printf_("   reducing matrix from g to f set\n");
@@ -122,5 +122,5 @@
 		inputs->Add("velocity",ug,numberofdofspernode,numberofnodes);
 		
-		PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type); 
+		PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type,sub_analysis_type); 
 
 		//Figure out if convergence is reached.
@@ -167,5 +167,5 @@
 		kflag=1; pflag=0; //stiffness generation only
 	
-		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); 
+		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
 	
 		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);
Index: issm/trunk/src/c/parallel/objectivefunctionC.cpp
===================================================================
--- issm/trunk/src/c/parallel/objectivefunctionC.cpp	(revision 464)
+++ issm/trunk/src/c/parallel/objectivefunctionC.cpp	(revision 465)
@@ -39,4 +39,5 @@
 	double* p_g_copy=NULL;
 	int     analysis_type;
+	int     sub_analysis_type;
 	Vec     u_g=NULL;
 	double* u_g_double=NULL;
@@ -59,4 +60,5 @@
 	femmodel->parameters->FindParam((void*)&fit,"fit");
 	femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type");
+	femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
 	femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
 
@@ -75,5 +77,5 @@
 
 	//Run diagnostic with updated parameters.
-	diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,analysis_type);
+	diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,analysis_type,sub_analysis_type);
 	VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
 
@@ -81,5 +83,5 @@
 	inputs->Add("fit",fit[n]);
 	Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, 
-		u_g_double,u_g_obs, inputs,analysis_type);
+		u_g_double,u_g_obs, inputs,analysis_type,sub_analysis_type);
 
 
Index: issm/trunk/src/c/parallel/parallel.h
===================================================================
--- issm/trunk/src/c/parallel/parallel.h	(revision 464)
+++ issm/trunk/src/c/parallel/parallel.h	(revision 465)
@@ -12,6 +12,6 @@
 
 void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs);
-void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type);
-void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int  analysis_type);
+void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
+void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int  analysis_type,int sub_analysis_type);
 int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel);
 
@@ -29,9 +29,9 @@
 //int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams);
 void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename);
-int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type);
+int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type,int sub_analysis_type);
 void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets);
 void WriteLockFile(char* filename);
 
-void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type);
+void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type);
 //int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename);
 
Index: issm/trunk/src/c/parallel/thermal.cpp
===================================================================
--- issm/trunk/src/c/parallel/thermal.cpp	(revision 465)
+++ issm/trunk/src/c/parallel/thermal.cpp	(revision 465)
@@ -0,0 +1,102 @@
+/*!\file:  thermal.cpp
+ * \brief: thermal solution
+ */ 
+
+#include "../issm.h"
+#include "./parallel.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "thermal"
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+int main(int argc,char* *argv){
+
+	/*I/O: */
+	FILE* fid=NULL;
+	char* inputfilename=NULL;
+	char* outputfilename=NULL;
+	char* lockname=NULL;
+	int   numberofnodes;
+
+	/*Fem models : */
+	FemModel femmodels[2];
+
+	/*initial velocity and pressure: */
+	Vec u_g=NULL;
+	Vec p_g=NULL;
+
+	
+	/*solution vectors: */
+	Vec t_g=NULL;
+	Vec m_g=NULL;
+
+	ParameterInputs* inputs=NULL;
+	int waitonlock=0;
+	
+	MODULEBOOT();
+
+	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
+	throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!");
+	#endif
+
+	PetscInitialize(&argc,&argv,(char *)0,"");  
+
+	/*Size and rank: */
+	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
+	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
+
+	_printf_("recover , input file name and output file name:\n");
+	inputfilename=argv[2];
+	outputfilename=argv[3];
+	lockname=argv[4];
+
+	/*Open handle to data on disk: */
+	fid=fopen(inputfilename,"rb");
+	if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading")); 
+
+	_printf_("read and create finite element model:\n");
+	CreateFemModel(&femmodels[0],fid,"thermal","");
+	CreateFemModel(&femmodels[1],fid,"melting","");
+
+	_printf_("initialize inputs:\n");
+	femmodels[0].parameters->FindParam((void*)&u_g,"u_g");
+	femmodels[0].parameters->FindParam((void*)&p_g,"p_g");
+	femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+
+	inputs=new ParameterInputs;
+	//inputs->Add("velocity",u_g_initial,3,numberofnodes);
+
+	//erase velocities: 
+	//param=(Param*)femmodels[0].parameters->FindParamObject("u_g");
+	//femmodels[0].parameters->DeleteObject((Object*)param);
+
+	_printf_("call computational core:\n");
+	diagnostic_core(&u_g,&p_g,femmodels,inputs);
+
+	_printf_("write results to disk:\n");
+	OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename);
+
+	_printf_("write lock file:\n");
+	femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
+	if (waitonlock){
+		WriteLockFile(lockname);
+	}
+		
+	_printf_("closing MPI and Petsc\n");
+	MPI_Barrier(MPI_COMM_WORLD);
+
+	/*Close MPI libraries: */
+	PetscFinalize(); 
+
+
+	/*end module: */
+	MODULEEND();
+	
+	return 0; //unix success return;
+}
+	
Index: issm/trunk/src/c/parallel/thermalsteady.c
===================================================================
--- issm/trunk/src/c/parallel/thermalsteady.c	(revision 464)
+++ 	(revision )
@@ -1,111 +1,0 @@
-/*
- * cielothermalsteady.c:
- */
-
-#include "../include/cielo.h"
-#include "../modules.h"
-#include "./parallel.h"
-
-#undef __FUNCT__ 
-#define __FUNCT__ "cielothermalsteady"
-
-int main(int argc,char* *argv){
-	
-	/*I/O: */
-	FILE* fid=NULL;
-	char* inputfilename=NULL;
-	char* outputfilename=NULL;
-	char* lockname=NULL;
-	char* analysis_type_t="thermalsteady";
-	char* analysis_type_m="melting";
-
-	/*Intermediary: */
-	
-	#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
-	FemModel femmodel_t;
-	FemModel femmodel_m;	Vec* t_g=NULL;
-	#endif
-	ParameterInputs* inputs=NULL;
-
-	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
-		_printf_("%s%s\n",__FUNCT__," error message: parallel executable was compiled without support of parallel libraries!");
-		return 1;
-	#else
-
-		/*Initialize MPI environment: */
-		PetscInitialize(&argc,&argv,(char *)0,"");  
-
-		/*Size and rank: */
-		MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
-		MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
-
-		/*Some checks on size of cluster*/
-		if (num_procs<=1){
-			_printf_("\nSize of MPI COMM WORLD is 1, needs to be at least 2. Include more nodes\n"); 
-			PetscFinalize(); 
-			return 0;
-		}
-
-		/*Recover dbdir, input file name and output file name: */
-		dbdir=argv[1];
-		inputfilename=argv[2];
-		outputfilename=argv[3];
-		lockname=argv[4];
-
-		
-		
-		/*Read and create thermal finite element model: */
-		fid=fopen(inputfilename,"rb");
-		if(fid==NULL){
-			_printf_("%s%s\n",__FUNCT__," error message: could not open file ",inputfilename," for binary reading"); 
-			return 0;
-		}
-		if(!CreateFemModel(&femmodel_t,fid,analysis_type_t)){
-			_printf_("%s\n",__FUNCT__," error message: could not read melting finite element model!\n");
-			return 0;
-		}
-		fclose(fid);
-
-		/*Read and create melting finite element model: */
-		fid=fopen(inputfilename,"rb");
-		if(!CreateFemModel(&femmodel_m,fid,analysis_type_m)){
-			_printf_("%s\n",__FUNCT__," error message: could not read thermal finite element model!\n");
-			return 0;
-		}
-		fclose(fid);
-
-		/*Initialize inputs: */
-		inputs=NewParameterInputs();
-			
-		ParameterInputsAddFromMat(inputs,WorkspaceParamsGetVelocity(femmodel_t.workspaceparams),femmodel_t.workspaceparams->gsize,"velocity");
-		ParameterInputsAddFromMat(inputs,WorkspaceParamsGetPressure(femmodel_t.workspaceparams),femmodel_t.workspaceparams->gsize,"pressure");
-
-		ParameterInputsAddFromMat(inputs,&femmodel_t.batchparams->dt,1,"dt");
-		
-		/*Call computational core: */
-		if(!cielothermal_core(&t_g,inputs,&femmodel_t)){
-			_printf_("%s\n",__FUNCT__," error message: could not run computational core!\n");
-			return 0;
-		}
-		
-		/*Write results to disk: */
-		OutputThermal(t_g,femmodel_t.tpartition,outputfilename,analysis_type_t);
-
-		/*Write lock file if requested: */
-		if (femmodel_t.batchparams->waitonlock){
-			WriteLockFile(lockname);
-		}
-					
-		cleanup_and_return:
-		_printf_("done.\n");
-		
-		/*Synchronize everyone before exiting: */
-		MPI_Barrier(MPI_COMM_WORLD);
-
-		/*Close MPI libraries: */
-		PetscFinalize(); 
-		
-		return 0; //unix success return;
-
-	#endif
-}
Index: issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp
===================================================================
--- issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp	(revision 465)
+++ issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp	(revision 465)
@@ -0,0 +1,47 @@
+/*!\file:  DistributeNumDofs.cpp
+ * \brief: figure out the maximum number of dofs per grid.
+ */ 
+
+#include "../../include/macros.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../shared.h"
+
+#undef __FUNCT__  
+#define __FUNCT__ "DistributeNumDofs"
+
+int DistributeNumDofs(int *pnumdofs,int analysis_type,int sub_analysis_type){
+
+	int numdofs=2; //default numdofs
+	int i;
+
+	/*ok, according to analysis type: */
+	if (analysis_type==ControlAnalysisEnum()){
+		numdofs=2;
+	}
+	else if (analysis_type==DiagnosticAnalysisEnum()){
+		if (sub_analysis_type==HorizAnalysisEnum()){
+			numdofs=2;
+		}
+		else if (sub_analysis_type==VertAnalysisEnum()){
+			numdofs=1;
+		}
+		else if (sub_analysis_type==StokesAnalysisEnum()){
+			numdofs=4;
+		}
+		else if (sub_analysis_type==HutterAnalysisEnum()){
+			numdofs=2;
+		}
+	}
+	else if (analysis_type==SlopeComputeAnalysisEnum()){
+		numdofs=1;
+	}
+	else if (analysis_type==ThermalAnalysisEnum()){
+		numdofs=1;
+	}
+	else throw ErrorException(__FUNCT__," analysis type not implemented yet!");
+
+	/*Assign output pointers:*/
+	*pnumdofs=numdofs;;
+
+	return 1;
+}
Index: issm/trunk/src/c/shared/Dofs/dofs.h
===================================================================
--- issm/trunk/src/c/shared/Dofs/dofs.h	(revision 464)
+++ issm/trunk/src/c/shared/Dofs/dofs.h	(revision 465)
@@ -7,4 +7,5 @@
 
 double* dofsetgen(int numdofs,int* doflist,int dofspernode,int totaldofs);
+int DistributeNumDofs(int *pnumdofs,int analysis_type,int sub_analysis_type);
 
 
Index: issm/trunk/src/c/shared/Numerics/OptFunc.cpp
===================================================================
--- issm/trunk/src/c/shared/Numerics/OptFunc.cpp	(revision 464)
+++ issm/trunk/src/c/shared/Numerics/OptFunc.cpp	(revision 465)
@@ -23,5 +23,5 @@
 	double J;
 
-	mxArray*   inputs[8];
+	mxArray*   inputs[9];
 	mxArray*   psearch_scalar=NULL;
 	mxArray*   mxJ=NULL;
@@ -30,5 +30,4 @@
 	psearch_scalar=mxCreateDoubleScalar(scalar);
 	inputs[0]=psearch_scalar;
-	
 	inputs[1]=optargs->m;
 	inputs[2]=optargs->inputs;
@@ -38,6 +37,7 @@
 	inputs[6]=optargs->n;
 	inputs[7]=optargs->analysis_type;
+	inputs[8]=optargs->sub_analysis_type;
 
-	mexCallMATLAB( 1, &mxJ, 8, (mxArray**)inputs, optargs->function_name);
+	mexCallMATLAB( 1, &mxJ, 9, (mxArray**)inputs, optargs->function_name);
 
 	/*extract misfit from mxArray*/
Index: issm/trunk/src/c/shared/Numerics/numerics.h
===================================================================
--- issm/trunk/src/c/shared/Numerics/numerics.h	(revision 464)
+++ issm/trunk/src/c/shared/Numerics/numerics.h	(revision 465)
@@ -15,5 +15,4 @@
 void cross(double* result,double* vector1,double* vector2);
 double norm(double* vector);
-int DistributeNumDofs(int *pnumdofs,int analysis_type);
 
 #endif //ifndef _NUMERICS_H_
Index: issm/trunk/src/m/classes/@model/model.m
===================================================================
--- issm/trunk/src/m/classes/@model/model.m	(revision 464)
+++ issm/trunk/src/m/classes/@model/model.m	(revision 465)
@@ -255,4 +255,5 @@
 	md.solverstring='';
 	md.analysis_type='';
+	md.sub_analysis_type='';
 
 	%management of large models
Index: issm/trunk/src/m/classes/public/BuildQueueingScript.m
===================================================================
--- issm/trunk/src/m/classes/public/BuildQueueingScript.m	(revision 464)
+++ issm/trunk/src/m/classes/public/BuildQueueingScript.m	(revision 465)
@@ -1,7 +1,7 @@
-function BuildQueueingScript(md,solutiontype,executionpath,codepath)
+function BuildQueueingScript(md,executionpath,codepath)
 %BUILDQUEUEINGSCRIPT - 
 %
 %   Usage:
-%      BuildQueueingScript(md,solutiontype,executionpath,codepath)
+%      BuildQueueingScript(md,executionpath,codepath)
 
 disp('building queueing script');
@@ -11,7 +11,7 @@
 if exist(function_name,'file'),
 	%Call this function:
-	eval([function_name '(md,solutiontype,executionpath,codepath);']);
+	eval([function_name '(md,executionpath,codepath);']);
 else
 	%Call the generic BuildQueueingScript:
-	BuildQueueingScriptGeneric(md,solutiontype,executionpath,codepath);
+	BuildQueueingScriptGeneric(md,executionpath,codepath);
 end
Index: issm/trunk/src/m/classes/public/BuildQueueingScriptGeneric.m
===================================================================
--- issm/trunk/src/m/classes/public/BuildQueueingScriptGeneric.m	(revision 464)
+++ issm/trunk/src/m/classes/public/BuildQueueingScriptGeneric.m	(revision 465)
@@ -1,3 +1,3 @@
-function BuildQueueingScriptGeneric(md,solutiontype,executionpath,codepath)
+function BuildQueueingScriptGeneric(md,executionpath,codepath)
 %BUILDQUEUEINGSCRIPTGENERIC - ...
 %
@@ -17,15 +17,14 @@
 fprintf(fid,'mpirun -np %i ',md.np);
 
-if strcmpi(solutiontype,'diagnostic_horiz') |  strcmpi(solutiontype,'diagnostic'),
+if strcmpi(md.analysis_type,'diagnostic'),
 	fprintf(fid,'%s/diagnostic.exe',codepath);
-elseif strcmpi(solutiontype,'control'),
+elseif strcmpi(md.analysis_type,'control'),
 	fprintf(fid,'%s/control.exe',codepath);
-elseif strcmpi(solutiontype,'thermalsteady'),
-	fprintf(fid,'%s/thermalsteady.exe',codepath);
+elseif strcmpi(md.analysis_type,'thermal'),
+	fprintf(fid,'%s/thermal.exe',codepath);
 else
 	error('BuildQueueingScriptGeneric error message: unsupported solution type!');
 end
 
-
 fprintf(fid,' %s %s.bin %s.outbin %s.lock 2> %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name);
 fclose(fid);
Index: issm/trunk/src/m/classes/public/BuildQueueingScriptcosmos.m
===================================================================
--- issm/trunk/src/m/classes/public/BuildQueueingScriptcosmos.m	(revision 464)
+++ issm/trunk/src/m/classes/public/BuildQueueingScriptcosmos.m	(revision 465)
@@ -2,5 +2,5 @@
 %BUILDQUEUEINGSCRIPTCOSMOS - build queueing script for batch system when running parallel
 %
-%   solutiontype is 'diagnostic','prognostic','transient','thermalsteady','thermaltransient','parameters' or 'control'
+%   solutiontype is 'diagnostic','prognostic','transient','thermal','parameters' or 'control'
 %   and varargin is an optional package name ('cielo', 'ice' or 'macayeal')
 %
@@ -38,12 +38,12 @@
 	fprintf(fid,'mpirun.lsf /home/larour/bin/alloc_cleanup.exe\n');
 end
-if strcmpi(solutiontype,'diagnostic_horiz') |  strcmpi(solutiontype,'diagnostic'),
+if strcmpi(solutiontype,'diagnostic'),
 	fprintf(fid,'mpirun.lsf %s/diagnostic.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name);
 elseif strcmpi(solutiontype,'control') ,
 	fprintf(fid,'mpirun.lsf %s/control.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name);
-elseif strcmpi(solutiontype,'thermalsteady') ,
-	fprintf(fid,'mpirun.lsf %s/thermalsteady.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name);
+elseif strcmpi(solutiontype,'thermal') ,
+	fprintf(fid,'mpirun.lsf %s/thermal.exe %s %s.bin %s.outbin %s.lock',codepath,executionpath,md.name,md.name,md.name);
 else
-	error('BuildQueueingScriptcosmso error message: unsupported solution type!');
+	error('BuildQueueingScriptcosmos error message: unsupported solution type!');
 end
 fclose(fid);
Index: issm/trunk/src/m/classes/public/BuildQueueingScriptgemini.m
===================================================================
--- issm/trunk/src/m/classes/public/BuildQueueingScriptgemini.m	(revision 464)
+++ issm/trunk/src/m/classes/public/BuildQueueingScriptgemini.m	(revision 465)
@@ -22,10 +22,10 @@
 fprintf(fid,'rm -rf %s.outlog %s.errlog %s.lock\n',md.name,md.name,md.name);
 
-if strcmpi(solutiontype,'diagnostic_horiz') |  strcmpi(solutiontype,'diagnostic'),
+if strcmpi(solutiontype,'diagnostic') ,
 	fprintf(fid,'mpirun -np %i %s/diagnostic.exe %s %s.bin %s.outbin %s.lock',md.np,codepath,executionpath,md.name,md.name,md.name);
 elseif strcmpi(solutiontype,'control') ,
 	fprintf(fid,'mpirun -np %i %s/control.exe %s %s.bin %s.outbin %s.lock',md.np,codepath,executionpath,md.name,md.name,md.name);
-elseif strcmpi(solutiontype,'thermalsteady') ,
-	fprintf(fid,'mpirun -np %i %s/thermalsteady.exe %s %s.bin %s.outbin %s.lock',md.np,codepath,executionpath,md.name,md.name,md.name);
+elseif strcmpi(solutiontype,'thermal') ,
+	fprintf(fid,'mpirun -np %i %s/thermal.exe %s %s.bin %s.outbin %s.lock',md.np,codepath,executionpath,md.name,md.name,md.name);
 else
 	error('BuildQueueingScriptgemini error message: unsupported solution type!');
Index: issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 464)
+++ issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 465)
@@ -1,11 +1,11 @@
-function bool=ismodelselfconsistent(md,solutiontype,package),
+function bool=ismodelselfconsistent(md,package),
 %ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
 %
 %   Usage:
-%      bool=ismodelselfconsistent(md,solutiontype,package),
+%      bool=ismodelselfconsistent(md,package),
 
 %tolerance we use in litmus tests for the consistency of the model
 tolerance=10^-12;
-if (nargin~=3  )
+if (nargin~=2  )
 	help ismodelselfconsistent
 	error('ismodelselfconsistent error message: wrong usage');
@@ -115,5 +115,5 @@
 
 %SIZE NUMBEROFGRIDS
-fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface','deadgrids'};
+fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
 for i=1:length(fields),
 	if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
@@ -159,5 +159,5 @@
 
 %DIAGNOSTIC
-if strcmpi(solutiontype,'diagnostic')
+if strcmpi(md.analysis_type,'diagnostic')
 
 	%HUTTER ON ICESHELF WARNING
@@ -189,5 +189,5 @@
 
 %PROGNOSTIC
-if strcmp(solutiontype,'prognostic'),
+if strcmp(md.analysis_type,'prognostic'),
 	%old testing
 	%	if isempty(find(md.gridondirichlet_diag)),
@@ -208,23 +208,25 @@
 
 %THERMAL TRANSIENT
-if strcmp(solutiontype,'thermaltransient')
-
-	%INITIAL TEMPERATURE, MELTING AND ACCUMULATION
-	if isempty(md.temperature),
-		disp(['An  initial temperature is needed for a transient thermal computation'])
-		bool=0;return;
-	end
-	if isstruct(md.temperature) |  isstruct(md.melting) | isstruct(md.accumulation),
-		disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
-		bool=0;return;
+if strcmp(md.analysis_type,'thermal')
+	if strcmp(md.sub_analysis_type,'transient')
+
+		%INITIAL TEMPERATURE, MELTING AND ACCUMULATION
+		if isempty(md.temperature),
+			disp(['An  initial temperature is needed for a transient thermal computation'])
+			bool=0;return;
+		end
+		if isstruct(md.temperature) |  isstruct(md.melting) | isstruct(md.accumulation),
+			disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
+			bool=0;return;
+		end
 	end
 end
 
 %THERMAL STEADY AND THERMAL TRANSIENT
-if strcmpi(solutiontype,'thermalsteady') |  strcmp(solutiontype,'thermaltransient'),
+if strcmpi(md.analysis_type,'thermal'),
 
 	%EXTRUSION
 	if strcmp(md.type,'2d'),
-		disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!'])
+		disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
 		bool=0;return;
 	end
@@ -242,5 +244,5 @@
 
 %THERMAL TRANSIENT AND TRANSIENT
-if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'transient'),
+if strcmp(md.analysis_type,'thermal'),
 
 	%DT and NDT
@@ -261,5 +263,5 @@
 
 %PARAMETERS
-if strcmp(solutiontype,'parameters')
+if strcmp(md.analysis_type,'parameters')
 
 	%OUTPUT
@@ -289,5 +291,5 @@
 
 %CONTROL
-if strcmpi(solutiontype,'control'),
+if strcmpi(md.analysis_type,'control'),
 
 	%CONTROL TYPE
@@ -320,5 +322,5 @@
 
 %QMU
-if strcmpi(solutiontype,'qmu'),
+if strcmpi(md.analysis_type,'qmu'),
 	if ~strcmpi(md.cluster,'none'),
 		if md.waitonlock==0,
@@ -330,5 +332,5 @@
 
 %MESH
-if strcmpi(solutiontype,'mesh'),
+if strcmpi(md.analysis_type,'mesh'),
 	%this solution is a little special. It should come right after the md=model;  operation. So a lot less checks!
 
@@ -338,5 +340,5 @@
 
 %MESH2GRID
-if strcmpi(solutiontype,'mesh2grid'),
+if strcmpi(md.analysis_type,'mesh2grid'),
 	if ~strcmpi(md.cluster,'none'),
 		disp(['model is not correctly configured: mesh2grid not supported in parallel yet!']);
Index: issm/trunk/src/m/classes/public/isresultconsistent.m
===================================================================
--- issm/trunk/src/m/classes/public/isresultconsistent.m	(revision 464)
+++ issm/trunk/src/m/classes/public/isresultconsistent.m	(revision 465)
@@ -1,3 +1,3 @@
-function bool=isresultconsistent(md,solutiontype),
+function bool=isresultconsistent(md)
 %ISRESULTCONSISTENT: check that results are consistent
 %
@@ -6,15 +6,15 @@
 %
 %   Usage:
-%      bool=IsModelSelfConsistent(model,solutiontype)
+%      bool=IsModelSelfConsistent(model)
 
 %tolerance we use in litmus tests for the consistency of the model
 tolerance=10^-2;
 
-if (nargin~=2  )
+if (nargin~=1  )
 	IsResultConsistentUsage();
 	error(' ');
 end
 %Check velocities i(no penalties for 2d meshes)
-if (strcmp(solutiontype,'diagnostic') || strcmp(solutiontype,'diagnostic_horiz') ),
+if strcmp(md.analysis_type,'diagnostic'),
 
 	if strcmpi(md.type,'3d'),
@@ -37,9 +37,9 @@
 end
 
-if strcmp(solutiontype,'thermalsteady')
+if strcmp(md.analysis_type,'thermal')
 
 	%check melting 
 	if any(abs(md.melting(md.numberofgrids2d+1:end))>tolerance)
-		disp(['''thermalsteady'' result not consistent: increase penalty_melting (negative melting)']);
+		disp(['''thermal'' result not consistent: increase penalty_melting (negative melting)']);
 		bool=0; return; 
 	end
@@ -47,5 +47,5 @@
 end
 
-if strcmp(solutiontype,'thermaltransient')
+if strcmp(md.analysis_type,'thermaltransient')
 
 	for iter=1:length(md.thermaltransient_results)
@@ -61,5 +61,5 @@
 
 
-if (strcmp(solutiontype,'transient') & strcmpi(md.type,'3d')),
+if (strcmp(md.analysis_type,'transient') & strcmpi(md.type,'3d')),
 	for iter=1:length(md.transient_results)
 
@@ -99,4 +99,4 @@
 function IsResultConsistentUsage()
 disp(' ');
-disp('   IsResultConsistent usage: flag=IsResultConsistent(model,solutiontype)');
+disp('   IsResultConsistent usage: flag=IsResultConsistent(model)');
 disp('         where model is an instance of the @model class, and flag a boolean');
Index: issm/trunk/src/m/classes/public/loadresultsfromcluster.m
===================================================================
--- issm/trunk/src/m/classes/public/loadresultsfromcluster.m	(revision 464)
+++ issm/trunk/src/m/classes/public/loadresultsfromcluster.m	(revision 465)
@@ -1,7 +1,7 @@
-function md=loadresultsfromcluster(md,solutiontype)
+function md=loadresultsfromcluster(md)
 %LOADRESULTSFROMCLUSTER - load results of solution sequence from cluster
 %
 %   Usage:
-%      md=loadresultsfromcluster(md,solutiontype);
+%      md=loadresultsfromcluster(md);
 
 %Get cluster.rc location
@@ -31,5 +31,5 @@
 
 %If we are here, no errors in the solution sequence, call loadresultsfromdisk.
-md=loadresultsfromdisk(md,[md.name '.outbin'],solutiontype);
+md=loadresultsfromdisk(md,[md.name '.outbin'],md.analysis_type);
 
 %erase the log and output files
Index: issm/trunk/src/m/classes/public/loadresultsfromdisk.m
===================================================================
--- issm/trunk/src/m/classes/public/loadresultsfromdisk.m	(revision 464)
+++ issm/trunk/src/m/classes/public/loadresultsfromdisk.m	(revision 465)
@@ -1,3 +1,3 @@
-function md=loadresultsfromdisk(md,filename,solutiontype)
+function md=loadresultsfromdisk(md,filename)
 %LOADRESULTSFROMDISK - load results of solution sequence from disk file "filename"            
 %
@@ -103,5 +103,5 @@
 disp(sprintf('%s\n','checking result consistency'));
 
-if ~isresultconsistent(md,solutiontype),
+if ~isresultconsistent(md),
 	disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
 end
Index: issm/trunk/src/m/classes/public/marshall.m
===================================================================
--- issm/trunk/src/m/classes/public/marshall.m	(revision 464)
+++ issm/trunk/src/m/classes/public/marshall.m	(revision 465)
@@ -1,11 +1,10 @@
-function marshall(md,solutiontype,package)
+function marshall(md)
 %MARSHALL - output Cielo compatible binary file from @model md, for certain solution type.
 %
 %   The routine creates a Cielo compatible binary file from @model md
-%   solutiontype can be 'diagnostic','prognostic' or 'control'. 
 %   This binary file will be used for parallel runs in cielo.
 %
 %   Usage:
-%      marshall(md,solutiontype,package)
+%      marshall(md)
 
 %some checks on list of arguments
@@ -23,9 +22,6 @@
 end
 
-if strcmp(solutiontype,'diagnostic'),
-	WriteData(fid,'diagnostic_horiz','String','analysis_type');
-else
-	WriteData(fid,solutiontype,'String','analysis_type');
-end
+WriteData(fid,md.analysis_type','String','analysis_type');
+WriteData(fid,md.sub_analysis_type','String','sub_analysis_type');
 WriteData(fid,md.type,'String','type');
 WriteData(fid,md.numberofgrids,'Integer','numberofgrids');
Index: issm/trunk/src/m/classes/public/process_solve_options.m
===================================================================
--- issm/trunk/src/m/classes/public/process_solve_options.m	(revision 465)
+++ issm/trunk/src/m/classes/public/process_solve_options.m	(revision 465)
@@ -0,0 +1,79 @@
+function outoptions=process_solve_options(options)
+%DEFAULT_SOLVE_OPTIONS - set up default options for solve phase
+%
+%   Usage:
+%      options=process_solve_options(options)
+%
+%   See also: SOLVE,RECOVER_SOLVE_OPTIONS
+
+
+
+
+%analysis_type: check on this option, error out otherwise
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'analysis_type'),
+		analysis_type=options{i,2};
+		found=1;
+	end
+end
+if ~found,
+	error('recover_solve_options error message: no ''analysis_type'' was provided');
+end
+
+%package: is there one? check to ''cielo''
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'package'),
+		package=options{i,2};
+		found=1;
+	end
+end
+if ~found,
+	disp('recover_solve_options info message: no ''package'' was provided, defaulting to ''cielo''');
+	options(end+1,:)={'package' 'cielo'};
+	package='cielo';
+end
+
+if ~ischar(package), 
+	error(['process_solve_options error message: package ' package ' not supported yet']);
+end
+
+%sub_analysis_type: check on it, not mandatory
+found=0;
+for i=1:size(options,1),
+	if strcmpi(options{i,1},'sub_analysis_type'),
+		sub_analysis_type=options{i,2};
+		found=1;
+	end
+end
+if ~found,
+	sub_analysis_type='';
+end
+
+%check solution type is supported
+if ~(strcmpi(analysis_type,'control') |  ...
+	strcmpi(analysis_type,'diagnostic') |  ...
+	strcmpi(analysis_type,'prognostic') |  ...
+	strcmpi(analysis_type,'thermal') |  ...
+	strcmpi(analysis_type,'mesh') |  ...
+	strcmpi(analysis_type,'mesh2grid') |  ...
+	strcmpi(analysis_type,'qmu') |  ...
+	strcmpi(analysis_type,'transient') ),
+	error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']);
+end
+
+if ~(strcmpi(sub_analysis_type,'') |  ...
+	strcmpi(sub_analysis_type,'steady_state') |  ...
+	strcmpi(sub_analysis_type,'horiz') |  ...
+	strcmpi(sub_analysis_type,'vert') |  ...
+	strcmpi(sub_analysis_type,'transient') ),
+	error(['process_solve_options error message: sub_analysis_type ' sub_analysis_type ' not supported yet!']);
+end
+
+
+
+%setup final options structure
+outoptions.analysis_type=analysis_type;
+outoptions.package=package;
+outoptions.sub_analysis_type=sub_analysis_type;
Index: issm/trunk/src/m/classes/public/recover_solve_options.m
===================================================================
--- issm/trunk/src/m/classes/public/recover_solve_options.m	(revision 465)
+++ issm/trunk/src/m/classes/public/recover_solve_options.m	(revision 465)
@@ -0,0 +1,22 @@
+function options=recover_solve_options(md,varargin)
+%RECOVER_SOLVE_OPTIONS - recover solution options
+%
+%   Usage:
+%      options=recover_solve_options(md,varargin);
+%
+%   See also: SOLVE
+
+%initialize options.
+options=cell(0,2);
+
+%make sure length(varargin) is even, ie options come in pairs.
+if mod(length(varargin),2),
+	error('recover_solve_options error message: an even number of options is necessary');
+end
+
+%go through varargin, extract options 
+for i=1:length(varargin)/2,
+
+	options(end+1,:)={varargin{2*i-1} varargin{2*i}};
+
+end
Index: issm/trunk/src/m/classes/public/solve.m
===================================================================
--- issm/trunk/src/m/classes/public/solve.m	(revision 464)
+++ issm/trunk/src/m/classes/public/solve.m	(revision 465)
@@ -1,64 +1,47 @@
-function md=solve(md,solutiontype,varargin)
+function md=solve(md,varargin)
 %SOLVE - apply solution sequence for this model
 %
-%   solutiontype is 'diagnostic','prognostic','transient','thermalsteady','thermaltransient','parameters' or 'control'
-%   and varargin is an optional package name ('cielo', 'ice' or 'macayeal')
-%
 %   Usage:
-%       md=solve(md,solutiontype,varargin)
+%       md=solve(md,varargin)
+%       where varargin is a lit of paired arguments. 
+%       arguments can be: 'analysis_type': 'diagnostic','thermal','prognostic','transient'
+%       arguments can be: 'sub_analysis_type': 'transient',''
+%       arguments can be: 'package': 'macayeal','ice','cielo'
 %
 %   Examples:
-%      md=solve(md,'diagnostic','macayeal');
-%      md=solve(md,'control','cielo');
+%      md=solve(md,'analysis_type','diagnostic','package','cielo');
+%      md=solve(md,'analysis_type','control','package','ice');
+%      md=solve(md,'analysis_type','thermal','sub_analysis_type','transient','package','ice');
+%      md=solve(md,'analysis_type','thermal','sub_analysis_type','','package','cielo');
 
 %some checks on list of arguments
-
 global ISSM_DIR
 
-solutions={'mesh2grid','mesh','thermaltransient','thermalsteady','qmu','diagnostic','diagnostic_horiz','prognostic','transient','parameters','control'};
-found=0;
-for i=1:length(solutions),
-	if strcmpi(solutiontype,solutions{i}),
-		found=1;
-		break;
-	end
-end
+%recover options
+options=recover_solve_options(md,varargin{:});
 
-if found==0,
-	error(['solve error message: solution type ' solutiontype ' not supported yet!']);
-end
+%add default options
+options=process_solve_options(options);
 
-%we are good with this solutiontype, put in the analysis_type field of md: 
-md.analysis_type=solutiontype;
-
-%Recover type of package being used: 
-if nargin==2,
-	package='Ice';
-else
-	package=varargin{1};
-end
-
-if ~ischar(package), 
-	error('Package specified in varargin can only be ''ice'', or ''cielo''');
-end
-
-if ~(strcmpi(package,'ice') || strcmpi(package,'cielo') || strcmpi(package,'macayeal'))
-	error('Package specified in varargin can only be ''ice'', ''macayeal'', or ''cielo''');
-end
+%recover some fields
+md.analysis_type=options.analysis_type;
+md.sub_analysis_type=options.sub_analysis_type;
+package=options.package;
 
 %Use package to set solution namespace
 usenamespace(package);
 
-if ~ismodelselfconsistent(md,solutiontype,package),
-		error(' '); %previous error messages should explain what is going on.
+if ~ismodelselfconsistent(md,package),
+	error(' '); %previous error messages should explain what is going on.
 end
 
 displaystring(md.debug,'\n%s\n','launching solution sequence');
 
+
 %If running in parallel, we have a different way of launching the solution
-%sequences.
-if ~strcmpi(solutiontype,'qmu'),
+%sequences. qmu is the only solution sequence that cannot run in parallel
+if ~strcmpi(md.analysis_type,'qmu'),
 	if ~strcmpi(md.cluster,'none'),
-		md=solveparallel(md,solutiontype,package);
+		md=solveparallel(md)
 		return;
 	end
@@ -66,30 +49,27 @@
 
 %Lauch correct solution sequence
-if strcmpi(solutiontype,'diagnostic'),
+if strcmpi(md.analysis_type,'diagnostic'),
 	md=diagnostic(md);
 
-elseif strcmpi(solutiontype,'mesh'),
+elseif strcmpi(md.analysis_type,'mesh'),
 	md=mesh(md);
 
-elseif strcmpi(solutiontype,'transient'),
+elseif strcmpi(md.analysis_type,'transient'),
 	md=transient(md);
 
-elseif strcmpi(solutiontype,'qmu'),
+elseif strcmpi(md.analysis_type,'qmu'),
 	md=qmu(md,package);
 
-elseif strcmpi(solutiontype,'parameters'),
-	md=parameters(md);;
-
-elseif strcmpi(solutiontype,'mesh2grid'),
+elseif strcmpi(md.analysis_type,'mesh2grid'),
 	md=mesh2grid(md);;
 
-elseif strcmpi(solutiontype,'prognostic'),
+elseif strcmpi(md.analysis_type,'prognostic'),
 	md=prognostic(md);;
 
-elseif strcmpi(solutiontype,'control'),
+elseif strcmpi(md.analysis_type,'control'),
 	md=control(md);
 
-elseif strcmpi(solutiontype,'thermalsteady') | strcmpi(solutiontype,'thermaltransient'),
-	md=thermal(md,solutiontype);
+elseif strcmpi(md.analysis_type,'thermal'),
+	md=thermal(md);
 else
 	error('solution type not supported for this model configuration.');
@@ -98,5 +78,5 @@
 %Check result is consistent
 displaystring(md.debug,'%s\n','checking result consistency');
-if ~isresultconsistent(md,solutiontype),
+if ~isresultconsistent(md),
 	disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
 end
@@ -110,5 +90,5 @@
 function solveusage();
 disp(' ');
-disp('   Solve usage: md=solve(md,solutiontype,package)');
-disp('         solutiontype can be ''diagnostic'',''transient'', ''thermalsteady'',''thermaltransient'',''parameters'',''mesh2grid''  or ''control''');
+disp('   Solve usage: md=solve(md,md.analysis_type,package)');
+disp('         md.analysis_type can be ''diagnostic'',''transient'', ''thermal'',''thermaltransient'',''parameters'',''mesh2grid''  or ''control''');
 disp('         package is either ''cielo'' or ''ice''');
Index: issm/trunk/src/m/classes/public/solveparallel.m
===================================================================
--- issm/trunk/src/m/classes/public/solveparallel.m	(revision 464)
+++ issm/trunk/src/m/classes/public/solveparallel.m	(revision 465)
@@ -1,22 +1,7 @@
-function md=solveparallel(md,solutiontype,varargin)
+function md=solveparallel(md)
 %SOLVEPARALLEL - solution sequence using a cluster in parallel mode
 %
 %   Usage:
-%      md=solveparallel(md,solutiontype,varargin)
-
-%Recover type of package being used: 
-if nargin==2,
-	package='Ice';
-else
-	package=varargin{1};
-end
-
-if ~ischar(package), 
-	error('Package specified in varargin can only be ''ice'', or ''cielo''');
-end
-
-if ~(strcmpi(package,'ice') || strcmpi(package,'cielo') || strcmpi(package,'macayeal'))
-	error('Package specified in varargin can only be ''ice'', ''macayeal'', or ''cielo''');
-end
+%      md=solveparallel(md);
 
 %Get cluster.rc location
@@ -27,8 +12,8 @@
 
 %Marshall model data into a binary file.
-marshall(md,solutiontype,package);
+marshall(md);
 
 %Now, we need to build the queuing script, used by the cluster to launch the job.
-BuildQueueingScript(md,solutiontype,executionpath,codepath);
+BuildQueueingScript(md,executionpath,codepath);
 
 %Now, launch the queueing script
@@ -40,5 +25,5 @@
 	waitonlock([executionpath '/' md.name '.lock']);
 	%load results
-	md=loadresultsfromcluster(md,solutiontype);
+	md=loadresultsfromcluster(md);
 else
 	return;
Index: issm/trunk/src/m/solutions/cielo/GradJCompute.m
===================================================================
--- issm/trunk/src/m/solutions/cielo/GradJCompute.m	(revision 464)
+++ issm/trunk/src/m/solutions/cielo/GradJCompute.m	(revision 465)
@@ -1,3 +1,3 @@
-function [u_g grad_g]=GradJCompute(m,inputs, u_g_obs,analysis_type);
+function [u_g grad_g]=GradJCompute(m,inputs, u_g_obs,analysis_type,sub_analysis_type);
 
 %Recover solution for this stiffness and right hand side: 
@@ -5,5 +5,5 @@
 	disp('         computing velocities...')
 end
-[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type);
+[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
 
 %Buid Du, difference between observed velocity and model velocity.
@@ -11,5 +11,5 @@
 	disp('         computing Du...')
 end
-[Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g,u_g_obs,inputs);
+[Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
 
 %Reduce adjoint load from g-set to f-set
@@ -26,3 +26,3 @@
 
 %Compute gradJ 
-grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, lambda_g,inputs);
+grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, lambda_g,inputs,analysis_type,sub_analysis_type);
Index: issm/trunk/src/m/solutions/cielo/control.m
===================================================================
--- issm/trunk/src/m/solutions/cielo/control.m	(revision 464)
+++ issm/trunk/src/m/solutions/cielo/control.m	(revision 465)
@@ -6,4 +6,5 @@
 	%Build all models requested for control simulation
 	md.analysis_type='control';
+	md.sub_analysis_type='';
 	m=CreateFemModel(md); 
 
@@ -40,5 +41,5 @@
 
 		disp('      computing gradJ...');
-		[u_g c(n).grad_g]=GradJCompute(m,inputs,u_g_obs,md.analysis_type);
+		[u_g c(n).grad_g]=GradJCompute(m,inputs,u_g_obs,md.analysis_type,md.sub_analysis_type);
 		inputs=add(inputs,'velocity',u_g,'doublevec',2,m.parameters.numberofnodes);
 		disp('      done.');
@@ -58,5 +59,5 @@
 		
 		disp('      optimizing along gradient direction...'); 
-		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,inputs,p_g,u_g_obs,c(n).grad_g,n,md.analysis_type);
+		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,inputs,p_g,u_g_obs,c(n).grad_g,n,md.analysis_type,md.sub_analysis_type);
 		disp('      done.');
 
@@ -76,5 +77,5 @@
 		%some temporary saving 
 		if(mod(n,5)==0),
-			solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type);
+			solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type,md.sub_analysis_type);
 			save temporary_control_results solution
 		end
@@ -85,5 +86,5 @@
 	%Create final solution
 	disp('      preparing final velocity solution...');
-	solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type);
+	solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type,md.sub_analysis_type);
 	disp('      done.');
 	
Index: issm/trunk/src/m/solutions/cielo/controlfinalsol.m
===================================================================
--- issm/trunk/src/m/solutions/cielo/controlfinalsol.m	(revision 464)
+++ issm/trunk/src/m/solutions/cielo/controlfinalsol.m	(revision 465)
@@ -1,7 +1,7 @@
-function solution=controlfinalsol(c,m,p_g,inputs,analysis_type);
+function solution=controlfinalsol(c,m,p_g,inputs,analysis_type,sub_analysis_type);
 
 %From parameters, build inputs for icediagnostic_core, using the final parameters
 inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes);
-u_g=diagnostic_core_nonlinear(m,inputs,analysis_type);
+u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
 
 %Build partitioning vectors to recover solution
Index: issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 464)
+++ issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 465)
@@ -10,17 +10,17 @@
 	%Build all models requested for diagnostic simulation
 	displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
-	md.analysis_type='diagnostic_horiz'; m_dh=CreateFemModel(md);
+	md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; m_dh=CreateFemModel(md);
 	
 	displaystring(md.debug,'\n%s',['reading diagnostic vert model data']);
-	md.analysis_type='diagnostic_vert'; m_dv=CreateFemModel(md);
+	md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; m_dv=CreateFemModel(md);
 	
 	displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']);
-	md.analysis_type='diagnostic_stokes'; m_ds=CreateFemModel(md);
+	md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; m_ds=CreateFemModel(md);
 
 	displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']);
-	md.analysis_type='diagnostic_hutter'; m_dhu=CreateFemModel(md);
+	md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; m_dhu=CreateFemModel(md);
 	
 	displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']);
-	md.analysis_type='slope_compute'; m_sl=CreateFemModel(md);
+	md.analysis_type='slope_compute'; md.sub_analysis_type=''; m_sl=CreateFemModel(md);
 	
 	% figure out number of dof: just for information purposes.
Index: issm/trunk/src/m/solutions/cielo/diagnostic_core.m
===================================================================
--- issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 464)
+++ issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 465)
@@ -16,6 +16,6 @@
 
 	displaystring(debug,'\n%s',['computing surface slope (x and y derivatives)...']);
-	slopex=diagnostic_core_linear(m_sl,inputs,'surface_slope_compute_x');
-	slopey=diagnostic_core_linear(m_sl,inputs,'surface_slope_compute_y');
+	slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','surfacex');
+	slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','surfacey');
 
 	if dim==3,
@@ -30,5 +30,5 @@
 
 	displaystring(debug,'\n%s',['computing hutter velocities...']);
-	u_g=diagnostic_core_linear(m_dhu,inputs,'diagnostic_hutter');
+	u_g=diagnostic_core_linear(m_dhu,inputs,'diagnostic','hutter');
 
 	displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']);
@@ -46,5 +46,5 @@
 
 	displaystring(debug,'\n%s',['computing horizontal velocities...']);
-	u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic_horiz');
+	u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic','horiz');
 
 	displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']);
@@ -59,5 +59,5 @@
 	displaystring(debug,'\n%s',['computing vertical velocities...']);
 	inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
-	u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic_vert');
+	u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic','vert');
 
 	displaystring(debug,'\n%s',['combining horizontal and vertical velocities...']);
@@ -75,6 +75,6 @@
 
 		displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
-		slopex=diagnostic_core_linear(m_sl,inputs,'bed_slope_compute_x');
-		slopey=diagnostic_core_linear(m_sl,inputs,'bed_slope_compute_y');
+		slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedx');
+		slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedy');
 		slopex=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex);
 		slopey=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey);
@@ -95,5 +95,5 @@
 
 		displaystring(debug,'\n%s',['computing stokes velocities and pressure ...']);
-		u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic_stokes');
+		u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic','stokes');
 	
 		%"decondition" pressure
Index: issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m
===================================================================
--- issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m	(revision 464)
+++ issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m	(revision 465)
@@ -1,3 +1,3 @@
-function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type)
+function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type)
 %INPUT function [ru_g varargout]=cielodiagnostic_core_nonlinear(m,inputs)
 	
@@ -30,8 +30,8 @@
 		
 		%system matrices 
-		[K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type);
+		[K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 		
 		%penalties
-		[K_gg , p_g]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type);
+		[K_gg , p_g]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 
 
@@ -60,5 +60,5 @@
 	
 		%penalty constraints
-		[loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs);
+		[loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 
 	%   Figure out if convergence is reached.
@@ -107,5 +107,5 @@
 		inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
 		m.parameters.kflag=1; m.parameters.pflag=0; 
-		[K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type);
+		[K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
 		varargout(1)={K_ff};
Index: issm/trunk/src/m/solutions/cielo/objectivefunctionC.m
===================================================================
--- issm/trunk/src/m/solutions/cielo/objectivefunctionC.m	(revision 464)
+++ issm/trunk/src/m/solutions/cielo/objectivefunctionC.m	(revision 465)
@@ -1,3 +1,3 @@
-function J =objectivefunctionC(search_scalar,m,inputs,p_g,u_g_obs,grad_g,n,analysis_type);
+function J =objectivefunctionC(search_scalar,m,inputs,p_g,u_g_obs,grad_g,n,analysis_type,sub_analysis_type);
         
 %recover some parameters
@@ -13,6 +13,6 @@
 
 %Run diagnostic with updated parameters.
-u_g=diagnostic_core_nonlinear(m,inputs,analysis_type);
+u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
 
 %Compute misfit for this velocity field. 
-J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, u_g_obs,inputs);
+J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, u_g_obs,inputs,analysis_type,sub_analysis_type);
Index: issm/trunk/src/m/solutions/cielo/thermal.m
===================================================================
--- issm/trunk/src/m/solutions/cielo/thermal.m	(revision 464)
+++ issm/trunk/src/m/solutions/cielo/thermal.m	(revision 465)
@@ -7,33 +7,32 @@
 	%timing
 	t1=clock;
-	
-	%md.analysis_type is set to thermalsteady or thermaltransient. save it.
-	solutiontype=md.analysis_type;
 
-	%Build all models requested for thermal and melting  simulation
-	md.analysis_type=solutiontype; m_t=CreateFemModel(md); 
-	md.analysis_type='melting' m_m=CreateFEMModel(md);
+	%Build all models requested for diagnostic simulation
+	displaystring(md.debug,'%s',['reading thermal model data']);
+	md.analysis_type='thermal'; m_t=CreateFemModel(md);
+
+	displaystring(md.debug,'%s',['reading melting model data']);
+	md.analysis_type='melting'; m_m=CreateFemModel(md);
 
 	% figure out number of dof: just for information purposes.
-	md.dof=m_t.nodesets.fsize; %biggest dof number
+	md.dof=m_t.nodesets.fsize;
 
-
-	if strcmpi(solutiontype,'thermalsteady'),
-
-		%initialize velocity and pressure 
-		u_g=m_t.parameters.u_g;
-		p_g=m_t.parameters.p_g;
-
-		inputs=struct('pressure',p_g,'velocity',u_g);
-
-		disp(sprintf('\n%s\n','computing temperature...'));
-		[t_g,m_t]=thermal_core(m_t,inputs);
-
-		disp(sprintf('\n%s\n','computing melting...'));
-		inputs=struct('pressure',p_g,'temperature',t_g);
+	%initialize inputs
+	displaystring(debug,'\n%s',['setup inputs...']);
+	inputs=inputlist;
+	inputs=add(inputs,'velocity',m_t.parameters.u_g,'doublevec',3,m_t.parameters.numberofnodes);
+	inputs=add(inputs,'pressure',m_t.parameters.p_g,'doublevec',1,m_t.parameters.numberofnodes);
+	inputs=add(inputs,'dt',m_t.parameters.dt,'double');
+	
+	if strcmpi(md.sub_analysis_type,'steady'),
+	
+		displaystring(debug,'\n%s',['computing temperatures...']);
+		[t_g loads_t melting_offset]=thermal_core(m_t,inputs,'thermal','steady');
 		
-		melting_g=melting_core(m_m,inputs);
-
-		%load results onto model
+		displaystring(debug,'\n%s',['computing melting...']);
+		inputs=add(inputs,'melting_offset',melting_offset,'double');
+		melting_g=melting_core(m_m,inputs,'melting','steady');
+		
+		displaystring(debug,'\n%s',['load results...']);
 		md.temperature=t_g;
 		md.melting=melting_g*md.yts; %from m/s to m/a
@@ -41,9 +40,9 @@
 	else
 
-		%initialize velocity, pressure  and temperature
-		u_g=m.parameters.u_g;
-		p_g=m.parameters.p_g;
-		t_g=m.parameters.t_g;
-		melting_g=m.parameters.melting_g;
+		%initialize velocity, pressure, temperature and melting
+		u_g=m_t.parameters.u_g;
+		p_g=m_t.parameters.p_g;
+		t_g=m_t.parameters.t_g;
+		melting_g=m_t.parameters.melting_g;
 
 		nsteps=m_t.parameters.ndt/m_t.parameters.dt;
@@ -55,9 +54,9 @@
 		for n=1:nsteps, 
 
-			disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt));
-			disp('   computing temperature...');
-		
-			inputs=struct('pressure',p_g,'temperature',soln(n).t_g,'velocity',u_g,'dt',m_t.parameters.dt);
-			[soln(n+1).t_g,m_t]=thermal_core(m_t,inputs);
+			displaystring(debug,'\n%s%i/%i\n','time step: ',n,nsteps);
+			
+			displaystring(debug,'\n%s',['    computing temperatures...']);
+			inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
+			[soln(n+1).t_g loads_t melting_offset]=thermal_core(m_t,inputs);
 			
 			disp('   computing melting...');
@@ -80,3 +79,6 @@
 		md.melting=solution_melting;
 	end
-	
+
+	%stop timing
+	t2=clock;
+	displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);	
Index: issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp
===================================================================
--- issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp	(revision 464)
+++ issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp	(revision 465)
@@ -51,4 +51,5 @@
 	optargs.n=STEP;
 	optargs.analysis_type=ANALYSIS;
+	optargs.sub_analysis_type=SUBANALYSIS;
 
 	optpars.xmin=xmin;
@@ -73,5 +74,5 @@
 {
 	_printf_("\n");
-	_printf_("   usage: [search_scalar J] = %s(function_name,xmin,xmax,options,m,inputs,p_g,u_g_obs,grad_g,step,analysis_type)\n",__FUNCT__);
+	_printf_("   usage: [search_scalar J] = %s(function_name,xmin,xmax,options,m,inputs,p_g,u_g_obs,grad_g,step,analysis_type,sub_analysis_type)\n",__FUNCT__);
 	_printf_("\n");
 }
Index: issm/trunk/src/mex/ControlOptimization/ControlOptimization.h
===================================================================
--- issm/trunk/src/mex/ControlOptimization/ControlOptimization.h	(revision 464)
+++ issm/trunk/src/mex/ControlOptimization/ControlOptimization.h	(revision 465)
@@ -28,4 +28,5 @@
 #define STEP (mxArray*)prhs[9]
 #define ANALYSIS (mxArray*)prhs[10]
+#define SUBANALYSIS (mxArray*)prhs[11]
 
 /* serial output macros: */
@@ -37,5 +38,5 @@
 #define NLHS  2
 #undef NRHS
-#define NRHS  11
+#define NRHS  12
 
 
Index: issm/trunk/src/mex/Du/Du.cpp
===================================================================
--- issm/trunk/src/mex/Du/Du.cpp	(revision 464)
+++ issm/trunk/src/mex/Du/Du.cpp	(revision 465)
@@ -15,8 +15,11 @@
 	DataSet* loads=NULL;
 	DataSet* materials=NULL;
-	int      analysis_type;
 	double*  u_g=NULL;
 	double*  u_g_obs=NULL;
 	ParameterInputs* inputs=NULL;
+	char*             analysis_type_string=NULL;
+	int               analysis_type;
+	char*             sub_analysis_type_string=NULL;
+	int               sub_analysis_type;
 
 	/* output datasets: */
@@ -35,7 +38,10 @@
 	FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
 	FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
-	FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL);
 	FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec");
 	FetchData((void**)&u_g_obs,NULL,NULL,UGOBS,"Vector","Vec");
+	FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
+	analysis_type=AnalysisTypeAsEnum(analysis_type_string);
+	FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
+	sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
 
 	/*Fetch inputs: */
@@ -44,5 +50,5 @@
 
 	/*!Call core code: */
-	Dux(&du_g, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type);
+	Dux(&du_g, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
 
 	/*write output : */
@@ -58,4 +64,6 @@
 	VecFree(&du_g);
 	delete inputs;
+	xfree((void**)&analysis_type_string);
+	xfree((void**)&sub_analysis_type_string);
 
 	/*end module: */
Index: issm/trunk/src/mex/Du/Du.h
===================================================================
--- issm/trunk/src/mex/Du/Du.h	(revision 464)
+++ issm/trunk/src/mex/Du/Du.h	(revision 465)
@@ -25,4 +25,6 @@
 #define UGOBS (mxArray*)prhs[6]
 #define INPUTS (mxArray*)prhs[7]
+#define ANALYSIS (mxArray*)prhs[8]
+#define SUBANALYSIS (mxArray*)prhs[9]
 
 /* serial output macros: */
@@ -33,5 +35,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  8
+#define NRHS  10
 
 
Index: issm/trunk/src/mex/Gradj/Gradj.cpp
===================================================================
--- issm/trunk/src/mex/Gradj/Gradj.cpp	(revision 464)
+++ issm/trunk/src/mex/Gradj/Gradj.cpp	(revision 465)
@@ -15,9 +15,12 @@
 	DataSet* loads=NULL;
 	DataSet* materials=NULL;
-	int      analysis_type;
 	char*    control_type=NULL;
 	double*  u_g=NULL;
 	double*  lambda_g=NULL;
 	ParameterInputs* inputs=NULL;
+	char*             analysis_type_string=NULL;
+	int               analysis_type;
+	char*             sub_analysis_type_string=NULL;
+	int               sub_analysis_type;
 
 	/* output datasets: */
@@ -36,8 +39,11 @@
 	FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
 	FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
-	FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL);
 	FetchData((void**)&control_type,NULL,NULL,mxGetField(PARAMETERS,0,"control_type"),"String",NULL);
 	FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec");
 	FetchData((void**)&lambda_g,NULL,NULL,LAMBDAG,"Vector","Vec");
+	FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
+	analysis_type=AnalysisTypeAsEnum(analysis_type_string);
+	FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
+	sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
 
 	/*Fetch inputs: */
@@ -46,5 +52,5 @@
 
 	/*!Call core code: */
-	Gradjx(&grad_g, elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,control_type);
+	Gradjx(&grad_g, elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type);
 
 	/*write output : */
@@ -61,4 +67,6 @@
 	delete inputs;
 	VecFree(&grad_g);
+	xfree((void**)&analysis_type_string);
+	xfree((void**)&sub_analysis_type_string);
 
 	/*end module: */
Index: issm/trunk/src/mex/Gradj/Gradj.h
===================================================================
--- issm/trunk/src/mex/Gradj/Gradj.h	(revision 464)
+++ issm/trunk/src/mex/Gradj/Gradj.h	(revision 465)
@@ -25,4 +25,6 @@
 #define LAMBDAG (mxArray*)prhs[6]
 #define INPUTS (mxArray*)prhs[7]
+#define ANALYSIS (mxArray*)prhs[8]
+#define SUBANALYSIS (mxArray*)prhs[9]
 
 /* serial output macros: */
@@ -33,5 +35,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  8
+#define NRHS  10
 
 
Index: issm/trunk/src/mex/Misfit/Misfit.cpp
===================================================================
--- issm/trunk/src/mex/Misfit/Misfit.cpp	(revision 464)
+++ issm/trunk/src/mex/Misfit/Misfit.cpp	(revision 465)
@@ -15,8 +15,11 @@
 	DataSet* loads=NULL;
 	DataSet* materials=NULL;
-	int      analysis_type;
 	double*  u_g=NULL;
 	double*  u_g_obs=NULL;
 	ParameterInputs* inputs=NULL;
+	char*             analysis_type_string=NULL;
+	int               analysis_type;
+	char*             sub_analysis_type_string=NULL;
+	int               sub_analysis_type;
 
 	/* output datasets: */
@@ -34,7 +37,10 @@
 	FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
 	FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
-	FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL);
 	FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec");
 	FetchData((void**)&u_g_obs,NULL,NULL,UGOBS,"Vector","Vec");
+	FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
+	analysis_type=AnalysisTypeAsEnum(analysis_type_string);
+	FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
+	sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
 
 	/*Fetch inputs: */
@@ -43,5 +49,5 @@
 
 	/*!Call core code: */
-	Misfitx(&J, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type);
+	Misfitx(&J, elements,nodes,loads,materials,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
 
 	/*write output : */
@@ -55,5 +61,7 @@
 	xfree((void**)&u_g);
 	xfree((void**)&u_g_obs);
-	delete inputs
+	delete inputs;
+	xfree((void**)&analysis_type_string);
+	xfree((void**)&sub_analysis_type_string);
 
 	/*end module: */
Index: issm/trunk/src/mex/Misfit/Misfit.h
===================================================================
--- issm/trunk/src/mex/Misfit/Misfit.h	(revision 464)
+++ issm/trunk/src/mex/Misfit/Misfit.h	(revision 465)
@@ -25,4 +25,6 @@
 #define UGOBS (mxArray*)prhs[6]
 #define INPUTS (mxArray*)prhs[7]
+#define ANALYSIS (mxArray*)prhs[8]
+#define SUBANALYSIS (mxArray*)prhs[9]
 
 /* serial output macros: */
@@ -33,5 +35,5 @@
 #define NLHS  1
 #undef NRHS
-#define NRHS  8
+#define NRHS  10
 
 
Index: issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp
===================================================================
--- issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp	(revision 464)
+++ issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp	(revision 465)
@@ -32,13 +32,5 @@
 
 	/*Create elements, nodes and materials: */
-	CreateElementsNodesAndMaterials(&elements,&nodes,&materials,model,MODEL);
-
-	/*Create constraints: */
-	CreateConstraints(&constraints,model,MODEL);
-
-	/*Create loads: */
-	CreateLoads(&loads,model,MODEL);
-	/*Create parameters: */
-	CreateParameters(&parameters,model,MODEL);
+	CreateDataSets(&elements,&nodes,&materials,&constraints, &loads, &parameters, model,MODEL);
 
 	/*Write output data: */
Index: issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp
===================================================================
--- issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp	(revision 464)
+++ issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp	(revision 465)
@@ -16,5 +16,8 @@
 	DataSet* materials=NULL;
 	ParameterInputs* inputs=NULL;
+	char*             analysis_type_string=NULL;
 	int               analysis_type;
+	char*             sub_analysis_type_string=NULL;
+	int               sub_analysis_type;
 
 	/*output: */
@@ -35,5 +38,9 @@
 	
 	/*parameters: */
-	FetchData((void**)&analysis_type,NULL,NULL,mxGetField(PARAMETERS,0,"analysis_type"),"Integer",NULL);
+	FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
+	analysis_type=AnalysisTypeAsEnum(analysis_type_string);
+
+	FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
+	sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
 
 	/*Fetch inputs: */
@@ -42,5 +49,5 @@
 
 	/*!Generate internal degree of freedom numbers: */
-	PenaltyConstraintsx(&converged, &num_unstable_constraints, elements,nodes,loads,materials,inputs,analysis_type); 
+	PenaltyConstraintsx(&converged, &num_unstable_constraints, elements,nodes,loads,materials,inputs,analysis_type,sub_analysis_type); 
 
 	/*write output datasets: */
@@ -56,4 +63,6 @@
 	delete materials;
 	delete inputs;
+	xfree((void**)&analysis_type_string);
+	xfree((void**)&sub_analysis_type_string);
 
 	/*end module: */
Index: issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.h
===================================================================
--- issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.h	(revision 464)
+++ issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.h	(revision 465)
@@ -23,4 +23,6 @@
 #define PARAMETERS (mxArray*)prhs[4]
 #define INPUTS (mxArray*)prhs[5]
+#define ANALYSIS (mxArray*)prhs[6]
+#define SUBANALYSIS (mxArray*)prhs[7]
 
 /* serial output macros: */
@@ -33,5 +35,5 @@
 #define NLHS  3
 #undef NRHS
-#define NRHS  6
+#define NRHS  8
 
 #endif  /* _PENALTYCONSTRAINTS_H */
Index: issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp
===================================================================
--- issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp	(revision 464)
+++ issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp	(revision 465)
@@ -21,4 +21,6 @@
 	int               analysis_type;
 	char*             analysis_type_string=NULL;
+	int               sub_analysis_type;
+	char*             sub_analysis_type_string=NULL;
 	
 	/*Boot module: */
@@ -42,4 +44,7 @@
 	analysis_type=AnalysisTypeAsEnum(analysis_type_string);
 
+	FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
+	sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
+
 
 	/*Fetch inputs: */
@@ -48,5 +53,5 @@
 
 	/*!Generate stiffnesses from penalties: */
-	PenaltySystemMatricesx(Kgg, pg,elements,nodes,loads,materials,kflag,pflag,inputs,analysis_type); 
+	PenaltySystemMatricesx(Kgg, pg,elements,nodes,loads,materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); 
 
 	/*write output datasets: */
@@ -63,4 +68,5 @@
 	VecFree(&pg);
 	xfree((void**)&analysis_type_string);
+	xfree((void**)&sub_analysis_type_string);
 
 	/*end module: */
Index: issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.h
===================================================================
--- issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.h	(revision 464)
+++ issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.h	(revision 465)
@@ -26,4 +26,5 @@
 #define INPUTS (mxArray*)prhs[7]
 #define ANALYSIS (mxArray*)prhs[8]
+#define SUBANALYSIS (mxArray*)prhs[9]
 
 /* serial output macros: */
@@ -35,5 +36,5 @@
 #define NLHS  2
 #undef NRHS
-#define NRHS  9
+#define NRHS  10
 
 
Index: issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp
===================================================================
--- issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp	(revision 464)
+++ issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp	(revision 465)
@@ -21,4 +21,6 @@
 	char*             analysis_type_string=NULL;
 	int               analysis_type;
+	char*             sub_analysis_type_string=NULL;
+	int               sub_analysis_type;
 	
 	/* output datasets: */
@@ -46,4 +48,7 @@
 	analysis_type=AnalysisTypeAsEnum(analysis_type_string);
 
+	FetchData((void**)&sub_analysis_type_string,NULL,NULL,SUBANALYSIS,"String",NULL);
+	sub_analysis_type=AnalysisTypeAsEnum(sub_analysis_type_string);
+
 	/*Fetch inputs: */
 	inputs=new ParameterInputs;
@@ -51,5 +56,5 @@
 
 	/*!Generate internal degree of freedom numbers: */
-	SystemMatricesx(&Kgg, &pg,elements,nodes,loads,materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type); 
+	SystemMatricesx(&Kgg, &pg,elements,nodes,loads,materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
 
 	/*write output datasets: */
@@ -60,4 +65,5 @@
 	/*Free ressources: */
 	xfree((void**)&analysis_type_string);
+	xfree((void**)&sub_analysis_type_string);
 	delete elements;
 	delete nodes;
Index: issm/trunk/src/mex/SystemMatrices/SystemMatrices.h
===================================================================
--- issm/trunk/src/mex/SystemMatrices/SystemMatrices.h	(revision 464)
+++ issm/trunk/src/mex/SystemMatrices/SystemMatrices.h	(revision 465)
@@ -24,4 +24,5 @@
 #define INPUTS (mxArray*)prhs[5]
 #define ANALYSIS (mxArray*)prhs[6]
+#define SUBANALYSIS (mxArray*)prhs[7]
 
 /* serial output macros: */
@@ -33,5 +34,5 @@
 #define NLHS  2
 #undef NRHS
-#define NRHS  7
+#define NRHS  8
 
 
