Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 9012)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 9013)
@@ -541,5 +541,6 @@
 	NpartEnum,
 	QmuMassFluxNumProfilesEnum,
-	PartEnum
+	PartEnum,
+	MaxSteadystateIterationsEnum
 };
 
Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp
===================================================================
--- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp	(revision 9012)
+++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp	(revision 9013)
@@ -483,4 +483,5 @@
 		case QmuMassFluxNumProfilesEnum : return "QmuMassFluxNumProfiles";
 		case PartEnum : return "Part";
+		case MaxSteadystateIterationsEnum : return "MaxSteadystateIterations";
 		default : return "unknown";
 
Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 9012)
+++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 9013)
@@ -39,4 +39,5 @@
 	parameters->AddObject(new DoubleParam(EpsAbsEnum,iomodel->eps_abs));
 	parameters->AddObject(new IntParam(MaxNonlinearIterationsEnum,(IssmInt)iomodel->max_nonlinear_iterations));
+	parameters->AddObject(new IntParam(MaxSteadystateIterationsEnum,(IssmInt)iomodel->max_steadystate_iterations));
 	parameters->AddObject(new DoubleParam(EpsvelEnum,iomodel->epsvel));
 	parameters->AddObject(new DoubleParam(YtsEnum,iomodel->yts));
Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp
===================================================================
--- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp	(revision 9012)
+++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp	(revision 9013)
@@ -481,4 +481,5 @@
 	else if (strcmp(name,"QmuMassFluxNumProfiles")==0) return QmuMassFluxNumProfilesEnum;
 	else if (strcmp(name,"Part")==0) return PartEnum;
+	else if (strcmp(name,"MaxSteadystateIterations")==0) return MaxSteadystateIterationsEnum;
 	else _error_("Enum %s not found",name);
 
Index: /issm/trunk/src/c/objects/IoModel.cpp
===================================================================
--- /issm/trunk/src/c/objects/IoModel.cpp	(revision 9012)
+++ /issm/trunk/src/c/objects/IoModel.cpp	(revision 9013)
@@ -185,4 +185,5 @@
 	IoModelFetchData(&this->eps_abs,iomodel_handle,EpsAbsEnum);
 	IoModelFetchData(&this->max_nonlinear_iterations,iomodel_handle,MaxNonlinearIterationsEnum);
+	IoModelFetchData(&this->max_steadystate_iterations,iomodel_handle,MaxSteadystateIterationsEnum);
 	IoModelFetchData(&this->dt,iomodel_handle,DtEnum);
 	IoModelFetchData(&this->ndt,iomodel_handle,NdtEnum);
@@ -355,4 +356,5 @@
 	this->eps_abs=0;
 	this->max_nonlinear_iterations=0;
+	this->max_steadystate_iterations=0;
 	this->dt=0;
 	this->ndt=0;
Index: /issm/trunk/src/c/objects/IoModel.h
===================================================================
--- /issm/trunk/src/c/objects/IoModel.h	(revision 9012)
+++ /issm/trunk/src/c/objects/IoModel.h	(revision 9013)
@@ -150,4 +150,5 @@
 		double  eps_abs;
 		double  max_nonlinear_iterations;
+		double  max_steadystate_iterations;
 		double  dt,ndt;
 		int     time_adapt;
Index: /issm/trunk/src/c/solutions/steadystate_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/steadystate_core.cpp	(revision 9012)
+++ /issm/trunk/src/c/solutions/steadystate_core.cpp	(revision 9013)
@@ -20,5 +20,5 @@
 	int dim;
 	int solution_type;
-	int max_its;
+	int max_steadystate_iterations;
 	bool control_analysis;
 	
@@ -27,5 +27,5 @@
 	femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-	femmodel->parameters->FindParam(&max_its,MaxNonlinearIterationsEnum);
+	femmodel->parameters->FindParam(&max_steadystate_iterations,MaxSteadystateIterationsEnum);
 
 	/*intialize counters: */
@@ -44,6 +44,6 @@
 			if(steadystateconvergence(femmodel)) break;
 		}
