Changeset 25915

Timestamp:

12/29/20 13:27:22 (4 years ago)

Author:

jdquinn

Message:

CHG: Initial commit of changes for example build; clean up

Location:

issm/trunk-jpl

Files:

• examples/AMR/runme.m (modified) (5 diffs)
• examples/EsaGRACE/runme.m (modified) (3 diffs)
• examples/EsaWahr/runme.m (modified) (3 diffs)
• examples/Functions/grace.m (modified) (1 diff)
• examples/Greenland/runme.m (modified) (15 diffs)
• examples/ISMIP/CheatyIsmipA.par (deleted)
• examples/ISMIP/CheatyIsmipF.par (deleted)
• examples/ISMIP/CheatyRunme.m (deleted)
• examples/ISMIP/runme.m (modified) (11 diffs)
• examples/IceBridge/Greenland.par (modified) (1 diff)
• examples/IceBridge/runme.m (modified) (7 diffs)
• examples/Inversion/runme.m (modified) (7 diffs)
• examples/Jakobshavn/runme.m (modified) (3 diffs)
• examples/LcurveAnalysis/runme.m (modified) (10 diffs)
• examples/Pig/runme.m (modified) (10 diffs)
• examples/Pig2/runme.m (modified) (10 diffs)
• examples/PigSensitivity/runme.m (modified) (6 diffs)
• examples/SlrFarrell/runme.m (modified) (3 diffs)
• examples/SlrGRACE/runme.m (modified) (5 diffs)
• examples/SlrGRACE_NIMS/runme.m (modified) (5 diffs)
• examples/UncertaintyQuantification/runme.m (modified) (13 diffs)
• examples/shakti/runme.m (modified) (5 diffs)
• jenkins/examples_tests.sh (added)
• jenkins/pine_island-mac-examples (copied) (copied from issm/trunk-jpl/jenkins/pine_island-mac-full ) (5 diffs)
• test/NightlyRun/test418.m (modified) (1 diff)

issm/trunk-jpl/examples/AMR/runme.m

@@ -1 +1 @@
 % Mismip3D experiment with AMR using BAMG
-steps=[1:3]; 
+steps=[1:3];
 
 if any(steps==1)
         disp('   Step 1: Coarse mesh');
-        
+
         %Generate an unstructured coarse mesh on the MISMIP domain with typical element edge length equal to 10,000 m
         md=bamg(model,'domain','./domain.exp','hmax',10000,'splitcorners',1);
 
         save AMRCoarseMesh md
-end 
+end
+
 if any(steps==2)
         disp('   Step 2: Parameterization');
@@ -16 +17 @@
 
         md=setmask(md,'','');
-        
+
         % Run parameterization script to set up geometry, inital velocity, material properties, etc.
         md=parameterize(md,'./mismip.par');
@@ -23 +24 @@
         % Here, we are refining around the grounding line
         % We impose the element resolution at grounding equal to 1000 m (1 km)
-        % The criterion used is the element distance to the grounding line.
+        % The criterion used is the element distance to the grounding line
         % The distance used here is 10000 m (10 km), used in both side around the grouding line (upstream and downstream)
         md.amr.groundingline_resolution=1000;
@@ -35 +36 @@
         save AMRParam md
 end
-if any(steps==3);
-   disp('   Step 3: Solve!');
-         
+
+if any(steps==3)
+        disp('   Step 3: Solve!');
+
         md=loadmodel('AMRParam');
 
@@ -43 +45 @@
         md.timestepping.time_step=1;
         md.timestepping.final_time=500; % here, as example, only 500 yr.
-        md.settings.output_frequency=10;% here, save results every 10 yr 
+        md.settings.output_frequency=10;% here, save results every 10 yr
         md.stressbalance.maxiter=30;
         md.stressbalance.abstol=NaN;
         md.stressbalance.restol=1;
+        md.settings.solver_residue_threshold=1e-2; % relaxing (the first stress balance solver iteration presents values higher than the original threshold. This probably happens because the initial velocity is set to one).
         md.verbose=verbose('convergence',false,'solution',true);
-        
+
         % Specify that you want to run the model on your current (local host) computer
         % Change the number of processors according to your machine (here np=2)

issm/trunk-jpl/examples/EsaGRACE/runme.m

@@ -1 @@
-
 clear all;
 addpath('../Data','../Functions');
 
-steps=[1]; % [1:5] 
+steps=[1]; % [1:5]
 
-if any(steps==1) 
+if any(steps==1) %{{{
         disp('   Step 1: Global mesh creation');
 
-        resolution=300;                 % [km] 
-        radius = 6.371012*10^3; % [km] 
+        resolution=300;                 % [km]
+        radius = 6.371012*10^3; % [km]
 
-        md=model; 
-        md.mask=maskpsl(); 
+        md=model;
         md.mesh=gmshplanet('radius',radius,'resolution',resolution);
 
-        md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 
+        md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
 
         save ./Models/EsaGRACE_Mesh md;
-        
+
         plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
+end %}}}
 
-end 
-
-if any(steps==2) 
+if any(steps==2) %{{{
         disp('   Step 2: Define loads in meters of ice height equivalent');
         md = loadmodel('./Models/EsaGRACE_Mesh');
 
         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)); 
-        
-        md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent 
+        time_range = [year_month year_month];
 
-        save ./Models/EsaGRACE_Loads md; 
-        
+        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
+
+        md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent
+
+        save ./Models/EsaGRACE_Loads md;
+
         plotmodel (md,'data',md.esa.deltathickness,...
                 'view',[90 -90],'caxis',[-.1 .1],...
                 'title','Ice height equivalent [m]');
+end %}}}
 
-end 
-
-if any(steps==3) 
+if any(steps==3) %{{{
         disp('   Step 3: Parameterization');
         md = loadmodel('./Models/EsaGRACE_Loads');
 
-        nlove=10001;    
-        md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[];
-        md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[];
+        love_numbers = lovenumbers('maxdeg',10000,'referenceframe','CF');
+        md.esa.love_h = love_numbers.h;
+        md.esa.love_l = love_numbers.l;
 
-        md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); 
-        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 
+        md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1);
+        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
         pos=find(md.esa.deltathickness~=0);
-        md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; 
-        md.mask.land_levelset = 1-md.mask.ocean_levelset; 
+        md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+        %md.mask.land_levelset = 1-md.mask.ocean_levelset;
 
         di=md.materials.rho_ice/md.materials.rho_water;
@@ -61 +57 @@
         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='EsaGRACE';
-        
-        save ./Models/EsaGRACE_Parameterization md; 
 
-end 
+        save ./Models/EsaGRACE_Parameterization md;
+end %}}}
 
-if any(steps==4) 
+if any(steps==4) %{{{
         disp('   Step 4: Solve Esa solver');
         md = loadmodel('./Models/EsaGRACE_Parameterization');
 
         md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'};
-        
-        md.cluster=generic('name',oshostname(),'np',3); 
+
+        md.cluster=generic('name',oshostname(),'np',3);
         md.verbose=verbose('111111111');
 
         md=solve(md,'Esa');
 
-        save ./Models/EsaGRACE_Solution md; 
+        save ./Models/EsaGRACE_Solution md;
+end %}}}
 
-end 
-
-if any(steps==5) 
+if any(steps==5) %{{{
         disp('   Step 5: Plot solutions');
         md = loadmodel('./Models/EsaGRACE_Solution');
 
-        sol1 = md.esa.deltathickness*100;                                       % [cm] 
-        sol2 = md.results.EsaSolution.EsaUmotion*1000;  % [mm] 
-        sol3 = md.results.EsaSolution.EsaNmotion*1000;  % [mm] 
-        sol4 = md.results.EsaSolution.EsaEmotion*1000;  % [mm] 
+        sol1 = md.esa.deltathickness*100;                               % [cm]
+        sol2 = md.results.EsaSolution.EsaUmotion*1000;  % [mm]
+        sol3 = md.results.EsaSolution.EsaNmotion*1000;  % [mm]
+        sol4 = md.results.EsaSolution.EsaEmotion*1000;  % [mm]
 
         sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',...
-                'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; 
-        fig_name={'Fig_dH.pdf','Fig_vert.pdf','Fig_horzNS.pdf','Fig_horzEW.pdf'}; 
+                'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'};
+        fig_name={'Fig_dH.pdf','Fig_vert.pdf','Fig_horzNS.pdf','Fig_horzEW.pdf'};
 
-        res = 1.0; % [degree] 
+        res = 1.0; % [degree]
 
         [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
-        sol_grid = zeros(size(lat_grid)); 
+        sol_grid = zeros(size(lat_grid));
 
-        for kk=1:4 
+        for kk=1:4
                 sol=eval(sprintf('sol%d',kk));
-        
-                if length(sol)==md.mesh.numberofelements 
+
+                if length(sol)==md.mesh.numberofelements
                         for jj=1:md.mesh.numberofelements
                                 ii=(jj-1)*3;
@@ -114 +108 @@
                         end
                         for jj=1:md.mesh.numberofvertices
-                                pos=ceil(find(pp==jj)/3); 
-                                temp(jj)=mean(sol(pos)); 
+                                pos=ceil(find(pp==jj)/3);
+                                temp(jj)=mean(sol(pos));
                         end
-                        sol=temp'; 
-                end 
+                        sol=temp';
+                end
 
-                F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
+                F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol);
                 F.Method = 'linear';
-                F.ExtrapolationMethod = '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; 
+                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; 
+                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'); 
-                end 
-                plot(long,lat,'k'); hold off; 
+                        geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
+                end
+                plot(long,lat,'k'); hold off;
                 c1=colorbar;
                 colormap('haxby');
-                caxis([-min(abs(min(sol)),abs(max(sol))) min(abs(min(sol)),abs(max(sol)))]); 
-                xlim([-180 180]); 
-                ylim([-90 90]); 
-                grid on; 
-                title(sol_name(kk)); 
+                caxis([-min(abs(min(sol)),abs(max(sol))) 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}); 
+                %export_fig(fig_name{kk});
         end
-
-end 
-
+end %}}}

issm/trunk-jpl/examples/EsaWahr/runme.m

@@ -1 @@
-
 clear all;
 addpath('../Functions');
 
-steps=[0]; % [0:6] 
+steps=[1]; % [1:7]
 
