Index: issm/trunk-jpl/src/m/classes/modellist.m
===================================================================
--- issm/trunk-jpl/src/m/classes/modellist.m	(revision 13011)
+++ issm/trunk-jpl/src/m/classes/modellist.m	(revision 13012)
@@ -292,2 +292,96 @@
 	end
 end
+
+function BuildMultipleQueueingScript(cluster,name,executionpath,codepath)% {{{
+%BUILDMULTIPLEQUEUEINGSCRIPT - 
+%
+%   Usage:
+%      BuildMultipleQueueingScript(executionpath,codepath)
+
+disp('building queueing script');
+
+%First try and figure out if there is a special script for this particular cluster
+function_name=['BuildMultipleQueueingScript' cluster];
+
+%some specific treatment of identical cluster, gemini, castor and pollux
+if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
+	function_name='BuildMultipleQueueingScriptgemini';
+end
+
+if exist(function_name,'file'),
+	%Call this function:
+	eval([function_name '(name,executionpath,codepath);']);
+else
+	%Call the generic BuildQueueingScript:
+	BuildMultipleQueueingScriptGeneric(name,executionpath,codepath);
+end
+end % }}}
+function BuildQueueingScriptgemini(name,executionpath,codepath)% {{{
+%BUILDQUEUEINGSCRIPTGEMINI - ...
+%
+%   Usage:
+%      BuildQueueingScriptgemini(md,executionpath,codepath)
+
+scriptname=[name '.queue'];
+
+fid=fopen(scriptname,'w');
+if fid==-1,
+	error(['BuildQueueingScriptgeminierror message: could not open ' scriptname ' file for ascii writing']);
+end
+
+fprintf(fid,'#!/bin/sh\n');
+fprintf(fid,'cd %s\n',executionpath);
+fprintf(fid,'mkdir %s\n',name);
+fprintf(fid,'cd %s\n',name);
+fprintf(fid,'mv ../ModelList.tar.gz ./\n');
+fprintf(fid,'tar -zxvf ModelList.tar.gz\n');
+fprintf(fid,'foreach i (%s-*vs*.queue)\n',name);
+fprintf(fid,'qsub $i\n');
+fprintf(fid,'end\n');
+fclose(fid);
+end% }}}
+function LaunchMultipleQueueJob(cluster,name,executionpath)% {{{
+%LAUNCHMULTIPLEQUEUEJOB - ...
+%
+%   Usage:
+%      LaunchMultipleQueueJob(executionpath)
+
+%First try and figure out if there is a special script for thie particular cluster
+function_name=['LaunchMultipleQueueJob' cluster];
+
+%some specific treatment of identical cluster, gemini, castor and pollux
+if strcmpi(cluster,'castor') || strcmpi(cluster,'pollux'),
+	function_name='LaunchMultipleQueueJobgemini';
+end
+
+if exist(function_name,'file'),
+	%Call this function:
+	eval([function_name '(cluster,name,executionpath);']);
+else
+	%Call the generic LaunchMultipleQueueJob:
+	LaunchMultipleQueueJobGeneric(cluster,name,executionpath);
+end
+end% }}}
+function md=LaunchMultipleQueueJobgemini(cluster,name,executionpath)% {{{
+%LAUNCHMULTIPLEQUEUEJOBGEMINI - Launch multiple queueing script on Gemini cluster
+%
+%   Usage:
+%      LaunchMultipleQueueJobgemini(cluster,name,executionpath)
+
+
+%first, check we have the binary file and the queueing script
+if ~exist([ name '.queue'],'file'),
+	error('LaunchMultipleQueueJobgemini error message: queueing script issing, cannot go forward');
+end
+
+if ~exist('ModelList.tar.gz','file'),
+	error('LaunchMultipleQueueJobgemini error message: inputs models file missing, cannot go forward');
+end
+
+%upload both files to cluster
+disp('uploading input file,  queueing script and variables script');
+eval(['!scp ModelList.tar.gz ' name '.queue '  cluster ':' executionpath]);
+
+disp('launching solution sequence on remote cluster');
+issmssh(cluster,login,['"cd ' executionpath ' && source ' name '.queue "']);
+end% }}}
Index: issm/trunk-jpl/src/m/contrib/bamg/BamgCall.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/bamg/BamgCall.m	(revision 13012)
+++ issm/trunk-jpl/src/m/contrib/bamg/BamgCall.m	(revision 13012)
@@ -0,0 +1,84 @@
+function md=BamgCall(md,field,hmin,hmax,gradation,epsilon),
+%BAMGCALL - call bam
+%
+%   build a metric using the Hessian of the given field
+%   call Bamg and the output mesh is plugged onto the model
+%   -hmin = minimum edge length (m)
+%   -hmax = maximum edge length (m)
+%   -gradation = maximum edge length gradation between 2 elements
+%   -epsilon = average error on each element (m/yr)
+%
+%   Usage:
+%      md=BamgCall(md,field,hmin,hmax,gradation,epsilon);
+%
+%   Example:
+%      md=BamgCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
+
+%2d geometric parameter (do not change)
+scale=2/9; 
+
+%Compute Hessian
+t1=clock; fprintf('%s','      computing Hessian...');
+hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%Compute metric
+t1=clock; fprintf('%s','      computing metric...');
+if length(md.nodeonwater)==md.mesh.numberofvertices,
+	pos=find(md.nodeonwater);
+else
+	pos=[];
+end
+metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%write files
+t1=clock; fprintf('%s','      writing initial mesh files...');
+fid=fopen('carre0.met','w');
+fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
+fprintf(fid,'%i %i %i\n',metric');
+fclose(fid);
+
+fid=fopen('carre0.mesh','w');
+
+%initialiation
+fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
+
+%dimension
+fprintf(fid,'\n%s\n%i\n','Dimension',2);
+
+%Vertices
+fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
+fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
+
+%Triangles
+fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
+fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
+numberofelements1=md.mesh.numberofelements;
+
+%close
+fclose(fid);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%call bamg
+fprintf('%s\n','      call Bamg...');
+system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
+
+%plug new mesh
+t1=clock; fprintf('\n%s','      reading final mesh files...');
+A=meshread('carre1.mesh');
+md.mesh.x=A.x;
+md.mesh.y=A.y;
+md.z=zeros(A.nods,1);
+md.mesh.elements=A.index;
+md.mesh.numberofvertices=A.nods;
+md.mesh.numberofelements=A.nels;
+numberofelements2=md.mesh.numberofelements;
+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%display number of elements
+fprintf('\n%s %i','      inital number of elements:',numberofelements1);
+fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
+
+%clean up:
+system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');
Index: issm/trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m	(revision 13012)
+++ issm/trunk-jpl/src/m/contrib/bamg/BamgCallFromMetric.m	(revision 13012)
@@ -0,0 +1,65 @@
+function md=BamgCallFromMetric(md,metric,gradation),
+%BAMGCALL - call bam
+%
+%   call Bamg and the output mesh is plugged onto the model
+%   -gradation = maximum edge length gradation between 2 elements
+%
+%   Usage:
+%      md=BamgCallFromMetric(md,metric,gradation);
+%
+%   Example:
+%      md=BamgCall(md,metric,1500,10^8,1.3,0.9);
+
+%2d geometric parameter (do not change)
+scale=2/9; 
+
+%write files
+t1=clock; fprintf('%s','      writing initial mesh files...');
+fid=fopen('carre0.met','w');
+fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
+fprintf(fid,'%i %i %i\n',metric');
+fclose(fid);
+
+fid=fopen('carre0.mesh','w');
+
+%initialiation
+fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
+
+%dimension
+fprintf(fid,'\n%s\n%i\n','Dimension',2);
+
+%Vertices
+fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
+fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
+
+%Triangles
+fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
+fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
+numberofelements1=md.mesh.numberofelements;
+
+%close
+fclose(fid);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%call bamg
+fprintf('%s\n','      call Bamg...');
+system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
+
+%plug new mesh
+t1=clock; fprintf('\n%s','      reading final mesh files...');
+A=meshread('carre1.mesh');
+md.mesh.x=A.x;
+md.mesh.y=A.y;
+md.z=zeros(A.nods,1);
+md.mesh.elements=A.index;
+md.mesh.numberofvertices=A.nods;
+md.mesh.numberofelements=A.nels;
+numberofelements2=md.mesh.numberofelements;
+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%display number of elements
+fprintf('\n%s %i','      inital number of elements:',numberofelements1);
+fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
+
+%clean up:
+system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');
Index: issm/trunk-jpl/src/m/contrib/bamg/YamsCall.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/bamg/YamsCall.m	(revision 13012)
+++ issm/trunk-jpl/src/m/contrib/bamg/YamsCall.m	(revision 13012)
@@ -0,0 +1,104 @@
+function md=YamsCall(md,field,hmin,hmax,gradation,epsilon),
+%YAMSCALL - call yams
+%
+%   build a metric using the Hessian of the given field
+%   call Yams and the output mesh is plugged onto the model
+%   -hmin = minimum edge length (m)
+%   -hmax = maximum edge length (m)
+%   -gradation = maximum edge length gradation between 2 elements
+%   -epsilon = average error on each element (m/yr)
+%
+%   Usage:
+%      md=YamsCall(md,field,hmin,hmax,gradation,epsilon);
+%
+%   Example:
+%      md=YamsCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
+
+%2d geometric parameter (do not change)
+scale=2/9; 
+
+%Compute Hessian
+t1=clock; fprintf('%s','      computing Hessian...');
+hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%Compute metric
+t1=clock; fprintf('%s','      computing metric...');
+if length(md.mask.vertexonwater)==md.mesh.numberofvertices,
+	pos=find(md.mask.vertexonwater);
+else
+	pos=[];
+end
+metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%write files
+t1=clock; fprintf('%s','      writing initial mesh files...');
+save -ascii carre0.met  metric
+
+fid=fopen('carre0.mesh','w');
+
+%initialiation
+fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1);
+
+%dimension
+fprintf(fid,'\n%s\n%i\n','Dimension',2);
+
+%Vertices
+fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
+fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]');
+
+%Triangles
+fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
+fprintf(fid,'%i %i %i %i\n',[md.mesh.elements zeros(md.mesh.numberofelements,1)]');
+numberofelements1=md.mesh.numberofelements;
+	
+%Deal with rifts
+if ~isnan(md.rifts.riftstruct),
+	
+	%we have the list of triangles that make up the rift. keep those triangles around during refinement.
+	triangles=[];
+	for i=1:size(md.rifts.riftstruct,1),
+		triangles=[triangles md.rifts(i).segments(:,3)'];
+	end
+
+	fprintf(fid,'\n\n%s\n%i\n\n','RequiredTriangles',length(triangles));
+	fprintf(fid,'%i\n',triangles);
+end
+
+%close
+fclose(fid);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%call yams
+fprintf('%s\n','      call Yams...');
+if ispc
+	%windows
+	system(['yams2-win -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
+elseif ismac
+	%Macosx
+	system(['yams2-osx -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
+else
+	%Linux
+	system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
+end
+
+%plug new mesh
+t1=clock; fprintf('\n%s','      reading final mesh files...');
+Tria=load('carre1.tria');
+Coor=load('carre1.coor');
+md.mesh.x=Coor(:,1);
+md.mesh.y=Coor(:,2);
+md.mesh.z=zeros(size(Coor,1),1);
+md.mesh.elements=Tria;
+md.mesh.numberofvertices=size(Coor,1);
+md.mesh.numberofelements=size(Tria,1);
+numberofelements2=md.mesh.numberofelements;
+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%display number of elements
+fprintf('\n%s %i','      inital number of elements:',numberofelements1);
+fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
+
+%clean up:
+system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb');
Index: issm/trunk-jpl/src/m/contrib/bamg/meshread.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/bamg/meshread.m	(revision 13012)
+++ issm/trunk-jpl/src/m/contrib/bamg/meshread.m	(revision 13012)
@@ -0,0 +1,41 @@
+function Struct=meshread(filename);
+
+%some checks
+if ~exist(filename),
+	error(['meshread error message: file ' filename ' not found!']);
+end
+
+fid=fopen(filename,'r');
+
+while (~feof(fid)),
+
+	A=fscanf(fid,'%s',1);
+
+	if strcmp(A,'MeshVersionFormatted');
+		Struct.Version=fscanf(fid,'%s',1);
+
+	elseif strcmp(A,'Dimension'),
+		Struct.Dimension=fscanf(fid,'%i',1);
+
+	elseif strcmp(A,'Vertices'),
+		Struct.nods=fscanf(fid,'%i',1);
+		A=fscanf(fid,'%f %f %f',[3 Struct.nods]);
+		Struct.x=A(1,:)';
+		Struct.y=A(2,:)';
+
+	elseif strcmp(A,'Triangles'),
+		Struct.nels=fscanf(fid,'%i',1);
+		A=fscanf(fid,'%i %i %i',[4 Struct.nels]);
+		Struct.index=A(1:3,:)';
+
+	elseif strcmp(A,'Quadrilaterals'),
+		Struct.nels=fscanf(fid,'%i',1);
+		A=fscanf(fid,'%i %i %i %i',[5 Struct.nels]);
+		Struct.index=A(1:4,:)';
+	else
+		%do nothing
+
+	end
+end
+
+fclose(fid);
Index: issm/trunk-jpl/src/m/contrib/bamg/yams.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/bamg/yams.m	(revision 13012)
+++ issm/trunk-jpl/src/m/contrib/bamg/yams.m	(revision 13012)
@@ -0,0 +1,171 @@
+function md=yams(md,varargin);
+%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
+%
+%   Usage:
+%      md=yams(md,varargin);
+%      where varargin is a lit of paired arguments. 
+%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
+%      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
+%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
+%                          this option is used to minimize the metric on water (no refinement)
+%      optional arguments: 'resolution': initial mesh resolution [m]
+%      optional arguments: 'nsteps': number of steps of mesh adaptation
+%      optional arguments: 'epsilon': average interpolation error wished [m/yr]
+%      optional arguments: 'hmin': minimum edge length
+%      optional arguments: 'hmanx': maximum edge
+%      optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
+%      
+%
+%   Examples:
+%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat');
+%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
+%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
+
+%recover options
+options=pairoptions(varargin{:});
+options=deleteduplicates(options,1);
+
+%recover some fields
+disp('MeshYams Options:')
+domainoutline=getfieldvalue(options,'domainoutline');
+disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
+riftoutline=getfieldvalue(options,'riftoutline','N/A');
+disp(sprintf('   %-15s: ''%s''','riftoutline',riftoutline));
+groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
+disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
+velocities=getfieldvalue(options,'velocities');
+disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
+resolution=getfieldvalue(options,'resolution',5000);
+disp(sprintf('   %-15s: %f','Resolution',resolution));
+nsteps=getfieldvalue(options,'nsteps',6);
+disp(sprintf('   %-15s: %i','nsteps',nsteps));
+gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
+disp(sprintf('   %-15s: %g','gradation',gradation(1)));
+epsilon=getfieldvalue(options,'epsilon',3);
+disp(sprintf('   %-15s: %f','epsilon',epsilon));
+hmin=getfieldvalue(options,'hmin',500);
+disp(sprintf('   %-15s: %f','hmin',hmin));
+hmax=getfieldvalue(options,'hmax',150*10^3);
+disp(sprintf('   %-15s: %f\n','hmax',hmax));
+
+%mesh with initial resolution
+disp('Initial mesh generation...');
+if strcmpi(riftoutline,'N/A');
+	md=setmesh(md,domainoutline,resolution);
+else
+	md=setmesh(md,domainoutline,riftoutline,resolution);
+	md=meshprocessrifts(md,domainoutline);
+end
+disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
+
+%load velocities 
+disp('loading velocities...');
+Names=VelFindVarNames(velocities);
+Vel=load(velocities);
+
+%start mesh adaptation
+for i=1:nsteps,
+	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
+
+	%interpolate velocities onto mesh
+	disp('   interpolating velocities...');
+	if strcmpi(Names.interp,'node'),
+		vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
+		vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
+	else
+		vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
+		vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
+	end
+	field=sqrt(vx_obs.^2+vy_obs.^2);
+
+	%set mask.vertexonwater  field
+	if ~strcmp(groundeddomain,'N/A'),
+		nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
+		md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
+		md.mask.vertexonwater(find(nodeground))=0;
+	else
+		md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
+	end
+
+	%adapt according to velocities
+	disp('   adapting...');
+	md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
+
+	%if we have rifts, we just messed them up, we need to recreate the segments that constitute those 
+	%rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
+	if md.rifts.numrifts, 
+		md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+		md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+		md.mesh.segments=findsegments(md);
+		md=yamsrecreateriftsegments(md);
+	end
+
+end
+	
+disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
+
+%Now, build the connectivity tables for this mesh.
+md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+
+%recreate segments
+md.mesh.segments=findsegments(md);
+md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+
+%Fill in rest of fields:
+md.mesh.z=zeros(md.mesh.numberofvertices,1);
+md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
+md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
+md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
+md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
+if ~strcmp(groundeddomain,'N/A'),
+	nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
+	md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
+	md.mask.vertexonwater(find(nodeground))=0;
+else
+	md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
+end
+if strcmpi(Names.interp,'node'),
+	md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
+	md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
+else
+	md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
+	md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
+end
+md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
+
+%deal with rifts 
+if md.rifts.numrifts,
+	%first, recreate rift segments
+	md=meshyamsrecreateriftsegments(md);
+
+	%using the segments, recreate the penaltypairs
+	for j=1:md.rifts.numrifts,
+		rift=md.rifts.riftstruct(j);
+
+		%build normals and lengths of segments:
+		lengths=sqrt((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))).^2 + (md.mesh.y(rift.segments(:,1))-md.mesh.y(rift.segments(:,2))).^2 );
+		normalsx=cos(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.y(rift.segments(:,1)))));
+		normalsy=sin(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.y(rift.segments(:,1)))));
+
+		%ok, build penaltypairs: 
+		numpenaltypairs=length(rift.segments)/2-1;
+		rift.penaltypairs=zeros(numpenaltypairs,7);
+
+		for i=1:numpenaltypairs,
+			rift.penaltypairs(i,1)=rift.segments(i,2);
+			rift.penaltypairs(i,2)=rift.segments(end-i,2);
+			rift.penaltypairs(i,3)=rift.segments(i,3);
+			rift.penaltypairs(i,4)=rift.segments(end-i,3);
+			rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
+			rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
+			rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
+		end
+		%renormalize norms: 
+		norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
+		rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
+		rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
+
+		md.rifts.riftstruct(j)=rift;
+	end
+end
Index: issm/trunk-jpl/src/m/geometry/GetAreas.m
===================================================================
--- issm/trunk-jpl/src/m/geometry/GetAreas.m	(revision 13012)
+++ issm/trunk-jpl/src/m/geometry/GetAreas.m	(revision 13012)
@@ -0,0 +1,51 @@
+function areas=GetAreas(index,x,y,varargin)
+%GETAREAS - compute areas or volumes of elements
+%
+%   compute areas of triangular elements or volumes 
+%   of pentahedrons
+%
+%   Usage:
+%      areas  =GetAreas(index,x,y);
+%      volumes=GetAreas(index,x,y,z);
+%
+%   Examples:
+%      areas  =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
+%      volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z);
+
+%get number of elements and number of nodes
+nels=size(index,1);
+nods=length(x);
+if nargin==4, z=varargin{1}; end
+
+%some checks
+if nargout~=1 | (nargin~=3 & nargin~=4),
+	help GetAreas
+	error('GetAreas error message: bad usage')
+end
+if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)),
+	error('GetAreas error message: x,y and z do not have the same length')
+end
+if max(index(:))>nods,
+	error(['GetAreas error message: index should not have values above ' num2str(nods) ])
+end
+if (nargin==3 & size(index,2)~=3),
+	error('GetAreas error message: index should have 3 columns for 2d meshes.')
+end
+if (nargin==4 & size(index,2)~=6),
+	error('GetAreas error message: index should have 6 columns for 3d meshes.')
+end
+
+%initialization
+areas=zeros(nels,1);
+x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3));
+y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
+
+%compute the volume of each element
+if nargin==3,
+	%compute the surface of the triangle
+	areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
+else
+	%V=area(triangle)*1/3(z1+z2+z3)
+	thickness=mean(z(index(:,4:6)),2)-mean(z(index(:,1:3)),2);
+	areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))).*thickness;
+end
Index: issm/trunk-jpl/src/m/mesh/BamgCall.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/BamgCall.m	(revision 13011)
+++ 	(revision )
@@ -1,84 +1,0 @@
-function md=BamgCall(md,field,hmin,hmax,gradation,epsilon),
-%BAMGCALL - call bam
-%
-%   build a metric using the Hessian of the given field
-%   call Bamg and the output mesh is plugged onto the model
-%   -hmin = minimum edge length (m)
-%   -hmax = maximum edge length (m)
-%   -gradation = maximum edge length gradation between 2 elements
-%   -epsilon = average error on each element (m/yr)
-%
-%   Usage:
-%      md=BamgCall(md,field,hmin,hmax,gradation,epsilon);
-%
-%   Example:
-%      md=BamgCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
-
-%2d geometric parameter (do not change)
-scale=2/9; 
-
-%Compute Hessian
-t1=clock; fprintf('%s','      computing Hessian...');
-hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%Compute metric
-t1=clock; fprintf('%s','      computing metric...');
-if length(md.nodeonwater)==md.mesh.numberofvertices,
-	pos=find(md.nodeonwater);
-else
-	pos=[];
-end
-metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%write files
-t1=clock; fprintf('%s','      writing initial mesh files...');
-fid=fopen('carre0.met','w');
-fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
-fprintf(fid,'%i %i %i\n',metric');
-fclose(fid);
-
-fid=fopen('carre0.mesh','w');
-
-%initialiation
-fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
-
-%dimension
-fprintf(fid,'\n%s\n%i\n','Dimension',2);
-
-%Vertices
-fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
-fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
-
-%Triangles
-fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
-fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
-numberofelements1=md.mesh.numberofelements;
-
-%close
-fclose(fid);
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%call bamg
-fprintf('%s\n','      call Bamg...');
-system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
-
-%plug new mesh
-t1=clock; fprintf('\n%s','      reading final mesh files...');
-A=meshread('carre1.mesh');
-md.mesh.x=A.x;
-md.mesh.y=A.y;
-md.z=zeros(A.nods,1);
-md.mesh.elements=A.index;
-md.mesh.numberofvertices=A.nods;
-md.mesh.numberofelements=A.nels;
-numberofelements2=md.mesh.numberofelements;
-t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%display number of elements
-fprintf('\n%s %i','      inital number of elements:',numberofelements1);
-fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
-
-%clean up:
-system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');
Index: issm/trunk-jpl/src/m/mesh/BamgCallFromMetric.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/BamgCallFromMetric.m	(revision 13011)
+++ 	(revision )
@@ -1,65 +1,0 @@
-function md=BamgCallFromMetric(md,metric,gradation),
-%BAMGCALL - call bam
-%
-%   call Bamg and the output mesh is plugged onto the model
-%   -gradation = maximum edge length gradation between 2 elements
-%
-%   Usage:
-%      md=BamgCallFromMetric(md,metric,gradation);
-%
-%   Example:
-%      md=BamgCall(md,metric,1500,10^8,1.3,0.9);
-
-%2d geometric parameter (do not change)
-scale=2/9; 
-
-%write files
-t1=clock; fprintf('%s','      writing initial mesh files...');
-fid=fopen('carre0.met','w');
-fprintf(fid,'%i %i\n',md.mesh.numberofvertices,3);
-fprintf(fid,'%i %i %i\n',metric');
-fclose(fid);
-
-fid=fopen('carre0.mesh','w');
-
-%initialiation
-fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
-
-%dimension
-fprintf(fid,'\n%s\n%i\n','Dimension',2);
-
-%Vertices
-fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
-fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)]');
-
-%Triangles
-fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
-fprintf(fid,'%i %i %i %i\n',[md.mesh.elements ones(md.mesh.numberofelements,1)]');
-numberofelements1=md.mesh.numberofelements;
-
-%close
-fclose(fid);
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%call bamg
-fprintf('%s\n','      call Bamg...');
-system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
-
-%plug new mesh
-t1=clock; fprintf('\n%s','      reading final mesh files...');
-A=meshread('carre1.mesh');
-md.mesh.x=A.x;
-md.mesh.y=A.y;
-md.z=zeros(A.nods,1);
-md.mesh.elements=A.index;
-md.mesh.numberofvertices=A.nods;
-md.mesh.numberofelements=A.nels;
-numberofelements2=md.mesh.numberofelements;
-t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%display number of elements
-fprintf('\n%s %i','      inital number of elements:',numberofelements1);
-fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
-
-%clean up:
-system('rm carre0.mesh carre0.met carre1.mesh carre1.mesh.gmsh');
Index: issm/trunk-jpl/src/m/mesh/FixMesh.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/FixMesh.m	(revision 13011)
+++ issm/trunk-jpl/src/m/mesh/FixMesh.m	(revision 13012)
@@ -1,10 +1,10 @@
 function  [index2 x2 y2 value2]=FixMesh(index,x,y,value)
-%FixMesh fix mesh with broken triangles, orphan vertices, etc ...
+% FIXMESH - FixMesh fix mesh with broken triangles, orphan vertices, etc ...
 %
-% Usage: 
-%            [index2 x2 y2 value2]=FixMesh(index,x,y,value)
-%            where index,x,y is a delaunay triangulation, 
-%                  value is a field on the input triangulation, with values at the vertices
-%                  index2,x2,y2,value2 is the repaired triangulation, with new values on new vertices
+%   Usage: 
+%      [index2 x2 y2 value2]=FixMesh(index,x,y,value)
+%      where index,x,y is a delaunay triangulation, 
+%      value is a field on the input triangulation, with values at the vertices
+%      index2,x2,y2,value2 is the repaired triangulation, with new values on new vertices
 %
 %
Index: issm/trunk-jpl/src/m/mesh/GetAreas.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/GetAreas.m	(revision 13011)
+++ 	(revision )
@@ -1,51 +1,0 @@
-function areas=GetAreas(index,x,y,varargin)
-%GETAREAS - compute areas or volumes of elements
-%
-%   compute areas of triangular elements or volumes 
-%   of pentahedrons
-%
-%   Usage:
-%      areas  =GetAreas(index,x,y);
-%      volumes=GetAreas(index,x,y,z);
-%
-%   Examples:
-%      areas  =GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
-%      volumes=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y,md.z);
-
-%get number of elements and number of nodes
-nels=size(index,1);
-nods=length(x);
-if nargin==4, z=varargin{1}; end
-
-%some checks
-if nargout~=1 | (nargin~=3 & nargin~=4),
-	help GetAreas
-	error('GetAreas error message: bad usage')
-end
-if ((length(y)~=nods) | (nargin==4 & length(z)~=nods)),
-	error('GetAreas error message: x,y and z do not have the same length')
-end
-if max(index(:))>nods,
-	error(['GetAreas error message: index should not have values above ' num2str(nods) ])
-end
-if (nargin==3 & size(index,2)~=3),
-	error('GetAreas error message: index should have 3 columns for 2d meshes.')
-end
-if (nargin==4 & size(index,2)~=6),
-	error('GetAreas error message: index should have 6 columns for 3d meshes.')
-end
-
-%initialization
-areas=zeros(nels,1);
-x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3));
-y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
-
-%compute the volume of each element
-if nargin==3,
-	%compute the surface of the triangle
-	areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
-else
-	%V=area(triangle)*1/3(z1+z2+z3)
-	thickness=mean(z(index(:,4:6)),2)-mean(z(index(:,1:3)),2);
-	areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1))).*thickness;
-end
Index: issm/trunk-jpl/src/m/mesh/GetCharacteristicLength.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/GetCharacteristicLength.m	(revision 13011)
+++ 	(revision )
@@ -1,45 +1,0 @@
-function length=GetCharacteristicLength(index,x,y,varargin)
-%GETCHARACTERISTICLENGTH - compute characteristic length for a mesh
-%
-%   compute characteristic lengths of every element of a mesh.
-%
-%   Usage:
-%      length  =GetCharacteristicLength(index,x,y);
-%      length  =GetCharacteristicLength(index,x,y,z);
-%
-%   Examples:
-%      length  =GetCharacteristicLength(md.mesh.elements,md.mesh.x,md.mesh.y);
-%      length  =GetCharacteristicLength(md.mesh.elements,md.mesh.x,md.mesh.y,md.z);
-
-
-%get number of elements and number of nodes
-nels=size(index,1);
-nods=numel(x);
-
-%some checks
-if nargout~=1 | (nargin~=3 & nargin~=4),
-	help GetCharacteristicLength
-	error('GetCharacteristicLength error message: bad usage')
-end
-if ((numel(y)~=nods) | (nargin==4 & numel(z)~=nods)),
-	error('GetCharacteristicLength error message: x,y and z do not have the same length')
-end
-if max(index(:))>nods,
-	error(['GetCharacteristicLength error message: index should not have values above ' num2str(nods) ])
-end
-if (nargin==3 & size(index,2)~=3),
-	error('GetCharacteristicLength error message: index should have 3 columns for 2d meshes.')
-end
-if (nargin==4 & size(index,2)~=6),
-	error('GetCharacteristicLength error message: index should have 6 columns for 3d meshes.')
-end
-
-%get areas or volumes.
-areas=GetAreas(index,x,y,varargin{:});
-
-%for a 2d mesh: 
-if nargin==3,
-	length=sqrt(2*areas);
-else
-	error('not supported yet');
-end
Index: issm/trunk-jpl/src/m/mesh/NodeInElement.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/NodeInElement.m	(revision 13011)
+++ issm/trunk-jpl/src/m/mesh/NodeInElement.m	(revision 13012)
@@ -1,6 +1,7 @@
 function node_in_element=NodeInElement(newx,newy,elements,x,y,nodeconnectivity);
