Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24688)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 24689)
@@ -70,4 +70,5 @@
 				IssmDouble time0     = Ta_input_tr->GetTimeByOffset(-1);
 				IssmDouble timeend   = Ta_input_tr->GetTimeByOffset(offsetend);
+				timeend   = ceil(timeend/dt)*dt;
 				if (time>time0 & timeend>time0){
 					delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
@@ -111,9 +112,9 @@
 	#ifndef _HAVE_AD_  //avoid the round operation check!
 	if (dgpTop != round(dgpTop)){ 
-		_error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop");
+		_error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop\n");
 	}
 	#endif
 	if(dzTop < 0.05){
-		_printf_("initial top grid cell length (dzTop) is < 0.05 m");
+		_printf_("initial top grid cell length (dzTop) is < 0.05 m\n");
 	}
 	gpTop=reCast<int,IssmDouble>(dgpTop);
@@ -282,6 +283,8 @@
 
 	// Take forward differences on left and right edges
-	dT[0] = (T[1] - T[0])/(zGPC[1]-zGPC[0]);
-	if(m>1) dT[m-1] = (T[m-1] - T[m-2])/(zGPC[m-1]-zGPC[m-2]);
+	if(m>1){
+		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]);
+	}
 
 	//Take centered differences on interior points
@@ -343,17 +346,19 @@
 			// higher than 15 degC m-1 (Brun et al., 1992)
 			//When wet-snow grains (class 6) are submitted to a
-			// temperature gradient higher than degC m-1, their sphericity
+			// temperature gradient higher than 5 degC m-1, their sphericity
 			// decreases according to Equations (4). When sphericity
 			// reaches 0, their size increases according to the functions
 			// determined by Marbouty. (Brun et al., 1992)
-			if(gsp[i]>0.0+Gdntol && gsp[i]<1.0-Gdntol && (dT[i]>15.0+Ttol || (T[i]>5.0+Ttol && W[i]>0.0+Wtol)) ){
-				//determine coeff0icients
+			if(gsp[i]>0.0+Gdntol && gsp[i]<1.0-Gdntol && (dT[i]>15.0+Ttol || (dT[i]>5.0+Ttol && W[i]>0.0+Wtol)) ){
+				//determine coefficients
 				IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4);
 				gsp[i] +=C;
 				if (gsp[i]<0.0+Gdntol)gsp[i]=0.0;
+				//determine new grain size (mm)
+				gsz[i] = 0.35 + (0.5-gsp[i])*0.1;
 			}
 			/*Dry snow metamorphism (Marbouty, 1980) grouped model coefficents
 			 *from Marbouty, 1980: Figure 9*/
-			else if(W[i]<=0.0+Wtol || (gsp[i]<=0.0+Gdntol && dT[i]>15.0+Ttol) || (gsp[i]<=0.0+Gdntol && T[i]>5.0+Ttol && W[i]>0.0+Wtol)){
+			else if(W[i]<=0.0+Wtol || (gsp[i]<=0.0+Gdntol && dT[i]>15.0+Ttol) || (gsp[i]<=0.0+Gdntol && dT[i]>5.0+Ttol && W[i]>0.0+Wtol)){
 				//_printf_("Nondendritic snow metamorphism\n");
 				Q = Marbouty(T[i],d[i],dT[i]);
@@ -366,5 +371,5 @@
 				//_printf_("Nondendritic wet snow metamorphism\n");
 				//wet rate of change coefficient
-				IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *dts);   // [mm^3 s^-1]
+				IssmDouble E = (1.28e-8 + 4.22e-10 * pow(lwc[i],3.0))* (dt *dts);   // [mm^3 s^-1]
 
 				// calculate change in grain volume and convert to grain size