-if any(steps==0) 
-        disp('   Step 0: Mesh creation');
+if any(steps==1) %{{{
+        disp('   Step 1: Mesh creation');
 
-        md=roundmesh(model,100000,10000);  % Domain radius and element size [m] 
+        md=roundmesh(model,100000,10000);  % Domain radius and element size [m]
 
         save ./Models/EsaWahr_Mesh md;
-        
+
         plotmodel(md,'data','mesh');
+end %}}}
 
-end 
+if any(steps==2) %{{{
+        disp('   Step 2: Anisotropic mesh creation');
 
-if any(steps==1) 
-        disp('   Step 1: Anisotropic mesh creation');
+        md=roundmesh(model,100000,1000);
 
-        md=roundmesh(model,100000,1000); 
+        disc_radius=20*1000;
+        rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2);       
+        field = abs(rad_dist-disc_radius);
 
-        disc_radius=20*1000; 
-        rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2);       
-        field = abs(rad_dist-disc_radius); 
-
-        md = bamg(md,'field',field,'err',50,'hmax',10000); 
+        md = bamg(md,'field',field,'err',50,'hmax',10000);
 
         save ./Models/EsaWahr_Mesh md;
-        
+
         plotmodel (md,'data','mesh');
+end %}}}
 
-end 
-
-if any(steps==2) 
-        disp('   Step 2: Define loads');
+if any(steps==3) %{{{
+        disp('   Step 3: Define loads');
         md = loadmodel('./Models/EsaWahr_Mesh');
 
-        rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; 
+        rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice;
 
         index=md.mesh.elements;         
-        x_cent=mean(md.mesh.x(index),2); 
-        y_cent=mean(md.mesh.y(index),2); 
+        x_cent=mean(md.mesh.x(index),2);
+        y_cent=mean(md.mesh.y(index),2);
 
-        md.esa.deltathickness = zeros(md.mesh.numberofelements,1); 
-        disc_radius=20; % [km] 
+        md.esa.deltathickness = zeros(md.mesh.numberofelements,1);
+        disc_radius=20; % [km]
         rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000;        
         md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i;
 
-        save ./Models/EsaWahr_Loads md; 
-        
+        save ./Models/EsaWahr_Loads md;
+
         plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]');
+end %}}}
 
-end 
-
-if any(steps==3) 
-        disp('   Step 3: Parameterization');
+if any(steps==4) %{{{
+        disp('   Step 4: Parameterization');
         md = loadmodel('./Models/EsaWahr_Loads');
 
-        nlove=10001;
-        md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[];
-        md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[];
+        love_numbers = lovenumbers('maxdeg',10000,'referenceframe','CF');
+        md.esa.love_h = love_numbers.h;
+        md.esa.love_l = love_numbers.l;
 
-        md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); 
-        md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1); 
+        md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1);
+        md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1);
 
         di=md.materials.rho_ice/md.materials.rho_water;
@@ -70 +66 @@
         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='EsaWahr';
-        
-        save ./Models/EsaWahr_Parameterization md; 
 
-end 
+        save ./Models/EsaWahr_Parameterization md;
+end %}}}
 
-if any(steps==4) 
-        disp('   Step 4: Solve Esa solver');
+if any(steps==5) %{{{
+        disp('   Step 5: Solve Esa solver');
         md = loadmodel('./Models/EsaWahr_Parameterization');
 
         md.esa.requested_outputs = {'EsaUmotion','EsaXmotion','EsaYmotion'};
-        
-        md.cluster=generic('name',oshostname(),'np',3); 
+
+        md.cluster=generic('name',oshostname(),'np',3);
         md.verbose=verbose('111111111');
 
         md=solve(md,'Esa');
 
-        save ./Models/EsaWahr_Solution md; 
+        save ./Models/EsaWahr_Solution md;
+end %}}}
 
-end 
-
-if any(steps==5) 
-        disp('   Step 5: Plot solutions');
+if any(steps==6) %{{{
+        disp('   Step 6: Plot solutions');
         md = loadmodel('./Models/EsaWahr_Solution');
 
-        vert = md.results.EsaSolution.EsaUmotion*1000;          % [mm] 
-        horz_n = md.results.EsaSolution.EsaYmotion*1000;        % [mm] 
-        horz_e = md.results.EsaSolution.EsaXmotion*1000;        % [mm] 
-        horz = sqrt(horz_n.^2+horz_e.^2);                                               % [mm]  
+        vert = md.results.EsaSolution.EsaUmotion*1000;          % [mm]
+        horz_n = md.results.EsaSolution.EsaYmotion*1000;        % [mm]
+        horz_e = md.results.EsaSolution.EsaXmotion*1000;        % [mm]
+        horz = sqrt(horz_n.^2+horz_e.^2);                                               % [mm] 
 
-        set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); 
+        set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6);
         figure('Position', [100, 100, 800, 600]);
         plotmodel(md,'data',vert,...
@@ -125 +119 @@
                 'caxis#3-4',[-0.5 0.5],...
                 'axispos',[0.505 0.02 0.47 0.47],...
-                'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); 
-        %export_fig('Fig5.pdf'); 
+                'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]');
+        %export_fig('Fig5.pdf');
+end %}}}
 
-end 
-
-if any(steps==6) 
-        disp('   Step 6: Compare results against Wahr semi-analytic solutions');
+if any(steps==7) %{{{
+        disp('   Step 7: Compare results against Wahr semi-analytic solutions');
         md = loadmodel('./Models/EsaWahr_Solution');
 
-        vert = md.results.EsaSolution.EsaUmotion*1000;          % [mm] 
-        horz_n = md.results.EsaSolution.EsaYmotion*1000;        % [mm] 
-        horz_e = md.results.EsaSolution.EsaXmotion*1000;        % [mm] 
-        horz = sqrt(horz_n.^2+horz_e.^2);                                               % [mm]  
-        
-        xi=[0:500:100000]; % grid points [m] 
-        yi=zeros(1,length(xi)); 
-        vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); 
-        horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); 
+        vert = md.results.EsaSolution.EsaUmotion*1000;          % [mm]
+        horz_n = md.results.EsaSolution.EsaYmotion*1000;        % [mm]
+        horz_e = md.results.EsaSolution.EsaXmotion*1000;        % [mm]
+        horz = sqrt(horz_n.^2+horz_e.^2);                                               % [mm] 
 
-        % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1) 
-        disc_radius = 20*1000; % [m] 
+        xi=[0:500:100000]; % grid points [m]
+        yi=zeros(1,length(xi));
+        vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear');
+        horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear');
+
+        % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1)
+        disc_radius = 20*1000; % [m]
         [vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l);
 
-        set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); 
+        set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6);
         figure1=figure('Position', [100, 100, 700, 400]);
-        ylabel_1={'0',' ','1','','2','','3',''}; 
+        ylabel_1={'0',' ','1','','2','','3',''};
         axes1 = axes('Layer','top','Position',[0.1 0.15 0.8 0.8],...
                 'XTick',[0:10:100],'xlim',[0 100],...
-                'ylim',[0 3.5],'Ytick',[0:0.5:3.5],'yticklabel',ylabel_1); 
-                box(axes1,'on'); hold(axes1,'all'); grid on; 
-                xlabel(axes1,'Radial distance [km]'); 
+                'ylim',[0 3.5],'Ytick',[0:0.5:3.5],'yticklabel',ylabel_1);
+                box(axes1,'on'); hold(axes1,'all'); grid on;
+                xlabel(axes1,'Radial distance [km]');
                 ylabel(axes1,'Displacement [mm]');
-                plot([20 20],[0 3.5],'-k','linewidth',2,'parent',axes1); 
-                % analytic soln 
-                h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1); 
-                h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1); 
-                % ISSM soln 
-                h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1); 
-                h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); 
+                plot([20 20],[0 3.5],'-k','linewidth',2,'parent',axes1);
+                % analytic soln
+                h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1);
+                h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1);
+                % ISSM soln
+                h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1);
+                h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1);
                 ag1 = gca;
                 leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)');
-                set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16); 
+                set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16);
                 set(gcf,'color','w');
-        %export_fig('Fig6.pdf'); 
-
-end 
-
+        %export_fig('Fig6.pdf');
+end %}}}

issm/trunk-jpl/examples/Functions/grace.m

@@ -43 +43 @@
 
         % 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]
+        filename=['../Data/GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'];
+        try
+                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]
+        catch e
+                disp('If dataset file exists at another location, modify the path in examples/Functions/grace.m')
+                rethrow(e)
+        end
+
 
         time_yr = 2002+time_0/365; % [yr] 

issm/trunk-jpl/examples/Greenland/runme.m

@@ -5 +5 @@
 ncdata='../Data/Greenland_5km_dev1.2.nc';
 
