Index: /issm/trunk-jpl/test/NightlyRun/test4003.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test4003.m	(revision 24248)
+++ /issm/trunk-jpl/test/NightlyRun/test4003.m	(revision 24249)
@@ -16,5 +16,5 @@
 presentdirectory=pwd;
 
-% {{{ Parameters: 
+% {{{ Parameters:
 if perform(org,'Parameters'),
 	Nx=3; %number of longitude cells
@@ -28,5 +28,5 @@
 	dLat=1; %latitude grid spacing
 	delZ=10; %thickness of vertical levels
-	icefront_position_ratio=.75; 
+	icefront_position_ratio=.75;
 	ice_thickness=1000;
 	rho_ice=917;
@@ -42,5 +42,5 @@
 	mlSalt   = 33.4; % open boundary salinity (PSU)
 	mlTheta  = -1.9; % open boundary potential temperature (deg C)
-	obcUvel  = -0.1; % open boundary velocity (m/s)        
+	obcUvel  = -0.1; % open boundary velocity (m/s)
 
 	MITgcmDeltaT=600; % MITgcm time step in seconds
@@ -69,5 +69,5 @@
 end
 % }}}
-% {{{ Bathymetry: 
+% {{{ Bathymetry:
 if perform(org,'Bathymetry'),
 
@@ -94,7 +94,7 @@
 end
 % }}}
-% {{{ IceSheetGeometry: 
+% {{{ IceSheetGeometry:
 if perform(org,'IceSheetGeometry'),
-	
+
 	loaddata(org,'Parameters');
 	loaddata(org,'Bathymetry');
@@ -102,9 +102,9 @@
 	latmax=max(lat(:));
 
-	%put ice_thickness constant layer of ice over the bathymetry, unless it floats: 
+	%put ice_thickness constant layer of ice over the bathymetry, unless it floats:
 	s=size(bathymetry);
 	thickness=ice_thickness*ones(s);
 
-	%figure out ice shelf: 
+	%figure out ice shelf:
 	pos=find(-di*thickness>bathymetry);
 	iceshelf_mask=zeros(s);
@@ -116,9 +116,9 @@
 	iceshelf_mask(pos)=0;
 
-	%compute draft of ice shelf: 
+	%compute draft of ice shelf:
 	draft=bathymetry;
 	pos=find(iceshelf_mask);
 	draft(pos)=-di*thickness(pos);
-	pos=find(~ice_mask); 
+	pos=find(~ice_mask);
 	draft(pos)=0;
 
@@ -128,10 +128,10 @@
 
 %Configure MITgcm
-% {{{ GetMITgcm: 
+% {{{ GetMITgcm:
 if perform(org,'GetMITgcm'),
   system([pwd '/../MITgcm/get_mitgcm_dngoldberg.sh']);
 end
 % }}}
-% {{{ BuildMITgcm: 
+% {{{ BuildMITgcm:
 if perform(org,'BuildMITgcm'),
     system(['../MITgcm/build_4003.sh generic ' pwd '/../MITgcm']);
@@ -139,5 +139,5 @@
 % }}}
 addpath(recursivepath([pwd '/../MITgcm']));
-% {{{ RunUncoupledMITgcm: 
+% {{{ RunUncoupledMITgcm:
 if perform(org,'RunUncoupledMITgcm'),
     % rename previous run directory and create new one
@@ -161,11 +161,11 @@
 
 %Configure ISSM
-% {{{ CreateMesh: 
+% {{{ CreateMesh:
 if perform(org,'CreateMesh'),
-	
-	%create model: 
+
+	%create model:
 	md=model();
-	
-	%Grab lat,long from MITgcm: 
+
+	%Grab lat,long from MITgcm:
 	long=readbin('run/XG.data',[3 200]);
 	long=[long long(:,end)]; long=[long; -105.1250*ones(1,size(long,2))];
@@ -173,5 +173,5 @@
 	lat=[lat -73.8832*ones(size(lat,1),1)]; lat=[lat; lat(end,:)];
 
-	%project lat,long: 
+	%project lat,long:
 	[x,y]=ll2xy(lat(:),long(:),-1);
 
@@ -179,7 +179,7 @@
 	index=[];
 	%  C  D
-	%  A  B 
+	%  A  B
 	for j=1:Ny-1,
-		for i=1:Nx-1, 
+		for i=1:Nx-1,
 			A=(j-1)*Nx+i;
 			B=(j-1)*Nx+i+1;
@@ -191,5 +191,5 @@
 	end
 
-	%fill mesh and model: 
+	%fill mesh and model:
 	md=meshconvert(md,index,x,y);
 	md.mesh.lat=lat(:);
@@ -200,10 +200,10 @@
 end
 % }}}
-% {{{ MeshGeometry: 
+% {{{ MeshGeometry:
 if perform(org,'MeshGeometry'),
-	
+
 	loaddata(org,'CreateMesh');
 
-	%transfer to vertices: 
+	%transfer to vertices:
 	bathymetry=readbin('run/bathy.box',[3 200],1,'real*8');
 	bathymetry=[bathymetry bathymetry(:,end)]; bathymetry=[bathymetry(1,:); bathymetry];
@@ -214,11 +214,11 @@
 	thickness=[thickness thickness(:,end)]; thickness=[thickness; thickness(end,:)];
 
-	%start filling some of the fields 
+	%start filling some of the fields
 	md.geometry.bed=bathymetry(:);
 	md.geometry.thickness=thickness(:);
-	md.geometry.base=-917/1028*md.geometry.thickness; 
+	md.geometry.base=-917/1028*md.geometry.thickness;
 	md.geometry.surface=md.geometry.base+md.geometry.thickness;
 
-	%nothing passes icefront: 
+	%nothing passes icefront:
 	pos=find((~ice_mask(:) & ice_mask(:)~=0) | thickness(:)==0);
 	md.geometry.thickness(pos)=1;
@@ -226,16 +226,16 @@
 	md.geometry.base(pos)=-di*md.geometry.thickness(pos);
 
-	%level sets: 
+	%level sets:
 	md.mask.groundedice_levelset=iceshelf_mask(:);
 	md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
 
-	pos=find(~ice_mask(:) & thickness(:)==0); md.mask.ice_levelset(pos)=1; 
+	pos=find(~ice_mask(:) & thickness(:)==0); md.mask.ice_levelset(pos)=1;
 
 	savemodel(org,md);
 end
 % }}}
-% {{{ ParameterizeIce: 
+% {{{ ParameterizeIce:
 if perform(org,'ParameterizeIce'),
-	
+
 	loaddata(org,'MeshGeometry');
 
@@ -243,10 +243,10 @@
 	md.miscellaneous.name='test4003';
 
-	%initial velocity: 
+	%initial velocity:
 	md.initialization.vx=zeros(md.mesh.numberofvertices,1);
 	md.initialization.vy=zeros(md.mesh.numberofvertices,1);
 	md.initialization.vz=zeros(md.mesh.numberofvertices,1);
 
-	%friction: 
+	%friction:
 	md.friction.coefficient=0*ones(md.mesh.numberofvertices,1);
 	pos=find(md.mask.groundedice_levelset>0);
@@ -260,10 +260,10 @@
 	md.smb.mass_balance = 0*ones(md.mesh.numberofvertices,1);
 
-	%Flow law 
+	%Flow law
 	md.materials.rheology_B=paterson(md.initialization.temperature);
 	md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
 	md.damage.D=zeros(md.mesh.numberofvertices,1);
 	md.damage.spcdamage=NaN*ones(md.mesh.numberofvertices,1);
-	
+
 	%the spcs going
 	md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
@@ -272,7 +272,7 @@
 	md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
 	md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
-	md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1); 
-
-	%get some flux at the ice divide: 
+	md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
+
+	%get some flux at the ice divide:
 	pos=find(md.mesh.lat==min(md.mesh.lat));
 	md.masstransport.spcthickness(pos)=md.geometry.thickness(pos);
@@ -280,5 +280,5 @@
 	md.stressbalance.spcvy(pos)=0;
 
-	%deal with boundaries, excluding icefront: 
+	%deal with boundaries, excluding icefront:
 	pos=find(md.mesh.long==min(md.mesh.long) | md.mesh.long==max(md.mesh.long));
 	md.stressbalance.spcvy(pos)=0;
@@ -295,5 +295,5 @@
 	md.basalforcings.geothermalflux=.064*ones(md.mesh.numberofvertices,1);
 
-	%flow equations: 
+	%flow equations:
 	md=setflowequation(md,'SSA','all');
 
@@ -301,11 +301,11 @@
 end
 % }}}
-% {{{ RunUncoupledISSM: 
+% {{{ RunUncoupledISSM:
 if perform(org,'RunUncoupledISSM'),
-	
+
 	loaddata(org,'Parameters');
 	loaddata(org,'ParameterizeIce');
 
-	%timestepping: 
+	%timestepping:
 	md.timestepping.final_time=100;
 	md.timestepping.time_step=0.5;
@@ -327,7 +327,7 @@
 
 %Run MITgcm/ISSM
-% {{{ RunCoupledMITgcmISSM: 
+% {{{ RunCoupledMITgcmISSM:
 if perform(org,'RunCoupledMITgcmISSM'),
-	
+
 	loaddata(org,'Parameters');
 	loaddata(org,'ParameterizeIce');
@@ -352,5 +352,5 @@
 
         % }}}
-        % {{{ prepare MITgcm 
+        % {{{ prepare MITgcm
         % rename previous run directory and create new one
         if exist ('run.old')
@@ -379,8 +379,8 @@
 end
 % }}}
-% {{{ RunCoupledMITgcmISSMMatlabCoupling: 
+% {{{ RunCoupledMITgcmISSMMatlabCoupling:
 if perform(org,'RunCoupledMITgcmISSMMatlabCoupling'),
 
-	%load data: 
+	%load data:
 	loaddata(org,'Parameters');
 	loaddata(org,'ParameterizeIce');
@@ -390,5 +390,5 @@
          floor(time_step*y2s*async_step_MITgcm_multiplier/MITgcmDeltaT));
 
-        % {{{ prepare MITgcm 
+        % {{{ prepare MITgcm
         % rename previous run directory and create new one
         if exist ('run.old')
@@ -449,5 +449,5 @@
             tline = fgetl(fidi);
             if ~ischar(tline), break, end
-            %do the change here: 
+            %do the change here:
             if strcmpi(tline,' OB_Iwest = 40*1,'),
                 fprintf(fido,'%s%i%s\n',' OB_Iwest = ',Ny,'*1,');
@@ -467,9 +467,9 @@
         writebin('run/bathymetry.bin',bathymetry);
         % }}}
-            
+
 	% {{{ ISSM settings:
 
-	setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib') 
-	%timestepping: 
+	setenv('DYLD_LIBRARY_PATH', '/usr/local/gfortran/lib')
+	%timestepping:
 	md.timestepping.start_time=start_time;
 	md.timestepping.final_time=final_time;
@@ -499,5 +499,5 @@
             writebin('run/icetopo.bin',draft);
 
-	    % {{{ generate MITgcm parameter file data 
+	    % {{{ generate MITgcm parameter file data
 	    fidi=fopen('../MITgcm/input/data','r');
 	    fido=fopen('run/data','w');
@@ -507,5 +507,5 @@
                 tline = fgetl(fidi);
                 if ~ischar(tline), break, end
-		%do the change here: 
+		%do the change here:
 		if strcmpi(tline,' xgOrigin = 0.0,'),
 		    fprintf(fido,'%s%i%s\n',' xgOrigin = ',xgOrigin,',');
@@ -578,5 +578,5 @@
                 tmp(find(tmp<0))=0;
                 [im jm]=find(tmp); % horizontal indices where there is melt
-                
+
                 % Extrapolate T/S to locations where ice shelf retreated
                 for i=1:length(im)
@@ -612,5 +612,5 @@
             % }}}
 
-            % {{{ system call to run MITgcm 
+            % {{{ system call to run MITgcm
             cd run
             eval(['!mpiexec -np ' int2str(nPx*nPy) ' ./mitgcmuv']);
@@ -638,5 +638,5 @@
 		 md.basalforcings.floatingice_melting_rate=melting_rate;
 
-	    % {{{ run ISSM and recover results 
+	    % {{{ run ISSM and recover results
 
 	    md.timestepping.start_time=t;
@@ -655,5 +655,5 @@
 		md.mask.groundedice_levelset=md.results.TransientSolution(end).MaskGroundediceLevelset;
 		md.results.TransientSolution(end).FloatingiceMeltingRate=md.basalforcings.floatingice_melting_rate;
-		
+
 		%save these results in the model, otherwise, they'll be wiped out
 		results(end+1)=md.results;
@@ -681,5 +681,5 @@
 	'Melting2','Vx3','Vy3','Thickness3','Base3','MaskGroundediceLevelset3','FloatingiceMeltingRate3',...
 	'Melting3','Vx4','Vy4','Thickness4','Base4','MaskGroundediceLevelset4','FloatingiceMeltingRate4','Melting4'};
-field_tolerances={2e-13,1e-13,5e-13,5e-13,1e-13,1e-13,1e-13,1e-12,...
+field_tolerances={2e-13,1e-13,6e-13,6e-13,1e-13,1e-13,1e-13,1e-12,...
 	1e-13, 1e-12, 1e-12, 1e-13, 1e-13, 1e-13, 1e-12,...
 	1e-13, 1e-12, 1e-12, 1e-13, 1e-13, 1e-13, 1e-12, 1e-13 };