Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24507)
+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24508)
@@ -493,7 +493,7 @@
 			for(int i=0;i<m;i++)a[i] -= d_a[i];                            // new albedo
 
 			// 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
 			if ( EC > 0.0 + Dtol && T[0] < 0.0-Ttol) P = P + (EC/dSnow) * 1000.0;  // add cond to precip [mm]
@@ -1339,6 +1339,8 @@
 
 	/*outputs:*/
 	IssmDouble  mAdd = 0.0;
+	IssmDouble  surplusE = 0.0;
+	IssmDouble  surplusT = 0.0;
 	IssmDouble dz_add = 0.0;
 	IssmDouble  Rsum = 0.0;
 	IssmDouble* T=*pT;
@@ -1432,15 +1434,25 @@
 
 			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]);
+				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];
 
-				surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
-				surpE[i+1] = surpT[i+1] * CI * m[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)
 				exsT[i] = LF/CI;
 				surpE[i] = 0.0;
@@ -1648,7 +1660,7 @@
 			W[i+1] = W[i+1] + W[i];                     // combine liquid water
 			m[i+1] = m_new;                             // combine top masses
 
-			// set cell to 99999 for deletion
+			// set cell to -99999 for deletion
 			m[i] = Delflag;
 		}
 	}
@@ -1655,7 +1667,7 @@
 
 	//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){
 				X2=i;
@@ -1678,7 +1690,7 @@
  		W[X1] = W[X1] + W[X2];                     // combine liquid water
 		m[X1] = m_new;                             // combine top masses
 
-		// set cell to 99999 for deletion
+		// set cell to -99999 for deletion
 		m[X2] = Delflag;
 	}
 
@@ -1703,7 +1715,7 @@
 
 	// 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;
 			break;
@@ -1714,21 +1726,21 @@
 	while(j<=X){
 		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++;
 	}
@@ -1816,7 +1828,7 @@
 	/*only in forward mode! avoid round in AD mode as it is not differentiable: */
 	#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");
 	#endif
@@ -1977,8 +1989,14 @@
 				H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
 				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);
 				c0 = M0*c0arth;
@@ -1991,10 +2009,15 @@
 				H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
 				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;
 				break;