Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 18611)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 18612)
@@ -207,9 +207,8 @@
 	InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum);
 
-	int solution_type;
-	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-	if(solution_type!=SteadystateSolutionEnum){
-		PostProcessing(femmodel);
-	}
+	IssmDouble dt;
+	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+	if(dt==0.) ComputeBasalMeltingrate(femmodel);
+	else PostProcessing(femmodel);
 
 }/*}}}*/
@@ -563,6 +562,9 @@
 	if(!element->IsOnBase() || element->IsFloating()) return NULL;
 
-	IssmDouble  dt,Jdet,enthalpy,pressure,watercolumn,geothermalflux,vx,vy,vz;
-	IssmDouble  enthalpyup,pressureup,alpha2,scalar,basalfriction,heatflux;
+	int i, state;
+	IssmDouble  dt,Jdet,scalar;
+	IssmDouble	enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
+	IssmDouble	vx,vy,vz;
+	IssmDouble  alpha2,basalfriction,geothermalflux,heatflux;
 	IssmDouble *xyz_list_base = NULL;
 
@@ -580,10 +582,10 @@
 	Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
 	Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
-	Input* enthalpy_input       = element->GetInput(EnthalpyPicardEnum);              _assert_(enthalpy_input);
-	Input* pressure_input       = element->GetInput(PressureEnum);                    _assert_(pressure_input);
+	Input* enthalpy_input		 = element->GetInput(EnthalpyPicardEnum);					 _assert_(enthalpy_input);
+	Input* pressure_input		 = element->GetInput(PressureEnum);							 _assert_(pressure_input);
+	Input* watercolumn_input	 = element->GetInput(WatercolumnEnum);							 _assert_(watercolumn_input);
+	Input* meltingrate_input	 = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);							 _assert_(meltingrate_input);
 	Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
-	Input* watercolumn_input    = element->GetInput(WatercolumnEnum);                 _assert_(watercolumn_input);
-	IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
-	IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+	IssmDouble  rho_ice			 = element->GetMaterialParameter(MaterialsRhoIceEnum);
 
 	/*Build friction element, needed later: */
@@ -591,6 +593,6 @@
 
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss   = element->NewGaussBase(2);
-	Gauss* gaussup = element->NewGaussTop(2);
+	GaussPenta* gauss=new GaussPenta(2,2);
+	GaussPenta* gaussup=new GaussPenta(2,2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
@@ -600,37 +602,39 @@
 
 		enthalpy_input->GetInputValue(&enthalpy,gauss);
+		enthalpy_input->GetInputValue(&enthalpyup,gaussup);
 		pressure_input->GetInputValue(&pressure,gauss);
+		pressure_input->GetInputValue(&pressureup,gaussup);
 		watercolumn_input->GetInputValue(&watercolumn,gauss);
-
-		if((dt==0.) || ((watercolumn<=0.) && (enthalpy<PureIceEnthalpy(element,pressure)))){
-			/* the above check is equivalent to 
-			 NOT [(watercolumn>0.) AND (enthalpy<PIE)] AND (enthalpy<PIE)*/
-			geothermalflux_input->GetInputValue(&geothermalflux,gauss);
-
-			friction->GetAlpha2(&alpha2,gauss);
-			vx_input->GetInputValue(&vx,gauss);
-			vy_input->GetInputValue(&vy,gauss);
-			vz_input->GetInputValue(&vz,gauss);
-			basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
-			heatflux      = (basalfriction+geothermalflux)/(rho_ice);
-
-			scalar = gauss->weight*Jdet*heatflux;
-			if(dt!=0.) scalar=dt*scalar;
-
-			for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
-		}
-		else if(enthalpy >= PureIceEnthalpy(element,pressure)){
-			/* check positive thickness of temperate basal ice layer */
-			enthalpy_input->GetInputValue(&enthalpyup,gaussup);
-			pressure_input->GetInputValue(&pressureup,gaussup);
-			if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
-				// do nothing, set grad enthalpy*n=0.
-			}
-			else{
-				// only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
-			}
-		}
-		else{
-			// base cold, but watercolumn positive. Set base to pressure melting point enthalpy
+		meltingrate_input->GetInputValue(&meltingrate,gauss);
+
+		state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
+		switch (state) {
+			case 0:
+				// cold, dry base: apply basal surface forcing
+				geothermalflux_input->GetInputValue(&geothermalflux,gauss);
+				friction->GetAlpha2(&alpha2,gauss);
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				vz_input->GetInputValue(&vz,gauss);
+				basalfriction=alpha2*(vx*vx+vy*vy+vz*vz);
+				heatflux=(basalfriction+geothermalflux)/(rho_ice);
+				scalar=gauss->weight*Jdet*heatflux;
+				if(dt!=0.) scalar=dt*scalar;
+				for(i=0;i<numnodes;i++) 
+					pe->values[i]+=scalar*basis[i];
+				break;
+			case 1:
+				// cold, wet base: keep at pressure melting point 
+			case 2:
+				// temperate, thin refreezing base: release spc
+			case 3:
+				// temperate, thin melting base: set spc
+			case 4:
+				// temperate, thick melting base: set grad H*n=0
+				for(i=0;i<numnodes;i++) 
+					pe->values[i]+=0.;
+				break;
+			default:
+				_printf0_("	unknown thermal basal state found!");
 		}
 	}
