Index: /issm/trunk/test/Exp/RoundFrontEISMINT.exp
===================================================================
--- /issm/trunk/test/Exp/RoundFrontEISMINT.exp	(revision 5117)
+++ /issm/trunk/test/Exp/RoundFrontEISMINT.exp	(revision 5117)
@@ -0,0 +1,47 @@
+## Name:
+## Icon:0
+# Points Count Value
+41 1
+# X pos Y pos
+1000000.000000 0.000000
+987688.340595 156434.465040
+951056.516295 309016.994375
+891006.524188 453990.499740
+809016.994375 587785.252292
+707106.781187 707106.781187
+587785.252292 809016.994375
+453990.499740 891006.524188
+309016.994375 951056.516295
+156434.465040 987688.340595
+0.000000 1000000.000000
+-156434.465040 987688.340595
+-309016.994375 951056.516295
+-453990.499740 891006.524188
+-587785.252292 809016.994375
+-707106.781187 707106.781187
+-809016.994375 587785.252292
+-891006.524188 453990.499740
+-951056.516295 309016.994375
+-987688.340595 156434.465040
+-1000000.000000 0.000000
+-987688.340595 -156434.465040
+-951056.516295 -309016.994375
+-891006.524188 -453990.499740
+-809016.994375 -587785.252292
+-707106.781187 -707106.781187
+-587785.252292 -809016.994375
+-453990.499740 -891006.524188
+-309016.994375 -951056.516295
+-156434.465040 -987688.340595
+-0.000000 -1000000.000000
+156434.465040 -987688.340595
+309016.994375 -951056.516295
+453990.499740 -891006.524188
+587785.252292 -809016.994375
+707106.781187 -707106.781187
+809016.994375 -587785.252292
+891006.524188 -453990.499740
+951056.516295 -309016.994375
+987688.340595 -156434.465040
+1000000.000000 0.000000
+
Index: /issm/trunk/test/NightlyRun/test1205.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1205.m	(revision 5117)
+++ /issm/trunk/test/NightlyRun/test1205.m	(revision 5117)
@@ -0,0 +1,87 @@
+%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
+numlayers=10;
+resolution=30000;
+
+%To begin with the numerical model
+md=model;
+md=roundmesh(md,750000,resolution);
+md=geography(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
+md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
+
+%Calculation of the analytical 2d velocity field
+constant=0.3;
+vx_obs=constant/2*md.x.*(md.thickness).^-1;
+vy_obs=constant/2*md.y.*(md.thickness).^-1;
+vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2));
+
+%We extrude the model to have a 3d model
+md=extrude(md,numlayers,1);
+md=setelementstype(md,'hutter','all');
+
+%Spc the grids on the bed
+pos=find(md.gridonbed);
+md.spcvelocity(pos,1:3)=1;
+
+%Now we can solve the problem 
+md=solve(md,'analysis_type',DiagnosticSolutionEnum);
+
+%Calculate the depth averaged velocity field (2d):
+vx=PatchToVec(md.results.DiagnosticSolution.Vx);
+vy=PatchToVec(md.results.DiagnosticSolution.Vy);
+vel=zeros(md.numberofgrids2d,1);
+
+grid_vel=0;
+for i=1:md.numberofgrids2d
+	for j=1:(md.numlayers-1)
+		grid_vel=grid_vel+1/(2*(md.numlayers-1))*(sqrt(vx(i+j*md.numberofgrids2d,1).^2+...
+			vy(i+j*md.numberofgrids2d,1).^2)+...
+			sqrt(vx(i+(j-1)*md.numberofgrids2d,1).^2+vy(i+(j-1)*md.numberofgrids2d,1).^2));
+	end
+	vel(i,1)=grid_vel;
+	grid_vel=0;
+end
+
+%Plot of the velocity from the exact and calculated solutions
+figure(1)
+subplot(2,2,1)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+vel,'FaceColor','interp','EdgeColor','none');
+title('Modelled velocity','FontSize',14,'FontWeight','bold')
+colorbar; 
+caxis([0 200]);
+   
+subplot(2,2,2)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+vel_obs,'FaceColor','interp','EdgeColor','none');
+title('Analytical velocity','FontSize',14,'FontWeight','bold')
+colorbar; 
+caxis([0 200]);
+
+subplot(2,2,3)
+hold on;
+plot(sqrt((md.x(1:md.numberofgrids2d)).^2+(md.y(1:md.numberofgrids2d)).^2),vel,'r.');
+plot(sqrt((md.x2d).^2+(md.y2d).^2),vel_obs,'b.');
+title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold');
+xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold');
+ylabel('velocity [m/yr]','FontSize',14,'FontWeight','bold');
+legend('calculated velocity','exact velocity');
+axis([0 750000 0 200]);
+hold off;
+
+subplot(2,2,4)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none');
+title('Relative misfit [%]','FontSize',14,'FontWeight','bold')
+colorbar;
+caxis([0 100]);
+
+%Fields and tolerances to track changes
+field_names     ={ ...
+	'Vx','Vy','Vel', ...
+};
+field_tolerances={...
+	1e-13,1e-13,1e-13, ...
+};
+field_values={
+	vx,vy,vel, ...
+};
Index: /issm/trunk/test/NightlyRun/test1206.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1206.m	(revision 5117)
+++ /issm/trunk/test/NightlyRun/test1206.m	(revision 5117)
@@ -0,0 +1,87 @@
+%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
+numlayers=10;
+resolution=30000;
+
+%To begin with the numerical model
+md=model;
+md=roundmesh(md,750000,resolution);
+md=geography(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
+md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
+
+%Calculation of the analytical 2d velocity field
+constant=0.3;
+vx_obs=constant/2*md.x.*(md.thickness).^-1;
+vy_obs=constant/2*md.y.*(md.thickness).^-1;
+vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2));
+
+%We extrude the model to have a 3d model
+md=extrude(md,numlayers,1);
+md=setelementstype(md,'pattyn','all');
+
+%Spc the grids on the bed
+pos=find(md.gridonbed);
+md.spcvelocity(pos,1:3)=1;
+
+%Now we can solve the problem 
+md=solve(md,'analysis_type',DiagnosticSolutionEnum);
+
+%Calculate the depth averaged velocity field (2d):
+vx=PatchToVec(md.results.DiagnosticSolution.Vx);
+vy=PatchToVec(md.results.DiagnosticSolution.Vy);
+vel=zeros(md.numberofgrids2d,1);
+
+grid_vel=0;
+for i=1:md.numberofgrids2d
+	for j=1:(md.numlayers-1)
+		grid_vel=grid_vel+1/(2*(md.numlayers-1))*(sqrt(vx(i+j*md.numberofgrids2d,1).^2+...
+			vy(i+j*md.numberofgrids2d,1).^2)+...
+			sqrt(vx(i+(j-1)*md.numberofgrids2d,1).^2+vy(i+(j-1)*md.numberofgrids2d,1).^2));
+	end
+	vel(i,1)=grid_vel;
+	grid_vel=0;
+end
+
+%Plot of the velocity from the exact and calculated solutions
+figure(1)
+subplot(2,2,1)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+vel,'FaceColor','interp','EdgeColor','none');
+title('Modelled velocity','FontSize',14,'FontWeight','bold')
+colorbar; 
+caxis([0 200]);
+   
+subplot(2,2,2)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+vel_obs,'FaceColor','interp','EdgeColor','none');
+title('Analytical velocity','FontSize',14,'FontWeight','bold')
+colorbar; 
+caxis([0 200]);
+
+subplot(2,2,3)
+hold on;
+plot(sqrt((md.x(1:md.numberofgrids2d)).^2+(md.y(1:md.numberofgrids2d)).^2),vel,'r.');
+plot(sqrt((md.x2d).^2+(md.y2d).^2),vel_obs,'b.');
+title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold');
+xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold');
+ylabel('velocity [m/yr]','FontSize',14,'FontWeight','bold');
+legend('calculated velocity','exact velocity');
+axis([0 750000 0 200]);
+hold off;
+
+subplot(2,2,4)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none');
+title('Relative misfit [%]','FontSize',14,'FontWeight','bold')
+colorbar;
+caxis([0 100]);
+
+%Fields and tolerances to track changes
+field_names     ={ ...
+	'Vx','Vy','Vel', ...
+};
+field_tolerances={...
+	1e-13,1e-13,1e-13, ...
+};
+field_values={
+	vx,vy,vel, ...
+};
Index: /issm/trunk/test/NightlyRun/test1207.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1207.m	(revision 5117)
+++ /issm/trunk/test/NightlyRun/test1207.m	(revision 5117)
@@ -0,0 +1,87 @@
+%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
+numlayers=10;
+resolution=30000;
+
+%To begin with the numerical model
+md=model;
+md=roundmesh(md,750000,resolution);
+md=geography(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
+md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
+
+%Calculation of the analytical 2d velocity field
+constant=0.3;
+vx_obs=constant/2*md.x.*(md.thickness).^-1;
+vy_obs=constant/2*md.y.*(md.thickness).^-1;
+vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2));
+
+%We extrude the model to have a 3d model
+md=extrude(md,numlayers,1);
+md=setelementstype(md,'pattyn','all','stokes','all');
+
+%Spc the grids on the bed
+pos=find(md.gridonbed);
+md.spcvelocity(pos,1:3)=1;
+
+%Now we can solve the problem 
+md=solve(md,'analysis_type',DiagnosticSolutionEnum);
+
+%Calculate the depth averaged velocity field (2d):
+vx=PatchToVec(md.results.DiagnosticSolution.Vx);
+vy=PatchToVec(md.results.DiagnosticSolution.Vy);
+vel=zeros(md.numberofgrids2d,1);
+
+grid_vel=0;
+for i=1:md.numberofgrids2d
+	for j=1:(md.numlayers-1)
+		grid_vel=grid_vel+1/(2*(md.numlayers-1))*(sqrt(vx(i+j*md.numberofgrids2d,1).^2+...
+			vy(i+j*md.numberofgrids2d,1).^2)+...
+			sqrt(vx(i+(j-1)*md.numberofgrids2d,1).^2+vy(i+(j-1)*md.numberofgrids2d,1).^2));
+	end
+	vel(i,1)=grid_vel;
+	grid_vel=0;
+end
+
+%Plot of the velocity from the exact and calculated solutions
+figure(1)
+subplot(2,2,1)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+vel,'FaceColor','interp','EdgeColor','none');
+title('Modelled velocity','FontSize',14,'FontWeight','bold')
+colorbar; 
+caxis([0 200]);
+   
+subplot(2,2,2)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+vel_obs,'FaceColor','interp','EdgeColor','none');
+title('Analytical velocity','FontSize',14,'FontWeight','bold')
+colorbar; 
+caxis([0 200]);
+
+subplot(2,2,3)
+hold on;
+plot(sqrt((md.x(1:md.numberofgrids2d)).^2+(md.y(1:md.numberofgrids2d)).^2),vel,'r.');
+plot(sqrt((md.x2d).^2+(md.y2d).^2),vel_obs,'b.');
+title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold');
+xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold');
+ylabel('velocity [m/yr]','FontSize',14,'FontWeight','bold');
+legend('calculated velocity','exact velocity');
+axis([0 750000 0 200]);
+hold off;
+
+subplot(2,2,4)
+p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
+abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none');
+title('Relative misfit [%]','FontSize',14,'FontWeight','bold')
+colorbar;
+caxis([0 100]);
+
+%Fields and tolerances to track changes
+field_names     ={ ...
+	'Vx','Vy','Vel', ...
+};
+field_tolerances={...
+	1e-13,1e-13,1e-13, ...
+};
+field_values={
+	vx,vy,vel, ...
+};
Index: /issm/trunk/test/Par/ISMIPC.par
===================================================================
--- /issm/trunk/test/Par/ISMIPC.par	(revision 5116)
+++ /issm/trunk/test/Par/ISMIPC.par	(revision 5117)
@@ -9,5 +9,5 @@
 disp('      creating drag');
 md.drag_type=2; %0 none 1 plastic 2 viscous
