Index: /issm/trunk-jpl/src/m/contrib/tsantos/remesh.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/tsantos/remesh.m	(revision 24514)
+++ /issm/trunk-jpl/src/m/contrib/tsantos/remesh.m	(revision 24515)
@@ -1,4 +1,7 @@
-function mdOut = remesh(md,parfile)
+function mdOut = remesh(md,parfile,time_step)
 %Set the new mesh (refined) into the model md
+if nargin<3,
+	time_step=length(md.results.TransientSolution);
+end
 
 NewModel = model;
@@ -6,7 +9,7 @@
 % geometry
 NewModel.mesh				= mesh2d();
-NewModel.mesh.x			= md.results.TransientSolution(end).MeshX;
-NewModel.mesh.y			= md.results.TransientSolution(end).MeshY;
-NewModel.mesh.elements	= md.results.TransientSolution(end).MeshElements;
+NewModel.mesh.x			= md.results.TransientSolution(time_step).MeshX;
+NewModel.mesh.y			= md.results.TransientSolution(time_step).MeshY;
+NewModel.mesh.elements	= md.results.TransientSolution(time_step).MeshElements;
 
 % build segments. CONVEX HULL: IT JUST WORKS FOR REGULAR MESHES, WITHOUT "BAYS"
@@ -28,5 +31,5 @@
 
 	if segments(s,3)==0,
-		error('Element not found!');
+		%error('Element not found!');
 	end
 end
@@ -51,14 +54,14 @@
 
 % Setting initialization
-NewModel.initialization.vx				= md.results.TransientSolution(end).Vx;
-NewModel.initialization.vy				= md.results.TransientSolution(end).Vy;
+NewModel.initialization.vx				= md.results.TransientSolution(time_step).Vx;
+NewModel.initialization.vy				= md.results.TransientSolution(time_step).Vy;
 NewModel.initialization.vz				= zeros(md.mesh.numberofvertices,1);
-NewModel.initialization.vel			= md.results.TransientSolution(end).Vel;
-NewModel.initialization.pressure    = md.results.TransientSolution(end).Pressure;
-NewModel.geometry.surface				= md.results.TransientSolution(end).Surface;
-NewModel.geometry.base					= md.results.TransientSolution(end).Base;
-NewModel.geometry.bed					= md.results.TransientSolution(end).Bed;%md.geometry.bed; %use from parameterize
-NewModel.geometry.thickness			= md.results.TransientSolution(end).Thickness;
-NewModel.mask.groundedice_levelset  = md.results.TransientSolution(end).MaskGroundediceLevelset;
+NewModel.initialization.vel			= md.results.TransientSolution(time_step).Vel;
+NewModel.initialization.pressure    = md.results.TransientSolution(time_step).Pressure;
+NewModel.geometry.surface				= md.results.TransientSolution(time_step).Surface;
+NewModel.geometry.base					= md.results.TransientSolution(time_step).Base;
+NewModel.geometry.bed					= md.results.TransientSolution(time_step).Bed;%md.geometry.bed; %use from parameterize
+NewModel.geometry.thickness			= md.results.TransientSolution(time_step).Thickness;
+NewModel.mask.groundedice_levelset  = md.results.TransientSolution(time_step).MaskGroundediceLevelset;
     
 %copy other data