@@ -653,5 +657,5 @@
 	if(!element->IsOnBase() || !element->IsFloating()) return NULL;
 
-	IssmDouble  h_pmp,dt,Jdet,scalar_ocean,pressure;
+	IssmDouble  Hpmp,dt,Jdet,scalar_ocean,pressure;
 	IssmDouble *xyz_list_base = NULL;
 
@@ -683,7 +687,7 @@
 
 		pressure_input->GetInputValue(&pressure,gauss);
-		h_pmp=element->PureIceEnthalpy(pressure);
-
-		scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*h_pmp/(heatcapacity*rho_ice);
+		Hpmp=element->PureIceEnthalpy(pressure);
+
+		scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*Hpmp/(heatcapacity*rho_ice);
 		if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
 
@@ -887,34 +891,23 @@
 
 	/*Intermediaries*/
-	int solution_type, i;
 	bool computebasalmeltingrates=true;
-	bool isdrainage=true;
-	bool updatebasalconstraints=true;
-
-	if(isdrainage){
-		/*Drain excess water fraction in ice column: */
-		for(i=0;i<femmodel->elements->Size();i++){
-			Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			DrainWaterfractionIcecolumn(element);
-		}
-	}
-
-	if(computebasalmeltingrates){
-		/*Compute basal melting rates: */
-		for(i=0;i<femmodel->elements->Size();i++){
-			Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			ComputeBasalMeltingrate(element);
-		}
-	}
-
-	if(updatebasalconstraints){
-		/*Update basal dirichlet BCs for enthalpy in transient runs: */
-		femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-		if(solution_type==TransientSolutionEnum){
-			for(i=0;i<femmodel->elements->Size();i++){
-				Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-				UpdateBasalConstraints(element);
-			}
-		}
+	bool drainicecolumn=true;
+	bool isdynamicbasalspc;
+	IssmDouble dt;
+
+	femmodel->parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+	//TODO: use dt to decide what to do
+	if(drainicecolumn)	DrainWaterfraction(femmodel);
+	if(computebasalmeltingrates)	ComputeBasalMeltingrate(femmodel);
+	if(isdynamicbasalspc)	UpdateBasalConstraints(femmodel);
+
+}/*}}}*/
+void EnthalpyAnalysis::ComputeBasalMeltingrate(FemModel* femmodel){/*{{{*/
+	/*Compute basal melting rates: */
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ComputeBasalMeltingrate(element);
 	}
 }/*}}}*/
@@ -931,11 +924,11 @@
 	/* Intermediaries */
 	const int   dim=3;
-	int         i,is,vertexdown,vertexup,numvertices,numsegments;
-	IssmDouble  heatflux;
-	IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim];
-	IssmDouble  basalfriction,alpha2;
+	int         i,is,state;
+	int			vertexdown,vertexup,numvertices,numsegments;
+	const int	enthalpy_enum=EnthalpyEnum;
+	IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim],d1pressure[dim];
+	IssmDouble  basalfriction,alpha2,geothermalflux,heatflux;
 	IssmDouble  dt,yts;
 	IssmDouble  melting_overshoot,lambda;
-	IssmDouble  geothermalflux;
 	IssmDouble  vx,vy,vz;
 	IssmDouble *xyz_list      = NULL;