-if any(steps==1)
+if any(steps==1) %{{{
         disp('   Step 1: Mesh creation');
 
@@ -22 +22 @@
         %Mesh Greenland
         md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel,'err',8);
-        
+
         %convert x,y coordinates (Polar stereo) to lat/lon
         [md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,+1,39,71);
@@ -29 +29 @@
 
         plotmodel (md,'data','mesh');
-end
-
-if any(steps==2)
+end %}}}
+
+if any(steps==2) %{{{
         disp('   Step 2: Parameterization');
         md = loadmodel('./Models/Greenland.Mesh_generation');
@@ -39 +39 @@
         md = setflowequation(md,'SSA','all');
 
-        save ./Models/Greenland.Parameterization md; 
-end
-
-if any(steps==3)
+        save ./Models/Greenland.Parameterization md;
+end %}}}
+
+if any(steps==3) %{{{
         disp('   Step 3: Control method friction');
         md = loadmodel('./Models/Greenland.Parameterization');
@@ -66 +66 @@
 
         %Additional parameters
-        md.stressbalance.restol=0.01; md.stressbalance.reltol=0.1; 
+        md.stressbalance.restol=0.01; md.stressbalance.reltol=0.1;
         md.stressbalance.abstol=NaN;
 
@@ -78 +78 @@
         md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
 
-        save ./Models/Greenland.Control_drag md; 
-end
-
-if any(steps==4)
+        save ./Models/Greenland.Control_drag md;
+end %}}}
+
+if any(steps==4) %{{{
         disp('   Step 4: Transient run');
         md = loadmodel('./Models/Greenland.Control_drag');
@@ -101 +101 @@
         %Additional options
         md.inversion.iscontrol=0;
-        md.transient.requested_outputs={'IceVolume','TotalSmb', ...
-                                             'SmbMassBalance'};
+        md.transient.requested_outputs={'IceVolume','TotalSmb','SmbMassBalance'};
         md.verbose=verbose('solution',true,'module',true,'convergence',true);
 
@@ -109 +108 @@
         md=solve(md,'Transient');
 
-        save ./Models/Greenland.Transient md; 
-end
-
-if any(steps==5)
-        disp('   Step 5: Plotting'); 
+        save ./Models/Greenland.Transient md;
+end %}}}
+
+if any(steps==5) %{{{
+        disp('   Step 5: Plotting');
         md = loadmodel('./Models/Greenland.Transient');
 
@@ -143 +142 @@
         subplot(3,1,3); plot([0.2:0.2:20],volume); title('Ice Volume');
         xlabel('years')
-end
-
-if any(steps==6)
+end %}}}
+
+if any(steps==6) %{{{
         disp('   Step 6: Extract Box SMB');
         md = loadmodel('./Models/Greenland.Transient');
@@ -172 +171 @@
 
         clear smbbox
-
-end
-
-if any(steps==7)
+end %}}}
+
+if any(steps==7) %{{{
         disp('   Step 7: Historical Relaxation run');
         md = loadmodel('./Models/Greenland.Control_drag');
@@ -198 +196 @@
         %Additional options
         md.inversion.iscontrol=0;
-        md.transient.requested_outputs={'IceVolume','TotalSmb', ...
-                'SmbMassBalance'};
+        md.transient.requested_outputs={'IceVolume','TotalSmb','SmbMassBalance'};
         md.verbose=verbose('solution',true,'module',true);
 
@@ -207 +204 @@
 
         save ./Models/Greenland.HistoricTransient md;
-end
-
-if any(steps==8)
+end % step 7 end
+
+if any(steps==8) %{{{
         disp('   Step 8: Plotting exercise');
 
         %Load historic transient model
 
-
-        %Create Line Plots of relaxation run.  Create a figure.
-
+        %Create Line Plots of relaxation run. Create a figure.
 
         %Save surface mass balance, by looping through 200 years (1000 steps)
-        % Note, the first output will always contain output from time step 1
-
+        %Note, the first output will always contain output from time step 1
 
         %Plot surface mass balance time series in first subplot
 
-        
         %Title this plot Mean surface mass balance
 
-
         %Save velocity by looping through 200 years
 
-
         %Plot velocity time series in second subplot
 
-        
         %Title this plot Mean Velocity
 
-
-        %Save Ice Volume by looping through 200 years 
-
+        %Save Ice Volume by looping through 200 years
 
         %Plot volume time series in third subplot
 
-        
         %Title this plot Mean Velocity and add an x label of years
 
-end
-
-if any(steps==9)
+end % step 8 end
+
+if any(steps==9) %{{{
         disp('   Step 9: Box Transient run');
         md = loadmodel('./Models/Greenland.HistoricTransient_200yr');
@@ -292 +279 @@
         %Additional options
         md.inversion.iscontrol=0;
-        md.transient.requested_outputs={'IceVolume','TotalSmb', ...
-                'SmbMassBalance'};
+        md.transient.requested_outputs={'IceVolume','TotalSmb','SmbMassBalance'};
         md.verbose=verbose('solution',true,'module',true);
 
@@ -301 +287 @@
 
         save ./Models/Greenland.BoxTransient md;
-end
-
-if any(steps==10)
+end %}}}
+
+if any(steps==10) %{{{
         disp('   Step 10: Plot Box Transient');
         md = loadmodel('./Models/Greenland.BoxTransient');
@@ -327 +313 @@
         subplot(3,1,3); plot(t,volume); title('Ice Volume');
         xlabel('years')
-end
+end %}}}

issm/trunk-jpl/examples/ISMIP/runme.m

@@ -1 @@
-%wich steps to perform steps are from 1 to 8 
+%which steps to perform; steps are from 1 to 8
 %step 7 is specific to ISMIPA
 %step 8 is specific to ISMIPF
@@ -11 +11 @@
 
 %Mesh Generation #1
-if any(steps==1) 
+if any(steps==1) %{{{
 
         %initialize md as a new model #help model
@@ -26 +26 @@
         %->
 
-end
+end %}}}
 
 %Masks #2
-if any(steps==2) 
+if any(steps==2) %{{{
 
         % load the preceding step #help loadmodel
@@ -45 +45 @@
         %->
 
-end
+end %}}}
 
 %Parameterization #3
-if any(steps==3) 
+if any(steps==3) %{{{
 
         % load the preceding step #help loadmodel
@@ -62 +62 @@
         %->
 
-end
+end %}}}
 
 %Extrusion #4
-if any(steps==4)
+if any(steps==4) %{{{
         
         % load the preceding step #help loadmodel
@@ -81 +81 @@
         %->
 
-end
+end %}}}
 
 %Set the flow computing method #5
-if any(steps==5)
+if any(steps==5) %{{{
 
         % load the preceding step #help loadmodel
@@ -97 +97 @@
         %->
 
-end
+end %}}}
 
 %Set Boundary Conditions #6
-if any(steps==6)
+if any(steps==6) %{{{
 
         % load the preceding step #help loadmodel
@@ -118 +118 @@
         %->
 
-  % set the sliding to zero on the bed (Vx and Vy)
+        % set the sliding to zero on the bed (Vx and Vy)
         %->
 
@@ -146 +146 @@
                 % if we are dealing with IsmipF the solution is in masstransport
                 md.masstransport.vertex_pairing=md.stressbalance.vertex_pairing;
-  end
-        % save the given model
-        %->
-
-end
+        end
+        % save the given model
+        %->
+
+end %}}}
 
 %Solving #7
-if any(steps==7)
+if any(steps==7) %{{{
         % load the preceding step #help loadmodel
         % path is given by the organizer with the name of the given step
@@ -175 +175 @@
         % plot the surface velocities #plotdoc
         %->
-end
+end %}}}
 
 %Solving #8
-if any(steps==8)
+if any(steps==8) %{{{
         % load the preceding step #help loadmodel
         % path is given by the organizer with the name of the given step
@@ -212 +212 @@
         % plot the surface velocities #plotdoc
         %->
-
-end
+end %}}}

issm/trunk-jpl/examples/IceBridge/Greenland.par

@@ -30 +30 @@
 %Reading IceBridge data for Jakobshavn
 disp('      reading IceBridge Jakobshavn bedrock');
-fid  = fopen('../Data/Jakobshavn_2008_2011_Composite_XYZGrid.txt');
+fid  = fopen('../Data/Jakobshavn_2008_2011_Composite/grids/Jakobshavn_2008_2011_Composite_XYZGrid.txt');
 titles = fgets(fid); 
 data = fscanf(fid,'%g,%g,%g,%g,%g',[5 266400])';

issm/trunk-jpl/examples/IceBridge/runme.m

@@ -1 +1 @@
 clear all;
-close all;
 steps=[1];
 
@@ -6 +5 @@
 ncdata='../Data/Greenland_5km_dev1.2.nc';
 
-if any(steps==1)
+if any(steps==1) %{{{
         disp('   Step 1: Mesh creation');
 
@@ -24 +23 @@
         md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel,'err',8);
         return;
-        
+
         %Refine mesh in the region of Jakobshavn (resolution = 3000 m)
         hmaxVertices=NaN*ones(md.mesh.numberofvertices,1);
@@ -37 +36 @@
 
         save ./Models/Greenland.Mesh_generation md;
-end
+end %}}}
 
-if any(steps==2)
+if any(steps==2) %{{{
         disp('   Step 2: Parameterization');
         md = loadmodel('./Models/Greenland.Mesh_generation');
 
         md = setmask(md,'','');
-        md = parameterize(md,'Greenland.par');
+        md = parameterize(md,'./Greenland.par');
         md = setflowequation(md,'SSA','all');
 
         save ./Models/Greenland.Parameterization2 md; 
-end
+end %}}}
 
-if any(steps==3)
+if any(steps==3) %{{{
         disp('   Step 3: Control method friction');
         md = loadmodel('./Models/Greenland.Parameterization2');
@@ -93 +92 @@
 
         save ./Models/Greenland.Control_drag md; 
-end
+end %}}}
 
-if any(steps==4)
+if any(steps==4) %{{{
         disp('   Step 4: Transient run');
         md = loadmodel('./Models/Greenland.Control_drag');
@@ -122 +121 @@
 
         save ./Models/Greenland.Transient md; 
-end
+end %}}}
 
-if any(steps==5)
+if any(steps==5) %{{{
         disp('   Step 5: Plotting'); 
         md = loadmodel('./Models/Greenland.Transient');
@@ -155 +154 @@
         subplot(3,1,3); plot(time_plot,volume); title('Ice Volume');
         xlabel('years')
-end
+end %}}}

issm/trunk-jpl/examples/Inversion/runme.m

@@ -1 @@
-step=1;
-if any(step==1)
+steps=[1];
+
+if any(steps==1) %{{{
         %Generate observations
         md = model;
@@ -9 +10 @@
         md.cluster = generic('np',2);
         md = solve(md,'Stressbalance');
-        plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'title','"True" B',...
+        plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[1.3 1.9]*10^8,'title','"True" B',...
                 'data',md.results.StressbalanceSolution.Vel,'title','"observed velocities"')
         save model1 md
-end
-if any(step==2)
+end %}}}
+
+if any(steps==2) %{{{
         %Modify rheology, now constant
         loadmodel('model1.mat');
@@ -20 +22 @@
         %results of previous run are taken as observations
         md.inversion=m1qn3inversion();
-        md.inversion.vx_obs  = md.results.StressbalanceSolution.Vx;
-        md.inversion.vy_obs  = md.results.StressbalanceSolution.Vy;
-        md.inversion.vel_obs = md.results.StressbalanceSolution.Vel;
+        md.inversion.vx_obs             = md.results.StressbalanceSolution.Vx;
+        md.inversion.vy_obs             = md.results.StressbalanceSolution.Vy;
+        md.inversion.vel_obs    = md.results.StressbalanceSolution.Vel;
 
         md = solve(md,'Stressbalance');
-        plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'title','B first guess',...
+        plotmodel(md,'axis#all','tight','data',md.materials.rheology_B,'caxis',[1.3 1.9]*10^8,'title','B first guess',...
                 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities')
         save model2 md
-end
-if any(step==3)
+end %}}}
+
+if any(steps==3) %{{{
         %invert for ice rigidity
         loadmodel('model2.mat');
@@ -40 +43 @@
         md.inversion.cost_functions = 101;
         md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1);
-        md.inversion.min_parameters    = cuffey(273)*ones(md.mesh.numberofvertices,1);
-        md.inversion.max_parameters    = cuffey(200)*ones(md.mesh.numberofvertices,1);
+        md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1);
+        md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1);
 
         %Go solve!
@@ -48 +51 @@
         plotmodel(md,'axis#all','tight','data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'title','inferred B',...
                 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities')
-end
-if any(step==4)
+end %}}}
+
+if any(steps==4) %{{{
         %invert for ice rigidity
         loadmodel('model2.mat');
@@ -59 +63 @@
         md.inversion.maxsteps = maxsteps;
         md.inversion.cost_functions = [101 502];
-        md.inversion.cost_functions_coefficients      = ones(md.mesh.numberofvertices,1);
-        md.inversion.cost_functions_coefficients(:,2) = 10^-16*ones(md.mesh.numberofvertices,1);
-        md.inversion.min_parameters    = cuffey(273)*ones(md.mesh.numberofvertices,1);
-        md.inversion.max_parameters    = cuffey(200)*ones(md.mesh.numberofvertices,1);
+        md.inversion.cost_functions_coefficients                = ones(md.mesh.numberofvertices,1);
+        md.inversion.cost_functions_coefficients(:,2)   = 10^-16*ones(md.mesh.numberofvertices,1);
+        md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1);
+        md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1);
 
         %Go solve!
@@ -69 +73 @@
         plotmodel(md,'axis#all','tight','data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'title','inferred B',...
                 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities')
-end
+end %}}}

issm/trunk-jpl/examples/Jakobshavn/runme.m

@@ -1 +1 @@
 steps=[1];
 
-if any(steps==1)
-        disp('  Step 1: Mesh creation'); 
+if any(steps==1) %{{{
+        disp('   Step 1: Mesh creation');
         md=triangle(model,'Domain.exp',2000);
 
         %Get observed velocity field on mesh nodes
         ncdata='../Data/Greenland_5km_dev1.2.nc';
-        if ~exist(ncdata,'file'), 
+        if ~exist(ncdata,'file'),
                 error('File Greenland_5km_dev1.2.nc not downloaded in Data Directory');
         end
-        x1   = ncread(ncdata,'x1');
-        y1   = ncread(ncdata,'y1');
-        velx = ncread(ncdata,'surfvelx');
-        vely = ncread(ncdata,'surfvely');
-        vx   = InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0);
-        vy   = InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0);
-        vel  = sqrt(vx.^2+vy.^2);
+        x1              = ncread(ncdata,'x1');
+        y1              = ncread(ncdata,'y1');
+        velx    = ncread(ncdata,'surfvelx');
+        vely    = ncread(ncdata,'surfvely');
+        vx              = InterpFromGridToMesh(x1,y1,velx',md.mesh.x,md.mesh.y,0);
+        vy              = InterpFromGridToMesh(x1,y1,vely',md.mesh.x,md.mesh.y,0);
+        vel             = sqrt(vx.^2+vy.^2);
 
         %refine mesh using surface velocities as metric
         md=bamg(md,'hmin',1200,'hmax',15000,'field',vel,'err',5);
         [md.mesh.lat,md.mesh.long]  = xy2ll(md.mesh.x,md.mesh.y,+1,39,71);
-        
+
         save JksMesh md
-end 
-if any(steps==2)
-        disp('  Step 2: Parameterization');
+end %}}}
+
+if any(steps==2) %{{{
+        disp('   Step 2: Parameterization');
         md=loadmodel('JksMesh');
-        
+
         md=setmask(md,'','');
-        md=parameterize(md,'Jks.par'); 
+        md=parameterize(md,'Jks.par');
 
         save JksPar md
-end 
-if any(steps==3)
-        disp('  Step 3: Control method friction');
+end %}}}
+
+if any(steps==3) %{{{
+        disp('   Step 3: Control method friction');
         md=loadmodel('JksPar');
 
@@ -66 +68 @@
         md.cluster=generic('name',oshostname,'np',4);
         md=solve(md,'Stressbalance');
-        
+
         save JksControl md
-end 
-if any(steps==4)
-        disp('  Plotting')
+end %}}}
+
+if any(steps==4) %{{{
+        disp('   Plotting')
         md=loadmodel('JksControl');
 
@@ -82 +85 @@
                 'title','Friction Coefficient',...
                 'colorbartitle#3','(m)');
-end
+end %}}}

issm/trunk-jpl/examples/LcurveAnalysis/runme.m

@@ -1 @@
-step=[1];
-if any(step==1)
+steps=[1];
+
+if any(steps==1) %{{{
         % Generate observations
         md = model;
@@ -12 +13 @@
                 'data',md.results.StressbalanceSolution.Vel,'title','"observed velocities"')
         save model1 md
-end
-if any(step==2)
+end %}}}
+
+if any(steps==2) %{{{
         % Modify rheology, now constant
         loadmodel('model1.mat');
@@ -28 +30 @@
                 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities')
         save model2 md
-end
-if any(step==3)
+end %}}}
+
+if any(steps==3) %{{{
         % Perform L-curve analysis for ice rigidity inversion
         loadmodel('model2.mat');
@@ -42 +45 @@
         md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1);
         md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1);
-   md.verbose = verbose('solution',false,'control',true);
-        
+        md.verbose = verbose('solution',false,'control',true);
+
         % Starting L-curve analysis:
         %
@@ -56 +59 @@
         % respective value of Jo and R (R versus Jo).
         %
-        min_alpha   = 1.e-20;
-        max_alpha   = 1.e-11;
-        nstep_alpha = 30;
-        log_step    = (log10(max_alpha)-log10(min_alpha))/nstep_alpha;
-        log_alphas  = [log10(min_alpha):log_step:log10(max_alpha)];
-        alphas      = 10.^log_alphas;
-        J           = zeros(length(alphas),length(md.inversion.cost_functions)+1);
+        min_alpha       = 1.e-20;
+        max_alpha       = 1.e-11;
+        nstep_alpha     = 30;
+        log_step        = (log10(max_alpha)-log10(min_alpha))/nstep_alpha;
+        log_alphas      = [log10(min_alpha):log_step:log10(max_alpha)];
+        alphas          = 10.^log_alphas;
+        J                       = zeros(length(alphas),length(md.inversion.cost_functions)+1);
         % Loop over the alphas
         for i=1:length(alphas),
@@ -78 +81 @@
                 Jo = Jo + J(:,i); % sum of the cost functions (no regularization term). In this example, only 101
         end
-        R  = J(:,end-1)./alphas(:); % only the regularization term
-        
+        R = J(:,end-1)./alphas(:); % only the regularization term
+
         % Tip:
         % A rescale in the axes may be useful to visualize the L-curve.
@@ -86 +89 @@
         %
         % Apply a linear transformation on the original axis (Jo, R): 
-   %
+        %
         % |   1       alpha | | Jo  |   | Jo + alpha*R |   |    J    |
         % |                 | |     | = |              | = |         |
@@ -100 +103 @@
         hoffset=0.1*Jo;
         text(Jo+hoffset,R+voffset,[repmat('\alpha = ',length(alphas),1) num2str(alphas(:),'%2.0e')],...
-            'FontSize',10,'HorizontalAlignment','left','VerticalAlignment','Middle')
+                'FontSize',10,'HorizontalAlignment','left','VerticalAlignment','Middle')
         xlabel('$\mathrm{log}(\mathcal{J}_0$)','Interpreter','latex')
         ylabel('$\mathrm{log}(\mathcal{R})$','Interpreter','latex')
