Changeset 421

Timestamp:

05/14/09 09:53:05 (16 years ago)

Author:

Mathieu Morlighem

Message:

SetMarineIceSheetBC can be done without Front.exp file

File:

• issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (modified) (1 diff)

issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m

@@ -1 @@
-function md=SetMarineIceSheetBC(md,icefrontfile)
+function md=SetMarineIceSheetBC(md,varargin)
 %SETICEMARINESHEETBC - Create the boundary conditions for diagnostic and thermal models for a  Marine Ice Sheet with Ice Front
 %
-%   Neumann BC are used on the ice front (an ANRGUS contour around the ice front
-%   must be given in input)
+%   Neumann BC are used on the ice front (an ARGUS contour around the ice front
+%   can be given in input, or it will be deduced as oniceshelf & onboundary)
 %   Dirichlet BC are used elsewhere for diagnostic
 %
 %   Usage:
 %      md=SetMarineIceSheetBC(md,icefrontfile)
+%      md=SetMarineIceSheetBC(md)
 %
 %   Example:
 %      md=SetMarineIceSheetBC(md,'Front.exp')
+%      md=SetMarineIceSheetBC(md)
 %
 %   See also: SETICESHELFBC, SETMARINEICESHEETBC
 
 %grid on Dirichlet (boundary and ~icefront)
-if ~exist(icefrontfile)
-        error(['SetMarineIceSheetBC error message: ice front file ' icefrontfile ' not found']);
+if nargin==2,
+        %User provided Front.exp, use it
+        icefrontfile=varargin{1};
+        if ~exist(icefrontfile)
+                error(['SetMarineIceSheetBC error message: ice front file ' icefrontfile ' not found']);
+        end
+        gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread(icefrontfile,1),'node',2);
+        gridonicefront=double(md.gridonboundary & gridinsideicefront);
+
+else
+        %Guess where the ice front is
+        pos=find(md.gridonboundary & md.gridoniceshelf);
+        md.gridondirichlet_diag=zeros(md.numberofgrids,1);
+        md.gridondirichlet_diag(pos)=1;
 end
-gridinsideicefront=ContourToMesh(md.elements,md.x,md.y,expread(icefrontfile,1),'node',2);
-gridonicefront=double(md.gridonboundary & gridinsideicefront);
 pos=find(md.gridonboundary & ~gridonicefront);
+if isempty(pos),
+        error('SetMarineIceSheetBC error message: ice front all around the glacier, no dirichlet found. Set BC manually')
+end
 md.gridondirichlet_diag=zeros(md.numberofgrids,1);
 md.gridondirichlet_diag(pos)=1;