Changeset 26142


Ignore:
Timestamp:
03/23/21 16:59:53 (4 years ago)
Author:
jdquinn
Message:

CHG: Sean’s updates

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r25986 r26142  
    44steps=[7];
    55
    6 if any(steps==1) 
     6if any(steps==1)
    77        disp('   Step 1: Global mesh creation');
    88
    9         numrefine=1;
    10         resolution=300*1e3;                     % inital resolution [m]
     9        numrefine=5;
     10        resolution=150*1e3;                     % inital resolution [m]
    1111        radius = 6.371012*10^6;         % mean radius of Earth, m
    1212        mindistance_coast=150*1e3;      % coastal resolution [m]
    … …  
    2323                distance=zeros(md.mesh.numberofvertices,1);
    2424
    25                 pos=find(~ocean_mask);  coaste.lat=md.mesh.lat(pos);    coaste.long=md.mesh.long(pos);
    26                 pos=find(ocean_mask);   coasto.lat=md.mesh.lat(pos);    coasto.long=md.mesh.long(pos);
    27 
     25                pos=find(~ocean_mask);  coaste.lat=md.mesh.lat(pos);    coaste.long=md.mesh.long(pos);          pos=find(ocean_mask);   coasto.lat=md.mesh.lat(pos);    coasto.long=md.mesh.long(pos);
    2826                for j=1:md.mesh.numberofvertices
    2927                        phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
    … …  
    4038                end
    4139                pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
    42                
     40
    4341                pos2=find(ocean_mask~=1 & distance>mindistance_land);
    4442                distance(pos2)=mindistance_land;
    … …  
    4947
    5048        ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
    51         pos = find(ocean_mask==0);
    52         md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
    53         md.mask.ocean_levelset(pos)=1;
    54 
     49        pos = find(ocean_mask==0);      md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);       md.mask.ocean_levelset(pos)=1;
    5550        save ./Models/SlrGRACE_Mesh md;
    5651
    5752        plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
    58 end
    59 
    60 if any(steps==2)
     53end
     54if any(steps==2)
    6155        disp('   Step 2: Define loads in meters of ice height equivalent');
    6256        md = loadmodel('./Models/SlrGRACE_Mesh');
    … …  
    7266
    7367        plotmodel (md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]');
    74 end
    75 
    76 if any(steps==3)
     68end
     69if any(steps==3)
    7770        disp('   Step 3: Parameterization');
    7871        md = loadmodel('./Models/SlrGRACE_Loads');
    … …  
    8073        md.solidearth.lovenumbers = lovenumbers('maxdeg',10000);
    8174
    82         md.mask.ice_levelset=-md.mask.ocean_levelset;
    83 
     75        md.mask.ice_levelset=-md.mask.ocean_levelset;
    8476        md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    8577        md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    … …  
    8981        md.solidearth.settings.ocean_area_scaling=1;
    9082
    91         % arbitary to pass consistency check.
    92         md.geometry.bed=-ones(md.mesh.numberofvertices,1);
     83        % arbitary to pass consistency check.   md.geometry.bed=-ones(md.mesh.numberofvertices,1);
    9384        md.geometry.surface=ones(md.mesh.numberofvertices,1);
    94         md.geometry.base=md.geometry.bed;
    95         md.geometry.thickness=md.geometry.surface-md.geometry.base;
    96 
     85        md.geometry.base=md.geometry.bed;       md.geometry.thickness=md.geometry.surface-md.geometry.base;
    9786        md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    9887        md.materials.rheology_B=paterson(md.initialization.temperature);
    … …  
    10291
    10392        save ./Models/SlrGRACE_Parameterization md;
    104 end
    105 
    106 if any(steps==4)
     93end
     94if any(steps==4)
    10795        disp('   Step 4: Solve Slr solver');
    10896        md = loadmodel('./Models/SlrGRACE_Parameterization');
    … …  
    120108
    121109        save ./Models/SlrGRACE_Solution md;
    122 end
    123 
    124 if any(steps==5)
     110end
     111if any(steps==5)
    125112        disp('   Step 5: Plot solutions');
    126113        md = loadmodel('./Models/SlrGRACE_Solution');
    127114
    128         sol1 = md.solidearth.surfaceload.icethicknesschange*100;                                                % [cm]
     115        sol1 = md.solidearth.surfaceload.icethicknesschange*100; % [cm]
    129116        sol2 = md.results.SealevelriseSolution.SealevelRSL*1000;        % [mm]
    130117
    … …  
    155142                F.Method = 'linear';
    156143                F.ExtrapolationMethod = 'linear';
    157                
     144
    158145                sol_grid = F(lat_grid, lon_grid);
    159146                sol_grid(isnan(sol_grid))=0;
    … …  
    162149                set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
    163150                figure1=figure('Position', [100, 100, 1000, 500]);
    164                 gcf; load coast; cla;
     151                gcf; load coastlines; cla;
    165152                pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
    166153                if (kk==1)
    167                         geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
     154                        geoshow(flipud(coastlat),flipud(coastlon),'DisplayType','polygon','FaceColor','white');
    168155                else
    169                         geoshow(lat,long,'DisplayType','polygon','FaceColor','white');
    170                 end
    171                 plot(long,lat,'k'); hold off;
     156                        geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor',[.5 1 .5]);
     157                end
     158                plot(coastlon, coastlat,'k'); hold off;
    172159                c1=colorbar;
    173160                colormap('haxby');
    … …  
    180167                %export_fig(fig_name{kk});
    181168        end
    182 end
    183 
    184 if any(steps==6)
     169end
     170if any(steps==6)
    185171        disp('   Step 6: Transient run');
    186172        md = loadmodel('./Models/SlrGRACE_Parameterization');
    … …  
    199185        md.transient.isstressbalance=0;
    200186        md.transient.isthermal=0;
    201         md.transient.isgia=1;
    202         md.transient.isslr=1;
     187        md.transient.isgia=1;   md.transient.isslr=1;
    203188
    204189        md.timestepping.start_time=-1;
    … …  
    212197
    213198        md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
    214        
     199
    215200        md=solve(md,'tr');
    216201
    217202        save ./Models/SlrGRACE_Transient md;
    218 end
    219 
    220 if any(steps==7)
     203end
     204if any(steps==7)
    221205        disp('   Step 7: Plot transient');
    222206        md = loadmodel('./Models/SlrGRACE_Transient');
    … …  
    230214        end
    231215        sol_name = {'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
    232         movie_name = {'Movie_dH.avi','Movie_slr.avi'};
     216        movie_name = {'Movie_dH.mp4','Movie_slr.mp4'};
    233217
    234218        res = 1.0;
    … …  
    254238                vidObj = VideoWriter(movie_name{kk});
    255239                vidObj.FrameRate=2; % frames per second
     240                vifObj.FileFormat='mp4';
    256241                open(vidObj);
    257242
    258243                for jj=1:length(time)
     244                        fprintf('creating frame %d...', jj);
    259245                        F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj));
    260246                        F.Method = 'linear';
    261247                        F.ExtrapolationMethod = 'linear';
    262                        
     248
    263249                        sol_grid = F(lat_grid, lon_grid);
    264250                        sol_grid(isnan(sol_grid))=0;
    … …  
    267253                        set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
    268254                        figure1=figure('Position', [100, 100, 1000, 500]);
    269                         gcf; load coast; cla;
     255                        gcf; load coastlines; cla;
    270256                        pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
    271257                        if (kk==1)
    272                                 geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
     258                                geoshow(flipud(coastlat),flipud(coastlon),'DisplayType','polygon','FaceColor','white');
    273259                        else
    274                                 geoshow(lat,long,'DisplayType','polygon','FaceColor','white');
    275                         end
    276                         plot(long,lat,'k'); hold off;
     260                                geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor','white');
     261                        end
     262                        plot(coastlon,coastlat,'k'); hold off;
    277263                        c1=colorbar;
    278264                        colormap('haxby');
    … …  
    285271                        writeVideo(vidObj,getframe(gcf));
    286272                        close
    287                 end
     273                        fprintf('done\n');
     274                end
     275                disp('closing vidObj...');
    288276                close(vidObj);
    289277        end
    290         !open *.avi;
     278        !open *.mp4;
    291279
    292280        % plot GMSL time series
    … …  
    296284        set(gcf,'color','w');
    297285        %export_fig('Fig7.pdf');
    298 end 
     286end
Note: See TracChangeset for help on using the changeset viewer.