Index: sm/workshop/2012/Talks/07_Ice_flow_models/Code/SetIceShelfBC.m
===================================================================
--- /issm/workshop/2012/Talks/07_Ice_flow_models/Code/SetIceShelfBC.m	(revision 14127)
+++ 	(revision )
@@ -1,100 +1,0 @@
-function md=SetIceShelfBC(md,varargin)
-%SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a  Ice Shelf with Ice Front
-%
-%   Neumann BC are used on the ice front (an ANRGUS contour around the ice front
-%   must be given in input)
-%   Dirichlet BC are used elsewhere for diagnostic
-%
-%   Usage:
-%      md=SetIceShelfBC(md,varargin)
-%
-%   Example:
-%      md=SetIceShelfBC(md);
-%      md=SetIceShelfBC(md,'Front.exp');
-%
-%   See also: SETICESHEETBC, SETMARINEICESHEETBC
-
-%node on Dirichlet (boundary and ~icefront)
-if nargin==2,
-	icefrontfile=varargin{1};
-	if ~exist(icefrontfile), error(['SetIceShelfBC error message: ice front file ' icefrontfile ' not found']); end
-	nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2);
-	nodeonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront);
-elseif nargin==1,
-	nodeonicefront=zeros(md.mesh.numberofvertices,1);
-else
-	help SetIceShelfBC
-	error('bad usage');
-end
-pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
-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);
-md.diagnostic.spcvx(pos)=0;
-md.diagnostic.spcvy(pos)=0;
-md.diagnostic.spcvz(pos)=0;
-md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);
-
-%Dirichlet Values
-if (length(md.inversion.vx_obs)==md.mesh.numberofvertices & length(md.inversion.vy_obs)==md.mesh.numberofvertices)
-	disp('      boundary conditions for diagnostic model: spc set as observed velocities');
-	md.diagnostic.spcvx(pos)=md.inversion.vx_obs(pos);
-	md.diagnostic.spcvy(pos)=md.inversion.vy_obs(pos);
-else
-	disp('      boundary conditions for diagnostic model: spc set as zero');
-end
-
-%segment on Ice Front
-%segment on Neumann (Ice Front)
-pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
-if (md.mesh.dimension==2)
-	pressureload=md.mesh.segments(pos,:);
-elseif md.mesh.dimension==3
-	pressureload_layer1=[md.mesh.segments(pos,1:2)  md.mesh.segments(pos,2)+md.mesh.numberofvertices2d  md.mesh.segments(pos,1)+md.mesh.numberofvertices2d  md.mesh.segments(pos,3)];
-	pressureload=[];
-	for i=1:md.mesh.numberoflayers-1,
-		pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
-	end
-end
-
-%Add water or air enum depending on the element
-pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
-
-%plug onto model
-md.diagnostic.icefront=pressureload;
-
-%Create zeros basalforcings.melting_rate, surfaceforcings.ablation_rate, surfaceforcings.accumulation_rate
-% and surfaceforcings.mass_balance if not specified
-if isnan(md.surfaceforcings.accumulation_rate),
-	md.surfaceforcings.accumulation_rate=zeros(md.mesh.numberofvertices,1);
-	disp('      no surfaceforcings.accumulation_rate specified: values set as zero');
-end
-if isnan(md.surfaceforcings.ablation_rate),
-	md.surfaceforcings.ablation_rate=zeros(md.mesh.numberofvertices,1);
-	disp('      no surfaceforcings.ablation_rate specified: values set as zero');
-end
-if isnan(md.surfaceforcings.mass_balance),
-	md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1);
-	disp('      no surfaceforcings.mass_balance specified: values set as zero');
-end
-if isnan(md.basalforcings.melting_rate),
-	md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
-	disp('      no basalforcings.melting_rate specified: values set as zero');
-end
-if isnan(md.balancethickness.thickening_rate),
-	md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1);
-	disp('      no balancethickness.thickening_rate specified: values set as zero');
-end
-
-md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
-md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
-
-if (length(md.initialization.temperature)==md.mesh.numberofvertices),
-	md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
-	pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); %impose observed temperature on surface
-	if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
-		md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1);
-	end
-else
-	disp('      no thermal boundary conditions created: no observed temperature found');
-end