Index: /issm/trunk/test/Validation/EISMINT/Diagnostic1/Square.par
===================================================================
--- /issm/trunk/test/Validation/EISMINT/Diagnostic1/Square.par	(revision 899)
+++ /issm/trunk/test/Validation/EISMINT/Diagnostic1/Square.par	(revision 900)
@@ -1,123 +1,32 @@
-
 %Ok, start defining model parameters here
 
-%material parameters
-	md.g=9.81;
-	md.rho_ice=917;
-	md.rho_water=1028;
-	di=md.rho_ice/md.rho_water;
-	md.yts=365*24*3600;
-	md.heatcapacity=2009;
-	md.thermalconductivity=2.2; %W/mK
-	md.beta=9.8*10^-8;
+disp('      creating thickness');
+ymin=min(md.y);
+ymax=max(md.y);
+md.thickness=500*ones(md.numberofgrids,1);
+md.firn_layer=0*ones(md.numberofgrids,1);
+md.bed=-di*md.thickness;
+md.surface=md.bed+md.thickness;
 
-%Solution parameters
-	%parallelization 
-	md.cluster='none';
-	md.np=2;
-	md.time=1;
-	md.exclusive=0;
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
 
-	%statics
-	md.lowmem=1;
-	md.eps_rel=0.01;
-	md.eps_abs=10;
-	md.penalty_offset=3;
-	md.penalty_melting=10^7;
-	if md.numberofgrids<1000000,
-	md.sparsity=.001;
-	else
-	md.sparsity=.0001;
-	end
+disp('      creating temperature');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
 
-	%dynamics
-	md.dt=1*md.yts; %1 year
-	md.ndt=md.dt*10; 
-	md.artificial_diffusivity=1;
+disp('      creating flow law paramter');
+md.B=1.7687*10^8*ones(md.numberofgrids,1);
+md.n=3*ones(md.numberofelements,1);
 
-	%control
-	md.control_type={'drag'}; %'drag', 'B'
-	md.nsteps=5;
-	md.tolx=10^-4;
-	md.maxiter=20;
-	md.optscal=10;
-	md.fit='logarithmic'; %'absolute','relative','logarithmic'
-	md.meanvel=1000/md.yts; %1000 meters/year
-	md.epsvel=eps;
+disp('      creating accumulation rates');
+md.accumulation=0.2*ones(md.numberofgrids,1); %0.2m/a
+md.melting=0*ones(md.numberofgrids,1); %0m/a
 
-
-	disp('      creating thickness');
-	ymin=min(md.y);
-	ymax=max(md.y);
-	md.thickness=500*ones(md.numberofgrids,1);
-	md.firn_layer=0*ones(md.numberofgrids,1);
-	md.bed=-di*md.thickness;
-	md.surface=md.bed+md.thickness;
-
-	disp('      creating velocities');
-	md.vx_obs=zeros(md.numberofgrids,1);
-	md.vy_obs=zeros(md.numberofgrids,1);
-	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-	
-	disp('      creating drag');
-	md.drag_type=2; %0 none 1 plastic 2 viscous
-	md.drag=200*ones(md.numberofgrids,1); %q=1.
-	%Take care of iceshelves: no basal drag
-	pos=find(md.elementoniceshelf);
-	md.drag(md.elements(pos,:))=0;
-	md.p=ones(md.numberofelements,1);
-	md.q=ones(md.numberofelements,1);
-
-	disp('      creating temperature');
-	md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
-
-	disp('      creating flow law paramter');
-	%md.B=paterson(md.observed_temperature);
-	md.B=1.7687*10^8*ones(md.numberofgrids,1);
-	md.n=3*ones(md.numberofelements,1);
-
-	disp('      creating accumulation rates');
-	md.accumulation=0.2*ones(md.numberofgrids,1); %10m/a
-	md.melting=0*ones(md.numberofgrids,1); %1m/a
-
-	%Deal with boundary conditions:
-	
-	disp('      boundary conditions for diagnostic model: ');
-	%Build gridonicefront, array of boundary grids belonging to the icefront:
-	gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
-	gridonicefront=double(md.gridonboundary & gridinsideicefront);
-
-	md.gridondirichlet_diag=zeros(md.numberofgrids,1);
-	pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
-	md.dirichletvalues_diag=zeros(md.numberofgrids,2);
-
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_diag=md.segments(pos,:);
-	md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
-
-	disp('      boundary conditions for prognostic model: ');
-	md.gridondirichlet_prog=zeros(md.numberofgrids,1);
-	md.dirichletvalues_prog=zeros(md.numberofgrids,1);
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_prog=md.segments(pos,:);
-	md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
-	md.neumannvalues_prog(:)=NaN; %free radiation
-	
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_prog2=md.segments(pos,:);
-	md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
-	md.neumannvalues_prog2(:)=NaN; %free radiation
-	
-	disp('      boundary conditions for thermal model: ');
-	md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
-	md.dirichletvalues_thermal=md.observed_temperature;
-	md.geothermalflux=zeros(md.numberofgrids,1); 
-	pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
-
-
-	
-
-% Some Cielo code, ignore.
-if strcmp(md.cluster,'yes')
-	ServerDisconnect;
-end   
+disp('      boundary conditions ');
+md=SetMarineIceSheetBC(md,'Front.exp');
Index: /issm/trunk/test/Validation/EISMINT/Diagnostic2/Square.par
===================================================================
--- /issm/trunk/test/Validation/EISMINT/Diagnostic2/Square.par	(revision 899)
+++ /issm/trunk/test/Validation/EISMINT/Diagnostic2/Square.par	(revision 900)
@@ -1,122 +1,34 @@
-
 %Ok, start defining model parameters here
 
