Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 17013)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 17014)
@@ -403,11 +403,15 @@
 
 	/*Intermediaries*/
-	int         stabilization;
+	int         i, stabilization;
 	IssmDouble  Jdet,phi,dt;
-	IssmDouble  enthalpy;
-	IssmDouble  kappa,tau_parameter,diameter;
+	IssmDouble  enthalpy, Hpmp;
+	IssmDouble  enthalpypicard, d1enthalpypicard[3];
+	IssmDouble  pressure, d1pressure[3], d2pressure;
+	IssmDouble  waterfractionpicard;
+	IssmDouble  kappa,tau_parameter,diameter,kappa_w;
 	IssmDouble  u,v,w;
-	IssmDouble  scalar_def,scalar_transient;
+	IssmDouble  scalar_def, scalar_sens ,scalar_transient;
 	IssmDouble* xyz_list = NULL;
+	IssmDouble  d1H_d1P, d1P2;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -423,5 +427,9 @@
 	element->GetVerticesCoordinates(&xyz_list);
 	IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
 	IssmDouble  thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
+	IssmDouble  temperateiceconductivity = element->GetMaterialParameter(MaterialsTemperateiceconductivityEnum);
+	IssmDouble  beta                = element->GetMaterialParameter(MaterialsBetaEnum);
+	IssmDouble  latentheat          = element->GetMaterialParameter(MaterialsLatentheatEnum);
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
 	element->FindParam(&stabilization,ThermalStabilizationEnum);
@@ -429,5 +437,7 @@
 	Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
 	Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input);
-	Input* enthalpy_input = NULL;
+	Input* enthalpypicard_input=element->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input);
+	Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
+	Input* enthalpy_input=NULL;
 	if(reCast<bool,IssmDouble>(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);}
 	if(stabilization==2){
@@ -443,4 +453,6 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
+		
+		/*viscous dissipation*/
 		element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
 
@@ -448,6 +460,25 @@
 		if(dt!=0.) scalar_def=scalar_def*dt;
 
-		/*TODO: add -beta*laplace T_m(p)*/
-		for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
+		for(i=0;i<numnodes;i++) pe->values[i]+=scalar_def*basis[i];
+
+		/*sensible heat flux in temperate ice*/
+		enthalpypicard_input->GetInputValue(&enthalpypicard,gauss);
+		pressure_input->GetInputValue(&pressure,gauss);
+		Hpmp=this->PureIceEnthalpy(element, pressure);
+
+		if(enthalpypicard>=Hpmp){
+			enthalpypicard_input->GetInputDerivativeValue(&d1enthalpypicard[0],xyz_list,gauss);
+			pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
+			d2pressure=0.; // for linear elements, 2nd derivative is zero
+			
+			d1H_d1P=0.;
+			for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
+			d1P2=0.;
+			for(i=0;i<3;i++) d1P2+=pow(d1pressure[i],2.);
+
+			scalar_sens=-beta*((temperateiceconductivity - thermalconductivity)/latentheat*(d1H_d1P + beta*heatcapacity*d1P2))/rho_ice;
+			if(dt!=0.) scalar_sens=scalar_sens*dt;
+			for(i=0;i<numnodes;i++) pe->values[i]+=scalar_sens*basis[i];
+		}		
 
 		/* Build transient now */
@@ -455,5 +486,5 @@
 			enthalpy_input->GetInputValue(&enthalpy, gauss);
 			scalar_transient=enthalpy*Jdet*gauss->weight;
-			for(int i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
+			for(i=0;i<numnodes;i++) pe->values[i]+=scalar_transient*basis[i];
 		}
 
@@ -466,8 +497,8 @@
 			tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa/rho_ice);
 
-			for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
+			for(i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
 
 			if(dt!=0.){
-				for(int i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
+				for(i=0;i<numnodes;i++) pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]);
 			}
 		}
@@ -484,5 +515,5 @@
 ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
 
-	/* Geothermal flux on ice sheet base and basal friction */
+	/* implementation of the basal condition decision chart of Aschwanden 2012, Fig.5 */
 	if(!element->IsOnBed() || element->IsFloating()) return NULL;
 