-		if(step>max_its){
-			_printf_(VerboseSolution(),"%s%i%s\n","   maximum number of iterations ",max_its," reached");
+		if(step>max_steadystate_iterations){
+			_printf_(VerboseSolution(),"%s%i%s\n","   maximum number steadystate iterations ",max_steadystate_iterations," reached");
 			break;
 		}
Index: /issm/trunk/src/m/classes/model.m
===================================================================
--- /issm/trunk/src/m/classes/model.m	(revision 9012)
+++ /issm/trunk/src/m/classes/model.m	(revision 9013)
@@ -193,4 +193,5 @@
 		 eps_abs                  = {0,true,'Double'};
 		 max_nonlinear_iterations = {0,true,'Double'};
+		 max_steadystate_iterations = {0,true,'Double'};
 		 sparsity                 = {0,true,'Double'};
 		 connectivity             = {0,true,'Integer'};
@@ -660,4 +661,7 @@
 			 %maximum of non-linear iterations.
 			 md.max_nonlinear_iterations=100;
+			 
+			 %maximum of steady state iterations
+			 md.max_steadystate_iterations=100;
 
 			 %sparsity
Index: /issm/trunk/src/m/enum/MaxSteadystateIterationsEnum.m
===================================================================
--- /issm/trunk/src/m/enum/MaxSteadystateIterationsEnum.m	(revision 9013)
+++ /issm/trunk/src/m/enum/MaxSteadystateIterationsEnum.m	(revision 9013)
@@ -0,0 +1,11 @@
+function macro=MaxSteadystateIterationsEnum()
+%MAXSTEADYSTATEITERATIONSENUM - Enum of MaxSteadystateIterations
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh
+%            Please read src/c/EnumDefinitions/README for more information
+%
+%   Usage:
+%      macro=MaxSteadystateIterationsEnum()
+
+macro=StringToEnum('MaxSteadystateIterations');
Index: /issm/trunk/src/m/model/extrude.m
===================================================================
--- /issm/trunk/src/m/model/extrude.m	(revision 9012)
+++ /issm/trunk/src/m/model/extrude.m	(revision 9013)
@@ -224,4 +224,5 @@
 md.surface=project3d(md,'vector',md.surface,'type','node');
 md.thickness=project3d(md,'vector',md.thickness,'type','node');
+md.thickness_coeff=project3d(md,'vector',md.thickness_coeff,'type','node');
 md.bed=project3d(md,'vector',md.bed,'type','node');
 md.nodeonboundary=project3d(md,'vector',md.nodeonboundary,'type','node');
Index: /issm/trunk/src/m/model/partition/AreaAverageOntoPartition.m
===================================================================
--- /issm/trunk/src/m/model/partition/AreaAverageOntoPartition.m	(revision 9012)
+++ /issm/trunk/src/m/model/partition/AreaAverageOntoPartition.m	(revision 9013)
@@ -1,10 +1,40 @@
-function partvector=AreaAverageOntoPartition(md,vector)
+function partvector=AreaAverageOntoPartition(md,vector,layer)
 %AREAAVERAGEONTOPARTITION  compute partition values for a certain vector expressed on the vertices of the mesh. Use area weighted average.
 %
 %   Usage: average=AreaAverageOntoPartition(md,vector)
+%           average=AreaAverageOntoPartition(md,vector,layer) %if in 3D, chose which layer is partitioned
 %
 
-%ok, first check that part is Matlab matlab indexed
+%some checks
+if md.dim==3,
+	if nargin~=3,
+		error('layer should be provided onto which Area Averaging occurs');
+	end
+	%save 3D model
+	md3d=md;
+	
+	md.elements=md.elements2d;
+	md.x=md.x2d;
+	md.y=md.y2d;
+	md.numberofnodes=md.numberofnodes2d;
+	md.numberofelements=md.numberofelements2d;
+	md.vwgt=[];
+	md.nodeconnectivity=[];
+
+	%run connectivity routine
+	md=adjacency(md);
+
+	%finally, project vector: 
+	vector=project2d(md3d,vector,layer);
+	md.part=project2d(md3d,md3d.part,layer);
+end
+
+%ok, first check that part is Matlab indexed
 part=md.part+1;
+
+%some check: 
+if md.npart~=max(part),
+	error('AreaAverageOntoPartition error message: ''npart'' should be equal to max(md.part)');
+end
 
 %initialize output
@@ -17,2 +47,7 @@
 	partvector(i)=sum(weightedvector(pos))/sum(md.vwgt(pos));
 end
+
+%in 3D, restore 3D model:
+if md.dim==3,
+	md=md3d;
+end
Index: /issm/trunk/src/m/model/partition/partitioner.m
===================================================================
--- /issm/trunk/src/m/model/partition/partitioner.m	(revision 9012)
+++ /issm/trunk/src/m/model/partition/partitioner.m	(revision 9013)
@@ -30,4 +30,17 @@
 npart=getfieldvalue(options,'npart');
 recomputeadjacency=getfieldvalue(options,'recomputeadjacency');
+
+if(md.dim==3),
+	%partitioning essentially happens in 2D. So partition in 2D, then 
+	%extrude the partition vector vertically. 
+	md3d=md; %save  for later
+	md.elements=md.elements2d;
+	md.x=md.x2d;
+	md.y=md.y2d;
+	md.numberofnodes=md.numberofnodes2d;
+	md.numberofelements=md.numberofelements2d;
+	md.vwgt=[];
+	md.nodeconnectivity=[];
+end
 
 %adjacency matrix if needed:
@@ -87,3 +100,9 @@
 end
 
+%extrude if we are in 3D:
+if md.dim==3,
+	md=md3d;
+	part=project3d(md,'vector',part','type','node');
+end
+
 md.part=part;
Index: /issm/trunk/src/m/model/tres.m
===================================================================
--- /issm/trunk/src/m/model/tres.m	(revision 9012)
+++ /issm/trunk/src/m/model/tres.m	(revision 9013)
@@ -85,8 +85,8 @@
 	md.temperature=PatchToVec(md.results.SteadystateSolution.Temperature);
 	md.basal_melting_rate=PatchToVec(md.results.SteadystateSolution.BasalMeltingRate);
-	
+
 	if md.control_analysis==1,
-		if control_type==md.control_type
-			md.(EnumToModelField(control_type))=PatchToVec(md.results.DiagnosticSolution.(EnumToString(control_type)));
+		for control_type=md.control_type
+			md.(EnumToModelField(control_type))=PatchToVec(md.results.SteadystateSolution.(EnumToString(control_type)));
 		end
 	end