Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 848)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 849)
@@ -21,6 +21,6 @@
 	if dim==3,
 		displaystring(debug,'\n%s',['extruding slopes in 3d...']);
-		slopex=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex);
-		slopey=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey);
+		slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex,"slopex",0);
+		slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey,"slopey",0);
 	end
 
@@ -55,5 +55,5 @@
 
 	displaystring(debug,'\n%s',['extruding horizontal velocities...']);
-	u_g_horiz=VelocityExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g);
+	u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g,"velocity",1);
 
 	displaystring(debug,'\n%s',['computing vertical velocities...']);
@@ -77,6 +77,6 @@
 		slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedx');
 		slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedy');
-		slopex=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex);
-		slopey=SlopeExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey);
+		slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex,"slopex",0);
+		slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey,"slopey",0);
 
 		inputs=add(inputs,'bedslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
Index: /issm/trunk/src/m/solutions/cielo/prognostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic_core.m	(revision 848)
+++ /issm/trunk/src/m/solutions/cielo/prognostic_core.m	(revision 849)
@@ -8,5 +8,5 @@
 	%Take only the first two dofs of m.parameters.u_g
 	u_g=get(inputs,'velocity',[1 1 0 0]);
-	u_g=VelocityDepthAverage(m.elements,m.nodes,m.loads,m.materials,u_g);
+	u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,u_g,"velocity");
 	inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes);
 
@@ -15,5 +15,5 @@
 
 	displaystring(m.parameters.debug,'\n%s',['extrude computed thickness on all layers:']);
-	h_g=ThicknessExtrude(m.elements,m.nodes,m.loads,m.materials,h_g);
+	h_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,h_g,"thickness",0);
 
 end %end function
Index: /issm/trunk/src/m/solutions/cielo/thermal_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 848)
+++ /issm/trunk/src/m/solutions/cielo/thermal_core.m	(revision 849)
@@ -45,5 +45,5 @@
 
 		disp('   computing melting...');
-		inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
+		inputs=add(inputs,'temperature',soln(n+1).t_g,'doublevec',1,m_t.parameters.numberofnodes);
 		inputs=add(inputs,'melting_offset',melting_offset,'double');
 		soln(n+1).m_g=diagnostic_core_linear(m_m,inputs,'melting','transient');
