Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24507)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24508)
@@ -494,5 +494,5 @@
 
 			// modification of albedo due to thin layer of snow or solid
-			// condensation (deposition) at the surface surface
+			// condensation (deposition) at the surface 
 
 			// check if condensation occurs & if it is deposited in solid phase
@@ -1340,4 +1340,6 @@
 	/*outputs:*/
 	IssmDouble  mAdd = 0.0;
+	IssmDouble  surplusE = 0.0;
+	IssmDouble  surplusT = 0.0;
 	IssmDouble dz_add = 0.0;
 	IssmDouble  Rsum = 0.0;
@@ -1433,12 +1435,22 @@
 			int i = 0;
 			while (cellsum(surpE,n) > 0.0+Ttol && i<n){
-				// 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] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
-				T[i+1] = min(CtoK, T[i+1]);
-
-				surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
-				surpE[i+1] = surpT[i+1] * CI * m[i+1];
+
+				if (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] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
+					T[i+1] = min(CtoK, T[i+1]);
+
+					surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
+					surpE[i+1] = surpT[i+1] * CI * m[i+1];
+				}
+				else{
+					surplusT=max(0.0, exsT[i] - LF/CI);
+					surplusE=surpE[i];
+					if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0){
+						_printf0_(" WARNING: surplus energy at the base of GEMB column\n");
+					}
+				}
 
 				// adjust current cell properties (again 159.1342 is the max T)
@@ -1649,5 +1661,5 @@
 			m[i+1] = m_new;                             // combine top masses
 
-			// set cell to 99999 for deletion
+			// set cell to -99999 for deletion
 			m[i] = Delflag;
 		}
@@ -1656,5 +1668,5 @@
 	//If last cell has to be merged
 	if(lastCellFlag){
-         //find closest cell to merge with
+      //find closest cell to merge with
 		for(int i=n-2;i>=0;i--){
 			if(m[i]!=Delflag){
@@ -1679,5 +1691,5 @@
 		m[X1] = m_new;                             // combine top masses
 
-		// set cell to 99999 for deletion
+		// set cell to -99999 for deletion
 		m[X2] = Delflag;
 	}
@@ -1704,5 +1716,5 @@
 	// check if any of the top 10 cell depths are too large
 	X=0;
-	for(int i=9;i>=0;i--){
+	for(int i=min(9,n-1);i>=0;i--){
 		if(dz[i]> 2.0*dzMin+Dtol){
 			X=i;
@@ -1715,19 +1727,19 @@
 		if (dz[j] > dzMin*2.0+Dtol){
 
-				// _printf_("dz > dzMin * 2");
-				// split in two
-				cellsplit(&dz, n, j,.5);
-				cellsplit(&W, n, j,.5);
-				cellsplit(&m, n, j,.5);
-				cellsplit(&T, n, j,1.0);
-				cellsplit(&d, n, j,1.0);
-				cellsplit(&a, n, j,1.0);
-				cellsplit(&EI, n, j,.5);
-				cellsplit(&EW, n, j,.5);
-				cellsplit(&re, n, j,1.0);
-				cellsplit(&gdn, n, j,1.0);
-				cellsplit(&gsp, n, j,1.0);
-				n++;
-				X=X+1;
+			// _printf_("dz > dzMin * 2");
+			// split in two
+			cellsplit(&dz, n, j,.5);
+			cellsplit(&W, n, j,.5);
+			cellsplit(&m, n, j,.5);
+			cellsplit(&T, n, j,1.0);
+			cellsplit(&d, n, j,1.0);
+			cellsplit(&a, n, j,1.0);
+			cellsplit(&EI, n, j,.5);
+			cellsplit(&EW, n, j,.5);
+			cellsplit(&re, n, j,1.0);
+			cellsplit(&gdn, n, j,1.0);
+			cellsplit(&gsp, n, j,1.0);
+			n++;
+			X=X+1;
 		}
 		else j++;
@@ -1817,5 +1829,5 @@
 	#ifndef _HAVE_AD_
 	dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
-	dE = round(sumE0 - sumE1 - sumER +  addE);
+	dE = round(sumE0 - sumE1 - sumER +  addE - surplusE);
 	if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
 			<< "dm: " << dm << " dE: " << dE << "\n");
@@ -1978,6 +1990,12 @@
 				c0arth = 0.07 * H;
 				c1arth = 0.03 * H;
+				//ERA-5
+				//M0 = max(2.3999 - (0.2610 * log(C)),0.25);
+				//M1 = max(2.7469 - (0.3228 * log(C)),0.25);
+				//RACMO
 				M0 = max(1.6599 - (0.1724 * log(C)),0.25);
 				M1 = max(2.0102 - (0.2458 * log(C)),0.25);
+				//From Ligtenberg
+				//H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
 				//M0 = max(1.435 - (0.151 * log(C)),0.25);
 				//M1 = max(2.366 - (0.293 * log(C)),0.25);
@@ -1992,8 +2010,13 @@
 				c0arth = 0.07 * H;
 				c1arth = 0.03 * H;
+				// ERA5
+				//M0 = max(1.8920 - (0.1569 * log(C)),0.25);
+				//M1 = max(2.5662 - (0.2274 * log(C)),0.25);
+				// RACMO
+				M0 = max(1.6201 - (0.1450 * log(C)),0.25);
+				M1 = max(2.5577 - (0.2899 * log(C)),0.25);
+				// From Kuipers Munneke
 				//M0 = max(1.042 - (0.0916 * log(C)),0.25);
 				//M1 = max(1.734 - (0.2039 * log(C)),0.25);
-				M0 = max(1.6201 - (0.1450 * log(C)),0.25);
-				M1 = max(2.5577 - (0.2899 * log(C)),0.25);
 				c0 = M0*c0arth;
 				c1 = M1*c1arth;