-end
-if any(step==4)
+end %}}}
+
+if any(steps==4) %{{{
         %invert for ice rigidity
         loadmodel('model2.mat');
@@ -114 +118 @@
         md.inversion.maxsteps = maxsteps;
         md.inversion.cost_functions = [101 502];
-        md.inversion.cost_functions_coefficients      = ones(md.mesh.numberofvertices,1);
-        md.inversion.cost_functions_coefficients(:,2) = 4.e-17*ones(md.mesh.numberofvertices,1); % here you can use the best value found for alpha
-        md.inversion.min_parameters    = cuffey(273)*ones(md.mesh.numberofvertices,1);
-        md.inversion.max_parameters    = cuffey(200)*ones(md.mesh.numberofvertices,1);
+        md.inversion.cost_functions_coefficients                = ones(md.mesh.numberofvertices,1);
+        md.inversion.cost_functions_coefficients(:,2)   = 4.e-17*ones(md.mesh.numberofvertices,1); % here you can use the best value found for alpha
+        md.inversion.min_parameters = cuffey(273)*ones(md.mesh.numberofvertices,1);
+        md.inversion.max_parameters = cuffey(200)*ones(md.mesh.numberofvertices,1);
 
         %Go solve!
@@ -124 +128 @@
         plotmodel(md,'axis#all','tight','data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'title','inferred B',...
                 'data',md.results.StressbalanceSolution.Vel,'title','modeled velocities')
-end
+end %}}}

issm/trunk-jpl/examples/Pig/runme.m

@@ -1 +1 @@
 steps=[1];
 
-if any(steps==1)   %Mesh Generation #1
+if any(steps==1) %Mesh Generation #1 %{{{
         %Mesh parameters
         domain =['./DomainOutline.exp'];
-        hinit=10000;   % element size for the initial mesh
-        hmax=40000;    % maximum element size of the final mesh
-        hmin=5000;     % minimum element size of the final mesh
-        gradation=1.7; % maximum size ratio between two neighboring elements
-        err=8;         % maximum error between interpolated and control field
+        hinit=10000;    % element size for the initial mesh
+        hmax=40000;             % maximum element size of the final mesh
+        hmin=5000;              % minimum element size of the final mesh
+        gradation=1.7;  % maximum size ratio between two neighboring elements
+        err=8;                  % maximum error between interpolated and control field
 
         % Generate an initial uniform mesh (resolution = hinit m)
@@ -14 +14 @@
 
         % Load Velocities
-        nsidc_vel='../Data/Antarctica_ice_velocity.nc';         
+        nsidc_vel='../Data/Antarctica_ice_velocity.nc';
 
         % Get necessary data to build up the velocity grid
