Index: /issm/trunk/src/m/model/plot/plot_BC.m
===================================================================
--- /issm/trunk/src/m/model/plot/plot_BC.m	(revision 10128)
+++ /issm/trunk/src/m/model/plot/plot_BC.m	(revision 10129)
@@ -7,7 +7,7 @@
 
 %plot dirichlets
-h1=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvx))),md.mesh.y(find(~isnan(md.diagnostic.spcvx))),'ro','MarkerSize',14,'MarkerFaceColor','r');
-h2=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvy))),md.mesh.y(find(~isnan(md.diagnostic.spcvy))),'bo','MarkerSize',10,'MarkerFaceColor','b');
-h3=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvz))),md.mesh.y(find(~isnan(md.diagnostic.spcvz))),'yo','MarkerSize',6 ,'MarkerFaceColor','y');
+h1=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvx(1:md.mesh.numberofvertices,1)))),md.mesh.y(find(~isnan(md.diagnostic.spcvx(1:md.mesh.numberofvertices,1)))),'ro','MarkerSize',14,'MarkerFaceColor','r');
+h2=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvy(1:md.mesh.numberofvertices,1)))),md.mesh.y(find(~isnan(md.diagnostic.spcvy(1:md.mesh.numberofvertices,1)))),'bo','MarkerSize',10,'MarkerFaceColor','b');
+h3=plot(md.mesh.x(find(~isnan(md.diagnostic.spcvz(1:md.mesh.numberofvertices,1)))),md.mesh.y(find(~isnan(md.diagnostic.spcvz(1:md.mesh.numberofvertices,1)))),'yo','MarkerSize',6 ,'MarkerFaceColor','y');
 
 %update legend
Index: /issm/trunk/src/m/model/plot/plot_transient_movie.m
===================================================================
--- /issm/trunk/src/m/model/plot/plot_transient_movie.m	(revision 10128)
+++ /issm/trunk/src/m/model/plot/plot_transient_movie.m	(revision 10129)
@@ -16,29 +16,40 @@
 	end
 
+	results=md.results.TransientSolution;
 	%loop over the time steps
+	if exist(options,'transient_movie_limit'),
+		limit=getfieldvalue(options,'transient_movie_limit');
+		steps=[limit(1):limit(end)];
+	else
+		steps=1:length(results);
+	end
 	results=md.results.TransientSolution;
-	for i=1:length(results)
+	nstep=1;
+	for i=steps
 
