Index: /issm/trunk-jpl/src/m/classes/qmustatistics.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/qmustatistics.py	(revision 25697)
+++ /issm/trunk-jpl/src/m/classes/qmustatistics.py	(revision 25698)
@@ -85,10 +85,10 @@
             for f in range(len(m['fields'])):
                 if not isinstance(m['fields'][f], str):
-                    raise Exception('qmustatistics consistency check error: qmu.statistics.method[{}][\'fields\'[{}] is not a string!'.format(i, f))
+                    raise Exception('qmustatistics consistency check error: qmu.statistics.method[{}][\'fields\'][{}] is not a string!'.format(i, f))
             for s in range(len(m['steps'])):
                 if m['steps'][s] <= 0:
-                    raise Exception('qmustatistics consistency check error: qmu.statistics.method[{}][\'steps\'[{}] should be > 0!'.format(i, s))
+                    raise Exception('qmustatistics consistency check error: qmu.statistics.method[{}][\'steps\'][{}] should be > 0!'.format(i, s))
                 if m['steps'][s] > md.mesh.numberofvertices:
-                    raise Exception('qmustatistics consistency check error: qmu.statistics.method[{}][\'steps\'[{}] should be < md.mesh.numberofvertices!'.format(i, s))
+                    raise Exception('qmustatistics consistency check error: qmu.statistics.method[{}][\'steps\'][{}] should be < md.mesh.numberofvertices!'.format(i, s))
     #}}}
 
Index: /issm/trunk-jpl/src/m/classes/results.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/results.py	(revision 25697)
+++ /issm/trunk-jpl/src/m/classes/results.py	(revision 25698)
@@ -2,59 +2,132 @@
 
 
-class results(object):
+class results(object): #{{{
     """RESULTS class definition
 
     Usage:
-        results = results()
-
-    TODO:
-    - Modify output so that it matches that of 
-
-        disp(md.results.<<solutionstring>>)
-
-    where <<solutionstring>> is one of the values from solve.m
+        md.results = results()
     """
 
-    def __init__(self, *args):  # {{{
+    def __init__(self, *args): #{{{
         pass
-    # }}}
+    #}}}
 
-    def __repr__(self):  # {{{
-        s = "   Model results:\n"
-
+    def __repr__(self): #{{{
+        s = ''
         if 'step' in self.__dict__:
-            s += "%s\n" % fielddisplay(self, 'step', "step number")
+            s += '{}\n'.format(fielddisplay(self, 'step', "step number"))
         if 'time' in self.__dict__:
-            s += "%s\n" % fielddisplay(self, 'time', "time value")
+            s += '{}\n'.format(fielddisplay(self, 'time', "time value"))
         if 'SolutionType' in self.__dict__:
-            s += "%s\n" % fielddisplay(self, 'SolutionType', "solution type")
+            s += '{}\n'.format(fielddisplay(self, 'SolutionType', "solution type"))
 
         for name in list(self.__dict__.keys()):
             if name not in ['step', 'time', 'SolutionType', 'errlog', 'outlog']:
                 if isinstance(getattr(self, name), list):
-                    s += "%s\n" % fielddisplay(self, name, "model results list")
+                    s += '{}\n'.format(fielddisplay(self, name, "model results list"))
                 elif isinstance(getattr(self, name), results):
-                    s += "%s\n" % fielddisplay(self, name, "model results case")
+                    s += '{}\n'.format(fielddisplay(self, name, "model results case"))
                 else:
-                    s += "%s\n" % fielddisplay(self, name, "")
+                    s += '{}\n'.format(fielddisplay(self, name, ""))
 
         if 'errlog' in self.__dict__:
-            s += "%s\n" % fielddisplay(self, 'errlog', "error log file")
+            s += '{}\n'.format(fielddisplay(self, 'errlog', "error log file"))
         if 'outlog' in self.__dict__:
-            s += "%s\n" % fielddisplay(self, 'outlog', "output log file")
+            s += '{}\n'.format(fielddisplay(self, 'outlog', "output log file"))
 
         return s
-    # }}}
+    #}}}
 
-    def setdefaultparameters(self):  # {{{
+    def setdefaultparameters(self): #{{{
         #do nothing
         return self
-    # }}}
+    #}}}
 
-    def checkconsistency(self, md, solution, analyses):  # {{{
+    def checkconsistency(self, md, solution, analyses): #{{{
         return md
-    # }}}
+    #}}}
 
-    def marshall(self, prefix, md, fid):  # {{{
+    def marshall(self, prefix, md, fid): #{{{
         pass
-    # }}}
+    #}}}
+#}}}
+
+class resultselement(object): #{{{
+    """RESULTSELEMENT class definition - Value of attribute of results
+
+    Usage:
+        setattr(md.results, '<solution_type>', resultselement())
+    """
+
+    def __init__(self, *args): #{{{
+        self.elements = [result()]
+    #}}}
+
+    def __repr__(self): #{{{
+        s = '  1x{} struct array with fields:\n'.format(len(self.elements))
+        s += '\n'
+
+        if len(self.elements) >= 1:
+            for name in list(self.elements[0].__dict__.keys()):
+                s += '    {}\n'.format(name)
+
+        return s
+    #}}}
+
+    def __len__(self): #{{{
+        return len(self.elements)
+    #}}}
+
+    def __getitem__(self, index): #{{{
+        while True:
+            try:
+                return self.elements[index]
+            except:
+                self.elements.append(result())
+    #}}}
+
+    def __setitem__(self, index, value): #{{{
+        while True:
+            try:
+                self.elements[index] = value
+                break
+            except:
+                self.elements.append(result())
+    #}}}
+
+    def setdefaultparameters(self): #{{{
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): #{{{
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        pass
+    #}}}
+#}}}
+
+class result(object): #{{{
+    """RESULT class definition - Element of resultselement::elements
+
+    Usage:
+        resultselement.elements.append(result())
+    """
+
+    def __init__(self, *args): #{{{
+        pass
+    #}}}
+
+    def setdefaultparameters(self): #{{{
+        return self
+    #}}}
+
+    def checkconsistency(self, md, solution, analyses): #{{{
+        return md
+    #}}}
+
+    def marshall(self, prefix, md, fid): #{{{
+        pass
+    #}}}
+#}}}
Index: /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.m
===================================================================
--- /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.m	(revision 25697)
+++ /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.m	(revision 25698)
@@ -65,5 +65,5 @@
 	end
 
