Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 673)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 674)
@@ -121,5 +121,4 @@
 	md.surface=NaN;
 	md.thickness=NaN;
-	md.new_thickness=NaN;
 	md.bed=NaN;
 
@@ -168,6 +167,4 @@
 	md.ndt=0;
 	md.artificial_diffusivity=0;
-	md.transient_results=NaN;
-	md.thermaltransient_results=NaN;
 	md.timestepping=0;
 	md.deltaH=0;
@@ -279,5 +276,5 @@
 
 	%results
-	md.results=NaN;
+	md.results=struct();
 
 	%Ice solver string
Index: /issm/trunk/src/m/solutions/cielo/loadresults.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/loadresults.m	(revision 673)
+++ /issm/trunk/src/m/solutions/cielo/loadresults.m	(revision 674)
@@ -15,23 +15,28 @@
 
 if dim==2,
-	
-	md.vx=u_g(1:2:end)*md.yts;
-	md.vy=u_g(2:2:end)*md.yts;
-	md.vel=sqrt(md.vx.^2+md.vy.^2);
-	md.pressure=p_g;
+	results.step=1;
+	results.time=1;
+	results.vx=u_g(1:2:end)*md.yts;
+	results.vy=u_g(2:2:end)*md.yts;
+	results.vel=sqrt(results.vx.^2+results.vy.^2);
+	results.pressure=p_g;
+	md.results.diagnostic=results;
 
 else
 	
+	results.step=1;
+	results.time=1;
 	if isstokes,
-		md.vx=u_g(1:4:end)*md.yts;
-		md.vy=u_g(2:4:end)*md.yts;
-		md.vz=u_g(3:4:end)*md.yts;
-		md.pressure=p_g;
+		results.vx=u_g(1:4:end)*md.yts;
+		results.vy=u_g(2:4:end)*md.yts;
+		results.vz=u_g(3:4:end)*md.yts;
+		results.pressure=p_g;
 	else
-		md.vx=u_g(1:3:end)*md.yts;
-		md.vy=u_g(2:3:end)*md.yts;
-		md.vz=u_g(3:3:end)*md.yts;
-		md.pressure=p_g;
+		results.vx=u_g(1:3:end)*md.yts;
+		results.vy=u_g(2:3:end)*md.yts;
+		results.vz=u_g(3:3:end)*md.yts;
+		results.pressure=p_g;
 	end
-	md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
+	results.vel=sqrt(results.vx.^2+results.vy.^2+results.vz.^2);
+	md.results.diagnostic=results;
 end
Index: /issm/trunk/src/m/solutions/cielo/prognostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 673)
+++ /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 674)
@@ -34,5 +34,7 @@
 	
 	displaystring(md.debug,'\n%s',['load results...']);
-	md.new_thickness=h_g;
+	md.results.prognostic.step=1;
+	md.results.prognostic.time=0;
+	md.results.prognostic.thickness=h_g;
 
 	%stop timing
Index: /issm/trunk/src/m/solutions/cielo/thermal.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 673)
+++ /issm/trunk/src/m/solutions/cielo/thermal.m	(revision 674)
@@ -36,6 +36,8 @@
 		
 		displaystring(md.debug,'\n%s',['load results...']);
-		md.temperature=t_g;
-		md.melting=m_g*md.yts; %from m/s to m/a
+		md.results.thermal.step=1;
+		md.results.thermal.time=0;
+		md.results.thermal.temperature=t_g;
+		md.results.thermal.melting=m_g*md.yts; %from m/s to m/a
 
 	else
@@ -73,6 +75,7 @@
 			solution(n).melting=soln(n).m_g*md.yts; %in m/year
 			solution(n).time=soln(n).time/md.yts;         %in years
+			solution(n).step=n;
 		end
-		md.thermaltransient_results=solution;
+		md.results.thermal=solution;
 	end
 
Index: /issm/trunk/src/m/solutions/ice/Loadresults.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/Loadresults.m	(revision 673)
+++ /issm/trunk/src/m/solutions/ice/Loadresults.m	(revision 674)
@@ -27,10 +27,13 @@
 	end
 
