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

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

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