Index: sm/trunk/src/m/model/BasinConstrain2.m
===================================================================
--- /issm/trunk/src/m/model/BasinConstrain2.m	(revision 9366)
+++ 	(revision )
@@ -1,64 +1,0 @@
-function md=BasinConstrain2(md,domain);
-%BASINCONSTRAIN - constrain basin
-%
-%   Constrain basin using a constraint domain outline, 
-%   to dirichlet boundary conditions.
-%   constraindomain is an Argus domain outline file enclosing 
-%   the geographical area of interest.
-%
-%   Usage: 
-%      md=BasinConstrain(md,constraindomain)
-%
-%   Example:
-%      md=BasinConstrain(md,'DomainOutline.exp');
-%      md=BasinConstrain(md,'~Iceshelves.exp');
-
-%now, flag nodes and elements outside the domain outline.
-if ischar(domain),
-	if isempty(domain),
-		elementondomain=zeros(md.numberofelements,1);
-		nodeondomain=zeros(md.numberofnodes,1);
-		invert=0;
-	elseif strcmpi(domain,'all')
-		elementondomain=ones(md.numberofelements,1);
-		nodeondomain=ones(md.numberofnodes,1);
-		invert=0;
-	else
-		%make sure that we actually don't want the elements outside the domain outline!
-		if strcmpi(domain(1),'~'),
-			domain=domain(2:end);
-			invert=1;
-		else
-			invert=0;
-		end
-		%ok, flag elements and nodes
-		[nodeondomain elementondomain]=ContourToMesh(md.elements(:,1:3),md.x,md.y,domain,'element and node',2);
-	end
-	if invert,
-		nodeondomain=~nodeondomain;
-		elementondomain=~elementondomain;
-	end
-else
-	error('BasinConstrain error message: domain type not supported yet');
-end
-
-%list of elements and nodes not on domain
-nodenotondomain=find(~nodeondomain);
-elementnotondomain=find(~elementondomain);
-
-%all elements outside the constraint domain are equivalent to water. all nodes outside are spc'd.
-md.spcvx(nodenotondomain)=md.vx_obs(nodenotondomain);
-md.spcvy(nodenotondomain)=md.vy_obs(nodenotondomain);
-md.elementonwater(elementnotondomain)=1;
-
-%now, make sure all elements on water have nodes that are spc'd, otherwise, we'll get a singular problem.
-pos=find(~md.elementonwater);
-numpos=unique(md.elements(pos,:));
-nodes=setdiff(1:1:md.numberofnodes,numpos);
-md.spcvx(nodes)=md.vx_obs(nodes);
-md.spcvy(nodes)=md.vy_obs(nodes);
-
-
-%make sure icefronts that are completely spc'd are taken out:
-free_segments=find((~isnan(md.spcvx(md.pressureload(:,1:2))) + ~isnan(md.spcvy(md.pressureload(:,1:2))))~=2);
-md.pressureload=md.pressureload(free_segments,:);
Index: sm/trunk/src/m/model/balvel.m
===================================================================
--- /issm/trunk/src/m/model/balvel.m	(revision 9366)
+++ 	(revision )
@@ -1,26 +1,0 @@
-function [velx,vely,vel]=balvel(md)
-%BALVEL - computation of balanced velocities
-%
-%   This routine uses the model of Hutter to compute the velocities
-%   of a 2d model using the surface slope
-%
-%   Usage:
-%      [velx,vely,vel]=balvel(md)
-
-if md.dim~=2,
-	error('Only 2d meshes are allowed to compute velocity balances');
-end
-
-%Get slope
-[sx,sy,s]=slope(md);
-
-%Average thickness and B over all elements.
-summer=[1;1;1];
-hel=md.thickness(md.elements)*summer/3;
-Bel=md.B(md.elements)*summer/3;
-
-Ael=Bel.^(-3);
-
-velx=-2*(md.rho_ice*md.g)^3*s.^2.*sx.*Ael/4.*hel.^4;
-vely=-2*(md.rho_ice*md.g)^3*s.^2.*sy.*Ael/4.*hel.^4;
-vel=sqrt(velx.^2+vely.^2);
Index: sm/trunk/src/m/model/export.m
===================================================================
--- /issm/trunk/src/m/model/export.m	(revision 9366)
+++ 	(revision )
@@ -1,16 +1,0 @@
-function export_results=export(md,results, xlower,ylower,xposting,yposting,xsize,ysize)
-%EXPORT - exports results from a mesh to a regular grid
-%
-%   this routine exports results from a mesh to a regular grid 
-%   of postings (xposting,yposint) with coordinates for lower left 
-%   corner (xlower,ylower), and of size (xsize,ysize) 
-%
-%   Usage:
-%      export_results=export(md,results, xlower,ylower,xposting,yposting,xsize,ysize)
-
-%Build x_m and y_m
-x_m=(xlower):xposting:(xlower+xposting*(xsize-1));
-y_m=(ylower):yposting:(ylower+yposting*(ysize-1));
-
-[X,Y]=meshgrid(x_m,y_m);
-export_results=griddata(md.x,md.y,results,X,Y);
Index: sm/trunk/src/m/model/plug.m
===================================================================
--- /issm/trunk/src/m/model/plug.m	(revision 9366)
+++ 	(revision )
@@ -1,38 +1,0 @@
-function md=plug(md,field,value)
-%PLUG - plugs the holes of a model field
-%
-%   This routine takes a model field interpolated on a mesh, and plugs the holes (made of value)
-%   to the nearest interpolation
-%
-%   Usage:
-%      md=plug(md,field,value)
-
-%First retrieve field
-fieldvalues=eval(['md.' field]);
-
-if isnan(value),
-	holevalues=find(isnan(fieldvalues));
-elseif isinf(value),
-	holevalues=find(isinf(fieldvalues));
-else
-	holevalues=find(fieldvalues==value);
-end
-
-for i=1:length(holevalues),
-
-	index=holevalues(i);
-
-	if ~mod(i,100),
-		%disp(['Step ' num2str(i) '/' num2str(length(holevalues))]);
-	end
-
-	distance=sqrt( (md.x(index)-md.x).^2 + (md.y(index)-md.y).^2 );
-	distance(holevalues)=10^10; %be sure we are not picking up the nan nodes!
-
-	index2=find(distance==min(distance));index2=index2(1);
-
-	fieldvalues(index)=fieldvalues(index2);
-end
-
-%Set fieldvalues in md.field
-eval(['md.' field '=fieldvalues;']);
Index: sm/trunk/src/m/model/plugvelocitiesraw.m
===================================================================
--- /issm/trunk/src/m/model/plugvelocitiesraw.m	(revision 9366)
+++ 	(revision )
@@ -1,98 +1,0 @@
-function md=plugvelocitiesraw(md,filename,default_value)
-%PLUGVELOCITIESRAW - load raw velocities on a model
-%
-%   load a matlab file (extension .mat) which holds 4 variables
-%   x,y,vx,vy to be plugged onto the model (or similar names)
-%   x and y must be vectors, vx, vy matrices
-%
-%   Usage:
-%      md=plugvelocitiesraw(md,filename,default_value)
-%
-%   Example:
-%      md=plugvelocitiesraw(md,'velocityfile.mat',0);
-%
-%   See also: PLUGVELOCITIES INTERPFROMFILE, GRIDDATA
-
-%some checks
-if nargin~=3 | nargout~=1
-	help plugvelocitiesraw
-	error('plugvelocitiesraw error message: bad usage');
-end
-if ~exist(filename)
-	error(['plugvelocitiesraw error message: file ' filename  ' does not exist']);
-end
-
-%Get variables
-A=whos('-file',filename);
-
-%find x,y,vx and vy
-xenum=NaN; yenum=NaN; vxenum=NaN; vyenum=NaN;
-if length(A)==4,
-	for i=1:4
-		if strcmpi(A(i).name(1),'x');
-			xenum=i;
-		elseif strcmpi(A(i).name(1),'y');
-			yenum=i;
-		else
-			if strcmpi(A(i).name(end),'x');
-				vxenum=i;
-			elseif strcmpi(A(i).name(end),'y');
-				vyenum=i;
-			end
-		end
-	end
-else
-	error(['plugvelocitiesraw error message: file ' filename  ' not supported yet (it should hold 4 variables x,y,vx and vy)']);
-end
-
-%assum that we have found at least vxenum and vyenum
-if ( isnan(vxenum) | isnan(vyenum))
-	error(['plugvelocitiesraw error message: file ' filename  ' not supported yet (it should hold 4 variables x,y,vx and vy)']);
-end
-
-%find x y
-if (isnan(xenum) | isnan(yenum))
-
-	%check the size
-	if A(vxenum).size(1)==A(vxenum).size(2),
-		error(['plugvelocitiesraw error message: file ' filename  ' not supported (velocities is a square matrix, save x and y with another name)']);
-	end
-	if ~(A(vxenum).size(1)==A(vyenum).size(1) & A(vxenum).size(2)==A(vyenum).size(2)),
-		error(['plugvelocitiesraw error message: file ' filename  ' not supported (vx and vy matrices do not have the same size)']);
-	end
-
-	%find xenum and yenum
-	for i=1:4
-		lengthi=max(A(i).size);
-		if ((i~=vxenum) & (lengthi==A(vxenum).size(1) | lengthi==A(vxenum).size(1)+1)),
-			yenum=i;
-		elseif ((i~=vxenum) & (lengthi==A(vxenum).size(2) | lengthi==A(vxenum).size(2)+1)),
-			xenum=i;
-		end
-	end
-
-	%last check
-	if (isnan(xenum) | isnan(yenum))
-		error(['plugdata error message: file ' filename  ' not supported yet']);
-	end
-end
-
-%create names:
-xname=A(xenum).name;
-yname=A(yenum).name;
-vxname=A(vxenum).name;
-vyname=A(vyenum).name;
-
-%load data
-B=load(filename);
-
-%get x y vx and vy
-eval(['x=B.' xname ';'])
-eval(['y=B.' yname ';'])
-eval(['vx=B.' vxname ';'])
-eval(['vy=B.' vyname ';'])
-
-%interpolate
-md.vx_obs_raw=InterpFromGridToMesh(x,y,vx,md.x,md.y,default_value);
-md.vy_obs_raw=InterpFromGridToMesh(x,y,vy,md.x,md.y,default_value);
-md.vel_obs_raw=sqrt(md.vx_obs_raw.^2+md.vy_obs_raw.^2);
Index: sm/trunk/src/m/model/removeholes.m
===================================================================
--- /issm/trunk/src/m/model/removeholes.m	(revision 9366)
+++ 	(revision )
@@ -1,46 +1,0 @@
-function md2=removeholes(md,field)
-%REMOVEHOLES - interpolate a field on the a mesh without any hole
-%
-%   as its name indicates, this routine takes a model, a field (value of some 
-%   physical quantity evaluated at the mesh nodes) of the model, and interpolates this field 
-%   on the model mesh, without any hole.
-%
-%   Usage:
-%      md=removeholes(md,field)
-
-%Check that model is complete
-if md.counter<3,
-	error('removeholes error message: model is incomplete ... exiting');
-end
-
-if nargin~=2,
-	removeholesusage;
-	error('removeholes error message');
-end
-
-if ~ischar(field), 
-	removeholesusage;
-	error('removeholes error message');
-end
-
-%Ok, retrieve and write domain outline without holes, to disk.
-domainoutline_string=md.domainoutline;
-name_slots=findstr(domainoutline_string,'## Name');
-domainoutline_string=domainoutline_string(1:(name_slots(2)-1)); %only keep first outline
-writefile('DomainOutlineTemp.exp',domainoutline_string);
-
-%Now create new model with mesh based on DomainOutlineTemp: 
-%get average resolution
-resolution=mean(sqrt(2*area(md.elements,md.x,md.y)));
-md2=model;
-md2=mesh(md2,'DomainOutlineTemp.exp',resolution);
-
-%Ok, now interpolate field onto this new mesh
-fieldvalue=getfield(md,field);
-md2=setfield(md2,field,griddata(md.x,md.y,fieldvalue,md2.x,md2.y));
-
-end %end of function
-
-function removeholesusage(),
-disp('usage: md2=removeholes(md,field)');
-end
Index: /issm/trunk/src/m/model/sia.m
===================================================================
--- /issm/trunk/src/m/model/sia.m	(revision 9367)
+++ /issm/trunk/src/m/model/sia.m	(revision 9367)
@@ -0,0 +1,26 @@
+function [velx,vely,vel]=sia(md)
+%BALVEL - computation of Shallow Ice velocities
+%
+%   This routine uses the model of Hutter to compute the velocities
+%   of a 2d model using the surface slope
+%
+%   Usage:
+%      [velx,vely,vel]=sia(md)
+
+if md.dim~=2,
+	error('Only 2d meshes are allowed to compute velocity balances');
+end
+
+%Get slope
+[sx,sy,s]=slope(md);
+
+%Average thickness and B over all elements.
+summer=[1;1;1];
+hel=md.thickness(md.elements)*summer/3;
+Bel=md.B(md.elements)*summer/3;
+
+Ael=Bel.^(-3);
+
+velx=-2*(md.rho_ice*md.g)^3*s.^2.*sx.*Ael/4.*hel.^4;
+vely=-2*(md.rho_ice*md.g)^3*s.^2.*sy.*Ael/4.*hel.^4;
+vel=sqrt(velx.^2+vely.^2);