-%NODEINELEMENT: find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to.
+% NODEINELEMENT - find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to.
 %
-%  Usage: node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity);
+%    Usage:
+%      node_in_element=NodeInElement(newx,newy,elements,x,y,md.mesh.vertexconnectivity);
 %
 %  See also Nodeconnectivity
Index: issm/trunk-jpl/src/m/mesh/YamsCall.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/YamsCall.m	(revision 13011)
+++ 	(revision )
@@ -1,104 +1,0 @@
-function md=YamsCall(md,field,hmin,hmax,gradation,epsilon),
-%YAMSCALL - call yams
-%
-%   build a metric using the Hessian of the given field
-%   call Yams and the output mesh is plugged onto the model
-%   -hmin = minimum edge length (m)
-%   -hmax = maximum edge length (m)
-%   -gradation = maximum edge length gradation between 2 elements
-%   -epsilon = average error on each element (m/yr)
-%
-%   Usage:
-%      md=YamsCall(md,field,hmin,hmax,gradation,epsilon);
-%
-%   Example:
-%      md=YamsCall(md,md.inversion.vel_obs,1500,10^8,1.3,0.9);
-
-%2d geometric parameter (do not change)
-scale=2/9; 
-
-%Compute Hessian
-t1=clock; fprintf('%s','      computing Hessian...');
-hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,field,'node');
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%Compute metric
-t1=clock; fprintf('%s','      computing metric...');
-if length(md.mask.vertexonwater)==md.mesh.numberofvertices,
-	pos=find(md.mask.vertexonwater);
-else
-	pos=[];
-end
-metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%write files
-t1=clock; fprintf('%s','      writing initial mesh files...');
-save -ascii carre0.met  metric
-
-fid=fopen('carre0.mesh','w');
-
-%initialiation
-fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1);
-
-%dimension
-fprintf(fid,'\n%s\n%i\n','Dimension',2);
-
-%Vertices
-fprintf(fid,'\n%s\n%i\n\n','Vertices',md.mesh.numberofvertices);
-fprintf(fid,'%8g %8g %i\n',[md.mesh.x md.mesh.y zeros(md.mesh.numberofvertices,1)]');
-
-%Triangles
-fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.mesh.numberofelements);
-fprintf(fid,'%i %i %i %i\n',[md.mesh.elements zeros(md.mesh.numberofelements,1)]');
-numberofelements1=md.mesh.numberofelements;
-	
-%Deal with rifts
-if ~isnan(md.rifts.riftstruct),
-	
-	%we have the list of triangles that make up the rift. keep those triangles around during refinement.
-	triangles=[];
-	for i=1:size(md.rifts.riftstruct,1),
-		triangles=[triangles md.rifts(i).segments(:,3)'];
-	end
-
-	fprintf(fid,'\n\n%s\n%i\n\n','RequiredTriangles',length(triangles));
-	fprintf(fid,'%i\n',triangles);
-end
-
-%close
-fclose(fid);
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%call yams
-fprintf('%s\n','      call Yams...');
-if ispc
-	%windows
-	system(['yams2-win -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
-elseif ismac
-	%Macosx
-	system(['yams2-osx -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
-else
-	%Linux
-	system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
-end
-
-%plug new mesh
-t1=clock; fprintf('\n%s','      reading final mesh files...');
-Tria=load('carre1.tria');
-Coor=load('carre1.coor');
-md.mesh.x=Coor(:,1);
-md.mesh.y=Coor(:,2);
-md.mesh.z=zeros(size(Coor,1),1);
-md.mesh.elements=Tria;
-md.mesh.numberofvertices=size(Coor,1);
-md.mesh.numberofelements=size(Tria,1);
-numberofelements2=md.mesh.numberofelements;
-t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-%display number of elements
-fprintf('\n%s %i','      inital number of elements:',numberofelements1);
-fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
-
-%clean up:
-system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb');
Index: issm/trunk-jpl/src/m/mesh/isconnected.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/isconnected.m	(revision 13011)
+++ 	(revision )
@@ -1,13 +1,0 @@
-function flag=isconnected(elements,A,B)
-%ISCONNECTED: are two nodes connected by a triangulation?
-%
-%   Usage: flag=isconnected(elements,A,B)
-%
-%
-
-elements=ElementsFromEdge(elements,A,B);
-if isempty(elements),
-	flag=0;
-else
-	flag=1;
-end
Index: issm/trunk-jpl/src/m/mesh/meshread.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/meshread.m	(revision 13011)
+++ 	(revision )
@@ -1,41 +1,0 @@
-function Struct=meshread(filename);
-
-%some checks
-if ~exist(filename),
-	error(['meshread error message: file ' filename ' not found!']);
-end
-
-fid=fopen(filename,'r');
-
-while (~feof(fid)),
-
-	A=fscanf(fid,'%s',1);
-
-	if strcmp(A,'MeshVersionFormatted');
-		Struct.Version=fscanf(fid,'%s',1);
-
-	elseif strcmp(A,'Dimension'),
-		Struct.Dimension=fscanf(fid,'%i',1);
-
-	elseif strcmp(A,'Vertices'),
-		Struct.nods=fscanf(fid,'%i',1);
-		A=fscanf(fid,'%f %f %f',[3 Struct.nods]);
-		Struct.x=A(1,:)';
-		Struct.y=A(2,:)';
-
-	elseif strcmp(A,'Triangles'),
-		Struct.nels=fscanf(fid,'%i',1);
-		A=fscanf(fid,'%i %i %i',[4 Struct.nels]);
-		Struct.index=A(1:3,:)';
-
-	elseif strcmp(A,'Quadrilaterals'),
-		Struct.nels=fscanf(fid,'%i',1);
-		A=fscanf(fid,'%i %i %i %i',[5 Struct.nels]);
-		Struct.index=A(1:4,:)';
-	else
-		%do nothing
-
-	end
-end
-
-fclose(fid);
Index: issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m	(revision 13011)
+++ issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m	(revision 13012)
@@ -93,2 +93,18 @@
 md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
 md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
+end
+
+function flag=isconnected(elements,A,B)% {{{
+	%ISCONNECTED: are two nodes connected by a triangulation?
+	%
+	%   Usage: flag=isconnected(elements,A,B)
+	%
+	%
+
+	elements=ElementsFromEdge(elements,A,B);
+	if isempty(elements),
+		flag=0;
+	else
+		flag=1;
+	end
+end % }}}
Index: issm/trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m	(revision 13012)
+++ issm/trunk-jpl/src/m/mesh/rifts/rifttipsonmesh.m	(revision 13012)
@@ -0,0 +1,26 @@
+function tips=rifttipsonmesh(md,riftoutline)
+%RIFTTIPSONMESH: identify, using a rift outline, the nodes that are tips of 
+%                rifts.
+
+%read rifts from outline file
+rifts=expread(riftoutline,1);
+
+tips=[];
+
+for i=1:length(rifts),
+	rift=rifts(i);
+
+	x_tip=rift.x(1);
+	y_tip=rift.y(1);
+	
+	index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
+	tips(end+1)=index;
+
+	x_tip=rift.x(end);
+	y_tip=rift.y(end);
+	
+	index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
+	tips(end+1)=index;
+
+end
+
Index: issm/trunk-jpl/src/m/mesh/rifttipsonmesh.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/rifttipsonmesh.m	(revision 13011)
+++ 	(revision )
@@ -1,26 +1,0 @@
-function tips=rifttipsonmesh(md,riftoutline)
-%RIFTTIPSONMESH: identify, using a rift outline, the nodes that are tips of 
-%                rifts.
-
-%read rifts from outline file
-rifts=expread(riftoutline,1);
-
-tips=[];
-
-for i=1:length(rifts),
-	rift=rifts(i);
-
-	x_tip=rift.x(1);
-	y_tip=rift.y(1);
-	
-	index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
-	tips(end+1)=index;
-
-	x_tip=rift.x(end);
-	y_tip=rift.y(end);
-	
-	index=find_point(md.mesh.x,md.mesh.y,x_tip,y_tip);
-	tips(end+1)=index;
-
-end
-
Index: issm/trunk-jpl/src/m/mesh/yams.m
===================================================================
--- issm/trunk-jpl/src/m/mesh/yams.m	(revision 13011)
+++ 	(revision )
@@ -1,171 +1,0 @@
-function md=yams(md,varargin);
-%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
-%
-%   Usage:
-%      md=yams(md,varargin);
-%      where varargin is a lit of paired arguments. 
-%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
-%      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
-%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
-%                          this option is used to minimize the metric on water (no refinement)
-%      optional arguments: 'resolution': initial mesh resolution [m]
-%      optional arguments: 'nsteps': number of steps of mesh adaptation
-%      optional arguments: 'epsilon': average interpolation error wished [m/yr]
-%      optional arguments: 'hmin': minimum edge length
-%      optional arguments: 'hmanx': maximum edge
-%      optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
-%      
-%
-%   Examples:
-%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat');
-%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
-%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
-
-%recover options
-options=pairoptions(varargin{:});
-options=deleteduplicates(options,1);
-
-%recover some fields
-disp('MeshYams Options:')
-domainoutline=getfieldvalue(options,'domainoutline');
-disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
-riftoutline=getfieldvalue(options,'riftoutline','N/A');
-disp(sprintf('   %-15s: ''%s''','riftoutline',riftoutline));
-groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
-disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
-velocities=getfieldvalue(options,'velocities');
-disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
-resolution=getfieldvalue(options,'resolution',5000);
-disp(sprintf('   %-15s: %f','Resolution',resolution));
-nsteps=getfieldvalue(options,'nsteps',6);
-disp(sprintf('   %-15s: %i','nsteps',nsteps));
-gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
-disp(sprintf('   %-15s: %g','gradation',gradation(1)));
-epsilon=getfieldvalue(options,'epsilon',3);
-disp(sprintf('   %-15s: %f','epsilon',epsilon));
-hmin=getfieldvalue(options,'hmin',500);
-disp(sprintf('   %-15s: %f','hmin',hmin));
-hmax=getfieldvalue(options,'hmax',150*10^3);
-disp(sprintf('   %-15s: %f\n','hmax',hmax));
-
-%mesh with initial resolution
-disp('Initial mesh generation...');
-if strcmpi(riftoutline,'N/A');
-	md=setmesh(md,domainoutline,resolution);
-else
-	md=setmesh(md,domainoutline,riftoutline,resolution);
-	md=meshprocessrifts(md,domainoutline);
-end
-disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
-
-%load velocities 
-disp('loading velocities...');
-Names=VelFindVarNames(velocities);
-Vel=load(velocities);
-
-%start mesh adaptation
-for i=1:nsteps,
-	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
-
-	%interpolate velocities onto mesh
-	disp('   interpolating velocities...');
-	if strcmpi(Names.interp,'node'),
-		vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
-		vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
-	else
-		vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
-		vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
-	end
-	field=sqrt(vx_obs.^2+vy_obs.^2);
-
-	%set mask.vertexonwater  field
-	if ~strcmp(groundeddomain,'N/A'),
-		nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
-		md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
-		md.mask.vertexonwater(find(nodeground))=0;
-	else
-		md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
-	end
-
-	%adapt according to velocities
-	disp('   adapting...');
-	md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
-
-	%if we have rifts, we just messed them up, we need to recreate the segments that constitute those 
-	%rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
-	if md.rifts.numrifts, 
-		md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
-		md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
-		md.mesh.segments=findsegments(md);
-		md=yamsrecreateriftsegments(md);
-	end
-
-end
-	
-disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
-
-%Now, build the connectivity tables for this mesh.
-md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
-md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
-
-%recreate segments
-md.mesh.segments=findsegments(md);
-md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-
-%Fill in rest of fields:
-md.mesh.z=zeros(md.mesh.numberofvertices,1);
-md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
-md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
-md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
-if ~strcmp(groundeddomain,'N/A'),
-	nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
-	md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
-	md.mask.vertexonwater(find(nodeground))=0;
-else
-	md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
-end
-if strcmpi(Names.interp,'node'),
-	md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
-	md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
-else
-	md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
-	md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
-end
-md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
-
-%deal with rifts 
-if md.rifts.numrifts,
-	%first, recreate rift segments
-	md=meshyamsrecreateriftsegments(md);
-
-	%using the segments, recreate the penaltypairs
-	for j=1:md.rifts.numrifts,
-		rift=md.rifts.riftstruct(j);
-
-		%build normals and lengths of segments:
-		lengths=sqrt((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))).^2 + (md.mesh.y(rift.segments(:,1))-md.mesh.y(rift.segments(:,2))).^2 );
-		normalsx=cos(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.y(rift.segments(:,1)))));
-		normalsy=sin(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.y(rift.segments(:,1)))));
-
-		%ok, build penaltypairs: 
-		numpenaltypairs=length(rift.segments)/2-1;
-		rift.penaltypairs=zeros(numpenaltypairs,7);
-
-		for i=1:numpenaltypairs,
-			rift.penaltypairs(i,1)=rift.segments(i,2);
-			rift.penaltypairs(i,2)=rift.segments(end-i,2);
-			rift.penaltypairs(i,3)=rift.segments(i,3);
-			rift.penaltypairs(i,4)=rift.segments(end-i,3);
-			rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
-			rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
-			rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
-		end
-		%renormalize norms: 
-		norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
-		rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
-		rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
-
-		md.rifts.riftstruct(j)=rift;
-	end
-end
Index: issm/trunk-jpl/src/m/miscellaneous/create_region.m
===================================================================
--- issm/trunk-jpl/src/m/miscellaneous/create_region.m	(revision 13011)
+++ 	(revision )
@@ -1,15 +1,0 @@
-function create_region(name)
-%CREATE_REGION - create region  ????
-%
-%   very temporary function.
-%   
-%   Usage: 
-%      create_region(name)
-
-eval(['mkdir ' name]);
-eval(['cd ' name ]);
-!mkdir Delivery Exp_Par Results
-cd Exp_Par
-!cp ../../RonneShelf/Exp_Par/* ./
-!rm -rf Hole*
-eval(['!mv Ronne.par ' name '.par']);
Index: issm/trunk-jpl/src/m/miscellaneous/findarg.m
===================================================================
--- issm/trunk-jpl/src/m/miscellaneous/findarg.m	(revision 13011)
+++ 	(revision )
@@ -1,51 +1,0 @@
-function  vals=findarg(arglist,field)
-%FINDARG - find argument associated to a field in a list
-%
-%   This function parses through an argument list (typically varargin in a routine)
-%   looking for a character array equal to field. Once this is found, we return the 
-%   next value in the varargin (if possible). 
-%   Because field might appear several times in the argument list, we return a structure 
-%   holding all these values. 
-%   Note that all comparisons to field value are case independent.
-%
-%   Usage:
-%      vals=findarg(arglist,field)
-%
-%   Example:
-%      routine foobar calls vals=findarg('Data',varargin)
-%      with varargin='Data',1,'Data','foo','Plot','velocity','Arrow',4
-%      findarg would return the following structure: vals(1).value=1, vals(2).value='foo'; 
-
-%some argument checking: 
-if ((nargin==0) | (nargout==0)),
-	help findarg;
-	error('findarg error message');
-end
-
-if ~ischar(field),
-	error('findarg error message: field should be a string');
-end
-
-if ~iscell(arglist),
-	error('findarg error message: argument list should be a cell array.');
-end
-
-%Recover data to plot
-founddata=0;
-
-for i=1:(length(arglist)-1), %data in arglist comes in pairs, hence the -1.
-	if ischar(arglist{i}),
-		if (strcmpi(arglist{i},field)),
-			founddata=founddata+1;
-			if founddata==1,
-				vals.value=arglist{i+1};
-			else
-				vals(end+1).value=arglist{i+1};
-			end
-		end
-	end
-end
-
-if founddata==0,
-	vals=[];
-end
Index: issm/trunk-jpl/src/m/miscellaneous/netcdf.m
===================================================================
--- issm/trunk-jpl/src/m/miscellaneous/netcdf.m	(revision 13011)
+++ 	(revision )
@@ -1,136 +1,0 @@
-function S = netcdf(File,varargin)
-% Function to read NetCDF files
-%   S = netcdf(File)
-% Input Arguments
-%   File = NetCDF file to read
-% Optional Input Arguments:
-%   'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)]
-%   'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs]
-% Output Arguments:
-%   S    = Structure of NetCDF data organised as per NetCDF definition
-% Notes:
-%   Only version 1, classic 32bit, NetCDF files are supported. By default
-% data are extracted into the S.VarArray().Data field for all variables.
-% To read the header only call S = netcdf(File,'Var',[]);
-%
-% SEE ALSO
-% ---------------------------------------------------------------------------
-S = [];
-
-try
-   if exist(File,'file') fp = fopen(File,'r','b');
-   else fp = []; error('File not found'); end
-   if fp == -1   error('Unable to open file'); end
-
-% Read header
-   Magic = fread(fp,4,'uint8=>char');
-   if strcmp(Magic(1:3),'CDF') error('Not a NetCDF file'); end
-   if uint8(Magic(4))~=1       error('Version not supported'); end
-   S.NumRecs  = fread(fp,1,'uint32=>uint32');
-   S.DimArray = DimArray(fp);
-   S.AttArray = AttArray(fp);
-   S.VarArray = VarArray(fp);
-
-% Setup indexing to arrays and records
-   Var = ones(1,length(S.VarArray));
-   Rec = ones(1,S.NumRecs);
-   for i = 1:2:length(varargin)
-      if     strcmp(upper(varargin{i}),'VAR') Var=Var*0; Var(varargin{i+1})=1;
-      elseif strcmp(upper(varargin{i}),'REC') Rec=Rec*0; Rec(varargin{i+1})=1;
-      else error('Optional input argument not recognised'); end
-   end
-   if sum(Var)==0 fclose(fp); return; end
-
-% Read non-record variables
-   Dim = double(cat(2,S.DimArray.Dim));
-   ID  = double(cat(2,S.VarArray.Type));
-
-   for i = 1:length(S.VarArray)
-      D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0);
-      if isempty(RecID{i})
-         if length(D)==0 D = [1,1]; N = 1; elseif length(D)==1 D=[D,1]; end
-         if Var(i)
-            S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D);
-            fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8');
-         else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
-      else S.VarArray(i).Data = []; end
-   end
-
-% Read record variables
-   for k = 1:S.NumRecs
-      for i = 1:length(S.VarArray)
-         if ~isempty(RecID{i})
-            D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D);
-            if length(D)==1 D=[D,1]; end
-            if Var(i) & Rec(k)
-               S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,...
-                  ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D));
-               if N > 1 fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); end
-            else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
-         end
-      end
-   end
-
-   fclose(fp);
-catch
-   Err = lasterror; fprintf('%s\n',Err.message);
-   if ~isempty(fp) && fp ~= -1 fclose(fp); end
-end
-
-% ---------------------------------------------------------------------------------------
-% Utility functions
-
-function S = Size(ID)
-% Size of NetCDF data type, ID, in bytes
-   S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}}));
-
-function T = Type(ID)
-% Matlab string for CDF data type, ID
-   T = subsref({'int8','char','int16','int32','single','double'},...
-               struct('type','{}','subs',{{ID}}));
-
-function N = Pad(Num,ID)
-% Number of elements to read after padding to 4 bytes for type ID
-   N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0);
-
-function S = String(fp)
-% Read a CDF string; Size,[String,[Padding]]
-   S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').';
-
-function A = ReOrder(A,S)
-% Rearrange CDF array A to size S with matlab ordering
-   A = permute(reshape(A,fliplr(S)),fliplr(1:length(S)));
-
-function S = DimArray(fp)
-% Read DimArray into structure
-   if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION
-      for i = 1:fread(fp,1,'uint32=>uint32')
-         S(i).Str = String(fp);
-         S(i).Dim = fread(fp,1,'uint32=>uint32');
-      end
-   else fread(fp,1,'uint32=>uint32'); S = []; end
-
-function S = AttArray(fp)
-% Read AttArray into structure
-   if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE
-      for i = 1:fread(fp,1,'uint32=>uint32')
-         S(i).Str = String(fp);
-         ID       = fread(fp,1,'uint32=>uint32');
-         Num      = fread(fp,1,'uint32=>uint32');
-         S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).';
-      end
-   else fread(fp,1,'uint32=>uint32'); S = []; end
-
-function S = VarArray(fp)
-% Read VarArray into structure
-   if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE
-      for i = 1:fread(fp,1,'uint32=>uint32')
-         S(i).Str      = String(fp);
-         Num           = double(fread(fp,1,'uint32=>uint32'));
-         S(i).DimID    = double(fread(fp,Num,'uint32=>uint32'));
-         S(i).AttArray = AttArray(fp);
-         S(i).Type     = fread(fp,1,'uint32=>uint32');
-         S(i).VSize    = fread(fp,1,'uint32=>uint32');
-         S(i).Begin    = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only
-      end
-   else fread(fp,1,'uint32=>uint32'); S = []; end
Index: issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m
===================================================================
--- issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m	(revision 13011)
+++ issm/trunk-jpl/src/m/miscellaneous/netcdf2struct.m	(revision 13012)
@@ -6,5 +6,5 @@
 
 %Read netcdf file
-data=netcdf(File);
+data=readnetcdf(File);
 
 %initialize output
@@ -26,2 +26,141 @@
 	S.(fieldname)=fieldvalue;
 end
+end
+
+function S = readnetcdf(File,varargin)
+% Function to read NetCDF files
+%   S = netcdf(File)
+% Input Arguments
+%   File = NetCDF file to read
+% Optional Input Arguments:
+%   'Var',Var - Read data for VarArray(Var), default [1:length(S.VarArray)]
+%   'Rec',Rec - Read data for Record(Rec), default [1:S.NumRecs]
+% Output Arguments:
+%   S    = Structure of NetCDF data organised as per NetCDF definition
+% Notes:
+%   Only version 1, classic 32bit, NetCDF files are supported. By default
+% data are extracted into the S.VarArray().Data field for all variables.
+% To read the header only call S = netcdf(File,'Var',[]);
+%
+% SEE ALSO
+% ---------------------------------------------------------------------------
+S = [];
+
+try
+   if exist(File,'file') fp = fopen(File,'r','b');
+   else fp = []; error('File not found'); end
+   if fp == -1   error('Unable to open file'); end
+
+% Read header
+   Magic = fread(fp,4,'uint8=>char');
+   if strcmp(Magic(1:3),'CDF') error('Not a NetCDF file'); end
+   if uint8(Magic(4))~=1       error('Version not supported'); end
+   S.NumRecs  = fread(fp,1,'uint32=>uint32');
+   S.DimArray = DimArray(fp);
+   S.AttArray = AttArray(fp);
+   S.VarArray = VarArray(fp);
+
+% Setup indexing to arrays and records
+   Var = ones(1,length(S.VarArray));
+   Rec = ones(1,S.NumRecs);
+   for i = 1:2:length(varargin)
+      if     strcmp(upper(varargin{i}),'VAR') Var=Var*0; Var(varargin{i+1})=1;
+      elseif strcmp(upper(varargin{i}),'REC') Rec=Rec*0; Rec(varargin{i+1})=1;
+      else error('Optional input argument not recognised'); end
+   end
+   if sum(Var)==0 fclose(fp); return; end
+
+% Read non-record variables
+   Dim = double(cat(2,S.DimArray.Dim));
+   ID  = double(cat(2,S.VarArray.Type));
+
+   for i = 1:length(S.VarArray)
+      D = Dim(S.VarArray(i).DimID+1); N = prod(D); RecID{i}=find(D==0);
+      if isempty(RecID{i})
+         if length(D)==0 D = [1,1]; N = 1; elseif length(D)==1 D=[D,1]; end
+         if Var(i)
+            S.VarArray(i).Data = ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D);
+            fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8');
+         else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
+      else S.VarArray(i).Data = []; end
+   end
+
+% Read record variables
+   for k = 1:S.NumRecs
+      for i = 1:length(S.VarArray)
+         if ~isempty(RecID{i})
+            D = Dim(S.VarArray(i).DimID+1); D(RecID{i}) = 1; N = prod(D);
+            if length(D)==1 D=[D,1]; end
+            if Var(i) & Rec(k)
+               S.VarArray(i).Data = cat(RecID{i},S.VarArray(i).Data,...
+                  ReOrder(fread(fp,N,[Type(ID(i)),'=>',Type(ID(i))]),D));
+               if N > 1 fread(fp,(Pad(N,ID(i))-N)*Size(ID(i)),'uint8=>uint8'); end
+            else fseek(fp,Pad(N,ID(i))*Size(ID(i)),'cof'); end
+         end
+      end
+   end
+
+   fclose(fp);
+catch
+   Err = lasterror; fprintf('%s\n',Err.message);
+   if ~isempty(fp) && fp ~= -1 fclose(fp); end
+end
+
+% ---------------------------------------------------------------------------------------
+% Utility functions
+
+function S = Size(ID)
+% Size of NetCDF data type, ID, in bytes
+   S = subsref([1,1,2,4,4,8],struct('type','()','subs',{{ID}}));
+
+function T = Type(ID)
+% Matlab string for CDF data type, ID
+   T = subsref({'int8','char','int16','int32','single','double'},...
+               struct('type','{}','subs',{{ID}}));
+
+function N = Pad(Num,ID)
+% Number of elements to read after padding to 4 bytes for type ID
+   N = (double(Num) + mod(4-double(Num)*Size(ID),4)/Size(ID)).*(Num~=0);
+
+function S = String(fp)
+% Read a CDF string; Size,[String,[Padding]]
+   S = fread(fp,Pad(fread(fp,1,'uint32=>uint32'),1),'uint8=>char').';
+
+function A = ReOrder(A,S)
+% Rearrange CDF array A to size S with matlab ordering
+   A = permute(reshape(A,fliplr(S)),fliplr(1:length(S)));
+
+function S = DimArray(fp)
+% Read DimArray into structure
+   if fread(fp,1,'uint32=>uint32') == 10 % NC_DIMENSION
+      for i = 1:fread(fp,1,'uint32=>uint32')
+         S(i).Str = String(fp);
+         S(i).Dim = fread(fp,1,'uint32=>uint32');
+      end
+   else fread(fp,1,'uint32=>uint32'); S = []; end
+
+function S = AttArray(fp)
+% Read AttArray into structure
+   if fread(fp,1,'uint32=>uint32') == 12 % NC_ATTRIBUTE
+      for i = 1:fread(fp,1,'uint32=>uint32')
+         S(i).Str = String(fp);
+         ID       = fread(fp,1,'uint32=>uint32');
+         Num      = fread(fp,1,'uint32=>uint32');
+         S(i).Val = fread(fp,Pad(Num,ID),[Type(ID),'=>',Type(ID)]).';
+      end
+   else fread(fp,1,'uint32=>uint32'); S = []; end
+
+function S = VarArray(fp)
+% Read VarArray into structure
+   if fread(fp,1,'uint32=>uint32') == 11 % NC_VARIABLE
+      for i = 1:fread(fp,1,'uint32=>uint32')
+         S(i).Str      = String(fp);
+         Num           = double(fread(fp,1,'uint32=>uint32'));
+         S(i).DimID    = double(fread(fp,Num,'uint32=>uint32'));
+         S(i).AttArray = AttArray(fp);
+         S(i).Type     = fread(fp,1,'uint32=>uint32');
+         S(i).VSize    = fread(fp,1,'uint32=>uint32');
+         S(i).Begin    = fread(fp,1,'uint32=>uint32'); % Classic 32 bit format only
+      end
+   else fread(fp,1,'uint32=>uint32'); S = []; end
+end
