1 | def SetIceShelfBC(md,icefrontfile='')
|
---|
2 | #SETICESHELFBC - Create the boundary conditions for diagnostic and thermal models for a Ice Shelf with Ice Front
|
---|
3 | #
|
---|
4 | # Neumann BC are used on the ice front (an ANRGUS contour around the ice front
|
---|
5 | # must be given in input)
|
---|
6 | # Dirichlet BC are used elsewhere for diagnostic
|
---|
7 | #
|
---|
8 | # Usage:
|
---|
9 | # md=SetIceShelfBC(md,varargin)
|
---|
10 | #
|
---|
11 | # Example:
|
---|
12 | # md=SetIceShelfBC(md);
|
---|
13 | # md=SetIceShelfBC(md,'Front.exp');
|
---|
14 | #
|
---|
15 | # See also: SETICESHEETBC, SETMARINEICESHEETBC
|
---|
16 |
|
---|
17 | #node on Dirichlet (boundary and ~icefront)
|
---|
18 | if not icefrontfile:
|
---|
19 | nodeonicefront=zeros(md.mesh.numberofvertices)
|
---|
20 | else:
|
---|
21 | if not os.path.isfile(icefrontfile):
|
---|
22 | print 'SetIceShelfBC error message: ice front file '+icefrontfile+ ' not found'
|
---|
23 | return []
|
---|
24 | nodeinsideicefront=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,icefrontfile,'node',2);
|
---|
25 | nodeonicefront=double(md.mesh.vertexonboundary & nodeinsideicefront);
|
---|
26 |
|
---|
27 | pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
|
---|
28 | md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1);
|
---|
29 | md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1);
|
---|
30 | md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1);
|
---|
31 | md.diagnostic.spcvx(pos)=0;
|
---|
32 | md.diagnostic.spcvy(pos)=0;
|
---|
33 | md.diagnostic.spcvz(pos)=0;
|
---|
34 | md.diagnostic.referential=NaN*ones(md.mesh.numberofvertices,6);
|
---|
35 |
|
---|
36 | #Dirichlet Values
|
---|
37 | if (length(md.inversion.vx_obs)==md.mesh.numberofvertices & length(md.inversion.vy_obs)==md.mesh.numberofvertices)
|
---|
38 | disp(' boundary conditions for diagnostic model: spc set as observed velocities');
|
---|
39 | md.diagnostic.spcvx(pos)=md.inversion.vx_obs(pos);
|
---|
40 | md.diagnostic.spcvy(pos)=md.inversion.vy_obs(pos);
|
---|
41 | else
|
---|
42 | disp(' boundary conditions for diagnostic model: spc set as zero');
|
---|
43 | end
|
---|
44 |
|
---|
45 | #segment on Ice Front
|
---|
46 | #segment on Neumann (Ice Front)
|
---|
47 | pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
|
---|
48 | if (md.mesh.dimension==2)
|
---|
49 | pressureload=md.mesh.segments(pos,:);
|
---|
50 | elseif md.mesh.dimension==3
|
---|
51 | 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)];
|
---|
52 | pressureload=[];
|
---|
53 | for i=1:md.mesh.numberoflayers-1,
|
---|
54 | pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
|
---|
55 | end
|
---|
56 | end
|
---|
57 |
|
---|
58 | #Add water or air enum depending on the element
|
---|
59 | pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
|
---|
60 |
|
---|
61 | #plug onto model
|
---|
62 | md.diagnostic.icefront=pressureload;
|
---|
63 |
|
---|
64 | #Create zeros basalforcings and surfaceforcings
|
---|
65 | if isnan(md.surfaceforcings.precipitation),
|
---|
66 | md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
|
---|
67 | disp(' no surfaceforcings.precipitation specified: values set as zero');
|
---|
68 | end
|
---|
69 | if isnan(md.surfaceforcings.mass_balance),
|
---|
70 | md.surfaceforcings.mass_balance=zeros(md.mesh.numberofvertices,1);
|
---|
71 | disp(' no surfaceforcings.mass_balance specified: values set as zero');
|
---|
72 | end
|
---|
73 | if isnan(md.basalforcings.melting_rate),
|
---|
74 | md.basalforcings.melting_rate=zeros(md.mesh.numberofvertices,1);
|
---|
75 | disp(' no basalforcings.melting_rate specified: values set as zero');
|
---|
76 | end
|
---|
77 | if isnan(md.balancethickness.thickening_rate),
|
---|
78 | md.balancethickness.thickening_rate=zeros(md.mesh.numberofvertices,1);
|
---|
79 | disp(' no balancethickness.thickening_rate specified: values set as zero');
|
---|
80 | end
|
---|
81 |
|
---|
82 | md.prognostic.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
|
---|
83 | md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
|
---|
84 |
|
---|
85 | if (length(md.initialization.temperature)==md.mesh.numberofvertices),
|
---|
86 | md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1);
|
---|
87 | pos=find(md.mesh.vertexonsurface); md.thermal.spctemperature(pos)=md.initialization.temperature(pos); #impose observed temperature on surface
|
---|
88 | if (length(md.basalforcings.geothermalflux)~=md.mesh.numberofvertices),
|
---|
89 | md.basalforcings.geothermalflux=zeros(md.mesh.numberofvertices,1);
|
---|
90 | end
|
---|
91 | else
|
---|
92 | disp(' no thermal boundary conditions created: no observed temperature found');
|
---|
93 | end
|
---|