-        xmin    = ncreadatt(nsidc_vel,'/','xmin');
-        ymax    = ncreadatt(nsidc_vel,'/','ymax');
-        spacing = ncreadatt(nsidc_vel,'/','spacing');
-        nx      = double(ncreadatt(nsidc_vel,'/','nx'));
-        ny      = double(ncreadatt(nsidc_vel,'/','ny'));
-        vx      = double(ncread(nsidc_vel,'vx'));
-        vy      = double(ncread(nsidc_vel,'vy'));
+        xmin    = ncreadatt(nsidc_vel,'/','xmin');
+        ymax    = ncreadatt(nsidc_vel,'/','ymax');
+        spacing = ncreadatt(nsidc_vel,'/','spacing');
+        nx              = double(ncreadatt(nsidc_vel,'/','nx'));
+        ny              = double(ncreadatt(nsidc_vel,'/','ny'));
+        vx              = double(ncread(nsidc_vel,'vx'));
+        vy              = double(ncread(nsidc_vel,'vy'));
 
         % Read coordinates
-        xmin = strtrim(xmin);  
-        xmin = str2num(xmin(1:end-2)); 
-        ymax = strtrim(ymax);  
-        ymax = str2num(ymax(1:end-2));  
+        xmin = strtrim(xmin);
+        xmin = str2num(xmin(1:end-2));
+        ymax = strtrim(ymax);
+        ymax = str2num(ymax(1:end-2));
         spacing = strtrim(spacing);
-        spacing = str2num(spacing(1:end-2));  
+        spacing = str2num(spacing(1:end-2));
 
         % Build the coordinates
         x=xmin+(0:1:nx)'*spacing;
         y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
-        
+
         % Interpolate velocities onto coarse mesh
         vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
@@ -45 +45 @@
         % Adapt the mesh to minimize error in velocity interpolation
         md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',vel_obs,'err',err);
-        
+
         %ploting
         plotmodel(md,'data','mesh')
@@ -51 +51 @@
         % Save model
         save ./Models/PIG_Mesh_generation md;
-end
+end %}}}
 
-if any(steps==2)  %Masks #2
+if any(steps==2) %Masks #2 %{{{
+        md = loadmodel('./Models/PIG_Mesh_generation');
 
-        md = loadmodel('./Models/PIG_Mesh_generation'); 
-
-        % Load SeaRISe dataset for Antarctica  
+        % Load SeaRISe dataset for Antarctica
         % http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica
         searise='../Data/Antarctica_5km_withshelves_v0.75.nc';
-        
+
         %read thickness mask from SeaRISE
         x1=double(ncread(searise,'x1'));
         y1=double(ncread(searise,'y1'));
         thkmask=double(ncread(searise,'thkmask'));
-        
+
         %interpolate onto our mesh vertices
         groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0));
@@ -77 +76 @@
         %ploting
         plotmodel(md,'data',md.mask.ocean_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice')
-        
+
         % Save model
         save ./Models/PIG_SetMask md;
-end
+end %}}}
 
-if any(steps==3)  %Parameterization #3
-
+if any(steps==3) %Parameterization #3 %{{{
         md = loadmodel('./Models/PIG_SetMask');
         md = parameterize(md,'./Pig.par');
@@ -89 +87 @@
         % Use a MacAyeal flow model
         md = setflowequation(md,'SSA','all');
-        
+
         % Save model
         save ./Models/PIG_Parameterization md;
-end
+end %}}}
 
-if any(steps==4)  %Control Method #4
-
+if any(steps==4) %Control Method #4 %{{{
         md = loadmodel('./Models/PIG_Parameterization');
 
@@ -135 +132 @@
         % Save model
         save ./Models/PIG_Control_drag md;
-end
+end %}}}
 
-if any(steps==5) %Plot #5
-
+if any(steps==5) %Plot #5 %{{{
         md = loadmodel('./Models/PIG_Control_drag');
 
@@ -149 +145 @@
                 'caxis#1-2',([1.5,4000]),...
                 'colorbartitle#3','(m)', 'log#1-2',10);
-end
+end %}}}
 
-if any(steps==6)  %Higher-Order #6
-
+if any(steps==6) %Higher-Order #6 %{{{
         % Load Model
 
@@ -165 +160 @@
         % Save Model
 
-end
+end % step 6 end
 
-if any(steps==7)  %Plot #7
-
+if any(steps==7) %Plot #7 %{{{
         mdHO = loadmodel('./Models/PIG_ModelHO');
         mdSSA = loadmodel('./Models/PIG_Control_drag');
@@ -176 +170 @@
 
         plotmodel(mdHO,'nlines',3,'ncols',2,'axis#all','equal',...
-                                                'data',mdHO.initialization.vel,'title','Observed velocity',...
-                                                'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',...
-                                                'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',...
-                                                'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',...
-                                                'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',...
-                                                'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',...
-                                                'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),...
-                                                'colorbar#all','on','view#all',2,...
-                                                'colorbartitle#all','(m/yr)',...
-                                                'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10);
-end
+                'data',mdHO.initialization.vel,'title','Observed velocity',...
+                'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',...
+                'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',...
+                'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',...
+                'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',...
+                'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',...
+                'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),...
+                'colorbar#all','on','view#all',2,...
+                'colorbartitle#all','(m/yr)',...
+                'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10);
+end %}}}

issm/trunk-jpl/examples/Pig2/runme.m

@@ -1 @@
-steps=1;
-
-if any(steps==1) %Mesh Generation #1
+steps=[1];
+
+if any(steps==1) %Mesh Generation #1 %{{{
 
         %Mesh parameters
         domain =['./DomainOutline.exp'];
-        hinit=5000;   % element size for the initial mesh
-        hmax=40000;    % maximum element size of the final mesh
-        hmin=4000;     % minimum element size of the final mesh
-        gradation=1.7; % maximum size ratio between two neighboring elements
-        err=8;         % maximum error between interpolated and control field
+        hinit=5000;             % element size for the initial mesh
+        hmax=40000;             % maximum element size of the final mesh
+        hmin=4000;              % minimum element size of the final mesh
+        gradation=1.7;  % maximum size ratio between two neighboring elements
+        err=8;                  % maximum error between interpolated and control field
 
         % Generate an initial uniform mesh (resolution = hinit m)
@@ -15 +15 @@
 
         % Get necessary data to build up the velocity grid
-        nsidc_vel='../Data/Antarctica_ice_velocity.nc';         
-        xmin    = strsplit(ncreadatt(nsidc_vel,'/','xmin'));      xmin    = str2num(xmin{2});
-        ymax    = strsplit(ncreadatt(nsidc_vel,'/','ymax'));      ymax    = str2num(ymax{2});
-        spacing = strsplit(ncreadatt(nsidc_vel,'/','spacing'));   spacing = str2num(spacing{2});
-        nx      = double(ncreadatt(nsidc_vel,'/','nx'));
-        ny      = double(ncreadatt(nsidc_vel,'/','ny'));
-        vx      = double(ncread(nsidc_vel,'vx'));
-        vy      = double(ncread(nsidc_vel,'vy'));
+        nsidc_vel       ='../Data/Antarctica_ice_velocity.nc';
+        xmin            = strsplit(ncreadatt(nsidc_vel,'/','xmin'));    xmin    = str2num(xmin{2});
+        ymax            = strsplit(ncreadatt(nsidc_vel,'/','ymax'));    ymax    = str2num(ymax{2});
+        spacing         = strsplit(ncreadatt(nsidc_vel,'/','spacing')); spacing = str2num(spacing{2});
+        nx                      = double(ncreadatt(nsidc_vel,'/','nx'));
+        ny                      = double(ncreadatt(nsidc_vel,'/','ny'));
+        vx                      = double(ncread(nsidc_vel,'vx'));
+        vy                      = double(ncread(nsidc_vel,'vy'));
 
         % Build the coordinates
         x=xmin+(0:1:nx)'*spacing;
         y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
-        
+
         % Interpolate velocities onto coarse mesh
         vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
@@ -40 +40 @@
         plotmodel(md,'data','mesh')
         save ./Models/PIG_Mesh_generation md;
-end
-
-if any(steps==2)  %Masks #2
-
-        md = loadmodel('./Models/PIG_Mesh_generation'); 
+end %}}}
+
+if any(steps==2) %Masks #2 %{{{
+        md = loadmodel('./Models/PIG_Mesh_generation');
 
         % Load SeaRISe dataset for Antarctica  http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica
         searise='../Data/Antarctica_5km_withshelves_v0.75.nc';
-        
+
         %read thickness mask from SeaRISE
         x1=double(ncread(searise,'x1'));
         y1=double(ncread(searise,'y1'));
         thkmask=double(ncread(searise,'thkmask'));
-        
+
         %interpolate onto our mesh vertices
         groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0));
@@ -64 +63 @@
 
         plotmodel(md,'data',md.mask.ocean_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice')
-        
+
         save ./Models/PIG_SetMask md;
-end
-
-if any(steps==3)  %Parameterization #3
-
+end %}}}
+
+if any(steps==3) %Parameterization #3 %{{{
         md = loadmodel('./Models/PIG_SetMask');
         md = setflowequation(md,'SSA','all');
         md = parameterize(md,'./Pig.par');
-        
+
         save ./Models/PIG_Parameterization md;
-end
-
-if any(steps==4)  %Rheology B inversion
-
+end %}}}
+
+if any(steps==4) %Rheology B inversion %{{{
         md = loadmodel('./Models/PIG_Parameterization');
 
@@ -120 +117 @@
         % Save model
         save ./Models/PIG_Control_B md;
-end
-
-if any(steps==5)  %drag inversion
-
+end %}}}
+
+if any(steps==5) %drag inversion %{{{
         md = loadmodel('./Models/PIG_Control_B');
 
@@ -153 +149 @@
                 'caxis#1-2',([1.5,4000]),...
                 'colorbartitle#3','(m)', 'log#1-2',10);
+
         save ./Models/PIG_Control_drag md;
-end
-
-if any(steps==6) %Transient Run #1
+end %}}}
+
+if any(steps==6) %Transient Run #1 %{{{
 
         md = loadmodel('./Models/PIG_Control_drag');    
@@ -183 +180 @@
         % Save model
         save ./Models/PIG_Transient md;
-end
+end %}}}
 
 if any(steps==7) %High Melt #2
@@ -196 +193 @@
 
         save ./Models/PIG_HighMelt md;
