source: issm/oecreview/Archive/25834-26739/ISSM-26257-26258.diff@ 26740

Last change on this file since 26740 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 13.0 KB
RevLine 
[26740]1Index: ../trunk-jpl/examples/Functions/grace.m
2===================================================================
3--- ../trunk-jpl/examples/Functions/grace.m (revision 26257)
4+++ ../trunk-jpl/examples/Functions/grace.m (revision 26258)
5@@ -1,8 +1,11 @@
6 % map grace loads in meters [m] of water equivalent height
7-function water_load = grace(index,lat,lon,tmin,tmax);
8+function water_load = grace(index,lat,lon,tmin,tmax,onvertex);
9
10 %compute centroids using (lat,lon) data {{{
11 ne = length(index); % number of elements
12+ nv = length(lat); % number of vertices
13+ % [lat,lon] \in [-90:90,-180,180];
14+ lat_vertex_0=lat; long_vertex_0=lon;
15 % lat -> [0,180]; long -> [0,360] to compute centroids
16 lat=90-lat;
17 lon(lon<0)=180+(180+lon(lon<0));
18@@ -115,10 +118,18 @@
19 time_yr=time_yr(pos1:pos2);
20 load_grace_plus=load_grace_plus(pos1:pos2,:);
21 num_yr=length(time_yr);
22- water_load_0=zeros(ne,num_yr);
23+ if onvertex
24+ water_load_0=zeros(nv,num_yr);
25+ else
26+ water_load_0=zeros(ne,num_yr);
27+ end
28
29 for jj=1:num_yr
30- water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element);
31+ if onvertex
32+ water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_vertex_0,lon);
33+ else
34+ water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element);
35+ end
36 disp([num2str(jj),' of ',num2str(num_yr),' months done!']);
37 end
38
39Index: ../trunk-jpl/examples/SlrFarrell/runme.m
40===================================================================
41--- ../trunk-jpl/examples/SlrFarrell/runme.m (revision 26257)
42+++ ../trunk-jpl/examples/SlrFarrell/runme.m (revision 26258)
43@@ -1,10 +1,7 @@
44 clear all;
45
46-steps=[1:5];
47+steps=[4];
48
49-try
50- % [1:5]
51-
52 if any(steps==1)
53 disp(' Step 1: Global mesh creation');
54
55@@ -19,16 +16,16 @@
56 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km]
57
58 for i=1:numrefine,
59- md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
60+ ocean_mask_levelset=gmtmask(md.mesh.lat,md.mesh.long);
61
62 distance=zeros(md.mesh.numberofvertices,1);
63
64- pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos);
65- pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos);
66+ pos=find(~ocean_mask_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos);
67+ pos=find(ocean_mask_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos);
68
69 for j=1:md.mesh.numberofvertices
70 phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
71- if md.mask.ocean_levelset(j),
72+ if ocean_mask_levelset(j),
73 phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi;
74 deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
75 d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
76@@ -41,7 +38,7 @@
77 end
78 pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
79
80- pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
81+ pos2=find(ocean_mask_levelset~=1 & distance>mindistance_land);
82 distance(pos2)=mindistance_land;
83
84 dist=min(maxdistance,distance);
85@@ -48,7 +45,10 @@
86 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist);
87 end
88
89- md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
90+ ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
91+ pos = find(ocean_mask==0);
92+ md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
93+ md.mask.ocean_levelset(pos)=1;
94
95 save ./Models/SlrFarrell_Mesh md;
96
97@@ -59,13 +59,19 @@
98 disp(' Step 2: Define source as in Farrell, 1972, Figure 1');
99 md = loadmodel('./Models/SlrFarrell_Mesh');
100
101- md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere
102- md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
103- md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
104+ md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofvertices,1);
105+ md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
106
107+ pos = find(md.mask.ocean_levelset==-1);
108+ md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
109+ md.solidearth.initialsealevel(pos)=1; % 1 m SLR everywhere
110+ md.dsl.global_average_thermosteric_sea_level_change=[0;0];
111+ md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
112+ md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
113+
114 save ./Models/SlrFarrell_Loads md;
115
116- plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]');
117+ plotmodel (md,'data',md.solidearth.initialsealevel,'view',[90 90],'title#all','Initial sea-level [m]');
118 end
119
120 if any(steps==3)
121@@ -74,18 +80,17 @@
122
123 md.solidearth.lovenumbers=lovenumbers('maxdeg',10000);
124
125- md.mask.land_levelset = 1-md.mask.ocean_levelset;
126- md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
127- md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1);
128- pos=find(md.mesh.lat <-80);
129- md.mask.ice_levelset(pos(1))=-1; % ice yes!
130- md.mask.ocean_levelset(pos(1))=1; % ice grounded!
131+ %md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
132+ %md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1);
133+ %pos=find(md.mesh.lat <-80);
134+ %md.mask.ice_levelset(pos(1))=-1; % ice yes!
135+ %md.mask.ocean_levelset(pos(1))=1; % ice grounded!
136
137- di=md.materials.rho_ice/md.materials.rho_water;
138- md.geometry.thickness=ones(md.mesh.numberofvertices,1);
139- md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
140- md.geometry.base=md.geometry.surface-md.geometry.thickness;
141- md.geometry.bed=md.geometry.base;
142+ % arbitary to pass consistency check.
143+ md.geometry.bed=-ones(md.mesh.numberofvertices,1);
144+ md.geometry.surface=ones(md.mesh.numberofvertices,1);
145+ md.geometry.base=md.geometry.bed;
146+ md.geometry.thickness=md.geometry.surface-md.geometry.base;
147
148 md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
149 md.materials.rheology_B=paterson(md.initialization.temperature);
150@@ -103,7 +108,7 @@
151 md.cluster=generic('name',oshostname(),'np',3);
152 md.verbose=verbose('111111111');
153
154- md.slr.reltol = 0.1/100; % percent change in solution
155+ md.solidearth.settings.reltol = 0.1/100; % percent change in solution
156
157 md=solve(md,'Slr');
158
159Index: ../trunk-jpl/examples/EsaGRACE/runme.m
160===================================================================
161--- ../trunk-jpl/examples/EsaGRACE/runme.m (revision 26257)
162+++ ../trunk-jpl/examples/EsaGRACE/runme.m (revision 26258)
163@@ -1,7 +1,7 @@
164 clear all;
165 addpath('../Data','../Functions');
166
167-steps=[1]; % [1:5]
168+steps=[5];
169
170 if any(steps==1)
171 disp(' Step 1: Global mesh creation');
172Index: ../trunk-jpl/examples/EsaWahr/runme.m
173===================================================================
174--- ../trunk-jpl/examples/EsaWahr/runme.m (revision 26257)
175+++ ../trunk-jpl/examples/EsaWahr/runme.m (revision 26258)
176@@ -1,7 +1,7 @@
177 clear all;
178 addpath('../Functions');
179
180-steps=[1]; % [1:7]
181+steps=[2];
182
183 if any(steps==1)
184 disp(' Step 1: Mesh creation');
185Index: ../trunk-jpl/examples/SlrGRACE/runme.m
186===================================================================
187--- ../trunk-jpl/examples/SlrGRACE/runme.m (revision 26257)
188+++ ../trunk-jpl/examples/SlrGRACE/runme.m (revision 26258)
189@@ -1,9 +1,9 @@
190 clear all;
191 addpath('../Data','../Functions');
192
193-steps=[7];
194+steps=[1:7];
195
196-if any(steps==1)
197+if any(steps==1) % {{{
198 disp(' Step 1: Global mesh creation');
199
200 numrefine=5;
201@@ -50,8 +50,8 @@
202 save ./Models/SlrGRACE_Mesh md;
203
204 plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
205-end
206-if any(steps==2)
207+end % }}}
208+if any(steps==2) % {{{
209 disp(' Step 2: Define loads in meters of ice height equivalent');
210 md = loadmodel('./Models/SlrGRACE_Mesh');
211
212@@ -58,40 +58,59 @@
213 year_month = 2007+15/365;
214 time_range = [year_month year_month];
215
216- water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
217+ onvertex = 1; % map data on vertex. If 0, it maps on the elemental centroid.
218+ water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2),onvertex);
219+ rho_water2ice = md.materials.rho_freshwater/md.materials.rho_ice;
220+ ice_load = water_load*rho_water2ice; % ice height equivalent.
221
222- md.solidearth.surfaceload.icethicknesschange = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
223+ %Geometry for the bed, arbitrary thickness of 100:
224+ md.geometry.bed=zeros(md.mesh.numberofvertices,1);
225+ md.geometry.base=md.geometry.bed;
226+ md.geometry.thickness = 100*ones(md.mesh.numberofvertices,1);
227+ md.geometry.surface=md.geometry.bed+md.geometry.thickness;
228
229+ md.masstransport.spcthickness = [md.geometry.thickness + ice_load; 0];
230+
231+ md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
232+
233 save ./Models/SlrGRACE_Loads md;
234
235- plotmodel (md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]');
236-end
237-if any(steps==3)
238+ plotmodel (md,'data',ice_load,...
239+ 'edgecolor','k','view',[45 45],'caxis',[-.1 .1],...
240+ 'title','Ice height equivalent [m]');
241+end % }}}
242+if any(steps==3) % {{{
243 disp(' Step 3: Parameterization');
244 md = loadmodel('./Models/SlrGRACE_Loads');
245-
246+
247+ md.mask.ice_levelset=-md.mask.ocean_levelset;
248+
249 md.solidearth.lovenumbers = lovenumbers('maxdeg',10000);
250+ md.solidearth.settings.reltol=NaN;
251+ md.solidearth.settings.abstol=1e-3;
252+ md.solidearth.settings.computesealevelchange=1;
253+ md.solidearth.settings.isgrd=1;
254+ md.solidearth.settings.grdmodel=1;
255+ md.solidearth.settings.maxiter=10;
256+
257+ %time stepping:
258+ md.timestepping.start_time=0;
259+ md.timestepping.time_step=1;
260+ md.timestepping.final_time=1;
261
262- md.mask.ice_levelset=-md.mask.ocean_levelset;
263- md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
264- md.dsl.global_average_thermosteric_sea_level_change=[0;0];
265- md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
266- md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
267+ %masstransport:
268+ md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
269+ md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
270+ md.initialization.vx=zeros(md.mesh.numberofvertices,1);
271+ md.initialization.vy=zeros(md.mesh.numberofvertices,1);
272+ md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
273+ md.initialization.str=0;
274
275- md.solidearth.settings.ocean_area_scaling=1;
276-
277- % arbitary to pass consistency check. md.geometry.bed=-ones(md.mesh.numberofvertices,1);
278- md.geometry.surface=ones(md.mesh.numberofvertices,1);
279- md.geometry.base=md.geometry.bed; md.geometry.thickness=md.geometry.surface-md.geometry.base;
280- md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
281- md.materials.rheology_B=paterson(md.initialization.temperature);
282- md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
283-
284 md.miscellaneous.name='SlrGRACE';
285
286 save ./Models/SlrGRACE_Parameterization md;
287-end
288-if any(steps==4)
289+end % }}}
290+if any(steps==4) % {{{
291 disp(' Step 4: Solve Slr solver');
292 md = loadmodel('./Models/SlrGRACE_Parameterization');
293
294@@ -99,21 +118,28 @@
295 md.solidearth.settings.elastic=1;
296 md.solidearth.settings.rotation=1;
297
298- md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
299+ %Physics:
300+ md.transient.issmb=0;
301+ md.transient.isstressbalance=0;
302+ md.transient.isthermal=0;
303+ md.transient.ismasstransport=1;
304+ md.transient.isslc=1;
305+
306+ md.solidearth.requested_outputs={'Sealevel','Bed'};
307
308 md.cluster=generic('name',oshostname(),'np',3);
309 md.verbose=verbose('111111111');
310
311- md=solve(md,'Slr');
312+ md=solve(md,'Transient');
313
314 save ./Models/SlrGRACE_Solution md;
315-end
316-if any(steps==5)
317+end % }}}
318+if any(steps==5) % {{{
319 disp(' Step 5: Plot solutions');
320 md = loadmodel('./Models/SlrGRACE_Solution');
321
322- sol1 = md.solidearth.surfaceload.icethicknesschange*100; % [cm]
323- sol2 = md.results.SealevelriseSolution.SealevelRSL*1000; % [mm]
324+ sol1 = (md.masstransport.spcthickness(1:end-1)-md.geometry.thickness)*100; % [cm]
325+ sol2 = (md.results.TransientSolution.Sealevel-md.results.TransientSolution.Bed)*1000; % [mm]
326
327 sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
328 fig_name={'Fig_dH.pdf','Fig_slr.pdf'};
329@@ -153,7 +179,7 @@
330 if (kk==1)
331 geoshow(flipud(coastlat),flipud(coastlon),'DisplayType','polygon','FaceColor','white');
332 else
333- geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor',[.5 1 .5]);
334+ geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor','none');
335 end
336 plot(coastlon, coastlat,'k'); hold off;
337 c1=colorbar;
338@@ -166,8 +192,8 @@
339 set(gcf,'color','w');
340 %export_fig(fig_name{kk});
341 end
342-end
343-if any(steps==6)
344+end % }}}
345+if any(steps==6) % {{{
346 disp(' Step 6: Transient run');
347 md = loadmodel('./Models/SlrGRACE_Parameterization');
348
349@@ -200,8 +226,8 @@
350 md=solve(md,'tr');
351
352 save ./Models/SlrGRACE_Transient md;
353-end
354-if any(steps==7)
355+end % }}}
356+if any(steps==7) % {{{
357 disp(' Step 7: Plot transient');
358 md = loadmodel('./Models/SlrGRACE_Transient');
359
360@@ -283,4 +309,5 @@
361 ylabel('GMSL [mm/yr]');
362 set(gcf,'color','w');
363 %export_fig('Fig7.pdf');
364-end
365+end % }}}
366+
Note: See TracBrowser for help on using the repository browser.