Changeset 25939


Ignore:
Timestamp:
01/15/21 13:23:34 (4 years ago)
Author:
adhikari
Message:

CHG: upgraded the script to comply with new solid earth formalism

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r25915 r25939  
    22addpath('../Data','../Functions');
    33
    4 steps=[1:8];
    5 
    6 try
    7  % [1:8]
     4steps=[8];
    85
    96if any(steps==1) %{{{
    … …  
    1613        md.mesh=gmshplanet('radius',radius,'resolution',resolution);
    1714
    18         md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
     15        ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
     16        pos = find(ocean_mask==0);
     17        md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
     18        md.mask.ocean_levelset(pos)=1;
    1919
    2020        save ./Models/SlrGRACE_Mesh md;
    … …  
    3939        for i=1:numrefine,
    4040
    41                 md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
     41                ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
    4242
    4343                distance=zeros(md.mesh.numberofvertices,1);
    4444
    45                 pos=find(~md.mask.ocean_levelset);      coaste.lat=md.mesh.lat(pos);    coaste.long=md.mesh.long(pos);
    46                 pos=find(md.mask.ocean_levelset);       coasto.lat=md.mesh.lat(pos);    coasto.long=md.mesh.long(pos);
     45                pos=find(~ocean_mask);  coaste.lat=md.mesh.lat(pos);    coaste.long=md.mesh.long(pos);
     46                pos=find(ocean_mask);   coasto.lat=md.mesh.lat(pos);    coasto.long=md.mesh.long(pos);
    4747
    4848                for j=1:md.mesh.numberofvertices
    4949                        phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
    50                         if md.mask.ocean_levelset(j),
     50                        if ocean_mask(j),
    5151                                phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi;
    5252                                deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
    … …  
    6161                pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
    6262               
    63                 pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
     63                pos2=find(ocean_mask~=1 & distance>mindistance_land);
    6464                distance(pos2)=mindistance_land;
    6565
    … …  
    7272        end
    7373
    74         md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
     74        ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
     75        pos = find(ocean_mask==0);
     76        md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
     77        md.mask.ocean_levelset(pos)=1;
    7578
    7679        save ./Models/SlrGRACE_Mesh_refined md;
    … …  
    127130        mask = domain_mask(lat_element_0,lon_element_0,domain);
    128131
    129         md.slr.deltathickness = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice;
     132        md.solidearth.surfaceload.icethicknesschange = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice;
    130133
    131134        save ./Models/SlrGRACE_Loads md;
    132135
    133         plotmodel(md,'data',md.slr.deltathickness,'view',[90 90],'title','Ice height equivalent [m]');
     136        plotmodel(md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 90],'title','Ice height equivalent [m]');
    134137end %}}}
    135138
    … …  
    140143        md.solidearth.lovenumbers = lovenumbers('maxdeg',10000);
    141144
    142         md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1);
    143         md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
    144         pos=find(md.slr.deltathickness~=0);
    145         md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
    146         md.mask.land_levelset = 1-md.mask.ocean_levelset;
    147 
    148         md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
    149         md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
    150         md.slr.ocean_area_scaling = 1;
    151 
    152         md.slr.spcthickness=NaN*zeros(md.mesh.numberofvertices,1);
    153         md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
    154         md.slr.Ugia=zeros(md.mesh.numberofvertices,1);
    155 
    156         md.slr.geodetic=1;
    157 
    158         di=md.materials.rho_ice/md.materials.rho_water;
    159         md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    160         md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
    161         md.geometry.base=md.geometry.surface-md.geometry.thickness;
    162         md.geometry.bed=md.geometry.base;
     145        md.mask.ice_levelset=-md.mask.ocean_levelset;
     146
     147        md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
     148        md.dsl.global_average_thermosteric_sea_level_change=[0;0];
     149        md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
     150        md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
     151
     152        md.solidearth.settings.ocean_area_scaling=1;
     153
     154        % arbitary to pass consistency check.
     155        md.geometry.bed=-ones(md.mesh.numberofvertices,1);
     156        md.geometry.surface=ones(md.mesh.numberofvertices,1);
     157        md.geometry.base=md.geometry.bed;
     158        md.geometry.thickness=md.geometry.surface-md.geometry.base;
    163159
    164160        md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    … …  
    175171        md = loadmodel('./Models/SlrGRACE_Parameterization');
    176172
    177         md.slr.rigid=1;
    178         md.slr.elastic=1;
    179         md.slr.rotation=1;
     173        md.solidearth.settings.rigid=1;
     174        md.solidearth.settings.elastic=1;
     175        md.solidearth.settings.rotation=1;
     176
     177        md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
    180178
    181179        md.cluster=generic('name',oshostname(),'np',3);
    … …  
    191189        md = loadmodel('./Models/SlrGRACE_Solution');
    192190
    193         sol1 = md.slr.deltathickness*100;                                               % [cm]
    194         sol2 = md.results.SealevelriseSolution.Sealevel*1000;   % [mm]
     191        sol1 = md.solidearth.surfaceload.icethicknesschange*100;                                                % [cm]
     192        sol2 = md.results.SealevelriseSolution.SealevelRSL*1000;        % [mm]
    195193
    196194        sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
    … …  
    252250
    253251        disp('Projecting  loads onto the mesh...');
    254         time_range = [2005+15/365 2014+350/365];
    255         %time_range = 2007 + [15 350]/365;
     252        %time_range = [2005+15/365 2014+350/365];
     253        time_range = 2007 + [15 350]/365;
    256254        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
    257255
    258256        [ne,nt]=size(water_load);
    259         md.slr.deltathickness = zeros(ne+1,nt);
    260         md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
    261         md.slr.deltathickness(ne+1,:)=[1:nt]; % times
     257        md.solidearth.surfaceload.icethicknesschange = zeros(ne+1,nt);
     258        md.solidearth.surfaceload.icethicknesschange(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
     259        md.solidearth.surfaceload.icethicknesschange(ne+1,:)=[1:nt]; % times
    262260
    263261        md.transient.issmb=0;
    … …  
    266264        md.transient.isthermal=0;
    267265        md.transient.isslr=1;
     266        md.transient.isgia=1;
    268267
    269268        md.timestepping.start_time=-1;
    … …  
    273272        md.settings.output_frequency=1;
    274273
     274        md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
     275       
    275276        md.cluster=generic('name',oshostname(),'np',3);
    276277        md.verbose=verbose('111111111');
    … …  
    285286        md = loadmodel('./Models/SlrGRACE_Transient');
    286287
    287         time = md.slr.deltathickness(end,:);
     288        time = md.solidearth.surfaceload.icethicknesschange(end,:);
    288289
    289290        tide_x = 37+29/60;      % degree N
    … …  
    291292
    292293        for tt=1:length(time)
    293                 gmsl(tt) = md.results.TransientSolution(tt).SealevelRSLEustatic*1000;
    294                 sol = (md.results.TransientSolution(tt+1).Sealevel-md.results.TransientSolution(tt).Sealevel)*1000;
     294                gmsl(tt) = md.results.TransientSolution(tt).Bslr*1000;
     295                sol = md.results.TransientSolution(tt+1).SealevelRSL*1000;
    295296                slr_incheon(tt) = griddata(md.mesh.lat,md.mesh.long,sol,tide_x,tide_y);
    296297        end
Note: See TracChangeset for help on using the changeset viewer.