@@ -943,5 +936,5 @@
 	int        *pairindices   = NULL;
 
-	/*Fetch parameters and inputs */
+	/*Fetch parameters*/
 	element->GetVerticesCoordinates(&xyz_list);
 	element->GetVerticesCoordinatesBase(&xyz_list_base);
@@ -949,10 +942,20 @@
 	IssmDouble rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
 	IssmDouble rho_water  = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
-	Input* enthalpy_input         = element->GetInput(EnthalpyEnum);                    _assert_(enthalpy_input);
+	IssmDouble beta		 = element->GetMaterialParameter(MaterialsBetaEnum);
+	IssmDouble kappa		 = EnthalpyDiffusionParameterVolume(element,EnthalpyEnum);     _assert_(kappa>=0.);
+	IssmDouble kappa_mix;
+
+	/*retrieve inputs*/
+	Input* enthalpy_input         = element->GetInput(enthalpy_enum);                    _assert_(enthalpy_input);
+	Input* pressure_input			= element->GetInput(PressureEnum);							 _assert_(pressure_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.);
+
+	/*Build friction element, needed later: */
+	Friction* friction=new Friction(element,dim);
+
+	/******** MELTING RATES  ************************************//*{{{*/
 	element->NormalBase(&normal_base[0],xyz_list_base);
 	element->VerticalSegmentIndices(&pairindices,&numsegments);
@@ -960,51 +963,48 @@
 	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,BasalforcingsGroundediceMeltingRateEnum);
+	IssmDouble* enthalpies = xNew<IssmDouble>(numvertices);
+	IssmDouble* pressures = xNew<IssmDouble>(numvertices);
+	IssmDouble* watercolumns = xNew<IssmDouble>(numvertices);
+	IssmDouble* basalmeltingrates = xNew<IssmDouble>(numvertices);
+	element->GetInputListOnVertices(enthalpies,enthalpy_enum);
+	element->GetInputListOnVertices(pressures,PressureEnum);
+	element->GetInputListOnVertices(watercolumns,WatercolumnEnum);
+	element->GetInputListOnVertices(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum);
 
 	Gauss* gauss=element->NewGauss();
