source: issm/trunk-jpl/test/SandBox/test2004.m@ 24809

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

CHG: converging towards a final version.
File size: 13.6 KB
Line 
1%Test Name: Earth_Antarctica_GIA
2
3testagainst2002=1;
4
5%Data paths: {{{
6modeldatapath='/Users/larour/ModelData';
7shppath='/Users/larour/issm-jpl/proj-group/qgis/NightlyRun';
8gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1-NightlyRun.shp'];
9%}}}
10
11%create sealevel model to hold our information:
12 sl=sealevelmodel();
13
14%Create basins using boundaries from shapefile:
15%some projections we'll rely on: %{{{
16proj4326=epsg2proj(4326);
17proj3031=epsg2proj(3031);
18%}}}
19%HemisphereWest: {{{
20sl.addbasin(basin('continent','hemispherewest','name','hemispherewest','proj',laea(0,-90),'boundaries',{... %Peru projection 3587
21 boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326,'orientation','reverse'),...
22 boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),...
23 boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031,'orientation','reverse'),...
24 boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),...
25 boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse'),...
26 boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),...
27 boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031,'orientation','reverse'),...
28 boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031)...
29 }));
30 %}}}
31 %HemisphereEast: {{{
32 sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long
33 boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),...
34 boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
35 boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),...
36 boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)...
37 }));
38%}}}
39%Antarctica excluding Ronne: {{{
40sl.addbasin(basin('continent','antarctica','name','antarctica-grounded','proj',proj3031,'boundaries',{...
41 boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),...
42 boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031),...
43 boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
44 boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031)...
45 boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031)...
46 boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031)...
47 boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031)...
48 boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031)...
49 }));
50 %}}}
51 %Ronne: {{{
52 sl.addbasin(basin('continent','antarctica','name','ronne','proj',proj3031,'boundaries',{...
53 boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),...
54 boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),...
55 boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),...
56 boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse')...
57 }));
58%}}}
59
60%Meshing
61%Go through basins and mesh: %{{{
62%meshing parameters: {{{
63hmin=2000; hmax=10000; hmin=hmin*1000; hmax=hmax*1000;
64tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes
65threshold=5;
66defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinment',1};
67alreadyloaded=0;
68%}}}
69for i=sl.basinindx('basin','all'),
70
71 bas=sl.basins{i};
72 disp(sprintf('Meshing basin %s\n',bas.name));
73
74 %recover basin domain:
75 domain=bas.contour();
76
77 %recover coastline inside basin, using GSHHS_c_L1. It's a lat,long file, hence epsg 4326
78 coastline=bas.shapefilecrop('shapefile',gshhsshapefile,'epsgshapefile',4326,'threshold',threshold);
79
80 %mesh:
81 md=bamg(model,'domain',domain,'subdomains',coastline,'hmin',hmin,'hmax',hmax,defaultoptions{:});
82 plotmodel(md,'data','mesh');pause(1);
83
84 %miscellaneous:
85 md.mesh.proj=bas.proj; md.miscellaneous.name=bas.name;
86
87 %recover mask where we have land:
88 md.private.bamg.landmask=double(md.private.bamg.mesh.Triangles(:,4)>=1);
89
90 %vertex connectivity:
91 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
92
93 %add model to sl icecaps:
94 sl.addicecap(md);
95end
96%}}}
97
98%Parameterization:
99%Parameterize ice sheets : {{{
100
101for ind=sl.basinindx('continent',{'antarctica'}),
102 disp(sprintf('Parameterizing basin %s\n', sl.icecaps{ind}.miscellaneous.name));
103
104 md=sl.icecaps{ind}; bas=sl.basins{ind};
105 %masks : %{{{
106 %new type of mask:
107 md.mask=maskpsl;
108
109 %land and ocean:
110 md.mask.land_levelset=ones(md.mesh.numberofvertices,1);
111 md.mask.ocean_levelset=zeros(md.mesh.numberofvertices,1);
112
113 %ice levelset from domain outlines:
114 md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
115 pos=find(md.mesh.vertexonboundary); md.mask.ice_levelset(pos)=0;
116
117 %ice front:
118 pos=find(md.mask.ice_levelset==0); md.mask.ocean_levelset(pos)=-1; md.mask.land_levelset(pos)=1;
119
120 %now, glaciers:
121 md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
122 %}}}
123 %initial grounding line position: {{{
124 if bas.isnameany('antarctica-grounded'),
125
126 md.mask.groundedice_levelset=ones(md.mesh.numberofvertices,1);
127 end
128 if bas.isnameany('ronne'),
129
130 md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
131
132 %correction to land and ocean levelset: ice shelf is not on land!
133 pos=find(md.mask.ice_levelset<=0 & md.mask.groundedice_levelset<=0);
134 md.mask.ocean_levelset(pos)=1;
135 md.mask.land_levelset(pos)=-1;
136 end
137 %}}}
138 %latlong: % {{{
139 [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326');
140 %}}}
141 %geometry: {{{
142 if bas.iscontinentany('antarctica'),
143 di=md.materials.rho_ice/md.materials.rho_water;
144
145 disp(' reading bedrock');
146 md.geometry.bed=interpBedmap2(md.mesh.x,md.mesh.y,'bed');
147 end % }}}
148 %Slr: {{{
149 if bas.iscontinentany('antarctica'),
150 if testagainst2002,
151 md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
152 %antarctica
153 late=sum(md.mesh.lat(md.mesh.elements),2)/3;
154 longe=sum(md.mesh.long(md.mesh.elements),2)/3;
155 pos=find(late <-80);
156 md.slr.deltathickness(pos)=-100;
157 else
158 delH=textread([modeldatapath '/AdhikariSciAdvTrends/AIS_delH_trend.txt']);
159 longAIS=delH(:,1);
160 latAIS=delH(:,2);
161 delHAIS=delH(:,3);
162 index=delaunay(latAIS,longAIS);
163
164 lat=md.mesh.lat;
165 long=md.mesh.long+360;
166 pos=find(long>360);long(pos)=long(pos)-360;
167
168 delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat);
169 northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0;
170
171 md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100;
172 end
173
174 md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
175 md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
176 md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
177 md.slr.Ugia=zeros(md.mesh.numberofvertices,1);
178
179 md.dsl.global_average_thermosteric_sea_level_change=[0;0];
180 md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
181 md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
182 md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1);
183
184 end %}}}
185 % material properties: {{{
186 md.materials=materials('hydro');
187 %}}}
188 %diverse: {{{
189 md.miscellaneous.name=bas.name;
190 % }}}
191
192 sl.icecaps{ind}=md;
193end
194%}}}
195% ParameterizeContinents {{{
196
197sl.basinindx('continent',{'hemisphereeast','hemispherewest'})
198
199for ind=sl.basinindx('continent',{'hemisphereeast','hemispherewest'}),
200 disp(sprintf('Masks for basin %s\n', sl.icecaps{ind}.miscellaneous.name));
201 md=sl.icecaps{ind}; bas=sl.basins{ind};
202
203 %recover lat,long:
204 [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326');
205
206 %mask: %{{{
207 %new type of mask:
208 md.mask=maskpsl;
209 %Figure out mask from initial mesh: deal with land and ocean masks (land include grounded ice). %{{{
210 %first, transform land element mask into vertex driven one.
211 land=md.private.bamg.landmask;
212 land_mask=-ones(md.mesh.numberofvertices,1);
213
214 landels=find(land);
215 land_mask(md.mesh.elements(landels,:))=1;
216
217 %gothrough edges of each land element
218 connectedels=md.mesh.elementconnectivity(landels,:);
219 connectedisonocean=~land(connectedels);
220 sumconnectedisonocean=sum(connectedisonocean,2);
221
222 %figure out which land elements are connected to the ocean:
223 landelsconocean=landels(find(sumconnectedisonocean));
224
225 ind1=[md.mesh.elements(landelsconocean,1); md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3)];
226 ind2=[md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3); md.mesh.elements(landelsconocean,1)];
227
228 h=waitbar(0,'Starting mask processing');
229
230 %edge ind1 and ind2:
231 for i=1:length(ind1),
232 if ~mod(i,100),
233 waitbar(i/length(ind1),h,sprintf('%.2i%% along...',i/length(ind1)*100));
234 end
235
236 els1=md.mesh.vertexconnectivity(ind1(i),1: md.mesh.vertexconnectivity(ind1(i),end));
237 els2=md.mesh.vertexconnectivity(ind2(i),1: md.mesh.vertexconnectivity(ind2(i),end));
238 els=intersect(els1,els2);
239
240 if length(find(land(els)))==1,
241 %this edge is on the beach, 0 the edge:
242 land_mask(ind1(i))=0;
243 land_mask(ind2(i))=0;
244 end
245 end
246 md.mask.land_levelset=land_mask;
247 close(h);
248 %}}}
249 %work on ocean, glaciers and ice: {{{
250 %ocean: opposite of land:
251 md.mask.ocean_levelset=-md.mask.land_levelset;
252
253 %initliaze:
254 md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
255
256 %grounded ice:
257 md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);
258
259 md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
260
261 % }}}
262 %}}}
263 %slr loading/calibration: {{{
264
265 if testagainst2002,
266 md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
267 %greenland
268 late=sum(md.mesh.lat(md.mesh.elements),2)/3;
269 longe=sum(md.mesh.long(md.mesh.elements),2)/3;
270 pos=find(late > 70 & late < 80 & longe>-60 & longe<-30);
271 md.slr.deltathickness(pos)=-100;
272
273 %correct mask:
274 md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
275 else
276 md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
277 end
278
279 md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
280 md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
281 md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
282 md.slr.Ugia=zeros(md.mesh.numberofvertices,1);
283
284 md.dsl.global_average_thermosteric_sea_level_change=[0;0];
285 %md.slr.steric_rate=(1.1+.38)*ones(md.mesh.numberofvertices,1); %steric + water storage.
286 md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
287 md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
288 md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1);
289
290 %love numbers:
291 md.slr.ocean_area_scaling=1;
292 %}}}
293 %geometry: {{{
294 di=md.materials.rho_ice/md.materials.rho_water;
295 md.geometry.bed=-ones(md.mesh.numberofvertices,1);
296 % }}}
297 %materials: {{{
298 md.materials=materials('hydro');
299 % }}}
300 sl.icecaps{ind}=md;
301end
302% }}}
303
304%%Assemble Earth in 3D {{{
305
306%parameters:
307plotting=1;
308tolerance=100;
309loneedgesdetect=0;
310
311%create earth model by concatenating all the icecaps in 3d:
312sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect);
313
314%figure out how each icecap's mesh connects to the larger earth mesh:
315sl.intersections('force',1);
316
317%figure out connectivity:
318disp('Mesh connectivity');
319sl.earth.mesh.vertexconnectivity=NodeConnectivity(sl.earth.mesh.elements,sl.earth.mesh.numberofvertices);
320
321%areas:
322disp('Mesh nodal areas');
323sl.earth.mesh.area=averaging(sl.earth,GetAreas3DTria(sl.earth.mesh.elements,sl.earth.mesh.x,sl.earth.mesh.y,sl.earth.mesh.z),4);
324
325%transfer a list of fields from each icecap and continent back to Earth:
326sl.transfer('mask.groundedice_levelset');
327sl.transfer('mask.ice_levelset');
328sl.transfer('mask.ocean_levelset');
329sl.transfer('mask.land_levelset');
330sl.transfer('mask.glacier_levelset');
331sl.transfer('geometry.bed');
332sl.transfer('mesh.lat');
333sl.transfer('mesh.long');
334sl.transfer('slr.deltathickness');
335sl.transfer('slr.spcthickness');
336sl.transfer('slr.Ngia');
337sl.transfer('slr.Ugia');
338sl.transfer('slr.hydro_rate');
339sl.transfer('slr.sealevel');
340sl.transfer('dsl.sea_surface_height_change_above_geoid');
341sl.transfer('dsl.sea_water_pressure_change_at_sea_floor');
342
343%radius:
344sl.earth.mesh.r=sqrt(sl.earth.mesh.x.^2+sl.earth.mesh.y.^2+sl.earth.mesh.z.^2);
345
346%check on the mesh transitions: {{{
347plotting=1;
348if plotting,
349 flags=ones(sl.earth.mesh.numberofelements,1);
350 for i=1:length(sl.eltransitions),
351 flags(sl.eltransitions{i})=i;
352 end
353 plotmodel(sl.earth,'data',flags,'shading','faceted','coastline','on','coast_color','g')
354end
355%}}}}
356
357% }}}
358
359%Solve Sea-level equation on Earth only: {{{
360md=sl.earth; %we don't do computations on ice sheets or land.
361
362%Materials:
363md.materials=materials('hydro');
364
365%elastic loading from love numbers:
366nlov=1001;
367md.slr.love_h = love_numbers('h'); md.slr.love_h(nlov+1:end)=[];
368md.slr.love_k = love_numbers('k'); md.slr.love_k(nlov+1:end)=[];
369md.slr.love_l = love_numbers('l'); md.slr.love_l(nlov+1:end)=[];
370md.slr.ocean_area_scaling = 0;
371
372%Miscellaneous
373md.miscellaneous.name='test2004';
374
375%New stuff
376md.dsl.global_average_thermosteric_sea_level_change=[0;0];
377
378%Solution parameters
379md.slr.reltol=NaN;
380md.slr.abstol=1e-3;
381md.slr.geodetic=1;
382
383% max number of iteration reverted back to 10 (i.e., the original default value)
384md.slr.maxiter=10;
385
386%eustatic run:
387md.slr.rigid=0; md.slr.elastic=0;md.slr.rotation=0;
388md=solve(md,'Sealevelrise');
389Seustatic=md.results.SealevelriseSolution.Sealevel;
390
391%eustatic + rigid run:
392md.slr.rigid=1; md.slr.elastic=0;md.slr.rotation=0;
393md=solve(md,'Sealevelrise');
394Srigid=md.results.SealevelriseSolution.Sealevel;
395
396%eustatic + rigid + elastic run:
397md.slr.rigid=1; md.slr.elastic=1;md.slr.rotation=0;
398md=solve(md,'Sealevelrise');
399Selastic=md.results.SealevelriseSolution.Sealevel;
400
401%eustatic + rigid + elastic + rotation run:
402md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1;
403md=solve(md,'Sealevelrise');
404Srotation=md.results.SealevelriseSolution.Sealevel;
405
406%}}}
Note: See TracBrowser for help on using the repository browser.