Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24219)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24220)
@@ -239,5 +239,5 @@
 	// Take forward differences on left and right edges
 	dT[0] = (T[1] - T[0])/(zGPC[1]-zGPC[0]);
-	dT[m-1] = (T[m-1] - T[m-2])/(zGPC[m-1]-zGPC[m-2]);
+	if(m>1) dT[m-1] = (T[m-1] - T[m-2])/(zGPC[m-1]-zGPC[m-2]);
 
 	//Take centered differences on interior points
@@ -667,6 +667,6 @@
 	KP = xNew<IssmDouble>(m);
 
-	KU[0] = UNDEF;
-	KD[m-1] = UNDEF;
+	KU[0] = Delflag;
+	KD[m-1] = Delflag;
 	for(int i=1;i<m;i++) KU[i]= K[i-1];
 	for(int i=0;i<m-1;i++) KD[i] = K[i+1];
@@ -676,6 +676,6 @@
 	dzU = xNew<IssmDouble>(m);
 	dzD = xNew<IssmDouble>(m);
-	dzU[0]=UNDEF;
-	dzD[m-1]=UNDEF;
+	dzU[0]=Delflag;
+	dzD[m-1]=Delflag;
 
 	for(int i=1;i<m;i++) dzU[i]= dz[i-1];
@@ -1401,5 +1401,5 @@
 		for(int i=0;i<n;i++) m[i] += dW[i];                                            // new mass [kg]
 		for(int i=0;i<n;i++) d[i] = m[i] / dz[i];                                    // density [kg m-3]   
-		for(int i=0;i<n;i++) T[i] = T[i] + (dW[i]*(LF+(CtoK - T[i])*CI)/(m[i]*CI));      // temperature [K]
+		for(int i=0;i<n;i++) if(m[i]>Wtol) T[i] = T[i] + (dW[i]*(LF+(CtoK - T[i])*CI)/(m[i]*CI));      // temperature [K]
 
 		// if pore water froze in ice then adjust d and dz thickness
@@ -1433,5 +1433,5 @@
 
 			int i = 0;
-			while (cellsum(surpE,n) > 0.0+Ttol){
+			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];
@@ -1451,5 +1451,9 @@
 
 		// convert temperature excess to melt [kg]
-		for(int i=0;i<n;i++) M[i] = min(exsT[i] * d[i] * dz[i] * CI / LF, m[i]);  // melt
+		double Mmax=0.0;
+		for(int i=0;i<n;i++){
+			Mmax=exsT[i] * d[i] * dz[i] * CI / LF;
+			M[i] = min(Mmax, m[i]);  // melt
+		}
 		sumM = cellsum(M,n);                                                       // total melt [kg]
 
@@ -1511,5 +1515,5 @@
 			// some or all meltwater refreezes
 			else{
-				// change in density density and temperature
+				// change in density and temperature
 				// _printf_("MELT REFREEZE");
 				//-----------------------melt water-----------------------------
@@ -1541,5 +1545,7 @@
 
 				flxDn[i+1] = inM - F1 - dW[i];     // meltwater out        
-				T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
+				if (m[i]>Wtol){
+					T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
+				}
 
 				// check if an ice layer forms 
@@ -2003,5 +2009,5 @@
 
 		// do not allow densities to exceed the density of ice
-		if(d[i] > dIce-Dtol) d[i]=dIce;
+		if(d[i] > dIce-Ptol) d[i]=dIce;
 
 		// calculate new grid cell length