Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 22735)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 22736)
@@ -1333,5 +1333,5 @@
 
 	// 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);
 	
@@ -1382,19 +1382,18 @@
             
 			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);  
 		}
 
@@ -1690,5 +1689,5 @@
 		// 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
@@ -1724,5 +1723,5 @@
 		// 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]);
         
@@ -1745,5 +1744,4 @@
 		n=D_size;
 		xDelete<int>(D);
-        
 	}
 
@@ -1763,5 +1761,5 @@
 	/*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"