Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15894)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15895)
@@ -5245,5 +5245,4 @@
 	IssmDouble  enthalpy, enthalpyup;
 	IssmDouble  pressure, pressureup;
-	IssmDouble  meltrate, watercolumn;
 	IssmDouble  temperature, waterfraction;
 	IssmDouble  latentheat;
@@ -5282,4 +5281,6 @@
 	SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type);
 	_assert_(numindices==numindicesup); 
+    
+	IssmDouble  meltrate[numindices], watercolumn[numindices];
 
 	/*Ok, get meltrates now from basal conditions*/
@@ -5287,18 +5288,20 @@
 	GaussPenta* gaussup=new GaussPenta();
 	for(ig=0;ig<numindices;ig++){
-		gauss->GaussNode(this->element_type,indices[ig]);
+		gauss->GaussVertex(this->element_type,indices[ig]);
 		gaussup->GaussNode(this->element_type,indicesup[ig]);
 
-		watercolumn_input->GetInputValue(&watercolumn, gauss);
+        // TODO: make sure that no node is computed twice (insert mask)
+
+		watercolumn_input->GetInputValue(&watercolumn[indices[ig]], gauss);
 		enthalpy_input->GetInputValue(&enthalpy, gauss);
 		pressure_input->GetInputValue(&pressure, gauss);
 
 		/*Calculate basal meltrate*/
-		if((watercolumn>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
+		if((watercolumn[indices[ig]]>0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
 			enthalpy=matpar->PureIceEnthalpy(pressure);
 		}
 		else if(enthalpy<matpar->PureIceEnthalpy(pressure)){
-			meltrate=0.;   //TODO: set zero meltrate and watercolumn in model
-			watercolumn=0.;
+			meltrate[indices[ig]]=0.;   //TODO: set zero meltrate and watercolumn in model
+			watercolumn[indices[ig]]=0.;
 			return;
 		}
@@ -5330,14 +5333,16 @@
 		vz_input->GetInputValue(&vz,gauss);
 		basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
-		meltrate=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
+		meltrate[indices[ig]]=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
 
 		/*Update water column*/
 		this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
 		if(reCast<bool,IssmDouble>(dt))
-		 watercolumn+=dt*meltrate;
+            watercolumn[indices[ig]]+=dt*meltrate[indices[ig]];
 		else
-		 watercolumn=meltrate;
-		// TODO: feed meltrate & watercolumn back to model
+            watercolumn[indices[ig]]=meltrate[indices[ig]];
 	}  
+    /*update meltrate and watercolumn*/
+    this->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumn, P1Enum));
+    this->inputs->AddInput(new PentaInput(BasalMeltrateEnum, meltrate, P1Enum));
 
 	/*Clean up and return*/
@@ -5354,11 +5359,11 @@
 
     /*Intermediaries*/
-    const int numdof=NDOF1*NUMVERTICES;
     int ig;
     bool isenthalpy;
-    IssmDouble waterfraction, temperature;
-    IssmDouble enthalpy, pressure; 
+    IssmDouble waterfraction_array[NUMVERTICES], temperature[NUMVERTICES];
+    IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; 
     IssmDouble latentheat, dt;
-
+    GaussPenta* gauss;
+    
     Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
 	Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
@@ -5368,22 +5373,25 @@
 	parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
 	if(!isenthalpy) return;       
+    
     this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-    
-    GaussPenta* gauss=new GaussPenta(2,3);
-    for(ig=gauss->begin();ig<gauss->end();ig++){ 
-        gauss->GaussPoint(ig);
-        
-        enthalpy_input->GetInputValue(&enthalpy, gauss);
-        pressure_input->GetInputValue(&pressure, gauss);
-        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy,pressure); 
+    latentheat=matpar->GetLatentHeat();
+    gauss=new GaussPenta(2,3);
+    for(ig=0;ig<NUMVERTICES;ig++){ 
+        gauss->GaussVertex(ig);
+        /*TODO: Check whether Vertex has been drained already*/
+        enthalpy_input->GetInputValue(&enthalpy[ig], gauss);
+        pressure_input->GetInputValue(&pressure[ig], gauss);
+        matpar->EnthalpyToThermal(&temperature[ig], &waterfraction[ig], enthalpy[ig],pressure[ig]); 
     
         /*drain water fraction*/
-        waterfraction-=dt*DrainageFunctionWaterfraction(waterfraction);
-        if(waterfraction<0) waterfraction=0.;
+        waterfraction[ig]-=dt*DrainageFunctionWaterfraction(waterfraction[ig]);
+        if(waterfraction[ig]<0) waterfraction[ig]=0.;
         /*update enthalpy*/
-        latentheat=matpar->GetLatentHeat();
-        matpar->ThermalToEnthalpy(&enthalpy, temperature, waterfraction, pressure);
-        //TODO feed result back into model
+        matpar->ThermalToEnthalpy(&enthalpy[ig], temperature[ig], waterfraction[ig], pressure[ig]);        
     }
+    /*feed updated results back into model*/
+    this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
+    this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
+    this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum));    
 }
 /*}}}*/