Index: /issm/trunk/src/m/solutions/cielo/UpdateGeometry.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/UpdateGeometry.m	(revision 738)
+++ /issm/trunk/src/m/solutions/cielo/UpdateGeometry.m	(revision 738)
@@ -0,0 +1,27 @@
+function [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,bed,surface)
+%UPDATEGEOMETRY - update the geometry during after a time step
+%
+%   Update geometry once given a new thickness, in the transient solutions
+%
+%   Usage:
+%      [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,bed,surface)
+
+%initialize bed
+new_bed=zeros(md.numberofgrids,1);
+new_surface=zeros(md.numberofgrids,1);
+
+%First, check that thickness is >0
+pos=find(new_thickness<0);
+new_thickness(pos)=10; %minimum thickness
+
+%For grids on ice sheet, the new bed remains the same, the new surface becomes bed+new_thickness. 
+pos=find(md.gridonicesheet);
+new_bed(pos)=bed(pos);
+new_surface(pos)=bed(pos)+new_thickness(pos);
+
+%For grids on ice shelt, we have hydrostatic equilibrium (for now)
+pos=find(md.gridoniceshelf);
+new_bed(pos)=-md.rho_ice/md.rho_water*new_thickness(pos);
+new_surface(pos)=(1-md.rho_ice/md.rho_water)*new_thickness(pos);
+
+%Deal with grounding line migration: later on.
Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 737)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic.m	(revision 738)
@@ -1,4 +1,5 @@
-function md=diagnostic_proto(md);
+function md=diagnostic(md);
 %DIAGNOSTIC - compute the velocity field of a model
+%
 %   Usage:
 %      md=diagnostic(md)
Index: /issm/trunk/src/m/solutions/cielo/prognostic.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 737)
+++ /issm/trunk/src/m/solutions/cielo/prognostic.m	(revision 738)
@@ -1,4 +1,4 @@
 function md=prognostic(md)
-%THERMAL - prognostic solution sequence.
+%PROGNOSITC - prognostic solution sequence.
 %
 %   Usage:
Index: /issm/trunk/src/m/solutions/cielo/transient.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/transient.m	(revision 738)
+++ /issm/trunk/src/m/solutions/cielo/transient.m	(revision 738)
@@ -0,0 +1,18 @@
+function md=transient(md);
+%TRANSIENT - transient solution
+%
+%   Usage:
+%      md=transient(md)
+
+%start timing
+t1=clock;
+
+if strcmpi(md.type,'2d'),
+	md=transient2d(md);
+else
+	md=transient3d(md);
+end
+
+%stop timing
+t2=clock;
+disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']));
Index: /issm/trunk/src/m/solutions/cielo/transient2d.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/transient2d.m	(revision 738)
+++ /issm/trunk/src/m/solutions/cielo/transient2d.m	(revision 738)
@@ -0,0 +1,111 @@
+function md=transient2d(md);
+%TRANSIENT2D - 2d transient solution
+%
+%   Usage:
+%      md=transient(md)
+%
+%   See also: TRANSIENT3D, TRANSIENT
+
+%Build all models requested
+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);
+
+%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=[];
+solution.h_g=m_p.parameters.h_g;
+solution.s_g=md.surface;
+solution.b_g=md.bed;
+
+%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',m_p.parameters.m_g,'doublevec',1,m_p.parameters.numberofnodes);
+inputs=add(inputs,'accumulation',m_p.parameters.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=md.dt;
+finaltime=md.ndt;
+time=dt;
+n=1; %counter
+
+while  time<finaltime+dt, %make sure we run up to finaltime.
+
+	disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time/md.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);
+
+	%Deal with velocities.
+
+	%Get horizontal solution. 
+	[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
+	disp(sprintf('%s','   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
+	disp(sprintf('%s','   updating geometry...'));
+	[new_bed,new_surface,new_thickness]=UpdateGeometry(md,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
+	%disp(sprintf('%s','   checking time stepping...'));
+	%[back,dt,time]=TimeStepping(md,solution,dt,time);
+	%if back,
+	%	continue;
+	%end
+
+	%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;
+	results(i).vx=solution(i).u_g(1:2:end);
+	results(i).vy=solution(i).u_g(2:2:end);
+	results(i).vel=sqrt(solution(i).u_g(1:2:end).^2+solution(i).u_g(2:2:end).^2)*m_dh.parameters.yts;
+	results(i).bed=solution(i).b_g;
+	results(i).surface=solution(i).s_g;
+	results(i).thickness=solution(i).h_g;
+end
+md.results.transient=results;