-	md.vx=u_g(indx)*md.yts;
-	md.vy=u_g(indy)*md.yts;
-	md.vel=sqrt(md.vx.^2+md.vy.^2);
-	md.pressure=u_g(indp);
+	md.results.diagnostic.step=1;
+	md.results.diagnostic.time=0;
+	md.results.diagnostic.vx=u_g(indx)*md.yts;
+	md.results.diagnostic.vy=u_g(indy)*md.yts;
+	md.results.diagnostic.vel=sqrt(md.vx.^2+md.vy.^2);
+	md.results.diagnostic.pressure=u_g(indp);
 
 else
+
 	m_dv=fem.m_dv;
 	indx=1:6:m_dv.gridset.gsize;
@@ -43,8 +46,10 @@
 	indp=indp(m_dv.tpart);
 
-	md.vx=u_g(indx)*md.yts;
-	md.vy=u_g(indy)*md.yts;
-	md.vz=u_g(indz)*md.yts;
-	md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
-	md.pressure=u_g(indp);
+	md.results.diagnostic.step=1;
+	md.results.diagnostic.time=0;
+	md.results.diagnostic.vx=u_g(indx)*md.yts;
+	md.results.diagnostic.vy=u_g(indy)*md.yts;
+	md.results.diagnostic.vz=u_g(indz)*md.yts;
+	md.results.diagnostic.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
+	md.results.diagnostic.pressure=u_g(indp);
 end
Index: /issm/trunk/src/m/solutions/ice/prognostic.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/prognostic.m	(revision 673)
+++ /issm/trunk/src/m/solutions/ice/prognostic.m	(revision 674)
@@ -59,8 +59,10 @@
 %project collapsed thickness onto 3d mesh
 indx=1:6:m.gridset.gsize; indx=indx(m.tpart);
+md.results.prognostic.step=1;
+md.results.prognostic.time=0;
 if strcmpi(md.type,'3d'),
-	md.new_thickness=project3d(md,project2d(md,h_g(indx),1),'node');
+	md.results.prognostic.thickness=project3d(md,project2d(md,h_g(indx),1),'node');
 else
-	md.new_thickness=h_g(indx);
+	md.results.prognostic.thickness=h_g(indx);
 end
 
Index: /issm/trunk/src/m/solutions/ice/thermal.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/thermal.m	(revision 673)
+++ /issm/trunk/src/m/solutions/ice/thermal.m	(revision 674)
@@ -45,6 +45,8 @@
 	indx=1:6:gridset.gsize;
 	indx=indx(m_t.tpart);
-	md.temperature=t_g(indx);
-	md.melting=melting_g(indx)*md.yts; %from m/s to m/a
+	md.results.thermal.step=1;
+	md.results.thermal.time=0;
+	md.results.thermal.temperature=t_g(indx);
+	md.results.thermal.melting=melting_g(indx); %from m/s to m/a
 else
 
@@ -86,9 +88,10 @@
 
 	for n=1:nsteps+1,
+		solution(n).step=n;
 		solution(n).temperature=soln(n).t_g(indx);
 		solution(n).melting=soln(n).melting_g(indx)*md.yts; %in m/a
 		solution(n).time=soln(n).time/md.yts; %in year
 	end
-	md.thermaltransient_results=solution;
+	md.results.thermal=solution;
 end
 
Index: /issm/trunk/src/m/solutions/ice/transient2d.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/transient2d.m	(revision 673)
+++ /issm/trunk/src/m/solutions/ice/transient2d.m	(revision 674)
@@ -108,4 +108,5 @@
 	solution2(i).bed=solution(i).bed(indx);
 	solution2(i).time=solution(i).time/md.yts;
+	solution2(i).step=i;
 end
-md.transient_results=solution2;
+md.results.transient=solution2;
Index: /issm/trunk/src/m/solutions/ice/transient3d.m
===================================================================
--- /issm/trunk/src/m/solutions/ice/transient3d.m	(revision 673)
+++ /issm/trunk/src/m/solutions/ice/transient3d.m	(revision 674)
@@ -162,5 +162,5 @@
 	solution2(i).bed=solution(i).bed(indx);
 	solution2(i).time=solution(i).time/md.yts;