-%material parameters
-	md.g=9.81;
-	md.rho_ice=917;
-	md.rho_water=1028;
-	di=md.rho_ice/md.rho_water;
-	md.yts=365*24*3600;
-	md.heatcapacity=2009;
-	md.thermalconductivity=2.2; %W/mK
-	md.beta=9.8*10^-8;
+disp('      creating thickness');
+ymin=min(md.y);
+ymax=max(md.y);
+md.thickness=500*ones(md.numberofgrids,1);
+md.firn_layer=0*ones(md.numberofgrids,1);
+md.bed=-di*md.thickness;
+md.surface=md.bed+md.thickness;
 
-%Solution parameters
-	%parallelization 
-	md.cluster='none';
-	md.np=2;
-	md.time=1;
-	md.exclusive=0;
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
 
-	%statics
-	md.lowmem=1;
-	md.eps_rel=0.01;
-	md.eps_abs=10;
-	md.penalty_offset=3;
-	md.penalty_melting=10^7;
-	if md.numberofgrids<1000000,
-	md.sparsity=.001;
-	else
-	md.sparsity=.0001;
-	end
+disp('      creating temperature');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
 
-	%dynamics
-	md.dt=1*md.yts; %1 year
-	md.ndt=md.dt*10; 
-	md.artificial_diffusivity=1;
+disp('      creating flow law paramter');
+md.B=1.7687*10^8*ones(md.numberofgrids,1);
+md.n=3*ones(md.numberofelements,1);
 
-	%control
-	md.control_type={'drag'}; %'drag', 'B'
-	md.nsteps=5;
-	md.tolx=10^-4;
-	md.maxiter=20;
-	md.optscal=10;
-	md.fit='logarithmic'; %'absolute','relative','logarithmic'
-	md.meanvel=1000/md.yts; %1000 meters/year
-	md.epsvel=eps;
+disp('      creating accumulation rates');
+md.accumulation=0.2*ones(md.numberofgrids,1); %0.2m/a
+md.melting=0*ones(md.numberofgrids,1); %0m/a
 
