Index: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPA.m
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPA.m	(revision 14124)
+++ 	(revision )
@@ -1,44 +1,0 @@
-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,'./ISMIPA.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: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPA.par
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPA.par	(revision 14124)
+++ 	(revision )
@@ -1,20 +1,0 @@
-%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: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPA_jgr.par
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPA_jgr.par	(revision 14124)
+++ 	(revision )
@@ -1,24 +1,0 @@
-%Ok, start defining model parameters here
-
-disp('      creating thickness');
-md.surface=-md.x*tan(0.5*pi/180);
-md.bed=md.surface-1000+500*sin(md.x*2*pi/max(md.x)).*sin(md.y*2*pi/max(md.x));
-md.thickness=md.surface-md.bed;
-md.firn_layer=0*ones(md.numberofgrids,1);
-
-disp('      creating drag');
-md.drag_type=2; %0 none 1 plastic 2 viscous
-md.drag_coefficient=200*ones(md.numberofgrids,1); %q=1.
-%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 flow law paramter');
-md.rheology_B=6.8067*10^7*ones(md.numberofgrids,1);
-md.rheology_n=3*ones(md.numberofelements,1);
-
-disp('      boundary conditions for diagnostic model');
-%Create grid on boundary fist (because we cannot use mesh)
-md=SetIceSheetBC(md);
Index: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPC_jgr.m
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPC_jgr.m	(revision 14124)
+++ 	(revision )
@@ -1,252 +1,0 @@
-%This test is a test from the ISMP-HOM Intercomparison project Pattyn and Payne 2006
-steps=[6];
-L=80000; % 5000,10000,20000,40000,80000,160000
-
-%Hard coded parameters
-clustername='pfe';
-interactive=1;
-printingflag=false;
-%Cluster parameters{{{1
-if strcmpi(clustername,'pfe');
-	cluster=pfe;
-	cluster.numnodes=40;
-	cluster.cpuspernode=12;
-	cluster.time=2*60;
-	cluster.interactive=interactive;
-	cluster.port=1024;
-	cluster.processor='wes';
-elseif strcmpi(clustername,'cosmos')
-	cluster=cosmos('np',100,'queue','shortq','time',1*60,'login','morlighe');
-elseif strcmpi(clustername,oshostname()),
-	cluster=eval([ oshostname() '(''np'',14);']);
-else
-	error('cluster not supported yet');
-end%}}}
-%Checks{{{
-if ~ismember(L,[5000,10000,20000,40000,80000,160000]),
-	error('L not supported (5000,10000,20000,40000,80000,160000)');
-end
-%}}}
-
-%Run Steps
-org=organizer('repository','./Models','prefix',['ISMIPC_' num2str(L) '.'],'steps',steps);
-if perform(org,'SetupHO'),% {{{1
-
-	nx=100; %number of nodes in x direction
-	ny=100;
-	md=squaremesh(model,L,L,nx,ny);
-	md=geography(md,'',''); %ice sheet test
-	md=parameterize(md,'./Par/ISMIPC.par');
-	md=extrude(md,20,1);
-	md=setelementstype(md,'pattyn','all');
-
-	%Create MPCs to have periodic boundary conditions
-	md.spcvelocity=zeros(md.numberofgrids,6);
-	posx=find(md.x==0 & md.mesh.y~=0 & md.mesh.y~=L);
-	posx2=find(md.x==L & md.mesh.y~=0 & md.mesh.y~=L);
-	posy=find(md.mesh.y==0 & md.x~=0 & md.x~=L); %Don't take the same grids two times
-	posy2=find(md.mesh.y==L & md.x~=0 & md.x~=L);
-	md.penalties=[posx,posx2;posy,posy2];
-
-	%Add spc on the corners
-	pos=find((md.x==0 | md.x==L) & (md.mesh.y==0 | md.mesh.y==L) & md.gridonbed);
-	md.spcvelocity(pos,1:2)=1;
-	if(L==5000),
-		md.spcvelocity(pos,4)=15.66;
-		md.spcvelocity(pos,5)=-0.1967;
-	elseif(L==10000),
-		md.spcvelocity(pos,4)=16.04;
-		md.spcvelocity(pos,5)=-0.1977;
-	elseif(L==20000),
-		md.spcvelocity(pos,4)=16.53;
-		md.spcvelocity(pos,5)=-1.27;
-	elseif(L==40000),
-		md.spcvelocity(pos,4)=17.23;
-		md.spcvelocity(pos,5)=-3.17;
-	elseif(L==80000),
-		md.spcvelocity(pos,4)=16.68;
-		md.spcvelocity(pos,5)=-2.69;
-	elseif(L==160000),
-		md.spcvelocity(pos,4)=16.03;
-		md.spcvelocity(pos,5)=-1.27;
-	end
-
-	%Spc the bed at zero for vz
-	pos=find(md.gridonbed);
-	md.spcvelocity(pos,3)=1;
-
-	savemodel(org,md);
-end %}}}
-if perform(org,'DiagHO'),% {{{1
-
-	md=loadmodel(org,'SetupHO');
-
-	md=setcluster(md,cluster);
-	md.verbose=verbose('convergence',true,'solution',true);
-	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
-
-	savemodel(org,md);
-end %}}}
-
-if perform(org,'SetupFS'),% {{{1
-
-	md=mesh(model,['./Exp/Square_' num2str(L) '.exp'],L/50);
-	md=geography(md,'',''); %ice sheet test
-	md=parameterize(md,'./Par/ISMIPC_stokes.par');
-	md=extrude(md,15,1);
-	md=setelementstype(md,'stokes','all');
-
-	%Create MPCs to have periodic boundary conditions
-	md.spcvelocity=zeros(md.numberofgrids,6);
-
-	%Add spc on the corners
-	pos=find((md.x==0 | md.x==max(md.x) | md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
-	md.spcvelocity(pos,1:3)=1;
-	if(L==5000),
-		md.spcvelocity(pos,4)=15.76;
-		md.spcvelocity(pos,5)=-0.1967;
-	elseif(L==10000),
-		md.spcvelocity(pos,4)=16.04;
-		md.spcvelocity(pos,5)=-0.1977;
-	elseif(L==20000),
-		md.spcvelocity(pos,4)=16.53;
-		md.spcvelocity(pos,5)=-1.27;
-	elseif(L==40000),
-		md.spcvelocity(pos,4)=17.23;
-		md.spcvelocity(pos,5)=-3.17;
-	elseif(L==80000),
-		md.spcvelocity(pos,4)=16.68;
-		md.spcvelocity(pos,5)=-2.69;
-	elseif(L==160000),
-		md.spcvelocity(pos,4)=16.03;
-		md.spcvelocity(pos,5)=-1.27;
-	end
-
-	savemodel(org,md);
-end %}}}
-if perform(org,'DiagFS'),% {{{1
-
-	md=loadmodel(org,'SetupFS');
-
-	md.eps_res=10^-5;
-	md=setcluster(md,cluster);
-	md.verbose=verbose('convergence',true,'solution',true);
-	md=solve(md,'analysis_type',DiagnosticSolutionEnum);
-
-	savemodel(org,md);
-end %}}}
-
-if perform(org,'ExtractResults'),% {{{1
-
-	%There are 250 points
-	A=expread('./Exp/ISMIP5000.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_5000.DiagHO');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx5HO=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-	A=expread('./Exp/ISMIP5000_large.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_5000.DiagFS');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx5FS=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-
-	A=expread('./Exp/ISMIP10000.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_10000.DiagHO');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx10HO=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-	A=expread('./Exp/ISMIP10000_large.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_10000.DiagFS');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx10FS=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-
-	A=expread('./Exp/ISMIP20000.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_20000.DiagHO');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx20HO=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-	A=expread('./Exp/ISMIP20000_large.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_20000.DiagFS');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx20FS=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-
-	A=expread('./Exp/ISMIP40000.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_40000.DiagHO');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx40HO=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-	A=expread('./Exp/ISMIP40000_large.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_40000.DiagFS');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx40FS=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-
-	A=expread('./Exp/ISMIP80000.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_80000.DiagHO');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx80HO=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-	A=expread('./Exp/ISMIP80000_large.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_80000.DiagFS');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx80FS=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-
-	A=expread('./Exp/ISMIP160000.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_160000.DiagHO');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx160HO=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-	A=expread('./Exp/ISMIP160000_large.exp');X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPC_160000.DiagFS');
-	vx=project2d(md,PatchToVec(md.results.DiagnosticSolution.Vx),md.numlayers);
-	md.x=md.x2d; md.mesh.y=md.mesh.y2d; md.elements=md.elements2d; md.dim=2;
-	vx160FS=InterpFromMeshToMesh2d(md.elements,md.x,md.mesh.y,vx,X,Y);
-
-	S=linspace(0,1,250)';
-	save ./Models/ISMIPC_results S vx5HO vx5FS vx10HO vx10FS vx20HO vx20FS vx40HO vx40FS vx80HO vx80FS vx160HO vx160FS;
-	clear vx* vy*
-
-end %}}}
-if perform(org,'Plot'),% {{{1
-
-	%Pattyn's standard plot
-	disp('Load and plot Pattyn''s results');
-	addpath([issmtier() '/test/IsmipHomResults/plotting_scripts/']);
-	plot_exp_c;
-	set(gcf,'Position',[500 500 1000 600],'color','w');
-	drawnow;
-
-	%We add our own contribution
-	disp('Load and plot ISSM''s results');
-	load('./Models/ISMIPC_results');
-	fsFillColor = [0.254 0.625 0.660];
-	nfsFillColor = [0.422 0.688 0.367];
-	subplot(2,3,1); set(gca,'Layer','top');
-	plot(S,vx5FS,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	plot(S,vx5HO,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	subplot(2,3,2); set(gca,'Layer','top');
-	plot(S,vx10FS,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	plot(S,vx10HO,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	subplot(2,3,3); set(gca,'Layer','top');
-	plot(S,vx20FS,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	plot(S,vx20HO,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	subplot(2,3,4); set(gca,'Layer','top');
-	plot(S,vx40FS,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	plot(S,vx40HO,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	subplot(2,3,5); set(gca,'Layer','top');
-	plot(S,vx80FS,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	plot(S,vx80HO,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	subplot(2,3,6); set(gca,'Layer','top');
-	p1=plot(S,vx160FS,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	p2=plot(S,vx160HO,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	legend([fh1 fh2 p1 p2],{'FS', 'NFS','FS ISSM','NFS ISSM'},'Location','West');
-
-	printingflag=true;
-	if printingflag, 
-		addpath([issmtier() '/externalpackages/export_fig']);
-		export_fig('../Article/figures/IsmipC.pdf');
-	end
-end %}}}
Index: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPC_jgr.par
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPC_jgr.par	(revision 14124)
+++ 	(revision )
@@ -1,25 +1,0 @@
-%Ok, start defining model parameters here
-
-disp('      creating thickness');
-md.surface=-md.x*tan(0.1*pi/180);
-md.bed=md.surface-1000;
-md.thickness=md.surface-md.bed;
-md.firn_layer=0*ones(md.numberofgrids,1);
-
-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/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))));
-%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=zeros(md.numberofelements,1);
-
-disp('      creating flow law paramter');
-md.rheology_B=6.8067*10^7*ones(md.numberofgrids,1);
-md.rheology_n=3*ones(md.numberofelements,1);
-
-disp('      boundary conditions for diagnostic model: ');
-%Create grid on boundary fist (because wi can not use mesh)
-md=SetIceSheetBC(md);
Index: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPF.m
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPF.m	(revision 14124)
+++ 	(revision )
@@ -1,41 +1,0 @@
-%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,'./ISMIPF.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: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPF.par
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPF.par	(revision 14124)
+++ 	(revision )
@@ -1,24 +1,0 @@
-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);
Index: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPF_jgr.m
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPF_jgr.m	(revision 14124)
+++ 	(revision )
@@ -1,309 +1,0 @@
-steps=[4];
-
-%Hard coded parameters
-clustername='pfe';
-interactive=1;
-printingflag=false;
-%Cluster parameters{{{1
-if strcmpi(clustername,'pfe');
-	cluster=pfe;
-	cluster.numnodes=30;
-	cluster.cpuspernode=6;
-	cluster.time=12*60;
-	cluster.interactive=interactive;
-	cluster.port=1024;
-	cluster.processor='wes';
-elseif strcmpi(clustername,'cosmos')
-	cluster=cosmos('np',100,'queue','shortq','time',1*60,'login','morlighe');
-elseif strcmpi(clustername,oshostname()),
-	cluster=eval([ oshostname() '(''np'',8);']);
-else
-	error('cluster not supported yet');
-end%}}}
-
-org=organizer('repository','./Models','prefix','ISMIPF.','steps',steps);
-if perform(org,'Pattyn_frozen'),% {{{1
-
-L=100000; %in m
-nx=50; %numberof nodes in x direction
-ny=50;
-md=model;
-md=squaremesh(md,L,L,nx,ny);
-%md=meshconvert(md);
-md=setmask(md,'',''); %ice sheet test
-md=parameterize(md,'./Par/ISMIPF.par');
-md=extrude(md,20,1);
-
-md=setflowequation(md,'pattyn','all');
-%md=setflowequation(md,'pattyn','all','stokes','all'); 
-
-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;
-md.thermal.spctemperature(:)=255;
-
-%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];
-
-md.timestepping.time_step=3;
-md.timestepping.final_time=3*800;
-md.settings.output_frequency=50;
-md.prognostic.stabilization=1;
-md.diagnostic.maxiter=1;
-md.settings.results_on_vertices=1;
-md.transient.isthermal=0;
-md.transient.requested_outputs=IceVolumeEnum;
-
-md.initialization.vx(:)=100;
-md.initialization.vy(:)=100;
-%Compute the diagnostic
-md=setcluster(md,cluster);
-md.verbose=verbose('convergence',true,'solution',true);
-md=solve(md,TransientSolutionEnum);
-plotmodel(md,'data',md.results.TransientSolution(end).Vx,'layer',md.mesh.numberoflayers,'sectionvalue',[issmtier() '/test/Exp/ISMIP100000.exp'],'figure',1)
-plotmodel(md,'data',md.results.TransientSolution(end).Surface-md.geometry.surface,'layer',md.mesh.numberoflayers,'sectionvalue',[issmtier() '/test/Exp/ISMIP100000.exp'],'figure',2);
-
-savemodel(org,md);
-end
-%}}}
-if perform(org,'Pattyn_sliding'),% {{{1
-L=100000; %in m
-nx=60; %numberof nodes in x direction
-ny=60;
-md=model;
-md=squaremesh(md,L,L,nx,ny);
-%md=meshconvert(md);
-md=setmask(md,'',''); %ice sheet test
-md=parameterize(md,'./Par/ISMIPF.par');
-md=extrude(md,10,1);
-
-md=setflowequation(md,'pattyn','all');
-%md=setflowequation(md,'pattyn','all','stokes','all'); 
-
-md.spcvelocity(:,:)=0;
-%Create dirichlet on the bed if no slip
-pos=find(md.gridonbed & (md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
-md.spcvelocity(pos,1:3)=1;
-md.spcvelocity(pos,4)=100; %because we need a dirichlet somewhere
-md.spcvelocity(pos,6)=0; %because we need a dirichlet somewhere
-md.spctemperature(:,2)=255;
-md.spctemperature(:,1)=1;
-
-%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.penalties=[posx,posx2;posy,posy2];
-
-md.dt=1;
-md.ndt=1*600;
-md.output_frequency=200;
-md.prognostic_DG=0;
-md.max_nonlinear_iterations=1;
-
-md.initialization.vx(:)=100;
-md.initialization.vy(:)=100;
-%Compute the diagnostic
-md=setcluster(md,cluster);
-md.verbose=verbose('convergence',true,'solution',true);
-md=solve(md,'analysis_type',Transient3DSolutionEnum);
-plotmodel(md,'data',PatchToVec(md.results.Transient3DSolution(end).Vx),'layer',md.numlayers,'sectionvalue',[issmtier() '/test/Exp/ISMIP100000.exp'],'figure',1)
-plotmodel(md,'data',PatchToVec(md.results.Transient3DSolution(end).Surface)-md.geometry.surface,'layer',md.numlayers,'sectionvalue',[issmtier() '/test/Exp/ISMIP100000.exp'],'figure',2);
-
-savemodel(org,md);
-end
-%}}}
-if perform(org,'Stokes_frozen'),% {{{1
-L=100000; %in m
-nx=60; %numberof nodes in x direction
-ny=60;
-md=model;
-md=squaremesh(md,L,L,nx,ny);
-%md=meshconvert(md);
-md=setmask(md,'',''); %ice sheet test
-md=parameterize(md,'./Par/ISMIPF.par');
-md=extrude(md,20,1);
-
-md=setflowequation(md,'stokes','all');
-%md=setflowequation(md,'pattyn','all','stokes','all'); 
-
-md.spcvelocity(:,:)=0;
-%Create dirichlet on the bed if no slip
-pos=find(md.gridonbed);
-md.spcvelocity(pos,1:3)=1;
-md.spctemperature(:,2)=255;
-md.spctemperature(:,1)=1;
-
-%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.penalties=[posx,posx2;posy,posy2];
-
-md.dt=2;
-md.ndt=2*800;
-md.output_frequency=50;
-md.prognostic_DG=0;
-md.max_nonlinear_iterations=1;
-
-md.initialization.vx(:)=100;
-md.initialization.vy(:)=100;
-%Compute the diagnostic
-md=setcluster(md,cluster);
-md.verbose=verbose('convergence',true,'solution',true);
-md=solve(md,'analysis_type',Transient3DSolutionEnum);
-plotmodel(md,'data',PatchToVec(md.results.Transient3DSolution(end).Vx),'layer',md.numlayers,'sectionvalue',[issmtier() '/test/Exp/ISMIP100000.exp'],'figure',1,'grid','on')
-plotmodel(md,'data',PatchToVec(md.results.Transient3DSolution(end).Surface)-md.geometry.surface,'layer',md.numlayers,'sectionvalue',[issmtier() '/test/Exp/ISMIP100000.exp'],'figure',2,'grid','on');
-
-savemodel(org,md);
-end
-%}}}
-if perform(org,'Stokes_sliding'),% {{{1
-L=100000; %in m
-nx=60; %numberof nodes in x direction
-ny=60;
-md=model;
-md=squaremesh(md,L,L,nx,ny);
-%md=meshconvert(md);
-md=setmask(md,'',''); %ice sheet test
-md=parameterize(md,'./Par/ISMIPF.par');
-md=extrude(md,4,1);
-
-md=setflowequation(md,'stokes','all');
-
-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.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
-md.diagnostic.spcvx(pos)=100; %because we need a dirichlet somewhere
-md.diagnostic.spcvy(pos)=0;
-md.diagnostic.spcvz(pos)=0;
-md.thermal.spctemperature(:)=255;
-
-%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];
-
-md.timestepping.time_step=1;
-md.timestepping.final_time=1*800;
-md.settings.output_frequency=100;
-md.prognostic.stabilization=1;
-md.diagnostic.maxiter=1;
-md.settings.results_on_vertices=1;
-md.transient.isthermal=0;
-md.transient.requested_outputs=IceVolumeEnum;
-
-md.initialization.vx(:)=100;
-md.initialization.vy(:)=100;
-%Compute the diagnostic
-md=setcluster(md,cluster);
-md.settings.results_on_vertices=1;
-md.transient.isthermal=0;
-md.verbose=verbose('convergence',true,'solution',true);
-md.diagnostic.stokes_basal_penalty=0;
-md=solve(md,TransientSolutionEnum);
-plotmodel(md,'data',(md.results.TransientSolution(end).Vx),'layer',md.numlayers,'sectionvalue',[issmtier() '/test/Exp/ISMIP100000.exp'],'figure',1)
-plotmodel(md,'data',(md.results.TransientSolution(end).Surface)-md.geometry.surface,'layer',md.numlayers,'sectionvalue',[issmtier() '/test/Exp/ISMIP100000.exp'],'figure',2);
-
-savemodel(org,md);
-end
-%}}}
-
-if perform(org,'ExtractResults'),% {{{1
-
-	%There are 250 points
-	A=expread([issmtier() '/test/Exp/ISMIP100000.exp']);X=linspace(A.x(1),A.x(2),250)';Y=linspace(A.y(1),A.y(2),250)';
-	md=loadmodel('./Models/ISMIPF.Pattyn_frozen');
-	vx=project2d(md,PatchToVec(md.results.Transient3DSolution(end).Vx),md.numlayers);
-	sf=project2d(md,PatchToVec(md.results.Transient3DSolution(end).Surface)-md.geometry.surface,md.numlayers);
-	md.mesh.x=md.mesh.x2d; md.mesh.y=md.mesh.y2d; md.mesh.elements=md.mesh.elements2d; md.dim=2;
-	vxHOf=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,vx,X,Y);
-	sfHOf=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,sf,X,Y);
-
-	md=loadmodel('./Models/ISMIPF.Pattyn_sliding');
-	vx=project2d(md,PatchToVec(md.results.Transient3DSolution(end).Vx),md.numlayers);
-	sf=project2d(md,PatchToVec(md.results.Transient3DSolution(end).Surface)-md.geometry.surface,md.numlayers);
-	md.mesh.x=md.mesh.x2d; md.mesh.y=md.mesh.y2d; md.mesh.elements=md.mesh.elements2d; md.dim=2;
-	vxHOs=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,vx,X,Y);
-	sfHOs=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,sf,X,Y);
-
-	md=loadmodel('./Models/ISMIPF.Stokes_frozen');
-	vx=project2d(md,PatchToVec(md.results.Transient3DSolution(end).Vx),md.numlayers);
-	sf=project2d(md,PatchToVec(md.results.Transient3DSolution(end).Surface)-md.geometry.surface,md.numlayers);
-	md.mesh.x=md.mesh.x2d; md.mesh.y=md.mesh.y2d; md.mesh.elements=md.mesh.elements2d; md.dim=2;
-	vxFSf=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,vx,X,Y);
-	sfFSf=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,sf,X,Y);
-
-	md=loadmodel('./Models/ISMIPF.Stokes_sliding');
-	vx=project2d(md,(md.results.TransientSolution(end).Vx),md.mesh.numberoflayers);
-	sf=project2d(md,(md.results.TransientSolution(end).Surface)-md.geometry.surface,md.mesh.numberoflayers);
-	md.mesh.x=md.mesh.x2d; md.mesh.y=md.mesh.y2d; md.mesh.elements=md.mesh.elements2d; md.mesh.dimension=2;
-	vxFSs=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,vx,X,Y);
-	sfFSs=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,sf,X,Y);
-
-	S=linspace(-50,50,250)';
-	save ./Models/ISMIPF_results S vxHOf sfHOf vxHOs sfHOs vxFSf sfFSf vxFSs sfFSs;
-	clear vx* sf*
-
-end %}}}
-if perform(org,'Plot'),% {{{1
-
-	%Pattyn's standard plot
-	disp('Load and plot Pattyn''s results');
-	addpath([issmtier() '/test/IsmipHomResults/plotting_scripts/']);
-	plot_exp_f;
-	plot_exp_fs;
-	set(gcf,'Position',[500 500 700 700],'color','w');
-	drawnow;
-
-	%We add our own contribution
-	disp('Load and plot ISSM''s results');
-	load('./Models/ISMIPF_results');
-	fsFillColor = [0.254 0.625 0.660];
-	nfsFillColor = [0.422 0.688 0.367];
-	disp('display ISSM''s results');
-	subplot(2,2,1); set(gca,'Layer','top');
-	plot(S,vxFSf,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	plot(S,vxHOf,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	subplot(2,2,3); set(gca,'Layer','top');
-	plot(S,vxFSs,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	plot(S,vxHOs,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	subplot(2,2,2); set(gca,'Layer','top');
-	plot(S,sfFSf,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	plot(S,sfHOf,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-	subplot(2,2,4);set(gca,'Layer','top');
-	p1=plot(S,sfFSs,'Color',fsFillColor,'LineWidth',2,'LineStyle','-');
-	p2=plot(S,sfHOs,'Color',nfsFillColor,'LineWidth',2,'LineStyle','-');
-
-	legend([fh1 fh2 p1 p2],{'FS', 'NFS','FS ISSM','NFS ISSM'},2);
-
-	if true, 
-		addpath([issmtier() '/externalpackages/export_fig']);
-		export_fig('../Article/figures/IsmipF.pdf');
-	end
-end %}}}
Index: sm/workshop/2012/Talks/11_ISMIP/Files/ISMIPF_jgr.par
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/ISMIPF_jgr.par	(revision 14124)
+++ 	(revision )
@@ -1,41 +1,0 @@
-%Ok, start defining model parameters here
-disp('      creating geometry');
-alpha=3*pi/180;
-md.geometry.surface=-md.mesh.x*tan(alpha);
-md.geometry.bed=md.geometry.surface-1000+100*exp(-((md.mesh.x-max(md.mesh.x)/2).^2+(md.mesh.y-max(md.mesh.y)/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);
-md.materials.rheology_law='None';
-
-disp('      boundary conditions for diagnostic model');
-%Create grid on boundary fist (because we cannot use mesh)
-md=SetIceSheetBC(md);
-%md.spcvelocity(:,4)=100*ones(md.mesh.numberofvertices,1);
-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.vel=zeros(md.mesh.numberofvertices,1);
-md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
-md.initialization.temperature=255*ones(md.mesh.numberofvertices,1);
-pos=find(md.mesh.x==min(md.mesh.x) | md.mesh.x==max(md.mesh.x) | md.mesh.y==min(md.mesh.y) | md.mesh.y==max(md.mesh.y));
-%md.balancethickness.spcthickness(pos,1)=1;
-%md.balancethickness.spcthickness(pos,2)=md.geometry.thickness(pos);
-md.thermal.spctemperature=[ones(md.mesh.numberofvertices,1) 255*ones(md.mesh.numberofvertices,1)];
-md.basalforcings.geothermalflux=0.4*ones(md.mesh.numberofvertices,1);
-
-%Parallel options
-md.mesh.average_vertex_connectivity=200;
-
-%Transient options
-md.timestepping.time_step=1;
-md.timestepping.final_time=10;
-md.prognostic.stabilization=1;
-md.transient.isthermal=0;
Index: sm/workshop/2012/Talks/11_ISMIP/Files/Profile.exp
===================================================================
--- /issm/workshop/2012/Talks/11_ISMIP/Files/Profile.exp	(revision 14124)
+++ 	(revision )
@@ -1,7 +1,0 @@
-## Name:icefront
-## Icon:0
-# Points Count  Value
-2 1.
-# X pos Y pos
-0 50000
-100000 50000