Index: issm/trunk-jpl/src/m/classes/modellist.m
===================================================================
--- issm/trunk-jpl/src/m/classes/modellist.m	(revision 13007)
+++ issm/trunk-jpl/src/m/classes/modellist.m	(revision 13008)
@@ -10,4 +10,116 @@
 	end
 	methods
+		function md_list=modelsextract(md,flags,minel,varargin) % {{{
+			%modelsextract - extract several self contained models according to a list of element flags.
+			%
+			%   The difference between this routine and the modelextract.m routine (without an 's') is that 
+			%   as many models are extracted as there are closed contours defined in area. 
+			%   This routine is needed for example when doing data assimilation of ice shelves in Antarctica. 
+			%   Many independent ice shelves are present, and we don't want data assimilation on one ice shelf 
+			%   to be hindered by another totally independent ice shelf.
+			%
+			%   Usage:
+			%      md_list=modelsextract(md,elementfalgs,minel);
+			%
+			%   Examples:
+			%      md_list=modelsextract(md,md.mask.elementonfloatingice,1000);
+			%
+			%   See also: EXTRUDE, COLLAPSE, MODELEXTRACT
+
+			disp('selecting pools of elements');
+			%go through flags and build as many independent element flags as there are groups of connected 1s
+			%in flags.
+
+			%2D or 3D?
+			if md.mesh.dimension==3,
+				numberofelements=md.mesh.numberofelements2d; %this will be forgotten when we get out.
+				flags=project2d(md,flags,1);
+			else
+				numberofelements=md.mesh.numberofelements;
+			end
+
+			%recover extra arguments: 
+			distance=0;
+			if nargin==4,
+				distance=varargin{1};
+			end
+
+			flag_list=cell(0,1);
+
+			for i=1:size(flags,1),
+
+				if (flags(i)),
+
+					%ok, we are sure element i is part of a new pool.
+					pool=zeros(numberofelements,1);
+					pool=PropagateFlagsFromConnectivity(md.mesh.elementconnectivity,pool,i,flags);
+					flag_list{end+1,1}=pool;
+
+					%speed up rest of computation by taking pool out of flags: 
+					pos=find(pool);flags(pos)=0;
+
+				end
+			end
+
+			%go through flag_list and discard any pool of less than minel elements: 
+			ex_pos=[];
+			for i=1:length(flag_list),
+				if length(find(flag_list{i}))<minel,
+					ex_pos=[ex_pos; i];
+				end
+			end
+			flag_list(ex_pos)=[];
+
+			%now, if distance was specified, expand the flag_list by distance km: 
+			if distance,
+				for i=1:length(flag_list),
+					flag_list{i}=PropagateFlagsUntilDistance(md,flag_list{i},distance);
+				end
+			end
+
+			%now, go use the pools of flags to extract models: 
+			disp(['extracting ' num2str(size(flag_list,1)) ' models']);
+			models=cell(0,1);
+
+			for i=1:size(flag_list,1),
+				disp(['   ' num2str(i) '/' num2str(size(flag_list,1))]);
+				if md.mesh.dimension==3,
+					flags2d=flag_list{i};
+					realflags=project3d(md,flags2d,'element');
+				else
+					realflags=flag_list{i};
+				end
+				models{end+1,1}=modelextract(md,realflags);
+			end
+
+			%return model list
+			md_list=modellist(models);
+
+		end %end of this function }}}
+		function md_list=modelsextractfromdomains(md,directory) % {{{
+			%modelsextractfromdomains- extract several self contained models according to a list of domains
+			%
+			%   Usage:
+			%      md_list=modelsextractfromdomains(md,'Basins/');
+			%
+			%   Examples:
+			%      md_list=modelsextract(md,'Basins/');
+			%
+			%   See also: MODELSEXTRACTS, MODELEXTRACT
+
+			%go into directory and get list of files.
+			cd(directory);
+			basins=listfiles;
+			cd ..
+
+			models=cell(0,1);
+			for i=1:length(basins),
+				models{end+1,1}=modelextract(md,[directory '/' basins{i}]);
+			end
+
+			%return model list: 
+			md_list=modellist(models);
+
+		end % }}}
 		function obj = modellist(varargin) % {{{
 
Index: issm/trunk-jpl/src/m/contrib/ecco/MeltingGroundingLines.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/ecco/MeltingGroundingLines.m	(revision 13008)
+++ issm/trunk-jpl/src/m/contrib/ecco/MeltingGroundingLines.m	(revision 13008)
@@ -0,0 +1,26 @@
+function md=MeltingGroundingLines(md,distance,value)
+%MELTINGGROUNDINGLINES - set melting near grounding lines to a constant value
+%
+%   Usage:
+%      md=MeltingGroundingLines(md,distance,value)
+%
+
+%get nodes on ice sheet and on ice shelf
+pos_shelf=find(~md.mask.vertexongroundedice);
+pos_GL=intersect(unique(md.mesh.elements(find(md.mask.elementongroundedice),:)),unique(md.mesh.elements(find(md.mask.elementonfloatingice),:)));
+
+for i=1:length(pos_shelf)
+
+	if (mod(i,100)==0),
+		fprintf('\b\b\b\b\b\b\b%5.2f%s',i/length(pos_shelf)*100,' %');
+	end
+
+	%search the node on ice sheet the closest to i
+	[d posd]=min(sqrt((md.mesh.x(pos_shelf(i))-md.mesh.x(pos_GL)).^2+(md.mesh.y(pos_shelf(i))-md.mesh.y(pos_GL)).^2));
+
+	if d<distance,
+
+		md.melting(pos_shelf(i))=value;
+
+	end
+end
Index: issm/trunk-jpl/src/m/contrib/ecco/PropagateFlagsUntilDistance.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/ecco/PropagateFlagsUntilDistance.m	(revision 13008)
+++ issm/trunk-jpl/src/m/contrib/ecco/PropagateFlagsUntilDistance.m	(revision 13008)
@@ -0,0 +1,60 @@
+function new_flags=PropagateFlagsUntilDistance(md,flags,distance)
+%PROPAGATEFLAGSUNTILDISTANCE
+%
+% Usage: 
+%              flags=PropagateFlagsUntilDistance(md,flags,distance)
+%
+%
+
+new_flags=flags;
+
+%make 3d work in 2d: 
+if md.mesh.dimension==3,
+	md.mesh.x=md.mesh.x2d;
+	md.mesh.y=md.mesh.y2d;
+	md.mesh.elements=md.mesh.elements2d;
+end
+
+%find elements that are at the border of flags: 
+flag_elements=find(flags);
+conn=md.mesh.elementconnectivity(flag_elements,:);
+pos=find(conn);conn(pos)=~flags(conn(pos));
+sum_conn=sum(conn,2);
+border_elements=flag_elements(find(sum_conn>=1));
+
+%average x and y over elements: 
+x_elem=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
+y_elem=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
+
+while 1,
+
+	%keep copy of new_flags for this loop: 
+	new_flags_bak=new_flags;
+
+	%extend new flags by connectivity
+	pos=find(new_flags);
+
+	connected_elements=md.mesh.elementconnectivity(pos,:);
+	connected_elements=connected_elements(find(connected_elements));
+	new_flags(connected_elements)=1;
+
+	%get new elements: 
+	new_elements=find(new_flags & ~new_flags_bak);
+	if ~length(new_elements),
+		%we are done!
+		break;
+	end
+
+	%check which of these new elements are more than distance away from the border elements
+	for i=1:length(new_elements),
+		dist=sqrt(     (x_elem(border_elements)-x_elem(new_elements(i))).^2 + (y_elem(border_elements)-y_elem(new_elements(i))).^2)-distance;
+		if ~any(dist<0)
+			%none of the border elements are within distance, this element is outside out area of interest.
+			%ensure this element never gets found again in the connectivity.
+			pos=find(md.mesh.elementconnectivity==new_elements(i));
+			md.mesh.elementconnectivity(pos)=0;
+			%exclude this new element from the new_flags!
+			new_flags(new_elements(i))=0;
+		end
+	end
+end
Index: issm/trunk-jpl/src/m/contrib/hack/thicknessevolution.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/hack/thicknessevolution.m	(revision 13008)
+++ issm/trunk-jpl/src/m/contrib/hack/thicknessevolution.m	(revision 13008)
@@ -0,0 +1,28 @@
+function dhdt=thicknessevolution(md)
+%THICKNESSEVOLUTION - compute the new thickness of a model after ∆t
+%
+%   This routine compute the new thickness of a model after a time step
+%   according to the following formula:
+%   dh/dt=-div(Hu)
+%
+%   Usage:
+%      dhdt=thicknessevolution(md)
+
+if (length(md.initialization.vx)~=md.mesh.numberofvertices)|(length(md.initialization.vy)~=md.mesh.numberofvertices)
+	error(['thicknessevolution error message: vx and vy should have a length of ' num2str(md.mesh.numberofvertices)])
+end
+
+%load some variables 
+H=md.geometry.thickness;
+vx=md.initialization.vx;
+vy=md.initialization.vy;
+index=md.mesh.elements;
+
+%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
+[alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); 
+
+%compute dhdt=div(Hu)
+summation=1/3*ones(3,1);
+dhdt=(vx(index)*summation).*sum( H(index).*alpha,2) + (vy(index)*summation).*sum(H(index).*beta,2) ...
+	+ ( H(index)*summation).*sum(vx(index).*alpha,2) + ( H(index)*summation).*sum(vy(index).*beta,2);
+dhdt=-dhdt;
Index: issm/trunk-jpl/src/m/contrib/hack/tres.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/hack/tres.m	(revision 13008)
+++ issm/trunk-jpl/src/m/contrib/hack/tres.m	(revision 13008)
@@ -0,0 +1,119 @@
+function md=tres(md,string)
+%TRES - transfer results results to corresponding model fields. 
+%
+%    Usage: md=tres(md,string)
+%
+%    Example: md=tres(md,'diagnostic');
+
+%check number of arguments
+
+if strcmpi(string,'diagnostic'),
+	if md.mesh.dimension==2,
+		md.initialization.vx=md.results.DiagnosticSolution.Vx;
+		md.initialization.vy=md.results.DiagnosticSolution.Vy;
+	else 
+		md.initialization.vx=md.results.DiagnosticSolution.Vx;
+		md.initialization.vy=md.results.DiagnosticSolution.Vy;
+		md.initialization.vz=md.results.DiagnosticSolution.Vz;
+	end
+	md.initialization.vel=md.results.DiagnosticSolution.Vel;
+
+	if isfield(md.results.DiagnosticSolution,'Pressure'),
+		md.initialization.pressure=md.results.DiagnosticSolution.Pressure;
+	end
+	if md.rifts.numrifts,
+		if isfield(md.results.DiagnosticSolution,'riftproperties'),
+			md.rifts.riftproperties=md.results.DiagnosticSolution.riftproperties;
+		end
+	end
+	if md.inversion.iscontrol==1,
+		for control_parameters=md.inversion.control_parameters
+			%Will need to be updated... good luck ;)
+			md.(EnumToModelField(control_parameters))=md.results.DiagnosticSolution.(EnumToString(control_parameters));
+		end
+	end
+
+elseif strcmpi(string,'dakota'),
+	md.qmu.results=md.results.dakota;
+
+elseif strcmpi(string,'flaim'),
+	md.flaim.solution=md.results.FlaimSolution.solution;
+	md.flaim.quality =md.results.FlaimSolution.quality;
+
+elseif strcmpi(string,'transient'),
+	results=md.results.TransientSolution;
+	results2.Vel=NaN;
+	count=1;
+	for i=1:length(results),
+		if ~isempty(md.results.TransientSolution(i).Vel),
+			results2(count).Vel=md.results.TransientSolution(i).Vel;
+			results2(count).Surface=md.results.TransientSolution(i).Surface;
+			results2(count).Thickness=md.results.TransientSolution(i).Thickness;
+			results2(count).Bed=md.results.TransientSolution(i).Bed;
+			results2(count).Vx=md.results.TransientSolution(i).Vx;
+			results2(count).Vy=md.results.TransientSolution(i).Vy;
+			results2(count).time=md.results.TransientSolution(i).time;
+			results2(count).step=md.results.TransientSolution(i).step;
+			if ~strcmpi(md.groundingline.migration,'None'),
+				results2(count).ElementOnIceShelf=md.results.TransientSolution(i).ElementOnIceShelf;
+			end
+			count=count+1;
+		end
+	end
+	md.results.TransientSolution=results2;
+	clear results,results2;
+elseif strcmpi(string,'steadystate'),
+	md.initialization.vx=md.results.SteadystateSolution.Vx;
+	md.initialization.vy=md.results.SteadystateSolution.Vy;
+	if isfield(md.results.SteadystateSolution,'Vz'),
+		md.initialization.vz=md.results.SteadystateSolution.Vz;
+	end
+
+	md.initialization.vel=md.results.SteadystateSolution.Vel;
+	md.initialization.pressure=md.results.SteadystateSolution.Pressure;
+	md.initialization.temperature=md.results.SteadystateSolution.Temperature;
+	md.basalforcings.melting_rate=md.results.SteadystateSolution.BasalforcingsMeltingRate;
+
+	if md.inversion.iscontrol==1,
+		for control_parameters=md.inversion.control_parameters
+			md.(EnumToModelField(control_parameters))=md.results.SteadystateSolution.(EnumToString(control_parameters));
+		end
+	end
+
+elseif strcmpi(string,'thermal'),
+	md.initialization.temperature=md.results.ThermalSolution.Temperature;
+	md.basalforcings.melting_rate=md.results.ThermalSolution.BasalMeltingRate;
+elseif strcmpi(string,'hydrology'),
+	md.initialization.watercolumn=md.results.HydrologySolution.Watercolumn;
+
+else 
+	error(['tres error message: analysis ' string ' not supported yet!']);
+end
+end 
+function string=EnumToModelField(enum) % {{{
+	%ENUMTOMODELFIELD - output string of model field associated to enum
+	%
+	%   Usage:
+	%      string=EnumToModelField(enum)
+
+	disp('Warning: EnumToModelField is deprecated, it cannot work with new model definition. This function will be removed in the future');
+
+	switch enum,
+
+		case ThicknessEnum(), string='thickness'; return
+		case FrictionCoefficientEnum(), string='drag_coefficient'; return
+		case MaterialsRheologyBEnum(), string='rheology_B'; return
+		case MaterialsRheologyBbarEnum(), string='rheology_B'; return
+		case BalancethicknessThickeningRateEnum: string='dhdt'; return
+		case VxEnum(), string='vx'; return
+		case InversionVxObsEnum(), string='vx_obs'; return
+		case VyEnum(), string='vy'; return
+		case InversionVyObsEnum(), string='vy_obs'; return
+		case BasalforcingsMeltingRateEnum(), string='basal_melting_rate'; return
+		case SurfaceforcingsAccumulationRateEnum(), string='surface_accumulation_rate'; return
+		case SurfaceforcingsAblationRateEnum(), string='surface_ablation_rate'; return
+		case SurfaceforcingsMassBalanceEnum(), string='surface_mass_balance'; return
+		otherwise, error(['Enum ' num2str(enum)  ' not found associated to any model field']);
+
+		end
+	end % }}}
Index: issm/trunk-jpl/src/m/contrib/massbalance/outflow.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/massbalance/outflow.m	(revision 13008)
+++ issm/trunk-jpl/src/m/contrib/massbalance/outflow.m	(revision 13008)
@@ -0,0 +1,18 @@
+function flag=outflow(md);
+%OUTFLOW - flag nodes on outflux boundary
+%
+%   Usage:
+%      flag=outflow(md);
+
+A=md.mesh.segments(:,1);
+B=md.mesh.segments(:,2);
+Nx=-(md.mesh.y(A)-md.mesh.y(B));
+Ny=  md.mesh.x(A)-md.mesh.x(B);
+Vx=(md.initialization.vx(A)+md.initialization.vx(B))/2;
+Vy=(md.initialization.vy(A)+md.initialization.vy(B))/2;
+
+%dot product
+VdotN=Vx.*Nx+Vy.*Ny;
+
+flag=zeros(md.mesh.numberofvertices,1);
+flag(A(find(VdotN>0)))=1;
Index: issm/trunk-jpl/src/m/geometry/ThicknessCorrection.m
===================================================================
--- issm/trunk-jpl/src/m/geometry/ThicknessCorrection.m	(revision 13008)
+++ issm/trunk-jpl/src/m/geometry/ThicknessCorrection.m	(revision 13008)
@@ -0,0 +1,77 @@
+function md=ThicknessCorrection(md,varargin)
+%THICKNESSCORRECTION - correct the thickness of the ice shelf near the grounding line
+%
+%   This routine corrects the thickness and the bed on the transition zone
+%   by forcing the hydrostatic equilibrium.
+%   the thickness is modified as follows:
+%      thickness = (1-coeff) * thickness_observation + coeff * thickness_hydrostatic
+%   where:
+%      coeff=(d/distance)^2;
+%      distance=10km by default but can be specified
+%
+%   Usage:
+%      md=ThicknessCorrection(md,varargin);
+%
+%   Example:
+%      md=ThicknessCorrection(md);
+%      md=ThicknessCorrection(md,15000);
+
+%initialize thickness with the observations, and get hydrostatic thickness from the dem
+thickness=md.geometry.thickness;
+thickness_hydro=md.geometry.surface/(1-md.materials.rho_ice/md.materials.rho_water);
+hydrostatic_ratio=zeros(size(md.geometry.thickness));
+
+%get nodes on ice sheet and on ice shelf
+pos_shelf=find(~md.mask.vertexongroundedice);
+pos_GL=intersect(unique(md.mesh.elements(find(md.mask.elementongroundedice),:)),unique(md.mesh.elements(find(md.mask.elementonfloatingice),:)));
+debug=(length(pos_shelf)>50000);
+
+%check that there is a GL
+if isempty(pos_GL)
+	error('ThicknessCorrection error message: no grounding line has been detected. Check the model mask');
+end
+
+%get distance
+if nargin==2,
+	distance=varargin{1};
+else
+	distance=10000;
+end
+
+%modify thickness
+if (debug), fprintf('%s','      correction progress:   0.00 %'); end
+for i=1:length(pos_shelf)
+
+	if (debug & mod(i,100)==0),
+		fprintf('\b\b\b\b\b\b\b%5.2f%s',i/length(pos_shelf)*100,' %');
+	end
+
+	%search the node on ice sheet the closest to i
+	[d posd]=min(sqrt((md.mesh.x(pos_shelf(i))-md.mesh.x(pos_GL)).^2+(md.mesh.y(pos_shelf(i))-md.mesh.y(pos_GL)).^2));
+
+	if d>distance,
+
+		%if d > 15km, hydrostatic equilibrium
+		hydrostatic_ratio(pos_shelf(i))=1;
+		thickness(pos_shelf(i))=thickness_hydro(pos_shelf(i));
+
+	else
+
+		%else: quadratic combination of hydrostatic equilibrium and observations
+		hydrostatic_ratio(pos_shelf(i))=(d/distance)^2;
+		thickness(pos_shelf(i))=(1-hydrostatic_ratio(pos_shelf(i)))*thickness(pos_shelf(i))+hydrostatic_ratio(pos_shelf(i))*thickness_hydro(pos_shelf(i));
+
+	end
+end
+if (debug), fprintf('\b\b\b\b\b\b\b%5.2f%s\n',100,' %'); end
+
+%check the computed thickness
+minth=1/(1-md.materials.rho_ice/md.materials.rho_water);
+pos=find(isnan(thickness) | (thickness<=0));
+thickness(pos)=minth;
+hydrostatic_ratio(pos)=-1;
+
+%change bed to take into account the changes in thickness
+md.geometry.thickness=thickness;
+md.geometry.hydrostatic_ratio=hydrostatic_ratio;
+md.geometry.bed=md.geometry.surface-md.geometry.thickness;
Index: issm/trunk-jpl/src/m/geometry/recover_areas.m
===================================================================
--- issm/trunk-jpl/src/m/geometry/recover_areas.m	(revision 13008)
+++ issm/trunk-jpl/src/m/geometry/recover_areas.m	(revision 13008)
@@ -0,0 +1,22 @@
+function [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
+%RECOVER_AREAS - flag the element depending on the physical model that is assigned to them
+%
+%   This routine is called by setelementstype, do not use
+%
+%   Usage:
+%      [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
+
+	%go through varargin, extract options and plug them into subtype options, by order of appearance
+	options=pairoptions(varargin{:});
+	options=deleteduplicates(options,1);
+
+	%recover elements distribution
+	hutterflag  =FlagElements(md,getfieldvalue(options,'hutter',''));
+	macayealflag=FlagElements(md,getfieldvalue(options,'macayeal',''));
+	pattynflag  =FlagElements(md,getfieldvalue(options,'pattyn',''));
+	stokesflag  =FlagElements(md,getfieldvalue(options,'stokes',''));
+	filltype    =getfieldvalue(options,'fill','none');
+
+end %end function
+
+
Index: issm/trunk-jpl/src/m/geometry/recover_areas.py
===================================================================
--- issm/trunk-jpl/src/m/geometry/recover_areas.py	(revision 13008)
+++ issm/trunk-jpl/src/m/geometry/recover_areas.py	(revision 13008)
@@ -0,0 +1,26 @@
+from pairoptions import *
+from FlagElements import *
+
+def recover_areas(md,*args):
+	"""
+	RECOVER_AREAS - flag the element depending on the physical model that is assigned to them
+
+	   This routine is called by setelementstype, do not use
+
+	   Usage:
+	      [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
+	"""
+
+	#go through varargin, extract options and plug them into subtype options, by order of appearance
+	options=pairoptions(*args)
+#	options=deleteduplicates(options,1);
+
+	#recover elements distribution
+	hutterflag  =FlagElements(md,options.getfieldvalue('hutter',''))
+	macayealflag=FlagElements(md,options.getfieldvalue('macayeal',''))
+	pattynflag  =FlagElements(md,options.getfieldvalue('pattyn',''))
+	stokesflag  =FlagElements(md,options.getfieldvalue('stokes',''))
+	filltype    =options.getfieldvalue('fill','none')
+
+	return hutterflag,macayealflag,pattynflag,stokesflag,filltype
+
Index: issm/trunk-jpl/src/m/geometry/slope.m
===================================================================
--- issm/trunk-jpl/src/m/geometry/slope.m	(revision 13008)
+++ issm/trunk-jpl/src/m/geometry/slope.m	(revision 13008)
@@ -0,0 +1,32 @@
+function [sx,sy,s]=slope(md)
+%SLOPE - compute the surface slope
+%
+%   Usage:
+%      [sx,sy,s]=slope(md)
+
+%load some variables (it is much faster if the variab;es are loaded from md once for all) 
+if (md.mesh.dimension==2),
+	numberofelements=md.mesh.numberofelements;
+	numberofnodes=md.mesh.numberofvertices;
+	index=md.mesh.elements;
+	x=md.mesh.x; y=md.mesh.y;
+else
+	numberofelements=md.mesh.numberofelements2d;
+	numberofnodes=md.mesh.numberofvertices2d;
+	index=md.mesh.elements2d;
+	x=md.mesh.x2d; y=md.mesh.y2d;
+end
+
+%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
+[alpha beta]=GetNodalFunctionsCoeff(index,x,y);
+
+summation=[1;1;1];
+sx=(md.geometry.surface(index).*alpha)*summation;
+sy=(md.geometry.surface(index).*beta)*summation;
+s=sqrt(sx.^2+sy.^2);
+
+if md.mesh.dimension==3,
+	sx=project3d(md,'vector',sx,'type','element');
+	sy=project3d(md,'vector',sy,'type','element');
+	s=sqrt(sx.^2+sy.^2);
+end
Index: issm/trunk-jpl/src/m/interp/ProfileValues.m
===================================================================
--- issm/trunk-jpl/src/m/interp/ProfileValues.m	(revision 13008)
+++ issm/trunk-jpl/src/m/interp/ProfileValues.m	(revision 13008)
@@ -0,0 +1,21 @@
+function [Z,data_interp]=ProfileValues(md,data,xprof,yprof,resolution)
+%PROFILEVALUES - compute the value of a field on a vertical profile
+%
+%   This routine gets the value of a given field of the model on points
+%   given by filname (Argus type file)
+%
+%   Usage:
+%      [z,data]=ProfileValues(md,data,filename,resolution)
+%      [z,data]=ProfileValues(md,data,profile_structure,resolution)
+
+%Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
+offset=10^-3;
+bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.bed,1),xprof,yprof)+offset;
+surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),xprof,yprof)-offset;
+
+%Some useful parameters
+layers=ceil(mean(md.geometry.thickness)/resolution);
+Z=(bed:resolution:surface)';
+X=xprof*ones(size(Z));
+Y=yprof*ones(size(Z));
+data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X,Y,Z,NaN);
Index: issm/trunk-jpl/src/m/interp/SectionValues.m
===================================================================
--- issm/trunk-jpl/src/m/interp/SectionValues.m	(revision 13008)
+++ issm/trunk-jpl/src/m/interp/SectionValues.m	(revision 13008)
@@ -0,0 +1,127 @@
+function [index,X,Y,Z,S,data_interp]=SectionValues(md,data,infile,resolution)
+%SECTIONVALUES - compute the value of a field on a section
+%
+%   This routine gets the value of a given field of the model on points
+%   given by filname (Argus type file)
+%
+%   Usage:
+%      [elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution)
+%      [elements,x,y,z,s,data]=SectionValues(md,data,profile_structure,resolution)
+
+%check what we have for profile as input
+if ischar(infile),
+	%read infile:
+	profile=expread(infile);
+	nods=profile.nods;
+	x=profile.x;
+	y=profile.y;
+else
+	%read infile:
+	nods=infile.nods;
+	x=infile.x;
+	y=infile.y;
+end
+
+
+%get the specified resolution
+if isnumeric(resolution(1))
+	res_h=resolution(1);
+else
+	error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]')
+end
+if md.mesh.dimension==3
+	if (length(resolution)==2 & isnumeric(resolution(2)))
+		res_v=resolution(2);
+	else
+		error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]')
+	end
+end
+
+%initialization
+X=[]; %X-coordinate
+Y=[]; %Y-coordinate
+S=0;  %curvilinear coordinate
+
+for i=1:nods-1
+
+	x_start=x(i);
+	x_end=x(i+1);
+	y_start=y(i);
+	y_end=y(i+1);
+	s_start=S(end);
+
+	length_segment=sqrt((x_end-x_start)^2+(y_end-y_start)^2);
+	portion=ceil(length_segment/res_h);
+
+	x_segment=zeros(portion,1);
+	y_segment=zeros(portion,1);
+	s_segment=zeros(portion,1);
+
+	for j=1:portion
+		x_segment(j)=x_start+(j-1)*(x_end-x_start)/portion;
+		y_segment(j)=y_start+(j-1)*(y_end-y_start)/portion;
+		s_segment(j)=s_start+j*length_segment/portion;
+	end
+
+	%plug into X and Y
+	X=[X;x_segment];
+	Y=[Y;y_segment];
+	S=[S;s_segment];
+end
+X(end+1)=x(nods);
+Y(end+1)=y(nods);
+
+%Number of nodes:
+numberofnodes=size(X,1);
+
+%Compute Z
+Z=zeros(numberofnodes,1);
+
+%New mesh and Data interpolation
+if (md.mesh.dimension==2)
+
+	%Interpolation of data on specified points
+	data_interp=InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y);
+	%data_interp=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y);
+	%data_interp=griddata(md.mesh.x,md.mesh.y,data,X,Y);
+
+	%Compute index
+	index=[1:1:(numberofnodes-1);2:1:numberofnodes]';
+
+else
+
+	%vertically extrude mesh
+
+	%Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
+	offset=10^-3;
+	bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.bed,1),X,Y)+offset;
+	surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),X,Y)-offset;
+
+	%Some useful parameters
+	layers=ceil(mean(md.geometry.thickness)/res_v);
+	nodesperlayer=numberofnodes;
+	nodestot=nodesperlayer*layers;
+	elementsperlayer=nodesperlayer-1;
+	elementstot=(nodesperlayer-1)*(layers-1);
+
+	%initialization
+	X3=zeros(nodesperlayer*layers,1); Y3=zeros(nodesperlayer*layers,1); Z3=zeros(nodesperlayer*layers,1); S3=zeros(nodesperlayer*layers,1); index3=zeros(elementstot,4);
+
+	%Get new coordinates in 3d
+	for i=1:layers
+		X3(i:layers:end)=X;
+		Y3(i:layers:end)=Y;
+		Z3(i:layers:end)=bed+(i-1)*(surface-bed)/(layers-1);
+		S3(i:layers:end)=S;
+
+		if i<layers %Build index3 with quads
+			index3((i-1)*elementsperlayer+1:i*elementsperlayer,:)=[i:layers:nodestot-layers; i+1:layers:nodestot-layers; i+layers+1:layers:nodestot; i+layers:layers:nodestot]';
+		end
+	end
+
+	%Interpolation of data on specified points
+	data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,NaN);
+
+	%build outputs
+	X=X3; Y=Y3; Z=Z3;  S=S3; index=index3;
+end
Index: issm/trunk-jpl/src/m/meca/shear2d.m
===================================================================
--- issm/trunk-jpl/src/m/meca/shear2d.m	(revision 13008)
+++ issm/trunk-jpl/src/m/meca/shear2d.m	(revision 13008)
@@ -0,0 +1,23 @@
+function [sx,sy,sxy,s]=shear2d(md)
+%SHEAR2D - computes 2d strain rate
+%
+%   This routine computes the strain rate of 2d models
+%
+%   Usage:
+%      [sx,sy,sxy,s]=shear2d(md);
+%      s=shear2d(md);
+
+[alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); 
+
+summation=[1;1;1];
+sx=(md.initialization.vx(md.mesh.elements).*alpha)*summation;
+uy=(md.initialization.vx(md.mesh.elements).*beta)*summation;
+vx=(md.initialization.vy(md.mesh.elements).*alpha)*summation;
+sy=(md.initialization.vy(md.mesh.elements).*beta)*summation;						
+sxy=(uy+vx)/2;
+s=sqrt(sx.^2+sy.^2+sxy.^2+sx.*sy);
+
+%if user requested only one output, it must be the norm
+if nargout==1,
+	sx=s;
+end
Index: issm/trunk-jpl/src/m/model/MeltingGroundingLines.m
===================================================================
--- issm/trunk-jpl/src/m/model/MeltingGroundingLines.m	(revision 13007)
+++ 	(revision )
@@ -1,26 +1,0 @@
-function md=MeltingGroundingLines(md,distance,value)
-%MELTINGGROUNDINGLINES - set melting near grounding lines to a constant value
-%
-%   Usage:
-%      md=MeltingGroundingLines(md,distance,value)
-%
-
-%get nodes on ice sheet and on ice shelf
-pos_shelf=find(~md.mask.vertexongroundedice);
-pos_GL=intersect(unique(md.mesh.elements(find(md.mask.elementongroundedice),:)),unique(md.mesh.elements(find(md.mask.elementonfloatingice),:)));
-
-for i=1:length(pos_shelf)
-
-	if (mod(i,100)==0),
-		fprintf('\b\b\b\b\b\b\b%5.2f%s',i/length(pos_shelf)*100,' %');
-	end
-
-	%search the node on ice sheet the closest to i
-	[d posd]=min(sqrt((md.mesh.x(pos_shelf(i))-md.mesh.x(pos_GL)).^2+(md.mesh.y(pos_shelf(i))-md.mesh.y(pos_GL)).^2));
-
-	if d<distance,
-
-		md.melting(pos_shelf(i))=value;
-
-	end
-end
Index: issm/trunk-jpl/src/m/model/ProfileValues.m
===================================================================
--- issm/trunk-jpl/src/m/model/ProfileValues.m	(revision 13007)
+++ 	(revision )
@@ -1,21 +1,0 @@
-function [Z,data_interp]=ProfileValues(md,data,xprof,yprof,resolution)
-%PROFILEVALUES - compute the value of a field on a vertical profile
-%
-%   This routine gets the value of a given field of the model on points
-%   given by filname (Argus type file)
-%
-%   Usage:
-%      [z,data]=ProfileValues(md,data,filename,resolution)
-%      [z,data]=ProfileValues(md,data,profile_structure,resolution)
-
-%Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
-offset=10^-3;
-bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.bed,1),xprof,yprof)+offset;
-surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),xprof,yprof)-offset;
-
-%Some useful parameters
-layers=ceil(mean(md.geometry.thickness)/resolution);
-Z=(bed:resolution:surface)';
-X=xprof*ones(size(Z));
-Y=yprof*ones(size(Z));
-data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X,Y,Z,NaN);
Index: issm/trunk-jpl/src/m/model/PropagateFlagsUntilDistance.m
===================================================================
--- issm/trunk-jpl/src/m/model/PropagateFlagsUntilDistance.m	(revision 13007)
+++ 	(revision )
@@ -1,64 +1,0 @@
-function new_flags=PropagateFlagsUntilDistance(md,flags,distance)
-%PROPAGATEFLAGSUNTILDISTANCE
-%
-% Usage: 
-%              flags=PropagateFlagsUntilDistance(md,flags,distance)
-%
-%
-
-
-	
-new_flags=flags;
-
-%make 3d work in 2d: 
-if md.mesh.dimension==3,
-	md.mesh.x=md.mesh.x2d;
-	md.mesh.y=md.mesh.y2d;
-	md.mesh.elements=md.mesh.elements2d;
-end
-
-%find elements that are at the border of flags: 
-flag_elements=find(flags);
-conn=md.mesh.elementconnectivity(flag_elements,:);
-pos=find(conn);conn(pos)=~flags(conn(pos));
-sum_conn=sum(conn,2);
-border_elements=flag_elements(find(sum_conn>=1));
-
-
-%average x and y over elements: 
-x_elem=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
-y_elem=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
-
-while 1,
-
-	%keep copy of new_flags for this loop: 
-	new_flags_bak=new_flags;
-
-	%extend new flags by connectivity
-	pos=find(new_flags);
-
-	connected_elements=md.mesh.elementconnectivity(pos,:);
-	connected_elements=connected_elements(find(connected_elements));
-	new_flags(connected_elements)=1;
-
-	%get new elements: 
-	new_elements=find(new_flags & ~new_flags_bak);
-	if ~length(new_elements),
-		%we are done!
-		break;
-	end
-
-
-	%check which of these new elements are more than distance away from the border elements
-	for i=1:length(new_elements),
-		dist=sqrt(     (x_elem(border_elements)-x_elem(new_elements(i))).^2 + (y_elem(border_elements)-y_elem(new_elements(i))).^2)-distance;
-		if ~any(dist<0)
-			%none of the border elements are within distance, this element is outside out area of interest.
-			%ensure this element never gets found again in the connectivity.
-			pos=find(md.mesh.elementconnectivity==new_elements(i));
-			md.mesh.elementconnectivity(pos)=0;
-			%exclude this new element from the new_flags!
-			new_flags(new_elements(i))=0;
-		end
-	end
-end
Index: issm/trunk-jpl/src/m/model/README
===================================================================
--- issm/trunk-jpl/src/m/model/README	(revision 13007)
+++ 	(revision )
@@ -1,11 +1,0 @@
-This directory is similar to @model, in that it only deals
-with methods proper to the @model class. But we need a different 
-directory from @model, to make those methods public, not private. 
-The advantage of public methods is that they will use subsref and 
-susasgn to access data in a @model object.
-
-For ex: 
-calling md.x in a "@model" routine just accesses the x field in the md 
-structure. But calling md.x in a "model" routine will access the x field 
-through the susref routine in @model. This ensures that we protect the 
-data in @model classes from improper use.
Index: issm/trunk-jpl/src/m/model/SectionValues.m
===================================================================
--- issm/trunk-jpl/src/m/model/SectionValues.m	(revision 13007)
+++ 	(revision )
@@ -1,127 +1,0 @@
-function [index,X,Y,Z,S,data_interp]=SectionValues(md,data,infile,resolution)
-%SECTIONVALUES - compute the value of a field on a section
-%
-%   This routine gets the value of a given field of the model on points
-%   given by filname (Argus type file)
-%
-%   Usage:
-%      [elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution)
-%      [elements,x,y,z,s,data]=SectionValues(md,data,profile_structure,resolution)
-
-%check what we have for profile as input
-if ischar(infile),
-	%read infile:
-	profile=expread(infile);
-	nods=profile.nods;
-	x=profile.x;
-	y=profile.y;
-else
-	%read infile:
-	nods=infile.nods;
-	x=infile.x;
-	y=infile.y;
-end
-
-
-%get the specified resolution
-if isnumeric(resolution(1))
-	res_h=resolution(1);
-else
-	error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]')
-end
-if md.mesh.dimension==3
-	if (length(resolution)==2 & isnumeric(resolution(2)))
-		res_v=resolution(2);
-	else
-		error('SectionValues error message: wrong resolution type. Resolution must be an array [horizontal_resolution vertical_resolution]')
-	end
-end
-
-%initialization
-X=[]; %X-coordinate
-Y=[]; %Y-coordinate
-S=0;  %curvilinear coordinate
-
-for i=1:nods-1
-
-	x_start=x(i);
-	x_end=x(i+1);
-	y_start=y(i);
-	y_end=y(i+1);
-	s_start=S(end);
-
-	length_segment=sqrt((x_end-x_start)^2+(y_end-y_start)^2);
-	portion=ceil(length_segment/res_h);
-
-	x_segment=zeros(portion,1);
-	y_segment=zeros(portion,1);
-	s_segment=zeros(portion,1);
-
-	for j=1:portion
-		x_segment(j)=x_start+(j-1)*(x_end-x_start)/portion;
-		y_segment(j)=y_start+(j-1)*(y_end-y_start)/portion;
-		s_segment(j)=s_start+j*length_segment/portion;
-	end
-
-	%plug into X and Y
-	X=[X;x_segment];
-	Y=[Y;y_segment];
-	S=[S;s_segment];
-end
-X(end+1)=x(nods);
-Y(end+1)=y(nods);
-
-%Number of nodes:
-numberofnodes=size(X,1);
-
-%Compute Z
-Z=zeros(numberofnodes,1);
-
-%New mesh and Data interpolation
-if (md.mesh.dimension==2)
-
-	%Interpolation of data on specified points
-	data_interp=InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y);
-	%data_interp=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y);
-	%data_interp=griddata(md.mesh.x,md.mesh.y,data,X,Y);
-
-	%Compute index
-	index=[1:1:(numberofnodes-1);2:1:numberofnodes]';
-
-else
-
-	%vertically extrude mesh
-
-	%Get bed and surface for each 2d point, offset to make sure that it is inside the glacier system
-	offset=10^-3;
-	bed=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.bed,1),X,Y)+offset;
-	surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),X,Y)-offset;
-
-	%Some useful parameters
-	layers=ceil(mean(md.geometry.thickness)/res_v);
-	nodesperlayer=numberofnodes;
-	nodestot=nodesperlayer*layers;
-	elementsperlayer=nodesperlayer-1;
-	elementstot=(nodesperlayer-1)*(layers-1);
-
-	%initialization
-	X3=zeros(nodesperlayer*layers,1); Y3=zeros(nodesperlayer*layers,1); Z3=zeros(nodesperlayer*layers,1); S3=zeros(nodesperlayer*layers,1); index3=zeros(elementstot,4);
-
-	%Get new coordinates in 3d
-	for i=1:layers
-		X3(i:layers:end)=X;
-		Y3(i:layers:end)=Y;
-		Z3(i:layers:end)=bed+(i-1)*(surface-bed)/(layers-1);
-		S3(i:layers:end)=S;
-
-		if i<layers %Build index3 with quads
-			index3((i-1)*elementsperlayer+1:i*elementsperlayer,:)=[i:layers:nodestot-layers; i+1:layers:nodestot-layers; i+layers+1:layers:nodestot; i+layers:layers:nodestot]';
-		end
-	end
-
-	%Interpolation of data on specified points
-	data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,NaN);
-
-	%build outputs
-	X=X3; Y=Y3; Z=Z3;  S=S3; index=index3;
-end
Index: issm/trunk-jpl/src/m/model/ThicknessCorrection.m
===================================================================
--- issm/trunk-jpl/src/m/model/ThicknessCorrection.m	(revision 13007)
+++ 	(revision )
@@ -1,77 +1,0 @@
-function md=ThicknessCorrection(md,varargin)
-%THICKNESSCORRECTION - correct the thickness of the ice shelf near the grounding line
-%
-%   This routine corrects the thickness and the bed on the transition zone
-%   by forcing the hydrostatic equilibrium.
-%   the thickness is modified as follows:
-%      thickness = (1-coeff) * thickness_observation + coeff * thickness_hydrostatic
-%   where:
-%      coeff=(d/distance)^2;
-%      distance=10km by default but can be specified
-%
-%   Usage:
-%      md=ThicknessCorrection(md,varargin);
-%
-%   Example:
-%      md=ThicknessCorrection(md);
-%      md=ThicknessCorrection(md,15000);
-
-%initialize thickness with the observations, and get hydrostatic thickness from the dem
-thickness=md.geometry.thickness;
-thickness_hydro=md.geometry.surface/(1-md.materials.rho_ice/md.materials.rho_water);
-hydrostatic_ratio=zeros(size(md.geometry.thickness));
-
-%get nodes on ice sheet and on ice shelf
-pos_shelf=find(~md.mask.vertexongroundedice);
-pos_GL=intersect(unique(md.mesh.elements(find(md.mask.elementongroundedice),:)),unique(md.mesh.elements(find(md.mask.elementonfloatingice),:)));
-debug=(length(pos_shelf)>50000);
-
-%check that there is a GL
-if isempty(pos_GL)
-	error('ThicknessCorrection error message: no grounding line has been detected. Check the model mask');
-end
-
-%get distance
-if nargin==2,
-	distance=varargin{1};
-else
-	distance=10000;
-end
-
-%modify thickness
-if (debug), fprintf('%s','      correction progress:   0.00 %'); end
-for i=1:length(pos_shelf)
-
-	if (debug & mod(i,100)==0),
-		fprintf('\b\b\b\b\b\b\b%5.2f%s',i/length(pos_shelf)*100,' %');
-	end
-
-	%search the node on ice sheet the closest to i
-	[d posd]=min(sqrt((md.mesh.x(pos_shelf(i))-md.mesh.x(pos_GL)).^2+(md.mesh.y(pos_shelf(i))-md.mesh.y(pos_GL)).^2));
-
-	if d>distance,
-
-		%if d > 15km, hydrostatic equilibrium
-		hydrostatic_ratio(pos_shelf(i))=1;
-		thickness(pos_shelf(i))=thickness_hydro(pos_shelf(i));
-
-	else
-
-		%else: quadratic combination of hydrostatic equilibrium and observations
-		hydrostatic_ratio(pos_shelf(i))=(d/distance)^2;
-		thickness(pos_shelf(i))=(1-hydrostatic_ratio(pos_shelf(i)))*thickness(pos_shelf(i))+hydrostatic_ratio(pos_shelf(i))*thickness_hydro(pos_shelf(i));
-
-	end
-end
-if (debug), fprintf('\b\b\b\b\b\b\b%5.2f%s\n',100,' %'); end
-
-%check the computed thickness
-minth=1/(1-md.materials.rho_ice/md.materials.rho_water);
-pos=find(isnan(thickness) | (thickness<=0));
-thickness(pos)=minth;
-hydrostatic_ratio(pos)=-1;
-
-%change bed to take into account the changes in thickness
-md.geometry.thickness=thickness;
-md.geometry.hydrostatic_ratio=hydrostatic_ratio;
-md.geometry.bed=md.geometry.surface-md.geometry.thickness;
Index: issm/trunk-jpl/src/m/model/inversions/misfit.m
===================================================================
--- issm/trunk-jpl/src/m/model/inversions/misfit.m	(revision 13008)
+++ issm/trunk-jpl/src/m/model/inversions/misfit.m	(revision 13008)
@@ -0,0 +1,37 @@
+function J=misfit(md)
+%MISFIT - compute misfit
+%
+%   Usage:
+%      J=misfit(md)
+%
+%   Example:
+%      J=misfit(md)
+%
+
+if md.mesh.dimension==2,
+	elements=md.mesh.elements;
+	x=md.mesh.x;
+	y=md.mesh.y;
+	vx=md.initialization.vx;
+	vy=md.initialization.vy;
+	vx_obs=md.inversion.vx_obs;
+	vy_obs=md.inversion.vy_obs;
+else
+	elements=md.mesh.elements2d;
+	x=md.mesh.x2d;
+	y=md.mesh.y2d;
+	vx=project2d(md,md.initialization.vx,md.mesh.numberoflayers);
+	vy=project2d(md,md.initialization.vy,md.mesh.numberoflayers);
+	vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
+	vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
+end
+
+%compute areas;
+areas=GetAreas(elements,x,y);
+
+%compute delta v on elements
+deltav=1/2*(   (vx-vx_obs).^2+(vy-vy_obs).^2)/md.constants.yts^2;
+deltav_elem=deltav(elements)*[1;1;1]/3;
+
+%compute misfit
+J=sum(deltav_elem.*areas);
Index: issm/trunk-jpl/src/m/model/mesh/bamg.m
===================================================================
--- issm/trunk-jpl/src/m/model/mesh/bamg.m	(revision 13007)
+++ issm/trunk-jpl/src/m/model/mesh/bamg.m	(revision 13008)
@@ -339,2 +339,145 @@
 	error('Output mesh has orphans. Decrease MaxCornerAngle to prevent outside points (ex: 0.01)');
 end
