Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 23845)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 23846)
@@ -3241,7 +3241,4 @@
 	int        swIdx=0;
 	IssmDouble cldFrac,t0wet, t0dry, K;
-	IssmDouble ulw=0.0;
-	IssmDouble netSW=0.0;
-	IssmDouble netLW=0.0;
 	IssmDouble lhf=0.0;
 	IssmDouble shf=0.0;
@@ -3274,4 +3271,10 @@
 	IssmDouble* gsp = NULL;
 	IssmDouble  EC = 0.0;
+	IssmDouble  ulw = 0.0;
+	IssmDouble  netSW=0.0;
+	IssmDouble  netLW=0.0;
+	IssmDouble  netULW=0.0;
+	IssmDouble  netLHF=0.0;
+	IssmDouble  netSHF=0.0;
 	IssmDouble* W = NULL;
 	IssmDouble* a = NULL;
@@ -3523,8 +3526,8 @@
 
 		/*Calculate net shortwave [W m-2]*/
-		netSW = cellsum(swf,m);
+		netSW = netSW + cellsum(swf,m);
 
 		/*Thermal profile computation:*/
-		if(isthermal)thermo(&EC, &T, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
+		if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
 
 		/*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
@@ -3544,8 +3547,9 @@
 		/*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
 		 * sub-time step in thermo equations*/
-		ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
+		//ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
 
 		/*Calculate net longwave [W m-2]*/
-		netLW = dlw - ulw;
+		netULW = netULW + ulw;
+		netLW = netLW + dlw - ulw;
 
 		/*Calculate turbulent heat fluxes [W m-2]*/
@@ -3564,8 +3568,12 @@
 						<< "gsp[" << cellsum(gsp,m)  << "] "
 						<< "swf[" << netSW << "] "
+						<< "lwf[" << netLW << "] "
 						<< "a[" << a << "] "
 						<< "te[" << teValue << "] "
 						<< "\n");
 		} /*}}}*/
+
+		netLHF = netLHF + lhf;
+		netSHF = netSHF + shf;
 
 		/*Sum component mass changes [kg m-2]*/
@@ -3620,4 +3628,9 @@
 	this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/dt/rho_ice));
 	this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/dt/rho_ice));
+	this->AddInput(new DoubleInput(SmbNetULWEnum,netULW));
+	this->AddInput(new DoubleInput(SmbNetLWEnum,netLW));
+	this->AddInput(new DoubleInput(SmbNetSWEnum,netSW));
+	this->AddInput(new DoubleInput(SmbNetLHFEnum,netLHF));
+	this->AddInput(new DoubleInput(SmbNetSHFEnum,netSHF));
 	this->AddInput(new DoubleInput(SmbDzAddEnum,sumdz_add));
 	this->AddInput(new DoubleInput(SmbMAddEnum,sumMassAdd/dt));
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 23845)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 23846)
@@ -509,5 +509,5 @@
 
 }  /*}}}*/
-void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid) { /*{{{*/
+void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid) { /*{{{*/
 
 	/* ENGLACIAL THERMODYNAMICS*/
@@ -536,4 +536,5 @@
 	// T: grid cell temperature [k]
 	// EC: evaporation/condensation [kg]
+	// ulwrf: upward longwave radiation flux [W m-2]
 
 	/*intermediary: */
@@ -583,4 +584,5 @@
 	IssmDouble EC=0.0;
 	IssmDouble* T=*pT;
+	IssmDouble ulwrf=0.0;
 
 	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   thermal module\n");
@@ -610,5 +612,6 @@
 
 	// Bulk-transfer coefficient for turbulent fluxes
-	An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
+	//An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
+	An =  pow(0.4,2) / (log(Tz/z0)*log(Vz/z0));     // Bulk-transfer coefficient
 	C = An * dAir * V;                        // shf & lhf common coefficient
 
@@ -820,5 +823,7 @@
 
 		// upward longwave contribution
-		ulw = - (SB * pow(Ts,4.0)* teValue) * dt ; //+20
+		ulw = - (SB * pow(Ts,4.0)* teValue) * dt; // - deltatest here
+		ulwrf = ulwrf - ulw;
+
 		dT_ulw = ulw / TCs;
 
@@ -875,4 +880,5 @@
 	*pEC=EC;
 	*pT=T;
+	*pulwrf=ulwrf;
 
 }  /*}}}*/
@@ -1939,5 +1945,6 @@
 			case 6: // Ligtenberg and others (2011) [semi-emperical], Antarctica 
 				// common variable
-				H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
+				// From literature: H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
+				H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
 				c0arth = 0.07 * H;
 				c1arth = 0.03 * H;
@@ -1952,5 +1959,6 @@
 			case 7: // Kuipers Munneke and others (2015) [semi-emperical], Greenland
 				// common variable
-				H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
+				// From literature: H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
+				H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
 				c0arth = 0.07 * H;
 				c1arth = 0.03 * H;
@@ -2068,5 +2076,6 @@
 
 	//// Bulk-transfer coefficient
-	An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
+	//An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
+	An =  pow(0.4,2) / (log(Tz/z0)*log(Vz/z0));     // Bulk-transfer coefficient
 	C = An * d_air * V;             // shf & lhf common coefficient
 
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 23845)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 23846)
@@ -31,5 +31,5 @@
 void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
 void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, 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 teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid);
+void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid);
 void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid);
 void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid);
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 23845)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 23846)
@@ -1135,4 +1135,9 @@
 	SmbDzAddEnum,
 	SmbFACEnum,
+	SmbNetULWEnum,
+	SmbNetLWEnum,
+	SmbNetSWEnum,
+	SmbNetLHFEnum,
+	SmbNetSHFEnum,
 	SMBforcingEnum,
 	SMBgcmEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 23845)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 23846)
@@ -1139,4 +1139,9 @@
 		case SmbDzAddEnum : return "SmbDzAdd";
 		case SmbFACEnum : return "SmbFAC";
+		case SmbNetULWEnum : return "SmbNetULW";
+		case SmbNetLWEnum : return "SmbNetLW";
+		case SmbNetSWEnum : return "SmbNetSW";
+		case SmbNetLHFEnum : return "SmbNetLHF";
+		case SmbNetSHFEnum : return "SmbNetSHF";
 		case SMBforcingEnum : return "SMBforcing";
 		case SMBgcmEnum : return "SMBgcm";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 23845)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 23846)
@@ -1166,4 +1166,9 @@
 	      else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
 	      else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
+	      else if (strcmp(name,"SmbNetULW")==0) return SmbNetULWEnum;
+	      else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
+	      else if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum;
+	      else if (strcmp(name,"SmbNetLHF")==0) return SmbNetLHFEnum;
+	      else if (strcmp(name,"SmbNetSHF")==0) return SmbNetSHFEnum;
 	      else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
 	      else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
@@ -1239,13 +1244,13 @@
 	      else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
 	      else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
-	      else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
 	      else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
 	      else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
 	      else if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum;
 	      else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
+	      else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
 	      else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
 	      else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;