source: issm/trunk-jpl/examples/SlrFarrell/runme.m@ 22818

Last change on this file since 22818 was 22818, checked in by adhikari, 7 years ago

CHG: several minor changes
File size: 5.1 KB
Line 
1
2clear all;
3
4steps=[1]; % [1:5]
5
6if any(steps==1)
7 disp(' Step 1: Global mesh creation');
8
9 numrefine=1;
10 resolution=150*1e3; % inital resolution [m]
11 radius = 6.371012*10^6; % mean radius of Earth, m
12 mindistance_coast=150*1e3; % coastal resolution [m]
13 mindistance_land=300*1e3; % resolution on the continents [m]
14 maxdistance=600*1e3; % max element size (on mid-oceans) [m]
15
16 md=model;
17 md.mask=maskpsl();
18 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km]
19
20 for i=1:numrefine,
21
22 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
23
24 distance=zeros(md.mesh.numberofvertices,1);
25
26 pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos);
27 pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos);
28
29 for j=1:md.mesh.numberofvertices
30 phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
31 if md.mask.ocean_levelset(j),
32 phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi;
33 deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
34 d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
35 else
36 phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi;
37 deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
38 d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
39 end
40 distance(j)=min(d);
41 end
42 pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
43
44 pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
45 distance(pos2)=mindistance_land;
46
47 dist=min(maxdistance,distance);
48 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist);
49
50 end
51
52 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
53
54 save ./Models/SlrFarrell_Mesh md;
55
56 plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
57
58end
59
60if any(steps==2)
61 disp(' Step 2: Define source as in Farrell, 1972, Figure 1');
62 md = loadmodel('./Models/SlrFarrell_Mesh');
63
64 md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere
65
66 md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
67 md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
68
69 save ./Models/SlrFarrell_Loads md;
70
71 plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]');
72
73end
74
75if any(steps==3)
76 disp(' Step 3: Parameterization');
77 md = loadmodel('./Models/SlrFarrell_Loads');
78
79 nlove=10001;
80 md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[];
81 md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[];
82 md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[];
83
84 md.mask.land_levelset = 1-md.mask.ocean_levelset;
85 md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
86 md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1);
87 pos=find(md.mesh.lat <-80);
88 md.mask.ice_levelset(pos(1))=-1; % ice yes!
89 md.mask.groundedice_levelset(pos(1))=1; % ice grounded!
90
91 di=md.materials.rho_ice/md.materials.rho_water;
92 md.geometry.thickness=ones(md.mesh.numberofvertices,1);
93 md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
94 md.geometry.base=md.geometry.surface-md.geometry.thickness;
95 md.geometry.bed=md.geometry.base;
96
97 md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
98 md.materials.rheology_B=paterson(md.initialization.temperature);
99 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
100
101 md.miscellaneous.name='SlrFarrell';
102
103 save ./Models/SlrFarrell_Parameterization md;
104
105end
106
107if any(steps==4)
108 disp(' Step 4: Solve Slr solver');
109 md = loadmodel('./Models/SlrFarrell_Parameterization');
110
111 md.cluster=generic('name',oshostname(),'np',3);
112 md.verbose=verbose('111111111');
113
114 md.slr.reltol = 0.1/100; % per cent change in solution
115
116 md=solve(md,'Slr');
117
118 save ./Models/SlrFarrell_Solution md;
119
120end
121
122if any(steps==5)
123 disp(' Step 5: Plot solutions');
124 md = loadmodel('./Models/SlrFarrell_Solution');
125
126 sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m)
127
128 res = 1; % [degree]
129
130 [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
131 sol_grid = zeros(size(lat_grid));
132
133 F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol);
134 F.Method = 'linear';
135 F.ExtrapolationMethod = 'linear';
136
137 sol_grid = F(lat_grid, lon_grid);
138 sol_grid(isnan(sol_grid))=0;
139 sol_grid(lat_grid>85 & sol_grid==0) =NaN;
140
141 set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
142 figure1=figure('Position', [100, 100, 1000, 500]);
143 gcf; load coast; cla;
144 pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
145 [C,h]=contour(lon_grid,lat_grid,sol_grid,[96 98 100 102 104 105],'-k','linewidth',2);
146 clabel(C,h,'FontSize',18,'Color','red','LabelSpacing',500);
147 geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]);
148 plot(long,lat,'k'); hold off;
149 c1=colorbar;
150 colormap(flipud(haxby));
151 caxis([96 105]);
152 xlim([-170 170]);
153 ylim([-85 85]);
154 grid on;
155 title('Relative sea-level [% of GMSL]');
156 set(gcf,'color','w');
157 %export_fig('Fig5.pdf');
158
159end
160
Note: See TracBrowser for help on using the repository browser.