-		%process data
-		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
-		[data datatype]=processdata(md,results(i).(field),options);
+		if ~isempty(results(i).(field)),
+			%process data
+			[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
+			[data datatype]=processdata(md,results(i).(field),options);
 
-		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
-		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
-		apply_options_movie(md,options,titlestring);
-		
-		if exist(options,'transient_movie_output'),
-			set(gcf,'Renderer','zbuffer','color','white'); %fixes a bug on Mac OS X (not needed in future Matlab version)
-			if i==1,
-				%initialize images and frame
-				frame=getframe(gcf);
-				[images,map]=rgb2ind(frame.cdata,256,'nodither');
-				images(1,1,1,length(results))=0;
+			titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
+			plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
+			apply_options_movie(md,options,titlestring);
+
+			if exist(options,'transient_movie_output'),
+				set(gcf,'Renderer','zbuffer','color','white'); %fixes a bug on Mac OS X (not needed in future Matlab version)
+				if nstep==1,
+					%initialize images and frame
+					frame=getframe(gcf);
+					[images,map]=rgb2ind(frame.cdata,256,'nodither');
+					images(1,1,1,length(steps))=0;
+				else
+					frame=getframe(gcf);
+					images(:,:,1,nstep) = rgb2ind(frame.cdata,map,'nodither');
+				end
 			else
-				frame=getframe(gcf);
-				images(:,:,1,i) = rgb2ind(frame.cdata,map,'nodither');
+				pause(0.1)
 			end
-		else
-			pause(0.1)
+			nstep=nstep+1;
 		end
 	end
Index: /issm/trunk/src/m/model/plot/plotmodel.m
===================================================================
--- /issm/trunk/src/m/model/plot/plotmodel.m	(revision 10128)
+++ /issm/trunk/src/m/model/plot/plotmodel.m	(revision 10129)
@@ -38,5 +38,5 @@
 		set(F,'Visible','Off');
 	else
-		figure(figurenumber),clf;
+		figure(figurenumber);clf;
 	end
 
Index: /issm/trunk/src/m/model/regionaltransient2d.m
===================================================================
--- /issm/trunk/src/m/model/regionaltransient2d.m	(revision 10129)
+++ /issm/trunk/src/m/model/regionaltransient2d.m	(revision 10129)
@@ -0,0 +1,149 @@
+function md2=regionaltransient2d(md1,area,hmin,hmax,err,stepres)
+%regionaltransient2d - extract a model according to an Argus contour or flag list and remesh
+%               at new resolution res
+%
+%   This routine extracts a submodel from a bigger model with respect to a given contour
+%   md must be followed by the corresponding exp domain file (argus type, .exp extension). 
+%   The model will be remeshed at high rsolution hmin and low resolution hmax.  The ice 
+%   boundary velocities will be spc'd to the transient velocities at saved transient steps
+%   at the resolution optionally provided for stepres.  A stepres of 2 means that you wish
+%   to skip every other saved transient step.  This is useful when extracting a long transient.
+%
+%   Usage:
+%      md2=regionaltransient2d(md1,area,hmin,hmax,err);
+%
+%   Examples:
+%      md2=regionaltransient2d(md,'Domain.exp',500,10000,[15 250]);
+%      md2=regionaltransient2d(md,'Domain.exp',3000,15000,[10 300],2);
+%
+%   See also: MODELEXTRACT, EXTRUDE, COLLAPSE
+
+%some checks
+if ((nargin~=5) & (nargin~=6)),
+	help regionaltransient2d 
+	error('regionaltransient2d error message: bad usage');
+end
+
+%get check option
+if (nargin==5),
+	stepres=1;
+end
+
+%take every fields from model
+mde=modelextract(md1,area);
+mde.private.bamg=[];
+mde.mesh.extractedvertices=nan;
+mde.mesh.extractedelements=nan;
+
+%remesh
+%md2=bamg(md2,'hmax',hmax,'hmin',hmin,'field',[md2.inversion.vel_obs md2.geometry.thickness],'splitcorner',1,'KeepVertices',0);
+md2=bamg(mde,'hmin',hmin,'hmax',hmax,'field',[mde.inversion.vel_obs mde.geometry.thickness],'splitcorner',1,'KeepVertices',0,'err',err);
+md2=setmask(md2,'','');
+
+%automatically modify fields
+
+	%loop over model fields
+	model_fields=fields(md1);
+	for i=1:length(model_fields),
+
+		%get field
+		field=md1.(model_fields{i});
+		fieldsize=size(field);
+
+		if isobject(field), %recursive call
+			object_fields=fields(md1.(model_fields{i}));
+			fname=['(model_fields{i}).(object_fields{j})'];
+		else
+			object_fields=field;
+			fname=['(model_fields{i})'];
+		end
+		for j=1:length(object_fields),
+			%get field
+			field=eval(['md2.' fname]);
+			fieldsize=size(field);
+
+			%size = number of nodes * n
+			for n=1:fieldsize(2)
+				if fieldsize(1)==mde.mesh.numberofvertices
+					if(sum(field(:,n) ~= field(1,n)) == 0)
+						eval(['md2.' fname '(1:md2.mesh.numberofvertices,n)=field(1,n)*ones(md2.mesh.numberofvertices,1);']);
+					else
+						eval(['md2.' fname '(1:md2.mesh.numberofvertices,n)=InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(:,n),md2.mesh.x,md2.mesh.y);']);
+					end
+					eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofvertices,n);']);
+				elseif fieldsize(1)==mde.mesh.numberofvertices+1
+					if(sum(field(1:end-1,n) ~= field(1,n)) == 0)
+						eval(['md2.' fname '(1:md2.mesh.numberofvertices+1,n)=[field(1,n)*ones(md2.mesh.numberofvertices,1); field(end,n)];']);
+					else
+						eval(['md2.' fname '(1:md2.mesh.numberofvertices+1,n)=[InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(1:end-1,n),md2.mesh.x,md2.mesh.y); field(end,n)];']);
+					end
+					eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofvertices+1,n)']);
+					%size = number of elements * n
+				elseif fieldsize(1)==mde.mesh.numberofelements
+					if(sum(field(1:end-1,n) ~= field(1,n)) == 0)
+						eval(['md2.' fname '(1:md2.mesh.numberofelements,n)=field(1,n)*ones(md2.mesh.numberofelements,1);']);
+					else
+						eval(['md2.' fname '(1:md2.mesh.numberofelements,n)=InterpFromMeshToMesh2d(mde.mesh.elements,mde.mesh.x,mde.mesh.y,field(:,n),md2.mesh.x,md2.mesh.y);']);
+					end
+					eval(['md2.' fname '(:,n)=md2.' fname '(1:md2.mesh.numberofelements,n);']);
+				end
+			end
+		end
+	end
+
+	spcx=[];
+	spcy=[];
+	steps=[];
+	nsteps=length(md1.results.TransientSolution);
+	count=0;
+	numElements=arrayfun(@(x) numel(x.step), md1.results.TransientSolution);
+	for t=find(numElements==1)
+		if ~isempty(md1.results.TransientSolution(t).Vel) & mod(count,stepres)==0,
+			vx=PatchToVec(md1.results.TransientSolution(t).Vx);
+			vy=PatchToVec(md1.results.TransientSolution(t).Vy);
+			spcx=[spcx InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vx,md2.mesh.x,md2.mesh.y)];
+			spcy=[spcy InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,vx,md2.mesh.x,md2.mesh.y)];
+			steps=[steps t*md1.timestepping.time_step];
+		end
+		count=count+1;
+	end
+
+	if nsteps > 0
+		md2.diagnostic.spcvx=md2.diagnostic.spcvx*ones(1,size(spcx,2));
+		md2.diagnostic.spcvy=md2.diagnostic.spcvy*ones(1,size(spcy,2));
+		md2.diagnostic.spcvz=md2.diagnostic.spcvz*ones(1,size(spcx,2));
+		md2.diagnostic.spcvx(find(md2.mesh.vertexonboundary),:)=spcx(find(md2.mesh.vertexonboundary),:);
+		md2.diagnostic.spcvy(find(md2.mesh.vertexonboundary),:)=spcy(find(md2.mesh.vertexonboundary),:);
+		md2.diagnostic.spcvz(find(md2.mesh.vertexonboundary),:)=0;
+		md2.diagnostic.spcvx=[md2.diagnostic.spcvx; steps];
+		md2.diagnostic.spcvy=[md2.diagnostic.spcvy; steps];
+		md2.diagnostic.spcvz=[md2.diagnostic.spcvz; steps];
+	end
+
+	%Diagnostic
+	if ~isnan(md2.diagnostic.icefront)
+		md1s=modelextract(md1,area);
+		%md2.diagnostic.icefront=[md2.mesh.segments 2];
+		e2=md2.mesh.segments(:,end);
+		e1=md1s.mesh.segments(:,end);
+
+		pload = nan*ones(size(md1s.mesh.elements,1),1);
+		pload(md1s.diagnostic.icefront(:,end-1))=md1s.diagnostic.icefront(:,end);
+
+		x2=mean(md2.mesh.x(md2.mesh.elements(e2,:)),2);
+      y2=mean(md2.mesh.y(md2.mesh.elements(e2,:)),2);
+		x1=mean(md1s.mesh.x(md1s.mesh.elements),2);
+      y1=mean(md1s.mesh.y(md1s.mesh.elements),2);
+
+		pload2=griddata(x1,y1,pload,x2,y2,'nearest');
+		md2.diagnostic.icefront=[md2.mesh.segments(~isnan(pload2),:) pload2(~isnan(pload2))];
+		md2.diagnostic.spcvx(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
+		md2.diagnostic.spcvy(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
+		md2.diagnostic.spcvz(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
+	end
+
+	%Clear results fields
+	if isstruct(md1.results),
+		md2.results=[];
+	end
+