Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 25992)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 25993)
@@ -3551,4 +3551,8 @@
 	IssmDouble teValue=1.0;
 	IssmDouble aValue=0.0;
+	IssmDouble szaValue=0.0;
+	IssmDouble cotValue=0.0;
+	IssmDouble ccsnowValue=0.0;
+	IssmDouble cciceValue=0.0;
 	IssmDouble dt,time,smb_dt;
 	int        aIdx=0;
@@ -3833,4 +3837,9 @@
 	teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
 	aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
+	//Not implemented as input
+	//szaValue;  // Solar Zenith Angle [degree]
+	//cotValue;  // Cloud Optical Thickness
+	//ccsnowValue;//concentration of light absorbing carbon for snow [ppm1]
+	//cciceValue;//concentration of light absorbing carbon for ice [ppm1]
 	//_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
 	/*}}}*/
@@ -3840,5 +3849,5 @@
 
 	/*Snow, firn and ice albedo:*/
-	if(isalbedo)albedo(&a,aIdx,re,dz,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,Msurf,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
+	if(isalbedo)albedo(&a,aIdx,re,dz,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,Msurf,ccsnowValue,cciceValue,szaValue,cotValue,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
 
 	/*Distribution of absorbed short wave radation with depth:*/
@@ -3927,5 +3936,6 @@
 	/*Check mass conservation:*/
 	if (dMass != 0.0){
-		_printf_("total system mass not conserved in MB function \n");
+		_printf_("Time: " << setprecision(8) << timeinputs/365.0/24.0/3600.0 << ",total system mass not conserved in MB function \n");
+		_printf_("sumMass: " << sumMass << " sumR: " << sumR << " sumW: " << sumW << " sumP: " << sumP << " sumEC: " << sumEC << " initMass: " << initMass << " sumMassAdd: " << sumMassAdd << " \n");
 	}
 	#endif
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 25992)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp	(revision 25993)
@@ -410,5 +410,5 @@
 
 }  /*}}}*/
-void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
+void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
 
 	//// Calculates Snow, firn and ice albedo as a function of:
@@ -432,4 +432,10 @@
 	// Methods 1 & 2
 	//   re      = surface effective grain radius [mm]
+	// Method 1, optional
+	//  clabSnow = concentration of light absorbing carbon  [ppm1], default 0
+	//  SZA      = solar zenith angle of the incident radiation [deg], default 0
+	//  COT      = cloud optical thickness, default 0
+	//  For TWO LAYER
+	//  clabIce  = concentration of light absorbing carbon of first ice layer [ppm1], default 0
 
 	// Methods 3
@@ -485,9 +491,75 @@
 			//function of effective grain radius
 
+			//gardnerAlb(S1, c1, SZA, t, z1, S2, c2)
+			//This is an implementation of the snow and ice broadband albedo
+			//  parameterization developed by Alex Gardner.
+			//Created By: Alex S. Gardner, Jet Propulsion Laboratory [alex.s.gardner@jpl.nasa.gov]
+			//  Last Modified: June, 2014
+			//Full Reference: Gardner, A. S., and Sharp, M. J.: A review of snow and
+			//  ice albedo and the development of a new physically based broadband albedo
+			//  parameterization, J. Geophys. Res., 115, F01009, 10.1029/2009jf001444,
+			//  2010.
+
+			//INPUTS
+			// ONE LAYER
+			//  - S1    : specific surface area of the snow or ice [cm^2 g-1]
+			//  - c1    : concentration of light absorbing carbon  [ppm1]
+			//  - SZA   : solar zenith angle of the incident radiation [deg]
+			//  - t     : cloud optical thickness
+			// TWO LAYER
+			//  - z1    : depth of snow suface layer [mm w.e.]
+			//  - S2    : specific surface area of bottom ice layer [cm^2 g-1]
+			//  - c2    : concentration of light absorbing carbon of bottom ice
+			//             layer [ppm1]
+			IssmDouble c1=clabSnow;
+			IssmDouble c2=clabIce;
+			IssmDouble t=COT;
+
+			//Single layer albedo parameterization
 			//convert effective radius to specific surface area [cm2 g-1]
-			IssmDouble S = 3.0 / (0.091 * re[0]);
-
-			//determine broadband albedo
-			a[0]= 1.48 - pow(S,-0.07);
+			IssmDouble S1 = 3.0 / (0.091 * re[0]);
+
+			//effective solar zenith angle
+			IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0);
+			IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0);
+
+			// pure snow albedo
+			IssmDouble as = 1.48 - pow(S1,-0.07);
+
+			//change in pure snow albedo due to soot loading
+			IssmDouble dac = max(0.04 - as, pow(-c1,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c1,0.6))*pow(S1,-0.25)));
+
+			//Two layer albedo parameterization
+			//  do two layer calculation if there is more than 1 layer
+			IssmDouble z1=0.0;
+			int lice=0;
+			for(int l=0;(l<m & d[l]<dPHC-Dtol);l++){
+				z1=z1+dz[l]*d[l]; //mm
+				lice=l+1;
+			}
+			if (m>0 & lice<m & z1 > Dtol){
+				// determine albedo values for bottom layer
+				IssmDouble S2 = 3.0 / (0.091 * re[lice]);
+
+				// pure snow albedo
+				IssmDouble as2 = 1.48 - pow(S2,-0.07);
+
+				// change in pure snow albedo due to soot loading
+				IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));
+
+				// determine the effective change due to finite depth and soot loading
+				IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13)));
+
+				dac =  (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2));
+			}
+
+			// change in albedo due to solar zenith angle
+			IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2);
+
+			// change in albedo due to cloud (apart from change in diffuse fraction)
+			IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as));
+
+			// Broadband albedo
+			a[0] = as + dac + dasz + dat;
 		}
 		else if(aIdx==2){
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 25992)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h	(revision 25993)
@@ -29,5 +29,5 @@
 IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
 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* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
+void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, 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* 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, bool isconstrainsurfaceT);