Changeset 22818
- Timestamp:
- 05/31/18 09:58:48 (7 years ago)
- Location:
- issm/trunk-jpl/examples
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/examples/EsaGRACE/runme.m
r22815 r22818 3 3 addpath('../Data','../Functions'); 4 4 5 steps=[ 0]; % [0:5];5 steps=[1]; % [1:5] 6 6 7 if any(steps==0) % Download GRACE land_mass data {{{ 8 disp(' Step 0: Download GRACE land_mass data'); 9 10 % Connect to ftp server and download 11 f = ftp('podaac-ftp.jpl.nasa.gov'); 12 %dir(f) 13 cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); 14 mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); 15 16 !mv *.nc '../Data/' 17 18 % display the content of the netcdf file. 19 %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') 20 21 end % }}} 22 23 if any(steps==1) % Global mesh creation {{{ 7 if any(steps==1) 24 8 disp(' Step 1: Global mesh creation'); 25 9 26 resolution=300; % km. uniform meshing...27 radius = 6.371012*10^3; % mean radius of Earth, km10 resolution=300; % [km] 11 radius = 6.371012*10^3; % [km] 28 12 29 %mesh earth:30 13 md=model; 31 md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset14 md.mask=maskpsl(); 32 15 md.mesh=gmshplanet('radius',radius,'resolution',resolution); 33 16 34 %figure out mask:35 17 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 36 18 37 save ./Models/EsaGRACE .Mesh md;19 save ./Models/EsaGRACE_Mesh md; 38 20 39 21 plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); 40 %export_fig('Fig1.pdf');41 22 42 end % }}}23 end 43 24 44 if any(steps==2) % Define loads {{{25 if any(steps==2) 45 26 disp(' Step 2: Define loads in meters of ice height equivalent'); 46 md = loadmodel('./Models/EsaGRACE .Mesh');27 md = loadmodel('./Models/EsaGRACE_Mesh'); 47 28 48 % define time interval for analysis49 29 year_month = 2007+15/365; 50 30 time_range = [year_month year_month]; 51 31 52 % map GRACE water load on to the mesh for the seleted month(s)53 32 water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); 54 33 55 34 md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent 56 35 57 save ./Models/EsaGRACE .Loads md;36 save ./Models/EsaGRACE_Loads md; 58 37 59 38 plotmodel (md,'data',md.esa.deltathickness,... 60 39 'view',[90 -90],'caxis',[-.1 .1],... 61 40 'title','Ice height equivalent [m]'); 62 %export_fig('Fig2.pdf');63 41 64 end % }}}42 end 65 43 66 if any(steps==3) % Parameterization {{{44 if any(steps==3) 67 45 disp(' Step 3: Parameterization'); 68 md = loadmodel('./Models/EsaGRACE .Loads');46 md = loadmodel('./Models/EsaGRACE_Loads'); 69 47 70 % Love numbers and reference frame: CF or CM (choose one!) 71 nlove=10001; % up to 10,000 degree 48 nlove=10001; 72 49 md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[]; 73 50 md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[]; 74 51 75 % Mask: for computational efficiency only those elements that have loads are convolved! 76 md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded 52 md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); 77 53 md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 78 54 pos=find(md.esa.deltathickness~=0); 79 md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads55 md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; 80 56 md.mask.land_levelset = 1-md.mask.ocean_levelset; 81 57 82 %% IGNORE BUT DO NOT DELETE %% {{{83 % Geometry: Important only when you want to couple with Ice Flow Model84 58 di=md.materials.rho_ice/md.materials.rho_water; 85 59 md.geometry.thickness=ones(md.mesh.numberofvertices,1); … … 87 61 md.geometry.base=md.geometry.surface-md.geometry.thickness; 88 62 md.geometry.bed=md.geometry.base; 89 % Materials:63 90 64 md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); 91 65 md.materials.rheology_B=paterson(md.initialization.temperature); 92 66 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); 93 % Miscellaneous:67 94 68 md.miscellaneous.name='EsaGRACE'; 95 %% IGNORE BUT DO NOT DELETE %% }}}96 69 97 save ./Models/EsaGRACE .Parameterization md;70 save ./Models/EsaGRACE_Parameterization md; 98 71 99 end % }}}72 end 100 73 101 if any(steps==4) % Solve {{{74 if any(steps==4) 102 75 disp(' Step 4: Solve Esa solver'); 103 md = loadmodel('./Models/EsaGRACE .Parameterization');76 md = loadmodel('./Models/EsaGRACE_Parameterization'); 104 77 105 % Request outputs106 78 md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'}; 107 79 108 % Cluster info109 80 md.cluster=generic('name',oshostname(),'np',3); 110 81 md.verbose=verbose('111111111'); 111 82 112 % Solve113 83 md=solve(md,'Esa'); 114 84 115 save ./Models/EsaGRACE .Solution md;85 save ./Models/EsaGRACE_Solution md; 116 86 117 end % }}}87 end 118 88 119 if any(steps==5) % Plot solutions {{{89 if any(steps==5) 120 90 disp(' Step 5: Plot solutions'); 121 md = loadmodel('./Models/EsaGRACE .Solution');91 md = loadmodel('./Models/EsaGRACE_Solution'); 122 92 123 % loads and solutions. 124 sol1 = md.esa.deltathickness*100; % WEH cm 125 sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] 126 sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] 127 sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] 93 sol1 = md.esa.deltathickness*100; % [cm] 94 sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] 95 sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] 96 sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] 128 97 129 98 sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',... … … 131 100 fig_name={'Fig_dH.pdf','Fig_vert.pdf','Fig_horzNS.pdf','Fig_horzEW.pdf'}; 132 101 133 res = 1.0; % degree102 res = 1.0; % [degree] 134 103 135 % Make a grid of lats and lons, based on the min and max of the original vectors136 104 [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); 137 105 sol_grid = zeros(size(lat_grid)); … … 140 108 sol=eval(sprintf('sol%d',kk)); 141 109 142 % if data are on elements, map those on to the vertices {{{143 110 if length(sol)==md.mesh.numberofelements 144 % map on to the vertices145 111 for jj=1:md.mesh.numberofelements 146 112 ii=(jj-1)*3; … … 152 118 end 153 119 sol=temp'; 154 end % }}}120 end 155 121 156 % Make a interpolation object157 122 F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 158 123 F.Method = 'linear'; 159 124 F.ExtrapolationMethod = 'linear'; 160 125 161 % Do the interpolation to get gridded solutions...162 126 sol_grid = F(lat_grid, lon_grid); 163 127 sol_grid(isnan(sol_grid))=0; 164 sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan128 sol_grid(lat_grid>85 & sol_grid==0)=NaN; 165 129 166 130 set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) … … 168 132 gcf; load coast; cla; 169 133 pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; 170 % mask out oceans {{{171 134 if (kk==1) 172 135 geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); 173 end % }}}136 end 174 137 plot(long,lat,'k'); hold off; 175 % colormap, xlim etc {{{176 138 c1=colorbar; 177 139 colormap('haxby'); … … 179 141 xlim([-180 180]); 180 142 ylim([-90 90]); 181 % }}}182 143 grid on; 183 144 title(sol_name(kk)); 184 145 set(gcf,'color','w'); 185 186 146 %export_fig(fig_name{kk}); 187 147 end 188 148 189 end % }}}149 end 190 150 -
issm/trunk-jpl/examples/EsaWahr/runme.m
r22815 r22818 3 3 addpath('../Functions'); 4 4 5 steps=[0]; % [0:6] ;5 steps=[0]; % [0:6] 6 6 7 if any(steps==0) % Simple mesh creation {{{7 if any(steps==0) 8 8 disp(' Step 0: Mesh creation'); 9 9 10 % Generate initial uniform mesh 11 md=roundmesh(model,100000,10000); % Domain radius and element size (meters) 10 md=roundmesh(model,100000,10000); % Domain radius and element size [m] 12 11 13 save ./Models/EsaWahr .Mesh md;12 save ./Models/EsaWahr_Mesh md; 14 13 15 14 plotmodel(md,'data','mesh'); 16 %export_fig('Fig0.pdf');17 15 18 end % }}}16 end 19 17 20 if any(steps==1) % {{{ Anisotropic mesh creation18 if any(steps==1) 21 19 disp(' Step 1: Anisotropic mesh creation'); 22 20 23 % Generate initial uniform mesh 24 md=roundmesh(model,100000,1000); % Domain radius and element size (meters) 25 %plotmodel(md,'data','mesh'); 21 md=roundmesh(model,100000,1000); 26 22 27 % Define a fake field to have anisotropic meshing 28 disc_radius=20*1000; % km => m 29 rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); % radial distance [m] 23 disc_radius=20*1000; 24 rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); 30 25 field = abs(rad_dist-disc_radius); 31 %plotmodel(md,'data',field);32 26 33 md = bamg(md,'field',field,'err',50,'hmax',10000); % error in field in m27 md = bamg(md,'field',field,'err',50,'hmax',10000); 34 28 35 save ./Models/EsaWahr .Mesh md;29 save ./Models/EsaWahr_Mesh md; 36 30 37 31 plotmodel (md,'data','mesh'); 38 %export_fig('Fig1.pdf');39 32 40 end % }}}33 end 41 34 42 if any(steps==2) % Define loads {{{35 if any(steps==2) 43 36 disp(' Step 2: Define loads'); 44 md = loadmodel('./Models/EsaWahr .Mesh');37 md = loadmodel('./Models/EsaWahr_Mesh'); 45 38 46 % primer47 39 rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; 48 40 49 % elemental centroids50 41 index=md.mesh.elements; 51 42 x_cent=mean(md.mesh.x(index),2); 52 43 y_cent=mean(md.mesh.y(index),2); 53 44 54 % Define a 20-km radius disc and unload it by 1 meter water height equivalent uniformly55 45 md.esa.deltathickness = zeros(md.mesh.numberofelements,1); 56 disc_radius=20; % km57 rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000; % radial distance in km58 md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i; % 1 m water withdrawl46 disc_radius=20; % [km] 47 rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000; 48 md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i; 59 49 60 save ./Models/EsaWahr .Loads md;50 save ./Models/EsaWahr_Loads md; 61 51 62 52 plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]'); 63 %export_fig('Fig2.pdf');64 53 65 end % }}}54 end 66 55 67 if any(steps==3) % Parameterization {{{56 if any(steps==3) 68 57 disp(' Step 3: Parameterization'); 69 md = loadmodel('./Models/EsaWahr .Loads');58 md = loadmodel('./Models/EsaWahr_Loads'); 70 59 71 % Love numbers and reference frame: CF or CM (choose one!) 72 nlove=10001; % up to 10,000 degree 60 nlove=10001; 73 61 md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[]; 74 62 md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[]; 75 63 76 % Mask: for computational efficiency only those elements that have loads are convolved! 77 md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); % -1 = ice loads 78 md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded 64 md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); 65 md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); 79 66 80 %% IGNORE BUT DO NOT DELETE %% {{{81 % Geometry: Important only when you want to couple with Ice Flow Model82 67 di=md.materials.rho_ice/md.materials.rho_water; 83 68 md.geometry.thickness=ones(md.mesh.numberofvertices,1); … … 85 70 md.geometry.base=md.geometry.surface-md.geometry.thickness; 86 71 md.geometry.bed=md.geometry.base; 87 % Materials:72 88 73 md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); 89 74 md.materials.rheology_B=paterson(md.initialization.temperature); 90 75 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); 91 % Miscellaneous:76 92 77 md.miscellaneous.name='EsaWahr'; 93 %% IGNORE BUT DO NOT DELETE %% }}}94 78 95 save ./Models/EsaWahr .Parameterization md;79 save ./Models/EsaWahr_Parameterization md; 96 80 97 end % }}}81 end 98 82 99 if any(steps==4) % Solve {{{83 if any(steps==4) 100 84 disp(' Step 4: Solve Esa solver'); 101 md = loadmodel('./Models/EsaWahr .Parameterization');85 md = loadmodel('./Models/EsaWahr_Parameterization'); 102 86 103 % Request outputs104 87 md.esa.requested_outputs = {'EsaUmotion','EsaXmotion','EsaYmotion'}; 105 88 106 % Cluster info107 89 md.cluster=generic('name',oshostname(),'np',3); 108 90 md.verbose=verbose('111111111'); 109 91 110 % Solve111 92 md=solve(md,'Esa'); 112 93 113 save ./Models/EsaWahr .Solution md;94 save ./Models/EsaWahr_Solution md; 114 95 115 end % }}}96 end 116 97 117 if any(steps==5) % Plot solutions {{{98 if any(steps==5) 118 99 disp(' Step 5: Plot solutions'); 119 md = loadmodel('./Models/EsaWahr .Solution');100 md = loadmodel('./Models/EsaWahr_Solution'); 120 101 121 % m => mm 122 vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 123 horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] 124 horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] 125 horz = sqrt(horz_n.^2+horz_e.^2); 102 vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 103 horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] 104 horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] 105 horz = sqrt(horz_n.^2+horz_e.^2); % [mm] 126 106 127 % plot128 107 set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); 129 108 figure('Position', [100, 100, 800, 600]); … … 147 126 'axispos',[0.505 0.02 0.47 0.47],... 148 127 'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); 149 150 128 %export_fig('Fig5.pdf'); 151 129 152 end % }}}130 end 153 131 154 if any(steps==6) % Compare results against semi-analytic solutions {{{132 if any(steps==6) 155 133 disp(' Step 6: Compare results against Wahr semi-analytic solutions'); 156 md = loadmodel('./Models/EsaWahr .Solution');134 md = loadmodel('./Models/EsaWahr_Solution'); 157 135 158 % numerical solutions 159 % m => mm 160 vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 161 horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] 162 horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] 163 horz = sqrt(horz_n.^2+horz_e.^2); 164 % grid data from the center 165 xi=[0:500:100000]; % 500 m resolution 136 vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 137 horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] 138 horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] 139 horz = sqrt(horz_n.^2+horz_e.^2); % [mm] 140 141 xi=[0:500:100000]; % grid points [m] 166 142 yi=zeros(1,length(xi)); 167 143 vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); … … 169 145 170 146 % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1) 171 disc_radius = 20*1000; % 20 km147 disc_radius = 20*1000; % [m] 172 148 [vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l); 173 149 174 % plot175 150 set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); 176 151 figure1=figure('Position', [100, 100, 700, 400]); … … 186 161 h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1); 187 162 h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1); 188 % ISSM 163 % ISSM soln 189 164 h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1); 190 165 h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); 191 % first box192 166 ag1 = gca; 193 167 leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)'); 194 168 set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16); 195 % 196 set(gcf,'color','w'); 197 169 set(gcf,'color','w'); 198 170 %export_fig('Fig6.pdf'); 199 171 200 end % }}}172 end 201 173 -
issm/trunk-jpl/examples/SlrFarrell/runme.m
r22815 r22818 2 2 clear all; 3 3 4 steps=[1]; % [1:5] ;4 steps=[1]; % [1:5] 5 5 6 if any(steps==1) % Global mesh creation {{{6 if any(steps==1) 7 7 disp(' Step 1: Global mesh creation'); 8 8 9 9 numrefine=1; 10 resolution=150*1e3; % inital resolution [m]. It determines, e.g., whether we capture small islands.11 radius = 6.371012*10^6; % mean radius of Earth, m12 mindistance_coast=150*1e3; 13 mindistance_land=300*1e3; 14 maxdistance=600*1e3; 10 resolution=150*1e3; % inital resolution [m] 11 radius = 6.371012*10^6; % mean radius of Earth, m 12 mindistance_coast=150*1e3; % coastal resolution [m] 13 mindistance_land=300*1e3; % resolution on the continents [m] 14 maxdistance=600*1e3; % max element size (on mid-oceans) [m] 15 15 16 %mesh earth:17 16 md=model; 18 md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset19 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km.17 md.mask=maskpsl(); 18 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] 20 19 21 20 for i=1:numrefine, 22 21 23 %figure out mask:24 22 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 25 23 26 %figure out distance to the coastline, in lat,long (not x,y,z):27 24 distance=zeros(md.mesh.numberofvertices,1); 28 25 … … 31 28 32 29 for j=1:md.mesh.numberofvertices 33 %figure out nearest coastline (using the great circle distance)34 30 phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; 35 31 if md.mask.ocean_levelset(j), … … 46 42 pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast; 47 43 48 % refine on the continents49 44 pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); 50 45 distance(pos2)=mindistance_land; 51 46 52 dist=min(maxdistance,distance); % max size 1000 km 53 %use distance to the coastline to refine mesh: 47 dist=min(maxdistance,distance); 54 48 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); 49 55 50 end 56 51 57 %figure out mask:58 52 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 59 53 60 save ./Models/SlrFarrell .Mesh md;54 save ./Models/SlrFarrell_Mesh md; 61 55 62 56 plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); 63 %export_fig('Fig1.pdf');64 57 65 end % }}}58 end 66 59 67 if any(steps==2) % Define source {{{60 if any(steps==2) 68 61 disp(' Step 2: Define source as in Farrell, 1972, Figure 1'); 69 md = loadmodel('./Models/SlrFarrell .Mesh');62 md = loadmodel('./Models/SlrFarrell_Mesh'); 70 63 71 % initial sea-level: 1 m RSL everywhere. 72 md.slr.sealevel=md.mask.ocean_levelset; 64 md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere 73 65 74 66 md.slr.deltathickness=zeros(md.mesh.numberofelements,1); 75 67 md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); 76 68 77 save ./Models/SlrFarrell .Loads md;69 save ./Models/SlrFarrell_Loads md; 78 70 79 plotmodel (md,'data',md.slr.sealevel,'view',[90 90],... 80 'title#all','Initial sea-level [m]'); 81 %export_fig('Fig2.pdf'); 71 plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]'); 82 72 83 end % }}}73 end 84 74 85 if any(steps==3) % Parameterization {{{75 if any(steps==3) 86 76 disp(' Step 3: Parameterization'); 87 md = loadmodel('./Models/SlrFarrell .Loads');77 md = loadmodel('./Models/SlrFarrell_Loads'); 88 78 89 % Love numbers and reference frame: CF or CM (choose one!) 90 nlove=10001; % up to 10,000 degree 79 nlove=10001; 91 80 md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[]; 92 81 md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[]; 93 82 md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[]; 94 83 95 % Mask: for computational efficiency only those elements that have loads are convolved!96 84 md.mask.land_levelset = 1-md.mask.ocean_levelset; 97 % fake ice load in one element! 98 md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); % no ice 99 md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1); % floated... 85 md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 86 md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1); 100 87 pos=find(md.mesh.lat <-80); 101 88 md.mask.ice_levelset(pos(1))=-1; % ice yes! 102 89 md.mask.groundedice_levelset(pos(1))=1; % ice grounded! 103 90 104 %% IGNORE BUT DO NOT DELETE %% {{{105 % Geometry: Important only when you want to couple with Ice Flow Model106 91 di=md.materials.rho_ice/md.materials.rho_water; 107 92 md.geometry.thickness=ones(md.mesh.numberofvertices,1); … … 109 94 md.geometry.base=md.geometry.surface-md.geometry.thickness; 110 95 md.geometry.bed=md.geometry.base; 111 % Materials:96 112 97 md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); 113 98 md.materials.rheology_B=paterson(md.initialization.temperature); 114 99 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); 115 % Miscellaneous:100 116 101 md.miscellaneous.name='SlrFarrell'; 117 %% IGNORE BUT DO NOT DELETE %% }}}118 102 119 save ./Models/SlrFarrell .Parameterization md;103 save ./Models/SlrFarrell_Parameterization md; 120 104 121 end % }}}105 end 122 106 123 if any(steps==4) % Solve {{{107 if any(steps==4) 124 108 disp(' Step 4: Solve Slr solver'); 125 md = loadmodel('./Models/SlrFarrell .Parameterization');109 md = loadmodel('./Models/SlrFarrell_Parameterization'); 126 110 127 % Cluster info128 111 md.cluster=generic('name',oshostname(),'np',3); 129 112 md.verbose=verbose('111111111'); 130 113 131 % Choose different convergence threshold. [10% 1% 0.1%] to match Farrell 3 panels in Fig. 1132 114 md.slr.reltol = 0.1/100; % per cent change in solution 133 115 134 % Solve135 116 md=solve(md,'Slr'); 136 117 137 save ./Models/SlrFarrell .Solution md;118 save ./Models/SlrFarrell_Solution md; 138 119 139 end % }}}120 end 140 121 141 if any(steps==5) % Plot solutions {{{122 if any(steps==5) 142 123 disp(' Step 5: Plot solutions'); 143 md = loadmodel('./Models/SlrFarrell .Solution');124 md = loadmodel('./Models/SlrFarrell_Solution'); 144 125 145 % solutions.146 126 sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) 147 127 148 res = 1; % degree128 res = 1; % [degree] 149 129 150 % Make a grid of lats and lons, based on the min and max of the original vectors151 130 [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); 152 131 sol_grid = zeros(size(lat_grid)); 153 132 154 % Make a interpolation object155 133 F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 156 F.Method = ' natural'; % for smooth contour157 F.ExtrapolationMethod = ' none';134 F.Method = 'linear'; 135 F.ExtrapolationMethod = 'linear'; 158 136 159 % Do the interpolation to get gridded solutions...160 137 sol_grid = F(lat_grid, lon_grid); 161 138 sol_grid(isnan(sol_grid))=0; 162 sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan139 sol_grid(lat_grid>85 & sol_grid==0) =NaN; 163 140 164 141 set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) … … 170 147 geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]); 171 148 plot(long,lat,'k'); hold off; 172 % define colormap, caxis, xlim etc {{{173 149 c1=colorbar; 174 150 colormap(flipud(haxby)); … … 176 152 xlim([-170 170]); 177 153 ylim([-85 85]); 178 % }}}179 154 grid on; 180 155 title('Relative sea-level [% of GMSL]'); 181 156 set(gcf,'color','w'); 182 183 157 %export_fig('Fig5.pdf'); 184 158 185 end % }}}159 end 186 160 -
issm/trunk-jpl/examples/SlrGRACE/runme.m
r22812 r22818 3 3 addpath('../Data','../Functions'); 4 4 5 steps=[0]; % [0:7]; 6 7 if any(steps==0) % Download GRACE land_mass data {{{ 8 disp(' Step 0: Download GRACE land_mass data'); 9 10 % Connect to ftp server and download 11 f = ftp('podaac-ftp.jpl.nasa.gov'); 12 %dir(f) 13 cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); 14 mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); 15 16 !mv *.nc '../Data/' 17 18 % display the content of the netcdf file. 19 %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') 20 21 end % }}} 22 23 if any(steps==1) % Global mesh creation {{{ 5 steps=[1]; % [1:7] 6 7 if any(steps==1) 24 8 disp(' Step 1: Global mesh creation'); 25 9 26 10 numrefine=1; 27 resolution=150*1e3; % inital resolution [m]. It determines, e.g., whether we capture small islands. 28 radius = 6.371012*10^6; % mean radius of Earth, m 29 mindistance_coast=150*1e3; % coastal resolution [m] 30 mindistance_land=300*1e3; % resolution on the continents [m] 31 maxdistance=600*1e3; % max element size (on mid-oceans) [m] 32 33 %mesh earth: 11 resolution=150*1e3; % inital resolution [m] 12 radius = 6.371012*10^6; % mean radius of Earth, m 13 mindistance_coast=150*1e3; % coastal resolution [m] 14 mindistance_land=300*1e3; % resolution on the continents [m] 15 maxdistance=600*1e3; % max element size (on mid-oceans) [m] 16 34 17 md=model; 35 md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset 36 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. 37 38 for i=1:numrefine; % refine mesh... {{{ 39 40 %figure out mask: 18 md.mask=maskpsl(); 19 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] 20 21 for i=1:numrefine, 22 41 23 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 42 24 43 %figure out distance to the coastline, in lat,long (not x,y,z):44 25 distance=zeros(md.mesh.numberofvertices,1); 45 26 … … 48 29 49 30 for j=1:md.mesh.numberofvertices 50 %figure out nearest coastline (using the great circle distance)51 31 phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; 52 32 if md.mask.ocean_levelset(j), … … 63 43 pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast; 64 44 65 % refine on the continents66 45 pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); 67 46 distance(pos2)=mindistance_land; 68 47 69 dist=min(maxdistance,distance); % max size 1000 km 70 %use distance to the coastline to refine mesh: 48 dist=min(maxdistance,distance); 71 49 md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); 72 end % }}} 73 74 %figure out mask: 50 51 end 52 75 53 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 76 54 77 save ./Models/SlrGRACE .Mesh md;55 save ./Models/SlrGRACE_Mesh md; 78 56 79 57 plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); 80 %export_fig('Fig1.pdf'); 81 82 end % }}} 83 84 if any(steps==2) % Define loads {{{ 58 59 end 60 61 if any(steps==2) 85 62 disp(' Step 2: Define loads in meters of ice height equivalent'); 86 md = loadmodel('./Models/SlrGRACE.Mesh'); 87 88 % define time interval for analysis 63 md = loadmodel('./Models/SlrGRACE_Mesh'); 64 89 65 year_month = 2007+15/365; 90 66 time_range = [year_month year_month]; 91 67 92 % map GRACE water load on to the mesh for the seleted month(s)93 68 water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); 94 69 95 md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent 96 97 save ./Models/SlrGRACE.Loads md; 98 99 plotmodel (md,'data',md.slr.deltathickness,... 100 'view',[90 -90],'caxis',[-.1 .1],... 101 'title','Ice height equivalent [m]'); 102 %export_fig('Fig2.pdf'); 103 104 end % }}} 105 106 if any(steps==3) % Parameterization {{{ 70 md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; 71 72 save ./Models/SlrGRACE_Loads md; 73 74 plotmodel (md,'data',md.slr.deltathickness,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]'); 75 76 end 77 78 if any(steps==3) 107 79 disp(' Step 3: Parameterization'); 108 md = loadmodel('./Models/SlrGRACE.Loads'); 109 110 % Love numbers and reference frame: CF or CM (choose one!) 111 nlove=10001; % up to 10,000 degree 80 md = loadmodel('./Models/SlrGRACE_Loads'); 81 82 nlove=10001; 112 83 md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[]; 113 84 md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[]; 114 85 md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[]; 115 86 116 % Mask: for computational efficiency only those elements that have loads are convolved! 117 md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded 87 md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); 118 88 md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 119 89 pos=find(md.slr.deltathickness~=0); 120 md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads90 md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; 121 91 md.mask.land_levelset = 1-md.mask.ocean_levelset; 122 92 123 % initialize124 93 md.slr.sealevel=zeros(md.mesh.numberofvertices,1); 125 % steric rate126 94 md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); 127 % solution is scaled according to "exact" area of the oceans128 95 md.slr.ocean_area_scaling = 1; 129 96 130 % convergence criteria131 %md.slr.reltol=NaN;132 %md.slr.abstol=1e-4;133 134 %% IGNORE BUT DO NOT DELETE %% {{{135 % Geometry: Important only when you want to couple with Ice Flow Model136 97 di=md.materials.rho_ice/md.materials.rho_water; 137 98 md.geometry.thickness=ones(md.mesh.numberofvertices,1); … … 139 100 md.geometry.base=md.geometry.surface-md.geometry.thickness; 140 101 md.geometry.bed=md.geometry.base; 141 % Materials:102 142 103 md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); 143 104 md.materials.rheology_B=paterson(md.initialization.temperature); 144 105 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); 145 % Miscellaneous:106 146 107 md.miscellaneous.name='SlrGRACE'; 147 %% IGNORE BUT DO NOT DELETE %% }}} 148 149 save ./Models/SlrGRACE.Parameterization md; 150 151 end % }}} 152 153 if any(steps==4) % Solve {{{ 108 109 save ./Models/SlrGRACE_Parameterization md; 110 111 end 112 113 if any(steps==4) 154 114 disp(' Step 4: Solve Slr solver'); 155 md = loadmodel('./Models/SlrGRACE.Parameterization'); 156 157 % What physics to capture? 115 md = loadmodel('./Models/SlrGRACE_Parameterization'); 116 158 117 md.slr.rigid=1; 159 118 md.slr.elastic=1; 160 119 md.slr.rotation=1; 161 120 162 % Cluster info163 121 md.cluster=generic('name',oshostname(),'np',3); 164 122 md.verbose=verbose('111111111'); 165 123 166 % Solve167 124 md=solve(md,'Slr'); 168 125 169 save ./Models/SlrGRACE .Solution md;170 171 end % }}}172 173 if any(steps==5) % Plot solutions {{{126 save ./Models/SlrGRACE_Solution md; 127 128 end 129 130 if any(steps==5) 174 131 disp(' Step 5: Plot solutions'); 175 md = loadmodel('./Models/SlrGRACE.Solution'); 176 177 % loads and solutions. 178 sol1 = md.slr.deltathickness*100; % WEH cm 132 md = loadmodel('./Models/SlrGRACE_Solution'); 133 134 sol1 = md.slr.deltathickness*100; % [cm] 179 135 sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] 180 136 … … 182 138 fig_name={'Fig_dH.pdf','Fig_slr.pdf'}; 183 139 184 res = 1.0; % degree 185 186 % Make a grid of lats and lons, based on the min and max of the original vectors 140 res = 1.0; 141 187 142 [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); 188 143 sol_grid = zeros(size(lat_grid)); … … 191 146 sol=eval(sprintf('sol%d',kk)); 192 147 193 % if data are on elements, map those on to the vertices {{{194 148 if length(sol)==md.mesh.numberofelements 195 % map on to the vertices196 149 for jj=1:md.mesh.numberofelements 197 150 ii=(jj-1)*3; … … 203 156 end 204 157 sol=temp'; 205 end % }}} 206 207 % Make a interpolation object 158 end 159 208 160 F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 209 161 F.Method = 'linear'; 210 162 F.ExtrapolationMethod = 'linear'; 211 163 212 % Do the interpolation to get gridded solutions...213 164 sol_grid = F(lat_grid, lon_grid); 214 165 sol_grid(isnan(sol_grid))=0; 215 sol_grid(lat_grid>85 & sol_grid==0) = NaN; % set polar unphysical 0s to Nan166 sol_grid(lat_grid>85 & sol_grid==0) = NaN; 216 167 217 168 set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) … … 219 170 gcf; load coast; cla; 220 171 pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; 221 % mask out land or oceans {{{222 172 if (kk==1) 223 173 geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); 224 174 else 225 175 geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); 226 end % }}}176 end 227 177 plot(long,lat,'k'); hold off; 228 % define colormap, caxis, xlim etc {{{229 178 c1=colorbar; 230 179 colormap('haxby'); … … 232 181 xlim([-180 180]); 233 182 ylim([-90 90]); 234 % }}}235 183 grid on; 236 184 title(sol_name(kk)); 237 185 set(gcf,'color','w'); 238 239 186 %export_fig(fig_name{kk}); 240 187 end 241 188 242 end % }}}243 244 if any(steps==6) % Transient {{{189 end 190 191 if any(steps==6) 245 192 disp(' Step 6: Transient run'); 246 md = loadmodel('./Models/SlrGRACE.Parameterization'); 247 248 % read GRACE data during the chosen time period << steps=2 >> 193 md = loadmodel('./Models/SlrGRACE_Parameterization'); 194 249 195 disp('Projecting loads onto the mesh...'); 250 196 time_range = 2007 + [15 350]/365; 251 197 water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); 252 198 253 % define ice load254 199 [ne,nt]=size(water_load); 255 200 md.slr.deltathickness = zeros(ne+1,nt); 256 md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent201 md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; 257 202 md.slr.deltathickness(ne+1,:)=[1:nt]; % times 258 203 259 % define transient run260 204 md.transient.issmb=0; 261 205 md.transient.ismasstransport=0; … … 264 208 md.transient.isslr=1; 265 209 266 % time stepping, output frequency etc.267 210 md.timestepping.start_time=-1; 268 211 md.timestepping.final_time=nt; … … 271 214 md.settings.output_frequency=1; 272 215 273 % Cluster info274 216 md.cluster=generic('name',oshostname(),'np',3); 275 217 md.verbose=verbose('111111111'); 276 218 277 % Solve278 219 md=solve(md,'tr'); 279 220 280 save ./Models/SlrGRACE .Transient md;281 282 end % }}}283 284 if any(steps==7) % Plot transient {{{221 save ./Models/SlrGRACE_Transient md; 222 223 end 224 225 if any(steps==7) 285 226 disp(' Step 7: Plot transient'); 286 md = loadmodel('./Models/SlrGRACE.Transient'); 287 288 time = md.slr.deltathickness(end,:); % year 289 290 % loads and solutions. 227 md = loadmodel('./Models/SlrGRACE_Transient'); 228 229 time = md.slr.deltathickness(end,:); 230 291 231 for tt=1:length(time) 292 232 gmsl(tt) = md.results.TransientSolution(tt).SealevelEustatic*1000; % GMSL rate mm/yr 293 233 sol1(:,tt) = md.slr.deltathickness(1:end-1,tt)*100; % ice equivalent height [cm/yr] 294 %% something weired happened here!!! first month solution is duplicated. Use offset of 1. Will fix it asap!295 234 sol2(:,tt) = md.results.TransientSolution(tt+1).Sealevel*1000; % mm/yr 296 235 end … … 298 237 movie_name = {'Movie_dH.avi','Movie_slr.avi'}; 299 238 300 res = 1.0; % degree 301 302 % Make a grid of lats and lons, based on the min and max of the original vectors 239 res = 1.0; 240 303 241 [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); 304 242 sol_grid = zeros(size(lat_grid)); … … 307 245 sol=eval(sprintf('sol%d',kk)); 308 246 309 % if data are on elements, map those on to the vertices {{{310 247 if length(sol)==md.mesh.numberofelements 311 % map on to the vertices312 248 for jj=1:md.mesh.numberofelements 313 249 ii=(jj-1)*3; … … 319 255 end 320 256 sol=temp2; 321 end % }}}257 end 322 258 323 259 vidObj = VideoWriter(movie_name{kk}); … … 326 262 327 263 for jj=1:length(time) 328 % Make a interpolation object329 264 F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); 330 265 F.Method = 'linear'; 331 266 F.ExtrapolationMethod = 'linear'; 332 267 333 % Do the interpolation to get gridded solutions...334 268 sol_grid = F(lat_grid, lon_grid); 335 269 sol_grid(isnan(sol_grid))=0; 336 sol_grid(lat_grid>85 & sol_grid==0) = NaN; % set polar unphysical 0s to Nan270 sol_grid(lat_grid>85 & sol_grid==0) = NaN; 337 271 338 272 set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) … … 340 274 gcf; load coast; cla; 341 275 pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; 342 % mask out land or oceans {{{343 276 if (kk==1) 344 277 geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); 345 278 else 346 279 geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); 347 end % }}}280 end 348 281 plot(long,lat,'k'); hold off; 349 % define colormap, caxis, xlim etc {{{350 282 c1=colorbar; 351 283 colormap('haxby'); … … 353 285 xlim([-180 180]); 354 286 ylim([-90 90]); 355 % }}}356 287 grid on; 357 288 title(sol_name(kk)); 358 289 set(gcf,'color','w'); 359 290 writeVideo(vidObj,getframe(gcf)); 360 close % close current figure291 close 361 292 end 362 363 % Close the file.364 293 close(vidObj); 365 294 end … … 371 300 ylabel('GMSL [mm/yr]'); 372 301 set(gcf,'color','w'); 373 374 302 %export_fig('Fig7.pdf'); 375 303 376 end % }}}377 304 end 305
Note:
See TracChangeset
for help on using the changeset viewer.