Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24651)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24652)
@@ -2986,5 +2986,5 @@
 	const int NUM_VERTICES_MONTHS_PER_YEAR	= NUM_VERTICES * 12;
 
-	int        	i;
+	int        	i,vertexlids[MAXVERTICES];;
 	IssmDouble* smb=xNew<IssmDouble>(NUM_VERTICES);		// surface mass balance
 	IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES);		// melting comp. of surface mass balance
@@ -3003,4 +3003,7 @@
 	IssmDouble inv_twelve=1./12.;								//factor for monthly average
 	IssmDouble time,yts,time_yr;
+
+	/*Get vertex Lids for later*/
+	this->GetVerticesLidList(&vertexlids[0]);
 
 	/*Get material parameters :*/
@@ -3089,19 +3092,46 @@
 					/*FIXME: look at other function Element::PositiveDegreeDay and propagate change! Just assert for now*/
 					PentaInput2* temp_input = xDynamicCast<PentaInput2*>(this->GetInput2(TemperatureEnum)); _assert_(temp_input);
-					_assert_(temp_input->GetInputInterpolationType()==P1Enum);
-
-					GetInputListOnVertices(&s[0],TemperatureEnum);
-					yearlytemperatures[0] = s[0];
-					yearlytemperatures[1] = s[1];
-					yearlytemperatures[2] = s[2];
-					this->AddInput2(TemperatureEnum,&yearlytemperatures[0],P1Enum);
+					switch(temp_input->GetInputInterpolationType()){
+						case P1Enum:
+							temp_input->element_values[3] = yearlytemperatures[3];
+							temp_input->element_values[4] = yearlytemperatures[4];
+							temp_input->element_values[5] = yearlytemperatures[5];
+							temp_input->SetInput(P1Enum,NUM_VERTICES,&vertexlids[0],temp_input->element_values);
+							break;
+						case P1DGEnum:
+						case P1xP2Enum:
+						case P1xP3Enum:
+						case P1xP4Enum:
+							temp_input->element_values[3] = yearlytemperatures[3];
+							temp_input->element_values[4] = yearlytemperatures[4];
+							temp_input->element_values[5] = yearlytemperatures[5];
+							temp_input->SetInput(temp_input->GetInputInterpolationType(),this->lid,this->GetNumberOfNodes(temp_input->GetInputInterpolationType()),temp_input->element_values);
+							break;
+						default:
+							_error_("Interpolation "<<EnumToStringx(temp_input->GetInputInterpolationType())<<" not supported yet");
+					}
+
 					if(isenthalpy){
 						/*Convert that to enthalpy for the enthalpy model*/
-						IssmDouble enthalpy[6];
-						GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
-						ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.);
-						ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[4],0.,0.);
-						ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[5],0.,0.);
-						this->AddInput2(EnthalpyEnum,&enthalpy[0],P1Enum);
+						PentaInput2* enth_input = xDynamicCast<PentaInput2*>(this->GetInput2(EnthalpyEnum)); _assert_(enth_input);
+						switch(enth_input->GetInputInterpolationType()){
+							case P1Enum:
+								ThermalToEnthalpy(&enth_input->element_values[3],yearlytemperatures[3],0.,0.);
+								ThermalToEnthalpy(&enth_input->element_values[4],yearlytemperatures[4],0.,0.);
+								ThermalToEnthalpy(&enth_input->element_values[5],yearlytemperatures[5],0.,0.);
+								enth_input->SetInput(P1Enum,NUM_VERTICES,&vertexlids[0],enth_input->element_values);
+								break;
+							case P1DGEnum:
+							case P1xP2Enum:
+							case P1xP3Enum:
+							case P1xP4Enum:
+								ThermalToEnthalpy(&enth_input->element_values[3],yearlytemperatures[3],0.,0.);
+								ThermalToEnthalpy(&enth_input->element_values[4],yearlytemperatures[4],0.,0.);
+								ThermalToEnthalpy(&enth_input->element_values[5],yearlytemperatures[5],0.,0.);
+								enth_input->SetInput(enth_input->GetInputInterpolationType(),this->lid,this->GetNumberOfNodes(enth_input->GetInputInterpolationType()),enth_input->element_values);
+								break;
+							default:
+								_error_("Interpolation "<<EnumToStringx(temp_input->GetInputInterpolationType())<<" not supported yet");
+						}
 					}
 				}