@@ -527,7 +558,7 @@
 		watercolumn_input->GetInputValue(&watercolumn,gauss);
 
-		if((watercolumn<=0.) && (enthalpy < PureIceEnthalpy(element,pressure))){
+		if((watercolumn<=0.) && (enthalpy<PureIceEnthalpy(element,pressure))){
 			/* the above check is equivalent to 
-			 NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
+			 NOT [(watercolumn>0.) AND (enthalpy<PIE)] AND (enthalpy<PIE)*/
 			geothermalflux_input->GetInputValue(&geothermalflux,gauss);
 
@@ -549,5 +580,5 @@
 			pressure_input->GetInputValue(&pressureup,gaussup);
 			if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
-				// TODO: temperate ice has positive thickness: grad enthalpy*n=0.
+				// do nothing, set grad enthalpy*n=0.
 			}
 			else{
@@ -556,5 +587,5 @@
 		}
 		else{
-			// base cold, but watercolumn positive. Set base to h_pmp.
+			// base cold, but watercolumn positive. Set base to pressure melting point enthalpy
 		}
 	}
@@ -810,13 +841,314 @@
 	}
 }/*}}}*/
+
+
+
+
+
+
 void EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
-
-}/*}}}*/
-void EnthalpyAnalysis::DrainWaterfraction(Element* element){/*{{{*/
-
-}/*}}}*/
+	/*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
+	/* melting rate is positive when melting, negative when refreezing*/
+
+	/* Intermediaries */
+	int         i,is,vertexdown,vertexup,dim=3,numvertices,numsegments;
+	IssmDouble  heatflux;
+	IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim];
+	IssmDouble  temperature, waterfraction;
+	IssmDouble  basalfriction,alpha2;
+	IssmDouble  dt,yts;
+	IssmDouble  melting_overshoot,lambda;
+	IssmDouble  geothermalflux;
+	IssmDouble  vx,vy,vz;
+	IssmDouble *xyz_list      = NULL;
+	IssmDouble *xyz_list_base = NULL;
+	int        *pairindices   = NULL;
+
+	/* Only compute melt rates at the base of grounded ice*/
+	if(!element->IsOnBed() || element->IsFloating()) return;
+
+	/*Fetch parameters and inputs */
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
+	IssmDouble rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	//Input* watercolumn_input      = inputs->GetInput(WatercolumnEnum);                 _assert_(watercolumn_input);
+	Input* enthalpy_input         = element->GetInput(EnthalpyEnum);                    _assert_(enthalpy_input);
+	//  Input* pressure_input         = inputs->GetInput(PressureEnum);                    _assert_(pressure_input);
+	//	Input* basalmeltingrate_input = inputs->GetInput(BasalforcingsMeltingRateEnum);    _assert_(basalmeltingrate_input);
+	Input* geothermalflux_input   = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
+	Input* vx_input               = element->GetInput(VxEnum);                          _assert_(vx_input);
+	Input* vy_input               = element->GetInput(VyEnum);                          _assert_(vy_input);
+	Input* vz_input               = element->GetInput(VzEnum);                          _assert_(vz_input);
+	IssmDouble kappa=EnthalpyDiffusionParameterVolume(element,EnthalpyEnum);     _assert_(kappa>0.);
+	element->NormalBase(&normal_base[0],xyz_list_base);
+	element->VerticalSegmentIndices(&pairindices,&numsegments);
+	IssmDouble* meltingrate_enthalpy = xNew<IssmDouble>(numsegments);
+	IssmDouble* heating = xNew<IssmDouble>(numsegments);	
+
+	/*Build friction element, needed later: */
+	Friction* friction=new Friction(element,dim);
+
+	/******** MELTING RATES  ************************************/
+	numvertices=element->GetNumberOfVertices();
+	IssmDouble* enthalpy = xNew<IssmDouble>(numvertices);
+	IssmDouble* pressure = xNew<IssmDouble>(numvertices);
+	IssmDouble* watercolumn = xNew<IssmDouble>(numvertices);
+	IssmDouble* basalmeltingrate = xNew<IssmDouble>(numvertices);
+	element->GetInputListOnVertices(enthalpy,EnthalpyEnum);
+	element->GetInputListOnVertices(pressure,PressureEnum);
+	element->GetInputListOnVertices(watercolumn,WatercolumnEnum);
+	element->GetInputListOnVertices(basalmeltingrate,BasalforcingsMeltingRateEnum);
+
+	Gauss* gauss=element->NewGauss();
+	
+	for(int is=0;is<numsegments;is++){
+		vertexdown = pairindices[is*2+0];
+		vertexup   = pairindices[is*2+1];
+		gauss->GaussVertex(vertexdown);
+		
+		bool checkpositivethickness=true;
+		_assert_(watercolumn>=0.);
+
+		/*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
+		meltingrate_enthalpy[is]=0.;
+		heating[is]=0.;
+		if((watercolumn[vertexdown]>0.) && (enthalpy[vertexdown]<PureIceEnthalpy(element,pressure[vertexdown]))){
+			/*ensure that no ice is at T<Tm(p), if water layer present*/
+			enthalpy[vertexdown]=element->PureIceEnthalpy(pressure[vertexdown]); 
+		}
+		else if(enthalpy[vertexdown]<element->PureIceEnthalpy(pressure[vertexdown])){
+			/*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
+			checkpositivethickness=false; // cold base, skip next test
+		}
+		else{/*we have a temperate base, go to next test*/}
+
+		if(checkpositivethickness){
+			/*From here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
+			bool istemperatelayer=false;
+			if(enthalpy[vertexup]>=element->PureIceEnthalpy(pressure[vertexup])) istemperatelayer=true;
+			if(istemperatelayer) for(i=0;i<dim;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp
+			else{
+				enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],xyz_list,gauss);
+				for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
+			}
+
+			/*heat flux along normal*/
+			heatflux=0.;
+			for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i];
+
+			/*basal friction*/
+			friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+			vx_input->GetInputValue(&vx,gauss);
+			vy_input->GetInputValue(&vy,gauss);
+			vz_input->GetInputValue(&vz,gauss);
+			basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
+
+			element->EnthalpyToThermal(&temperature,&waterfraction,enthalpy[vertexdown],pressure[vertexdown]);
+			geothermalflux_input->GetInputValue(&geothermalflux,gauss);
+			// -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
+			heating[is]=(heatflux+basalfriction+geothermalflux);
+			meltingrate_enthalpy[is]=heating[is]/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent //????
+		}
+	}
+	// enthalpy might have been changed, update 
+	//element->AddInput(EnthalpyEnum,enthalpy,P1Enum);
+
+	/******** DRAINAGE *****************************************/
+	IssmDouble* drainrate_column = xNew<IssmDouble>(numsegments); //TODO: xDelete?
+	IssmDouble* drainrate_element = xNew<IssmDouble>(numsegments);
+	for(is=0;is<numsegments;is++)	drainrate_column[is]=0.;
+	Element* elementi = element;
+	for(;;){
+		for(is=0;is<numsegments;is++)	drainrate_element[is]=0.;
+		DrainWaterfraction(elementi,drainrate_element); // TODO: make sure every vertex is only drained once
+		for(is=0;is<numsegments;is++)	drainrate_column[is]+=drainrate_element[is];
+
+		if(elementi->IsOnSurface()) break;
+		elementi=elementi->GetUpperElement();			
+	}
+	// add drained water to melting rate
+	for(is=0;is<numsegments;is++) meltingrate_enthalpy[is]+=drainrate_column[is];
+
+	/******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	for(is=0;is<numsegments;is++){
+		vertexdown = pairindices[is*2+0];
+		vertexup   = pairindices[is*2+1];
+		if(reCast<bool,IssmDouble>(dt)){
+			if(watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt<0.){	// prevent too much freeze on			
+				melting_overshoot=watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt;
+				lambda=melting_overshoot/(meltingrate_enthalpy[is]*dt); _assert_(lambda>0); _assert_(lambda<1);
+				basalmeltingrate[vertexdown]=(1.-lambda)*meltingrate_enthalpy[is];
+				watercolumn[vertexdown]=0.;
+				yts=365*24*60*60;
+				enthalpy[vertexdown]+=dt/yts*lambda*heating[is];
+			}
+			else{
+				basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
+				watercolumn[vertexdown]+=dt*meltingrate_enthalpy[is]; 
+			}
+		}
+		else{
+			basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
+			watercolumn[vertexdown]+=meltingrate_enthalpy[is];
+		}		
+		_assert_(watercolumn[vertexdown]>=0.);
+	}
+
+	/*feed updated variables back into model*/
+	element->AddInput(EnthalpyEnum,enthalpy,P1Enum);
+	element->AddInput(WatercolumnEnum,watercolumn,P1Enum);
+	element->AddInput(BasalforcingsMeltingRateEnum,basalmeltingrate,P1Enum);
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
+  	xDelete<IssmDouble>(enthalpy);
+  	xDelete<IssmDouble>(pressure);
+  	xDelete<IssmDouble>(watercolumn);
+  	xDelete<IssmDouble>(basalmeltingrate);
+	xDelete<IssmDouble>(meltingrate_enthalpy);
+	xDelete<IssmDouble>(heating);
+	xDelete<IssmDouble>(drainrate_column);
+	xDelete<IssmDouble>(drainrate_element);
+}/*}}}*/
+
+void EnthalpyAnalysis::DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element){/*{{{*/
+
+	/*Intermediaries*/
+	int iv,is,vertexdown,vertexup,numsegments;	
+	IssmDouble dt, height_element;
+	IssmDouble rho_water, rho_ice;
+	int numvertices = element->GetNumberOfVertices();
+
+	IssmDouble* xyz_list = NULL;
+	IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
+	IssmDouble* pressures = xNew<IssmDouble>(numvertices);
+	IssmDouble* temperatures = xNew<IssmDouble>(numvertices);
+	IssmDouble* waterfractions = xNew<IssmDouble>(numvertices);
+	IssmDouble* deltawaterfractions = xNew<IssmDouble>(numvertices);
+	int        *pairindices   = NULL;
+	
+	rho_ice=element->GetMaterialParameter(MaterialsRhoIceEnum);
+	rho_water=element->GetMaterialParameter(MaterialsRhoWaterEnum);
+
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetInputListOnVertices(enthalpies,EnthalpyEnum);
+	element->GetInputListOnVertices(pressures,PressureEnum);
+
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	for(iv=0;iv<numvertices;iv++){ 
+		element->EnthalpyToThermal(&temperatures[iv],&waterfractions[iv], enthalpies[iv],pressures[iv]); 
+		deltawaterfractions[iv]=DrainageFunctionWaterfraction(waterfractions[iv], dt);
+	}
+	
+	/*drain waterfraction, feed updated variables back into model*/
+	for(iv=0;iv<numvertices;iv++){
+		if(reCast<bool,IssmDouble>(dt))
+			waterfractions[iv]-=deltawaterfractions[iv]*dt;
+		else
+			waterfractions[iv]-=deltawaterfractions[iv];
+		element->ThermalToEnthalpy(&enthalpies[iv], temperatures[iv], waterfractions[iv], pressures[iv]);
+	}
+	element->AddInput(EnthalpyEnum,enthalpies,P1Enum);
+  	element->AddInput(WaterfractionEnum,waterfractions,P1Enum);
+
+	/*return meltwater column equivalent to drained water*/
+	element->VerticalSegmentIndices(&pairindices,&numsegments);
+	for(is=0;is<numsegments;is++){
+		vertexdown = pairindices[is*2+0];
+		vertexup   = pairindices[is*2+1];
+		height_element=fabs(xyz_list[vertexup*3+2]-xyz_list[vertexdown*3+2]);
+		pdrainrate_element[is]=(deltawaterfractions[vertexdown]+deltawaterfractions[vertexup])/2.*rho_water/rho_ice*height_element;
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(enthalpies);
+	xDelete<IssmDouble>(pressures);
+	xDelete<IssmDouble>(temperatures);
+	xDelete<IssmDouble>(waterfractions);
+	xDelete<IssmDouble>(deltawaterfractions);
+}/*}}}*/
+
+
+
+
+
 void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
 
-}/*}}}*/
+	/* /\*Intermediary*\/ */
+	/* bool        isdynamicbasalspc,setspc; */
+	/* int         numindices, numindicesup; */
+	/* IssmDouble  pressure, pressureup; */
+	/* IssmDouble  h_pmp, enthalpy, enthalpyup; */
+	/* IssmDouble  watercolumn; */
+	/* int        *indices = NULL, *indicesup = NULL; */
+	
+	/* /\* Only update Constraints at the base of grounded ice*\/ */
+	/* if(!IsOnBed() || IsFloating()) return; */
+
+	/* /\*Check wether dynamic basal boundary conditions are activated *\/ */
+	/* element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum); */
+	/* if(!isdynamicbasalspc) return; */
+
+	/* /\*Fetch indices of basal & surface nodes for this finite element*\/ */
+	/* // BasalNodeIndices(&numindices,&indices,this->element_type); */
+	/* // SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type); */
+	/* //	_assert_(numindices==numindicesup); */
+
+	/* /\*Get parameters and inputs: *\/ */
+	/* Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); */
+	/* Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); */
+	/* Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); */
+
+	/* /\*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*\/ */
+	/* GaussPenta* gauss=new GaussPenta(); */
+	/* GaussPenta* gaussup=new GaussPenta(); */
+	/* for(int i=0;i<numindices;i++){ */
+	/* 	gauss->GaussNode(this->element_type,indices[i]); */
+	/* 	gaussup->GaussNode(this->element_type,indicesup[i]); */
+
+	/* 	/\*Check wether there is a temperate layer at the base or not *\/ */
+	/* 	/\*check if node is temperate, if not, continue*\/ */
+	/* 	enthalpy_input->GetInputValue(&enthalpy, gauss); */
+	/* 	pressure_input->GetInputValue(&pressure, gauss); */
+	/* 	watercolumn_input->GetInputValue(&watercolumn,gauss); */
+	/* 	setspc = false; */
+	/* 	// TODO: add case H<Hpmp && watercolumn>0.; */
+	/* 	if (enthalpy>=element->PureIceEnthalpy(pressure)){ */
+	/* 		/\*check if upper node is temperate, too. */
+	/* 			if yes, then we have a temperate layer of positive thickness and reset the spc. */
+	/* 			if not, apply dirichlet BC.*\/ */
+	/* 		enthalpy_input->GetInputValue(&enthalpyup, gaussup); */
+	/* 		pressure_input->GetInputValue(&pressureup, gaussup); */
+	/* 		setspc=((enthalpyup<element->PureIceEnthalpy(pressureup)) && (watercolumn>=0.))?true:false; */
+	/* 	} */
+
+	/* 	if (setspc) { */
+	/* 		/\*Calculate enthalpy at pressure melting point *\/ */
+	/* 		h_pmp=matpar->PureIceEnthalpy(pressure); */
+	/* 		/\*Apply Dirichlet condition (dof = 0 here, since there is only one degree of freedom per node)*\/ */
+	/* 		//element->ApplyConstraint(indices[i],1,h_pmp); */
+	/* 		//nodes[indices[i]]->ApplyConstraint(1,h_pmp); */
+	/* 	} */
+	/* 	else { */
+	/* 		/\*remove spc*\/ */
+	/* 		//element->RemoveConstraint(indices[i],0); */
+	/* 		//nodes[indices[i]]->DofInFSet(0); */
+	/* 	} */
+	/* } */
+
+	/* /\*Free ressources:*\/ */
+	/* xDelete<int>(indices); */
+	/* xDelete<int>(indicesup); */
+	/* delete gauss; */
+	/* delete gaussup; */
+}/*}}}*/
+
+
+
+
 
 /*Intermediaries*/
Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 17013)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 17014)
@@ -40,12 +40,12 @@
 		static void PostProcessing(FemModel* femmodel);
 		static void ComputeBasalMeltingrate(Element* element);