-	
-	for(int is=0;is<numsegments;is++){
+	for(is=0;is<numsegments;is++){
 		vertexdown = pairindices[is*2+0];
 		vertexup   = pairindices[is*2+1];
 		gauss->GaussVertex(vertexdown);
-		
-		bool checkpositivethickness=true;
-		_assert_(watercolumn[vertexdown]>=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{
+
+		state=GetThermalBasalCondition(element, enthalpies[vertexdown], enthalpies[vertexup], pressures[vertexdown], pressures[vertexup], watercolumns[vertexdown], basalmeltingrates[vertexdown]);
+		switch (state) {
+			case 0:
+				// cold, dry base: apply basal surface forcing
+				for(i=0;i<3;i++) vec_heatflux[i]=0.;
+				break;
+			case 1:
+				// cold, wet base: keep at pressure melting point 
+			case 2:
+				// temperate, thin refreezing base: release spc
+
+			case 3:
+				// temperate, thin melting base: set spc
+				// enthalpies[vertexdown]=element->PureIceEnthalpy(pressures[vertexdown]); 
 				enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],xyz_list,gauss);
 				for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
-			}
-
+				break;
+			case 4:
+				// temperate, thick melting base: set grad H*n=0
+				kappa_mix=GetWetIceConductivity(element, enthalpies[vertexdown], pressures[vertexdown]);
+				pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
+				for(i=0;i<3;i++) vec_heatflux[i]=kappa_mix*beta*d1pressure[i];
+				break;
+			default:
+				_printf0_("	unknown thermal basal state found!");
+		}
+		if(state==0) meltingrate_enthalpy[is]=0.;
+		else{
 			/*heat flux along normal*/
 			heatflux=0.;
@@ -1013,9 +1013,6 @@
 			/*basal friction*/
 			friction->GetAlpha2(&alpha2,gauss);
-			vx_input->GetInputValue(&vx,gauss);
-			vy_input->GetInputValue(&vy,gauss);
-			vz_input->GetInputValue(&vz,gauss);
+			vx_input->GetInputValue(&vx,gauss);		vy_input->GetInputValue(&vy,gauss);		vz_input->GetInputValue(&vz,gauss);
 			basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
-
 			geothermalflux_input->GetInputValue(&geothermalflux,gauss);
 			/* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/
@@ -1023,37 +1020,39 @@
 			meltingrate_enthalpy[is]=heating[is]/(latentheat*rho_ice); // m/s water equivalent
 		}
-	}
-	/******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
+	}/*}}}*/
+
+	/******** UPDATE MELTINGRATES AND WATERCOLUMN **************//*{{{*/
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	element->FindParam(&yts, ConstantsYtsEnum);
 	for(is=0;is<numsegments;is++){
 		vertexdown = pairindices[is*2+0];
 		vertexup   = pairindices[is*2+1];
 		if(dt!=0.){
-			if(watercolumn[vertexdown]+meltingrate_enthalpy[is]*dt<0.){	// prevent too much freeze on			
-				lambda = -watercolumn[vertexdown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
-				watercolumn[vertexdown]=0.;
-				basalmeltingrate[vertexdown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
-				enthalpy[vertexdown]+=(1.-lambda)*meltingrate_enthalpy[is]*dt*latentheat; // use rest of energy to cool down base
+			if(watercolumns[vertexdown]+meltingrate_enthalpy[is]*dt<0.){	// prevent too much freeze on			
+				lambda = -watercolumns[vertexdown]/(dt*meltingrate_enthalpy[is]); _assert_(lambda>=0.); _assert_(lambda<1.);
+				watercolumns[vertexdown]=0.;
+				basalmeltingrates[vertexdown]=lambda*meltingrate_enthalpy[is]; // restrict freeze on only to size of watercolumn
+				enthalpies[vertexdown]+=(1.-lambda)*dt/yts*meltingrate_enthalpy[is]*latentheat*rho_ice; // use rest of energy to cool down base: dE=L*m, m=(1-lambda)*meltingrate*rho_ice
 			}
 			else{
-				basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
-				watercolumn[vertexdown]+=dt*meltingrate_enthalpy[is]; 
+				basalmeltingrates[vertexdown]=meltingrate_enthalpy[is];
+				watercolumns[vertexdown]+=dt*meltingrate_enthalpy[is]; 
 			}
 		}
 		else{
-			basalmeltingrate[vertexdown]=meltingrate_enthalpy[is];
-			if(watercolumn[vertexdown]+meltingrate_enthalpy[is]<0.)
-				watercolumn[vertexdown]=0.;
+			basalmeltingrates[vertexdown]=meltingrate_enthalpy[is];
+			if(watercolumns[vertexdown]+meltingrate_enthalpy[is]<0.)
+				watercolumns[vertexdown]=0.;
 			else
-				watercolumn[vertexdown]+=meltingrate_enthalpy[is];
+				watercolumns[vertexdown]+=meltingrate_enthalpy[is];
 		}	
-		basalmeltingrate[vertexdown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
-		_assert_(watercolumn[vertexdown]>=0.);
-	}
+		basalmeltingrates[vertexdown]*=rho_water/rho_ice; // convert meltingrate from water to ice equivalent
+		_assert_(watercolumns[vertexdown]>=0.);
+	}/*}}}*/
 
 	/*feed updated variables back into model*/
-	element->AddInput(EnthalpyEnum,enthalpy,P1Enum);
-	element->AddInput(WatercolumnEnum,watercolumn,P1Enum);
-	element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrate,P1Enum);
+	element->AddInput(EnthalpyEnum,enthalpies,P1Enum); //TODO: distinguis for steadystate and transient run
+	element->AddInput(WatercolumnEnum,watercolumns,P1Enum);
+	element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,P1Enum);
 
 	/*Clean up and return*/
@@ -1061,12 +1060,19 @@
 	delete friction;
 	xDelete<int>(pairindices);
-	xDelete<IssmDouble>(enthalpy);
-	xDelete<IssmDouble>(pressure);
-	xDelete<IssmDouble>(watercolumn);
-	xDelete<IssmDouble>(basalmeltingrate);
+	xDelete<IssmDouble>(enthalpies);
+	xDelete<IssmDouble>(pressures);
+	xDelete<IssmDouble>(watercolumns);
+	xDelete<IssmDouble>(basalmeltingrates);
 	xDelete<IssmDouble>(meltingrate_enthalpy);
 	xDelete<IssmDouble>(heating);
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(xyz_list_base);
+}/*}}}*/
+void EnthalpyAnalysis::DrainWaterfraction(FemModel* femmodel){/*{{{*/
+	/*Drain excess water fraction in ice column: */
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		DrainWaterfractionIcecolumn(element);
+	}
 }/*}}}*/
 void EnthalpyAnalysis::DrainWaterfractionIcecolumn(Element* element){/*{{{*/
@@ -1173,4 +1179,11 @@
 	xDelete<IssmDouble>(deltawaterfractions);
 }/*}}}*/
+void EnthalpyAnalysis::UpdateBasalConstraints(FemModel* femmodel){/*{{{*/
+	/*Update basal dirichlet BCs for enthalpy: */
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		UpdateBasalConstraints(element);
+	}
+}/*}}}*/
 void EnthalpyAnalysis::UpdateBasalConstraints(Element* element){/*{{{*/
 
@@ -1182,11 +1195,6 @@
 
 	/*Intermediary*/
-	bool        isdynamicbasalspc,setspc;
-	int         numindices, numindicesup;
-	IssmDouble  pressure, pressureup;
-	IssmDouble  h_pmp, enthalpy, enthalpyup;
-	IssmDouble  watercolumn;
-	int        *indices = NULL, *indicesup = NULL;
-	Node*       node = NULL;
+	bool        isdynamicbasalspc;
+	IssmDouble	dt;
 
 	/*Check wether dynamic basal boundary conditions are activated */
@@ -1194,16 +1202,127 @@
 	if(!isdynamicbasalspc) return;
 
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	if(dt==0.){
+		UpdateBasalConstraintsSteadystate(element);
+	}
+	else{
+		UpdateBasalConstraintsTransient(element);
+	}
+}/*}}}*/
+void EnthalpyAnalysis::UpdateBasalConstraintsTransient(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return;
+
+	/* Only update Constraints at the base of grounded ice*/
+	if(!(element->IsOnBase()) || element->IsFloating()) return;
+
+	/*Intermediary*/
+	bool        isdynamicbasalspc,setspc;
+	int         numindices, numindicesup, state;
+	int        *indices = NULL, *indicesup = NULL;
+	Node*       node = NULL;
+	IssmDouble	enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
+
+	/*Check wether dynamic basal boundary conditions are activated */
+	element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+	if(!isdynamicbasalspc) return;
+
+	/*Get parameters and inputs: */
+	Input* enthalpy_input       = element->GetInput(EnthalpyEnum);                    _assert_(enthalpy_input); //TODO: check EnthalpyPicard?
+	Input* pressure_input		 = element->GetInput(PressureEnum);							 _assert_(pressure_input);
+	Input* watercolumn_input	 = element->GetInput(WatercolumnEnum);							 _assert_(watercolumn_input);
+	Input* meltingrate_input	 = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);							 _assert_(meltingrate_input);
+
 	/*Fetch indices of basal & surface nodes for this finite element*/
 	Penta *penta =  (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
 	penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
-	penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType());
-	_assert_(numindices==numindicesup);
+	penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType());	_assert_(numindices==numindicesup);
+
+	GaussPenta* gauss=new GaussPenta();
+	GaussPenta* gaussup=new GaussPenta();
+
+	for(int i=0;i<numindices;i++){
+		gauss->GaussNode(element->GetElementType(),indices[i]);
+		gaussup->GaussNode(element->GetElementType(),indicesup[i]);
+		
+		enthalpy_input->GetInputValue(&enthalpy,gauss);
+		enthalpy_input->GetInputValue(&enthalpyup,gaussup);
+		pressure_input->GetInputValue(&pressure,gauss);
+		pressure_input->GetInputValue(&pressureup,gaussup);
+		watercolumn_input->GetInputValue(&watercolumn,gauss);
+		meltingrate_input->GetInputValue(&meltingrate,gauss);
+
+		state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
+
+		setspc=false;
+		switch (state) {
+			case 0:
+				// cold, dry base: apply basal surface forcing
+				break;
+			case 1:
+				// cold, wet base: keep at pressure melting point 
+				setspc=true;
+				break;
+			case 2:
+				// temperate, thin refreezing base: release spc
+				break;
+			case 3:
+				// temperate, thin melting base: set spc
+				setspc=true;
+				break;
+			case 4:
+				// temperate, thick melting base: set grad H*n=0
+				break;
+			default:
+				_printf0_("	unknown thermal basal state found!");
+		}
+
+		/*apply or release spc*/
+		node=element->GetNode(indices[i]);
+		if(setspc){
+			pressure_input->GetInputValue(&pressure, gauss);
+			node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
+		}
+		else			
+			node->DofInFSet(0);
+	}
+
+	/*Free ressources:*/
+	xDelete<int>(indices);
+	xDelete<int>(indicesup);
+	delete gauss;
+	delete gaussup;
+}/*}}}*/
+void EnthalpyAnalysis::UpdateBasalConstraintsSteadystate(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return;
+
+	/* Only update Constraints at the base of grounded ice*/
+	if(!(element->IsOnBase()) || element->IsFloating()) return;
+
+	/*Intermediary*/
+	bool        isdynamicbasalspc,setspc;
+	int         numindices, numindicesup, state;
+	int        *indices = NULL, *indicesup = NULL;
+	Node*       node = NULL;
+	IssmDouble	enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
+
+	/*Check wether dynamic basal boundary conditions are activated */
+	element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+	if(!isdynamicbasalspc) return;
 
 	/*Get parameters and inputs: */
-	Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
-	Input* enthalpy_input=element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
-	Input* watercolumn_input=element->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
-
-	/*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/
+	Input* enthalpy_input		 = element->GetInput(EnthalpyPicardEnum);					 _assert_(enthalpy_input);
+	Input* pressure_input		 = element->GetInput(PressureEnum);							 _assert_(pressure_input);
+	Input* watercolumn_input	 = element->GetInput(WatercolumnEnum);							 _assert_(watercolumn_input);
+	Input* meltingrate_input	 = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);							 _assert_(meltingrate_input);
+
+	/*Fetch indices of basal & surface nodes for this finite element*/
+	Penta *penta =  (Penta *) element; // TODO: add Basal-/SurfaceNodeIndices to element.h, and change this to Element*
+	penta->BasalNodeIndices(&numindices,&indices,element->GetElementType());
+	penta->SurfaceNodeIndices(&numindicesup,&indicesup,element->GetElementType());	_assert_(numindices==numindicesup);
+
 	GaussPenta* gauss=new GaussPenta();
 	GaussPenta* gaussup=new GaussPenta();
@@ -1212,26 +1331,44 @@
 		gaussup->GaussNode(element->GetElementType(),indicesup[i]);
 
-		/*Check wether there is a temperate layer at the base or not */
-		/*check if node is temperate, else continue*/
-		enthalpy_input->GetInputValue(&enthalpy, gauss);
-		pressure_input->GetInputValue(&pressure, gauss);
+		enthalpy_input->GetInputValue(&enthalpy,gauss);
+		enthalpy_input->GetInputValue(&enthalpyup,gaussup);
+		pressure_input->GetInputValue(&pressure,gauss);
+		pressure_input->GetInputValue(&pressureup,gaussup);
 		watercolumn_input->GetInputValue(&watercolumn,gauss);
-		h_pmp=PureIceEnthalpy(element,pressure);
-		if (enthalpy>=h_pmp){
-			/*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<PureIceEnthalpy(element,pressureup)) && (watercolumn>=0.))?true:false;
-		}
-		else
-			setspc = false;
-
+		meltingrate_input->GetInputValue(&meltingrate,gauss);
+
+		state=GetThermalBasalCondition(element, enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate);
+		setspc=false;
+		switch (state) {
+			case 0:
+				// cold, dry base: apply basal surface forcing
+				break;
+			case 1:
+				// cold, wet base: keep at pressure melting point 
+				setspc=true;
+				break;
+			case 2:
+				// temperate, thin refreezing base: release spc
+				break;
+			case 3:
+				// temperate, thin melting base: set spc
+				setspc=true;
+				break;
+			case 4:
+				// temperate, thick melting base: s
+				setspc=true;
+				break;
+			default:
+				_printf0_("	unknown thermal basal state found!");
+		}
+
+		/*apply or release spc*/
 		node=element->GetNode(indices[i]);
-		if(setspc) 
-			node->ApplyConstraint(0,h_pmp); /*apply spc*/ 
+		if(setspc){
+			pressure_input->GetInputValue(&pressure, gauss);
+			node->ApplyConstraint(0,PureIceEnthalpy(element,pressure));
+		}
 		else			
-			node->DofInFSet(0); /*remove spc*/ 
+			node->DofInFSet(0);
 	}
 
@@ -1241,4 +1378,56 @@
 	delete gauss;
 	delete gaussup;
+}/*}}}*/
+int EnthalpyAnalysis::GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpyup, IssmDouble pressure, IssmDouble pressureup, IssmDouble watercolumn, IssmDouble meltingrate){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return -1;
+
+	/* Only update Constraints at the base of grounded ice*/
+	if(!(element->IsOnBase())) return -1;
+
+	/*Intermediary*/
+	int state=-1;
+	IssmDouble	dt;
+
+	/*Get parameters and inputs: */
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+
+	if(dt==0.){ // steadystate case
+		state=0; //TODO: add consistent steadystate basal condition scheme
+// 		if(enthalpy<PureIceEnthalpy(element,pressure)){ /*is base cold?*/
+// 			if(watercolumn<=0.) state=0; // cold, dry base
+// 			else state=1; // cold, wet base (refreezing)
+// 		}
+// 		else{ /*base is temperate, check if upper node is temperate, too.*/
+// 			if(enthalpyup<PureIceEnthalpy(element,pressureup))
+// 				if(meltingrate<0.)	state=2; // refreezing temperate base (non-physical, only for steadystate solver)
+// 				else						state=3; // melting temperate base with no temperate layer
+// 			else
+// 				state=4; // melting temperate base with temperate layer of positive thickness
+// 		}
+	}
+	else{ // transient case
+		if(enthalpy<PureIceEnthalpy(element,pressure)){
+			if(watercolumn<=0.) state=0; // cold, dry base
+			else state=1; // cold, wet base (refreezing)
+		}
+		else{
+			if(enthalpyup<PureIceEnthalpy(element,pressureup))	state=3; // temperate base, but no temperate layer
+			else state=4; // temperate layer with positive thickness
+		}
+	}
+
+	_assert_(state>=0);
+	return state;
+}/*}}}*/
+IssmDouble EnthalpyAnalysis::GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure){/*{{{*/
+
+	IssmDouble temperature, waterfraction;
+	IssmDouble kappa_w = 0.6; // thermal conductivity of water (in W/m/K)
+	IssmDouble kappa_i = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
+	element->EnthalpyToThermal(&temperature, &waterfraction, enthalpy, pressure);
+
+	return (1.-waterfraction)*kappa_i + waterfraction*kappa_w;
 }/*}}}*/
 
Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 18611)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 18612)
@@ -8,4 +8,5 @@
 /*Headers*/
 #include "./Analysis.h"
