Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 20271)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 20272)
@@ -821,6 +821,7 @@
 	if(!element->IsOnBase() || element->IsFloating()) return NULL;
 
-	bool isdynamicbasalspc;
+	bool converged, isdynamicbasalspc;
 	int i, state;
+	int enthalpy_enum;
 	IssmDouble  dt,Jdet,scalar;
 	IssmDouble	enthalpy, enthalpyup, pressure, pressureup, watercolumn, meltingrate;
@@ -840,8 +841,11 @@
 	element->FindParam(&dt,TimesteppingTimeStepEnum);
 	element->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
+	element->GetInputValue(&converged,ConvergedEnum);
+	if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum; // use enthalpy from last iteration
+	else enthalpy_enum=EnthalpyEnum; // use enthalpy from last time step
 	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);
-	Input* enthalpy_input		 = element->GetInput(EnthalpyPicardEnum);					 _assert_(enthalpy_input);
+	Input* enthalpy_input		 = element->GetInput(enthalpy_enum);					 _assert_(enthalpy_input);
 	Input* pressure_input		 = element->GetInput(PressureEnum);							 _assert_(pressure_input);
 	Input* watercolumn_input	 = element->GetInput(WatercolumnEnum);							 _assert_(watercolumn_input);
@@ -876,6 +880,9 @@
 
 		switch (state) {
-			case 0:
-				// cold, dry base: apply basal surface forcing
+			case 0: case 1: case 2: case 3:
+				// cold, dry base; cold, wet base; refreezing temperate base; thin temperate base: 
+				// Apply basal surface forcing.
+				// Interpolated values of enthalpy on gauss nodes may indicate cold base, 
+				// although one node might have become temperate. So keep heat flux switched on.
 				geothermalflux_input->GetInputValue(&geothermalflux,gauss);
 				friction->GetAlpha2(&alpha2,gauss);
@@ -890,13 +897,4 @@
 					pe->values[i]+=scalar*basis[i];
 				break;
-			case 1:
-				// cold, wet base: keep at pressure melting point 
-				break;
-			case 2:
-				// temperate, thin refreezing base: release spc
-				break;
-			case 3:
-				// temperate, thin melting base: set spc
-				break;
 			case 4:
 				// temperate, thick melting base: set grad H*n=0