function md=geography(md,iceshelfname,icesheetname)
%GEOGRAPHY - establish boundaries between grounded and floating ice.
%
%   By default, ice is considered grounded. The contour iceshelfname defines grids 
%   for which ice is floting. The contour icesheetname defines grids inside an iceshelf, 
%   that are grounded (ie: ice rises, islands, etc ...)
%   All input files are in the Argus format (extension .exp).
%
%   Usage:
%      md=geography(md,iceshelfname,icesheetname)
%
%   Examples:
%      md=geography(md,'all','');
%      md=geography(md,'Iceshelves.exp','Islands.exp');

%some checks on list of arguments
if ((nargin~=3) | (nargout~=1)),
	geographyusage();
	error('geography error message');
end
if ((~ischar(iceshelfname)) | (~ischar(icesheetname))  )
	geographyusage();
	error('geography error message');
end

if md.counter>=2,
	choice=input('This model already has a geometry. Are you sure you want to go ahead? (y/n)','s');
	if ~strcmp(choice,'y')
		error('no geometry done ... exiting');
	end
else
	if (md.counter~=1)
		error('geography error message: you need to run mesh.m first on this model');
	else
		md.counter=2;
	end
end

%Get assigned fields
x=md.x;
y=md.y;
elements=md.elements;

%Fist recover arrays of ice shelf grids and elements, and ice sheet grids and elements.
if strcmp(iceshelfname,''), %no iceshelf contour file, we are dealing with a pure ice sheet.
	gridoniceshelf=zeros(md.numberofgrids,1);
	elementoniceshelf=zeros(md.numberofelements,1);
elseif strcmp(iceshelfname,'all'), %we are dealing with a pure ice shelf.
	gridoniceshelf=ones(md.numberofgrids,1);
	elementoniceshelf=ones(md.numberofelements,1);
else
	[gridoniceshelf,elementoniceshelf]=ArgusContourToMesh(elements,x,y,expread(iceshelfname,1),'element and node',2);
end

if strcmp(icesheetname,''), %no icesheet contour file, we are dealing with a pure ice shelf.
	gridonicesheet=zeros(md.numberofgrids,1);
	elementonicesheet=zeros(md.numberofelements,1);
else
	[gridonicesheet,elementonicesheet]=ArgusContourToMesh(elements,x,y,expread(icesheetname,1),'element and node',2);
end

%Because icesheet grids and elements can be included into an iceshelf, we need to update. Remember, all the previous 
%arrays come from domain outlines that can intersect one another: 
gridoniceshelf=double((gridoniceshelf & ~gridonicesheet));
elementoniceshelf=double((elementoniceshelf & ~elementonicesheet));
gridonicesheet=double(~gridoniceshelf);
elementonicesheet=double(~elementoniceshelf);

%Return: 
md.elementoniceshelf=elementoniceshelf;
md.gridoniceshelf=gridoniceshelf;

md.elementonicesheet=elementonicesheet;
md.gridonicesheet=gridonicesheet;

%Keep track of input files
if strcmp(iceshelfname,''), md.iceshelfoutline=''; elseif strcmp(iceshelfname,'all'), md.iceshelfoutline='all'; else md.iceshelfoutline=readfile(iceshelfname); end
if strcmp(icesheetname,''), md.icesheetoutline=''; elseif strcmp(icesheetname,'all'), md.icesheetoutline='all'; else md.icesheetoutline=readfile(icesheetname); end

function geographyusage(),
disp('INPUT md=geography(md,iceshelfname,icesheetname)');
disp('geography:  establish boundaries between grounded and floating ice.');
disp('     By default, ice is considered grounded. The contour iceshelfname defines grids ');
disp('     for which ice is floting. The contour icesheetname defines grids inside an iceshelf, ');
disp('     that are grounded (ie: ice rises, islands, etc ...)');
disp('     All input files are in the Argus format (extension .exp');
