Index: /issm/trunk-jpl/examples/Pig/CheatSheet.m
===================================================================
--- /issm/trunk-jpl/examples/Pig/CheatSheet.m	(revision 20732)
+++ /issm/trunk-jpl/examples/Pig/CheatSheet.m	(revision 20733)
@@ -1,4 +1,3 @@
-if any(steps==6)
-	% Load Model
+if step==6
 	md = loadmodel('./Models/PIG.Control_drag');
 	md.inversion.iscontrol=0;
@@ -6,13 +5,11 @@
 	disp('   Extruding mesh')
 	number_of_layers=3;
-	md=extrude(md,number_of_layers,0.9);
+	md=extrude(md,number_of_layers,1);
 
 	disp('   Using HO Ice Flow Model')
 	md=setflowequation(md, 'HO', 'all');
 
-	% Solve
 	md=solve(md,StressbalanceSolutionEnum);
 
-	% Save Model
 	save ./Models/PIG.ModelHO md;
 end
Index: /issm/trunk-jpl/examples/Pig/Pig.par
===================================================================
--- /issm/trunk-jpl/examples/Pig/Pig.par	(revision 20732)
+++ /issm/trunk-jpl/examples/Pig/Pig.par	(revision 20733)
@@ -17,34 +17,4 @@
 smb   = ncread(ncdata,'presprcp')';
 gflux = ncread(ncdata,'bheatflx_fox')';
-
-disp('   Loading velocities data from NetCDF');
-nsidc_vel='../Data/Antarctica_ice_velocity.nc';
-
-xmin = ncreadatt(nsidc_vel,'/','xmin');
-xmin = strtrim(xmin);  %this is a string, and we need to recover the double value
-xmin = xmin(1:end-2);  %get rid of the unit
-xmin = str2num(xmin);  %convert to double
-
-ymax = ncreadatt(nsidc_vel,'/','ymax'); 
-ymax = strtrim(ymax);  
-ymax = ymax(1:end-2);  
-ymax = str2num(ymax); 
-
-nx = ncreadatt(nsidc_vel,'/','nx');
-ny = ncreadatt(nsidc_vel,'/','ny');
-
-spacing = ncreadatt(nsidc_vel,'/','spacing'); 
-spacing = strtrim(spacing);  
-spacing = spacing(1:end-2);  
-spacing = str2num(spacing); 
-
-velx = double(ncread(nsidc_vel,'vx'));
-vely = double(ncread(nsidc_vel,'vy'));
-
-x2=xmin+(0:1:nx)'*spacing; 
-x2=double(x2);
-
-y2=(ymax-ny*spacing)+(0:1:ny)'*spacing; 
-y2=double(y2);
 
 %Geometry
@@ -70,5 +40,5 @@
 
 %Set min thickness to 1 meter
-pos0=find(md.geometry.thickness<=0);
+pos0=find(md.geometry.thickness<=1);
 md.geometry.thickness(pos0)=1;
 md.geometry.surface=md.geometry.thickness+md.geometry.base;
@@ -79,14 +49,30 @@
 clear temp;
 
+disp('   Loading velocities data from NetCDF');
+nsidc_vel='../Data/Antarctica_ice_velocity.nc';
+xmin    = ncreadatt(nsidc_vel,'/','xmin');
+ymax    = ncreadatt(nsidc_vel,'/','ymax');
+spacing = ncreadatt(nsidc_vel,'/','spacing');
+nx      = double(ncreadatt(nsidc_vel,'/','nx'));
+ny      = double(ncreadatt(nsidc_vel,'/','ny'));
+velx    = double(ncread(nsidc_vel,'vx'));
+vely    = double(ncread(nsidc_vel,'vy'));
+% Read coordinates
+xmin = strtrim(xmin);  
+xmin = str2num(xmin(1:end-2)); 
+ymax = strtrim(ymax);  
+ymax = str2num(ymax(1:end-2));  
+spacing = strtrim(spacing);
+spacing = str2num(spacing(1:end-2));  
+% Build the coordinates
+x2=xmin+(0:1:nx)'*spacing;
+y2=(ymax-ny*spacing)+(0:1:ny)'*spacing;
+
 disp('   Set observed velocities')
