Index: /issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m (revision 25938) +++ /issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m (revision 25939) @@ -2,8 +2,5 @@ addpath('../Data','../Functions'); -steps=[1:8]; - -try - % [1:8] +steps=[8]; if any(steps==1) %{{{ @@ -16,5 +13,8 @@ md.mesh=gmshplanet('radius',radius,'resolution',resolution); - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + ocean_mask=gmtmask(md.mesh.lat,md.mesh.long); + pos = find(ocean_mask==0); + md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); + md.mask.ocean_levelset(pos)=1; save ./Models/SlrGRACE_Mesh md; @@ -39,14 +39,14 @@ for i=1:numrefine, - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + ocean_mask=gmtmask(md.mesh.lat,md.mesh.long); distance=zeros(md.mesh.numberofvertices,1); - pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); - pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); + 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); for j=1:md.mesh.numberofvertices phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; - if md.mask.ocean_levelset(j), + if ocean_mask(j), phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); @@ -61,5 +61,5 @@ pos=find(distancemindistance_land); + pos2=find(ocean_mask~=1 & distance>mindistance_land); distance(pos2)=mindistance_land; @@ -72,5 +72,8 @@ end - md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + ocean_mask=gmtmask(md.mesh.lat,md.mesh.long); + pos = find(ocean_mask==0); + md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); + md.mask.ocean_levelset(pos)=1; save ./Models/SlrGRACE_Mesh_refined md; @@ -127,9 +130,9 @@ mask = domain_mask(lat_element_0,lon_element_0,domain); - md.slr.deltathickness = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice; + md.solidearth.surfaceload.icethicknesschange = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice; save ./Models/SlrGRACE_Loads md; - plotmodel(md,'data',md.slr.deltathickness,'view',[90 90],'title','Ice height equivalent [m]'); + plotmodel(md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 90],'title','Ice height equivalent [m]'); end %}}} @@ -140,25 +143,18 @@ md.solidearth.lovenumbers = lovenumbers('maxdeg',10000); - md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); - md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); - pos=find(md.slr.deltathickness~=0); - md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; - md.mask.land_levelset = 1-md.mask.ocean_levelset; - - md.slr.sealevel=zeros(md.mesh.numberofvertices,1); - md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); - md.slr.ocean_area_scaling = 1; - - md.slr.spcthickness=NaN*zeros(md.mesh.numberofvertices,1); - md.slr.Ngia=zeros(md.mesh.numberofvertices,1); - md.slr.Ugia=zeros(md.mesh.numberofvertices,1); - - md.slr.geodetic=1; - - di=md.materials.rho_ice/md.materials.rho_water; - md.geometry.thickness=ones(md.mesh.numberofvertices,1); - md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); - md.geometry.base=md.geometry.surface-md.geometry.thickness; - md.geometry.bed=md.geometry.base; + md.mask.ice_levelset=-md.mask.ocean_levelset; + + md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); + md.dsl.global_average_thermosteric_sea_level_change=[0;0]; + md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); + md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); + + md.solidearth.settings.ocean_area_scaling=1; + + % arbitary to pass consistency check. + md.geometry.bed=-ones(md.mesh.numberofvertices,1); + md.geometry.surface=ones(md.mesh.numberofvertices,1); + md.geometry.base=md.geometry.bed; + md.geometry.thickness=md.geometry.surface-md.geometry.base; md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); @@ -175,7 +171,9 @@ md = loadmodel('./Models/SlrGRACE_Parameterization'); - md.slr.rigid=1; - md.slr.elastic=1; - md.slr.rotation=1; + md.solidearth.settings.rigid=1; + md.solidearth.settings.elastic=1; + md.solidearth.settings.rotation=1; + + md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'}; md.cluster=generic('name',oshostname(),'np',3); @@ -191,6 +189,6 @@ md = loadmodel('./Models/SlrGRACE_Solution'); - sol1 = md.slr.deltathickness*100; % [cm] - sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] + sol1 = md.solidearth.surfaceload.icethicknesschange*100; % [cm] + sol2 = md.results.SealevelriseSolution.SealevelRSL*1000; % [mm] sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; @@ -252,12 +250,12 @@ disp('Projecting loads onto the mesh...'); - time_range = [2005+15/365 2014+350/365]; - %time_range = 2007 + [15 350]/365; + %time_range = [2005+15/365 2014+350/365]; + time_range = 2007 + [15 350]/365; water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); [ne,nt]=size(water_load); - md.slr.deltathickness = zeros(ne+1,nt); - md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; - md.slr.deltathickness(ne+1,:)=[1:nt]; % times + md.solidearth.surfaceload.icethicknesschange = zeros(ne+1,nt); + md.solidearth.surfaceload.icethicknesschange(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + md.solidearth.surfaceload.icethicknesschange(ne+1,:)=[1:nt]; % times md.transient.issmb=0; @@ -266,4 +264,5 @@ md.transient.isthermal=0; md.transient.isslr=1; + md.transient.isgia=1; md.timestepping.start_time=-1; @@ -273,4 +272,6 @@ md.settings.output_frequency=1; + md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'}; + md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('111111111'); @@ -285,5 +286,5 @@ md = loadmodel('./Models/SlrGRACE_Transient'); - time = md.slr.deltathickness(end,:); + time = md.solidearth.surfaceload.icethicknesschange(end,:); tide_x = 37+29/60; % degree N @@ -291,6 +292,6 @@ for tt=1:length(time) - gmsl(tt) = md.results.TransientSolution(tt).SealevelRSLEustatic*1000; - sol = (md.results.TransientSolution(tt+1).Sealevel-md.results.TransientSolution(tt).Sealevel)*1000; + gmsl(tt) = md.results.TransientSolution(tt).Bslr*1000; + sol = md.results.TransientSolution(tt+1).SealevelRSL*1000; slr_incheon(tt) = griddata(md.mesh.lat,md.mesh.long,sol,tide_x,tide_y); end