-md.drag_coefficient=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x/3)).*sin(md.y*2*pi/max(md.x/3)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
+md.drag_coefficient=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x/2)).*sin(md.y*2*pi/max(md.x/2)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
 %md.drag_coefficient=sqrt(md.yts.*(1000+1000*sin(md.x*2*pi/max(md.x)).*sin(md.y*2*pi/max(md.x)))./(md.g*(md.rho_ice*md.thickness+md.rho_water*md.bed)));
 %Take care of iceshelves: no basal drag
Index: /issm/trunk/test/Par/RoundSheetStaticEISMINT.par
===================================================================
--- /issm/trunk/test/Par/RoundSheetStaticEISMINT.par	(revision 5117)
+++ /issm/trunk/test/Par/RoundSheetStaticEISMINT.par	(revision 5117)
@@ -0,0 +1,60 @@
+disp('      creating thickness');
+hmin=0.01;
+hmax=2756.7;
+radius=(sqrt((md.x).^2+(md.y).^2));
+radiusmax=max(radius);
+md.thickness=hmin*ones(size(md.x,1),1)+hmax*(4*((1/2)^(4/3)*ones(size(md.x,1),1)-((radius)./(2*radiusmax)).^(4/3))).^(3/8);
+md.firn_layer=10*ones(md.numberofgrids,1);
+md.bed=0*md.thickness;
+md.surface=md.bed+md.thickness;
+
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag_coefficient=20*ones(md.numberofgrids,1); %q=1. %no drag is specified in the analytical solution
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag_coefficient(md.elements(pos,:))=0;
+md.drag_p=ones(md.numberofelements,1);
+md.drag_q=ones(md.numberofelements,1);
+
+disp('      creating temperatures');
+tmin=238.15; %K
+st=1.67*10^-2/1000; %k/m;
+md.observed_temperature=(tmin+st*radius);
+md.temperature=md.observed_temperature;
+md.geothermalflux=4.2*10^-2*ones(md.numberofgrids,1);
+
+disp('      creating flow law paramter');
+md.rheology_B=6.81*10^(7)*ones(md.numberofgrids,1); %to have the same B as the analytical solution 
+md.rheology_n=3*ones(md.numberofelements,1);
+
+disp('      creating accumulation rates');
+acc_max=0.5; %m/yr
+sb=10^-2/1000; %m/yr/m
+rel=450*1000; %m
+md.accumulation_rate=min(acc_max,sb*(rel-radius));
+
+disp('      creating velocities');
+constant=0.3;
+md.vx_obs=constant/2*md.x.*(md.thickness).^-1;
+md.vy_obs=constant/2*md.y.*(md.thickness).^-1;
+md.vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2));
+md.vx=zeros(md.numberofgrids,1);
+md.vy=zeros(md.numberofgrids,1);
+md.vz=zeros(md.numberofgrids,1);
+md.pressure=zeros(md.numberofgrids,1);
+
+%Deal with boundary conditions:
+disp('      boundary conditions for diagnostic model: ');
+md=SetMarineIceSheetBC(md,'../Exp/RoundFrontEISMINT.exp');
+
+radius=sqrt((md.x).*md.x+(md.y).*md.y);
+pos=find(radius==min(radius));
+md.x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center
+
+md.spcvelocity(pos,1:3)=1;
+md.spcvelocity(pos,4:6)=0;
+
+%parallel options
+md.np=6;
+md.cluster=oshostname();