+	solution2(i).step=i;
 end
-
-md.transient_results=solution2;
+md.results.transient=solution2;
Index: /issm/trunk/src/m/solutions/macayeal/control.m
===================================================================
--- /issm/trunk/src/m/solutions/macayeal/control.m	(revision 673)
+++ /issm/trunk/src/m/solutions/macayeal/control.m	(revision 674)
@@ -123,6 +123,4 @@
 end
 
-
-
 %check that the model is not a pure ice shelf (no friction on ice shelves)
 if strcmp(md.control_type,'drag') & (length(find(drag_coeff))==0),
@@ -168,5 +166,4 @@
        direction=direction/sqrt(direction'*direction);
 
-
        % rough orthagonalization
        direction=direction-(direction'*old_direction)*old_direction; 
@@ -242,10 +239,12 @@
     end
 
-    %Load results onto md:
-    md.cont_J=J;
-	 md.cont_parameter=B;
-	 md.cont_vx=u;
-	 md.cont_vy=v;
-	 md.cont_vel=sqrt(u.^2+v.^2);
+	 %Load results onto md:
+	 md.results.control.step=1;
+	 md.results.control.time=0;
+    md.results.control.J=J;
+	 md.results.control.parameter=B;
+	 md.results.control.vx=u;
+	 md.results.control.vy=v;
+	 md.results.control.vel=sqrt(u.^2+v.^2);
 
 elseif strcmp(md.control_type,'drag'),
Index: /issm/trunk/src/m/solutions/macayeal/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/macayeal/diagnostic.m	(revision 673)
+++ /issm/trunk/src/m/solutions/macayeal/diagnostic.m	(revision 674)
@@ -468,7 +468,9 @@
 
 %Load results onto md:
-md.vx=u*yts;
-md.vy=v*yts;
-md.vel=sqrt(md.vx.^2+md.vy.^2);
+md.results.diagnostic.step=1;
+md.results.diagnostic.time=0;
+md.results.diagnostic.vx=u*yts;
+md.results.diagnostic.vy=v*yts;
+md.results.diagnostic.vel=sqrt(u.^2+v.^2)*yts;
 
 %stop timing
Index: /issm/trunk/src/m/utils/Nightly/testsgetfields.m
===================================================================
--- /issm/trunk/src/m/utils/Nightly/testsgetfields.m	(revision 673)
+++ /issm/trunk/src/m/utils/Nightly/testsgetfields.m	(revision 674)
@@ -15,27 +15,27 @@
 if strcmpi(string,'diagnostic'),
 	if strcmpi(mdtype,'3d')
-		fields={'vy','vz'};
+		fields={'diagnostic.vy','diagnostic.vz'};
 	else
-		fields={'vel'};
+		fields={'diagnostic.vel'};
 	end
 
 elseif strcmpi(string,'thermalsteady'),
-	fields={'temperature','melting'};
+	fields={'thermal.temperature','thermal.melting'};
 
 elseif strcmpi(string,'thermaltransient'),
-	fields={'thermaltransient_results(end).temperature','thermaltransient_results(end).melting'};
+	fields={'thermal(end).temperature','thermal(end).melting'};
 
 elseif strcmpi(string,'prognostic'),
-	fields={'new_thickness'};
+	fields={'prognostic.thickness'};
 
 elseif strcmpi(string,'transient'), 
 	if strcmpi(mdtype,'3d')
-		fields={'transient_results(end).vy','transient_results(end).vz','transient_results(end).temperature','transient_results(end).melting','transient_results(end).thickness'};
+		fields={'transient(end).vy','transient(end).vz','transient(end).temperature','transient(end).melting','transient(end).thickness'};
 	else
-		fields={'transient_results(end).vel','transient_results(end).thickness'};
+		fields={'transient(end).vel','transient(end).thickness'};
 	end
 
 elseif strcmpi(string,'control'),
-	fields={'cont_vel','cont_J','cont_parameter'};
+	fields={'control.vel','control.J','control.parameter'};
 
 else