+end 
+
+function geom=processgeometry(geom,tol,outline); % {{{
+
+%Deal with edges
+disp('Checking Edge crossing...');
+i=0;
+while (i<size(geom.Edges,1)),
+
+	%edge counter
+	i=i+1;
+
+	%Get coordinates
+	x1=geom.Vertices(geom.Edges(i,1),1);
+	y1=geom.Vertices(geom.Edges(i,1),2);
+	x2=geom.Vertices(geom.Edges(i,2),1);
+	y2=geom.Vertices(geom.Edges(i,2),2);
+	color1=geom.Edges(i,3);
+
+	j=i; %test edges located AFTER i only
+	while (j<size(geom.Edges,1)),
+
+		%edge counter
+		j=j+1;
+
+		%Skip if the two edges already have a vertex in common
+		if any(ismember(geom.Edges(i,1:2),geom.Edges(j,1:2))),
+			continue
+		end
+
+		%Get coordinates
+		x3=geom.Vertices(geom.Edges(j,1),1);
+		y3=geom.Vertices(geom.Edges(j,1),2);
+		x4=geom.Vertices(geom.Edges(j,2),1);
+		y4=geom.Vertices(geom.Edges(j,2),2);
+		color2=geom.Edges(j,3);
+
+		%Check if the two edges are crossing one another
+		if SegIntersect([x1 y1; x2 y2],[x3 y3; x4 y4]),
+
+			%Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
+			x=det([det([x1 y1; x2 y2])  x1-x2;det([x3 y3; x4 y4])  x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+			y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+
+			%Add vertex to the list of vertices
+			geom.Vertices(end+1,:)=[x y min(color1,color2)];
+			id=size(geom.Vertices,1);
+
+			%Update edges i and j
+			edgei=geom.Edges(i,:);
+			edgej=geom.Edges(j,:);
+			geom.Edges(i,:)    =[edgei(1) id       edgei(3)];
+			geom.Edges(end+1,:)=[id       edgei(2) edgei(3)];
+			geom.Edges(j,:)    =[edgej(1) id       edgej(3)];
+			geom.Edges(end+1,:)=[id       edgej(2) edgej(3)];
+
+			%update current edge second tip
+			x2=x; y2=y;
+		end
+	end
+
+end
+
+%Check point outside
+disp('Checking for points outside the domain...');
+i=0;
+num=0;
+while (i<size(geom.Vertices,1)),
+
+	%vertex counter
+	i=i+1;
+
+	%Get coordinates
+	x=geom.Vertices(i,1);
+	y=geom.Vertices(i,2);
+	color=geom.Vertices(i,3);
+
+	%Check that the point is inside the domain
+	if (color~=1 & ~ContourToNodes(x,y,outline(1),1)),
+
+		%Remove points from list of Vertices
+		num=num+1;
+		geom.Vertices(i,:)=[];
+
+		%update edges
+		[posedges dummy]=find(geom.Edges==i);
+		geom.Edges(posedges,:)=[];
+		posedges=find(geom.Edges>i);
+		geom.Edges(posedges)=geom.Edges(posedges)-1;
+
+		%update counter
+		i=i-1;
+	end
+end
+if num,
+	disp(['WARNING: ' num2str(num) ' points outside the domain outline have been removed']);
+end
+
+%Check point spacing
+if ~isnan(tol),
+	disp('Checking point spacing...');
+	i=0;
+	while (i<size(geom.Vertices,1)),
+
+		%vertex counter
+		i=i+1;
+
+		%Get coordinates
+		x1=geom.Vertices(i,1);
+		y1=geom.Vertices(i,2);
+
+		j=i; %test edges located AFTER i only
+		while (j<size(geom.Vertices,1)),
+
+			%vertex counter
+			j=j+1;
+
+			%Get coordinates
+			x2=geom.Vertices(j,1);
+			y2=geom.Vertices(j,2);
+
+			%Check whether the two vertices are too close
+			if ((x2-x1)^2+(y2-y1)^2<tol^2)
+
+				%Remove points from list of Vertices
+				geom.Vertices(j,:)=[];
+
+				%update edges
+				posedges=find(ismember(geom.Edges,j));
+				geom.Edges(posedges)=i;
+				posedges=find(geom.Edges>j);
+				geom.Edges(posedges)=geom.Edges(posedges)-1;
+
+				%update counter
+				j=j-1;
+
+			end
+		end
+	end
+end
+%remove empty edges
+geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[];
+end % }}}
Index: issm/trunk-jpl/src/m/model/misfit.m
===================================================================
--- issm/trunk-jpl/src/m/model/misfit.m	(revision 13007)
+++ 	(revision )
@@ -1,37 +1,0 @@
-function J=misfit(md)
-%MISFIT - compute misfit
-%
-%   Usage:
-%      J=misfit(md)
-%
-%   Example:
-%      J=misfit(md)
-%
-
-if md.mesh.dimension==2,
-	elements=md.mesh.elements;
-	x=md.mesh.x;
-	y=md.mesh.y;
-	vx=md.initialization.vx;
-	vy=md.initialization.vy;
-	vx_obs=md.inversion.vx_obs;
-	vy_obs=md.inversion.vy_obs;
-else
-	elements=md.mesh.elements2d;
-	x=md.mesh.x2d;
-	y=md.mesh.y2d;
-	vx=project2d(md,md.initialization.vx,md.mesh.numberoflayers);
-	vy=project2d(md,md.initialization.vy,md.mesh.numberoflayers);
-	vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
-	vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
-end
-
-%compute areas;
-areas=GetAreas(elements,x,y);
-
-%compute delta v on elements
-deltav=1/2*(   (vx-vx_obs).^2+(vy-vy_obs).^2)/md.constants.yts^2;
-deltav_elem=deltav(elements)*[1;1;1]/3;
-
-%compute misfit
-J=sum(deltav_elem.*areas);
Index: issm/trunk-jpl/src/m/model/modelsextract.m
===================================================================
--- issm/trunk-jpl/src/m/model/modelsextract.m	(revision 13007)
+++ 	(revision )
@@ -1,87 +1,0 @@
-function md_list=modelsextract(md,flags,minel,varargin)
-%modelsextract - extract several self contained models according to a list of element flags.
-%
-%   The difference between this routine and the modelextract.m routine (without an 's') is that 
-%   as many models are extracted as there are closed contours defined in area. 
-%   This routine is needed for example when doing data assimilation of ice shelves in Antarctica. 
-%   Many independent ice shelves are present, and we don't want data assimilation on one ice shelf 
-%   to be hindered by another totally independent ice shelf.
-%
-%   Usage:
-%      md_list=modelsextract(md,elementfalgs,minel);
-%
-%   Examples:
-%      md_list=modelsextract(md,md.mask.elementonfloatingice,1000);
-%
-%   See also: EXTRUDE, COLLAPSE, MODELEXTRACT
-
-disp('selecting pools of elements');
-%go through flags and build as many independent element flags as there are groups of connected 1s
-%in flags.
-
-%2D or 3D?
-if md.mesh.dimension==3,
-	numberofelements=md.mesh.numberofelements2d; %this will be forgotten when we get out.
-	flags=project2d(md,flags,1);
-else
-	numberofelements=md.mesh.numberofelements;
-end
-
-%recover extra arguments: 
-distance=0;
-if nargin==4,
-	distance=varargin{1};
-end
-
-flag_list=cell(0,1);
-
-for i=1:size(flags,1),
-
-	if (flags(i)),
-
-		%ok, we are sure element i is part of a new pool.
-		pool=zeros(numberofelements,1);
-		pool=PropagateFlagsFromConnectivity(md.mesh.elementconnectivity,pool,i,flags);
-		flag_list{end+1,1}=pool;
-		
-		%speed up rest of computation by taking pool out of flags: 
-		pos=find(pool);flags(pos)=0;
-
-	end
-end
-
-%go through flag_list and discard any pool of less than minel elements: 
-ex_pos=[];
-for i=1:length(flag_list),
-	if length(find(flag_list{i}))<minel,
-		ex_pos=[ex_pos; i];
-	end
-end
-flag_list(ex_pos)=[];
-
-%now, if distance was specified, expand the flag_list by distance km: 
-if distance,
-	for i=1:length(flag_list),
-		flag_list{i}=PropagateFlagsUntilDistance(md,flag_list{i},distance);
-	end
-end
-
-%now, go use the pools of flags to extract models: 
-disp(['extracting ' num2str(size(flag_list,1)) ' models']);
-models=cell(0,1);
-
-for i=1:size(flag_list,1),
-	disp(['   ' num2str(i) '/' num2str(size(flag_list,1))]);
-	if md.mesh.dimension==3,
-		flags2d=flag_list{i};
-		realflags=project3d(md,flags2d,'element');
-	else
-		realflags=flag_list{i};
-	end
-	models{end+1,1}=modelextract(md,realflags);
-end
-
-%return model list
-md_list=modellist(models);
-
-end %end of this function
Index: issm/trunk-jpl/src/m/model/modelsextractfromdomains.m
===================================================================
--- issm/trunk-jpl/src/m/model/modelsextractfromdomains.m	(revision 13007)
+++ 	(revision )
@@ -1,25 +1,0 @@
-function md_list=modelsextractfromdomains(md,directory)
-%modelsextractfromdomains- extract several self contained models according to a list of domains
-%
-%   Usage:
-%      md_list=modelsextractfromdomains(md,'Basins/');
-%
-%   Examples:
-%      md_list=modelsextract(md,'Basins/');
-%
-%   See also: MODELSEXTRACTS, MODELEXTRACT
-
-%go into directory and get list of files.
-cd(directory);
-basins=listfiles;
-cd ..
-
-models=cell(0,1);
-for i=1:length(basins),
-	models{end+1,1}=modelextract(md,[directory '/' basins{i}]);
-end
-
-%return model list: 
-md_list=modellist(models);
-
-end 
Index: issm/trunk-jpl/src/m/model/modis.m
===================================================================
--- issm/trunk-jpl/src/m/model/modis.m	(revision 13007)
+++ 	(revision )
@@ -1,35 +1,0 @@
-function [xm,ym,modis]=modis(modisgeotif,xlim,ylim)
-%MODIS - from modis geotiff, return image
-%
-%   Usage:
-%      [xm,ym,modis]=modis(modisgeotif,xlim,ylim)
-%
-
-%find gdal coordinates
-x0=min(xlim);
-x1=max(xlim);
-
-y0=min(ylim);
-y1=max(ylim);
-
-%Get path  to gdal binaries
-path_gdal=[issmdir() '/externalpackages/gdal/install/bin/'];
-
-%Was gdal compiled? 
-if ~exist([path_gdal 'gdal_translate']),
-	error(['modis error message: GDAL library needs to be compiled to use this routine. Compile GDAL in ' issmdir() '/externalpackages/gdal to use this routine.']);
-end
-
-inputname='./temp.tif';
-system([path_gdal 'gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' modisgeotif ' ' inputname ]);
-
-%Read in temp.tif:
-modis=double(flipud(imread('temp.tif','TIFF')));
-xm=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
-ym=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
-
-%Erase image
-system('rm -rf ./temp.tif');
-
-
-end
Index: issm/trunk-jpl/src/m/model/outflow.m
===================================================================
--- issm/trunk-jpl/src/m/model/outflow.m	(revision 13007)
+++ 	(revision )
@@ -1,18 +1,0 @@
-function flag=outflow(md);
-%OUTFLOW - flag nodes on outflux boundary
-%
-%   Usage:
-%      flag=outflow(md);
-
-A=md.mesh.segments(:,1);
-B=md.mesh.segments(:,2);
-Nx=-(md.mesh.y(A)-md.mesh.y(B));
-Ny=  md.mesh.x(A)-md.mesh.x(B);
-Vx=(md.initialization.vx(A)+md.initialization.vx(B))/2;
-Vy=(md.initialization.vy(A)+md.initialization.vy(B))/2;
-
-%dot product
-VdotN=Vx.*Nx+Vy.*Ny;
-
-flag=zeros(md.mesh.numberofvertices,1);
-flag(A(find(VdotN>0)))=1;
Index: issm/trunk-jpl/src/m/model/plot/radarpower.m
===================================================================
--- issm/trunk-jpl/src/m/model/plot/radarpower.m	(revision 13008)
+++ issm/trunk-jpl/src/m/model/plot/radarpower.m	(revision 13008)
@@ -0,0 +1,118 @@
+function md=radarpower(md,varargin)
+%RADARPOWER - overlay a power radar image on an existing mesh
+%
+%   This routine will overlay a power radar image on an existing mesh.
+%   The power amplitude will be output to vel for now.
+%   In the future, think about a field to hold this value.
+%
+%   Usage:
+%      md=radarpower(md,options);
+%      md=radarpower(md)
+
+%If gdal does not work, uncomment the following line
+%setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/');
+%Parse inputs
+if nargin==1,
+	options=pairoptions;
+else
+	options=varargin{:};
+	if ~isa(options,'pairoptions'),
+		options=pairoptions(varargin{:});
+	end
+end
+
+highres=getfieldvalue(options,'highres',0);
+xlim=getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
+ylim=getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]);
+posting=getfieldvalue(options,'posting',0); % 0 -> image posting default
+
+%find gdal coordinates
+x0=min(xlim); x1=max(xlim);
+y0=min(ylim); y1=max(ylim);
+
+%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one.
+if ~exist(options,'overlay_image'),
+	if strcmpi(md.mesh.hemisphere,'n'),
+		if ~exist([jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg']),
+			error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']);
+		end
+		name = 'mog150_greenland_map';
+		%name = 'mog100_hp1_v10';
+		%name = 'mog500_hp1_v10';
+		jpgim=[jplsvn() '/projects/ModelData/MOG/' name '.jpg'];
+		geom=load([jplsvn() '/projects/ModelData/MOG/' name '.jpgw'],'ascii');
+
+		%geom:   xposting nbcols nbrows yposting xmin ymax
+		xmin=max(geom(5),x0);
+		xmax=min(geom(5)+geom(1)*geom(2),x1);
+		ymin=max(geom(6)-geom(3)*geom(4),y0);
+		ymax=min(geom(6),y1);
+
+		firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min
+		firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max
+		numcols=floor((xmax-xmin)/geom(1)); % x posting
+		numrows=floor((ymax-ymin)/geom(4)); % y posting
+		pixelskip=max(1,ceil(posting/geom(1)));
+
+		%Read and crop file
+		disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
+		im=imread(jpgim);
+		im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1);
+		md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
+		md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
+		md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);
+
+	elseif strcmpi(md.mesh.hemisphere,'s'),
+		if highres,
+			if ~exist([jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']),
+				error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']);
+			end
+			geotiff_name=[jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
+		else
+			if ~exist([jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),
+				error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);
+			end
+			geotiff_name=[jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
+		end
+
+		%Name of image
+		inputname='./temp.tif';
+		eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
+
+		%Read in temp.tif:
+		im=imread('temp.tif','TIFF');
+		pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
+		md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
+		md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
+		md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
+
+		%Erase image
+		system('rm -rf ./temp.tif');
+
+	else
+		error('field hemisphere should either be ''n'' or ''s''');
+	end
+else
+	%ok, user provided an image. check we also have overlay_xlim and overlay_ylim  options, to know what range of coordinates the image covers.
+	if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')),
+		error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option');
+	end
+	overlay_image=getfieldvalue(options,'overlay_image');
+	overlay_xlim=getfieldvalue(options,'overlay_xlim');
+	overlay_ylim=getfieldvalue(options,'overlay_ylim');
+	overlay_xposting=getfieldvalue(options,'overlay_xposting');
+	overlay_yposting=getfieldvalue(options,'overlay_yposting');
+
+	sizex=floor((x1-x0)/overlay_xposting);
+	sizey=floor((y1-y0)/overlay_yposting);
+	topleftx=floor((x0-overlay_xlim(1))/overlay_xposting); % x min
+	toplefty=floor((overlay_ylim(2)-y1)/overlay_yposting); % y max
+
+	%Read and crop file
+	disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
+	im=imread(overlay_image);
+	im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
+	md.radaroverlay.pwr=double(flipud(im));
+	md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
+	md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
+end
Index: issm/trunk-jpl/src/m/model/processgeometry.m
===================================================================
--- issm/trunk-jpl/src/m/model/processgeometry.m	(revision 13007)
+++ 	(revision )
@@ -1,140 +1,0 @@
-function geom=processgeometry(geom,tol,outline);
-
-%Deal with edges
-disp('Checking Edge crossing...');
-i=0;
-while (i<size(geom.Edges,1)),
-
-	%edge counter
-	i=i+1;
-
-	%Get coordinates
-	x1=geom.Vertices(geom.Edges(i,1),1);
-	y1=geom.Vertices(geom.Edges(i,1),2);
-	x2=geom.Vertices(geom.Edges(i,2),1);
-	y2=geom.Vertices(geom.Edges(i,2),2);
-	color1=geom.Edges(i,3);
-
-	j=i; %test edges located AFTER i only
-	while (j<size(geom.Edges,1)),
-
-		%edge counter
-		j=j+1;
-
-		%Skip if the two edges already have a vertex in common
-		if any(ismember(geom.Edges(i,1:2),geom.Edges(j,1:2))),
-			continue
-		end
-
-		%Get coordinates
-		x3=geom.Vertices(geom.Edges(j,1),1);
-		y3=geom.Vertices(geom.Edges(j,1),2);
-		x4=geom.Vertices(geom.Edges(j,2),1);
-		y4=geom.Vertices(geom.Edges(j,2),2);
-		color2=geom.Edges(j,3);
-
-		%Check if the two edges are crossing one another
-		if SegIntersect([x1 y1; x2 y2],[x3 y3; x4 y4]),
-
-			%Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
-			x=det([det([x1 y1; x2 y2])  x1-x2;det([x3 y3; x4 y4])  x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
-			y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
-
-			%Add vertex to the list of vertices
-			geom.Vertices(end+1,:)=[x y min(color1,color2)];
-			id=size(geom.Vertices,1);
-
-			%Update edges i and j
-			edgei=geom.Edges(i,:);
-			edgej=geom.Edges(j,:);
-			geom.Edges(i,:)    =[edgei(1) id       edgei(3)];
-			geom.Edges(end+1,:)=[id       edgei(2) edgei(3)];
-			geom.Edges(j,:)    =[edgej(1) id       edgej(3)];
-			geom.Edges(end+1,:)=[id       edgej(2) edgej(3)];
-
-			%update current edge second tip
-			x2=x; y2=y;
-		end
-	end
-
-end
-
-%Check point outside
-disp('Checking for points outside the domain...');
-i=0;
-num=0;
-while (i<size(geom.Vertices,1)),
-
-	%vertex counter
-	i=i+1;
-
-	%Get coordinates
-	x=geom.Vertices(i,1);
-	y=geom.Vertices(i,2);
-	color=geom.Vertices(i,3);
-
-	%Check that the point is inside the domain
-	if (color~=1 & ~ContourToNodes(x,y,outline(1),1)),
-
-		%Remove points from list of Vertices
-		num=num+1;
-		geom.Vertices(i,:)=[];
-
-		%update edges
-		[posedges dummy]=find(geom.Edges==i);
-		geom.Edges(posedges,:)=[];
-		posedges=find(geom.Edges>i);
-		geom.Edges(posedges)=geom.Edges(posedges)-1;
-
-		%update counter
-		i=i-1;
-	end
-end
-if num,
-	disp(['WARNING: ' num2str(num) ' points outside the domain outline have been removed']);
-end
-
-%Check point spacing
-if ~isnan(tol),
-	disp('Checking point spacing...');
-	i=0;
-	while (i<size(geom.Vertices,1)),
-
-		%vertex counter
-		i=i+1;
-
-		%Get coordinates
-		x1=geom.Vertices(i,1);
-		y1=geom.Vertices(i,2);
-
-		j=i; %test edges located AFTER i only
-		while (j<size(geom.Vertices,1)),
-
-			%vertex counter
-			j=j+1;
-
-			%Get coordinates
-			x2=geom.Vertices(j,1);
-			y2=geom.Vertices(j,2);
-
-			%Check whether the two vertices are too close
-			if ((x2-x1)^2+(y2-y1)^2<tol^2)
-
-				%Remove points from list of Vertices
-				geom.Vertices(j,:)=[];
-
-				%update edges
-				posedges=find(ismember(geom.Edges,j));
-				geom.Edges(posedges)=i;
-				posedges=find(geom.Edges>j);
-				geom.Edges(posedges)=geom.Edges(posedges)-1;
-
-				%update counter
-				j=j-1;
-
-			end
-		end
-	end
-end
-%remove empty edges
-geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[];
Index: issm/trunk-jpl/src/m/model/qstat.m
===================================================================
--- issm/trunk-jpl/src/m/model/qstat.m	(revision 13007)
+++ 	(revision )
@@ -1,8 +1,0 @@
-function qstat(md)
-%QSTAT - check job status on remote cluster
-%
-%   Usage:
-%      qstat(md)
-
-%run qstat command on remote cluster
-issmssh(md.cluster,['qstat -a']);
Index: issm/trunk-jpl/src/m/model/radarpower.m
===================================================================
--- issm/trunk-jpl/src/m/model/radarpower.m	(revision 13007)
+++ 	(revision )
@@ -1,118 +1,0 @@
-function md=radarpower(md,varargin)
-%RADARPOWER - overlay a power radar image on an existing mesh
-%
-%   This routine will overlay a power radar image on an existing mesh.
-%   The power amplitude will be output to vel for now.
-%   In the future, think about a field to hold this value.
-%
-%   Usage:
-%      md=radarpower(md,options);
-%      md=radarpower(md)
-
-%If gdal does not work, uncomment the following line
-%setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/');
-%Parse inputs
-if nargin==1,
-	options=pairoptions;
-else
-	options=varargin{:};
-	if ~isa(options,'pairoptions'),
-		options=pairoptions(varargin{:});
-	end
-end
-
-highres=getfieldvalue(options,'highres',0);
-xlim=getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
-ylim=getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]);
-posting=getfieldvalue(options,'posting',0); % 0 -> image posting default
-
-%find gdal coordinates
-x0=min(xlim); x1=max(xlim);
-y0=min(ylim); y1=max(ylim);
-
-%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one.
-if ~exist(options,'overlay_image'),
-	if strcmpi(md.mesh.hemisphere,'n'),
-		if ~exist([jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg']),
-			error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']);
-		end
-		name = 'mog150_greenland_map';
-		%name = 'mog100_hp1_v10';
-		%name = 'mog500_hp1_v10';
-		jpgim=[jplsvn() '/projects/ModelData/MOG/' name '.jpg'];
-		geom=load([jplsvn() '/projects/ModelData/MOG/' name '.jpgw'],'ascii');
-
-		%geom:   xposting nbcols nbrows yposting xmin ymax
-		xmin=max(geom(5),x0);
-		xmax=min(geom(5)+geom(1)*geom(2),x1);
-		ymin=max(geom(6)-geom(3)*geom(4),y0);
-		ymax=min(geom(6),y1);
-
-		firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min
-		firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max
-		numcols=floor((xmax-xmin)/geom(1)); % x posting
-		numrows=floor((ymax-ymin)/geom(4)); % y posting
-		pixelskip=max(1,ceil(posting/geom(1)));
-
-		%Read and crop file
-		disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
-		im=imread(jpgim);
-		im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1);
-		md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
-		md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
-		md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);
-
-	elseif strcmpi(md.mesh.hemisphere,'s'),
-		if highres,
-			if ~exist([jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']),
-				error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']);
-			end
-			geotiff_name=[jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
-		else
-			if ~exist([jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),
-				error(['radarpower error message: file ' jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);
-			end
-			geotiff_name=[jplsvn() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
-		end
-
-		%Name of image
-		inputname='./temp.tif';
-		eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);
-
-		%Read in temp.tif:
-		im=imread('temp.tif','TIFF');
-		pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
-		md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
-		md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
-		md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
-
-		%Erase image
-		system('rm -rf ./temp.tif');
-
-	else
-		error('field hemisphere should either be ''n'' or ''s''');
-	end
-else
-	%ok, user provided an image. check we also have overlay_xlim and overlay_ylim  options, to know what range of coordinates the image covers.
-	if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')),
-		error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option');
-	end
-	overlay_image=getfieldvalue(options,'overlay_image');
-	overlay_xlim=getfieldvalue(options,'overlay_xlim');
-	overlay_ylim=getfieldvalue(options,'overlay_ylim');
-	overlay_xposting=getfieldvalue(options,'overlay_xposting');
-	overlay_yposting=getfieldvalue(options,'overlay_yposting');
-
-	sizex=floor((x1-x0)/overlay_xposting);
-	sizey=floor((y1-y0)/overlay_yposting);
-	topleftx=floor((x0-overlay_xlim(1))/overlay_xposting); % x min
-	toplefty=floor((overlay_ylim(2)-y1)/overlay_yposting); % y max
-
-	%Read and crop file
-	disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
-	im=imread(overlay_image);
-	im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
-	md.radaroverlay.pwr=double(flipud(im));
-	md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
-	md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
-end
Index: issm/trunk-jpl/src/m/model/recover_areas.m
===================================================================
--- issm/trunk-jpl/src/m/model/recover_areas.m	(revision 13007)
+++ 	(revision )
@@ -1,22 +1,0 @@
-function [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
-%RECOVER_AREAS - flag the element depending on the physical model that is assigned to them
-%
-%   This routine is called by setelementstype, do not use
-%
-%   Usage:
-%      [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
-
-	%go through varargin, extract options and plug them into subtype options, by order of appearance
-	options=pairoptions(varargin{:});
-	options=deleteduplicates(options,1);
-
-	%recover elements distribution
-	hutterflag  =FlagElements(md,getfieldvalue(options,'hutter',''));
-	macayealflag=FlagElements(md,getfieldvalue(options,'macayeal',''));
-	pattynflag  =FlagElements(md,getfieldvalue(options,'pattyn',''));
-	stokesflag  =FlagElements(md,getfieldvalue(options,'stokes',''));
-	filltype    =getfieldvalue(options,'fill','none');
-
-end %end function
-
-
Index: issm/trunk-jpl/src/m/model/recover_areas.py
===================================================================
--- issm/trunk-jpl/src/m/model/recover_areas.py	(revision 13007)
+++ 	(revision )
@@ -1,26 +1,0 @@
-from pairoptions import *
-from FlagElements import *
-
-def recover_areas(md,*args):
-	"""
-	RECOVER_AREAS - flag the element depending on the physical model that is assigned to them
-
-	   This routine is called by setelementstype, do not use
-
-	   Usage:
-	      [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
-	"""
-
-	#go through varargin, extract options and plug them into subtype options, by order of appearance
-	options=pairoptions(*args)
-#	options=deleteduplicates(options,1);
-
-	#recover elements distribution
-	hutterflag  =FlagElements(md,options.getfieldvalue('hutter',''))
-	macayealflag=FlagElements(md,options.getfieldvalue('macayeal',''))
-	pattynflag  =FlagElements(md,options.getfieldvalue('pattyn',''))
-	stokesflag  =FlagElements(md,options.getfieldvalue('stokes',''))
-	filltype    =options.getfieldvalue('fill','none')
-
-	return hutterflag,macayealflag,pattynflag,stokesflag,filltype
-
Index: issm/trunk-jpl/src/m/model/regional/regionaltransient2d.m
===================================================================
--- issm/trunk-jpl/src/m/model/regional/regionaltransient2d.m	(revision 13008)
+++ issm/trunk-jpl/src/m/model/regional/regionaltransient2d.m	(revision 13008)
@@ -0,0 +1,158 @@
+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(mde,'hmin',hmin,'hmax',hmax,'field',[mde.inversion.vel_obs mde.geometry.surface],'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);
+
+		%copy field, interpolated to new mesh
+		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
+
+	%Read transient velocities and thickness, looping through only the populated times
+	spcx=[];
+	spcy=[];
+	spct=[];
+	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);
+			thickness=PatchToVec(md1.results.TransientSolution(t).Thickness);
+			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,vy,md2.mesh.x,md2.mesh.y)];
+			spct=[spct InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,thickness,md2.mesh.x,md2.mesh.y)];
+			steps=[steps t*md1.timestepping.time_step];
+		end
+		count=count+1;
+	end
+
+	%As long as there are recorded time steps, spc the boundaries with velocities
+	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.prognostic.spcthickness=md2.prognostic.spcthickness*ones(1,size(spct,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.prognostic.spcthickness(find(md2.mesh.vertexonboundary),:)=spct(find(md2.mesh.vertexonboundary),:);
+		md2.diagnostic.spcvx=[md2.diagnostic.spcvx; steps];
+		md2.diagnostic.spcvy=[md2.diagnostic.spcvy; steps];
+		md2.diagnostic.spcvz=[md2.diagnostic.spcvz; steps];
+		md2.prognostic.spcthickness=[md2.prognostic.spcthickness; steps];
+	end
+
+	%Diagnostic.  Don't spc the icefront vertices.
+	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;
+		md2.prognostic.spcthickness(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
+	end
+
+	%Clear results fields
+	if isstruct(md1.results),
+		md2.results=[];
+	end
+
Index: issm/trunk-jpl/src/m/model/regionaltransient2d.m
===================================================================
--- issm/trunk-jpl/src/m/model/regionaltransient2d.m	(revision 13007)
+++ 	(revision )
@@ -1,158 +1,0 @@
-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(mde,'hmin',hmin,'hmax',hmax,'field',[mde.inversion.vel_obs mde.geometry.surface],'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);
-
-		%copy field, interpolated to new mesh
-		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
-
-	%Read transient velocities and thickness, looping through only the populated times
-	spcx=[];
-	spcy=[];
-	spct=[];
-	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);
-			thickness=PatchToVec(md1.results.TransientSolution(t).Thickness);
-			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,vy,md2.mesh.x,md2.mesh.y)];
-			spct=[spct InterpFromMeshToMesh2d(md1.mesh.elements,md1.mesh.x,md1.mesh.y,thickness,md2.mesh.x,md2.mesh.y)];
-			steps=[steps t*md1.timestepping.time_step];
-		end
-		count=count+1;
-	end
-
-	%As long as there are recorded time steps, spc the boundaries with velocities
-	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.prognostic.spcthickness=md2.prognostic.spcthickness*ones(1,size(spct,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.prognostic.spcthickness(find(md2.mesh.vertexonboundary),:)=spct(find(md2.mesh.vertexonboundary),:);
-		md2.diagnostic.spcvx=[md2.diagnostic.spcvx; steps];
-		md2.diagnostic.spcvy=[md2.diagnostic.spcvy; steps];
-		md2.diagnostic.spcvz=[md2.diagnostic.spcvz; steps];
-		md2.prognostic.spcthickness=[md2.prognostic.spcthickness; steps];
-	end
-
-	%Diagnostic.  Don't spc the icefront vertices.
-	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;
-		md2.prognostic.spcthickness(unique(md2.diagnostic.icefront(:,1:2)),:)=nan;
-	end
-
-	%Clear results fields
-	if isstruct(md1.results),
-		md2.results=[];
-	end
-
Index: issm/trunk-jpl/src/m/model/shear2d.m
===================================================================
--- issm/trunk-jpl/src/m/model/shear2d.m	(revision 13007)
+++ 	(revision )
@@ -1,23 +1,0 @@
-function [sx,sy,sxy,s]=shear2d(md)
-%SHEAR2D - computes 2d strain rate
-%
-%   This routine computes the strain rate of 2d models
-%
-%   Usage:
-%      [sx,sy,sxy,s]=shear2d(md);
-%      s=shear2d(md);
-
-[alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); 
-
-summation=[1;1;1];
-sx=(md.initialization.vx(md.mesh.elements).*alpha)*summation;
-uy=(md.initialization.vx(md.mesh.elements).*beta)*summation;
-vx=(md.initialization.vy(md.mesh.elements).*alpha)*summation;
-sy=(md.initialization.vy(md.mesh.elements).*beta)*summation;						
-sxy=(uy+vx)/2;
-s=sqrt(sx.^2+sy.^2+sxy.^2+sx.*sy);
-
-%if user requested only one output, it must be the norm
-if nargout==1,
-	sx=s;
-end
Index: issm/trunk-jpl/src/m/model/sia.m
===================================================================
--- issm/trunk-jpl/src/m/model/sia.m	(revision 13007)
+++ 	(revision )
@@ -1,26 +1,0 @@
-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.mesh.dimension~=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.geometry.thickness(md.mesh.elements)*summer/3;
-Bel=md.B(md.mesh.elements)*summer/3;
-
-Ael=Bel.^(-3);
-
-velx=-2*(md.materials.rho_ice*md.constants.g)^3*s.^2.*sx.*Ael/4.*hel.^4;
-vely=-2*(md.materials.rho_ice*md.constants.g)^3*s.^2.*sy.*Ael/4.*hel.^4;
-vel=sqrt(velx.^2+vely.^2);
Index: issm/trunk-jpl/src/m/model/slope.m
===================================================================
--- issm/trunk-jpl/src/m/model/slope.m	(revision 13007)
+++ 	(revision )
@@ -1,32 +1,0 @@
-function [sx,sy,s]=slope(md)
-%SLOPE - compute the surface slope
-%
-%   Usage:
-%      [sx,sy,s]=slope(md)
-
-%load some variables (it is much faster if the variab;es are loaded from md once for all) 
-if (md.mesh.dimension==2),
-	numberofelements=md.mesh.numberofelements;
-	numberofnodes=md.mesh.numberofvertices;
-	index=md.mesh.elements;
-	x=md.mesh.x; y=md.mesh.y;
-else
-	numberofelements=md.mesh.numberofelements2d;
-	numberofnodes=md.mesh.numberofvertices2d;
-	index=md.mesh.elements2d;
-	x=md.mesh.x2d; y=md.mesh.y2d;
-end
-
-%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
-[alpha beta]=GetNodalFunctionsCoeff(index,x,y);
-
-summation=[1;1;1];
-sx=(md.geometry.surface(index).*alpha)*summation;
-sy=(md.geometry.surface(index).*beta)*summation;
-s=sqrt(sx.^2+sy.^2);
-
-if md.mesh.dimension==3,
-	sx=project3d(md,'vector',sx,'type','element');
-	sy=project3d(md,'vector',sy,'type','element');
-	s=sqrt(sx.^2+sy.^2);
-end
Index: issm/trunk-jpl/src/m/model/thicknessevolution.m
===================================================================
--- issm/trunk-jpl/src/m/model/thicknessevolution.m	(revision 13007)
+++ 	(revision )
@@ -1,28 +1,0 @@
-function dhdt=thicknessevolution(md)
-%THICKNESSEVOLUTION - compute the new thickness of a model after ∆t
-%
-%   This routine compute the new thickness of a model after a time step
-%   according to the following formula:
-%   dh/dt=-div(Hu)
-%
-%   Usage:
-%      dhdt=thicknessevolution(md)
-
-if (length(md.initialization.vx)~=md.mesh.numberofvertices)|(length(md.initialization.vy)~=md.mesh.numberofvertices)
-	error(['thicknessevolution error message: vx and vy should have a length of ' num2str(md.mesh.numberofvertices)])
-end
-
-%load some variables 
-H=md.geometry.thickness;
-vx=md.initialization.vx;
-vy=md.initialization.vy;
-index=md.mesh.elements;
-
-%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
-[alpha beta]=GetNodalFunctionsCoeff(md.mesh.elements,md.mesh.x,md.mesh.y); 
-
-%compute dhdt=div(Hu)
-summation=1/3*ones(3,1);
-dhdt=(vx(index)*summation).*sum( H(index).*alpha,2) + (vy(index)*summation).*sum(H(index).*beta,2) ...
-	+ ( H(index)*summation).*sum(vx(index).*alpha,2) + ( H(index)*summation).*sum(vy(index).*beta,2);
-dhdt=-dhdt;
Index: issm/trunk-jpl/src/m/model/tres.m
===================================================================
--- issm/trunk-jpl/src/m/model/tres.m	(revision 13007)
+++ 	(revision )
@@ -1,119 +1,0 @@
-function md=tres(md,string)
-%TRES - transfer results results to corresponding model fields. 
-%
-%    Usage: md=tres(md,string)
-%
-%    Example: md=tres(md,'diagnostic');
-
-%check number of arguments
-
-if strcmpi(string,'diagnostic'),
-	if md.mesh.dimension==2,
-		md.initialization.vx=md.results.DiagnosticSolution.Vx;
-		md.initialization.vy=md.results.DiagnosticSolution.Vy;
-	else 
-		md.initialization.vx=md.results.DiagnosticSolution.Vx;
-		md.initialization.vy=md.results.DiagnosticSolution.Vy;
-		md.initialization.vz=md.results.DiagnosticSolution.Vz;
-	end
-	md.initialization.vel=md.results.DiagnosticSolution.Vel;
-
-	if isfield(md.results.DiagnosticSolution,'Pressure'),
-		md.initialization.pressure=md.results.DiagnosticSolution.Pressure;
-	end
-	if md.rifts.numrifts,
-		if isfield(md.results.DiagnosticSolution,'riftproperties'),
-			md.rifts.riftproperties=md.results.DiagnosticSolution.riftproperties;
-		end
-	end
-	if md.inversion.iscontrol==1,
-		for control_parameters=md.inversion.control_parameters
-			%Will need to be updated... good luck ;)
-			md.(EnumToModelField(control_parameters))=md.results.DiagnosticSolution.(EnumToString(control_parameters));
-		end
-	end
-
-elseif strcmpi(string,'dakota'),
-	md.qmu.results=md.results.dakota;
-
-elseif strcmpi(string,'flaim'),
-	md.flaim.solution=md.results.FlaimSolution.solution;
-	md.flaim.quality =md.results.FlaimSolution.quality;
-
-elseif strcmpi(string,'transient'),
-	results=md.results.TransientSolution;
-	results2.Vel=NaN;
-	count=1;
-	for i=1:length(results),
-		if ~isempty(md.results.TransientSolution(i).Vel),
-			results2(count).Vel=md.results.TransientSolution(i).Vel;
-			results2(count).Surface=md.results.TransientSolution(i).Surface;
-			results2(count).Thickness=md.results.TransientSolution(i).Thickness;
-			results2(count).Bed=md.results.TransientSolution(i).Bed;
-			results2(count).Vx=md.results.TransientSolution(i).Vx;
-			results2(count).Vy=md.results.TransientSolution(i).Vy;
-			results2(count).time=md.results.TransientSolution(i).time;
-			results2(count).step=md.results.TransientSolution(i).step;
-			if ~strcmpi(md.groundingline.migration,'None'),
-				results2(count).ElementOnIceShelf=md.results.TransientSolution(i).ElementOnIceShelf;
-			end
-			count=count+1;
-		end
-	end
-	md.results.TransientSolution=results2;
-	clear results,results2;
-elseif strcmpi(string,'steadystate'),
-	md.initialization.vx=md.results.SteadystateSolution.Vx;
-	md.initialization.vy=md.results.SteadystateSolution.Vy;
-	if isfield(md.results.SteadystateSolution,'Vz'),
-		md.initialization.vz=md.results.SteadystateSolution.Vz;
-	end
-
-	md.initialization.vel=md.results.SteadystateSolution.Vel;
-	md.initialization.pressure=md.results.SteadystateSolution.Pressure;
-	md.initialization.temperature=md.results.SteadystateSolution.Temperature;
-	md.basalforcings.melting_rate=md.results.SteadystateSolution.BasalforcingsMeltingRate;
-
-	if md.inversion.iscontrol==1,
-		for control_parameters=md.inversion.control_parameters
-			md.(EnumToModelField(control_parameters))=md.results.SteadystateSolution.(EnumToString(control_parameters));
-		end
-	end
-
-elseif strcmpi(string,'thermal'),
-	md.initialization.temperature=md.results.ThermalSolution.Temperature;
-	md.basalforcings.melting_rate=md.results.ThermalSolution.BasalMeltingRate;
-elseif strcmpi(string,'hydrology'),
-	md.initialization.watercolumn=md.results.HydrologySolution.Watercolumn;
-
-else 
-	error(['tres error message: analysis ' string ' not supported yet!']);
-end
-end 
-function string=EnumToModelField(enum) % {{{
-	%ENUMTOMODELFIELD - output string of model field associated to enum
-	%
-	%   Usage:
-	%      string=EnumToModelField(enum)
-
-	disp('Warning: EnumToModelField is deprecated, it cannot work with new model definition. This function will be removed in the future');
-
-	switch enum,
-
-		case ThicknessEnum(), string='thickness'; return
-		case FrictionCoefficientEnum(), string='drag_coefficient'; return
-		case MaterialsRheologyBEnum(), string='rheology_B'; return
-		case MaterialsRheologyBbarEnum(), string='rheology_B'; return
-		case BalancethicknessThickeningRateEnum: string='dhdt'; return
-		case VxEnum(), string='vx'; return
-		case InversionVxObsEnum(), string='vx_obs'; return
-		case VyEnum(), string='vy'; return
-		case InversionVyObsEnum(), string='vy_obs'; return
-		case BasalforcingsMeltingRateEnum(), string='basal_melting_rate'; return
-		case SurfaceforcingsAccumulationRateEnum(), string='surface_accumulation_rate'; return
-		case SurfaceforcingsAblationRateEnum(), string='surface_ablation_rate'; return
-		case SurfaceforcingsMassBalanceEnum(), string='surface_mass_balance'; return
-		otherwise, error(['Enum ' num2str(enum)  ' not found associated to any model field']);
-
-		end
-	end % }}}
