Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 19565)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 19566)
@@ -129,5 +129,5 @@
 
 } /*}}}*/
-void grainGrowth(IssmDouble* re, IssmDouble* gdn, IssmDouble* gsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx){ /*{{{*/
+void grainGrowth(IssmDouble* re, IssmDouble* gdn, IssmDouble* gsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx,int sid){ /*{{{*/
 
 	/*Created by: Alex S. Gardner, University of Alberta
@@ -173,4 +173,6 @@
 	IssmDouble  Q=0;
 
+	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   grain growth module\n");
+	
 	/*only when aIdx = 1 or 2 do we run grainGrowth: */
 	if(aIdx!=1 & aIdx!=2){
@@ -302,5 +304,5 @@
 
 }  /*}}}*/
-void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, int m) { /*{{{*/
+void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, int m,int sid) { /*{{{*/
 
 	//// Calculates Snow, firn and ice albedo as a function of:
@@ -345,4 +347,6 @@
 	// a = albedo(4, [], [], [], 0.48, 0.85, [0.8 0.5 ... 0.48], ...
 	//   [273 272.5 ... 265], [0 0.001 ... 0], 0, 0.01, 15, 15, 7, 3600)
+
+	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   albedo module\n");
         
 	//some constants:
@@ -451,5 +455,5 @@
 	else if (isnan(a[0])) _error_("albedo == NAN\n");
 }  /*}}}*/
-void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz) { /*{{{*/
+void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz,int sid) { /*{{{*/
 
 	/* ENGLACIAL THERMODYNAMICS*/
@@ -527,4 +531,6 @@
 	IssmDouble EC;
 
+	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   thermal module\n");
+
 	// INITIALIZE
 	CI = 2102;          // heat capacity of snow/ice (J kg-1 k-1)
@@ -621,10 +627,10 @@
 	// must go evenly into one hour or the data frequency if it is smaller
 
-	// all integer factors of the number of second in a day (8600 [s])
+	// all integer factors of the number of second in a day (86400 [s])
 	int f[45] = {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 30, 36, 40, 45, 48, 50, 60,
     72, 75, 80, 90, 100, 120, 144, 150, 180, 200, 225, 240, 300, 360, 400, 450, 600, 720, 900, 1200, 1800, 3600};
 
 	// return the min integer factor that is < dt
-	max_fdt=0;
+	max_fdt=f[0];
 	for(int i=0;i<45;i++){
 		if (f[i]<dt)if(f[i]>=max_fdt)max_fdt=f[i];
@@ -773,13 +779,12 @@
     
 		//SW penetrates surface
-		for(int i=0;i<m;i++) T[i] = T[i] + dT_sw[i];
+		for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
 		T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
 		
 		// temperature diffusion
-		for(int i=0;i<m;i++)T0[1+i]=T[i];
-		for(int i=0;i<m;i++) Tu[i] = T0[i];
-		for(int i=0;i<m;i++) Td[i] = T0[2+i];
-    
-		for(int i=0;i<m;i++) T[i] = (Np[i] * T[i]) + (Nu[i] * Tu[i]) + (Nd[i] * Td[i]);
+		for(int j=0;j<m;j++)T0[1+j]=T[j];
+		for(int j=0;j<m;j++) Tu[j] = T0[j];
+		for(int j=0;j<m;j++) Td[j] = T0[2+j];
+		for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);
 		
 		// calculate cumulative evaporation (+)/condensation(-)
@@ -826,5 +831,5 @@
 
 }  /*}}}*/
-void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m){ /*{{{*/
+void shortwave(IssmDouble* swf, 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
@@ -848,11 +853,7 @@
 	//   swf     = absorbed shortwave radiation [W m-2]
 
-	/*outputs:*/
-	IssmDouble* swf = NULL;
+	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   shortwave module\n");
 
 	// SHORTWAVE FUNCTION
-	// initialize variables
-	swf = xNewZeroInit<IssmDouble>(m);
-
 	if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell
         
@@ -1010,5 +1011,5 @@
 
 			// add flux absorbed at surface
-			swf[0] = swf[0] + swf_s;
+			swf[0] += swf_s;
 
 			/*Free ressources:*/
@@ -1019,8 +1020,6 @@
 		}
 	}
-	/*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){ /*{{{*/
+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){ /*{{{*/
 
 	// Adds precipitation and deposition to the model grid
@@ -1062,4 +1061,6 @@
 	IssmDouble* gsp=NULL;
 	int         m;
+
+	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   accumulation module\n");
 
 	/*Recover pointers: */
@@ -1169,5 +1170,5 @@
 	*pm=m;
 } /*}}}*/
-void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin){ /*{{{*/
+void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, int sid){ /*{{{*/
 
 	//// MELT ROUTINE
@@ -1220,4 +1221,6 @@
 	int         n=*pn;
 	
+	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
+
 	//// INITIALIZATION
 
@@ -1564,5 +1567,5 @@
 
 } /*}}}*/ 
-void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,int m){ /*{{{*/
+void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble dIce, int m, int sid){ /*{{{*/
 
 	//// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPATION IS COMPNSATED FOR BY TRACES OF SNOW???]
@@ -1606,5 +1609,4 @@
 	//// MAIN FUNCTION
 	// specify constants
-	const IssmDouble dIce    = 910;         // density of ice [kg m-3]
 	dt      = dt / 86400;  // convert from [s] to [d]
 	// R     = 8.314        // gas constant [mol-1 K-1]
@@ -1615,4 +1617,6 @@
 	/*intermediary: */
 	IssmDouble c0,c1,H;
+	
+	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   densification module\n");
 
 	// initial mass
@@ -1679,4 +1683,5 @@
 		dz[i] = mass_init[i] / d[i];
 	}
+
 	/*Free ressources:*/
 	xDelete<IssmDouble>(mass_init);
@@ -1685,5 +1690,5 @@
 
 } /*}}}*/
-void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz){ /*{{{*/
+void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, int sid){ /*{{{*/
 
 	//// TURBULENT HEAT FLUX
@@ -1725,4 +1730,6 @@
 	/*output: */
 	IssmDouble shf, lhf, EC;
+	
+	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   turbulentFlux module\n");
 
 	// calculated air density [kg/m3]