Changeset 22803

Timestamp:

05/24/18 10:28:17 (7 years ago)

Author:

adhikari

Message:

CHG: added plotting scripts

File:

• issm/trunk-jpl/examples/EsaGRACE/runme.m (modified) (4 diffs)

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

@@ -1 +1 @@
 
 clear all;
-steps=[5]; % [1:6]; 
+steps=[0]; % [1:5]; 
 
 if any(steps==0) % Download GRACE land_mass data {{{
@@ -42 +42 @@
         md = loadmodel('./Models/EsaGRACE.Mesh');
 
+        % define time interval for analysis 
         year_month = 2007+15/365;
-        time_min=year_month; 
-        time_max=year_month; 
-
+        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_min,time_max); 
+        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 
@@ -117 +117 @@
         md = loadmodel('./Models/EsaGRACE.Solution');
 
-        % solutions. 
-        ur = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
-        un = md.results.EsaSolution.EsaNmotion*1000; % [mm] 
-        ue = md.results.EsaSolution.EsaEmotion*1000; % [mm] 
+        % 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] 
+        sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',...
+                'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; 
 
         res = 1.0; % degree 
@@ -128 +131 @@
         sol_grid = zeros(size(lat_grid)); 
 
-        sol = ue; 
+        for kk=1:4 
+                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;
+                                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 % }}}
 
-        % Make a interpolation object
-        F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
-        F.Method = 'linear';
-        F.ExtrapolationMethod = 'linear'; 
+                % 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; 
+                % 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 polar unphysical 0s to Nan 
-        sol_grid(lat_grid>85 & sol_grid==0) =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');
 
-        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; 
-        %caxis([-0.3 0.3])
-        hold on
-        plot(long,lat,'k');
-        %geoshow(lat,long,'DisplayType','polygon','FaceColor','white');         
-        hold off; 
-        c1=colorbar;
-        colormap(jet); 
-        xlim([-180 180]); 
-        ylim([-90 90]); 
-        grid on; 
-        title(['Average change in relative sea-level [mm/yr]']);
-        set(gcf,'color','w');
-
-        %export_fig('Fig5.pdf'); 
+                %export_fig('Fig5.pdf'); 
+        end
 
 end % }}} 
 
-if any(steps==6) % {{{ Transient 
-        disp('   Step 6: Transient run');
-
-end % }}} 
-