-		static void DrainWaterfraction(Element* element);
+		static void DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element);
 		static void UpdateBasalConstraints(Element* element);
 
 		/*Intermediaries*/
-		IssmDouble EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure);
-		IssmDouble EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum);
-		IssmDouble PureIceEnthalpy(Element* element,IssmDouble pressure);
-		IssmDouble TMeltingPoint(Element* element,IssmDouble pressure);
+		static IssmDouble EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure);
+		static IssmDouble EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum);
+		static IssmDouble PureIceEnthalpy(Element* element,IssmDouble pressure);
+		static IssmDouble TMeltingPoint(Element* element,IssmDouble pressure);
 };
 #endif
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17013)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17014)
@@ -105,4 +105,5 @@
 		virtual void       SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
 		virtual void   ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
+		virtual void   ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure)=0;
 		virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
 		virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0;
@@ -123,4 +124,7 @@
 		virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0;
 
+		virtual Element* GetUpperElement(void)=0;
+		virtual Element* GetLowerElement(void)=0;
+		virtual Element* GetSurfaceElement(void)=0;
 		virtual Element* GetBasalElement(void)=0;
 		virtual void	GetDofList(int** pdoflist,int approximation_enum,int setenum)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17013)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17014)
@@ -134,5 +134,5 @@
 				penta->inputs->AddInput(new PentaInput(input_enum,&extrudedvalues[0],P1Enum));
 				if (penta->IsOnSurface()) break;
