Index: /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPAback.m
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPAback.m	(revision 14143)
+++ /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPAback.m	(revision 14143)
@@ -0,0 +1,44 @@
+L=80000; % 5000,10000,20000,40000,80000,160000
+
+%Create model
+md=model;
+
+%Create mesh
+nx=20; %number of nodes in x direction
+ny=20;
+md=squaremesh(model,L,L,nx,ny);
+
+%Define mask
+md=setmask(md,'',''); %ice sheet test
+
+%Parameterize model
+md=parameterize(md,'./ISMIPAback.par');
+
+%Extrude to have a 3d model
+md=extrude(md,5,1);
+
+%Choose flow approximation
+md=setflowequation(md,'pattyn','all');
+
+%Create dirichlet on the bed only
+md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+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 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same grids two times
+posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x));
+md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
+
+%Set cluster and print messages
+md.cluster=generic('name',oshostname(),'np',2);
+md.verbose=verbose('convergence',true,'solution',true);
+
+%Solve 
+md=solve(md,DiagnosticSolutionEnum);
Index: /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPAback.par
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPAback.par	(revision 14143)
+++ /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPAback.par	(revision 14143)
@@ -0,0 +1,20 @@
+%Create geometry: bed, thickness and surface
+disp('      creating thickness');
+md.geometry.surface=-md.mesh.x*tan(0.5*pi/180);
+md.geometry.bed=md.geometry.surface-1000+500*sin(md.mesh.x*2*pi/max(md.mesh.x)).*sin(md.mesh.y*2*pi/max(md.mesh.x));
+md.geometry.thickness=md.geometry.surface-md.geometry.bed;
+
+%Create drag parameters (won't be used as ice is frozen on the bed)
+disp('      creating drag');
+md.friction.coefficient=200*ones(md.mesh.numberofvertices,1); %q=1.
+md.friction.p=ones(md.mesh.numberofelements,1);
+md.friction.q=ones(md.mesh.numberofelements,1);
+
+%Create materials parameters: rheology parameters
+disp('      creating flow law paramter');
+md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+
+%Create generic boundary conditions
+disp('      boundary conditions for diagnostic model');
+md=SetIceSheetBC(md);
Index: /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPFback.m
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPFback.m	(revision 14143)
+++ /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPFback.m	(revision 14143)
@@ -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,'./ISMIPFback.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/2012/Talks/11_ISMIP/Files/ISMIPFback.par
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPFback.par	(revision 14143)
+++ /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPFback.par	(revision 14143)
@@ -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);