Index: /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp	(revision 28102)
+++ /issm/trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp	(revision 28103)
@@ -10,5 +10,5 @@
 
 IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec,
-				 IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu, 
+				 IssmDouble* pdds, IssmDouble* pds, IssmDouble* melt, IssmDouble* accu,  
 				 IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac,
 				 IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps,IssmDouble rlapslgm,
@@ -70,5 +70,7 @@
   IssmDouble pddtj, hmx2;
   IssmDouble pddsnowfac0=4.3, pddicefac0=8.3;
-  IssmDouble snowfac, icefac;
+  IssmDouble snowfac, icefac; 
+  
+  IssmDouble snow;
 
   sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
@@ -159,6 +161,7 @@
   else if(Tsum< 10){
 	  snwmf = (0.15*(Tsum+1) + (2.65+snowfac-pddsnowfac0))*0.001;
-	  smf = (((17.22-icefac)/(pow(11,3)))*pow((10.-Tsum),3) + pddicefac0)*0.001;
-	  //JC,smf = (((icefac-pddicefac0)/(Tsum+1))*pow((10.-Tsum),3) + pddicefac0)*0.001;
+	  smf = (((17.22-pddicefac0)/(pow(11,3)))*pow((10.-Tsum),3) + icefac)*0.001;//icefac)*0.001;//pddicefac0)*0.001;
+	  //2024: Prior Calculation SMF:  smf = (((17.22-icefac)/(pow(11,3)))*pow((10.-Tsum),3) + pddicefac0)*0.001;
+	  //Original: Tarasov and Peltier: smf = (((icefac-pddicefac0)/(Tsum+1))*pow((10.-Tsum),3) + pddicefac0)*0.001;
   }
   else{
@@ -246,4 +249,10 @@
   pddtj=pddt;
 
+  //cout << sizeof(snwm);
+  //snw = mean(snwm);
+
+  //std::cout << snwm << std::endl;
+  //printf("snowmelt %f\n",snow);
+  
   return B;
 }