Index: /issm/trunk-jpl/src/m/contrib/hydrology/effectivepressure.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/hydrology/effectivepressure.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/contrib/hydrology/effectivepressure.m	(revision 13007)
@@ -0,0 +1,12 @@
+function Neff=effectivepressure(md)
+%EFFECTIVEPRESSURE - compute effective pressure
+%
+%   Usage:
+%      Neff=effectivepressure(md)
+%
+%   Example:
+%      Neff=effectivepressure(md)
+
+Neff=md.materials.rho_ice*md.constants.g*md.geometry.thickness+md.materials.rho_ice*md.constants.g*md.geometry.bed;
+pos=find(Neff<0);
+Neff(pos)=0;
Index: /issm/trunk-jpl/src/m/contrib/massbalance/divergence.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/massbalance/divergence.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/contrib/massbalance/divergence.m	(revision 13007)
@@ -0,0 +1,27 @@
+function div=divergence(md,a,b)
+%DIVERGENCE - divergence of [a;b] vector, using model's triangulation.
+%
+%   Usage:
+%      div=divergence(md,a,b)
+
+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; z=md.mesh.z;
+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];
+dx=(a(index).*alpha)*summation;
+dy=(b(index).*beta)*summation;
+
+div=dx+dy;
+div=averaging(md,div,1);
Index: /issm/trunk-jpl/src/m/kml/kmlimagesc.m
===================================================================
--- /issm/trunk-jpl/src/m/kml/kmlimagesc.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/kml/kmlimagesc.m	(revision 13007)
@@ -0,0 +1,67 @@
+function kmlimagesc(md,fieldname,varargin)
+%KMLIMAGESC - create lat,long kml image
+%
+%   Usage:
+%      kmlimagesc(md,field,options);
+%
+%   Options: 
+%      'hemisphere': default +1;
+%      'central_meridian: 45 for Greenland and 0 for Antarctica
+%      'standard_parallel: 70 for Greenland and 71 for Antarctica
+%      'posting': default .1 degree
+%
+
+%process varargin for options: 
+options=pairoptions(varargin{:});
+
+%recover field: 
+field=md.(fieldname);
+
+%recover some options, and set defaults
+fontsize=getfieldvalue(options,'fontsize',12);
+posting=getfieldvalue(options,'posting',.1);
+minlong=getfieldvalue(options,'minlong',min(md.mesh.long));
+maxlong=getfieldvalue(options,'maxlong',max(md.mesh.long));
+minlat=getfieldvalue(options,'minlat',min(md.mesh.lat));
+maxlat=getfieldvalue(options,'maxlat',max(md.mesh.lat));
+minfield=getfieldvalue(options,'minfield',min(field));
+maxfield=getfieldvalue(options,'maxfield',max(field));
+
+%do we have hemisphere setup?:
+if ~isstr(md.mesh.hemisphere),
+	error('md.mesh.hemisphere should be ''s'' or ''n''');
+end
+
+if strcmpi(md.mesh.hemisphere,'s'),
+	hemisphere=1;
+	central_meridian=getfieldvalue(options,'central_meridian',45);
+	standard_parallel=getfieldvalue(options,'standard_parallel',70);
+elseif strcmpi(md.mesh.hemisphere,'n'),
+	hemisphere=-1;
+	central_meridian=getfieldvalue(options,'central_meridian',0);
+	standard_parallel=getfieldvalue(options,'standard_parallel',71);
+else
+	error('md.mesh.hemisphere should be ''s'' or ''n''');
+end
+
+%figure out nlines and ncols in our image
+nlines=(maxlat-minlat)/posting;
+ncols=(maxlong-minlong)/posting;
+
+%regrid to lat,long grid
+[x_m,y_m,field]=InterpFromMeshToGrid(md.mesh.elements,md.mesh.long,md.mesh.lat,field,minlong,maxlat,posting,posting,nlines,ncols,NaN);
+field=flipud(field);
+
+%massage  and log:
+pos=find(field<minfield); field(pos)=minfield;
+pos=find(field>maxfield);field(pos)=maxfield;
+
+%create google earth kml file out of this regridded dataset:
+imagestr=ge_imagesc(x_m,y_m,field,'imgURL',[fieldname '.png'],'name',fieldname);
+imagestr=ge_folder(fieldname,imagestr);
+colorbarstr=ge_colorbar((min(x_m)+max(x_m))/2,(min(y_m)+max(y_m))/2,field,'name',fieldname);
+colorbarstr=ge_folder('Colorbar',colorbarstr);
+ge_output([fieldname '.kml'],[imagestr colorbarstr]);
+
+%now, create kmz file:
+system(['mv ' [fieldname '.kml'] ' doc.kml && zip ' [fieldname '.kmz'] ' doc.kml ' fieldname '.png && rm -rf doc.kml ' [fieldname '.png'] ]);
Index: /issm/trunk-jpl/src/m/meca/drivingstress.m
===================================================================
--- /issm/trunk-jpl/src/m/meca/drivingstress.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/meca/drivingstress.m	(revision 13007)
@@ -0,0 +1,18 @@
+function [px,py,pmag]=drivingstress(md)
+%DRIVINGSTRESS -  evaluates the driving stress
+%
+%   The driving stress is computed according to the following formula: 
+%   driving stress= rho_ice*g*H*slope
+%
+%   Usage:
+%      [Fx,Fy,Fmag]=drivingstress(md)
+
+%Get slope
+[sx,sy,s]=slope(md);
+
+%Average thickness over elements
+thickness_bar=(md.geometry.thickness(md.mesh.elements(:,1))+md.geometry.thickness(md.mesh.elements(:,2))+md.geometry.thickness(md.mesh.elements(:,3)))/3;
+
+px=md.materials.rho_ice*md.constants.g*thickness_bar.*sx;
+py=md.materials.rho_ice*md.constants.g*thickness_bar.*sy;
+pmag=sqrt(px.^2+py.^2);
Index: /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py
===================================================================
--- /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py	(revision 13007)
+++ /issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py	(revision 13007)
@@ -0,0 +1,52 @@
+def oshostname():
+	import socket
+
+	return socket.gethostname().lower().split('.')[0]
+
+def strcmp(s1,s2):
+
+	if s1 == s2:
+		return True
+	else:
+		return False
+
+def strncmp(s1,s2,n):
+
+	if s1[0:n] == s2[0:n]:
+		return True
+	else:
+		return False
+
+def strcmpi(s1,s2):
+
+	if s1.lower() == s2.lower():
+		return True
+	else:
+		return False
+
+def strncmpi(s1,s2,n):
+
+	if s1.lower()[0:n] == s2.lower()[0:n]:
+		return True
+	else:
+		return False
+
+def ismember(a,s):
+	import numpy
+
+	if not isinstance(s,(tuple,list,dict,numpy.ndarray)):
+		s=[s]
+
+	if not isinstance(a,(tuple,list,dict,numpy.ndarray)):
+		a=[a]
+
+	if not isinstance(a,numpy.ndarray):
+		b=[item in s for item in a]
+
+	else:
+		b=numpy.empty_like(a)
+		for i,item in enumerate(a.flat):
+			b.flat[i]=item in s
+
+	return b
+
Index: sm/trunk-jpl/src/m/model/DepthAverage.m
===================================================================
--- /issm/trunk-jpl/src/m/model/DepthAverage.m	(revision 13006)
+++ 	(revision )
@@ -1,33 +1,0 @@
-function  vector_average=DepthAverage(md,vector);
-%DEPTHAVERAGE - computes depth average of 3d vector, and return value on 2d mesh. 
-%
-%   Usage:
-%      vector_average=DepthAverage(md,vector);
-%
-%   Example:
-%      vel_bar=DepthAverage(md,md.initialization.vel);
-
-%check that the model given in input is 3d
-if ~md.mesh.dimension==3;
-	error('DepthAverage error message: the model given in input must be 3d')
-end
-
-%nods data
-if (length(vector)==md.mesh.numberofvertices),
-	vector_average=zeros(md.mesh.numberofvertices2d,1);
-	for i=1:md.mesh.numberoflayers-1,
-		vector_average=vector_average+(project2d(md,vector,i)+project2d(md,vector,i+1))/2.*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i));
-	end
-	vector_average=vector_average./project2d(md,md.geometry.thickness,1);
-
-%element data
-elseif (length(vector)==md.mesh.numberofelements),
-	vector_average=zeros(md.mesh.numberofelements2d,1);
-	for i=1:md.mesh.numberoflayers-1,
-		vector_average=vector_average+project2d(md,vector,i).*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i));
-	end
-	vector_average=vector_average./project2d(md,md.geometry.thickness,1);
-
-else
-	error('vector size not supported yet');
-end
Index: sm/trunk-jpl/src/m/model/EnumToModelField.m
===================================================================
--- /issm/trunk-jpl/src/m/model/EnumToModelField.m	(revision 13006)
+++ 	(revision )
@@ -1,26 +1,0 @@
-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
Index: sm/trunk-jpl/src/m/model/MatlabFuncs.py
===================================================================
--- /issm/trunk-jpl/src/m/model/MatlabFuncs.py	(revision 13006)
+++ 	(revision )
@@ -1,52 +1,0 @@
-def oshostname():
-	import socket
-
-	return socket.gethostname().lower().split('.')[0]
-
-def strcmp(s1,s2):
-
-	if s1 == s2:
-		return True
-	else:
-		return False
-
-def strncmp(s1,s2,n):
-
-	if s1[0:n] == s2[0:n]:
-		return True
-	else:
-		return False
-
-def strcmpi(s1,s2):
-
-	if s1.lower() == s2.lower():
-		return True
-	else:
-		return False
-
-def strncmpi(s1,s2,n):
-
-	if s1.lower()[0:n] == s2.lower()[0:n]:
-		return True
-	else:
-		return False
-
-def ismember(a,s):
-	import numpy
-
-	if not isinstance(s,(tuple,list,dict,numpy.ndarray)):
-		s=[s]
-
-	if not isinstance(a,(tuple,list,dict,numpy.ndarray)):
-		a=[a]
-
-	if not isinstance(a,numpy.ndarray):
-		b=[item in s for item in a]
-
-	else:
-		b=numpy.empty_like(a)
-		for i,item in enumerate(a.flat):
-			b.flat[i]=item in s
-
-	return b
-
Index: sm/trunk-jpl/src/m/model/MatlabProcessPatch.m
===================================================================
--- /issm/trunk-jpl/src/m/model/MatlabProcessPatch.m	(revision 13006)
+++ 	(revision )
@@ -1,65 +1,0 @@
-function structure=MatlabProcessPatch(structure);
-%PROCESSPATCH - create a structure from a patch
-%
-%   Usage:
-%      Result=ProcessPatch(Result);
-
-%return if there is no field Patch
-if (~isfield(structure,'Patch')),
-	return;
-end
-
-%loop over steps
-for i=1:length(structure),
-
-	%Get Patch for current step
-	Patch=structure(i).Patch;
-	numvertices=structure(i).PatchVertices;
-
-	%check that Patch is not empty
-	if length(Patch)==0 continue; end
-
-	%Get number of fields;
-	fields=unique(Patch(:,1));
-	steps=unique(Patch(:,2));
-
-	%parse steps
-	for j=1:length(steps),
-
-		posstep=find(Patch(:,2)==steps(j));
-
-		%Take all the lines of the Patch for this timestep
-		temporarypatch=Patch(posstep,:);
-		time=temporarypatch(1,3);
-		step=temporarypatch(1,2);
-
-		%parse fields
-		for i=1:length(fields),
-
-			%get name
-			fieldname=EnumToString(fields(i));
-
-			%get line positions
-			pos=find(temporarypatch(:,1)==fields(i));
-
-			%Fill Result structure
-			structure(step).steps=step;
-			structure(step).time=time;
-			structure(step).(fieldname).element=temporarypatch(pos,4);
-			structure(step).(fieldname).interpolation=temporarypatch(pos,5);
-			structure(step).(fieldname).index=temporarypatch(pos,6:5+numvertices);
-			if structure(step).(fieldname).interpolation==P1Enum,
-				structure(step).(fieldname).value=temporarypatch(pos,6+numvertices:end);
-			end
-			if structure(step).(fieldname).interpolation==P0Enum,
-				structure(step).(fieldname).value=temporarypatch(pos,6+numvertices);
-			end
-
-		end
-	end
-end
-
-%remove fields
-structure=rmfield(structure,'Patch');
-structure=rmfield(structure,'PatchVertices');
-structure=rmfield(structure,'PatchNodes');
Index: sm/trunk-jpl/src/m/model/MatlabProcessPatch.py
===================================================================
--- /issm/trunk-jpl/src/m/model/MatlabProcessPatch.py	(revision 13006)
+++ 	(revision )
@@ -1,19 +1,0 @@
-def MatlabProcessPatch(structure):
-	"""
-	PROCESSPATCH - create a structure from a patch
- 
-	   Usage:
-	      Result=ProcessPatch(Result);
-	"""
-
-	#loop over steps
-	for structurei in structure.itervalues():
-
-		#return if there is no field Patch
-		if not 'Patch' in structurei:
-			continue
-
-		raise SystemError("MatlabProcessPatch not implemented in Python.")
-
-	return structure
-
Index: sm/trunk-jpl/src/m/model/divergence.m
===================================================================
--- /issm/trunk-jpl/src/m/model/divergence.m	(revision 13006)
+++ 	(revision )
@@ -1,27 +1,0 @@
-function div=divergence(md,a,b)
-%DIVERGENCE - divergence of [a;b] vector, using model's triangulation.
-%
-%   Usage:
-%      div=divergence(md,a,b)
-
-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; z=md.mesh.z;
-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];
-dx=(a(index).*alpha)*summation;
-dy=(b(index).*beta)*summation;
-
-div=dx+dy;
-div=averaging(md,div,1);
Index: sm/trunk-jpl/src/m/model/drivingstress.m
===================================================================
--- /issm/trunk-jpl/src/m/model/drivingstress.m	(revision 13006)
+++ 	(revision )
@@ -1,18 +1,0 @@
-function [px,py,pmag]=drivingstress(md)
-%DRIVINGSTRESS -  evaluates the driving stress
-%
-%   The driving stress is computed according to the following formula: 
-%   driving stress= rho_ice*g*H*slope
-%
-%   Usage:
-%      [Fx,Fy,Fmag]=drivingstress(md)
-
-%Get slope
-[sx,sy,s]=slope(md);
-
-%Average thickness over elements
-thickness_bar=(md.geometry.thickness(md.mesh.elements(:,1))+md.geometry.thickness(md.mesh.elements(:,2))+md.geometry.thickness(md.mesh.elements(:,3)))/3;
-
-px=md.materials.rho_ice*md.constants.g*thickness_bar.*sx;
-py=md.materials.rho_ice*md.constants.g*thickness_bar.*sy;
-pmag=sqrt(px.^2+py.^2);
Index: sm/trunk-jpl/src/m/model/effectivepressure.m
===================================================================
--- /issm/trunk-jpl/src/m/model/effectivepressure.m	(revision 13006)
+++ 	(revision )
@@ -1,12 +1,0 @@
-function Neff=effectivepressure(md)
-%EFFECTIVEPRESSURE - compute effective pressure
-%
-%   Usage:
-%      Neff=effectivepressure(md)
-%
-%   Example:
-%      Neff=effectivepressure(md)
-
-Neff=md.materials.rho_ice*md.constants.g*md.geometry.thickness+md.materials.rho_ice*md.constants.g*md.geometry.bed;
-pos=find(Neff<0);
-Neff(pos)=0;
Index: /issm/trunk-jpl/src/m/model/extrusion/DepthAverage.m
===================================================================
--- /issm/trunk-jpl/src/m/model/extrusion/DepthAverage.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/model/extrusion/DepthAverage.m	(revision 13007)
@@ -0,0 +1,33 @@
+function  vector_average=DepthAverage(md,vector);
+%DEPTHAVERAGE - computes depth average of 3d vector, and return value on 2d mesh. 
+%
+%   Usage:
+%      vector_average=DepthAverage(md,vector);
+%
+%   Example:
+%      vel_bar=DepthAverage(md,md.initialization.vel);
+
+%check that the model given in input is 3d
+if ~md.mesh.dimension==3;
+	error('DepthAverage error message: the model given in input must be 3d')
+end
+
+%nods data
+if (length(vector)==md.mesh.numberofvertices),
+	vector_average=zeros(md.mesh.numberofvertices2d,1);
+	for i=1:md.mesh.numberoflayers-1,
+		vector_average=vector_average+(project2d(md,vector,i)+project2d(md,vector,i+1))/2.*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i));
+	end
+	vector_average=vector_average./project2d(md,md.geometry.thickness,1);
+
+%element data
+elseif (length(vector)==md.mesh.numberofelements),
+	vector_average=zeros(md.mesh.numberofelements2d,1);
+	for i=1:md.mesh.numberoflayers-1,
+		vector_average=vector_average+project2d(md,vector,i).*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i));
+	end
+	vector_average=vector_average./project2d(md,md.geometry.thickness,1);
+
+else
+	error('vector size not supported yet');
+end
Index: /issm/trunk-jpl/src/m/model/extrusion/project2d.m
===================================================================
--- /issm/trunk-jpl/src/m/model/extrusion/project2d.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/model/extrusion/project2d.m	(revision 13007)
@@ -0,0 +1,36 @@
+function projection_value=project2d(md3d,value,layer)
+%PROJECT2D - returns the value of a field for a given layer of the mesh
+%
+%   project 'value' vector taken at layer 'layer' from extruded 2d-3d mesh onto 2d mesh 
+%   used to do the extrusion. This routine is used to compare values between a 2d-3d mesh
+%   at a certain layer, and the equivalent value (if it exists), on the original 2d mesh. 
+%   This routine relies heavily on projections (contained in 3d model md) recored during 
+%   the extrude operation.
+%
+%   Usage:
+%      projection_value=project2d(md3d,value,layer)
+%
+%   Example:
+%      vel2=project2d(md3d,md3d.vel,2);
+
+%some checks on list of arguments
+if ((nargin~=3) ),
+	help project2d
+	error('project2d error message');
+end
+
+if (md3d.mesh.dimension~=3),
+	error('wrong model type ... should be ''3d''');
+end
+
+if ((layer<1) | (layer>md3d.mesh.numberoflayers)),
+	error(['layer must be between 1 and ' num2str(md3d.mesh.numberoflayers)]);
+end
+
+if size(value,1)==md3d.mesh.numberofvertices,
+	projection_value=value((layer-1)*md3d.mesh.numberofvertices2d+1:layer*md3d.mesh.numberofvertices2d,:);
+elseif size(value,1)==md3d.mesh.numberofvertices+1,
+	projection_value=[value((layer-1)*md3d.mesh.numberofvertices2d+1:layer*md3d.mesh.numberofvertices2d,:); value(end,:)];
+else
+	projection_value=value((layer-1)*md3d.mesh.numberofelements2d+1:layer*md3d.mesh.numberofelements2d,:);
+end
Index: /issm/trunk-jpl/src/m/model/extrusion/project3d.m
===================================================================
--- /issm/trunk-jpl/src/m/model/extrusion/project3d.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/model/extrusion/project3d.m	(revision 13007)
@@ -0,0 +1,82 @@
+function projected_vector=project3d(md,varargin);
+%PROJECT3D - vertically project a vector from 2d mesh
+%
+%   vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh.
+%   This vector can be a node vector of size (md.mesh.numberofvertices2d,N/A) or an 
+%   element vector of size (md.mesh.numberofelements2d,N/A). 
+%   arguments: 
+%      'vector': 2d vector
+%      'type': 'element' or 'node'. 
+%   options: 
+%      'layer' a layer number where vector should keep its values. If not specified, all layers adopt the 
+%             value of the 2d vector.
+%      'padding': default to 0 (value adopted by other 3d layers not being projected0
+%
+%   Egs:
+%      extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',NaN);
+%      extruded_vector=project3d(md,'vector',vector2d,'type','element','padding',0);
+%      extruded_vector=project3d(md,'vector',vector2d,'type','node');
+
+%some regular checks
+if nargin==0,
+	help project3d
+	error('bad usage');
+end
+if md.mesh.dimension~=3
+	error('input model is not 3d');
+end
+
+%retrieve parameters from options.
+options      = pairoptions(varargin{:});
+vector2d     = getfieldvalue(options,'vector');     %mandatory
+type         = getfieldvalue(options,'type');       %mandatory
+layer        = getfieldvalue(options,'layer',0);    %optional (do all layers otherwise)
+paddingvalue = getfieldvalue(options,'padding',0);  %0 by default
+
+if length(vector2d)==1,
+	projected_vector=vector2d;
+elseif strcmpi(type,'node'),
+
+	%Initialize 3d vector
+	if size(vector2d,1)==md.mesh.numberofvertices2d
+		projected_vector=paddingvalue*ones(md.mesh.numberofvertices,  size(vector2d,2));
+	elseif size(vector2d,1)==md.mesh.numberofvertices2d+1
+		projected_vector=paddingvalue*ones(md.mesh.numberofvertices+1,size(vector2d,2));
+		projected_vector(end,:)=vector2d(end,:);
+		vector2d=vector2d(1:end-1,:);
+	else
+		error('vector length not supported')
+	end
+
+	%Fill in
+	if layer==0,
+		for i=1:md.mesh.numberoflayers,
+			projected_vector(((i-1)*md.mesh.numberofvertices2d+1):(i*md.mesh.numberofvertices2d),:)=vector2d;
+		end
+	else
+		projected_vector(((layer-1)*md.mesh.numberofvertices2d+1):(layer*md.mesh.numberofvertices2d),:)=vector2d;
+	end
+elseif strcmpi(type,'element'),
+
+	%Initialize 3d vector
+	if size(vector2d,1)==md.mesh.numberofelements2d
+		projected_vector=paddingvalue*ones(md.mesh.numberofelements,  size(vector2d,2));
+	elseif size(vector2d,1)==md.mesh.numberofelements2d+1
+		projected_vector=paddingvalue*ones(md.mesh.numberofelements+1,size(vector2d,2));
+		projected_vector(end,:)=vector2d(end,:);
+		vector2d=vector2d(1:end-1,:);
+	else
+		error('vector length not supported')
+	end
+
+	if layer==0,
+		for i=1:(md.mesh.numberoflayers-1),
+			projected_vector( ((i-1)*md.mesh.numberofelements2d+1):(i*md.mesh.numberofelements2d),:)=vector2d;
+		end
+
+	else
+		projected_vector( ((layer-1)*md.mesh.numberofelements2d+1):(layer*md.mesh.numberofelements2d),:)=vector2d;
+	end
+else
+	error('project3d error message: unknown projection type');
+end
Index: sm/trunk-jpl/src/m/model/graddetection.m
===================================================================
--- /issm/trunk-jpl/src/m/model/graddetection.m	(revision 13006)
+++ 	(revision )
@@ -1,33 +1,0 @@
-function [direction,direction2]=graddetection(md)
-%GRADDETECTION detect gradient of control method between steps nsteps+1 and nsteps
-%
-% Usage: direction=graddetection(md);
-
-%keep copy of md: 
-md2=md;
-
-%solve first batch of control methods, with given settings.
-md2=solve(md2,'solution_type','DiagnosticAnalysis');
-
-%record final optimized parameter. 
-parameter1=md2.results.DiagnosticAnalysis.parameter;
-
-%plug optimized parameter in model. 
-md2.(EnumToModelField(md2.results.DiagnosticAnalysis.inversion.control_parameters))=parameter1;
-
-%put nsteps to 1: 
-md2.nsteps=1;
-md2.optscal=md2.optscal(end)*ones(md2.nsteps,1);
-md2.fit=md2.fit(end)*ones(md2.nsteps,1);
-md2.cm_jump=md2.cm_jump(end)*ones(md2.nsteps,1);
-md2.maxiter=md2.maxiter(end)*ones(md2.nsteps,1);
-
-%rerun control method with optimized parameter, only for 1 more step.
-md2=solve(md2,'solution_type','DiagnosticAnalysis');
-
-%get optimized parameter after 1 more step. 
-parameter2=md2.results.DiagnosticAnalysis.parameter;
-
-%return relative  difference between nsteps+1 and nsteps;
-direction=(parameter2-parameter1)./parameter1;
-direction2=(parameter2-parameter1);
Index: /issm/trunk-jpl/src/m/model/io/loadmodel.m
===================================================================
--- /issm/trunk-jpl/src/m/model/io/loadmodel.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/model/io/loadmodel.m	(revision 13007)
@@ -0,0 +1,46 @@
+function varargout=loadmodel(path)
+%LOADMODEL - load a model using built-in load module
+%
+%   check that model prototype has not changed. if so, adapt to new model prototype.
+%
+%   Usage:
+%      md=loadmodel(path)
+%      loadmodel path
+
+%check nargout
+if nargout>1,
+	error('loadmodel usage error: md=loadmodel(path)');
+end
+
+%check existence
+if exist(path,'file')
+	%do nothing
+elseif exist([path '.mat'],'file')
+	%add extension
+	path = [path '.mat'];
+else
+	error(['loadmodel error message: file ' path ' does not exist']);
+end
+
+try,
+	%recover model on file and name it md
+	warning off MATLAB:unknownElementsNowStruc;
+	warning off MATLAB:load:classNotFound
+	struc=load(path,'-mat');
+	warning on MATLAB:unknownElementsNowStruc;
+	warning on MATLAB:load:classNotFound
+
+	name=char(fieldnames(struc));
+	if size(name,1)>1,
+		error(['loadmodel error message: file ' path ' contains several variables. Only one model should be present.']); 
+	end
+	md=struc.(name);
+	if nargout,
+		varargout{1}=md;
+	else
+		assignin('caller',name,md);
+	end
+catch me
+	disp(getReport(me))
+	error(['could not load model ' path]);
+end
Index: /issm/trunk-jpl/src/m/model/io/loadmodellist.m
===================================================================
--- /issm/trunk-jpl/src/m/model/io/loadmodellist.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/model/io/loadmodellist.m	(revision 13007)
@@ -0,0 +1,50 @@
+function varargout=loadmodellist(path)
+%LOADMODELLIST- load a model using built-in load module
+%
+%   check that modellist prototype has not changed. if so, adapt to new modellist prototype.
+%
+%   Usage:
+%      mds=loadmodellist(path)
+%      loadmodellist path
+
+%check nargout
+if nargout>1,
+	error('loadmodellist usage error: mds=loadmodellist(path)');
+end
+%check existence
+if ~exist(path)
+	error(['loadmodellist error message: file ' path ' does not exist']);
+end
+
+%check that the file is readable
+[stat,mess]=fileattrib(path);
+if( stat==0 | mess.UserRead~=1),
+	error(['loadmodellist error message: file ' path ' is not readable (permission dinied).']);
+end
+
+%check number of variables
+if length(whos('-file',path))>1,
+	error(['loadmodellist error message: file ' path ' contains several variables. Only one model should be present.']);
+end
+
+try,
+	struc=load(path,'-mat');
+
+	%get name of model variable
+	fieldname=char(fieldnames(struc));
+	mds=eval(['struc.' fieldname]);
+	if ~strcmpi(class(mds),'model'),
+		mds2=modellist;
+		mds2=structtomodel(mds2,mds);
+		mds=mds2;
+		clear mds2;
+	end
+	if nargout,
+		varargout{1}=mds;
+	else
+		assignin('caller',fieldname,mds);
+	end
+catch me
+	disp(getReport(me))
+	error(['could not load model ' path]);
+end
Index: sm/trunk-jpl/src/m/model/kmlimagesc.m
===================================================================
--- /issm/trunk-jpl/src/m/model/kmlimagesc.m	(revision 13006)
+++ 	(revision )
@@ -1,67 +1,0 @@
-function kmlimagesc(md,fieldname,varargin)
-%KMLIMAGESC - create lat,long kml image
-%
-%   Usage:
-%      kmlimagesc(md,field,options);
-%   Options: 
-%      'hemisphere': default +1;
-%      'central_meridian: 45 for Greenland and 0 for Antarctica
-%      'standard_parallel: 70 for Greenland and 71 for Antarctica
-%      'posting': default .1 degree
-%
-
-%process varargin for options: 
-options=pairoptions(varargin{:});
-
-%recover field: 
-field=md.(fieldname);
-
-%recover some options, and set defaults
-fontsize=getfieldvalue(options,'fontsize',12);
-posting=getfieldvalue(options,'posting',.1);
-minlong=getfieldvalue(options,'minlong',min(md.mesh.long));
-maxlong=getfieldvalue(options,'maxlong',max(md.mesh.long));
-minlat=getfieldvalue(options,'minlat',min(md.mesh.lat));
-maxlat=getfieldvalue(options,'maxlat',max(md.mesh.lat));
-minfield=getfieldvalue(options,'minfield',min(field));
-maxfield=getfieldvalue(options,'maxfield',max(field));
-
-%do we have hemisphere setup?:
-if ~isstr(md.mesh.hemisphere),
-	error('md.mesh.hemisphere should be ''s'' or ''n''');
-end
-
-if strcmpi(md.mesh.hemisphere,'s'),
-	hemisphere=1;
-	central_meridian=getfieldvalue(options,'central_meridian',45);
-	standard_parallel=getfieldvalue(options,'standard_parallel',70);
-elseif strcmpi(md.mesh.hemisphere,'n'),
-	hemisphere=-1;
-	central_meridian=getfieldvalue(options,'central_meridian',0);
-	standard_parallel=getfieldvalue(options,'standard_parallel',71);
-else
-	error('md.mesh.hemisphere should be ''s'' or ''n''');
-end
-
-%figure out nlines and ncols in our image
-nlines=(maxlat-minlat)/posting;
-ncols=(maxlong-minlong)/posting;
-
-%regrid to lat,long grid
-[x_m,y_m,field]=InterpFromMeshToGrid(md.mesh.elements,md.mesh.long,md.mesh.lat,field,minlong,maxlat,posting,posting,nlines,ncols,NaN);
-field=flipud(field);
-
-%massage  and log:
-pos=find(field<minfield); field(pos)=minfield;
-pos=find(field>maxfield);field(pos)=maxfield;
-
-%create google earth kml file out of this regridded dataset:
-imagestr=ge_imagesc(x_m,y_m,field,'imgURL',[fieldname '.png'],'name',fieldname);
-imagestr=ge_folder(fieldname,imagestr);
-colorbarstr=ge_colorbar((min(x_m)+max(x_m))/2,(min(y_m)+max(y_m))/2,field,'name',fieldname);
-colorbarstr=ge_folder('Colorbar',colorbarstr);
-ge_output([fieldname '.kml'],[imagestr colorbarstr]);
-
-
-%now, create kmz file:
-system(['mv ' [fieldname '.kml'] ' doc.kml && zip ' [fieldname '.kmz'] ' doc.kml ' fieldname '.png && rm -rf doc.kml ' [fieldname '.png'] ]);
Index: sm/trunk-jpl/src/m/model/loadmodel.m
===================================================================
--- /issm/trunk-jpl/src/m/model/loadmodel.m	(revision 13006)
+++ 	(revision )
@@ -1,46 +1,0 @@
-function varargout=loadmodel(path)
-%LOADMODEL - load a model using built-in load module
-%
-%   check that model prototype has not changed. if so, adapt to new model prototype.
-%
-%   Usage:
-%      md=loadmodel(path)
-%      loadmodel path
-
-%check nargout
-if nargout>1,
-	error('loadmodel usage error: md=loadmodel(path)');
-end
-
-%check existence
-if exist(path,'file')
-	%do nothing
-elseif exist([path '.mat'],'file')
-	%add extension
-	path = [path '.mat'];
-else
-	error(['loadmodel error message: file ' path ' does not exist']);
-end
-
-try,
-	%recover model on file and name it md
-	warning off MATLAB:unknownElementsNowStruc;
-	warning off MATLAB:load:classNotFound
-	struc=load(path,'-mat');
-	warning on MATLAB:unknownElementsNowStruc;
-	warning on MATLAB:load:classNotFound
-
-	name=char(fieldnames(struc));
-	if size(name,1)>1,
-		error(['loadmodel error message: file ' path ' contains several variables. Only one model should be present.']); 
-	end
-	md=struc.(name);
-	if nargout,
-		varargout{1}=md;
-	else
-		assignin('caller',name,md);
-	end
-catch me
-	disp(getReport(me))
-	error(['could not load model ' path]);
-end
Index: sm/trunk-jpl/src/m/model/loadmodellist.m
===================================================================
--- /issm/trunk-jpl/src/m/model/loadmodellist.m	(revision 13006)
+++ 	(revision )
@@ -1,50 +1,0 @@
-function varargout=loadmodellist(path)
-%LOADMODELLIST- load a model using built-in load module
-%
-%   check that modellist prototype has not changed. if so, adapt to new modellist prototype.
-%
-%   Usage:
-%      mds=loadmodellist(path)
-%      loadmodellist path
-
-%check nargout
-if nargout>1,
-	error('loadmodellist usage error: mds=loadmodellist(path)');
-end
-%check existence
-if ~exist(path)
-	error(['loadmodellist error message: file ' path ' does not exist']);
-end
-
-%check that the file is readable
-[stat,mess]=fileattrib(path);
-if( stat==0 | mess.UserRead~=1),
-	error(['loadmodellist error message: file ' path ' is not readable (permission dinied).']);
-end
-
-%check number of variables
-if length(whos('-file',path))>1,
-	error(['loadmodellist error message: file ' path ' contains several variables. Only one model should be present.']);
-end
-
-try,
-	struc=load(path,'-mat');
-
-	%get name of model variable
-	fieldname=char(fieldnames(struc));
-	mds=eval(['struc.' fieldname]);
-	if ~strcmpi(class(mds),'model'),
-		mds2=modellist;
-		mds2=structtomodel(mds2,mds);
-		mds=mds2;
-		clear mds2;
-	end
-	if nargout,
-		varargout{1}=mds;
-	else
-		assignin('caller',fieldname,mds);
-	end
-catch me
-	disp(getReport(me))
-	error(['could not load model ' path]);
-end
Index: sm/trunk-jpl/src/m/model/project2d.m
===================================================================
--- /issm/trunk-jpl/src/m/model/project2d.m	(revision 13006)
+++ 	(revision )
@@ -1,36 +1,0 @@
-function projection_value=project2d(md3d,value,layer)
-%PROJECT2D - returns the value of a field for a given layer of the mesh
-%
-%   project 'value' vector taken at layer 'layer' from extruded 2d-3d mesh onto 2d mesh 
-%   used to do the extrusion. This routine is used to compare values between a 2d-3d mesh
-%   at a certain layer, and the equivalent value (if it exists), on the original 2d mesh. 
-%   This routine relies heavily on projections (contained in 3d model md) recored during 
-%   the extrude operation.
-%
-%   Usage:
-%      projection_value=project2d(md3d,value,layer)
-%
-%   Example:
-%      vel2=project2d(md3d,md3d.vel,2);
-
-%some checks on list of arguments
-if ((nargin~=3) ),
-	help project2d
-	error('project2d error message');
-end
-
-if (md3d.mesh.dimension~=3),
-	error('wrong model type ... should be ''3d''');
-end
-
-if ((layer<1) | (layer>md3d.mesh.numberoflayers)),
-	error(['layer must be between 1 and ' num2str(md3d.mesh.numberoflayers)]);
-end
-
-if size(value,1)==md3d.mesh.numberofvertices,
-	projection_value=value((layer-1)*md3d.mesh.numberofvertices2d+1:layer*md3d.mesh.numberofvertices2d,:);
-elseif size(value,1)==md3d.mesh.numberofvertices+1,
-	projection_value=[value((layer-1)*md3d.mesh.numberofvertices2d+1:layer*md3d.mesh.numberofvertices2d,:); value(end,:)];
-else
-	projection_value=value((layer-1)*md3d.mesh.numberofelements2d+1:layer*md3d.mesh.numberofelements2d,:);
-end
Index: sm/trunk-jpl/src/m/model/project3d.m
===================================================================
--- /issm/trunk-jpl/src/m/model/project3d.m	(revision 13006)
+++ 	(revision )
@@ -1,82 +1,0 @@
-function projected_vector=project3d(md,varargin);
-%PROJECT3D - vertically project a vector from 2d mesh
-%
-%   vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh.
-%   This vector can be a node vector of size (md.mesh.numberofvertices2d,N/A) or an 
-%   element vector of size (md.mesh.numberofelements2d,N/A). 
-%   arguments: 
-%      'vector': 2d vector
-%      'type': 'element' or 'node'. 
-%   options: 
-%      'layer' a layer number where vector should keep its values. If not specified, all layers adopt the 
-%             value of the 2d vector.
-%      'padding': default to 0 (value adopted by other 3d layers not being projected0
-%
-%   Egs:
-%      extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',NaN);
-%      extruded_vector=project3d(md,'vector',vector2d,'type','element','padding',0);
-%      extruded_vector=project3d(md,'vector',vector2d,'type','node');
-
-%some regular checks
-if nargin==0,
-	help project3d
-	error('bad usage');
-end
-if md.mesh.dimension~=3
-	error('input model is not 3d');
-end
-
-%retrieve parameters from options.
-options      = pairoptions(varargin{:});
-vector2d     = getfieldvalue(options,'vector');     %mandatory
-type         = getfieldvalue(options,'type');       %mandatory
-layer        = getfieldvalue(options,'layer',0);    %optional (do all layers otherwise)
-paddingvalue = getfieldvalue(options,'padding',0);  %0 by default
-
-if length(vector2d)==1,
-	projected_vector=vector2d;
-elseif strcmpi(type,'node'),
-
-	%Initialize 3d vector
-	if size(vector2d,1)==md.mesh.numberofvertices2d
-		projected_vector=paddingvalue*ones(md.mesh.numberofvertices,  size(vector2d,2));
-	elseif size(vector2d,1)==md.mesh.numberofvertices2d+1
-		projected_vector=paddingvalue*ones(md.mesh.numberofvertices+1,size(vector2d,2));
-		projected_vector(end,:)=vector2d(end,:);
-		vector2d=vector2d(1:end-1,:);
-	else
-		error('vector length not supported')
-	end
-
-	%Fill in
-	if layer==0,
-		for i=1:md.mesh.numberoflayers,
-			projected_vector(((i-1)*md.mesh.numberofvertices2d+1):(i*md.mesh.numberofvertices2d),:)=vector2d;
-		end
-	else
-		projected_vector(((layer-1)*md.mesh.numberofvertices2d+1):(layer*md.mesh.numberofvertices2d),:)=vector2d;
-	end
-elseif strcmpi(type,'element'),
-
-	%Initialize 3d vector
-	if size(vector2d,1)==md.mesh.numberofelements2d
-		projected_vector=paddingvalue*ones(md.mesh.numberofelements,  size(vector2d,2));
-	elseif size(vector2d,1)==md.mesh.numberofelements2d+1
-		projected_vector=paddingvalue*ones(md.mesh.numberofelements+1,size(vector2d,2));
-		projected_vector(end,:)=vector2d(end,:);
-		vector2d=vector2d(1:end-1,:);
-	else
-		error('vector length not supported')
-	end
-
-	if layer==0,
-		for i=1:(md.mesh.numberoflayers-1),
-			projected_vector( ((i-1)*md.mesh.numberofelements2d+1):(i*md.mesh.numberofelements2d),:)=vector2d;
-		end
-
-	else
-		projected_vector( ((layer-1)*md.mesh.numberofelements2d+1):(layer*md.mesh.numberofelements2d),:)=vector2d;
-	end
-else
-	error('project3d error message: unknown projection type');
-end
Index: /issm/trunk-jpl/src/m/model/solve/MatlabProcessPatch.m
===================================================================
--- /issm/trunk-jpl/src/m/model/solve/MatlabProcessPatch.m	(revision 13007)
+++ /issm/trunk-jpl/src/m/model/solve/MatlabProcessPatch.m	(revision 13007)
@@ -0,0 +1,65 @@
+function structure=MatlabProcessPatch(structure);
+%PROCESSPATCH - create a structure from a patch
+%
+%   Usage:
+%      Result=ProcessPatch(Result);
+
+%return if there is no field Patch
+if (~isfield(structure,'Patch')),
+	return;
+end
+
+%loop over steps
+for i=1:length(structure),
+
+	%Get Patch for current step
+	Patch=structure(i).Patch;
+	numvertices=structure(i).PatchVertices;
+
+	%check that Patch is not empty
+	if length(Patch)==0 continue; end
+
+	%Get number of fields;
+	fields=unique(Patch(:,1));
+	steps=unique(Patch(:,2));
+
+	%parse steps
+	for j=1:length(steps),
+
+		posstep=find(Patch(:,2)==steps(j));
+
+		%Take all the lines of the Patch for this timestep
+		temporarypatch=Patch(posstep,:);
+		time=temporarypatch(1,3);
+		step=temporarypatch(1,2);
+
+		%parse fields
+		for i=1:length(fields),
+
+			%get name
+			fieldname=EnumToString(fields(i));
+
+			%get line positions
+			pos=find(temporarypatch(:,1)==fields(i));
+
+			%Fill Result structure
+			structure(step).steps=step;
+			structure(step).time=time;
+			structure(step).(fieldname).element=temporarypatch(pos,4);
+			structure(step).(fieldname).interpolation=temporarypatch(pos,5);
+			structure(step).(fieldname).index=temporarypatch(pos,6:5+numvertices);
+			if structure(step).(fieldname).interpolation==P1Enum,
+				structure(step).(fieldname).value=temporarypatch(pos,6+numvertices:end);
+			end
+			if structure(step).(fieldname).interpolation==P0Enum,
+				structure(step).(fieldname).value=temporarypatch(pos,6+numvertices);
+			end
+
+		end
+	end
+end
+
+%remove fields
+structure=rmfield(structure,'Patch');
+structure=rmfield(structure,'PatchVertices');
+structure=rmfield(structure,'PatchNodes');
Index: /issm/trunk-jpl/src/m/model/solve/MatlabProcessPatch.py
===================================================================
--- /issm/trunk-jpl/src/m/model/solve/MatlabProcessPatch.py	(revision 13007)
+++ /issm/trunk-jpl/src/m/model/solve/MatlabProcessPatch.py	(revision 13007)
@@ -0,0 +1,19 @@
+def MatlabProcessPatch(structure):
+	"""
+	PROCESSPATCH - create a structure from a patch
+ 
+	   Usage:
+	      Result=ProcessPatch(Result);
+	"""
+
+	#loop over steps
+	for structurei in structure.itervalues():
+
+		#return if there is no field Patch
+		if not 'Patch' in structurei:
+			continue
+
+		raise SystemError("MatlabProcessPatch not implemented in Python.")
+
+	return structure
+
Index: /issm/trunk-jpl/src/m/model/tres.m
===================================================================
--- /issm/trunk-jpl/src/m/model/tres.m	(revision 13006)
+++ /issm/trunk-jpl/src/m/model/tres.m	(revision 13007)
@@ -90,2 +90,30 @@
 	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 % }}}