-
-	disp('      creating thickness');
-	ymin=min(md.y);
-	ymax=max(md.y);
-	md.thickness=500*ones(md.numberofgrids,1);
-	md.firn_layer=0*ones(md.numberofgrids,1);
-	md.bed=-di*md.thickness;
-	md.surface=md.bed+md.thickness;
-
-	disp('      creating velocities');
-	md.vx_obs=zeros(md.numberofgrids,1);
-	md.vy_obs=zeros(md.numberofgrids,1);
-	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-	
-	disp('      creating drag');
-	md.drag_type=2; %0 none 1 plastic 2 viscous
-	md.drag=200*ones(md.numberofgrids,1); %q=1.
-	%Take care of iceshelves: no basal drag
-	pos=find(md.elementoniceshelf);
-	md.drag(md.elements(pos,:))=0;
-	md.p=ones(md.numberofelements,1);
-	md.q=ones(md.numberofelements,1);
-
-	disp('      creating temperature');
-	md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
-
-	disp('      creating flow law paramter');
-	%md.B=paterson(md.observed_temperature);
-	md.B=1.7687*10^8*ones(md.numberofgrids,1);
-	md.n=3*ones(md.numberofelements,1);
-
-	disp('      creating accumulation rates');
-	md.accumulation=0.2*ones(md.numberofgrids,1); %10m/a
-	md.melting=0*ones(md.numberofgrids,1); %1m/a
-
-	%Deal with boundary conditions:
-	
-	disp('      boundary conditions for diagnostic model: ');
-	%Build gridonicefront, array of boundary grids belonging to the icefront:
-	gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
-	gridonicefront=double(md.gridonboundary & gridinsideicefront);
-
-	md.gridondirichlet_diag=zeros(md.numberofgrids,1);
-	pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
-	md.dirichletvalues_diag=zeros(md.numberofgrids,2);
-
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_diag=md.segments(pos,:);
-	md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
-
-	disp('      boundary conditions for prognostic model: ');
-	md.gridondirichlet_prog=zeros(md.numberofgrids,1);
-	md.dirichletvalues_prog=zeros(md.numberofgrids,1);
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	pos=find(md.y==max(md.y)); %grids on the upper boundary condition
-	md.dirichletvalues_diag(pos,2)=400*(((md.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);
-	
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_prog2=md.segments(pos,:);
-	md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
-	md.neumannvalues_prog2(:)=NaN; %free radiation
-	
-	disp('      boundary conditions for thermal model: ');
-	md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
-	md.dirichletvalues_thermal=md.observed_temperature;
-	md.geothermalflux=zeros(md.numberofgrids,1); 
-	pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
-
-
-	
-
-% Some Cielo code, ignore.
-if strcmp(md.cluster,'yes')
-	ServerDisconnect;
-end   
+disp('      boundary conditions ');
+md=SetMarineIceSheetBC(md,'Front.exp');
+pos=find(md.y==max(md.y)); %grids on the upper boundary condition
+md.dirichletvalues_diag(pos,2)=400*(((md.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);
Index: /issm/trunk/test/Validation/EISMINT/MassConservation/Square.par
===================================================================
--- /issm/trunk/test/Validation/EISMINT/MassConservation/Square.par	(revision 899)
+++ /issm/trunk/test/Validation/EISMINT/MassConservation/Square.par	(revision 900)
@@ -1,134 +1,44 @@
-
 %Ok, start defining model parameters here
 
-%material parameters
-	md.g=9.81;
-	md.rho_ice=917;
-	md.rho_water=1028;
-	di=md.rho_ice/md.rho_water;
-	md.yts=365*24*3600;
-	md.heatcapacity=2009;
-	md.thermalconductivity=2.2; %W/mK
-	md.beta=9.8*10^-8;
+disp('      creating thickness');
+ymin=min(md.y);
+ymax=max(md.y);
+md.thickness=500*ones(md.numberofgrids,1);
+md.firn_layer=0*ones(md.numberofgrids,1);
+md.bed=-di*md.thickness;
+md.surface=md.bed+md.thickness;
 
-%Solution parameters
-	%parallelization 
-	md.cluster='none';
-	md.np=2;
-	md.time=1;
-	md.exclusive=0;
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
 
-	%statics
-	md.lowmem=1;
-	md.eps_rel=0.01;
-	md.eps_abs=10;
-	md.penalty_offset=3;
-	md.penalty_melting=10^7;
-	if md.numberofgrids<1000000,
-	md.sparsity=.001;
-	else
-	md.sparsity=.0001;
-	end
+disp('      creating temperature');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
 
-	%dynamics
-	md.dt=1*md.yts; %1 year
-	md.ndt=md.dt*10; 
-	md.artificial_diffusivity=1;
+disp('      creating flow law paramter');
+%md.B=paterson(md.observed_temperature);
+md.B=1.7687*10^8*ones(md.numberofgrids,1);
+md.n=3*ones(md.numberofelements,1);
 
-	%control
-	md.control_type={'drag'}; %'drag', 'B'
-	md.nsteps=5;
-	md.tolx=10^-4;
-	md.maxiter=20;
-	md.optscal=10;
-	md.fit='logarithmic'; %'absolute','relative','logarithmic'
-	md.meanvel=1000/md.yts; %1000 meters/year
-	md.epsvel=eps;
+disp('      creating accumulation rates');
+md.accumulation=0.0*ones(md.numberofgrids,1); %0m/a
+md.melting=0*ones(md.numberofgrids,1); %0m/a
 
+disp('      boundary conditions ');
+md=SetMarineIceSheetBC(md,'Front.exp');
 
-	disp('      creating thickness');
-	ymin=min(md.y);
-	ymax=max(md.y);
-	md.thickness=500*ones(md.numberofgrids,1);
-	md.firn_layer=0*ones(md.numberofgrids,1);
-	md.bed=-di*md.thickness;
-	md.surface=md.bed+md.thickness;
+%analytical test
+pos=find(md.y==200000); %grids on the upper boundary condition
+md.gridondirichlet_diag=ones(md.numberofgrids,1);
+md.dirichletvalues_diag(:,1)=0;
+md.dirichletvalues_diag(:,2)=-400;
 
-	disp('      creating velocities');
-	md.vx_obs=zeros(md.numberofgrids,1);
-	md.vy_obs=zeros(md.numberofgrids,1);
-	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-	
-	disp('      creating drag');
-	md.drag_type=2; %0 none 1 plastic 2 viscous
-	md.drag=200*ones(md.numberofgrids,1); %q=1.
-	%Take care of iceshelves: no basal drag
-	pos=find(md.elementoniceshelf);
-	md.drag(md.elements(pos,:))=0;
-	md.p=ones(md.numberofelements,1);
-	md.q=ones(md.numberofelements,1);
-
-	disp('      creating temperature');
-	md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
-
-	disp('      creating flow law paramter');
-	%md.B=paterson(md.observed_temperature);
-	md.B=1.7687*10^8*ones(md.numberofgrids,1);
-	md.n=3*ones(md.numberofelements,1);
-
-	disp('      creating accumulation rates');
-	md.accumulation=0.0*ones(md.numberofgrids,1); %10m/a
-	md.melting=0*ones(md.numberofgrids,1); %1m/a
-
-	%Deal with boundary conditions:
-	
-	disp('      boundary conditions for diagnostic model: ');
-	%Build gridonicefront, array of boundary grids belonging to the icefront:
-	gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
-	gridonicefront=double(md.gridonboundary & gridinsideicefront);
-
-	md.gridondirichlet_diag=zeros(md.numberofgrids,1);
-	pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
-	md.dirichletvalues_diag=zeros(md.numberofgrids,2);
-
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_diag=md.segments(pos,:);
-	md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
-
-	%analytical test
-	pos=find(md.y==200000); %grids on the upper boundary condition
-	md.gridondirichlet_diag=ones(md.numberofgrids,1);
-	md.dirichletvalues_diag(:,1)=0;
-	md.dirichletvalues_diag(:,2)=-400;
-
-	disp('      boundary conditions for prognostic model: ');
-	md.gridondirichlet_prog=zeros(md.numberofgrids,1);
-	md.dirichletvalues_prog=zeros(md.numberofgrids,1);
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_prog=md.segments(pos,:);
-	md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
-	md.neumannvalues_prog(:)=NaN; %free radiation
-	
-	%analytical test
-	pos=find(md.y==200000); %grids on the upper boundary condition
-	md.gridondirichlet_prog(pos)=1;
-	md.dirichletvalues_prog(pos)=500;
-
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_prog2=md.segments(pos,:);
-	md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
-	md.neumannvalues_prog2(:)=NaN; %free radiation
-	
-	disp('      boundary conditions for thermal model: ');
-	md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
-	md.dirichletvalues_thermal=md.observed_temperature;
-	md.geothermalflux=zeros(md.numberofgrids,1); 
-	pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
-
-
-	
-
-% Some Cielo code, ignore.
-if strcmp(md.cluster,'yes')
-	ServerDisconnect;
-end   
+%analytical test
+pos=find(md.y==200000); %grids on the upper boundary condition
+md.gridondirichlet_prog(pos)=1;
+md.dirichletvalues_prog(pos)=500;
Index: /issm/trunk/test/Validation/EISMINT/Transient/Square.par
===================================================================
--- /issm/trunk/test/Validation/EISMINT/Transient/Square.par	(revision 899)
+++ /issm/trunk/test/Validation/EISMINT/Transient/Square.par	(revision 900)
@@ -1,128 +1,45 @@
-
 %Ok, start defining model parameters here
 
-%material parameters
-	md.g=9.81;
-	md.rho_ice=917;
-	md.rho_water=1028;
-	di=md.rho_ice/md.rho_water;
-	md.yts=365*24*3600;
-	md.heatcapacity=2009;
-	md.thermalconductivity=2.2; %W/mK
-	md.beta=9.8*10^-8;
+disp('      creating thickness');
+ymin=min(md.y);
+ymax=max(md.y);
+md.thickness=500*ones(md.numberofgrids,1);
+md.firn_layer=0*ones(md.numberofgrids,1);
+md.bed=-di*md.thickness;
+md.surface=md.bed+md.thickness;
 
-%Solution parameters
-	%parallelization 
-	md.cluster='none';
-	md.np=2;
-	md.time=1;
-	md.exclusive=0;
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
 
-	%statics
-	md.lowmem=1;
-	md.eps_rel=0.01;
-	md.eps_abs=10;
-	md.penalty_offset=3;
-	md.penalty_melting=10^7;
-	if md.numberofgrids<1000000,
-	md.sparsity=.001;
-	else
-	md.sparsity=.0001;
-	end
+disp('      creating temperature');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
 
-	%dynamics
-	md.dt=1*md.yts; %1 year
-	md.ndt=md.dt*10; 
-	md.artificial_diffusivity=1;
+disp('      creating flow law paramter');
+md.B=1.7687*10^8*ones(md.numberofgrids,1);
+md.n=3*ones(md.numberofelements,1);
 
-	%control
-	md.control_type={'drag'}; %'drag', 'B'
-	md.nsteps=5;
-	md.tolx=10^-4;
-	md.maxiter=20;
-	md.optscal=10;
-	md.fit='logarithmic'; %'absolute','relative','logarithmic'
-	md.meanvel=1000/md.yts; %1000 meters/year
-	md.epsvel=eps;
+disp('      creating accumulation rates');
+md.accumulation=0.2*ones(md.numberofgrids,1); %0.2m/a
+md.melting=0*ones(md.numberofgrids,1); %0m/a
 
+disp('      boundary conditions ');
+md=SetMarineIceSheetBC(md,'Front.exp');
 
-	disp('      creating thickness');
-	ymin=min(md.y);
-	ymax=max(md.y);
-	md.thickness=500*ones(md.numberofgrids,1);
-	md.firn_layer=0*ones(md.numberofgrids,1);
-	md.bed=-di*md.thickness;
-	md.surface=md.bed+md.thickness;
+disp('      boundary conditions for prognostic model: ');
+md.gridondirichlet_prog=zeros(md.numberofgrids,1);
+md.dirichletvalues_prog=zeros(md.numberofgrids,1);
 
-	disp('      creating velocities');
-	md.vx_obs=zeros(md.numberofgrids,1);
-	md.vy_obs=zeros(md.numberofgrids,1);
-	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-	
-	disp('      creating drag');
-	md.drag_type=2; %0 none 1 plastic 2 viscous
-	md.drag=200*ones(md.numberofgrids,1); %q=1.
-	%Take care of iceshelves: no basal drag
-	pos=find(md.elementoniceshelf);
-	md.drag(md.elements(pos,:))=0;
-	md.p=ones(md.numberofelements,1);
-	md.q=ones(md.numberofelements,1);
+pos=find(md.y==max(md.y)); %grids on the upper boundary condition
+md.dirichletvalues_diag(pos,2)=400*(((md.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);
+md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
+md.neumannvalues_prog(:)=NaN; %free radiation
 
-	disp('      creating temperature');
-	md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
-
-	disp('      creating flow law paramter');
-	%md.B=paterson(md.observed_temperature);
-	md.B=1.7687*10^8*ones(md.numberofgrids,1);
-	md.n=3*ones(md.numberofelements,1);
-
-	disp('      creating accumulation rates');
-	md.accumulation=0.2*ones(md.numberofgrids,1); %10m/a
-	md.melting=0*ones(md.numberofgrids,1); %1m/a
-
-	%Deal with boundary conditions:
-	
-	disp('      boundary conditions for diagnostic model: ');
-	%Build gridonicefront, array of boundary grids belonging to the icefront:
-	gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
-	gridonicefront=double(md.gridonboundary & gridinsideicefront);
-
-	md.gridondirichlet_diag=zeros(md.numberofgrids,1);
-	pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
-	md.dirichletvalues_diag=zeros(md.numberofgrids,2);
-
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_diag=md.segments(pos,:);
-	md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
-
-	disp('      boundary conditions for prognostic model: ');
-	md.gridondirichlet_prog=zeros(md.numberofgrids,1);
-	md.dirichletvalues_prog=zeros(md.numberofgrids,1);
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	pos=find(md.y==max(md.y)); %grids on the upper boundary condition
-	md.dirichletvalues_diag(pos,2)=400*(((md.x(pos)-100000)/25000).^2-ones(size(pos,1),1)).*heaviside((1+eps)*ones(size(pos,1),1)-((md.x(pos)-100000)/25000).^2);
-	md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
-	md.neumannvalues_prog(:)=NaN; %free radiation
-	
-	pos=find(md.y==max(md.y)); %grids on the upper boundary condition
-	md.gridondirichlet_prog(pos)=1;
-	md.dirichletvalues_prog(pos)=500;
-
-	pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-	md.segmentonneumann_prog2=md.segments(pos,:);
-	md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
-	md.neumannvalues_prog2(:)=NaN; %free radiation
-	
-	disp('      boundary conditions for thermal model: ');
-	md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
-	md.dirichletvalues_thermal=md.observed_temperature;
-	md.geothermalflux=zeros(md.numberofgrids,1); 
-	pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
-
-
-	
-
-% Some Cielo code, ignore.
-if strcmp(md.cluster,'yes')
-	ServerDisconnect;
-end   
+pos=find(md.y==max(md.y)); %grids on the upper boundary condition
+md.gridondirichlet_prog(pos)=1;
+md.dirichletvalues_prog(pos)=500;
Index: /issm/trunk/test/Validation/ThermalTests/Melting/Square.par
===================================================================
--- /issm/trunk/test/Validation/ThermalTests/Melting/Square.par	(revision 899)
+++ /issm/trunk/test/Validation/ThermalTests/Melting/Square.par	(revision 900)
@@ -1,47 +1,47 @@
 %Ok, start defining model parameters here
 
-	%parallelization 
-	md.cluster='none';
+%parallelization 
+md.cluster='none';
 
-	disp('      creating thickness');
-	h=1000;
-	md.thickness=h*ones(md.numberofgrids,1);
-	md.firn_layer=10*ones(md.numberofgrids,1);
-	md.bed=-1000*ones(md.numberofgrids,1);
-	md.surface=md.bed+md.thickness;
+disp('      creating thickness');
+h=1000;
+md.thickness=h*ones(md.numberofgrids,1);
+md.firn_layer=10*ones(md.numberofgrids,1);
+md.bed=-1000*ones(md.numberofgrids,1);
+md.surface=md.bed+md.thickness;
 
-	disp('      creating velocities');
-	md.vx_obs=zeros(md.numberofgrids,1);
-	md.vy_obs=zeros(md.numberofgrids,1);
-	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-	
-	disp('      creating drag');
-	md.drag_type=2; %0 none 1 plastic 2 viscous
-	md.drag=200*ones(md.numberofgrids,1); %q=1.
-	%Take care of iceshelves: no basal drag
-	pos=find(md.elementoniceshelf);
-	md.drag(md.elements(pos,:))=0;
-	md.p=ones(md.numberofelements,1);
-	md.q=ones(md.numberofelements,1);
+disp('      creating velocities');
+md.vx_obs=zeros(md.numberofgrids,1);
+md.vy_obs=zeros(md.numberofgrids,1);
+md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
 
-	disp('      creating temperatures');
-	md.observed_temperature=273.15*ones(md.numberofgrids,1);
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
 
-	disp('      creating flow law paramter');
-	md.B=paterson(md.observed_temperature);
-	md.n=3*ones(md.numberofelements,1);
-	
-	disp('      creating accumulation rates');
-	md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
-	md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
+disp('      creating temperatures');
+md.observed_temperature=273.15*ones(md.numberofgrids,1);
 
-	%Deal with boundary conditions:
-	
-	disp('      boundary conditions for diagnostic model');
-	md=SetIceShelfBC(md,'Front.exp');
-	
-	disp('      boundary conditions for thermal model');
-	md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
-	md.dirichletvalues_thermal=md.observed_temperature;
-	md.geothermalflux=zeros(md.numberofgrids,1); 
-	pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=1*10^-3; %1 mW/m^2
+disp('      creating flow law paramter');
+md.B=paterson(md.observed_temperature);
+md.n=3*ones(md.numberofelements,1);
+
+disp('      creating accumulation rates');
+md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
+md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
+
+%Deal with boundary conditions:
+
+disp('      boundary conditions for diagnostic model');
+md=SetIceShelfBC(md,'Front.exp');
+
+disp('      boundary conditions for thermal model');
+md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
+md.dirichletvalues_thermal=md.observed_temperature;
+md.geothermalflux=zeros(md.numberofgrids,1); 
+pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=1*10^-3; %1 mW/m^2
Index: /issm/trunk/test/Validation/ThermalTests/Simpleadvection/Square.par
===================================================================
--- /issm/trunk/test/Validation/ThermalTests/Simpleadvection/Square.par	(revision 899)
+++ /issm/trunk/test/Validation/ThermalTests/Simpleadvection/Square.par	(revision 900)
@@ -1,47 +1,47 @@
 %Ok, start defining model parameters here
 
-	%parallelization 
-	md.cluster='none';
+%parallelization 
+md.cluster='none';
 
-	disp('      creating thickness');
-	h=1000;
-	md.thickness=h*ones(md.numberofgrids,1);
-	md.firn_layer=10*ones(md.numberofgrids,1);
-	md.bed=-1000*ones(md.numberofgrids,1);
-	md.surface=md.bed+md.thickness;
+disp('      creating thickness');
+h=1000;
+md.thickness=h*ones(md.numberofgrids,1);
+md.firn_layer=10*ones(md.numberofgrids,1);
+md.bed=-1000*ones(md.numberofgrids,1);
+md.surface=md.bed+md.thickness;
 
-	disp('      creating velocities');
-	md.vx_obs=zeros(md.numberofgrids,1);
-	md.vy_obs=zeros(md.numberofgrids,1);
-	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-	
-	disp('      creating drag');
-	md.drag_type=2; %0 none 1 plastic 2 viscous
-	md.drag=200*ones(md.numberofgrids,1); %q=1.
-	%Take care of iceshelves: no basal drag
-	pos=find(md.elementoniceshelf);
-	md.drag(md.elements(pos,:))=0;
-	md.p=ones(md.numberofelements,1);
-	md.q=ones(md.numberofelements,1);
+disp('      creating velocities');
+md.vx_obs=zeros(md.numberofgrids,1);
+md.vy_obs=zeros(md.numberofgrids,1);
+md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
 
-	disp('      creating temperatures');
-	md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
 
-	disp('      creating flow law paramter');
-	md.B=paterson(md.observed_temperature);
-	md.n=3*ones(md.numberofelements,1);
-	
-	disp('      creating accumulation rates');
-	md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
-	md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
+disp('      creating temperatures');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
 
-	%Deal with boundary conditions:
-	
-	disp('      boundary conditions for diagnostic model');
-	md=SetIceShelfBC(md,'Front.exp');
-	
-	disp('      boundary conditions for thermal model');
-	md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
-	md.dirichletvalues_thermal=md.observed_temperature;
-	md.geothermalflux=zeros(md.numberofgrids,1); 
-	pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
+disp('      creating flow law paramter');
+md.B=paterson(md.observed_temperature);
+md.n=3*ones(md.numberofelements,1);
+
+disp('      creating accumulation rates');
+md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
+md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
+
+%Deal with boundary conditions:
+
+disp('      boundary conditions for diagnostic model');
+md=SetIceShelfBC(md,'Front.exp');
+
+disp('      boundary conditions for thermal model');
+md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
+md.dirichletvalues_thermal=md.observed_temperature;
+md.geothermalflux=zeros(md.numberofgrids,1); 
+pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
Index: /issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux/Square.par
===================================================================
--- /issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux/Square.par	(revision 899)
+++ /issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux/Square.par	(revision 900)
@@ -1,47 +1,47 @@
 %Ok, start defining model parameters here
 
-	%parallelization 
-	md.cluster='none';
+%parallelization 
+md.cluster='none';
 
-	disp('      creating thickness');
-	h=1000;
-	md.thickness=h*ones(md.numberofgrids,1);
-	md.firn_layer=10*ones(md.numberofgrids,1);
-	md.bed=zeros(md.numberofgrids,1);
-	md.surface=md.thickness;
+disp('      creating thickness');
+h=1000;
+md.thickness=h*ones(md.numberofgrids,1);
+md.firn_layer=10*ones(md.numberofgrids,1);
+md.bed=zeros(md.numberofgrids,1);
+md.surface=md.thickness;
 
-	disp('      creating velocities');
-	md.vx_obs=zeros(md.numberofgrids,1);
-	md.vy_obs=zeros(md.numberofgrids,1);
-	md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-	
-	disp('      creating drag');
-	md.drag_type=2; %0 none 1 plastic 2 viscous
-	md.drag=200*ones(md.numberofgrids,1); %q=1.
-	%Take care of iceshelves: no basal drag
-	pos=find(md.elementoniceshelf);
-	md.drag(md.elements(pos,:))=0;
-	md.p=ones(md.numberofelements,1);
-	md.q=ones(md.numberofelements,1);
+disp('      creating velocities');
+md.vx_obs=zeros(md.numberofgrids,1);
+md.vy_obs=zeros(md.numberofgrids,1);
+md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
 
-	disp('      creating temperatures');
-	md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=200*ones(md.numberofgrids,1); %q=1.
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
 
-	disp('      creating flow law paramter');
-	md.B=paterson(md.observed_temperature);
-	md.n=3*ones(md.numberofelements,1);
-	
-	disp('      creating accumulation rates');
-	md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
-	md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
+disp('      creating temperatures');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
 
-	%Deal with boundary conditions:
-	
-	disp('      boundary conditions for diagnostic model');
-	md=SetIceShelfBC(md,'Front.exp');
-	
-	disp('      boundary conditions for thermal model');
-	md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
-	md.dirichletvalues_thermal=md.observed_temperature;
-	md.geothermalflux=zeros(md.numberofgrids,1); 
-	pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=0.1; %100 mW/m^2
+disp('      creating flow law paramter');
+md.B=paterson(md.observed_temperature);
+md.n=3*ones(md.numberofelements,1);
+
+disp('      creating accumulation rates');
+md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
+md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
+
+%Deal with boundary conditions:
+
+disp('      boundary conditions for diagnostic model');
+md=SetIceShelfBC(md,'Front.exp');
+
+disp('      boundary conditions for thermal model');
+md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
+md.dirichletvalues_thermal=md.observed_temperature;
+md.geothermalflux=zeros(md.numberofgrids,1); 
+pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=0.1; %100 mW/m^2