Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15892)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15893)
@@ -5236,5 +5236,5 @@
 	/* Intermediaries */
 	bool        isenthalpy;
-	int         i,j,analysis_type,numindices,numindicesup;
+	int         i,j,ig,analysis_type,numindices,numindicesup;
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
@@ -5283,10 +5283,10 @@
 	_assert_(numindices==numindicesup); 
 
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	/*Ok, get meltrates now from basal conditions*/
 	GaussPenta* gauss=new GaussPenta();
 	GaussPenta* gaussup=new GaussPenta();
-	for(i=0;i<numindices;i++){
-		gauss->GaussNode(this->element_type,indices[i]);
-		gaussup->GaussNode(this->element_type,indicesup[i]);
+	for(ig=0;ig<numindices;ig++){
+		gauss->GaussNode(this->element_type,indices[ig]);
+		gaussup->GaussNode(this->element_type,indicesup[ig]);
 
 		watercolumn_input->GetInputValue(&watercolumn, gauss);
@@ -5308,5 +5308,5 @@
 		pressure_input->GetInputValue(&pressureup, gaussup);    
 		if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
-		 for(i=0;i<3;i++) vec_heatflux[i]=0.;
+            for(i=0;i<3;i++) vec_heatflux[i]=0.;
 		}
 		else{
@@ -5346,4 +5346,44 @@
 	delete gauss;
 	delete gaussup;
+}
+/*}}}*/
+/*FUNCTION Penta::DrainWaterfraction{{{*/
+void Penta::DrainWaterfraction(void){
+    
+// TODO: create vector to mark whether node has been drained already i.o. to not drain nodes multiple times
+
+    /*Intermediaries*/
+    const int numdof=NDOF1*NUMVERTICES;
+    int ig;
+    bool isenthalpy;
+    IssmDouble waterfraction, temperature;
+    IssmDouble enthalpy, pressure; 
+    IssmDouble latentheat, dt;
+
+    Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
+	Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
+    Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
+    
+    /*Check wether enthalpy is activated*/
+	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); 
+    
+        /*drain water fraction*/
+        waterfraction-=dt*DrainageFunctionWaterfraction(waterfraction);
+        if(waterfraction<0) waterfraction=0.;
+        /*update enthalpy*/
+        latentheat=matpar->GetLatentHeat();
+        matpar->ThermalToEnthalpy(&enthalpy, temperature, waterfraction, pressure);
+        //TODO feed result back into model
+    }
 }
 /*}}}*/