-vx_obs=InterpFromGridToMesh(x2,y2,flipud(velx'),md.mesh.x,md.mesh.y,0);
-vy_obs=InterpFromGridToMesh(x2,y2,flipud(vely'),md.mesh.x,md.mesh.y,0);
+md.initialization.vx=InterpFromGridToMesh(x2,y2,flipud(velx'),md.mesh.x,md.mesh.y,0);
+md.initialization.vy=InterpFromGridToMesh(x2,y2,flipud(vely'),md.mesh.x,md.mesh.y,0);
+md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+md.initialization.vel=sqrt(md.initialization.vx.^2+md.initialization.vy.^2);
 clear velx vely;
-
-vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
-md.initialization.vx=vx_obs;
-md.initialization.vy=vy_obs;
-md.initialization.vz=zeros(md.mesh.numberofvertices,1);
-md.initialization.vel=vel_obs;
 
 disp('   Set Pressure');
@@ -99,6 +85,6 @@
 %Forcings
 disp('   Interpolating surface mass balance');
-md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
-md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
+mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
+md.smb.mass_balance=mass_balance*md.materials.rho_water/md.materials.rho_ice;
 clear smb;
 
@@ -118,7 +104,7 @@
 
 md.inversion=m1qn3inversion();
-md.inversion.vx_obs=vx_obs;
-md.inversion.vy_obs=vy_obs;
-md.inversion.vel_obs=vel_obs;
+md.inversion.vx_obs=md.initialization.vx;
+md.inversion.vy_obs=md.initialization.vy;
+md.inversion.vel_obs=md.initialization.vel;
 
 disp('   Set boundary conditions');
@@ -126,4 +112,4 @@
 md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
 md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
-md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
-md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);
+md.thermal.spctemperature                 = md.initialization.temperature;
+md.masstransport.spcthickness             = NaN*ones(md.mesh.numberofvertices,1);
Index: /issm/trunk-jpl/examples/Pig/runme.m
===================================================================
--- /issm/trunk-jpl/examples/Pig/runme.m	(revision 20732)
+++ /issm/trunk-jpl/examples/Pig/runme.m	(revision 20733)
@@ -1,22 +1,15 @@
-%Which steps to be performed
-steps=[1] ;
+step=[7];
 
-%Run Steps
-
-%Mesh Generation #1
-if any(steps==1)
-
-	md.miscellaneous.name='PIG.Mesh_generation';
-
+if step==1   %Mesh Generation #1
 	%Mesh parameters
 	domain =['./DomainOutline.exp'];
+	hinit=10000;   % element size for the initial mesh
 	hmax=40000;    % maximum element size of the final mesh
 	hmin=5000;     % minimum element size of the final mesh
-	hinit=10000;   % element size for the initial mesh
 	gradation=1.7; % maximum size ratio between two neighboring elements
 	err=8;         % maximum error between interpolated and control field
 
 	% Generate an initial uniform mesh (resolution = hinit m)
-	md=bamg(model,'domain',domain,'hmax',hinit,'MaxCornerAngle',1);
+	md=bamg(model,'domain',domain,'hmax',hinit);
 
 	%ploting
@@ -24,35 +17,27 @@
 
 	% Load Velocities
-	% http://nsidc.org/data/nsidc-0484.html
 	nsidc_vel='../Data/Antarctica_ice_velocity.nc'; 	
 
 	% Get necessary data to build up the velocity grid
-	xmin = ncreadatt(nsidc_vel,'/','xmin');
-	xmin = strtrim(xmin);  % this is a string, and we need to recover the double value
-	xmin = xmin(1:end-2);  % get rid of the unit
-	xmin = str2num(xmin);  % convert to double
+	xmin    = ncreadatt(nsidc_vel,'/','xmin');
+	ymax    = ncreadatt(nsidc_vel,'/','ymax');
+	spacing = ncreadatt(nsidc_vel,'/','spacing');
+	nx      = double(ncreadatt(nsidc_vel,'/','nx'));
+	ny      = double(ncreadatt(nsidc_vel,'/','ny'));
+	vx      = double(ncread(nsidc_vel,'vx'));
+	vy      = double(ncread(nsidc_vel,'vy'));
+
+	% Read coordinates
+	xmin = strtrim(xmin);  
+	xmin = str2num(xmin(1:end-2)); 
+	ymax = strtrim(ymax);  
+	ymax = str2num(ymax(1:end-2));  
+	spacing = strtrim(spacing);
+	spacing = str2num(spacing(1:end-2));  
+
+	% Build the coordinates
+	x=xmin+(0:1:nx)'*spacing;
+	y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
 	
-	ymax = ncreadatt(nsidc_vel,'/','ymax'); 
-	ymax = strtrim(ymax);  
-	ymax = ymax(1:end-2);  
-	ymax = str2num(ymax); 
-	
-	nx = ncreadatt(nsidc_vel,'/','nx');
-	ny = ncreadatt(nsidc_vel,'/','ny');
-	
-	spacing = ncreadatt(nsidc_vel,'/','spacing'); 
-	spacing = strtrim(spacing);
-	spacing = spacing(1:end-2);  
-	spacing = str2num(spacing); 
-	
-	% Get velocities (Note: You can use ncdisp('file') to see an ncdump)
-	vx = double(ncread(nsidc_vel,'vx'));
-	vy = double(ncread(nsidc_vel,'vy'));
-	
-	x=xmin+(0:1:nx)'*spacing; 
-	x=double(x);
-	y=(ymax-ny*spacing)+(0:1:ny)'*spacing; 
-	y=double(y);
-
 	% Interpolate velocities onto coarse mesh
 	vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
@@ -67,13 +52,9 @@
 	plotmodel(md,'data','mesh')
 
-	% Convert x,y coordinates (Polar stereo) to lat/lon
-	[md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,-1);
-	
 	% Save model
 	save ./Models/PIG.Mesh_generation md;
 end
 
-%Masks #2
-if any(steps==2) 
+if step==2  %Masks #2
 
 	md = loadmodel('./Models/PIG.Mesh_generation');	
@@ -104,6 +85,5 @@
 end
 
-%Parameterization #3
-if any(steps==3) 
+if step==3  %Parameterization #3
 
 	md = loadmodel('./Models/PIG.SetMask');
@@ -117,6 +97,5 @@
 end
 
-%Control Method #4
-if any(steps==4)
+if step==4  %Control Method #4
 
 	md = loadmodel('./Models/PIG.Parameterization');
@@ -161,6 +140,5 @@
 end
 
-%Plot #5
-if any(steps==5)
+if step==5 %Plot #5
 
 	md = loadmodel('./Models/PIG.Control_drag');
@@ -179,6 +157,5 @@
 end
 
-%HO #6
-if any(steps==6)
+if step==6  %Higher-Order #6
 
 	% Load Model
@@ -196,6 +173,5 @@
 end
 
-%Plot #7
-if any(steps==7)
+if step==7  %Plot #7
 
 	mdHO = loadmodel('./Models/PIG.ModelHO');