-end
-
-if any(steps==8) %Ice Front retreat
+end %}}}
+
+if any(steps==8) %Ice Front retreat %{{{
         md = loadmodel('./Models/PIG_Transient');       
 
@@ -209 +206 @@
 
         save ./Models/PIG_FrontRetreat md;
-end
-
-if any(steps==9) %High surface mass balance #3
+end %}}}
+
+if any(steps==9) %High surface mass balance #3 %{{{
         %Load model from PIG_Transient
         %...
@@ -229 +226 @@
         %Save model
         save ./Models/PIG_SMB md;
-end
+end %}}}

issm/trunk-jpl/examples/PigSensitivity/runme.m

@@ -1 @@
-steps=[1];
+steps=[1:4];
 
-if any(steps==1) %Transient Run #1
+try
 
-        md = loadmodel('../Pig/Models/PIG_Control_drag');       
+
+if any(steps==1) %Transient Run #1 %{{{
+
+        md = loadmodel('../Pig/Models/PIG_Control_drag');
 
         md.inversion.iscontrol=0;
@@ -11 +14 @@
         md.transient.ismovingfront=0;
         md.transient.isthermal=0;
-        
+
         pos=find(md.mask.ocean_levelset<0);
         md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
         md.basalforcings.floatingice_melting_rate=25*ones(md.mesh.numberofvertices,1);
-        
+
         md.timestepping.time_step=0.1;
         md.timestepping.final_time=10;
@@ -24 +27 @@
         % Save model
         save ./Models/PIG_Transient md;
-end
+end %}}}
 
-if any(steps==2) %High Melt #2
-        md = loadmodel('./Models/PIG_Transient');       
+if any(steps==2) %High Melt #2 %{{{
+        md = loadmodel('./Models/PIG_Transient');
 
         md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
         md.basalforcings.floatingice_melting_rate=60*ones(md.mesh.numberofvertices,1);
-        
+
         md.timestepping.time_step=0.1;
         md.timestepping.final_time=10;
@@ -39 +42 @@
 
         save ./Models/PIG_HighMelt md;
-end
+end %}}}
 
-if any(steps==3) %Ice Front retreat
-        md = loadmodel('./Models/PIG_Transient');       
+if any(steps==3) %Ice Front retreat %{{{
+        md = loadmodel('./Models/PIG_Transient');
 
         md2=extract(md,'~FrontRetreat.exp');
@@ -58 +61 @@
 
         save ./Models/PIG_FrontRetreat md2;
-end
+end %}}}
 
-if any(steps==4) %High surface mass balance #3
+if any(steps==4)
+ disp('Needs work!') %High surface mass balance #3 %{{{
         %Load model
 
         %Change external forcing basal melting rate and surface mass balance)
-        
+
         %Refine time steps and time span of the simulation
 
@@ -73 +77 @@
         %Save model
 
-end
+end %}}}

issm/trunk-jpl/examples/SlrFarrell/runme.m

@@ -1 @@
-
 clear all;
 
-steps=[1]; % [1:5] 
+steps=[1:5];
 
-if any(steps==1) 
+try
+ % [1:5]
+
+if any(steps==1) %{{{
         disp('   Step 1: Global mesh creation');
 
         numrefine=1;
-        resolution=150*1e3;                     % inital resolution [m] 
+        resolution=150*1e3;                     % inital resolution [m]
         radius = 6.371012*10^6;         % mean radius of Earth, m
-        mindistance_coast=150*1e3;      % coastal resolution [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]
+        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] 
+        md=model;
+        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); 
+                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);  
+                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; 
+                        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; 
+                                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; 
+                                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));
@@ -41 +41 @@
                 end
                 pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
-                
-                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); 
-                distance(pos2)=mindistance_land; 
 
-                dist=min(maxdistance,distance); 
+                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
+                distance(pos2)=mindistance_land;
+
+                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); 
+        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');
+end %}}}
 
-end 
-
-if any(steps==2) 
+if any(steps==2) %{{{
         disp('   Step 2: Define source as in Farrell, 1972, Figure 1');
         md = loadmodel('./Models/SlrFarrell_Mesh');
 
-        md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere 
-        
+        md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere
         md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
         md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
 
-        save ./Models/SlrFarrell_Loads md; 
-        
+        save ./Models/SlrFarrell_Loads md;
+
         plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]');
+end %}}}
 
-end 
-
-if any(steps==3) 
+if any(steps==3) %{{{
         disp('   Step 3: Parameterization');
         md = loadmodel('./Models/SlrFarrell_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.solidearth.lovenumbers=lovenumbers('maxdeg',10000);
 
-        md.mask.land_levelset = 1-md.mask.ocean_levelset; 
-        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 
-        md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1); 
+        md.mask.land_levelset = 1-md.mask.ocean_levelset;
+        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
+        md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1);
         pos=find(md.mesh.lat <-80);
-        md.mask.ice_levelset(pos(1))=-1; % ice yes!  
-        md.mask.ocean_levelset(pos(1))=1; % ice grounded!  
+        md.mask.ice_levelset(pos(1))=-1; % ice yes!
+        md.mask.ocean_levelset(pos(1))=1; % ice grounded!
 
         di=md.materials.rho_ice/md.materials.rho_water;
@@ -94 +87 @@
         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='SlrFarrell';
-        
-        save ./Models/SlrFarrell_Parameterization md; 
 
-end 
+        save ./Models/SlrFarrell_Parameterization md;
+end %}}}
 
-if any(steps==4) 
+if any(steps==4) %{{{
         disp('   Step 4: Solve Slr solver');
         md = loadmodel('./Models/SlrFarrell_Parameterization');
 
-        md.cluster=generic('name',oshostname(),'np',3); 
+        md.cluster=generic('name',oshostname(),'np',3);
         md.verbose=verbose('111111111');
 
-        md.slr.reltol = 0.1/100; % per cent change in solution 
+        md.slr.reltol = 0.1/100; % percent change in solution
 
         md=solve(md,'Slr');
 
-        save ./Models/SlrFarrell_Solution md; 
+        save ./Models/SlrFarrell_Solution md;
+end %}}}
 
-end 
-
-if any(steps==5) 
+if any(steps==5) %{{{
         disp('   Step 5: Plot solutions');
         md = loadmodel('./Models/SlrFarrell_Solution');
 
-        sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m)  
+        sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m)
 
         res = 1; % [degree]
 
         [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
-        sol_grid = zeros(size(lat_grid)); 
+        sol_grid = zeros(size(lat_grid));
 
-        F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
-        F.Method = 'linear'; 
-        F.ExtrapolationMethod = 'linear'; 
-        
+        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; 
+        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; 
+        figure1=figure('Position', [100, 100, 1000, 500]);
+        gcf; load coast; cla;
         pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on;
         [C,h]=contour(lon_grid,lat_grid,sol_grid,[96 98 100 102 104 105],'-k','linewidth',2);
-        clabel(C,h,'FontSize',18,'Color','red','LabelSpacing',500); 
-        geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]); 
-        plot(long,lat,'k'); hold off; 
+        clabel(C,h,'FontSize',18,'Color','red','LabelSpacing',500);
+        geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]);
+        plot(long,lat,'k'); hold off;
         c1=colorbar;
-        colormap(flipud(haxby)); 
-        caxis([96 105]); 
-        xlim([-170 170]); 
-        ylim([-85 85]); 
-        grid on; 
-        title('Relative sea-level [% of GMSL]'); 
+        colormap(flipud(haxby));
+        caxis([96 105]);
+        xlim([-170 170]);
+        ylim([-85 85]);
+        grid on;
+        title('Relative sea-level [% of GMSL]');
         set(gcf,'color','w');
-        %export_fig('Fig5.pdf'); 
-
-end 
-
+        %export_fig('Fig5.pdf');
+end %}}}

issm/trunk-jpl/examples/SlrGRACE/runme.m

@@ -1 @@
-
 clear all;
 addpath('../Data','../Functions');
 
-steps=[1]; % [1:7]
-
-if any(steps==1) 
+steps=[1:7];
+
+try
+ % [1:7]
+
+if any(steps==1) %{{{
         disp('   Step 1: Global mesh creation');
 
         numrefine=1;
-        resolution=150*1e3;                     % inital resolution [m] 
+        resolution=150*1e3;                     % inital resolution [m]
         radius = 6.371012*10^6;         % mean radius of Earth, m
-        mindistance_coast=150*1e3;      % coastal resolution [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]
-
-        md=model; 
-        md.mask=maskpsl(); 
-        md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] 
+        maxdistance=600*1e3;            % max element size (on mid-oceans) [m]
+
+        md=model;
+        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); 
+                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);  
+                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; 
+                        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; 
+                                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; 
+                                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));
@@ -43 +44 @@
                 pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
                 
