Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 19581)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 19582)
@@ -7,5 +7,5 @@
 #include "../../toolkits/toolkits.h"
 
-const double Pi =3.141592653589793238462643383279502884197169399375105820974944592308;
+const double Pi = 3.141592653589793;
 
 void Gembx(FemModel* femmodel){  /*{{{*/
@@ -100,18 +100,18 @@
 
 	// initialize
-	IssmDouble F = 0, H=0, G=0, E;
-
+	IssmDouble F = 0, H=0, G=0;
+	const IssmDouble E = 0.09;        //[mm d-1] model time growth constant E
 	// convert T from K to ºC
 	T = T - 273.15;
 
 	// temperature coefficient F
-	if(T>-6.0) F =  0.7 + ((T/-6) * 0.3);
-	if(T<=-6.0 && T>-22.0) F =  1 - ((T+6)/-16 * 0.8);
-	if(T<=-22.0 && T>=-40.0) F =  0.2 - ((T+22)/-18 * 0.2);
+	if(T>-6.0) F =  0.7 + ((T/-6.0) * 0.3);
+	if(T<=-6.0 && T>-22.0) F =  1 - ((T+6.0)/-16.0 * 0.8);
+	if(T<=-22.0 && T>-40.0) F =  0.2 - ((T+22.0)/-18.0 * 0.2);
 
 	// density coefficient F
-	if(d<150) H=1.0;
-
-	if(d>=150 & d<400) H = 1 - ((d-150)/250);
+	if(d<150.0) H=1.0;
+
+	if(d>=150.0 & d<400.0) H = 1 - ((d-150.0)/250.0);
 
 	// temperature gradient coefficient G
@@ -121,7 +121,4 @@
 	if(dT >= 0.5  && dT < 0.7)  G = 0.9 + (((dT - 0.5)/0.2) * 0.1);
 	if(dT >= .7              )  G = 1.0;
-
-	// model time growth constant E
-	E = 0.09;        //[mm d-1]
 
 	// grouped coefficient Q
@@ -182,8 +179,8 @@
 
 	/*Figure out grain size from effective grain radius: */
-	gsz=xNew<IssmDouble>(m); for(int i=0;i<m;i++)gsz[i]=re[i]*2;
+	gsz=xNew<IssmDouble>(m); for(int i=0;i<m;i++)gsz[i]=re[i]*2.0;
 
 	/*Convert dt from seconds to day: */
-	dt=smb_dt/86400;
+	dt=smb_dt/86400.0;
 
 	/*Determine liquid-water content in percentage: */
@@ -191,5 +188,5 @@
 
 	//set maximum water content by mass to 9 percent (Brun, 1980)
-	for(int i=0;i<m;i++)if(lwc[i]>9) lwc[i]=9.0;
+	for(int i=0;i<m;i++)if(lwc[i]>9.0) lwc[i]=9.0;
 
 
@@ -224,5 +221,5 @@
 
 			//index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1
-			if(dT[i]<5){
+			if(dT[i]<5.0){
 				//determine coefficients
 				IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt;
@@ -242,10 +239,10 @@
 
 			// wet snow metamorphism
-			if(W[i]>0){
+			if(W[i]>0.0){
 
 				//_printf_("D}ritic wet snow metamorphism\n");
 
 				//determine coefficient
-				IssmDouble D = (1/16) * pow(lwc[i],3.0) * dt;
+				IssmDouble D = (1.0/16.0) * pow(lwc[i],3.0) * dt;
 
 				// new dendricity and sphericity for wet snow
@@ -255,8 +252,8 @@
 
 			// dendricity and sphericity can not be > 1 or < 0
-			if (gdn[i]<0)gdn[i]=0.0;
-			if (gsp[i]<0)gsp[i]=0.0;
-			if (gdn[i]>1)gdn[i]=1.0;
-			if (gsp[i]>1)gsp[i]=1.0;
+			if (gdn[i]<0.0)gdn[i]=0.0;
+			if (gsp[i]<0.0)gsp[i]=0.0;
+			if (gdn[i]>1.0)gdn[i]=1.0;
+			if (gsp[i]>1.0)gsp[i]=1.0;
 
 			// determine new grain size (mm)
@@ -276,16 +273,16 @@
 
 			//Wet snow metamorphism (Brun, 1989)
-			if(W[i]>0){
+			if(W[i]>0.0){
 				//_printf_("Nond}ritic wet snow metamorphism\n");
 				//wet rate of change coefficient
-				IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *86400);   // [mm^3 s^-1]
+				IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *86400.0);   // [mm^3 s^-1]
 
 				// calculate change in grain volume and convert to grain size
-				gsz[i] = 2 * pow(3/(Pi * 4)*((4 / 3)*Pi*pow(gsz[i]/2,3.0) + E),1.0/3.0);
+				gsz[i] = 2.0 * pow(3.0/(Pi * 4.0)*((4.0/ 3.0)*Pi*pow(gsz[i]/2.0,3.0) + E),1.0/3.0);
 
 			}
 
 			// grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992)
-			if (gsp[i]==1.0 && gsz[i]>2) gsz[i]=2.0;
+			if (gsp[i]==1.0 && gsz[i]>2.0) gsz[i]=2.0;
 
 			// grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992)
@@ -294,5 +291,5 @@
 
 		//convert grain size back to effective grain radius: 
-		re[i]=gsz[i]/2;
+		re[i]=gsz[i]/2.0;
 	}
 	
@@ -505,8 +502,8 @@
 	IssmDouble  dt;
 	IssmDouble max_fdt=0;
-	IssmDouble  Ts;
+	IssmDouble  Ts=0;
 	IssmDouble  L;
 	IssmDouble  eS;
-	IssmDouble  Ri;
+	IssmDouble  Ri=0;
 	IssmDouble  coefM;
 	IssmDouble  coefH;
@@ -690,5 +687,5 @@
 
 	// PREALLOCATE ARRAYS BEFORE LOOP FOR IMPROVED PERFORMANCE
-	T0 = xNew<IssmDouble>(m+2);
+	T0 = xNewZeroInit<IssmDouble>(m+2);
 	Tu=xNew<IssmDouble>(m);
 	Td=xNew<IssmDouble>(m);
@@ -831,5 +828,5 @@
 
 }  /*}}}*/
-void shortwave(IssmDouble* swf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid){ /*{{{*/
+void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid){ /*{{{*/
 
 	// DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE
@@ -852,6 +849,14 @@
 	// Outputs
 	//   swf     = absorbed shortwave radiation [W m-2]
+	//
+	
+	/*outputs: */
+	IssmDouble* swf=NULL;
 
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   shortwave module\n");
+
+	/*Initialize and allocate: */
+	swf=xNewZeroInit<IssmDouble>(m);
+
 
 	// SHORTWAVE FUNCTION
@@ -1020,4 +1025,7 @@
 		}
 	}
+	/*Assign output pointers: */
+	*pswf=swf;
+
 } /*}}}*/ 
 void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid){ /*{{{*/
@@ -1075,11 +1083,11 @@
 	m=*pm;
 
-	/*Allocate: */
+	// determine initial mass
 	mInit=xNew<IssmDouble>(m);
+	for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
+	massinit=0; for(int i=0;i<m;i++)massinit+=mInit[i];
 
 	if (P > 0){
-		// determine initial mass
-		for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
-	
+			
 
 		if (T_air <= 273.15){ // if snow
@@ -1148,5 +1156,4 @@
 		// check for conservation of mass
 		mass=0; for(int i=0;i<m;i++)mass+=d[i]*dz[i]; 
-		massinit=0; for(int i=0;i<m;i++)massinit+=mInit[i];
 
 		mass_diff = mass - massinit - P;
@@ -1188,6 +1195,5 @@
 	IssmDouble* F=NULL;
 	IssmDouble* flxDn=NULL;
-	IssmDouble* R=NULL;
-	IssmDouble* ER=NULL;
+	IssmDouble  ER=0;
 	IssmDouble* EI=NULL;
 	IssmDouble* EW=NULL;
@@ -1207,4 +1213,5 @@
 	IssmDouble Wi;
 	int        D_size;
+	int         i;
 
 	/*outputs:*/
@@ -1220,4 +1227,5 @@
 	IssmDouble* gsp=*pgsp;
 	int         n=*pn;
+	IssmDouble* R=0;
 	
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
@@ -1246,10 +1254,10 @@
 
 	// initialize melt and runoff scalars
-	R = 0;          // runoff [kg]
+	Rsum = 0;       // runoff [kg]
 	sumM = 0;       // total melt [kg]
 	mAdd = 0;       // mass added/removed to/from base of model [kg]
 	addE = 0;       // energy added/removed to/from base of model [J]
 
-	// calculate temperature excess above 0 °C
+	// calculate temperature excess above 0 deg C
 	exsT=xNewZeroInit<IssmDouble>(n);
 	for(int i=0;i<n;i++) exsT[i]= fmax(0, T[i] - CtoK);        // [K] to [°C]
@@ -1264,5 +1272,5 @@
 	// check if any pore water
 	if (cellsum(W,n) > 0){
-		// _printf_("PORE WATER REFREEZE");
+		if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("      pore water refreeze\n");
 		// calculate maximum freeze amount, maxF [kg]
 		for(int i=0;i<n;i++) maxF[i] = fmax(0, -((T[i] - CtoK) * m[i] * CI) / LF);
@@ -1276,9 +1284,11 @@
 
 		// if pore water froze in ice then adjust d and dz thickness
-		for(int i=0;i<n;i++)if(d[i]>dIce)d[i]=m[i]/d[i];
+		for(int i=0;i<n;i++)if(d[i]>dIce)d[i]=dIce;
+		for(int i=0;i<n;i++) dz[i]= m[i]/d[i];
 	}
 
 	// squeeze water from snow pack
-	exsW=xNew<IssmDouble>(n); for(int i=0;i<n;i++){
+	exsW=xNew<IssmDouble>(n); 
+	for(int i=0;i<n;i++){
 		Wi= (910 - d[i]) * Swi * (m[i] / d[i]);        // irreducible water content [kg]
 		exsW[i] = fmax(0, W[i] - Wi);                  // water "squeezed" from snow [kg]
@@ -1305,5 +1315,5 @@
 			while (cellsum(surpE,n) > 0){
 				// use surplus energy to increase the temperature of lower cell
-				T[i+i] = surpE[i] * m[i+1]/CI + T[i+i];
+				T[i+1] = surpE[i] * m[i+1]/CI + T[i+1];
 				surpT[i+1] = fmax(0, (T[i+1] - CtoK - 159.1342));
 				surpE[i+1] = surpT[i+1] * CI / m[i+1];
@@ -1326,6 +1336,7 @@
 
 		// initialize refreeze, runoff, flxDn and dW vectors [kg]
-		F = xNewZeroInit<IssmDouble>(n); 
-		R = xNewZeroInit<IssmDouble>(n);
+		IssmDouble* F = xNewZeroInit<IssmDouble>(n); 
+		IssmDouble* R=xNewZeroInit<IssmDouble>(n); 
+
 		for(int i=0;i<n;i++)dW[i] = 0;
 		flxDn=xNewZeroInit<IssmDouble>(n+1); for(int i=0;i<n;i++)flxDn[i+1]=F[i];
@@ -1436,8 +1447,15 @@
 		celldelete(&gdn,n,D,D_size);
 		celldelete(&gsp,n,D,D_size);
+		celldelete(&EI,n,D,D_size);
+		celldelete(&EW,n,D,D_size);
+		celldelete(&R,n,D,D_size);
 		n=D_size;
 	
 		// calculate new grid lengths
 		for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
+
+		//calculate Rsum: 
+		Rsum=cellsum(R,n);
+	
 	}
 
@@ -1476,4 +1494,49 @@
 	xDelete<int>(D); D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)D_size++; D=xNew<int>(D_size); 
 	D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999){ D[D_size] = i; D_size++;}
+
+	celldelete(&m,n,D,D_size);
+	celldelete(&W,n,D,D_size);
+	celldelete(&dz,n,D,D_size);
+	celldelete(&d,n,D,D_size);
+	celldelete(&T,n,D,D_size);
+	celldelete(&a,n,D,D_size);
+	celldelete(&re,n,D,D_size);
+	celldelete(&gdn,n,D,D_size);
+	celldelete(&gsp,n,D,D_size);
+	celldelete(&EI,n,D,D_size);
+	celldelete(&EW,n,D,D_size);
+	n=D_size;
+
+	// check if any of the top 10 cell depths are too large
+	X=0;
+	for(int i=9;i>=0;i--){
+		if(dz[i]> 2* dzMin){
+			X=i;
+			break;
+		}
+	}
+	
+	i=0;
+	while(i<=X){
+		if (dz [i] > dzMin *2){
+
+			// _printf_("dz > dzMin * 2");
+			// split in two
+			cellsplit(&dz, n, i,.5);
+			cellsplit(&W, n, i,.5);
+			cellsplit(&m, n, i,.5);
+			cellsplit(&T, n, i,1.0);
+			cellsplit(&d, n, i,1.0);
+			cellsplit(&a, n, i,1.0);
+			cellsplit(&EI, n, i,1.0);
+			cellsplit(&EW, n, i,1.0);
+			cellsplit(&re, n, i,1.0);
+			cellsplit(&gdn, n, i,1.0);
+			cellsplit(&gsp, n, i,1.0);
+			n++;
+			X=X+1;
+		}
+		else i++;
+	}
 
 	//// CORRECT FOR TOTAL MODEL DEPTH
@@ -1516,27 +1579,24 @@
 
 	// calculate final mass [kg] and energy [J]
-	ER=xNew<IssmDouble>(n);
-	Rsum = cellsum(R,n);
+	sumER = Rsum * (LF + CtoK * CI);
 	for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
-	for(int i=0;i<n;i++)ER[i] = R[i] * (LF + CtoK * CI);
 	for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);
 
 	mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
 	sumE1 = cellsum(EI,n) + cellsum(EW,n);
-	sumER = cellsum(ER,n);
 
 	dm = round(mSum0 - mSum1 + mAdd);
 	dE = round(sumE0 - sumE1 - sumER +  addE);
-
-	if (dm !=0  && dE !=0) _printf_("mass and energy are not conserved in melt equations");
-	else if (dm != 0) _printf_("mass is not conserved in melt equations");
-	else if (dE != 0) _printf_("energy is not conserved in melt equations");
-
-	// W = round(W * 10000)/10000;
-	for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations");
-
+	
+	for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations\n");
+
+	if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n" 
+			<< "dm: " << dm << " dE: " << dE << "\n");
+
+	
 	/*Free ressources:*/
 	if(m)xDelete<IssmDouble>(m);
 	if(EI)xDelete<IssmDouble>(EI);
+	if(EW)xDelete<IssmDouble>(EW);
 	if(maxF)xDelete<IssmDouble>(maxF);
 	if(dW)xDelete<IssmDouble>(dW);
@@ -1545,8 +1605,6 @@
 	if(surpT)xDelete<IssmDouble>(surpT);
 	if(surpE)xDelete<IssmDouble>(surpE);
-	if(F)xDelete<IssmDouble>(F);
 	if(flxDn)xDelete<IssmDouble>(flxDn);
 	if(D)xDelete<int>(D);
-	if(R)xDelete<IssmDouble>(R);
 	if(M)xDelete<IssmDouble>(M);
 	
@@ -1630,5 +1688,5 @@
 	IssmDouble* obp = xNew<IssmDouble>(m);
 	obp[0]=0;
-	for(int i=1;i<m;i++)obp[i]=cumdz[i]*d[i];
+	for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
 	
 	// calculate new snow/firn density for:
@@ -1655,6 +1713,6 @@
 				// common variable
 				H = exp((-60/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000,2.0);
-				c0 = 9.2E-9 * H;
-				c1 = 3.7E-9 * H;
+				c0 = 9.2e-9 * H;
+				c1 = 3.7e-9 * H;
 				break;
 
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 19581)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 19582)
@@ -25,5 +25,5 @@
 void grainGrowth(IssmDouble* pre, IssmDouble* pgdn, IssmDouble* pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid);
 void albedo(IssmDouble* a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt,int m, int sid);
-void shortwave(IssmDouble* swf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid);
+void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid);
 void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, int sid);
 void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid);