Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 21541)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 21542)
@@ -98,6 +98,9 @@
 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");
 }/*}}}*/
@@ -121,4 +124,8 @@
 	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;
@@ -126,5 +133,5 @@
 		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: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 21541)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 21542)
@@ -9,5 +9,7 @@
 
 	/*Intermediary*/
-	int finiteelement = P1Enum;
+	int finiteelement;
+	iomodel->FindConstant(&finiteelement,"md.thermal.fe");
+	_assert_(finiteelement==P1Enum); 
 
 	/*Only 3d mesh supported*/
@@ -62,5 +64,7 @@
 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");
@@ -79,5 +83,7 @@
 
 	/*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++){
Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 21541)
+++ /issm/trunk-jpl/src/c/classes/IoModel.cpp	(revision 21542)
@@ -656,4 +656,6 @@
 						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){
@@ -797,4 +799,6 @@
 						} 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));
Index: /issm/trunk-jpl/src/m/classes/thermal.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/thermal.m	(revision 21541)
+++ /issm/trunk-jpl/src/m/classes/thermal.m	(revision 21542)
@@ -15,4 +15,5 @@
 		isenthalpy        = 0;
 		isdynamicbasalspc = 0;
+		fe                = 'P1';
 		requested_outputs = {};
 	end
@@ -66,4 +67,7 @@
 			self.isdynamicbasalspc=0;
 
+			%Linear elements by default
+			self.fe='P1';
+
 			%default output
 			self.requested_outputs={'default'};
@@ -76,4 +80,5 @@
 			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),
 
@@ -110,4 +115,5 @@
 			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');
 
@@ -122,4 +128,5 @@
 			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');
 
Index: /issm/trunk-jpl/src/m/classes/thermal.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/thermal.py	(revision 21541)
+++ /issm/trunk-jpl/src/m/classes/thermal.py	(revision 21542)
@@ -23,4 +23,5 @@
 		self.isenthalpy        = 0
 		self.isdynamicbasalspc = 0;
+		self.fe                = 'P1';
 		self.requested_outputs = []
 
@@ -79,5 +80,8 @@
 
 		#will basal boundary conditions be set dynamically
-		self.isdynamicbasalspc=0;
+		self.isdynamicbasalspc=0
+
+		#Finite element interpolation
+		self.fe='P1'
 
 		#default output
@@ -124,4 +128,5 @@
 		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');