Index: /issm/trunk-jpl/examples/EsaGRACE/runme.m =================================================================== --- /issm/trunk-jpl/examples/EsaGRACE/runme.m (revision 22817) +++ /issm/trunk-jpl/examples/EsaGRACE/runme.m (revision 22818) @@ -3,83 +3,57 @@ addpath('../Data','../Functions'); -steps=[0]; % [0:5]; +steps=[1]; % [1:5] -if any(steps==0) % Download GRACE land_mass data {{{ - disp(' Step 0: Download GRACE land_mass data'); - - % Connect to ftp server and download - f = ftp('podaac-ftp.jpl.nasa.gov'); - %dir(f) - cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); - mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); - - !mv *.nc '../Data/' - - % display the content of the netcdf file. - %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') - -end % }}} - -if any(steps==1) % Global mesh creation {{{ +if any(steps==1) disp(' Step 1: Global mesh creation'); - resolution=300; % km. uniform meshing... - radius = 6.371012*10^3; % mean radius of Earth, km + resolution=300; % [km] + radius = 6.371012*10^3; % [km] - %mesh earth: md=model; - md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset + md.mask=maskpsl(); md.mesh=gmshplanet('radius',radius,'resolution',resolution); - %figure out mask: md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); - save ./Models/EsaGRACE.Mesh md; + save ./Models/EsaGRACE_Mesh md; plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); - %export_fig('Fig1.pdf'); -end % }}} +end -if any(steps==2) % Define loads {{{ +if any(steps==2) disp(' Step 2: Define loads in meters of ice height equivalent'); - md = loadmodel('./Models/EsaGRACE.Mesh'); + md = loadmodel('./Models/EsaGRACE_Mesh'); - % define time interval for analysis year_month = 2007+15/365; time_range = [year_month year_month]; - % map GRACE water load on to the mesh for the seleted month(s) water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent - save ./Models/EsaGRACE.Loads md; + save ./Models/EsaGRACE_Loads md; plotmodel (md,'data',md.esa.deltathickness,... 'view',[90 -90],'caxis',[-.1 .1],... 'title','Ice height equivalent [m]'); - %export_fig('Fig2.pdf'); -end % }}} +end -if any(steps==3) % Parameterization {{{ +if any(steps==3) disp(' Step 3: Parameterization'); - md = loadmodel('./Models/EsaGRACE.Loads'); + md = loadmodel('./Models/EsaGRACE_Loads'); - % Love numbers and reference frame: CF or CM (choose one!) - nlove=10001; % up to 10,000 degree + nlove=10001; md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[]; md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[]; - % Mask: for computational efficiency only those elements that have loads are convolved! - md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded + md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); pos=find(md.esa.deltathickness~=0); - md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; md.mask.land_levelset = 1-md.mask.ocean_levelset; - %% 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); @@ -87,43 +61,38 @@ 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='EsaGRACE'; - %% IGNORE BUT DO NOT DELETE %% }}} - save ./Models/EsaGRACE.Parameterization md; + save ./Models/EsaGRACE_Parameterization md; -end % }}} +end -if any(steps==4) % Solve {{{ +if any(steps==4) disp(' Step 4: Solve Esa solver'); - md = loadmodel('./Models/EsaGRACE.Parameterization'); + md = loadmodel('./Models/EsaGRACE_Parameterization'); - % Request outputs md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'}; - % Cluster info md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); - % Solve md=solve(md,'Esa'); - save ./Models/EsaGRACE.Solution md; + save ./Models/EsaGRACE_Solution md; -end % }}} +end -if any(steps==5) % Plot solutions {{{ +if any(steps==5) disp(' Step 5: Plot solutions'); - md = loadmodel('./Models/EsaGRACE.Solution'); + md = loadmodel('./Models/EsaGRACE_Solution'); - % loads and solutions. - sol1 = md.esa.deltathickness*100; % WEH cm - sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] - sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] - sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] + sol1 = md.esa.deltathickness*100; % [cm] + sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] + sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] + sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',... @@ -131,7 +100,6 @@ fig_name={'Fig_dH.pdf','Fig_vert.pdf','Fig_horzNS.pdf','Fig_horzEW.pdf'}; - res = 1.0; % degree + 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)); @@ -140,7 +108,5 @@ sol=eval(sprintf('sol%d',kk)); - % if data are on elements, map those on to the vertices {{{ if length(sol)==md.mesh.numberofelements - % map on to the vertices for jj=1:md.mesh.numberofelements ii=(jj-1)*3; @@ -152,15 +118,13 @@ end sol=temp'; - end % }}} + end - % 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 + sol_grid(lat_grid>85 & sol_grid==0)=NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) @@ -168,10 +132,8 @@ gcf; load coast; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; - % mask out oceans {{{ if (kk==1) geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); - end % }}} + end plot(long,lat,'k'); hold off; - % colormap, xlim etc {{{ c1=colorbar; colormap('haxby'); @@ -179,12 +141,10 @@ xlim([-180 180]); ylim([-90 90]); - % }}} grid on; title(sol_name(kk)); set(gcf,'color','w'); - %export_fig(fig_name{kk}); end -end % }}} +end Index: /issm/trunk-jpl/examples/EsaWahr/runme.m =================================================================== --- /issm/trunk-jpl/examples/EsaWahr/runme.m (revision 22817) +++ /issm/trunk-jpl/examples/EsaWahr/runme.m (revision 22818) @@ -3,81 +3,66 @@ addpath('../Functions'); -steps=[0]; % [0:6]; +steps=[0]; % [0:6] -if any(steps==0) % Simple mesh creation {{{ +if any(steps==0) disp(' Step 0: Mesh creation'); - % Generate initial uniform mesh - md=roundmesh(model,100000,10000); % Domain radius and element size (meters) + md=roundmesh(model,100000,10000); % Domain radius and element size [m] - save ./Models/EsaWahr.Mesh md; + save ./Models/EsaWahr_Mesh md; plotmodel(md,'data','mesh'); - %export_fig('Fig0.pdf'); -end % }}} +end -if any(steps==1) % {{{ Anisotropic mesh creation +if any(steps==1) disp(' Step 1: Anisotropic mesh creation'); - % Generate initial uniform mesh - md=roundmesh(model,100000,1000); % Domain radius and element size (meters) - %plotmodel(md,'data','mesh'); + md=roundmesh(model,100000,1000); - % Define a fake field to have anisotropic meshing - disc_radius=20*1000; % km => m - rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); % radial distance [m] + disc_radius=20*1000; + rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); field = abs(rad_dist-disc_radius); - %plotmodel(md,'data',field); - md = bamg(md,'field',field,'err',50,'hmax',10000); % error in field in m + md = bamg(md,'field',field,'err',50,'hmax',10000); - save ./Models/EsaWahr.Mesh md; + save ./Models/EsaWahr_Mesh md; plotmodel (md,'data','mesh'); - %export_fig('Fig1.pdf'); -end % }}} +end -if any(steps==2) % Define loads {{{ +if any(steps==2) disp(' Step 2: Define loads'); - md = loadmodel('./Models/EsaWahr.Mesh'); + md = loadmodel('./Models/EsaWahr_Mesh'); - % primer rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; - % elemental centroids index=md.mesh.elements; x_cent=mean(md.mesh.x(index),2); y_cent=mean(md.mesh.y(index),2); - % Define a 20-km radius disc and unload it by 1 meter water height equivalent uniformly md.esa.deltathickness = zeros(md.mesh.numberofelements,1); - disc_radius=20; % km - rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000; % radial distance in km - md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i; % 1 m water withdrawl + disc_radius=20; % [km] + rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000; + md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i; - save ./Models/EsaWahr.Loads md; + save ./Models/EsaWahr_Loads md; plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]'); - %export_fig('Fig2.pdf'); -end % }}} +end -if any(steps==3) % Parameterization {{{ +if any(steps==3) disp(' Step 3: Parameterization'); - md = loadmodel('./Models/EsaWahr.Loads'); + md = loadmodel('./Models/EsaWahr_Loads'); - % Love numbers and reference frame: CF or CM (choose one!) - nlove=10001; % up to 10,000 degree + nlove=10001; md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[]; md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[]; - % Mask: for computational efficiency only those elements that have loads are convolved! - md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); % -1 = ice loads - md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded + md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); + md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); - %% 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); @@ -85,45 +70,39 @@ 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='EsaWahr'; - %% IGNORE BUT DO NOT DELETE %% }}} - save ./Models/EsaWahr.Parameterization md; + save ./Models/EsaWahr_Parameterization md; -end % }}} +end -if any(steps==4) % Solve {{{ +if any(steps==4) disp(' Step 4: Solve Esa solver'); - md = loadmodel('./Models/EsaWahr.Parameterization'); + md = loadmodel('./Models/EsaWahr_Parameterization'); - % Request outputs md.esa.requested_outputs = {'EsaUmotion','EsaXmotion','EsaYmotion'}; - % Cluster info md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); - % Solve md=solve(md,'Esa'); - save ./Models/EsaWahr.Solution md; + save ./Models/EsaWahr_Solution md; -end % }}} +end -if any(steps==5) % Plot solutions {{{ +if any(steps==5) disp(' Step 5: Plot solutions'); - md = loadmodel('./Models/EsaWahr.Solution'); + md = loadmodel('./Models/EsaWahr_Solution'); - % m => mm - vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] - horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] - horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] - horz = sqrt(horz_n.^2+horz_e.^2); + vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] + horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] + horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] + horz = sqrt(horz_n.^2+horz_e.^2); % [mm] - % plot set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); figure('Position', [100, 100, 800, 600]); @@ -147,21 +126,18 @@ 'axispos',[0.505 0.02 0.47 0.47],... 'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); - %export_fig('Fig5.pdf'); -end % }}} +end -if any(steps==6) % Compare results against semi-analytic solutions {{{ +if any(steps==6) disp(' Step 6: Compare results against Wahr semi-analytic solutions'); - md = loadmodel('./Models/EsaWahr.Solution'); + md = loadmodel('./Models/EsaWahr_Solution'); - % numerical solutions - % m => mm - vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] - horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] - horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] - horz = sqrt(horz_n.^2+horz_e.^2); - % grid data from the center - xi=[0:500:100000]; % 500 m resolution + vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] + horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] + horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] + horz = sqrt(horz_n.^2+horz_e.^2); % [mm] + + xi=[0:500:100000]; % grid points [m] yi=zeros(1,length(xi)); vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); @@ -169,8 +145,7 @@ % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1) - disc_radius = 20*1000; % 20 km + disc_radius = 20*1000; % [m] [vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l); - % plot set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); figure1=figure('Position', [100, 100, 700, 400]); @@ -186,16 +161,13 @@ h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1); h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1); - % ISSM + % ISSM soln h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1); h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); - % first box ag1 = gca; leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)'); set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16); - % - set(gcf,'color','w'); - + set(gcf,'color','w'); %export_fig('Fig6.pdf'); -end % }}} +end Index: /issm/trunk-jpl/examples/SlrFarrell/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrFarrell/runme.m (revision 22817) +++ /issm/trunk-jpl/examples/SlrFarrell/runme.m (revision 22818) @@ -2,27 +2,24 @@ clear all; -steps=[1]; % [1:5]; +steps=[1]; % [1:5] -if any(steps==1) % Global mesh creation {{{ +if any(steps==1) 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] + resolution=150*1e3; % inital resolution [m] + 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. + md.mask=maskpsl(); + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes 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); @@ -31,5 +28,4 @@ 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), @@ -46,62 +42,51 @@ pos=find(distancemindistance_land); distance(pos2)=mindistance_land; - dist=min(maxdistance,distance); % max size 1000 km - %use distance to the coastline to refine mesh: + dist=min(maxdistance,distance); 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; + save ./Models/SlrFarrell_Mesh md; plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); - %export_fig('Fig1.pdf'); -end % }}} +end -if any(steps==2) % Define source {{{ +if any(steps==2) disp(' Step 2: Define source as in Farrell, 1972, Figure 1'); - md = loadmodel('./Models/SlrFarrell.Mesh'); + md = loadmodel('./Models/SlrFarrell_Mesh'); - % initial sea-level: 1 m RSL everywhere. - md.slr.sealevel=md.mask.ocean_levelset; + md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere md.slr.deltathickness=zeros(md.mesh.numberofelements,1); md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); - save ./Models/SlrFarrell.Loads md; + save ./Models/SlrFarrell_Loads md; - plotmodel (md,'data',md.slr.sealevel,'view',[90 90],... - 'title#all','Initial sea-level [m]'); - %export_fig('Fig2.pdf'); + plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]'); -end % }}} +end -if any(steps==3) % Parameterization {{{ +if any(steps==3) disp(' Step 3: Parameterization'); - md = loadmodel('./Models/SlrFarrell.Loads'); + md = loadmodel('./Models/SlrFarrell_Loads'); - % Love numbers and reference frame: CF or CM (choose one!) - nlove=10001; % up to 10,000 degree + nlove=10001; 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... + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1); 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); @@ -109,56 +94,48 @@ 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; + save ./Models/SlrFarrell_Parameterization md; -end % }}} +end -if any(steps==4) % Solve {{{ +if any(steps==4) disp(' Step 4: Solve Slr solver'); - md = loadmodel('./Models/SlrFarrell.Parameterization'); + md = loadmodel('./Models/SlrFarrell_Parameterization'); - % Cluster info md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); - % Choose different convergence threshold. [10% 1% 0.1%] to match Farrell 3 panels in Fig. 1 md.slr.reltol = 0.1/100; % per cent change in solution - % Solve md=solve(md,'Slr'); - save ./Models/SlrFarrell.Solution md; + save ./Models/SlrFarrell_Solution md; -end % }}} +end -if any(steps==5) % Plot solutions {{{ +if any(steps==5) disp(' Step 5: Plot solutions'); - md = loadmodel('./Models/SlrFarrell.Solution'); + md = loadmodel('./Models/SlrFarrell_Solution'); - % solutions. sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) - res = 1; % degree + res = 1; % [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 = 'natural'; % for smooth contour - F.ExtrapolationMethod = 'none'; + 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 + sol_grid(lat_grid>85 & sol_grid==0) =NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) @@ -170,5 +147,4 @@ geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]); plot(long,lat,'k'); hold off; - % define colormap, caxis, xlim etc {{{ c1=colorbar; colormap(flipud(haxby)); @@ -176,11 +152,9 @@ xlim([-170 170]); ylim([-85 85]); - % }}} grid on; title('Relative sea-level [% of GMSL]'); set(gcf,'color','w'); - %export_fig('Fig5.pdf'); -end % }}} +end Index: /issm/trunk-jpl/examples/SlrGRACE/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrGRACE/runme.m (revision 22817) +++ /issm/trunk-jpl/examples/SlrGRACE/runme.m (revision 22818) @@ -3,43 +3,24 @@ addpath('../Data','../Functions'); -steps=[0]; % [0:7]; - -if any(steps==0) % Download GRACE land_mass data {{{ - disp(' Step 0: Download GRACE land_mass data'); - - % Connect to ftp server and download - f = ftp('podaac-ftp.jpl.nasa.gov'); - %dir(f) - cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); - mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); - - !mv *.nc '../Data/' - - % display the content of the netcdf file. - %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') - -end % }}} - -if any(steps==1) % Global mesh creation {{{ +steps=[1]; % [1:7] + +if any(steps==1) 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: + resolution=150*1e3; % inital resolution [m] + 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] + 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; % refine mesh... {{{ - - %figure out mask: + md.mask=maskpsl(); + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] + + for i=1:numrefine, + 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); @@ -48,5 +29,4 @@ 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), @@ -63,75 +43,56 @@ pos=find(distancemindistance_land); distance(pos2)=mindistance_land; - dist=min(maxdistance,distance); % max size 1000 km - %use distance to the coastline to refine mesh: + dist=min(maxdistance,distance); md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); - end % }}} - - %figure out mask: + + end + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); - save ./Models/SlrGRACE.Mesh md; + save ./Models/SlrGRACE_Mesh md; plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); - %export_fig('Fig1.pdf'); - -end % }}} - -if any(steps==2) % Define loads {{{ + +end + +if any(steps==2) disp(' Step 2: Define loads in meters of ice height equivalent'); - md = loadmodel('./Models/SlrGRACE.Mesh'); - - % define time interval for analysis + md = loadmodel('./Models/SlrGRACE_Mesh'); + year_month = 2007+15/365; time_range = [year_month year_month]; - % map GRACE water load on to the mesh for the seleted month(s) water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); - md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent - - save ./Models/SlrGRACE.Loads md; - - plotmodel (md,'data',md.slr.deltathickness,... - 'view',[90 -90],'caxis',[-.1 .1],... - 'title','Ice height equivalent [m]'); - %export_fig('Fig2.pdf'); - -end % }}} - -if any(steps==3) % Parameterization {{{ + md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + + save ./Models/SlrGRACE_Loads md; + + plotmodel (md,'data',md.slr.deltathickness,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]'); + +end + +if any(steps==3) disp(' Step 3: Parameterization'); - md = loadmodel('./Models/SlrGRACE.Loads'); - - % Love numbers and reference frame: CF or CM (choose one!) - nlove=10001; % up to 10,000 degree + md = loadmodel('./Models/SlrGRACE_Loads'); + + nlove=10001; 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.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded + md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); pos=find(md.slr.deltathickness~=0); - md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; md.mask.land_levelset = 1-md.mask.ocean_levelset; - % initialize md.slr.sealevel=zeros(md.mesh.numberofvertices,1); - % steric rate md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); - % solution is scaled according to "exact" area of the oceans md.slr.ocean_area_scaling = 1; - % convergence criteria - %md.slr.reltol=NaN; - %md.slr.abstol=1e-4; - - %% 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); @@ -139,42 +100,37 @@ 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='SlrGRACE'; - %% IGNORE BUT DO NOT DELETE %% }}} - - save ./Models/SlrGRACE.Parameterization md; - -end % }}} - -if any(steps==4) % Solve {{{ + + save ./Models/SlrGRACE_Parameterization md; + +end + +if any(steps==4) disp(' Step 4: Solve Slr solver'); - md = loadmodel('./Models/SlrGRACE.Parameterization'); - - % What physics to capture? + md = loadmodel('./Models/SlrGRACE_Parameterization'); + md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1; - % Cluster info md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); - % Solve md=solve(md,'Slr'); - save ./Models/SlrGRACE.Solution md; - -end % }}} - -if any(steps==5) % Plot solutions {{{ + save ./Models/SlrGRACE_Solution md; + +end + +if any(steps==5) disp(' Step 5: Plot solutions'); - md = loadmodel('./Models/SlrGRACE.Solution'); - - % loads and solutions. - sol1 = md.slr.deltathickness*100; % WEH cm + md = loadmodel('./Models/SlrGRACE_Solution'); + + sol1 = md.slr.deltathickness*100; % [cm] sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] @@ -182,7 +138,6 @@ fig_name={'Fig_dH.pdf','Fig_slr.pdf'}; - res = 1.0; % degree - - % Make a grid of lats and lons, based on the min and max of the original vectors + res = 1.0; + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); sol_grid = zeros(size(lat_grid)); @@ -191,7 +146,5 @@ sol=eval(sprintf('sol%d',kk)); - % if data are on elements, map those on to the vertices {{{ if length(sol)==md.mesh.numberofelements - % map on to the vertices for jj=1:md.mesh.numberofelements ii=(jj-1)*3; @@ -203,15 +156,13 @@ end sol=temp'; - end % }}} - - % Make a interpolation object + end + 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 + sol_grid(lat_grid>85 & sol_grid==0) = NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) @@ -219,12 +170,10 @@ gcf; load coast; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; - % mask out land or oceans {{{ if (kk==1) geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); else geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); - end % }}} + end plot(long,lat,'k'); hold off; - % define colormap, caxis, xlim etc {{{ c1=colorbar; colormap('haxby'); @@ -232,30 +181,25 @@ xlim([-180 180]); ylim([-90 90]); - % }}} grid on; title(sol_name(kk)); set(gcf,'color','w'); - %export_fig(fig_name{kk}); end -end % }}} - -if any(steps==6) % Transient {{{ +end + +if any(steps==6) disp(' Step 6: Transient run'); - md = loadmodel('./Models/SlrGRACE.Parameterization'); - - % read GRACE data during the chosen time period << steps=2 >> + md = loadmodel('./Models/SlrGRACE_Parameterization'); + disp('Projecting loads onto the mesh...'); time_range = 2007 + [15 350]/365; water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); - % define ice load [ne,nt]=size(water_load); md.slr.deltathickness = zeros(ne+1,nt); - md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent + md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; md.slr.deltathickness(ne+1,:)=[1:nt]; % times - % define transient run md.transient.issmb=0; md.transient.ismasstransport=0; @@ -264,5 +208,4 @@ md.transient.isslr=1; - % time stepping, output frequency etc. md.timestepping.start_time=-1; md.timestepping.final_time=nt; @@ -271,26 +214,22 @@ md.settings.output_frequency=1; - % Cluster info md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); - % Solve md=solve(md,'tr'); - save ./Models/SlrGRACE.Transient md; - -end % }}} - -if any(steps==7) % Plot transient {{{ + save ./Models/SlrGRACE_Transient md; + +end + +if any(steps==7) disp(' Step 7: Plot transient'); - md = loadmodel('./Models/SlrGRACE.Transient'); - - time = md.slr.deltathickness(end,:); % year - - % loads and solutions. + md = loadmodel('./Models/SlrGRACE_Transient'); + + time = md.slr.deltathickness(end,:); + for tt=1:length(time) gmsl(tt) = md.results.TransientSolution(tt).SealevelEustatic*1000; % GMSL rate mm/yr sol1(:,tt) = md.slr.deltathickness(1:end-1,tt)*100; % ice equivalent height [cm/yr] - %% something weired happened here!!! first month solution is duplicated. Use offset of 1. Will fix it asap! sol2(:,tt) = md.results.TransientSolution(tt+1).Sealevel*1000; % mm/yr end @@ -298,7 +237,6 @@ movie_name = {'Movie_dH.avi','Movie_slr.avi'}; - res = 1.0; % degree - - % Make a grid of lats and lons, based on the min and max of the original vectors + res = 1.0; + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); sol_grid = zeros(size(lat_grid)); @@ -307,7 +245,5 @@ sol=eval(sprintf('sol%d',kk)); - % if data are on elements, map those on to the vertices {{{ if length(sol)==md.mesh.numberofelements - % map on to the vertices for jj=1:md.mesh.numberofelements ii=(jj-1)*3; @@ -319,5 +255,5 @@ end sol=temp2; - end % }}} + end vidObj = VideoWriter(movie_name{kk}); @@ -326,13 +262,11 @@ for jj=1:length(time) - % Make a interpolation object F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); 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 + sol_grid(lat_grid>85 & sol_grid==0) = NaN; set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) @@ -340,12 +274,10 @@ gcf; load coast; cla; pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; - % mask out land or oceans {{{ if (kk==1) geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); else geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); - end % }}} + end plot(long,lat,'k'); hold off; - % define colormap, caxis, xlim etc {{{ c1=colorbar; colormap('haxby'); @@ -353,13 +285,10 @@ xlim([-180 180]); ylim([-90 90]); - % }}} grid on; title(sol_name(kk)); set(gcf,'color','w'); writeVideo(vidObj,getframe(gcf)); - close % close current figure + close end - - % Close the file. close(vidObj); end @@ -371,7 +300,6 @@ ylabel('GMSL [mm/yr]'); set(gcf,'color','w'); - %export_fig('Fig7.pdf'); -end % }}} - +end +