-				penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
+				penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
 			}
 		}
@@ -148,5 +148,4 @@
 }
 /*}}}*/
-
 /*FUNCTION Penta::BasalFrictionCreateInput {{{*/
 void Penta::BasalFrictionCreateInput(void){
@@ -479,4 +478,9 @@
 }
 /*}}}*/
+/*FUNCTION Penta::ThermalToEnthalpy{{{*/
+void Penta::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){
+	matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
+}
+/*}}}*/
 /*FUNCTION Penta::EnthalpyToThermal{{{*/
 void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){
@@ -542,11 +546,52 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetBasalElement{{{*/
-Element* Penta::GetBasalElement(void){
+/*FUNCTION Penta::GetUpperPenta{{{*/
+Penta* Penta::GetUpperPenta(void){
+
+	Penta* upper_penta=NULL;
+
+	upper_penta=(Penta*)verticalneighbors[1]; //first one (0) under, second one (1) above
+
+	return upper_penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetLowerPenta{{{*/
+Penta* Penta::GetLowerPenta(void){
+
+	Penta* lower_penta=NULL;
+
+	lower_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
+
+	return lower_penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetSurfacePenta{{{*/
+Penta* Penta::GetSurfacePenta(void){
 
 	/*Output*/
 	Penta* penta=NULL;
 
-	/*Go through all elements till the bed is reached*/
+	/*Go through all pentas till the surface is reached*/
+	penta=this;
+	for(;;){
+		/*Stop if we have reached the surface, else, take upper penta*/
+		if (penta->IsOnSurface()) break;
+
+		/* get upper Penta*/
+		penta=penta->GetUpperPenta();
+		_assert_(penta->Id()!=this->id);
+	}
+
+	/*return output*/
+	return penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetBasalPenta{{{*/
+Penta* Penta::GetBasalPenta(void){
+
+	/*Output*/
+	Penta* penta=NULL;
+
+	/*Go through all pentas till the bed is reached*/
 	penta=this;
 	for(;;){
@@ -555,5 +600,5 @@
 
 		/* get lower Penta*/
-		penta=penta->GetLowerElement();
+		penta=penta->GetLowerPenta();
 		_assert_(penta->Id()!=this->id);
 	}
@@ -561,4 +606,36 @@
 	/*return output*/
 	return penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetUpperElement{{{*/
+Element* Penta::GetUpperElement(void){
+
+	/*Output*/
+	Element* upper_element=this->GetUpperPenta();
+	return upper_element;
+}
+/*}}}*/
+/*FUNCTION Penta::GetLowerElement{{{*/
+Element* Penta::GetLowerElement(void){
+
+	/*Output*/
+	Element* lower_element=this->GetLowerPenta();
+	return lower_element;
+}
+/*}}}*/
+/*FUNCTION Penta::GetSurfaceElement{{{*/
+Element* Penta::GetSurfaceElement(void){
+
+	/*Output*/
+	Element* element=this->GetSurfacePenta();
+	return element;
+}
+/*}}}*/
+/*FUNCTION Penta::GetBasalElement{{{*/
+Element* Penta::GetBasalElement(void){
+
+	/*Output*/
+	Element* element=this->GetBasalPenta();
+	return element;
 }
 /*}}}*/
@@ -833,14 +910,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetLowerElement{{{*/
-Penta* Penta::GetLowerElement(void){
-
-	Penta* upper_penta=NULL;
-
-	upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
-
-	return upper_penta;
-}
-/*}}}*/
 /*FUNCTION Penta::GetNodeIndex {{{*/
 int Penta::GetNodeIndex(Node* node){
@@ -1121,14 +1188,4 @@
 
 	return tau_parameter;
-}
-/*}}}*/
-/*FUNCTION Penta::GetUpperElement{{{*/
-Penta* Penta::GetUpperElement(void){
-
-	Penta* upper_penta=NULL;
-
-	upper_penta=(Penta*)verticalneighbors[1]; //first one under, second one above
-
-	return upper_penta;
 }
 /*}}}*/
@@ -1486,5 +1543,5 @@
 
 		/* get upper Penta*/
-		penta=penta->GetUpperElement();
+		penta=penta->GetUpperPenta();
 		_assert_(penta->Id()!=this->id);
 
@@ -1555,5 +1612,5 @@
 	for(;;){
 		/* get upper Penta*/
-		penta=penta->GetUpperElement();
+		penta=penta->GetUpperPenta();
 		_assert_(penta->Id()!=this->id);
 
@@ -1778,5 +1835,5 @@
 
 		/* get upper Penta*/
-		penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
+		penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
 	}
 
@@ -3416,5 +3473,5 @@
 
 		if(penta->IsOnSurface()) break;
-		penta=penta->GetUpperElement();			
+		penta=penta->GetUpperPenta();			
 	}
 	// add drained water to melting rate
@@ -3460,5 +3517,5 @@
 void Penta::DrainWaterfraction(IssmDouble* drainrate_element){
 
-    /*Intermediaries*/
+  /*Intermediaries*/
 	bool isenthalpy;
 	int iv, index0;
@@ -4609,5 +4666,5 @@
 			
 			/* get upper Penta*/
-			penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
+			penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
 		}
 	}
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17013)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17014)
@@ -72,5 +72,13 @@
 		void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
 		void   Delta18oParameterization(void);
+		void   ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
 		void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
+		Penta* GetUpperPenta(void);
+		Penta* GetLowerPenta(void);
+		Penta* GetSurfacePenta(void);
+		Penta* GetBasalPenta(void);
+		Element* GetUpperElement(void);
+		Element* GetLowerElement(void);
+		Element* GetSurfaceElement(void);
 		Element* GetBasalElement(void);
 		void	 GetDofList(int** pdoflist,int approximation_enum,int setenum);
@@ -211,6 +219,4 @@
 		void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		Node*          GetNode(int node_number);
-		Penta*         GetUpperElement(void);
-		Penta*         GetLowerElement(void);
 		void           InputChangeName(int input_enum, int enum_type_old);
 		void	         InputExtrude(int enum_type,int object_type);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17013)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17014)
@@ -71,8 +71,12 @@
 		void        Delta18oParameterization(void){_error_("not implemented yet");};
 		void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
+		void        ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
 		void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
 		IssmDouble  EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
 		IssmDouble  EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
 		int         FiniteElement(void);
+		Element*    GetUpperElement(void){_error_("not implemented yet");};
+	  	Element*    GetLowerElement(void){_error_("not implemented yet");};
+	  	Element*    GetSurfaceElement(void){_error_("not implemented yet");};
 		Element*    GetBasalElement(void){_error_("not implemented yet");};
 		void        GetDofList(int** pdoflist,int approximation_enum,int setenum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17013)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17014)
@@ -69,6 +69,10 @@
 		void        Delta18oParameterization(void);
 		void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
+		void        ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};
 		void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
 		int         FiniteElement(void);
+		Element*    GetUpperElement(void){_error_("not implemented yet");};
+	  	Element*    GetLowerElement(void){_error_("not implemented yet");};
+	  	Element*    GetSurfaceElement(void){_error_("not implemented yet");};
 		Element*    GetBasalElement(void){_error_("not implemented yet");};
 		void	      GetDofList(int** pdoflist,int approximation_enum,int setenum);