-                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); 
-                distance(pos2)=mindistance_land; 
-
-                dist=min(maxdistance,distance); 
+                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
+                distance(pos2)=mindistance_land;
+
+                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); 
+        end
+
+        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) 
+end %}}}
+
+if any(steps==2) %{{{
         disp('   Step 2: Define loads in meters of ice height equivalent');
         md = loadmodel('./Models/SlrGRACE_Mesh');
 
         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)); 
-        
-        md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; 
-
-        save ./Models/SlrGRACE_Loads md; 
-        
+        time_range = [year_month year_month];
+
+        water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
+
+        md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
+
+        save ./Models/SlrGRACE_Loads md;
+
         plotmodel (md,'data',md.slr.deltathickness,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]');
-
-end 
-
-if any(steps==3) 
+end %}}}
+
+if any(steps==3) %{{{
         disp('   Step 3: 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.ocean_levelset = ones(md.mesh.numberofvertices,1); 
-        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 
+        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.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);
@@ -100 +95 @@
         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==4) 
+
+        save ./Models/SlrGRACE_Parameterization md;
+end %}}}
+
+if any(steps==4) %{{{
         disp('   Step 4: Solve Slr solver');
         md = loadmodel('./Models/SlrGRACE_Parameterization');
 
-        md.slr.rigid=1; 
-        md.slr.elastic=1; 
+        md.slr.rigid=1;
+        md.slr.elastic=1;
         md.slr.rotation=1;
 
-        md.cluster=generic('name',oshostname(),'np',3); 
+        md.cluster=generic('name',oshostname(),'np',3);
         md.verbose=verbose('111111111');
 
         md=solve(md,'Slr');
 
-        save ./Models/SlrGRACE_Solution md; 
-
-end 
-
-if any(steps==5) 
+        save ./Models/SlrGRACE_Solution md;
+end %}}}
+
+if any(steps==5) %{{{
         disp('   Step 5: 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; 
+        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)); 
+        sol_grid = zeros(size(lat_grid));
 
         for kk=1:2
                 sol=eval(sprintf('sol%d',kk));
-        
-                if length(sol)==md.mesh.numberofelements 
+
+                if length(sol)==md.mesh.numberofelements
                         for jj=1:md.mesh.numberofelements
                                 ii=(jj-1)*3;
@@ -152 +145 @@
                         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); 
+                                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'; 
+                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; 
+                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; 
+                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'); 
+                        geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
                 else
-                        geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); 
-                end 
-                plot(long,lat,'k'); hold off; 
+                        geoshow(lat,long,'DisplayType','polygon','FaceColor','white');
+                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)); 
+                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==6) 
+                %export_fig(fig_name{kk});
+        end
+end %}}}
+
+if any(steps==6) %{{{
         disp('   Step 6: Transient run');
         md = loadmodel('./Models/SlrGRACE_Parameterization');
 
         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)); 
-
-        [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 
-        
+        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.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.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.cluster=generic('name',oshostname(),'np',3);
         md.verbose=verbose('111111111');
 
         md=solve(md,'tr');
 
-        save ./Models/SlrGRACE_Transient md; 
-
-end 
-
-if any(steps==7) 
+        save ./Models/SlrGRACE_Transient md;
+end %}}}
+
+if any(steps==7) %{{{
         disp('   Step 7: Plot transient');
         md = loadmodel('./Models/SlrGRACE_Transient');
 
-        time = md.slr.deltathickness(end,:); 
+        time = md.slr.deltathickness(end,:);
 
         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] 
-           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_dH.avi','Movie_slr.avi'}; 
-
-        res = 1.0; 
+                sol1(:,tt) = md.slr.deltathickness(1:end-1,tt)*100;                                     % ice equivalent height [cm/yr]
+                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_dH.avi','Movie_slr.avi'};
+
+        res = 1.0;
 
         [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
-        sol_grid = zeros(size(lat_grid)); 
+        sol_grid = zeros(size(lat_grid));
 
         for kk=1:2
                 sol=eval(sprintf('sol%d',kk));
-        
-                if length(sol)==md.mesh.numberofelements 
+
+                if length(sol)==md.mesh.numberofelements
                         for jj=1:md.mesh.numberofelements
                                 ii=(jj-1)*3;
@@ -251 +242 @@
                         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 
+                                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)
-                        F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); 
+                        F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj));
                         F.Method = 'linear';
-                        F.ExtrapolationMethod = '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; 
+                        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; 
+                        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'); 
+                                geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white');
                         else
-                                geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); 
-                        end 
-                        plot(long,lat,'k'); hold off; 
+                                geoshow(lat,long,'DisplayType','polygon','FaceColor','white');
+                        end
+                        plot(long,lat,'k'); hold off;
                         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)); 
+                        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 
+                        writeVideo(vidObj,getframe(gcf));
+                        close
                 end
                 close(vidObj);
         end
-        !open *.avi; 
-        
-        % plot GMSL time series 
-        plot(time,gmsl,'-*','linewidth',3); grid on; 
-        xlabel('# month'); 
+        !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 
-
+end %}}}

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

@@ -1 @@
-
 clear all;
 addpath('../Data','../Functions');
 
-steps=[1]; % [1:8]
-
-if any(steps==1) 
+steps=[1:8];
+
+try
+ % [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(); 
+        resolution=300;                 % [km]
+        radius = 6.371012*10^3; % [km]
+
+        md=model;
         md.mesh=gmshplanet('radius',radius,'resolution',resolution);
 
-        md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 
+        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) 
+end %}}}
+
+if any(steps==2) %{{{
         disp('   Step 2: Global mesh creation: refined coastlines');
 
         numrefine=1;
-        resolution=150*1e3;                     % inital resolution [m] 
+        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_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] 
+        maxdistance=600*1e3;            % max element size (on mid-oceans) [m]
+
+        md=model;
+        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); 
+                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);  
+                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; 
+                        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; 
+                                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; 
+                                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));
@@ -62 +61 @@
                 pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
                 
-                pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); 
-                distance(pos2)=mindistance_land; 
+                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);  
+                mask = domain_mask(md.mesh.lat,md.mesh.long,domain); 
                 distance(mask>0)=mindistance_source;
 
-                dist=min(maxdistance,distance); 
+                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); 
+        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) 
+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)); 
-        
+        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 
+                        end
                         if by_0(ii)==0
-                                by_0(ii)=360; 
-                        end 
-                        if cy_0(ii)==0 
-                                cy_0(ii)=360; 
-                        end
-                end 
+                                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; 
+
+        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;
+        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) 
+
+        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.ocean_levelset = ones(md.mesh.numberofvertices,1); 
-        md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 
+        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.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);
@@ -161 +154 @@
         md.slr.Ugia=zeros(md.mesh.numberofvertices,1);
 
-        md.slr.geodetic=1; 
+        md.slr.geodetic=1;
 
         di=md.materials.rho_ice/md.materials.rho_water;
@@ -168 +161 @@
         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) 
+
+        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.rigid=1;
+        md.slr.elastic=1;
         md.slr.rotation=1;
 
-        md.cluster=generic('name',oshostname(),'np',3); 
+        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) 
+        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; 
+        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)); 
+        sol_grid = zeros(size(lat_grid));
 
         for kk=1:2
                 sol=eval(sprintf('sol%d',kk));
-        
-                if length(sol)==md.mesh.numberofelements 
+
+                if length(sol)==md.mesh.numberofelements
                         for jj=1:md.mesh.numberofelements
                                 ii=(jj-1)*3;
@@ -220 +211 @@
                         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); 
+                                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'; 
+                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; 
+                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; 
+                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'); 
+                        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; 
+                        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)); 
+                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}); 
+                %export_fig(fig_name{kk});
         end
-end 
-
-if any(steps==7) 
+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 
-        
+        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.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.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.cluster=generic('name',oshostname(),'np',3);
         md.verbose=verbose('111111111');
 
         md=solve(md,'tr');
 
-        save ./Models/SlrGRACE_Transient md; 
-
-end 
-
-if any(steps==8) 
+        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
+        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; 
+                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'); 
+        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 
-
+        legend('GMSL','Incheon');
+end %}}}

issm/trunk-jpl/examples/UncertaintyQuantification/runme.m

@@ -1 +1 @@
 %PIG Uncertainty Quantification Application
-steps=[1]; 
-
-if any(steps==1) 
-        disp('   Step 1: plot flux gates'); 
-        
+steps=[1];
+
+if any(steps==1) %{{{
+        disp('   Step 1: plot flux gates');
+
         md = loadmodel('../Pig/Models/PIG_Control_drag');
 
@@ -28 +28 @@
                 '8','9','10','11','12','13'},...
                 'textposition',textpositions);
-
-end
-if any(steps==2) 
-        disp('   Step 2: compute cross overs from CRESIS'); 
-        
+end %}}}
+
+if any(steps==2) %{{{
+        disp('   Step 2: compute cross overs from CRESIS');
+
         md = loadmodel('../Pig/Models/PIG_Control_drag');
 
         %load cross overs: CRESIS McCord Antarctica, 2009 (courtesy of John Paden)
         load('../Data/CrossOvers2009.mat');
-        
+
         %interpolate cross over errors over our mesh vertices
         DeltaHH=InterpFromMeshToMesh2d(index,x,y,dhh,md.mesh.x,md.mesh.y);
@@ -46 +46 @@
         %filter out unrealistic error ranges
         flags=ContourToNodes(md.mesh.x,md.mesh.y,'ErrorContour.exp',1);
-        pos=find(~flags); DeltaHH(pos)=0; 
+        pos=find(~flags); DeltaHH(pos)=0;
 
         %avoid large unrealistic values
         pos=find(DeltaHH>1); DeltaHH(pos)=1;
         pos=find(DeltaHH<-1); DeltaHH(pos)=-1;
-        
+
         %transform into absolute errors and setup a minimum error everywhere.
         DeltaHH=abs(DeltaHH);
-        pos=find(DeltaHH==0); pos2=find(DeltaHH~=0); 
+        pos=find(DeltaHH==0); pos2=find(DeltaHH~=0);
         DeltaHH(pos)=min(DeltaHH(pos2));
 
         save Models/PIG.CrossOvers DeltaHH
-
-end
-if any(steps==3) 
-        disp('   Step 3: sampling analysis'); 
-        
+end %}}}
+
+if any(steps==3) %{{{
+        disp('   Step 3: sampling analysis');
+
         %load model and cross over errors
         md = loadmodel('../Pig/Models/PIG_Control_drag');
@@ -68 +68 @@
 
         %partition the mesh
-        md.qmu.numberofpartitions=50;
-        md=partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,...
-        'weighting','on');
-        md.qmu.partition=md.qmu.partition-1; %switch partition to c-indexing
+        npart=50;
+        partition=partitioner(md,'package','chaco','npart',npart,'weighting','on')-1;
 
         %make DeltaHH into our 3 sigma deviation
         DeltaHH=DeltaHH/6; %2 (to transform DeltaHH into a radius) x 3 (for 3 sigma)
-        DeltaHH_on_partition=AreaAverageOntoPartition(md,DeltaHH);
-        DeltaHH_on_grids=DeltaHH_on_partition(md.qmu.partition+1); %just to check in case
-        
-        md.qmu.variables.thickness=normal_uncertain('scaled_Thickness',1,1);
-        md.qmu.variables.thickness.stddev=DeltaHH_on_partition;
+        DeltaHH_on_partition=AreaAverageOntoPartition(md,DeltaHH,partition);
+        DeltaHH_on_grids=DeltaHH_on_partition(partition+1); %just to check in case
+
+        md.qmu.variables.thickness=normal_uncertain('descriptor','scaled_Thickness',...
+                'mean',ones(npart,1),...
+                'stddev',DeltaHH_on_partition,...
+                'partition',partition);
 
         %responses
@@ -110 +110 @@
 
         %mass flux profiles
-        md.qmu.mass_flux_profiles={'MassFlux1.exp',...
-                                   'MassFlux2.exp',...
-                                   'MassFlux3.exp',...
-                                   'MassFlux4.exp',...
-                                   'MassFlux5.exp',...
-                                   'MassFlux6.exp',...
-                                   'MassFlux7.exp',...
-                                   'MassFlux8.exp',...
-                                   'MassFlux9.exp',...
-                                   'MassFlux10.exp',...
-                                   'MassFlux11.exp',...
-                                   'MassFlux12.exp',...
-                                   'MassFlux13.exp'};
+        md.qmu.mass_flux_profiles={...
+                'MassFlux1.exp',...
+                'MassFlux2.exp',...
+                'MassFlux3.exp',...
+                'MassFlux4.exp',...
+                'MassFlux5.exp',...
+                'MassFlux6.exp',...
+                'MassFlux7.exp',...
+                'MassFlux8.exp',...
+                'MassFlux9.exp',...
+                'MassFlux10.exp',...
+                'MassFlux11.exp',...
+                'MassFlux12.exp',...
+                'MassFlux13.exp'...
+        };
         md.qmu.mass_flux_profile_directory='../MassFluxes/';
 
         %%  sampling analysis
-        md.qmu.method     =dakota_method('nond_samp');
-        md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),...
-        'seed',1234,...
-        'samples',30,...
-        'sample_type','lhs'); %random or lhs
+        md.qmu.method           =dakota_method('nond_samp');
+        md.qmu.method(end)      =dmeth_params_set(...
+                md.qmu.method(end),...
+                'seed',1234,...
+                'samples',30,...
+                'sample_type','lhs'...
+        ); %random or lhs
 
         %%  a variety of parameters
