Index: /issm/trunk-jpl/examples/SlrGRACE/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrGRACE/runme.m (revision 26141) +++ /issm/trunk-jpl/examples/SlrGRACE/runme.m (revision 26142) @@ -4,9 +4,9 @@ steps=[7]; -if any(steps==1) +if any(steps==1) disp(' Step 1: Global mesh creation'); - numrefine=1; - resolution=300*1e3; % inital resolution [m] + numrefine=5; + resolution=150*1e3; % inital resolution [m] radius = 6.371012*10^6; % mean radius of Earth, m mindistance_coast=150*1e3; % coastal resolution [m] @@ -23,7 +23,5 @@ distance=zeros(md.mesh.numberofvertices,1); - pos=find(~ocean_mask); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); - pos=find(ocean_mask); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); - + pos=find(~ocean_mask); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); pos=find(ocean_mask); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); for j=1:md.mesh.numberofvertices phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; @@ -40,5 +38,5 @@ end pos=find(distancemindistance_land); distance(pos2)=mindistance_land; @@ -49,14 +47,10 @@ ocean_mask=gmtmask(md.mesh.lat,md.mesh.long); - pos = find(ocean_mask==0); - md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); - md.mask.ocean_levelset(pos)=1; - + pos = find(ocean_mask==0); md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); md.mask.ocean_levelset(pos)=1; save ./Models/SlrGRACE_Mesh md; plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); -end - -if any(steps==2) +end +if any(steps==2) disp(' Step 2: Define loads in meters of ice height equivalent'); md = loadmodel('./Models/SlrGRACE_Mesh'); @@ -72,7 +66,6 @@ plotmodel (md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]'); -end - -if any(steps==3) +end +if any(steps==3) disp(' Step 3: Parameterization'); md = loadmodel('./Models/SlrGRACE_Loads'); @@ -80,6 +73,5 @@ md.solidearth.lovenumbers = lovenumbers('maxdeg',10000); - md.mask.ice_levelset=-md.mask.ocean_levelset; - + md.mask.ice_levelset=-md.mask.ocean_levelset; md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); md.dsl.global_average_thermosteric_sea_level_change=[0;0]; @@ -89,10 +81,7 @@ md.solidearth.settings.ocean_area_scaling=1; - % arbitary to pass consistency check. - md.geometry.bed=-ones(md.mesh.numberofvertices,1); + % arbitary to pass consistency check. md.geometry.bed=-ones(md.mesh.numberofvertices,1); md.geometry.surface=ones(md.mesh.numberofvertices,1); - md.geometry.base=md.geometry.bed; - md.geometry.thickness=md.geometry.surface-md.geometry.base; - + md.geometry.base=md.geometry.bed; md.geometry.thickness=md.geometry.surface-md.geometry.base; md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); @@ -102,7 +91,6 @@ save ./Models/SlrGRACE_Parameterization md; -end - -if any(steps==4) +end +if any(steps==4) disp(' Step 4: Solve Slr solver'); md = loadmodel('./Models/SlrGRACE_Parameterization'); @@ -120,11 +108,10 @@ save ./Models/SlrGRACE_Solution md; -end - -if any(steps==5) +end +if any(steps==5) disp(' Step 5: Plot solutions'); md = loadmodel('./Models/SlrGRACE_Solution'); - sol1 = md.solidearth.surfaceload.icethicknesschange*100; % [cm] + sol1 = md.solidearth.surfaceload.icethicknesschange*100; % [cm] sol2 = md.results.SealevelriseSolution.SealevelRSL*1000; % [mm] @@ -155,5 +142,5 @@ F.Method = 'linear'; F.ExtrapolationMethod = 'linear'; - + sol_grid = F(lat_grid, lon_grid); sol_grid(isnan(sol_grid))=0; @@ -162,12 +149,12 @@ set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) figure1=figure('Position', [100, 100, 1000, 500]); - gcf; load coast; cla; + gcf; load coastlines; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; if (kk==1) - geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + geoshow(flipud(coastlat),flipud(coastlon),'DisplayType','polygon','FaceColor','white'); else - geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); - end - plot(long,lat,'k'); hold off; + geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor',[.5 1 .5]); + end + plot(coastlon, coastlat,'k'); hold off; c1=colorbar; colormap('haxby'); @@ -180,7 +167,6 @@ %export_fig(fig_name{kk}); end -end - -if any(steps==6) +end +if any(steps==6) disp(' Step 6: Transient run'); md = loadmodel('./Models/SlrGRACE_Parameterization'); @@ -199,6 +185,5 @@ md.transient.isstressbalance=0; md.transient.isthermal=0; - md.transient.isgia=1; - md.transient.isslr=1; + md.transient.isgia=1; md.transient.isslr=1; md.timestepping.start_time=-1; @@ -212,11 +197,10 @@ md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'}; - + md=solve(md,'tr'); save ./Models/SlrGRACE_Transient md; -end - -if any(steps==7) +end +if any(steps==7) disp(' Step 7: Plot transient'); md = loadmodel('./Models/SlrGRACE_Transient'); @@ -230,5 +214,5 @@ end sol_name = {'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; - movie_name = {'Movie_dH.avi','Movie_slr.avi'}; + movie_name = {'Movie_dH.mp4','Movie_slr.mp4'}; res = 1.0; @@ -254,11 +238,13 @@ vidObj = VideoWriter(movie_name{kk}); vidObj.FrameRate=2; % frames per second + vifObj.FileFormat='mp4'; open(vidObj); for jj=1:length(time) + fprintf('creating frame %d...', jj); F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); F.Method = 'linear'; F.ExtrapolationMethod = 'linear'; - + sol_grid = F(lat_grid, lon_grid); sol_grid(isnan(sol_grid))=0; @@ -267,12 +253,12 @@ set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) figure1=figure('Position', [100, 100, 1000, 500]); - gcf; load coast; cla; + gcf; load coastlines; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; if (kk==1) - geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + geoshow(flipud(coastlat),flipud(coastlon),'DisplayType','polygon','FaceColor','white'); else - geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); - end - plot(long,lat,'k'); hold off; + geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor','white'); + end + plot(coastlon,coastlat,'k'); hold off; c1=colorbar; colormap('haxby'); @@ -285,8 +271,10 @@ writeVideo(vidObj,getframe(gcf)); close - end + fprintf('done\n'); + end + disp('closing vidObj...'); close(vidObj); end - !open *.avi; + !open *.mp4; % plot GMSL time series @@ -296,3 +284,3 @@ set(gcf,'color','w'); %export_fig('Fig7.pdf'); -end +end