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

Last change on this file since 12127 was 12127, checked in by Eric.Larour, 13 years ago

Starting conversion of SetIceShelfBC
Needs ContourToMesh python module, so converted this module to python interface.
Also simplified Contours* treatment.

File size: 4.0 KB
Line 
1def 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
Note: See TracBrowser for help on using the repository browser.