Changeset 22809

Timestamp:

05/29/18 14:07:24 (7 years ago)

Author:

adhikari

Message:

CHG: functions are located in a seperate directory etc

Location:

issm/trunk-jpl/examples

Files:

• EsaGRACE/runme.m (modified) (1 diff)
• Functions (added)
• Functions/grace.m (added)
• Functions/wahr.m (added)
• SlrGRACE/runme.m (modified) (10 diffs)

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

@@ -1 +1 @@
 
 clear all;
-steps=[0]; % [1:5]; 
+steps=[0:5]; % [1:5]; 
 
 if any(steps==0) % Download GRACE land_mass data {{{

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

@@ -1 +1 @@
 
 clear all;
-steps=[5]; % [1:6]; 
+addpath('../Data','../Functions');
+
+steps=[7]; % [0:7]; 
 
 if any(steps==0) % Download GRACE land_mass data {{{
@@ -11 +13 @@
         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. 
@@ -32 +36 @@
         md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. 
 
-        for i=1:numrefine,
+        for i=1:numrefine; % refine mesh... {{{
 
                 %figure out mask: 
@@ -66 +70 @@
                 %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
+        end % }}} 
 
         %figure out mask: 
@@ -156 +160 @@
         md.slr.rotation=1;
 
-        % Request outputs 
-%       md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'};
-        
         % Cluster info 
         md.cluster=generic('name',oshostname(),'np',3); 
@@ -206 +207 @@
                 F.Method = 'linear';
                 F.ExtrapolationMethod = 'linear'; 
-
+                
                 % Do the interpolation to get gridded solutions... 
                 sol_grid = F(lat_grid, lon_grid);
@@ -214 +215 @@
                 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; 
+                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 % }}}
                 plot(long,lat,'k'); hold off; 
+                % define colormap, caxis, xlim etc {{{
                 c1=colorbar;
-                colormap(jet); 
+                colormap('haxby'); 
+                caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))])
                 xlim([-180 180]); 
                 ylim([-90 90]); 
+                % }}} 
                 grid on; 
                 title(sol_name(kk)); 
@@ -237 +245 @@
 
         % read GRACE data during the chosen time period << steps=2 >>
-        disp('Projecting grace loads onto the mesh...');
+        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)); 
@@ -261 +269 @@
         md.settings.output_frequency=1;
 
-        % Request outputs 
-        md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'};
-        
         % Cluster info 
         md.cluster=generic('name',oshostname(),'np',3); 
@@ -271 +276 @@
         md=solve(md,'tr');
 
-end % }}} 
-
+        save ./Models/SlrGRACE.Transient md; 
+
+end % }}} 
+
+if any(steps==7) % Plot transient {{{ 
+        disp('   Step 7: Plot transient');
+        md = loadmodel('./Models/SlrGRACE.Transient');
+
+        time = md.slr.deltathickness(end,:); % year 
+
+        % loads and solutions. 
+        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
+        sol_name = {'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; 
+        movie_name = {'movie_loads.avi','movie_fingerprints.avi'}; 
+
+        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)); 
+
+        for kk=1:2
+                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); 
+                                temp2(jj,:)=mean(sol(pos,:)); 
+                        end
+                        sol=temp2; 
+                end % }}}
+
+                vidObj = VideoWriter(movie_name{kk}); 
+                vidObj.FrameRate=2; % frames per second 
+                open(vidObj);
+
+                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 
+
+                        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;
+                        % 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 % }}}
+                        plot(long,lat,'k'); hold off; 
+                        % define colormap, caxis, xlim etc {{{
+                        c1=colorbar;
+                        colormap('haxby'); 
+                        caxis([-floor(min(abs(min(min(sol))),abs(max(max(sol))))) floor(min(abs(min(min(sol))),abs(max(max(sol)))))])
+                        xlim([-180 180]); 
+                        ylim([-90 90]); 
+                        % }}} 
+                        grid on; 
+                        title(sol_name(kk)); 
+                        set(gcf,'color','w');
+                        writeVideo(vidObj,getframe(gcf)); 
+                        close % close current figure 
+                end
+
+                % Close the file.
+                close(vidObj);
+        end
+        !open *.avi; 
+        
+        % plot GMSL time series 
+        plot(time,gmsl,'-*','linewidth',3); grid on; 
+        xlabel('# month'); 
+        ylabel('GMSL [mm/yr]');
+        set(gcf,'color','w');
+
+        %export_fig('Fig7.pdf');
+
+end % }}} 
+