source: issm/trunk-jpl/examples/SlrFarrell/runme.m@ 22805

clear all;
steps=[2]; % 

if any(steps==1) % Global mesh creation {{{ 
        disp('   Step 1: Global mesh creation');

        numrefine=1;
        resolution=150*1e3;             % inital resolution [m]. It determines, e.g., whether we capture small islands. 
        radius = 6.371012*10^6; % mean radius of Earth, m
        mindistance_coast=150*1e3;              % coastal resolution [m] 
        mindistance_land=300*1e3; % resolution on the continents [m]
        maxdistance=600*1e3;     % max element size (on mid-oceans) [m]

        %mesh earth: 
        md=model; 
        md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset 
        md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. 

        for i=1:numrefine,

                %figure out mask: 
                md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 

                %figure out distance to the coastline, in lat,long (not x,y,z): 
                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
                        %figure out nearest coastline (using the great circle distance)
                        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;
                
                % refine on the continents
                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); 
                distance(pos2)=mindistance_land; 

                dist=min(maxdistance,distance); % max size 1000 km 
                %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

        %figure out mask: 
        md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 

        save ./Models/SlrFarrell.Mesh md;
        
        plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
        %export_fig('Fig1.pdf'); 

end % }}} 

if any(steps==2) % Define source {{{ 
        disp('   Step 2: Define source as in Farrell, 1972, Figure 1');
        md = loadmodel('./Models/SlrFarrell.Mesh');

        % initial sea-level: 1 m RSL everywhere. 
        md.slr.sealevel=md.mask.ocean_levelset; 

        save ./Models/SlrFarrell.Loads md; 
        
        plotmodel (md,'data',md.slr.sealevel,'view',[90 90],...
                'title#all','Initial sea-level [m]');
        %export_fig('Fig2.pdf'); 

end % }}} 

if any(steps==3) % Parameterization {{{ 
        disp('   Step 3: Parameterization');
        md = loadmodel('./Models/SlrFarrell.Loads');

        % Love numbers and reference frame: CF or CM (choose one!)  
        nlove=10001;    % up to 10,000 degree 
        md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[];
        md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[];
        md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[];

        % Mask: for computational efficiency only those elements that have loads are convolved! 
        md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded 
        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 
        pos=find(md.slr.deltathickness~=0);
        md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads 
        md.mask.land_levelset = 1-md.mask.ocean_levelset; 

        %% IGNORE BUT DO NOT DELETE %% {{{
        % Geometry: Important only when you want to couple with Ice Flow Model 
        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;
        % Materials: 
        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);
        % Miscellaneous: 
        md.miscellaneous.name='SlrFarrell';
        %% IGNORE BUT DO NOT DELETE %% }}}  
        
        save ./Models/SlrFarrell.Parameterization md; 

end % }}} 

if any(steps==4) % Solve {{{ 
        disp('   Step 4: Solve Slr solver');
        md = loadmodel('./Models/SlrFarrell.Parameterization');

        % Request outputs 
        md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'};
        
        % Cluster info 
        md.cluster=generic('name',oshostname(),'np',3); 
        md.verbose=verbose('111111111');

        % Solve 
        md=solve(md,'Slr');

        save ./Models/SlrFarrell.Solution md; 

end % }}} 

if any(steps==5) % Plot solutions {{{ 
        disp('   Step 5: Plot solutions');
        md = loadmodel('./Models/SlrFarrell.Solution');

        % loads and solutions. 
        sol1 = md.slr.deltathickness*100; % WEH cm 
        sol2 = md.results.SlrSolution.SlrUmotion*1000; % [mm] 
        sol3 = md.results.SlrSolution.SlrNmotion*1000; % [mm] 
        sol4 = md.results.SlrSolution.SlrEmotion*1000; % [mm] 
        sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',...
                'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; 

        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:4 
                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); 
                                temp(jj)=mean(sol(pos)); 
                        end
                        sol=temp'; 
                end % }}}

                % Make a interpolation object
                F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
                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; 
                plot(long,lat,'k'); hold off; 
                c1=colorbar;
                colormap(jet); 
                xlim([-180 180]); 
                ylim([-90 90]); 
                grid on; 
                title(sol_name(kk)); 
                set(gcf,'color','w');

                %export_fig('Fig5.pdf'); 
        end

end % }}}