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

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

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