@@ -143 +147 @@
         %Turn off verbosity
         md.verbose=verbose(0);
-        
+
         %Here, we choose to run with 4 processors, 3 for DAKOTA
         % while one serves as the master
@@ -157 +161 @@
         md=solve(md,'Stressbalance','overwrite','y');
 
-        save ./Models/PIG.Sampling md; 
-end
-if any(steps==4) 
-        disp('   Step 4: sensitivity analysis'); 
-        
+        save ./Models/PIG.Sampling md;
+end %}}}
+
+if any(steps==4) %{{{
+        disp('   Step 4: sensitivity analysis');
+
         %load model
         md = loadmodel('../Pig/Models/PIG_Control_drag');
@@ -195 +200 @@
 
         %mass flux profiles
-        md.qmu.mass_flux_profiles={'MassFlux1.exp',...
-                                   'MassFlux2.exp',...
-                                   'MassFlux3.exp',...
-                                   'MassFlux4.exp',...
-                                   'MassFlux5.exp',...
-                                   'MassFlux6.exp',...
-                                   'MassFlux7.exp',...
-                                   'MassFlux8.exp',...
-                                   'MassFlux9.exp',...
-                                   'MassFlux10.exp',...
-                                   'MassFlux11.exp',...
-                                   'MassFlux12.exp',...
-                                   'MassFlux13.exp'};
+        md.qmu.mass_flux_profiles={...
+                'MassFlux1.exp',...
+                'MassFlux2.exp',...
+                'MassFlux3.exp',...
+                'MassFlux4.exp',...
+                'MassFlux5.exp',...
+                'MassFlux6.exp',...
+                'MassFlux7.exp',...
+                'MassFlux8.exp',...
+                'MassFlux9.exp',...
+                'MassFlux10.exp',...
+                'MassFlux11.exp',...
+                'MassFlux12.exp',...
+                'MassFlux13.exp'...
+        };
         md.qmu.mass_flux_profile_directory='../MassFluxes/';
 
@@ -227 +234 @@
         % while one serves as the master
         md.cluster=generic('name',oshostname,'np',2);
-        
+
         %Dakota runs in parallel with a master/slave configuration.
         % At least 2 cpu's are needed to run the UQ
@@ -238 +245 @@
         md.verbose=verbose('qmu',true);
         md=solve(md,'Stressbalance','overwrite','y');
-        
-        save ./Models/PIG.Sensitivity md; 
-end
-if any(steps==5) 
+
+        save ./Models/PIG.Sensitivity md;
+end %}}}
+
+if any(steps==5) %{{{
         disp('   Step 5: plot partition');
-        
+
         %load model
         md = loadmodel('./Models/PIG.Sampling');
-        
+
         plotmodel(md,'data','mesh','partitionedges','on',...
         'linewidth',2, 'axis#all','image','unit','km','colorbar','off',...
         'title','Partition Edges on ISSM mesh','grid','on');
-
-end 
-if any(steps==6) 
-        disp('   Step 6: plot histogram'); 
-        
+end %}}}
+
+if any(steps==6) %{{{
+        disp('   Step 6: plot histogram');
+
         %load model
         md = loadmodel('./Models/PIG.Sampling');
@@ -270 +278 @@
         'none','EdgeColor','red');
 
-end
-if any(steps==7) 
-        disp('   Step 7: plot sensitivity'); 
-        
+end %}}}
+
+if any(steps==7) %{{{
+        disp('   Step 7: plot sensitivity');
+
         %load model
         md = loadmodel('./Models/PIG.Sensitivity');
@@ -280 +289 @@
 
         %which profile are we looking at?
-        index=1;        
+        index=1;
 
         %To plot sensitivities
@@ -309 +318 @@
                 'colorbartitle#3','If_{B}','unit#all','km','figure',2,...
                 'title','Importance Factors: H, \alpha, B');
-
-end
+end %}}}

issm/trunk-jpl/examples/shakti/runme.m

@@ -1 +1 @@
 steps=[1:3];
 
-if any(steps==1)
+if any(steps==1) %{{{
         disp('  Step 1: Mesh');
 
@@ -8 +8 @@
 
         save MoulinMesh md
-end
-if any(steps==2)
+end %}}}
+
+if any(steps==2) %{{{
         disp('  Step 2: Parameterization');
         md=loadmodel('MoulinMesh');
@@ -21 +22 @@
         % Change hydrology class to Sommers' SHaKTI model
         md.hydrology=hydrologyshakti();
-        
+
         % Define initial water head such that water pressure is 50% of ice overburden pressure
         md.hydrology.head = 0.5*md.materials.rho_ice/md.materials.rho_freshwater*md.geometry.thickness + md.geometry.base;
@@ -48 +49 @@
 
         save MoulinParam md;
-end
-if any(steps==3);
+end %}}}
+
+if any(steps==3) %{{{
         disp('  Step 3: Solve!');
         md=loadmodel('MoulinParam');
@@ -80 +82 @@
 
         save MoulinTransient md
-end
+end %}}}

issm/trunk-jpl/jenkins/pine_island-mac-examples

@@ -1 @@
-# NOTE: This configuration adds solid earth and Dakota capabilities to the 
-                basic build.
+# NOTE: Same configuration as pine_island-mac-full
 
 #--------------------#
@@ -6 +5 @@
 #--------------------#
 
-MATLAB_PATH="/Applications/MATLAB_R2018a.app"
+MATLAB_PATH="/Applications/MATLAB_R2019b.app"
 
 ISSM_CONFIG='\
@@ -19 +18 @@
         --with-fortran-lib="-L/usr/local/Cellar/gcc/10.2.0/lib/gcc/10 -lgfortran" \
         --with-mpi-include=${ISSM_DIR}/externalpackages/petsc/install/include \
-        --with-mpi-libflags="-L${ISSM_DIR}/externalpackages/petsc/install/lib -lmpi -lmpicxx -lmpifort" \
-        --with-blas-lapack-dir=${ISSM_DIR}/externalpackages/petsc/install \
-        --with-metis-dir=${ISSM_DIR}/externalpackages/petsc/install \
-        --with-scalapack-dir=${ISSM_DIR}/externalpackages/petsc/install \
-        --with-mumps-dir=${ISSM_DIR}/externalpackages/petsc/install \
-        --with-hdf5-dir=${ISSM_DIR}/externalpackages/petsc/install \
-        --with-petsc-dir=${ISSM_DIR}/externalpackages/petsc/install \
-        --with-gsl-dir=${ISSM_DIR}/externalpackages/gsl/install \
-        --with-boost-dir=${ISSM_DIR}/externalpackages/boost/install \
-        --with-dakota-dir=${ISSM_DIR}/externalpackages/dakota/install \
-        --with-triangle-dir=${ISSM_DIR}/externalpackages/triangle/install \
-        --with-chaco-dir=${ISSM_DIR}/externalpackages/chaco/install \
-        --with-m1qn3-dir=${ISSM_DIR}/externalpackages/m1qn3/install \
-        --with-semic-dir=${ISSM_DIR}/externalpackages/semic/install \
+        --with-mpi-libflags="-L${ISSM_EXT_SHARED_DIR}/petsc/install/lib -lmpi -lmpicxx -lmpifort" \
+        --with-blas-lapack-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \
+        --with-metis-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \
+        --with-scalapack-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \
+        --with-mumps-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \
+        --with-hdf5-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \
+        --with-petsc-dir=${ISSM_EXT_SHARED_DIR}/petsc/install \
+        --with-gsl-dir=${ISSM_EXT_SHARED_DIR}/gsl/install \
+        --with-boost-dir=${ISSM_EXT_SHARED_DIR}/boost/install \
+        --with-dakota-dir=${ISSM_EXT_SHARED_DIR}/dakota/install \
+        --with-triangle-dir=${ISSM_EXT_SHARED_DIR}/triangle/install \
+        --with-chaco-dir=${ISSM_EXT_DIR}/chaco/install \
+        --with-m1qn3-dir=${ISSM_EXT_DIR}/m1qn3/install \
+        --with-semic-dir=${ISSM_EXT_DIR}/semic/install \
 '
 
@@ -49 +48 @@
         netcdf          install-4.7-parallel.sh
         proj            install-6.2.sh
-        gdal            install-3-python-netcdf.sh
+        gdal            install-3-python.sh
         gshhg           install.sh
-        gmt                     install-6.0-mac.sh
+        gmt                     install-6-mac.sh
         gmsh            install-4.sh
         triangle        install-mac.sh
@@ -68 +67 @@
 PYTHON_TEST=0
 JAVASCRIPT_TEST=0
-EXAMPLES_TEST=0
+EXAMPLES_TEST=1
 
 # Number of CPUs used in ISSM compilation

issm/trunk-jpl/test/NightlyRun/test418.m

@@ -9 +9 @@
 %partitioning
 npart=100;
-[partition, md]=partitioner(md,'package','chaco','npart',npart);
-partition=partition-1;
+partition=partitioner(md,'package','chaco','npart',npart)-1;
 
 vector=(1:1:md.mesh.numberofvertices)';