-%post processes qmu results if necessary
+%postprocess qmu results if necessary
 else
 	md=postqmu(md);
Index: /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py	(revision 25697)
+++ /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.py	(revision 25698)
@@ -3,5 +3,5 @@
 import numpy as np # Remove: used temporarily for debugging
 
-from helpers import fieldnames
+from helpers import fieldnames, struct
 from parseresultsfromdisk import parseresultsfromdisk
 from postqmu import postqmu
@@ -34,6 +34,6 @@
             raise OSError(err_msg)
 
-        # Initialize md.results if not a structure yet
-        if not isinstance(md.results, results):
+        # # Initialize md.results if not a structure yet
+        if not isinstance(md.results, struct):
             md.results = results()
 
@@ -73,5 +73,5 @@
             setattr(md.results, structure[0].SolutionType, structure[0])
 
-    # Post processes QMU results if necessary
+    # Postprocess QMU results if necessary
     else:
         md = postqmu(md)
Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py
===================================================================
--- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py	(revision 25697)
+++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py	(revision 25698)
@@ -4,5 +4,5 @@
 import numpy as np
 
-import results as resultsclass
+from results import resultselement
 
 
@@ -25,5 +25,5 @@
     saveres = []
 
-    #if we have done split I / O, ie, we have results that are fragmented across patches,
+    #if we have done split I/O, ie, we have results that are fragmented across patches,
     #do a first pass, and figure out the structure of results
     loadres = ReadDataDimensions(fid)
@@ -83,5 +83,5 @@
     except IOError as e:
         raise IOError("parseresultsfromdisk error message: could not open {} for binary reading".format(filename))
-        
+
     # Collect all results in a list
     allresults = []
@@ -108,7 +108,5 @@
 
     # Ok, now construct structure
-    results = []
-    for i in range(len(allsteps)):
-        results.append(resultsclass.results())
+    results = resultselement()
     for i in range(numresults):
         result = allresults[i]
Index: /issm/trunk-jpl/test/NightlyRun/test2006.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2006.m	(revision 25697)
+++ /issm/trunk-jpl/test/NightlyRun/test2006.m	(revision 25698)
@@ -182,8 +182,6 @@
 md=solve(md,'Transient');
 
-disp(mds.results.StatisticsSolution)
-
 %compare statistics with our own here:
-svalues=mds.results.StatisticsSolution(end).SealevelSamples; %all values at locations. 
+svalues=mds.results.StatisticsSolution(end).SealevelSamples; %all values at locations.
 
 dvalues=zeros(md.qmu.method.params.samples,length(locations));
@@ -191,7 +189,4 @@
 	dvalues(i,:)=md.results.dakota.modelresults{i}.TransientSolution(end).Sealevel(locations);
 end
-
-disp(dvalues)
-disp(size(dvalues))
 
 samplesnorm=norm(dvalues-svalues,'fro');
Index: /issm/trunk-jpl/test/NightlyRun/test2006.py
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2006.py	(revision 25697)
+++ /issm/trunk-jpl/test/NightlyRun/test2006.py	(revision 25698)
@@ -135,5 +135,5 @@
 md.qmu.variables.surfaceload = qmuvar.surfaceload
 
-locations = [1, 5, 10, 15, 20]
+locations = np.array([1, 5, 10, 15, 20])
 
 # Responses #{{{
@@ -208,5 +208,5 @@
 md.qmu.statistics.method[2]['fields'] = ['Sealevel', 'BslrIce']
 md.qmu.statistics.method[2]['steps'] = np.arange(1, 10 + 1).reshape(1, -1)
-md.qmu.statistics.method[2]['indices'] = np.asarray(locations).reshape(1, -1)
+md.qmu.statistics.method[2]['indices'] = locations.reshape(1, -1)
 #}}}
 
@@ -219,13 +219,9 @@
 
 # Compare statistics with our own here
-#
-# TODO: SealevelSamples should be an attribute of mds.results.StatisticsSolution[-1] (see test2006.m)
-#
-svalues = mds.results.StatisticsSolution[0].SealevelSamples # all values at locations
+svalues = mds.results.StatisticsSolution[-1].SealevelSamples # all values at locations
 
 dvalues = np.zeros((md.qmu.method.params.samples, len(locations)))
-
 for i in range(md.qmu.method.params.samples):
-    dvalues[i, :] = md.results.dakota.modelresults[i].TransientSolution[-1].Sealevel[locations].flatten()
+    dvalues[i, :] = md.results.dakota.modelresults[i].TransientSolution[-1].Sealevel[locations - 1].flatten()
 
 samplesnorm = np.linalg.norm(dvalues - svalues, 'fro')