source: issm/oecreview/Archive/25834-26739/ISSM-26257-26258.diff

Last change on this file was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 13.0 KB

  • ../trunk-jpl/examples/Functions/grace.m

     
    11% map grace loads in meters [m] of water equivalent height
    2 function water_load = grace(index,lat,lon,tmin,tmax);
     2function water_load = grace(index,lat,lon,tmin,tmax,onvertex);
    33
    44        %compute centroids using (lat,lon) data {{{
    55        ne = length(index); % number of elements
     6        nv = length(lat); % number of vertices
     7        % [lat,lon] \in [-90:90,-180,180];
     8        lat_vertex_0=lat; long_vertex_0=lon;
    69        % lat -> [0,180]; long -> [0,360] to compute centroids
    710        lat=90-lat;
    811        lon(lon<0)=180+(180+lon(lon<0));
    … …  
    115118        time_yr=time_yr(pos1:pos2);
    116119        load_grace_plus=load_grace_plus(pos1:pos2,:);
    117120        num_yr=length(time_yr);
    118         water_load_0=zeros(ne,num_yr);
     121        if onvertex
     122                water_load_0=zeros(nv,num_yr);
     123        else
     124                water_load_0=zeros(ne,num_yr);
     125        end
    119126
    120127        for jj=1:num_yr
    121                 water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element);
     128                if onvertex
     129                        water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_vertex_0,lon);
     130                else
     131                        water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element);
     132                end
    122133                disp([num2str(jj),' of ',num2str(num_yr),' months done!']);
    123134        end
    124135

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

     
    11clear all;
    22
    3 steps=[1:5];
     3steps=[4];
    44
    5 try
    6  % [1:5]
    7 
    85if any(steps==1)
    96        disp('   Step 1: Global mesh creation');
    107
    … …  
    1916        md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km]
    2017
    2118        for i=1:numrefine,
    22                 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
     19                ocean_mask_levelset=gmtmask(md.mesh.lat,md.mesh.long);
    2320
    2421                distance=zeros(md.mesh.numberofvertices,1);
    2522
    26                 pos=find(~md.mask.ocean_levelset);      coaste.lat=md.mesh.lat(pos);    coaste.long=md.mesh.long(pos);
    27                 pos=find(md.mask.ocean_levelset);       coasto.lat=md.mesh.lat(pos);    coasto.long=md.mesh.long(pos);
     23                pos=find(~ocean_mask_levelset); coaste.lat=md.mesh.lat(pos);    coaste.long=md.mesh.long(pos);
     24                pos=find(ocean_mask_levelset);  coasto.lat=md.mesh.lat(pos);    coasto.long=md.mesh.long(pos);
    2825
    2926                for j=1:md.mesh.numberofvertices
    3027                        phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
    31                         if md.mask.ocean_levelset(j),
     28                        if ocean_mask_levelset(j),
    3229                                phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi;
    3330                                deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
    3431                                d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
    … …  
    4138                end
    4239                pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
    4340
    44                 pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
     41                pos2=find(ocean_mask_levelset~=1 & distance>mindistance_land);
    4542                distance(pos2)=mindistance_land;
    4643
    4744                dist=min(maxdistance,distance);
    … …  
    4845                md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist);
    4946        end
    5047
    51         md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
     48        ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
     49        pos = find(ocean_mask==0);
     50        md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
     51        md.mask.ocean_levelset(pos)=1;
    5252
    5353        save ./Models/SlrFarrell_Mesh md;
    5454
    … …  
    5959        disp('   Step 2: Define source as in Farrell, 1972, Figure 1');
    6060        md = loadmodel('./Models/SlrFarrell_Mesh');
    6161
    62         md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere
    63         md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
    64         md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
     62        md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofvertices,1);
     63        md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
    6564
     65        pos = find(md.mask.ocean_levelset==-1);
     66        md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
     67        md.solidearth.initialsealevel(pos)=1; % 1 m SLR everywhere
     68        md.dsl.global_average_thermosteric_sea_level_change=[0;0];
     69        md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
     70        md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
     71
    6672        save ./Models/SlrFarrell_Loads md;
    6773
    68         plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]');
     74        plotmodel (md,'data',md.solidearth.initialsealevel,'view',[90 90],'title#all','Initial sea-level [m]');
    6975end
    7076
    7177if any(steps==3)
    … …  
    7480
    7581        md.solidearth.lovenumbers=lovenumbers('maxdeg',10000);
    7682
    77         md.mask.land_levelset = 1-md.mask.ocean_levelset;
    78         md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
    79         md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1);
    80         pos=find(md.mesh.lat <-80);
    81         md.mask.ice_levelset(pos(1))=-1; % ice yes!
    82         md.mask.ocean_levelset(pos(1))=1; % ice grounded!
     83        %md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
     84        %md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1);
     85        %pos=find(md.mesh.lat <-80);
     86        %md.mask.ice_levelset(pos(1))=-1; % ice yes!
     87        %md.mask.ocean_levelset(pos(1))=1; % ice grounded!
    8388
    84         di=md.materials.rho_ice/md.materials.rho_water;
    85         md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    86         md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
    87         md.geometry.base=md.geometry.surface-md.geometry.thickness;
    88         md.geometry.bed=md.geometry.base;
     89        % arbitary to pass consistency check.
     90        md.geometry.bed=-ones(md.mesh.numberofvertices,1);
     91        md.geometry.surface=ones(md.mesh.numberofvertices,1);
     92        md.geometry.base=md.geometry.bed;
     93        md.geometry.thickness=md.geometry.surface-md.geometry.base;
    8994
    9095        md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    9196        md.materials.rheology_B=paterson(md.initialization.temperature);
    … …  
    103108        md.cluster=generic('name',oshostname(),'np',3);
    104109        md.verbose=verbose('111111111');
    105110
    106         md.slr.reltol = 0.1/100; % percent change in solution
     111        md.solidearth.settings.reltol = 0.1/100; % percent change in solution
    107112
    108113        md=solve(md,'Slr');
    109114

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

     
    11clear all;
    22addpath('../Data','../Functions');
    33
    4 steps=[1]; % [1:5]
     4steps=[5];
    55
    66if any(steps==1)
    77        disp('   Step 1: Global mesh creation');

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

     
    11clear all;
    22addpath('../Functions');
    33
    4 steps=[1]; % [1:7]
     4steps=[2];
    55
    66if any(steps==1)
    77        disp('   Step 1: Mesh creation');

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

     
    11clear all;
    22addpath('../Data','../Functions');
    33
    4 steps=[7];
     4steps=[1:7];
    55
    6 if any(steps==1)
     6if any(steps==1) % {{{
    77        disp('   Step 1: Global mesh creation');
    88
    99        numrefine=5;
    … …  
    5050        save ./Models/SlrGRACE_Mesh md;
    5151
    5252        plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
    53 end
    54 if any(steps==2)
     53end % }}}
     54if any(steps==2) % {{{
    5555        disp('   Step 2: Define loads in meters of ice height equivalent');
    5656        md = loadmodel('./Models/SlrGRACE_Mesh');
    5757
    … …  
    5858        year_month = 2007+15/365;
    5959        time_range = [year_month year_month];
    6060
    61         water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
     61        onvertex = 1; % map data on vertex. If 0, it maps on the elemental centroid.
     62        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2),onvertex);
     63        rho_water2ice = md.materials.rho_freshwater/md.materials.rho_ice;
     64        ice_load = water_load*rho_water2ice; % ice height equivalent.
    6265
    63         md.solidearth.surfaceload.icethicknesschange = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
     66        %Geometry for the bed, arbitrary thickness of 100:
     67        md.geometry.bed=zeros(md.mesh.numberofvertices,1);
     68        md.geometry.base=md.geometry.bed;
     69        md.geometry.thickness = 100*ones(md.mesh.numberofvertices,1);
     70        md.geometry.surface=md.geometry.bed+md.geometry.thickness;
    6471
     72        md.masstransport.spcthickness = [md.geometry.thickness + ice_load; 0];
     73
     74        md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
     75
    6576        save ./Models/SlrGRACE_Loads md;
    6677
    67         plotmodel (md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]');
    68 end
    69 if any(steps==3)
     78        plotmodel (md,'data',ice_load,...
     79                'edgecolor','k','view',[45 45],'caxis',[-.1 .1],...
     80                'title','Ice height equivalent [m]');
     81end % }}}
     82if any(steps==3) % {{{
    7083        disp('   Step 3: Parameterization');
    7184        md = loadmodel('./Models/SlrGRACE_Loads');
    72 
     85       
     86        md.mask.ice_levelset=-md.mask.ocean_levelset;
     87       
    7388        md.solidearth.lovenumbers = lovenumbers('maxdeg',10000);
     89        md.solidearth.settings.reltol=NaN;
     90        md.solidearth.settings.abstol=1e-3;
     91        md.solidearth.settings.computesealevelchange=1;
     92        md.solidearth.settings.isgrd=1;
     93        md.solidearth.settings.grdmodel=1;
     94        md.solidearth.settings.maxiter=10;
     95       
     96        %time stepping:
     97        md.timestepping.start_time=0;
     98        md.timestepping.time_step=1;
     99        md.timestepping.final_time=1;
    74100
    75         md.mask.ice_levelset=-md.mask.ocean_levelset;
    76         md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
    77         md.dsl.global_average_thermosteric_sea_level_change=[0;0];
    78         md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
    79         md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
     101        %masstransport:
     102        md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     103        md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
     104        md.initialization.vx=zeros(md.mesh.numberofvertices,1);
     105        md.initialization.vy=zeros(md.mesh.numberofvertices,1);
     106        md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
     107        md.initialization.str=0;
    80108
    81         md.solidearth.settings.ocean_area_scaling=1;
    82 
    83         % arbitary to pass consistency check.   md.geometry.bed=-ones(md.mesh.numberofvertices,1);
    84         md.geometry.surface=ones(md.mesh.numberofvertices,1);
    85         md.geometry.base=md.geometry.bed;       md.geometry.thickness=md.geometry.surface-md.geometry.base;
    86         md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    87         md.materials.rheology_B=paterson(md.initialization.temperature);
    88         md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    89 
    90109        md.miscellaneous.name='SlrGRACE';
    91110
    92111        save ./Models/SlrGRACE_Parameterization md;
    93 end
    94 if any(steps==4)
     112end % }}}
     113if any(steps==4) % {{{
    95114        disp('   Step 4: Solve Slr solver');
    96115        md = loadmodel('./Models/SlrGRACE_Parameterization');
    97116
    … …  
    99118        md.solidearth.settings.elastic=1;
    100119        md.solidearth.settings.rotation=1;
    101120
    102         md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
     121        %Physics:
     122        md.transient.issmb=0;
     123        md.transient.isstressbalance=0;
     124        md.transient.isthermal=0;
     125        md.transient.ismasstransport=1;
     126        md.transient.isslc=1;
     127       
     128        md.solidearth.requested_outputs={'Sealevel','Bed'};
    103129
    104130        md.cluster=generic('name',oshostname(),'np',3);
    105131        md.verbose=verbose('111111111');
    106132
    107         md=solve(md,'Slr');
     133        md=solve(md,'Transient');
    108134
    109135        save ./Models/SlrGRACE_Solution md;
    110 end
    111 if any(steps==5)
     136end % }}}
     137if any(steps==5) % {{{
    112138        disp('   Step 5: Plot solutions');
    113139        md = loadmodel('./Models/SlrGRACE_Solution');
    114140
    115         sol1 = md.solidearth.surfaceload.icethicknesschange*100; % [cm]
    116         sol2 = md.results.SealevelriseSolution.SealevelRSL*1000;        % [mm]
     141        sol1 = (md.masstransport.spcthickness(1:end-1)-md.geometry.thickness)*100; % [cm]
     142        sol2 = (md.results.TransientSolution.Sealevel-md.results.TransientSolution.Bed)*1000;   % [mm]
    117143
    118144        sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
    119145        fig_name={'Fig_dH.pdf','Fig_slr.pdf'};
    … …  
    153179                if (kk==1)
    154180                        geoshow(flipud(coastlat),flipud(coastlon),'DisplayType','polygon','FaceColor','white');
    155181                else
    156                         geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor',[.5 1 .5]);
     182                        geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor','none');
    157183                end
    158184                plot(coastlon, coastlat,'k'); hold off;
    159185                c1=colorbar;
    … …  
    166192                set(gcf,'color','w');
    167193                %export_fig(fig_name{kk});
    168194        end
    169 end
    170 if any(steps==6)
     195end % }}}
     196if any(steps==6) % {{{
    171197        disp('   Step 6: Transient run');
    172198        md = loadmodel('./Models/SlrGRACE_Parameterization');
    173199
    … …  
    200226        md=solve(md,'tr');
    201227
    202228        save ./Models/SlrGRACE_Transient md;
    203 end
    204 if any(steps==7)
     229end % }}}
     230if any(steps==7) % {{{
    205231        disp('   Step 7: Plot transient');
    206232        md = loadmodel('./Models/SlrGRACE_Transient');
    207233
    … …  
    283309        ylabel('GMSL [mm/yr]');
    284310        set(gcf,'color','w');
    285311        %export_fig('Fig7.pdf');
    286 end
     312end % }}}
     313
Note: See TracBrowser for help on using the repository browser.