Index: /issm/trunk/src/m/solutions/cielo/transient3d.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/transient3d.m	(revision 849)
+++ /issm/trunk/src/m/solutions/cielo/transient3d.m	(revision 849)
@@ -0,0 +1,147 @@
+function md=transient3d(md);
+%TRANSIENT2D - 3d transient solution
+%
+%   Usage:
+%      md=transient3d(md)
+%
+%   See also: TRANSIENT2D, TRANSIENT
+
+%Build all models related to diagnostic
+displaystring(md.debug,'%s',['reading diagnostic horiz model data']);
+md.analysis_type='diagnostic'; md.sub_analysis_type='horiz'; m_dh=CreateFemModel(md);
+
+displaystring(md.debug,'\n%s',['reading diagnostic vert model data']);
+md.analysis_type='diagnostic'; md.sub_analysis_type='vert'; m_dv=CreateFemModel(md);
+
+displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']);
+md.analysis_type='diagnostic'; md.sub_analysis_type='stokes'; m_ds=CreateFemModel(md);
+
+displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']);
+md.analysis_type='diagnostic'; md.sub_analysis_type='hutter'; m_dhu=CreateFemModel(md);
+
+displaystring(md.debug,'\n%s',['reading surface and bed slope computation model data']);
+md.analysis_type='slope_compute'; md.sub_analysis_type=''; m_sl=CreateFemModel(md);
+
+displaystring(md.debug,'%s',['reading prognostic model data']);
+md.analysis_type='prognostic'; m_p=CreateFemModel(md);
+
+%Build all models related to thermal
+displaystring(md.debug,'%s',['reading thermal model data']);
+md.analysis_type='thermal'; m_t=CreateFemModel(md);
+
+displaystring(md.debug,'%s',['reading melting model data']);
+md.analysis_type='melting'; m_m=CreateFemModel(md);
+
+
+%initialize solution
+solution=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
+solution.step=1;
+solution.time=0;
+solution.u_g=m_dh.parameters.u_g(dofsetgen([1,2],3,length(m_dh.parameters.u_g)));
+solution.p_g=m_p.parameters.p_g;
+solution.h_g=m_p.parameters.h_g;
+solution.s_g=m_p.parameters.s_g;
+solution.b_g=m_p.parameters.b_g;
+solution.t_g=m_p.parameters.t_g;
+solution.m_g=m_p.parameters.m_g;
+
+%initialize inputs
+displaystring(md.debug,'\n%s',['setup inputs...']);
+inputs=inputlist;
+inputs=add(inputs,'velocity',solution.u_g,'doublevec',2,m_p.parameters.numberofnodes);
+inputs=add(inputs,'melting',solution.m_g,'doublevec',1,m_p.parameters.numberofnodes);
+inputs=add(inputs,'accumulation',solution.a_g,'doublevec',1,m_p.parameters.numberofnodes);
+inputs=add(inputs,'dt',m_p.parameters.dt,'double');
+
+% figure out number of dof: just for information purposes.
+md.dof=modelsize(m_dh,m_dv,m_ds,m_dhu,m_sl);
+
+%first time step is given by model. 
+dt=m_p.parameters.dt;
+finaltime=m_p.parameters.ndt;
+time=dt;
+yts=m_p.parameters.yts;
+n=1; %counter
+
+while  time<finaltime+dt, %make sure we run up to finaltime.
+
+	displaystring(md.debug,'\n%s%g%s%g%s%g\n','time [yr]: ',time/yts,'    iteration number: ',n,'/',floor(finaltime/dt));
+
+	solution(n+1).step=n+1;
+	solution(n+1).time=time;
+
+	%update inputs
+	inputs=add(inputs,'thickness',solution(n).h_g,'doublevec',1,m_p.parameters.numberofnodes);
+	inputs=add(inputs,'surface',solution(n).s_g,'doublevec',1,m_p.parameters.numberofnodes);
+	inputs=add(inputs,'bed',solution(n).b_g,'doublevec',1,m_p.parameters.numberofnodes);
+	inputs=add(inputs,'velocity',solution(n).u_g,'doublevec',m_p.parameters.numberofdofspernodes,m_p.parameters.numberofnodes);
+	inputs=add(inputs,'pressure',solution(n).p_g,'doublevec',1,m_p.parameters.numberofnodes);
+	inputs=add(inputs,'temperature',solution(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
+
+	%Deal with temperature first 
+	displaystring(md.debug,'\n%s',['    computing temperatures...']);
+	[solution(n+1).t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,'thermal','transient');
+	
+	displaystring(md.debug,'\n%s',['    computing melting...']);
+	inputs=add(inputs,'melting_offset',melting_offset,'double');
+	solution(n+1).m_g=diagnostic_core_linear(m_m,inputs,'melting','transient');
+
+	%Compute depth averaged temperature
+	temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,solution(n+1).t_g,'temperature');
+	inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,m_t.parameters.numberofnodes);
+
+	%Deal with velocities.
+	[solution(n+1).u_g,solution(n+1).p_g]=diagnostic_core(m_dh,m_dhu,m_dv,m_ds,m_sl,inputs);
+
+	%compute new thickness
+	displaystring(md.debug,'   computing new thickness...'));
+	inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',2,m_p.parameters.numberofnodes);
+	new_thickness=prognostic_core(m_p,inputs,'prognostic','');
+
+	%update surface and bed using the new thickness
+	displaystring(md.debug,'   updating geometry...'));
+	[new_thickness,new_bed,new_surface]=UpdateGeometry(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,m_p.parameters,new_thickness,solution(n).b_g,solution(n).s_g);
+
+	%Record bed surface and thickness in the solution
+	solution(n+1).h_g=new_thickness;
+	solution(n+1).b_g=new_bed;
+	solution(n+1).s_g=new_surface;
+
+	%figure out if time stepping is good
+	%displaystring(md.debug,'   checking time stepping...'));
+	%[back,dt,time]=TimeStepping(md,solution,dt,time);
+	%if back,
+	%	continue;
+	%end
+
+	%update nodes
+	displaystring(md.debug,'   updating node positions...'));
+	nodes_dh=UpdateNodePositions(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,new_bed,new_thickness);
+	nodes_dv=UpdateNodePositions(m_dv.elements,m_dv.nodes,m_dv.loads,m_dv.materials,new_bed,new_thickness);
+	nodes_ds=UpdateNodePositions(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,new_bed,new_thickness);
+	nodes_sl=UpdateNodePositions(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,new_bed,new_thickness);
+	nodes_p=UpdateNodePositions(m_p.elements,m_p.nodes,m_p.loads,m_p.materials,new_bed,new_thickness);
+	nodes_t=UpdateNodePositions(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,new_bed,new_thickness);
+	nodes_m=UpdateNodePositions(m_m.elements,m_m.nodes,m_m.loads,m_m.materials,new_bed,new_mhickness);
+
+	%update time and counter
+	time=time+dt;
+	n=n+1;
+end
+
+%Load results onto model
+results=struct('step',[],'time',[],'vx',[],'vy',[],'vel',[],'pressure',[],'thickness',[],'surface',[],'bed',[]);
+for i=1:length(solution),
+	results(i).step=solution(i).step;
+	results(i).time=solution(i).time/yts;
+	results(i).vx=solution(i).u_g(1:2:end)*yts;
+	results(i).vy=solution(i).u_g(2:2:end)*yts;
+	results(i).vel=sqrt(solution(i).u_g(1:2:end).^2+solution(i).u_g(2:2:end).^2)*yts;
+	results(i).bed=solution(i).b_g;
+	results(i).surface=solution(i).s_g;
+	results(i).thickness=solution(i).h_g;
+	results(i).temperature=solution(i).t_g;
+	results(i).melting=solution(i).m_g;
+end
+if isnan(md.results), md.results=struct(); end
+md.results.transient=results;