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