[issm-svn] r23407 - in issm/trunk/examples: . Functions SlrGRACE_NIMS
morlighe at issm.ess.uci.edu
morlighe at issm.ess.uci.edu
Wed Oct 10 04:57:44 PDT 2018
Author: morlighe
Date: 2018-10-10 04:57:44 -0700 (Wed, 10 Oct 2018)
New Revision: 23407
Added:
issm/trunk/examples/Functions/domain_mask.m
issm/trunk/examples/SlrGRACE_NIMS/
issm/trunk/examples/SlrGRACE_NIMS/Models/
issm/trunk/examples/SlrGRACE_NIMS/runme.m
Modified:
issm/trunk/examples/Functions/grace.m
Log:
CHG: sync'ed files for workshop
Added: issm/trunk/examples/Functions/domain_mask.m
===================================================================
--- issm/trunk/examples/Functions/domain_mask.m (rev 0)
+++ issm/trunk/examples/Functions/domain_mask.m 2018-10-10 11:57:44 UTC (rev 23407)
@@ -0,0 +1,65 @@
+% mask cryospheric domains: GrIS, AIS, RGI regions
+function mask = domain_mask(lat,lon,varargin);
+
+ nv = length(lat);
+ mask=zeros(nv,1);
+
+ if (strcmp(varargin{:},'Antarctica')==1)
+ % antarctica
+ rad=30*pi/180;
+ phi0=-90*pi/180; lambda0=0*pi/180;
+ for j=1:nv
+ phi=lat(j)/180*pi; lambda=lon(j)/180*pi;
+ delPhi=abs(phi-phi0); delLambda=abs(lambda-lambda0);
+ dist0(j)=2*asin(sqrt(sin(delPhi/2).^2+cos(phi).*cos(phi0).*sin(delLambda/2).^2));
+ end
+ pos=find(dist0<rad);
+ mask(pos)=1;
+ ocean_levelset=gmtmask(lat,lon);
+ mask(ocean_levelset==1)=0;
+
+ elseif (strcmp(varargin{:},'Greenland')==1)
+
+ load('Mask_GrIS_Tundra_Glaciers.mat');
+ mask_0=griddata(lat_element,lon_element,mask_tundra1_gris2_glaciers3,lat,lon);
+ mask(mask_0>1.5 & mask_0<2.5)=1; % that would be GrIS
+ mask(isnan(mask))=0; % there are not many glaciers there
+
+ elseif (strcmp(varargin{:},'Alaska')==1)
+
+ gla_data=dlmread('rgi_glacier_fraction_0125X0125.txt');
+ gla_lat=90-gla_data(:,1);
+ gla_lon=gla_data(:,2);
+ gla_lon(gla_lon>180)=gla_lon(gla_lon>180)-360;
+ mask_reg=gla_data(:,3); % Alaska
+ mask_reg(mask_reg>0)=1;
+ mask=griddata(gla_lat,gla_lon,mask_reg,lat,lon);
+ mask(mask>0)=1;
+ mask(isnan(mask))=0; % there are not many glaciers there
+
+ elseif (strcmp(varargin{:},'HMA')==1)
+
+ gla_data=dlmread('rgi_glacier_fraction_0125X0125.txt');
+ gla_lat=90-gla_data(:,1);
+ gla_lon=gla_data(:,2);
+ gla_lon(gla_lon>180)=gla_lon(gla_lon>180)-360;
+ mask_reg=sum(gla_data(:,15:17),2); % HMA regions 13-15
+ mask_reg(mask_reg>0)=1;
+ mask=griddata(gla_lat,gla_lon,mask_reg,lat,lon);
+ mask(mask>0)=1;
+ mask(isnan(mask))=0; % there are not many glaciers there
+
+ elseif (strcmp(varargin{:},'Glaciers')==1)
+
+ gla_data=dlmread('rgi_glacier_fraction_0125X0125.txt');
+ gla_lat=90-gla_data(:,1);
+ gla_lon=gla_data(:,2);
+ gla_lon(gla_lon>180)=gla_lon(gla_lon>180)-360;
+ mask_reg=sum(gla_data(:,3:20),2); % all RGI regions excluding AIS peripheral glaciers (#19)
+ mask_reg(mask_reg>0)=1;
+ mask=griddata(gla_lat,gla_lon,mask_reg,lat,lon);
+ mask(mask>0)=1;
+ mask(isnan(mask))=0; % there are not many glaciers there
+
+ end
+
Property changes on: issm/trunk/examples/Functions/domain_mask.m
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
Modified: issm/trunk/examples/Functions/grace.m
===================================================================
--- issm/trunk/examples/Functions/grace.m 2018-10-10 11:53:15 UTC (rev 23406)
+++ issm/trunk/examples/Functions/grace.m 2018-10-10 11:57:44 UTC (rev 23407)
@@ -65,7 +65,7 @@
grace_monthly(grace_monthly<-1344.6016 | grace_monthly>858.5046)=0;
% detrend over the entire time period
- grace_monthly = detrend(grace_monthly); % 159X64800 => remove trends from each column!
+ %grace_monthly = detrend(grace_monthly); % 159X64800 => remove trends from each column!
% fill out the blanks {{{
Added: issm/trunk/examples/SlrGRACE_NIMS/runme.m
===================================================================
--- issm/trunk/examples/SlrGRACE_NIMS/runme.m (rev 0)
+++ issm/trunk/examples/SlrGRACE_NIMS/runme.m 2018-10-10 11:57:44 UTC (rev 23407)
@@ -0,0 +1,316 @@
+
+clear all;
+addpath('../Data','../Functions');
+
+steps=[1]; % [1:8]
+
+if any(steps==1)
+ disp(' Step 1: Global mesh creation');
+
+ resolution=300; % [km]
+ radius = 6.371012*10^3; % [km]
+
+ md=model;
+ md.mask=maskpsl();
+ md.mesh=gmshplanet('radius',radius,'resolution',resolution);
+
+ md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
+
+ save ./Models/SlrGRACE_Mesh md;
+
+ plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
+
+end
+
+if any(steps==2)
+ disp(' Step 2: Global mesh creation: refined coastlines');
+
+ numrefine=1;
+ resolution=150*1e3; % inital resolution [m]
+ radius = 6.371012*10^6; % mean radius of Earth, m
+ mindistance_coast=150*1e3; % coastal resolution [m]
+ mindistance_source=75*1e3; % source resolution [m]
+ mindistance_land=300*1e3; % resolution on the continents [m]
+ maxdistance=600*1e3; % max element size (on mid-oceans) [m]
+
+ md=model;
+ md.mask=maskpsl();
+ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km]
+
+ for i=1:numrefine,
+
+ md.mask.ocean_levelset=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);
+
+ 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),
+ phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi;
+ deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
+ d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
+ else
+ phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi;
+ deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
+ d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
+ end
+ distance(j)=min(d);
+ end
+ pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
+
+ pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
+ distance(pos2)=mindistance_land;
+
+ domain = 'Greenland'; %'Antarctica'; %'HMA'; %'Alaska'; %'Glaciers'
+ mask = domain_mask(md.mesh.lat,md.mesh.long,domain);
+ distance(mask>0)=mindistance_source;
+
+ dist=min(maxdistance,distance);
+ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist);
+
+ end
+
+ md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
+
+ save ./Models/SlrGRACE_Mesh_refined md;
+
+ plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
+
+end
+
+if any(steps==3)
+ disp(' Step 3: Define loads in meters of ice height equivalent');
+ md = loadmodel('./Models/SlrGRACE_Mesh_refined');
+
+ year_month = 2007+15/365;
+ time_range = [year_month year_month];
+ water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
+
+ index = md.mesh.elements;
+ lat=90-md.mesh.lat;
+ lon=md.mesh.long;
+ lon(lon<0)=180+(180+lon(lon<0));
+
+ ax_0=lat(index(:,1)); ay_0=lon(index(:,1));
+ bx_0=lat(index(:,2)); by_0=lon(index(:,2));
+ cx_0=lat(index(:,3)); cy_0=lon(index(:,3));
+
+ for ii=1:md.mesh.numberofelements
+ if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180)
+ if ay_0(ii)==0
+ ay_0(ii)=360;
+ end
+ if by_0(ii)==0
+ by_0(ii)=360;
+ end
+ if cy_0(ii)==0
+ cy_0(ii)=360;
+ end
+ end
+ end
+
+ ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2;
+ by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2;
+ cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2;
+
+ ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2;
+ by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2;
+ cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2;
+
+ lat_element=(ax_0+bx_0+cx_0)/3;
+ lon_element=(ay_0+by_0+cy_0)/3;
+
+ lat_element_0 = 90-lat_element; lon_element_0 = lon_element;
+ lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180;
+
+ domain = 'Greenland'; %'HMA'; %'Alaska'; %'Antarctica'; %'Glaciers';
+ mask = domain_mask(lat_element_0,lon_element_0,domain);
+
+ md.slr.deltathickness = 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]');
+
+end
+
+if any(steps==4)
+ disp(' Step 4: Parameterization');
+ md = loadmodel('./Models/SlrGRACE_Loads');
+
+ nlove=10001;
+ md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[];
+ md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[];
+ md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[];
+
+ md.mask.groundedice_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.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
+ md.materials.rheology_B=paterson(md.initialization.temperature);
+ md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+
+ md.miscellaneous.name='SlrGRACE';
+
+ save ./Models/SlrGRACE_Parameterization md;
+
+end
+
+if any(steps==5)
+ disp(' Step 5: Solve Slr solver');
+ md = loadmodel('./Models/SlrGRACE_Parameterization');
+
+ md.slr.rigid=1;
+ md.slr.elastic=1;
+ md.slr.rotation=1;
+
+ md.cluster=generic('name',oshostname(),'np',3);
+ md.verbose=verbose('111111111');
+
+ md=solve(md,'Slr');
+
+ save ./Models/SlrGRACE_Solution md;
+
+end
+
+if any(steps==6)
+ disp(' Step 6: Plot solutions');
+ md = loadmodel('./Models/SlrGRACE_Solution');
+
+ sol1 = md.slr.deltathickness*100; % [cm]
+ sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm]
+
+ sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
+ fig_name={'Fig_dH.pdf','Fig_slr.pdf'};
+
+ res = 1.0;
+
+ [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
+ sol_grid = zeros(size(lat_grid));
+
+ for kk=1:2
+ sol=eval(sprintf('sol%d',kk));
+
+ if length(sol)==md.mesh.numberofelements
+ for jj=1:md.mesh.numberofelements
+ ii=(jj-1)*3;
+ pp(ii+1:ii+3)=md.mesh.elements(jj,:);
+ end
+ for jj=1:md.mesh.numberofvertices
+ pos=ceil(find(pp==jj)/3);
+ temp(jj)=mean(sol(pos));
+ end
+ sol=temp';
+ end
+
+ F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol);
+ F.Method = 'linear';
+ F.ExtrapolationMethod = 'linear';
+
+ sol_grid = F(lat_grid, lon_grid);
+ sol_grid(isnan(sol_grid))=0;
+ sol_grid(lat_grid>85 & sol_grid==0) = NaN;
+
+ set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
+ figure1=figure('Position', [100, 100, 1000, 500]);
+ gcf; load coast; cla;
+ pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
+ if (kk==1)
+ geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
+ else
+ geoshow(lat,long,'DisplayType','polygon','FaceColor',[0.75 0.75 0.75]);
+ end
+ plot(long,lat,'k'); hold off;
+ c1=colorbar;
+ colormap('haxby');
+ caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]);
+ xlim([-180 180]);
+ ylim([-90 90]);
+ grid on;
+ title(sol_name(kk));
+ set(gcf,'color','w');
+ %export_fig(fig_name{kk});
+ end
+end
+
+if any(steps==7)
+ disp(' Step 7: Transient run');
+ md = loadmodel('./Models/SlrGRACE_Parameterization');
+
+ disp('Projecting loads onto the mesh...');
+ 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.transient.issmb=0;
+ md.transient.ismasstransport=0;
+ md.transient.isstressbalance=0;
+ md.transient.isthermal=0;
+ md.transient.isslr=1;
+
+ md.timestepping.start_time=-1;
+ md.timestepping.final_time=nt;
+ md.timestepping.time_step=1;
+ md.timestepping.interp_forcings=0;
+ md.settings.output_frequency=1;
+
+ md.cluster=generic('name',oshostname(),'np',3);
+ md.verbose=verbose('111111111');
+
+ md=solve(md,'tr');
+
+ save ./Models/SlrGRACE_Transient md;
+
+end
+
+if any(steps==8)
+ disp(' Step 8: Plot transient');
+ md = loadmodel('./Models/SlrGRACE_Transient');
+
+ time = md.slr.deltathickness(end,:);
+
+ tide_x = 37+29/60; % degree N
+ tide_y = 126+20/60; % degree E
+
+ 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;
+ slr_incheon(tt) = griddata(md.mesh.lat,md.mesh.long,sol,tide_x,tide_y);
+ end
+
+ plot(2005+time/12,gmsl-gmsl(1),'-*','linewidth',3); grid on; hold on;
+ plot(2005+time/12,slr_incheon-slr_incheon(1),'-*r','linewidth',3); hold off;
+ xlabel('Year');
+ ylabel('GMSL [mm]');
+ set(gcf,'color','w');
+ legend('GMSL','Incheon');
+
+end
+
More information about the issm-svn
mailing list