Index: /issm/workshop/16_ISMIP/test2.m
===================================================================
--- /issm/workshop/16_ISMIP/test2.m	(revision 11052)
+++ /issm/workshop/16_ISMIP/test2.m	(revision 11052)
@@ -0,0 +1,41 @@
+%Create model
+md=model;
+
+%Create mesh
+L=100000; %in m
+nx=20; %numberof nodes in x direction
+ny=20;
+md=squaremesh(md,L,L,nx,ny);
+
+md=setmask(md,'',''); %ice sheet test
+md=parameterize(md,'./test2.par');
+md=extrude(md,5,1);
+md=setflowequation(md,'pattyn','all');
+
+%Boundary conditions
+md.diagnostic.spcvx(:)=NaN;
+md.diagnostic.spcvy(:)=NaN;
+md.diagnostic.spcvz(:)=NaN;
+%Create dirichlet on the bed if no slip
+pos=find(md.mesh.vertexonbed);
+md.diagnostic.spcvx(pos)=0;
+md.diagnostic.spcvy(pos)=0;
+md.diagnostic.spcvz(pos)=0;
+
+%Create MPCs to have periodic boundary conditions
+posx=find(md.mesh.x==0);
+posx2=find(md.mesh.x==max(md.mesh.x));
+posy=find(md.mesh.y==0); %Don't take the same grids two times
+posy2=find(md.mesh.y==max(md.mesh.y));
+md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
+md.prognostic.vertex_pairing=[posx,posx2;posy,posy2];
+
+%Transient parameters
+md.timestepping.time_step=4;
+md.timestepping.final_time=4*20;
+md.transient.isthermal=0;
+
+%Compute the diagnostic
+md.cluster=generic('name',oshostname(),'np',2);
+md.verbose=verbose('convergence',true,'solution',true);
+md=solve(md,TransientSolutionEnum);
Index: /issm/workshop/16_ISMIP/test2.par
===================================================================
--- /issm/workshop/16_ISMIP/test2.par	(revision 11052)
+++ /issm/workshop/16_ISMIP/test2.par	(revision 11052)
@@ -0,0 +1,24 @@
+disp('      creating geometry');
+L=100000;
+alpha=3*pi/180;
+md.geometry.surface=-md.mesh.x*tan(alpha);
+md.geometry.bed=md.geometry.surface-1000+100*exp(-((md.mesh.x-L/2).^2+(md.mesh.y-L/2).^2)/(10000^2)); 
+md.geometry.thickness=md.geometry.surface-md.geometry.bed;
+
+disp('      creating drag');
+md.friction.coefficient=sqrt(365.25*24*3600/(1000*2.140373*10^-7))*ones(md.mesh.numberofvertices,1);
+md.friction.p=ones(md.mesh.numberofelements,1);
+md.friction.q=zeros(md.mesh.numberofelements,1);
+
+disp('      creating flow law paramter');
+md.materials.rheology_B=(1/(2.140373*10^-7/(365.25*24*3600)))*ones(md.mesh.numberofvertices,1);
+md.materials.rheology_n=1*ones(md.mesh.numberofelements,1);
+
+disp('      boundary conditions for diagnostic model');
+md=SetIceSheetBC(md);
+
+%Field initialization
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+md.initialization.pressure=zeros(md.mesh.numberofvertices,1);