Index: /issm/trunk-jpl/examples/EsaGRACE/runme.m =================================================================== --- /issm/trunk-jpl/examples/EsaGRACE/runme.m (revision 22808) +++ /issm/trunk-jpl/examples/EsaGRACE/runme.m (revision 22809) @@ -1,5 +1,5 @@ clear all; -steps=[0]; % [1:5]; +steps=[0:5]; % [1:5]; if any(steps==0) % Download GRACE land_mass data {{{ Index: /issm/trunk-jpl/examples/Functions/grace.m =================================================================== --- /issm/trunk-jpl/examples/Functions/grace.m (revision 22809) +++ /issm/trunk-jpl/examples/Functions/grace.m (revision 22809) @@ -0,0 +1,121 @@ +% map grace loads in meters [m] of water equivalent height +function water_load = grace(index,lat,lon,tmin,tmax); + + %compute centroids using (lat,lon) data {{{ + ne = length(index); % number of elements + % lat -> [0,180]; long -> [0,360] to compute centroids + lat=90-lat; + 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)); + % find whether long is 0 or 360! This is important to compute centroids as well as elemental area + for ii=1:ne + 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 + % correction at the north pole + 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; + % correction at the south pole + 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,lon] \in [-90:90,-180,180]; + lat_element_0 = 90-lat_element; lon_element_0 = lon_element; + lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180; + % }}} + + % Monthly GRACE data + filename=['GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc']; + time_0=ncread(filename,'time'); % days since 2002-01-01 00:00:00 UTC + long_0=ncread(filename,'lon'); % longitudes: 0.27-359.75 + lati_0=ncread(filename,'lat'); % latitudes: -89.75:89.75 + rec=ncread(filename,'lwe_thickness');% rec_ensemble_mean [cm] + + time_yr = 2002+time_0/365; % [yr] + + [nn_0,mm_0] = size(squeeze(rec(:,:,1))); + for jj=1:mm_0 % chose a latitude + for kk=1:nn_0 + ii=nn_0*(jj-1)+kk; + lat_0(ii)=lati_0(jj); + lon_0(ii)=long_0(kk); + tws_monthly(:,ii) = rec(kk,jj,:); + end + end + % + %% translate variables as grace-related variables -- so I do not need to do too much editing + grace_monthly=tws_monthly; + 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! + + % fill out the blanks {{{ + + lat_grace=lat_0; + lon_grace=lon_0; + num_org=length(lon_grace); + + qq=1; mm=1; + for jj=2:num_org-1 + if (lat_grace(jj)~=lat_grace(jj+1)) + lat_new(qq)=lat_grace(jj); + lon_new(qq)=lon_grace(jj)+(lon_grace(jj)-lon_grace(jj-1)); + load_new(:,qq)=grace_monthly(:,mm); + lat_new(qq+1)=lat_grace(jj); + lon_new(qq+1)=lon_grace(mm)-(lon_grace(jj)-lon_grace(jj-1)); + load_new(:,qq+1)=grace_monthly(:,jj); + qq=qq+2; + mm=jj+1; % to find out the value for monthly data + end + end + + num_add=length(lat_new); + num_plus=num_org+num_add; + + lat_grace_plus=zeros(num_plus,1); + lon_grace_plus=zeros(num_plus,1); + load_grace_plus=zeros(length(time_0),num_plus); + + lat_grace_plus(1:num_org)=lat_grace; + lat_grace_plus(1+num_org:num_plus)=lat_new; + lon_grace_plus(1:num_org)=lon_grace; + lon_grace_plus(1+num_org:num_plus)=lon_new; + load_grace_plus(:,1:num_org)=grace_monthly; + load_grace_plus(:,1+num_org:num_plus)=load_new; %(:,:); + % }}} + + % choose selected months ONLY + [diff1,pos1] = min(abs(tmin-time_yr)); + [diff2,pos2] = min(abs(tmax-time_yr)); + + time_yr=time_yr(pos1:pos2); + load_grace_plus=load_grace_plus(pos1:pos2,:); + num_yr=length(time_yr); + water_load_0=zeros(ne,num_yr); + + for jj=1:num_yr + water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element); + disp([num2str(jj),' of ',num2str(num_yr),' months done!']); + end + + water_load=water_load_0/100; % cm -> meters of water + water_load(isnan(water_load))=0; + Index: /issm/trunk-jpl/examples/Functions/wahr.m =================================================================== --- /issm/trunk-jpl/examples/Functions/wahr.m (revision 22809) +++ /issm/trunk-jpl/examples/Functions/wahr.m (revision 22809) @@ -0,0 +1,46 @@ +% compute semi-analytic solutions for a disc load +function [vert, horz] = wahr(disc_rad,xi,love_h,love_l); + + disc_rad = disc_rad/1000; % km + % compute P(x), dP(x)/dx, d2P(x)/dx2 + %--------------------------------------------------------------------- + % compute p_value + theta=km2deg(xi/1000)'; + ang = theta/180*pi; + alpha=cos(ang); + m=length(alpha); + n=length(love_h)-1; + p_value = p_polynomial_value(m,n,alpha); + p_prime = p_polynomial_prime(m,n,alpha); + %--------------------------------------------------------------------- + nn=[0:n]; + nn_plus_1=nn+1; + + % disc radius in degree + disc_rad = km2deg(disc_rad)/180*pi; + tau=zeros(size(love_h)); + tau(1) = 0.5*(1-cos(disc_rad)); % tau_0 + p_value_disc = p_polynomial_value(1,n+1,cos(disc_rad)); + p_prime_disc = p_polynomial_prime(1,n,cos(disc_rad)); + for jj=2:n+1 + nnn = jj-1; + tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1)); + end + + const=zeros(size(love_h)); + for jj=1:n+1 + const(jj) = 1/(2*(jj-1)+1); + end + + disc=sum(bsxfun(@times,p_value,tau'),2); + + g1 = -sum(bsxfun(@times,p_value,(tau.*love_h.*const)'),2); + g5 = -sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,(tau.*love_l.*const)')),2); + + % coeff + coeff = 1000*4*pi*(6.67408*10^-11)*(6.3781*10^6)/9.81; + + % vertical and horizontal solutions in mm + vert = g1*coeff*1000; % mm + horz = g5*coeff*1000; % mm + Index: /issm/trunk-jpl/examples/SlrGRACE/runme.m =================================================================== --- /issm/trunk-jpl/examples/SlrGRACE/runme.m (revision 22808) +++ /issm/trunk-jpl/examples/SlrGRACE/runme.m (revision 22809) @@ -1,5 +1,7 @@ clear all; -steps=[5]; % [1:6]; +addpath('../Data','../Functions'); + +steps=[7]; % [0:7]; if any(steps==0) % Download GRACE land_mass data {{{ @@ -11,4 +13,6 @@ cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); + + !mv *.nc '../Data/' % display the content of the netcdf file. @@ -32,5 +36,5 @@ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. - for i=1:numrefine, + for i=1:numrefine; % refine mesh... {{{ %figure out mask: @@ -66,5 +70,5 @@ %use distance to the coastline to refine mesh: md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); - end + end % }}} %figure out mask: @@ -156,7 +160,4 @@ md.slr.rotation=1; - % Request outputs -% md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; - % Cluster info md.cluster=generic('name',oshostname(),'np',3); @@ -206,5 +207,5 @@ F.Method = 'linear'; F.ExtrapolationMethod = 'linear'; - + % Do the interpolation to get gridded solutions... sol_grid = F(lat_grid, lon_grid); @@ -214,13 +215,20 @@ 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; + gcf; load coast; cla; + pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; + % mask out land or oceans {{{ + if (kk==1) + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + else + geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); + end % }}} plot(long,lat,'k'); hold off; + % define colormap, caxis, xlim etc {{{ c1=colorbar; - colormap(jet); + 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)); @@ -237,5 +245,5 @@ % read GRACE data during the chosen time period << steps=2 >> - disp('Projecting grace loads onto the mesh...'); + disp('Projecting loads onto the mesh...'); time_range = 2007 + [15 350]/365; water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); @@ -261,7 +269,4 @@ md.settings.output_frequency=1; - % Request outputs - md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; - % Cluster info md.cluster=generic('name',oshostname(),'np',3); @@ -271,4 +276,100 @@ md=solve(md,'tr'); -end % }}} - + save ./Models/SlrGRACE.Transient md; + +end % }}} + +if any(steps==7) % Plot transient {{{ + disp(' Step 7: Plot transient'); + md = loadmodel('./Models/SlrGRACE.Transient'); + + time = md.slr.deltathickness(end,:); % year + + % loads and solutions. + for tt=1:length(time) + gmsl(tt) = md.results.TransientSolution(tt).SealevelEustatic*1000; % GMSL rate mm/yr + sol1(:,tt) = md.slr.deltathickness(1:end-1,tt)*100; % ice equivalent height [cm/yr] + %% something weired happened here!!! first month solution is duplicated. Use offset of 1. Will fix it asap! + sol2(:,tt) = md.results.TransientSolution(tt+1).Sealevel*1000; % mm/yr + end + sol_name = {'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; + movie_name = {'movie_loads.avi','movie_fingerprints.avi'}; + + res = 1.0; % degree + + % Make a grid of lats and lons, based on the min and max of the original vectors + [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 data are on elements, map those on to the vertices {{{ + if length(sol)==md.mesh.numberofelements + % map on to the vertices + 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); + temp2(jj,:)=mean(sol(pos,:)); + end + sol=temp2; + end % }}} + + vidObj = VideoWriter(movie_name{kk}); + vidObj.FrameRate=2; % frames per second + open(vidObj); + + for jj=1:length(time) + % Make a interpolation object + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + + % Do the interpolation to get gridded solutions... + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to 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; + % mask out land or oceans {{{ + if (kk==1) + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + else + geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); + end % }}} + plot(long,lat,'k'); hold off; + % define colormap, caxis, xlim etc {{{ + c1=colorbar; + colormap('haxby'); + caxis([-floor(min(abs(min(min(sol))),abs(max(max(sol))))) floor(min(abs(min(min(sol))),abs(max(max(sol)))))]) + xlim([-180 180]); + ylim([-90 90]); + % }}} + grid on; + title(sol_name(kk)); + set(gcf,'color','w'); + writeVideo(vidObj,getframe(gcf)); + close % close current figure + end + + % Close the file. + close(vidObj); + end + !open *.avi; + + % plot GMSL time series + plot(time,gmsl,'-*','linewidth',3); grid on; + xlabel('# month'); + ylabel('GMSL [mm/yr]'); + set(gcf,'color','w'); + + %export_fig('Fig7.pdf'); + +end % }}} +