Changeset 22818


Ignore:
Timestamp:
05/31/18 09:58:48 (7 years ago)
Author:
adhikari
Message:

CHG: several minor changes

Location:
issm/trunk-jpl/examples
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r22815 r22818  
    33addpath('../Data','../Functions');
    44
    5 steps=[0]; % [0:5];
     5steps=[1]; % [1:5]
    66
    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 {{{
     7if any(steps==1)
    248        disp('   Step 1: Global mesh creation');
    259
    26         resolution=300; % km. uniform meshing...
    27         radius = 6.371012*10^3; % mean radius of Earth, km
     10        resolution=300;                 % [km]
     11        radius = 6.371012*10^3; % [km]
    2812
    29         %mesh earth:
    3013        md=model;
    31         md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset
     14        md.mask=maskpsl();
    3215        md.mesh=gmshplanet('radius',radius,'resolution',resolution);
    3316
    34         %figure out mask:
    3517        md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
    3618
    37         save ./Models/EsaGRACE.Mesh md;
     19        save ./Models/EsaGRACE_Mesh md;
    3820       
    3921        plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
    40         %export_fig('Fig1.pdf');
    4122
    42 end % }}}
     23end
    4324
    44 if any(steps==2) % Define loads {{{
     25if any(steps==2)
    4526        disp('   Step 2: Define loads in meters of ice height equivalent');
    46         md = loadmodel('./Models/EsaGRACE.Mesh');
     27        md = loadmodel('./Models/EsaGRACE_Mesh');
    4728
    48         % define time interval for analysis
    4929        year_month = 2007+15/365;
    5030        time_range = [year_month year_month];
    5131       
    52         % map GRACE water load on to the mesh for the seleted month(s)
    5332        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
    5433       
    5534        md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent
    5635
    57         save ./Models/EsaGRACE.Loads md;
     36        save ./Models/EsaGRACE_Loads md;
    5837       
    5938        plotmodel (md,'data',md.esa.deltathickness,...
    6039                'view',[90 -90],'caxis',[-.1 .1],...
    6140                'title','Ice height equivalent [m]');
    62         %export_fig('Fig2.pdf');
    6341
    64 end % }}}
     42end
    6543
    66 if any(steps==3) % Parameterization {{{
     44if any(steps==3)
    6745        disp('   Step 3: Parameterization');
    68         md = loadmodel('./Models/EsaGRACE.Loads');
     46        md = loadmodel('./Models/EsaGRACE_Loads');
    6947
    70         % Love numbers and reference frame: CF or CM (choose one!) 
    71         nlove=10001;    % up to 10,000 degree
     48        nlove=10001;   
    7249        md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[];
    7350        md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[];
    7451
    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);
    7753        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
    7854        pos=find(md.esa.deltathickness~=0);
    79         md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads
     55        md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
    8056        md.mask.land_levelset = 1-md.mask.ocean_levelset;
    8157
    82         %% IGNORE BUT DO NOT DELETE %% {{{
    83         % Geometry: Important only when you want to couple with Ice Flow Model
    8458        di=md.materials.rho_ice/md.materials.rho_water;
    8559        md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    … …  
    8761        md.geometry.base=md.geometry.surface-md.geometry.thickness;
    8862        md.geometry.bed=md.geometry.base;
    89         % Materials:
     63       
    9064        md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    9165        md.materials.rheology_B=paterson(md.initialization.temperature);
    9266        md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    93         % Miscellaneous:
     67       
    9468        md.miscellaneous.name='EsaGRACE';
    95         %% IGNORE BUT DO NOT DELETE %% }}} 
    9669       
    97         save ./Models/EsaGRACE.Parameterization md;
     70        save ./Models/EsaGRACE_Parameterization md;
    9871
    99 end % }}}
     72end
    10073
    101 if any(steps==4) % Solve {{{
     74if any(steps==4)
    10275        disp('   Step 4: Solve Esa solver');
    103         md = loadmodel('./Models/EsaGRACE.Parameterization');
     76        md = loadmodel('./Models/EsaGRACE_Parameterization');
    10477
    105         % Request outputs
    10678        md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'};
    10779       
    108         % Cluster info
    10980        md.cluster=generic('name',oshostname(),'np',3);
    11081        md.verbose=verbose('111111111');
    11182
    112         % Solve
    11383        md=solve(md,'Esa');
    11484
    115         save ./Models/EsaGRACE.Solution md;
     85        save ./Models/EsaGRACE_Solution md;
    11686
    117 end % }}}
     87end
    11888
    119 if any(steps==5) % Plot solutions {{{
     89if any(steps==5)
    12090        disp('   Step 5: Plot solutions');
    121         md = loadmodel('./Models/EsaGRACE.Solution');
     91        md = loadmodel('./Models/EsaGRACE_Solution');
    12292
    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]
    12897
    12998        sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',...
    … …  
    131100        fig_name={'Fig_dH.pdf','Fig_vert.pdf','Fig_horzNS.pdf','Fig_horzEW.pdf'};
    132101
    133         res = 1.0; % degree
     102        res = 1.0; % [degree]
    134103
    135         % Make a grid of lats and lons, based on the min and max of the original vectors
    136104        [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
    137105        sol_grid = zeros(size(lat_grid));
    … …  
    140108                sol=eval(sprintf('sol%d',kk));
    141109       
    142                 % if data are on elements, map those on to the vertices {{{
    143110                if length(sol)==md.mesh.numberofelements
    144                         % map on to the vertices
    145111                        for jj=1:md.mesh.numberofelements
    146112                                ii=(jj-1)*3;
    … …  
    152118                        end
    153119                        sol=temp';
    154                 end % }}}
     120                end
    155121
    156                 % Make a interpolation object
    157122                F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol);
    158123                F.Method = 'linear';
    159124                F.ExtrapolationMethod = 'linear';
    160125
    161                 % Do the interpolation to get gridded solutions...
    162126                sol_grid = F(lat_grid, lon_grid);
    163127                sol_grid(isnan(sol_grid))=0;
    164                 sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan
     128                sol_grid(lat_grid>85 & sol_grid==0)=NaN;
    165129
    166130                set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
    … …  
    168132                gcf; load coast; cla;
    169133                pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
    170                 % mask out oceans {{{
    171134                if (kk==1)
    172135                        geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
    173                 end % }}}
     136                end
    174137                plot(long,lat,'k'); hold off;
    175                 % colormap, xlim etc {{{
    176138                c1=colorbar;
    177139                colormap('haxby');
    … …  
    179141                xlim([-180 180]);
    180142                ylim([-90 90]);
    181                 % }}}
    182143                grid on;
    183144                title(sol_name(kk));
    184145                set(gcf,'color','w');
    185                
    186146                %export_fig(fig_name{kk});
    187147        end
    188148
    189 end % }}}
     149end
    190150

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

    r22815 r22818  
    33addpath('../Functions');
    44
    5 steps=[0]; % [0:6];
     5steps=[0]; % [0:6]
    66
    7 if any(steps==0) % Simple mesh creation {{{
     7if any(steps==0)
    88        disp('   Step 0: Mesh creation');
    99
    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]
    1211
    13         save ./Models/EsaWahr.Mesh md;
     12        save ./Models/EsaWahr_Mesh md;
    1413       
    1514        plotmodel(md,'data','mesh');
    16         %export_fig('Fig0.pdf');
    1715
    18 end % }}}
     16end
    1917
    20 if any(steps==1) % {{{ Anisotropic mesh creation 
     18if any(steps==1)
    2119        disp('   Step 1: Anisotropic mesh creation');
    2220
    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);
    2622
    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);       
    3025        field = abs(rad_dist-disc_radius);
    31         %plotmodel(md,'data',field);
    3226
    33         md = bamg(md,'field',field,'err',50,'hmax',10000); % error in field in m
     27        md = bamg(md,'field',field,'err',50,'hmax',10000);
    3428
    35         save ./Models/EsaWahr.Mesh md;
     29        save ./Models/EsaWahr_Mesh md;
    3630       
    3731        plotmodel (md,'data','mesh');
    38         %export_fig('Fig1.pdf');
    3932
    40 end % }}}
     33end
    4134
    42 if any(steps==2) % Define loads {{{
     35if any(steps==2)
    4336        disp('   Step 2: Define loads');
    44         md = loadmodel('./Models/EsaWahr.Mesh');
     37        md = loadmodel('./Models/EsaWahr_Mesh');
    4538
    46         % primer
    4739        rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice;
    4840
    49         % elemental centroids
    5041        index=md.mesh.elements;         
    5142        x_cent=mean(md.mesh.x(index),2);
    5243        y_cent=mean(md.mesh.y(index),2);
    5344
    54         % Define a 20-km radius disc and unload it by 1 meter water height equivalent uniformly 
    5545        md.esa.deltathickness = zeros(md.mesh.numberofelements,1);
    56         disc_radius=20; % km
    57         rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000;        % radial distance in km
    58         md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i;    % 1 m water withdrawl
     46        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;
    5949
    60         save ./Models/EsaWahr.Loads md;
     50        save ./Models/EsaWahr_Loads md;
    6151       
    6252        plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]');
    63         %export_fig('Fig2.pdf');
    6453
    65 end % }}}
     54end
    6655
    67 if any(steps==3) % Parameterization {{{
     56if any(steps==3)
    6857        disp('   Step 3: Parameterization');
    69         md = loadmodel('./Models/EsaWahr.Loads');
     58        md = loadmodel('./Models/EsaWahr_Loads');
    7059
    71         % Love numbers and reference frame: CF or CM (choose one!) 
    72         nlove=10001;    % up to 10,000 degree
     60        nlove=10001;
    7361        md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[];
    7462        md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[];
    7563
    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);
    7966
    80         %% IGNORE BUT DO NOT DELETE %% {{{
    81         % Geometry: Important only when you want to couple with Ice Flow Model
    8267        di=md.materials.rho_ice/md.materials.rho_water;
    8368        md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    … …  
    8570        md.geometry.base=md.geometry.surface-md.geometry.thickness;
    8671        md.geometry.bed=md.geometry.base;
    87         % Materials:
     72       
    8873        md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    8974        md.materials.rheology_B=paterson(md.initialization.temperature);
    9075        md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    91         % Miscellaneous:
     76       
    9277        md.miscellaneous.name='EsaWahr';
    93         %% IGNORE BUT DO NOT DELETE %% }}} 
    9478       
    95         save ./Models/EsaWahr.Parameterization md;
     79        save ./Models/EsaWahr_Parameterization md;
    9680
    97 end % }}}
     81end
    9882
    99 if any(steps==4) % Solve {{{
     83if any(steps==4)
    10084        disp('   Step 4: Solve Esa solver');
    101         md = loadmodel('./Models/EsaWahr.Parameterization');
     85        md = loadmodel('./Models/EsaWahr_Parameterization');
    10286
    103         % Request outputs
    10487        md.esa.requested_outputs = {'EsaUmotion','EsaXmotion','EsaYmotion'};
    10588       
    106         % Cluster info
    10789        md.cluster=generic('name',oshostname(),'np',3);
    10890        md.verbose=verbose('111111111');
    10991
    110         % Solve
    11192        md=solve(md,'Esa');
    11293
    113         save ./Models/EsaWahr.Solution md;
     94        save ./Models/EsaWahr_Solution md;
    11495
    115 end % }}}
     96end
    11697
    117 if any(steps==5) % Plot solutions {{{
     98if any(steps==5)
    11899        disp('   Step 5: Plot solutions');
    119         md = loadmodel('./Models/EsaWahr.Solution');
     100        md = loadmodel('./Models/EsaWahr_Solution');
    120101
    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] 
    126106
    127         % plot
    128107        set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6);
    129108        figure('Position', [100, 100, 800, 600]);
    … …  
    147126                'axispos',[0.505 0.02 0.47 0.47],...
    148127                'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]');
    149 
    150128        %export_fig('Fig5.pdf');
    151129
    152 end % }}}
     130end
    153131
    154 if any(steps==6) % Compare results against semi-analytic solutions {{{
     132if any(steps==6)
    155133        disp('   Step 6: Compare results against Wahr semi-analytic solutions');
    156         md = loadmodel('./Models/EsaWahr.Solution');
     134        md = loadmodel('./Models/EsaWahr_Solution');
    157135
    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]
    166142        yi=zeros(1,length(xi));
    167143        vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear');
    … …  
    169145
    170146        % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1)
    171         disc_radius = 20*1000; % 20 km
     147        disc_radius = 20*1000; % [m]
    172148        [vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l);
    173149
    174         % plot
    175150        set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6);
    176151        figure1=figure('Position', [100, 100, 700, 400]);
    … …  
    186161                h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1);
    187162                h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1);
    188                 % ISSM
     163                % ISSM soln
    189164                h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1);
    190165                h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1);
    191                 % first box
    192166                ag1 = gca;
    193167                leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)');
    194168                set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16);
    195                 %
    196         set(gcf,'color','w');
    197        
     169                set(gcf,'color','w');
    198170        %export_fig('Fig6.pdf');
    199171
    200 end % }}}
     172end
    201173

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

    r22815 r22818  
    22clear all;
    33
    4 steps=[1]; % [1:5];
     4steps=[1]; % [1:5]
    55
    6 if any(steps==1) % Global mesh creation {{{
     6if any(steps==1)
    77        disp('   Step 1: Global mesh creation');
    88
    99        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, 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]
     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]
    1515
    16         %mesh earth:
    1716        md=model;
    18         md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset
    19         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]
    2019
    2120        for i=1:numrefine,
    2221
    23                 %figure out mask:
    2422                md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
    2523
    26                 %figure out distance to the coastline, in lat,long (not x,y,z):
    2724                distance=zeros(md.mesh.numberofvertices,1);
    2825
    … …  
    3128
    3229                for j=1:md.mesh.numberofvertices
    33                         %figure out nearest coastline (using the great circle distance)
    3430                        phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
    3531                        if md.mask.ocean_levelset(j),
    … …  
    4642                pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
    4743               
    48                 % refine on the continents
    4944                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
    5045                distance(pos2)=mindistance_land;
    5146
    52                 dist=min(maxdistance,distance); % max size 1000 km
    53                 %use distance to the coastline to refine mesh:
     47                dist=min(maxdistance,distance);
    5448                md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist);
     49
    5550        end
    5651
    57         %figure out mask:
    5852        md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
    5953
    60         save ./Models/SlrFarrell.Mesh md;
     54        save ./Models/SlrFarrell_Mesh md;
    6155       
    6256        plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
    63         %export_fig('Fig1.pdf');
    6457
    65 end % }}}
     58end
    6659
    67 if any(steps==2) % Define source {{{
     60if any(steps==2)
    6861        disp('   Step 2: Define source as in Farrell, 1972, Figure 1');
    69         md = loadmodel('./Models/SlrFarrell.Mesh');
     62        md = loadmodel('./Models/SlrFarrell_Mesh');
    7063
    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
    7365       
    7466        md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
    7567        md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
    7668
    77         save ./Models/SlrFarrell.Loads md;
     69        save ./Models/SlrFarrell_Loads md;
    7870       
    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]');
    8272
    83 end % }}}
     73end
    8474
    85 if any(steps==3) % Parameterization {{{
     75if any(steps==3)
    8676        disp('   Step 3: Parameterization');
    87         md = loadmodel('./Models/SlrFarrell.Loads');
     77        md = loadmodel('./Models/SlrFarrell_Loads');
    8878
    89         % Love numbers and reference frame: CF or CM (choose one!) 
    90         nlove=10001;    % up to 10,000 degree
     79        nlove=10001;   
    9180        md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[];
    9281        md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[];
    9382        md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[];
    9483
    95         % Mask: for computational efficiency only those elements that have loads are convolved!
    9684        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);
    10087        pos=find(md.mesh.lat <-80);
    10188        md.mask.ice_levelset(pos(1))=-1; % ice yes! 
    10289        md.mask.groundedice_levelset(pos(1))=1; % ice grounded! 
    10390
    104         %% IGNORE BUT DO NOT DELETE %% {{{
    105         % Geometry: Important only when you want to couple with Ice Flow Model
    10691        di=md.materials.rho_ice/md.materials.rho_water;
    10792        md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    … …  
    10994        md.geometry.base=md.geometry.surface-md.geometry.thickness;
    11095        md.geometry.bed=md.geometry.base;
    111         % Materials:
     96       
    11297        md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    11398        md.materials.rheology_B=paterson(md.initialization.temperature);
    11499        md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    115         % Miscellaneous:
     100       
    116101        md.miscellaneous.name='SlrFarrell';
    117         %% IGNORE BUT DO NOT DELETE %% }}} 
    118102       
    119         save ./Models/SlrFarrell.Parameterization md;
     103        save ./Models/SlrFarrell_Parameterization md;
    120104
    121 end % }}}
     105end
    122106
    123 if any(steps==4) % Solve {{{
     107if any(steps==4)
    124108        disp('   Step 4: Solve Slr solver');
    125         md = loadmodel('./Models/SlrFarrell.Parameterization');
     109        md = loadmodel('./Models/SlrFarrell_Parameterization');
    126110
    127         % Cluster info
    128111        md.cluster=generic('name',oshostname(),'np',3);
    129112        md.verbose=verbose('111111111');
    130113
    131         % Choose different convergence threshold. [10% 1% 0.1%] to match Farrell 3 panels in Fig. 1
    132114        md.slr.reltol = 0.1/100; % per cent change in solution
    133115
    134         % Solve
    135116        md=solve(md,'Slr');
    136117
    137         save ./Models/SlrFarrell.Solution md;
     118        save ./Models/SlrFarrell_Solution md;
    138119
    139 end % }}}
     120end
    140121
    141 if any(steps==5) % Plot solutions {{{
     122if any(steps==5)
    142123        disp('   Step 5: Plot solutions');
    143         md = loadmodel('./Models/SlrFarrell.Solution');
     124        md = loadmodel('./Models/SlrFarrell_Solution');
    144125
    145         % solutions.
    146126        sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) 
    147127
    148         res = 1; % degree
     128        res = 1; % [degree]
    149129
    150         % Make a grid of lats and lons, based on the min and max of the original vectors
    151130        [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
    152131        sol_grid = zeros(size(lat_grid));
    153132
    154         % Make a interpolation object
    155133        F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol);
    156         F.Method = 'natural'; % for smooth contour
    157         F.ExtrapolationMethod = 'none';
     134        F.Method = 'linear';
     135        F.ExtrapolationMethod = 'linear';
    158136       
    159         % Do the interpolation to get gridded solutions...
    160137        sol_grid = F(lat_grid, lon_grid);
    161138        sol_grid(isnan(sol_grid))=0;
    162         sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan
     139        sol_grid(lat_grid>85 & sol_grid==0) =NaN;
    163140
    164141        set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
    … …  
    170147        geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]);
    171148        plot(long,lat,'k'); hold off;
    172         % define colormap, caxis, xlim etc {{{
    173149        c1=colorbar;
    174150        colormap(flipud(haxby));
    … …  
    176152        xlim([-170 170]);
    177153        ylim([-85 85]);
    178         % }}}
    179154        grid on;
    180155        title('Relative sea-level [% of GMSL]');
    181156        set(gcf,'color','w');
    182 
    183157        %export_fig('Fig5.pdf');
    184158
    185 end % }}}
     159end
    186160

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

    r22812 r22818  
    33addpath('../Data','../Functions');
    44
    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 {{{
     5steps=[1]; % [1:7]
     6
     7if any(steps==1)
    248        disp('   Step 1: Global mesh creation');
    259
    2610        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
    3417        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
    4123                md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
    4224
    43                 %figure out distance to the coastline, in lat,long (not x,y,z):
    4425                distance=zeros(md.mesh.numberofvertices,1);
    4526
    … …  
    4829
    4930                for j=1:md.mesh.numberofvertices
    50                         %figure out nearest coastline (using the great circle distance)
    5131                        phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
    5232                        if md.mask.ocean_levelset(j),
    … …  
    6343                pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
    6444               
    65                 % refine on the continents
    6645                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
    6746                distance(pos2)=mindistance_land;
    6847
    69                 dist=min(maxdistance,distance); % max size 1000 km
    70                 %use distance to the coastline to refine mesh:
     48                dist=min(maxdistance,distance);
    7149                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
    7553        md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
    7654
    77         save ./Models/SlrGRACE.Mesh md;
     55        save ./Models/SlrGRACE_Mesh md;
    7856       
    7957        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
     59end
     60
     61if any(steps==2)
    8562        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
    8965        year_month = 2007+15/365;
    9066        time_range = [year_month year_month];
    9167       
    92         % map GRACE water load on to the mesh for the seleted month(s)
    9368        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
    9469       
    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
     76end
     77
     78if any(steps==3)
    10779        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;
    11283        md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[];
    11384        md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[];
    11485        md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[];
    11586
    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);
    11888        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
    11989        pos=find(md.slr.deltathickness~=0);
    120         md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads
     90        md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
    12191        md.mask.land_levelset = 1-md.mask.ocean_levelset;
    12292
    123         % initialize
    12493        md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
    125         % steric rate
    12694        md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
    127         % solution is scaled according to "exact" area of the oceans
    12895        md.slr.ocean_area_scaling = 1;
    12996
    130         % convergence criteria
    131         %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 Model
    13697        di=md.materials.rho_ice/md.materials.rho_water;
    13798        md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    … …  
    139100        md.geometry.base=md.geometry.surface-md.geometry.thickness;
    140101        md.geometry.bed=md.geometry.base;
    141         % Materials:
     102       
    142103        md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    143104        md.materials.rheology_B=paterson(md.initialization.temperature);
    144105        md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    145         % Miscellaneous:
     106       
    146107        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
     111end
     112
     113if any(steps==4)
    154114        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
    158117        md.slr.rigid=1;
    159118        md.slr.elastic=1;
    160119        md.slr.rotation=1;
    161120
    162         % Cluster info
    163121        md.cluster=generic('name',oshostname(),'np',3);
    164122        md.verbose=verbose('111111111');
    165123
    166         % Solve
    167124        md=solve(md,'Slr');
    168125
    169         save ./Models/SlrGRACE.Solution md;
    170 
    171 end % }}}
    172 
    173 if any(steps==5) % Plot solutions {{{
     126        save ./Models/SlrGRACE_Solution md;
     127
     128end
     129
     130if any(steps==5)
    174131        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]
    179135        sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm]
    180136       
    … …  
    182138        fig_name={'Fig_dH.pdf','Fig_slr.pdf'};
    183139
    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
    187142        [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
    188143        sol_grid = zeros(size(lat_grid));
    … …  
    191146                sol=eval(sprintf('sol%d',kk));
    192147       
    193                 % if data are on elements, map those on to the vertices {{{
    194148                if length(sol)==md.mesh.numberofelements
    195                         % map on to the vertices
    196149                        for jj=1:md.mesh.numberofelements
    197150                                ii=(jj-1)*3;
    … …  
    203156                        end
    204157                        sol=temp';
    205                 end % }}}
    206 
    207                 % Make a interpolation object
     158                end
     159
    208160                F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol);
    209161                F.Method = 'linear';
    210162                F.ExtrapolationMethod = 'linear';
    211163               
    212                 % Do the interpolation to get gridded solutions...
    213164                sol_grid = F(lat_grid, lon_grid);
    214165                sol_grid(isnan(sol_grid))=0;
    215                 sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan
     166                sol_grid(lat_grid>85 & sol_grid==0) = NaN;
    216167
    217168                set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
    … …  
    219170                gcf; load coast; cla;
    220171                pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
    221                 % mask out land or oceans {{{
    222172                if (kk==1)
    223173                        geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
    224174                else
    225175                        geoshow(lat,long,'DisplayType','polygon','FaceColor','white');
    226                 end % }}}
     176                end
    227177                plot(long,lat,'k'); hold off;
    228                 % define colormap, caxis, xlim etc {{{
    229178                c1=colorbar;
    230179                colormap('haxby');
    … …  
    232181                xlim([-180 180]);
    233182                ylim([-90 90]);
    234                 % }}}
    235183                grid on;
    236184                title(sol_name(kk));
    237185                set(gcf,'color','w');
    238 
    239186                %export_fig(fig_name{kk});
    240187        end
    241188
    242 end % }}}
    243 
    244 if any(steps==6) % Transient {{{
     189end
     190
     191if any(steps==6)
    245192        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
    249195        disp('Projecting  loads onto the mesh...');
    250196        time_range = 2007 + [15 350]/365;
    251197        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
    252198
    253         % define ice load
    254199        [ne,nt]=size(water_load);
    255200   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 equivalent
     201        md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
    257202        md.slr.deltathickness(ne+1,:)=[1:nt]; % times
    258203       
    259         % define transient run
    260204        md.transient.issmb=0;
    261205        md.transient.ismasstransport=0;
    … …  
    264208        md.transient.isslr=1;
    265209       
    266         % time stepping, output frequency etc.
    267210        md.timestepping.start_time=-1;
    268211   md.timestepping.final_time=nt;
    … …  
    271214        md.settings.output_frequency=1;
    272215
    273         % Cluster info
    274216        md.cluster=generic('name',oshostname(),'np',3);
    275217        md.verbose=verbose('111111111');
    276218
    277         % Solve
    278219        md=solve(md,'tr');
    279220
    280         save ./Models/SlrGRACE.Transient md;
    281 
    282 end % }}}
    283 
    284 if any(steps==7) % Plot transient {{{
     221        save ./Models/SlrGRACE_Transient md;
     222
     223end
     224
     225if any(steps==7)
    285226        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
    291231        for tt=1:length(time)
    292232                gmsl(tt) = md.results.TransientSolution(tt).SealevelEustatic*1000;      % GMSL rate mm/yr
    293233                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!
    295234           sol2(:,tt) = md.results.TransientSolution(tt+1).Sealevel*1000;                       % mm/yr
    296235        end
    … …  
    298237        movie_name = {'Movie_dH.avi','Movie_slr.avi'};
    299238
    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
    303241        [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
    304242        sol_grid = zeros(size(lat_grid));
    … …  
    307245                sol=eval(sprintf('sol%d',kk));
    308246       
    309                 % if data are on elements, map those on to the vertices {{{
    310247                if length(sol)==md.mesh.numberofelements
    311                         % map on to the vertices
    312248                        for jj=1:md.mesh.numberofelements
    313249                                ii=(jj-1)*3;
    … …  
    319255                        end
    320256                        sol=temp2;
    321                 end % }}}
     257                end
    322258
    323259                vidObj = VideoWriter(movie_name{kk});
    … …  
    326262
    327263                for jj=1:length(time)
    328                         % Make a interpolation object
    329264                        F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj));
    330265                        F.Method = 'linear';
    331266                        F.ExtrapolationMethod = 'linear';
    332267                       
    333                         % Do the interpolation to get gridded solutions...
    334268                        sol_grid = F(lat_grid, lon_grid);
    335269                        sol_grid(isnan(sol_grid))=0;
    336                         sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan
     270                        sol_grid(lat_grid>85 & sol_grid==0) = NaN;
    337271
    338272                        set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
    … …  
    340274                        gcf; load coast; cla;
    341275                        pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
    342                         % mask out land or oceans {{{
    343276                        if (kk==1)
    344277                                geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
    345278                        else
    346279                                geoshow(lat,long,'DisplayType','polygon','FaceColor','white');
    347                         end % }}}
     280                        end
    348281                        plot(long,lat,'k'); hold off;
    349                         % define colormap, caxis, xlim etc {{{
    350282                        c1=colorbar;
    351283                        colormap('haxby');
    … …  
    353285                        xlim([-180 180]);
    354286                        ylim([-90 90]);
    355                         % }}}
    356287                        grid on;
    357288                        title(sol_name(kk));
    358289                        set(gcf,'color','w');
    359290                        writeVideo(vidObj,getframe(gcf));
    360                         close % close current figure
     291                        close
    361292                end
    362 
    363                 % Close the file.
    364293                close(vidObj);
    365294        end
    … …  
    371300        ylabel('GMSL [mm/yr]');
    372301        set(gcf,'color','w');
    373 
    374302        %export_fig('Fig7.pdf');
    375303
    376 end % }}}
    377 
     304end
     305
Note: See TracChangeset for help on using the changeset viewer.