clear all; steps=[5]; % if any(steps==1) % Global mesh creation {{{ disp(' Step 1: Global mesh creation'); numrefine=1; resolution=150*1e3; % inital resolution [m]. It determines, e.g., whether we capture small islands. radius = 6.371012*10^6; % mean radius of Earth, m mindistance_coast=150*1e3; % coastal resolution [m] mindistance_land=300*1e3; % resolution on the continents [m] maxdistance=600*1e3; % max element size (on mid-oceans) [m] %mesh earth: md=model; md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. for i=1:numrefine, %figure out mask: md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); %figure out distance to the coastline, in lat,long (not x,y,z): distance=zeros(md.mesh.numberofvertices,1); pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); for j=1:md.mesh.numberofvertices %figure out nearest coastline (using the great circle distance) phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; if md.mask.ocean_levelset(j), phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); else phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); end distance(j)=min(d); end pos=find(distancemindistance_land); distance(pos2)=mindistance_land; dist=min(maxdistance,distance); % max size 1000 km %use distance to the coastline to refine mesh: md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); end %figure out mask: md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); save ./Models/SlrFarrell.Mesh md; plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); %export_fig('Fig1.pdf'); end % }}} if any(steps==2) % Define source {{{ disp(' Step 2: Define source as in Farrell, 1972, Figure 1'); md = loadmodel('./Models/SlrFarrell.Mesh'); % initial sea-level: 1 m RSL everywhere. md.slr.sealevel=md.mask.ocean_levelset; md.slr.deltathickness=zeros(md.mesh.numberofelements,1); md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); save ./Models/SlrFarrell.Loads md; plotmodel (md,'data',md.slr.sealevel,'view',[90 90],... 'title#all','Initial sea-level [m]'); %export_fig('Fig2.pdf'); end % }}} if any(steps==3) % Parameterization {{{ disp(' Step 3: Parameterization'); md = loadmodel('./Models/SlrFarrell.Loads'); % Love numbers and reference frame: CF or CM (choose one!) nlove=10001; % up to 10,000 degree md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[]; md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[]; md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[]; % Mask: for computational efficiency only those elements that have loads are convolved! md.mask.land_levelset = 1-md.mask.ocean_levelset; % fake ice load in one element! md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); % no ice md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1); % floated... pos=find(md.mesh.lat <-80); md.mask.ice_levelset(pos(1))=-1; % ice yes! md.mask.groundedice_levelset(pos(1))=1; % ice grounded! %% IGNORE BUT DO NOT DELETE %% {{{ % Geometry: Important only when you want to couple with Ice Flow Model di=md.materials.rho_ice/md.materials.rho_water; md.geometry.thickness=ones(md.mesh.numberofvertices,1); md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); md.geometry.base=md.geometry.surface-md.geometry.thickness; md.geometry.bed=md.geometry.base; % Materials: md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); md.materials.rheology_B=paterson(md.initialization.temperature); md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); % Miscellaneous: md.miscellaneous.name='SlrFarrell'; %% IGNORE BUT DO NOT DELETE %% }}} save ./Models/SlrFarrell.Parameterization md; end % }}} if any(steps==4) % Solve {{{ disp(' Step 4: Solve Slr solver'); md = loadmodel('./Models/SlrFarrell.Parameterization'); % Cluster info md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); % Solve md=solve(md,'Slr'); save ./Models/SlrFarrell.Solution md; end % }}} if any(steps==5) % Plot solutions {{{ disp(' Step 5: Plot solutions'); md = loadmodel('./Models/SlrFarrell.Solution'); % solutions. sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) res = 1.0; % degree % Make a grid of lats and lons, based on the min and max of the original vectors [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); sol_grid = zeros(size(lat_grid)); % Make a interpolation object F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); F.Method = 'linear'; F.ExtrapolationMethod = 'linear'; % Do the interpolation to get gridded solutions... sol_grid = F(lat_grid, lon_grid); sol_grid(isnan(sol_grid))=0; sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) figure1=figure('Position', [100, 100, 1000, 500]); gcf; load coast; cla; %pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; contourf(lon_grid,lat_grid,sol_grid,[96 98 100 102 104 105]); shading flat; hold on; geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); plot(long,lat,'k'); hold off; % define colormap, caxis, xlim etc {{{ c1=colorbar; %colormap('haxby'); %caxis([96 104]); xlim([-180 180]); ylim([-90 90]); % }}} grid on; title('Relative sea-level [mm]'); set(gcf,'color','w'); %export_fig('Fig5.pdf'); end % }}}