Index: ../trunk-jpl/src/m/classes/thermal.py
===================================================================
--- ../trunk-jpl/src/m/classes/thermal.py	(revision 21541)
+++ ../trunk-jpl/src/m/classes/thermal.py	(revision 21542)
@@ -22,6 +22,7 @@
 		self.penalty_factor    = 0
 		self.isenthalpy        = 0
 		self.isdynamicbasalspc = 0;
+		self.fe                = 'P1';
 		self.requested_outputs = []
 
 		#set defaults
@@ -78,8 +79,11 @@
 		self.isenthalpy=0
 
 		#will basal boundary conditions be set dynamically
-		self.isdynamicbasalspc=0;
+		self.isdynamicbasalspc=0
 
+		#Finite element interpolation
+		self.fe='P1'
+
 		#default output
 		self.requested_outputs=['default']
 		return self
@@ -123,6 +127,7 @@
 		WriteData(fid,prefix,'object',self,'fieldname','penalty_lock','format','Integer')
 		WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double')
 		WriteData(fid,prefix,'object',self,'fieldname','isenthalpy','format','Boolean')
+		WriteData(fid,prefix,'object',self,'fieldname','fe','format','String');
 		WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean');
 
 		#process requested outputs
Index: ../trunk-jpl/src/m/classes/thermal.m
===================================================================
--- ../trunk-jpl/src/m/classes/thermal.m	(revision 21541)
+++ ../trunk-jpl/src/m/classes/thermal.m	(revision 21542)
@@ -14,6 +14,7 @@
 		penalty_factor    = 0;
 		isenthalpy        = 0;
 		isdynamicbasalspc = 0;
+		fe                = 'P1';
 		requested_outputs = {};
 	end
 	methods
@@ -65,6 +66,9 @@
 			%will basal boundary conditions be set dynamically
 			self.isdynamicbasalspc=0;
 
+			%Linear elements by default
+			self.fe='P1';
+
 			%default output
 			self.requested_outputs={'default'};
 		end % }}}
@@ -75,6 +79,7 @@
 
 			md = checkfield(md,'fieldname','thermal.stabilization','numel',[1],'values',[0 1 2]);
 			md = checkfield(md,'fieldname','thermal.spctemperature','Inf',1,'timeseries',1,'>=',0);
+			md = checkfield(md,'fieldname','thermal.fe','values',{'P1','P1xP2'});
 			if (ismember('EnthalpyAnalysis',analyses) & md.thermal.isenthalpy & dimension(md.mesh)==3),
 
 				%Make sure the spc are less than melting point
@@ -109,6 +114,7 @@
 			fielddisplay(self,'penalty_factor','scaling exponent (default is 3)');
 			fielddisplay(self,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
 			fielddisplay(self,'isdynamicbasalspc',['enable dynamic setting of basal forcing. required for enthalpy formulation (default is 0)']);
+			fielddisplay(self,'fe','Finite Element type: ''P1'' (default), ''P1xP2''');
 			fielddisplay(self,'requested_outputs','additional outputs requested');
 
 		end % }}}
@@ -121,6 +127,7 @@
 			WriteData(fid,prefix,'object',self,'fieldname','penalty_lock','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','penalty_factor','format','Double');
 			WriteData(fid,prefix,'object',self,'fieldname','isenthalpy','format','Boolean');
+			WriteData(fid,prefix,'object',self,'fieldname','fe','format','String');
 			WriteData(fid,prefix,'object',self,'fieldname','isdynamicbasalspc','format','Boolean');
 
 			%process requested outputs
Index: ../trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/IoModel.cpp	(revision 21541)
+++ ../trunk-jpl/src/c/classes/IoModel.cpp	(revision 21542)
@@ -656,6 +656,8 @@
 						this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
 					} else if(strcmp(record_name,"md.flowequation.fe_FS")==0){
 						this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
+					} else if(strcmp(record_name,"md.thermal.fe")==0){
+						this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
 					} else if(strcmp(record_name,"md.groundingline.migration")==0){
 						this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
 					} else if(strcmp(record_name,"md.masstransport.hydrostatic_adjustment")==0){
@@ -796,6 +798,8 @@
 							this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
 						} else if(strcmp(record_name,"md.flowequation.fe_FS")==0){
 							this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
+						} else if(strcmp(record_name,"md.thermal.fe")==0){
+							this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
 						} else if(strcmp(record_name,"md.groundingline.migration")==0){
 							this->AddConstant(new IoConstant(StringToEnumx(string),record_name));
 						} else if(strcmp(record_name,"md.masstransport.hydrostatic_adjustment")==0){
Index: ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 21541)
+++ ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 21542)
@@ -97,8 +97,11 @@
 }/*}}}*/
 void EnthalpyAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+
 	if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
-	::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,P1Enum);
+	::CreateNodes(nodes,iomodel,EnthalpyAnalysisEnum,finiteelement);
 	iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
 }/*}}}*/
 int  EnthalpyAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
@@ -120,12 +123,16 @@
 	/*Fetch data needed: */
 	iomodel->FetchData(3,"md.initialization.temperature","md.initialization.waterfraction","md.initialization.pressure");
 
+	/*Finite element type*/
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+
 	/*Update elements: */
 	int counter=0;
 	for(int i=0;i<iomodel->numberofelements;i++){
 		if(iomodel->my_elements[i]){
 			Element* element=(Element*)elements->GetObjectByOffset(counter);
-			element->Update(i,iomodel,analysis_counter,analysis_type,P1Enum);
+			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
 			counter++;
 		}
 	}
Index: ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 21541)
+++ ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 21542)
@@ -8,7 +8,9 @@
 void ThermalAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
 	/*Intermediary*/
-	int finiteelement = P1Enum;
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+	_assert_(finiteelement==P1Enum); 
 
 	/*Only 3d mesh supported*/
 	if(iomodel->domaintype!=Domain3DEnum) _error_("not supported yet");
@@ -61,7 +63,9 @@
 }/*}}}*/
 void ThermalAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
-	int finiteelement = P1Enum;
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+	_assert_(finiteelement==P1Enum); 
 	
 	if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
 	::CreateNodes(nodes,iomodel,ThermalAnalysisEnum,finiteelement);
@@ -78,7 +82,9 @@
 	if(iomodel->domaintype==Domain2DhorizontalEnum)return;
 
 	/*Update elements: */
-	int finiteelement = P1Enum;
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+	_assert_(finiteelement==P1Enum); 
 	int counter=0;
 	for(int i=0;i<iomodel->numberofelements;i++){
 		if(iomodel->my_elements[i]){