Index: /issm/trunk-jpl/src/m/classes/SMBautoregression.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/SMBautoregression.py	(revision 26820)
+++ /issm/trunk-jpl/src/m/classes/SMBautoregression.py	(revision 26821)
@@ -121,11 +121,11 @@
         temprefelevation  = np.copy(md.smb.refelevation)
         if(np.any(np.isnan(md.smb.lapserate_pos))):
-            templapserate_pos = np.zeros((md.smb.num_basins))
+            templapserate_pos = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins)
             print('      smb.lapserate_pos not specified: set to 0')
         if(np.any(np.isnan(md.smb.lapserate_neg))):
-            templapserate_neg = np.zeros((md.smb.num_basins))
+            templapserate_neg = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins)
             print('      smb.lapserate_neg not specified: set to 0')
         if(np.any(np.isnan(md.smb.refelevation))):
-            temprefelevation = np.zeros((md.smb.num_basins))
+            temprefelevation = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins)
             areas = GetAreas(md.mesh.elements, md.mesh.x, md.mesh.y)
             for ii,bid in enumerate(np.unique(md.smb.basin_id)):
@@ -133,13 +133,6 @@
                 elemsh  = np.zeros((len(indices)))
                 for jj in range(len(indices)):
-                    #elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj],:]-1])
-                    sumval = 0
-                    cntval = 0
-                    elemindices = md.mesh.elements[indices[jj],:]-1
-                    for kk in range(len(elemindices)):
-                        sumval += md.geometry.surface[elemindices[kk]]
-                        cntval += 1
-                    elemsh[jj] = sumval/cntval
-                temprefelevation[ii] = np.sum(areas[indices]*elemsh)/np.sum(areas[indices])
+                    elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj],:]-1])
+                temprefelevation[0,ii] = np.sum(areas[indices]*elemsh)/np.sum(areas[indices])
             if(np.any(templapserate_pos!=0) or np.any(templapserate_neg!=0)):
                 print('      smb.refelevation not specified: Reference elevations set to mean surface elevation of basins')