source: issm/trunk-jpl/src/m/utils/BC/SetIceShelfBC.py@ 12127

def SetIceShelfBC(md,icefrontfile='')
        #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 not icefrontfile:
                nodeonicefront=zeros(md.mesh.numberofvertices)
        else:
                if not os.path.isfile(icefrontfile):
                        print 'SetIceShelfBC error message: ice front file '+icefrontfile+ ' not found'
                        return []
                nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2);
                nodeonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront);
        
        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 and surfaceforcings
        if isnan(md.surfaceforcings.precipitation),
                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.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