Changeset 22813

Timestamp:

05/29/18 23:51:25 (7 years ago)

Author:

adhikari

Message:

CHG: Farrell solutions recovered

File:

• issm/trunk-jpl/examples/SlrFarrell/runme.m (modified) (6 diffs)

issm/trunk-jpl/examples/SlrFarrell/runme.m

@@ -1 +1 @@
 
 clear all;
-steps=[2]; % 
+steps=[5]; % 
 
 if any(steps==1) % Global mesh creation {{{ 
@@ -70 +70 @@
         % 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; 
@@ -85 +88 @@
         % 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(nlov+1:end)=[];
-        md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[];
-        md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[];
+        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.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.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 %% {{{
@@ -119 +124 @@
         md = loadmodel('./Models/SlrFarrell.Parameterization');
 
-        % Request outputs 
-        md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'};
-        
         % Cluster info 
         md.cluster=generic('name',oshostname(),'np',3); 
@@ -137 +139 @@
         md = loadmodel('./Models/SlrFarrell.Solution');
 
-        % loads and solutions. 
-        sol1 = md.slr.deltathickness*100; % WEH cm 
-        sol2 = md.results.SlrSolution.SlrUmotion*1000; % [mm] 
-        sol3 = md.results.SlrSolution.SlrNmotion*1000; % [mm] 
-        sol4 = md.results.SlrSolution.SlrEmotion*1000; % [mm] 
-        sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',...
-                'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; 
-
+        % solutions. 
+        sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m)  
+        
         res = 1.0; % degree 
 
@@ -151 +148 @@
         sol_grid = zeros(size(lat_grid)); 
 
-        for kk=1:4 
-                sol=eval(sprintf('sol%d',kk));
+        % Make a interpolation object
+        F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
+        F.Method = 'linear';
+        F.ExtrapolationMethod = 'linear'; 
         
-                % 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;
-                                pp(ii+1:ii+3)=md.mesh.elements(jj,:);
-                        end
-                        for jj=1:md.mesh.numberofvertices
-                                pos=ceil(find(pp==jj)/3); 
-                                temp(jj)=mean(sol(pos)); 
-                        end
-                        sol=temp'; 
-                end % }}}
+        % 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 
 
-                % Make a interpolation object
-                F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
-                F.Method = 'linear';
-                F.ExtrapolationMethod = 'linear'; 
+        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');
 
-                % 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; 
-                plot(long,lat,'k'); hold off; 
-                c1=colorbar;
-                colormap(jet); 
-                xlim([-180 180]); 
-                ylim([-90 90]); 
-                grid on; 
-                title(sol_name(kk)); 
-                set(gcf,'color','w');
-
-                %export_fig('Fig5.pdf'); 
-        end
+        %export_fig('Fig5.pdf'); 
 
 end % }}} 
 
-