+#include "../classes/classes.h"
 
 class EnthalpyAnalysis: public Analysis{
@@ -41,8 +42,18 @@
 		/*Modules*/
 		static void PostProcessing(FemModel* femmodel);
+		static void ComputeBasalMeltingrate(FemModel* femmodel);
 		static void ComputeBasalMeltingrate(Element* element);
+		static void DrainWaterfraction(FemModel* femmodel);
 		static void DrainWaterfractionIcecolumn(Element* element);
 		static void DrainWaterfraction(Element* element, IssmDouble* pdrainrate_element);
+		static void UpdateBasalConstraints(FemModel* femmodel);
 		static void UpdateBasalConstraints(Element* element);
+		static void UpdateBasalConstraintsTransient(Element* element);
+		static void UpdateBasalConstraintsSteadystate(Element* element);
+// 		static int GetThermalBasalCondition(Element* element, Gauss* gauss, Gauss* gaussup, int enthalpy_enum);
+// 		static int GetThermalBasalCondition(Element* element, GaussPenta* gauss, GaussPenta* gaussup, int enthalpy_enum);
+		static int GetThermalBasalCondition(Element* element, IssmDouble enthalpy, IssmDouble enthalpy_up, IssmDouble pressure, IssmDouble pressure_up, IssmDouble watercolumn, IssmDouble meltingrate);
+		static IssmDouble GetWetIceConductivity(Element* element, IssmDouble enthalpy, IssmDouble pressure);
+
 
 		/*Intermediaries*/