Index: /issm/workshop/12_Inverse_methods/runme.m
===================================================================
--- /issm/workshop/12_Inverse_methods/runme.m	(revision 11049)
+++ /issm/workshop/12_Inverse_methods/runme.m	(revision 11050)
@@ -1,3 +1,3 @@
-step=1;
+step=3;
 if step==1
 	%Generate observation
Index: /issm/workshop/16_ISMIP/ISMIPA.m
===================================================================
--- /issm/workshop/16_ISMIP/ISMIPA.m	(revision 11050)
+++ /issm/workshop/16_ISMIP/ISMIPA.m	(revision 11050)
@@ -0,0 +1,42 @@
+L=80000; % 5000,10000,20000,40000,80000,160000
+
+%Create model
+md=model;
+
+%Create mesh with square mesh
+......
+
+%Define mask 
+......
+
+%Parameterize model
+......
+
+%Extrude to have a 3d model
+......
+
+%Choose flow approximation
+......
+
+%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 
+......
Index: /issm/workshop/16_ISMIP/ISMIPA.par
===================================================================
--- /issm/workshop/16_ISMIP/ISMIPA.par	(revision 11050)
+++ /issm/workshop/16_ISMIP/ISMIPA.par	(revision 11050)
@@ -0,0 +1,20 @@
+%Create geometry: bed, thickness and surface
+disp('      creating geometry: thickness, bed and surface');
+md.geometry.surface= ......
+md.geometry.bed= ......
+md.geometry.thickness= ......
+
+%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/16_ISMIP/ISMIPF.m
===================================================================
--- /issm/workshop/16_ISMIP/ISMIPF.m	(revision 11050)
+++ /issm/workshop/16_ISMIP/ISMIPF.m	(revision 11050)
@@ -0,0 +1,47 @@
+%Create model
+md=model;
+
+%Create mesh with square mesh
+......
+
+%Set mask for grounded ice
+......
+
+%Parameterize model
+......
+
+%Extrude to create 3d model
+......
+
+%Set approximation used
+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];
+
+%Create transient parameters: timestepping and options for transient
+......
+......
+......
+
+%Compute the diagnostic
+md.cluster=generic('name',oshostname(),'np',2);
+md.verbose=verbose('convergence',true,'solution',true);
+
+%Solve model
+......
Index: /issm/workshop/16_ISMIP/ISMIPF.par
===================================================================
--- /issm/workshop/16_ISMIP/ISMIPF.par	(revision 11050)
+++ /issm/workshop/16_ISMIP/ISMIPF.par	(revision 11050)
@@ -0,0 +1,24 @@
+disp('      creating geometry: thickness, bed, surface');
+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);