source: issm/branches/trunk-jpl-damage/src/m/boundaryconditions/SetMarineIceSheetBC.m@ 13101

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 ARGUS contour around the ice front
%   can be given in input, or it will be deduced as onfloatingice & 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

%node on Dirichlet (boundary and ~icefront)
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
        nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2);
        vertexonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront);
else
        %Guess where the ice front is
        vertexonfloatingice=zeros(md.mesh.numberofvertices,1);
        vertexonfloatingice(md.mesh.elements(find(md.mask.elementonfloatingice),:))=1;
        vertexonicefront=double(md.mesh.vertexonboundary & vertexonfloatingice);
end
pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
if isempty(pos),
        warning('SetMarineIceSheetBC warning: ice front all around the glacier, no dirichlet found. Dirichlet must be added manually')
end
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

md.hydrology.spcwatercolumn=zeros(md.mesh.numberofvertices,2);
pos=find(md.mesh.vertexonboundary); 
md.hydrology.spcwatercolumn(pos,1)=1;

%segment on Neumann (Ice Front)
pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(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))+ 0*md.mask.elementongroundedice(pressureload(:,end))];

%plug onto model
md.diagnostic.icefront=pressureload;


%Create zeros basalforcings and surfaceforcings
if (isnan(md.surfaceforcings.precipitation) & (md.surfaceforcings.ispdd==1)),
        md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
        disp('      no surfaceforcings.precipitation specified: values set as zero');
end
if (isnan(md.surfaceforcings.mass_balance) & (md.surfaceforcings.ispdd==0)),
        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);
                md.basalforcings.geothermalflux(find(md.mask.vertexongroundedice))=50*10^-3; %50mW/m2
        end
else
        disp('      no thermal boundary conditions created: no observed temperature found');
end