Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 19582)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 19583)
@@ -2164,4 +2164,5 @@
 	int         m;
 	IssmDouble  SmbMassBalance=0;
+	int         count=0;
 	/*}}}*/
 
@@ -2246,5 +2247,4 @@
 		a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aSnow;         //set albedo equal to fresh snow [fraction] 
 		T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
-		swf = xNewZeroInit<IssmDouble>(m); 
 
 		//fixed lower temperatuer bounday condition - T is fixed
@@ -2265,5 +2265,4 @@
 		DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input);
 		DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input);
-		DoubleArrayInput* swf_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbSwfEnum));_assert_(swf_input);
 		
 		/*Recover arrays: */
@@ -2277,5 +2276,7 @@
 		a_input->GetValues(&a,&m);
 		T_input->GetValues(&T,&m);
-		swf_input->GetValues(&swf,&m);
+		
+		//fixed lower temperatuer bounday condition - T is fixed
+		T_bottom=T[m-1];
 
 	} /*}}}*/
@@ -2292,7 +2293,8 @@
 
 	/*Start loop: */
+	count=1;
 	for (t=time;t<=time+dt;t=t+smb_dt){
 
-		if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << "\n");
+		if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n");
 
 		/*extract daily data:{{{*/
@@ -2315,6 +2317,9 @@
 					
 		/*Distribution of absorbed short wave radation with depth:*/
-        if(isshortwave)shortwave(swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid());
+        if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid());
 		
+		/*Calculate net shortwave [W m-2]*/
+        netSW = cellsum(swf,m);
+
 		/*Thermal profile computation:*/
         if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,this->Sid());     
@@ -2338,10 +2343,24 @@
         ulw = 5.67E-8 * pow(T[0],4.0);
 
-		/*Calculate net shortwave and longwave [W m-2]*/
-        netSW = cellsum(swf,m);
+		/*Calculate net longwave [W m-2]*/
         netLW = dlw - ulw;
 		
 		/*Calculate turbulent heat fluxes [W m-2]*/
         if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,this->Sid());
+		
+		/*Verbose some resuls in debug mode: {{{*/
+		if(VerboseSmb() && 0){ 
+			_printf_("smb log: count[" << count << "] m[" << m << "] " 
+				<< setprecision(16)   << "T[" << cellsum(T,m)  << "] " 
+					                  << "d[" << cellsum(d,m)  << "] "
+					                  << "dz[" << cellsum(dz,m)  << "] "
+					                  << "a[" << cellsum(a,m)  << "] "
+					                  << "W[" << cellsum(W,m)  << "] "
+					                  << "re[" << cellsum(re,m)  << "] "
+					                  << "gdn[" << cellsum(gdn,m)  << "] "
+					                  << "gsp[" << cellsum(gsp,m)  << "] "
+					                  << "swf[" << netSW << "] "
+									  << "\n");
+		} /*}}}*/
 		
 		/*Sum component mass changes [kg m-2]*/
@@ -2359,11 +2378,12 @@
 
 		/*Check mass conservation:*/
-        if(!VerboseSmb())if (dMass != 0.0) _printf_("total system mass not conserved in MB function");
+        if (dMass != 0.0) _printf_("total system mass not conserved in MB function");
 		
 		/*Check bottom grid cell T is unchanged:*/
-        if(!VerboseSmb())if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
+        if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
 		
 		SmbMassBalance += sumP + sumEC - sumR; //increment SMB for the entire time span of ice-flow dynamics.
 
+		count++;
 
 	} //for (t=time;t<=time+dt;t=t+smb_dt)
@@ -2392,4 +2412,6 @@
 	xDelete<IssmDouble>(a);
 	xDelete<IssmDouble>(T);
+	xDelete<IssmDouble>(swf);
+	delete gauss;
 	/*}}}*/
 }
