Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 22735)
+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 22736)
@@ -1332,7 +1332,7 @@
 	for(int i=0;i<n;i++) exsT[i]= fmax(0.0, T[i] - CtoK);        // [K] to [degC]
 
 	// new grid point center temperature, T [K]
-	//for(int i=0;i<n;i++) T[i]-=exsT[i];
+	// for(int i=0;i<n;i++) T[i]-=exsT[i];
 	for(int i=0;i<n;i++) T[i]=fmin(T[i],CtoK);
 	
 	// specify irreducible water content saturation [fraction]
@@ -1381,21 +1381,20 @@
 			surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI * m[i];
             
 			int i = 0;
-			while (cellsum(surpE,n) > 0.0 + Ttol){
+			while (cellsum(surpE,n) > 0.0 + Ttol && i<n-1){
+
 				// use surplus energy to increase the temperature of lower cell
 				T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
-                
-				exsT[i+1] = fmax(0.0, T[i+1] - CtoK) + exsT[i+1];
-				T[i+1] = fmin(CtoK, T[i+1]);
-                
-				surpT[i+1] = fmax(0.0, exsT[i+1] - LF/CI);
+				surpT[i+1] = fmax(0.0, (T[i+1] - CtoK - LF/CI));
 				surpE[i+1] = surpT[i+1] * CI * m[i+1];
-                
+
 				// adjust current cell properties (again 159.1342 is the max T)
-				exsT[i] = LF/CI;
+				T[i] = T[i]-surpE[i]/m[i+1]/CI;
 				surpE[i] = 0.0;
 				i = i + 1;
 			}
+			// recalculate temperature excess above 0 deg C 
+			for(int i=0;i<n;i++) exsT[i] = fmax(0.0, T[i] - CtoK);  
 		}
 
 		// convert temperature excess to melt [kg]
@@ -1689,7 +1688,7 @@
 		// printf("Total depth < zMin %f \n", Ztot);
 		// mass and energy to be added
 		mAdd = m[n-1]+W[n-1];
-		addE = T[n-1]*m[n-1]*CI;
+		addE = T[n-1]*m[n-1]*CI + W[n-1]*(LF+CtoK*CI);
         
 		// add a grid cell of the same size and temperature to the bottom
 		dz_bot=dz[n-1];
@@ -1723,7 +1722,7 @@
 		// printf("Total depth > zMax %f \n", Ztot);
 		// mass and energy loss
 		mAdd = -(m[n-1]+W[n-1]);
-		addE = -(T[n-1]*m[n-1]*CI);
+		addE = -(T[n-1]*m[n-1]*CI) - (W[n-1]*(LF+CtoK*CI));
 		dz_add=-(dz[n-1]);
         
 		// remove a grid cell from the bottom
@@ -1744,7 +1743,6 @@
 		celldelete(&EW,n,D,D_size);
 		n=D_size;
 		xDelete<int>(D);
-        
 	}
 
 	//// CHECK FOR MASS AND ENERGY CONSERVATION
@@ -1762,7 +1760,7 @@
 
 	/*only in forward mode! avoid round in AD mode as it is not differentiable: */
 	#ifndef _HAVE_ADOLC_
-	dm = round(mSum0 - mSum1 + mAdd);
+	dm = round((mSum0 - mSum1 + mAdd)*1000);
 	dE = round(sumE0 - sumE1 - sumER +  addE);
 	if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
 			<< "dm: " << dm << " dE: " << dE << "\n");