diff --git a/aux/addanchors.m b/aux/addanchors.m
new file mode 100644
index 0000000..5225b24
--- /dev/null
+++ b/aux/addanchors.m
@@ -0,0 +1,116 @@
+%% ADDANCHORS
+% Add anchors to horizontal or vertical segments of a polygon to ensure that the segments are not too long.
+% This is useful when trying to plot the polygon on a map, as the segments are not displayed as great circle arcs.
+%
+% Syntax
+% p = addanchors(p)
+% lonlat = addanchors(lonlat)
+% [lon, lat] = addanchors(lon, lat)
+%
+% Input arguments
+% p - A polyshape object
+% lonlat - A 2-column matrix of longitudes and latitudes
+% lon, lat - Vectors of longitudes and latitudes
+%
+% Output arguments
+% p - A polyshape object
+% lonlat - A 2-column matrix of longitudes and latitudes
+% lon, lat - Vectors of longitudes and latitudes
+%
+% Notes
+% This function was originally part of the slepian_ulmo package.
+%
+% Last modified by
+% 2024/10/15, En-Chi Lee (williameclee@arizona.edu)
+
+function varargout = addanchors(varargin)
+ isPoly = false;
+
+ if nargin == 1
+
+ if isa(varargin{1}, 'polyshape')
+ isPoly = true;
+ varargin{1} = varargin{1}.Vertices;
+ varargin{1} = closecoastline(varargin{1});
+ elseif size(varargin{1}, 1) == 2
+ varargin{1} = varargin{1}';
+ elseif size(varargin{1}, 2) ~= 2
+ error('Invalid input argument')
+ end
+
+ lon = varargin{1}(:, 1);
+ lat = varargin{1}(:, 2);
+ elseif nargin == 2
+ lon = varargin{1};
+ lat = varargin{2};
+ else
+ error('Invalid number of input arguments')
+ end
+
+ lon = lon(:);
+ lat = lat(:);
+
+ maxSep = 5;
+ tol = 1;
+
+ diffLat = abs(diff(lat));
+ diffLon = abs(diff(lon));
+
+ needsAnchor = find( ...
+ (diffLat > maxSep & diffLon <= tol) | ...
+ (diffLon > maxSep & diffLat <= tol));
+ difLength = max(diffLat(needsAnchor), diffLon(needsAnchor));
+
+ lata = lat;
+ lona = lon;
+
+ if isempty(needsAnchor)
+
+ if nargout == 1
+
+ if isPoly
+ varargout = {p};
+ else
+ varargout = {[lon, lat]};
+ end
+
+ else
+ varargout = {lon, lat};
+ end
+
+ return
+ end
+
+ for i = length(needsAnchor):-1:1
+ latInterp = interp1([0, 1], ...
+ [lata(needsAnchor(i)), lata(needsAnchor(i) + 1)], ...
+ linspace(0, 1, round(difLength(i) / tol)));
+ lata = [lata(1:needsAnchor(i) - 1); latInterp'; lata(needsAnchor(i) + 1:end)];
+ lonInterp = interp1([0, 1], ...
+ [lona(needsAnchor(i)), lona(needsAnchor(i) + 1)], ...
+ linspace(0, 1, round(difLength(i) / tol)));
+ lona = [lona(1:needsAnchor(i) - 1); lonInterp'; lona(needsAnchor(i) + 1:end)];
+ end
+
+ if nargout == 0
+ figName = sprintf('Polygon with more anchors (%s)', upper(mfilename));
+ figure(998)
+ set(gcf, 'Name', figName, 'NumberTitle', 'off')
+ clf
+ plotqdm(lona, lata)
+ hold on
+ plot(lona, lata, 'k.-')
+ hold off
+ elseif nargout == 1
+
+ if isPoly
+ varargout = {polyshape(lona, lata, "Simplify", false)};
+ else
+ varargout = {[lona, lata]};
+ end
+
+ else
+ varargout = {lona, lata};
+ end
+
+end
diff --git a/oceandomains/alloceans.m b/domains/alloceans.m
similarity index 100%
rename from oceandomains/alloceans.m
rename to domains/alloceans.m
diff --git a/oceandomains/arctic.m b/domains/arctic.m
similarity index 100%
rename from oceandomains/arctic.m
rename to domains/arctic.m
diff --git a/oceandomains/atlantic.m b/domains/atlantic.m
similarity index 100%
rename from oceandomains/atlantic.m
rename to domains/atlantic.m
diff --git a/oceandomains/earthquakes.m b/domains/earthquakes.m
similarity index 100%
rename from oceandomains/earthquakes.m
rename to domains/earthquakes.m
diff --git a/oceandomains/gshhscoastline.m b/domains/gshhscoastline.m
similarity index 100%
rename from oceandomains/gshhscoastline.m
rename to domains/gshhscoastline.m
diff --git a/oceandomains/gulf.m b/domains/gulf.m
similarity index 100%
rename from oceandomains/gulf.m
rename to domains/gulf.m
diff --git a/domains/icesheetPoly.m b/domains/icesheetPoly.m
new file mode 100644
index 0000000..5f23f13
--- /dev/null
+++ b/domains/icesheetPoly.m
@@ -0,0 +1,55 @@
+%% ICESHEETPOLY - Gets an approximate ice sheet polygon
+%
+% Syntax
+% [icePoly, iceLonlat] = ICESHEETPOLY(icesheet)
+%
+% Input arguments
+% icesheet - Name of the ice sheet
+% The following ice sheets are recognised (case-insensitive):
+% - "Fennoscandian", "Fennoscandian Ice Sheet", "FIS"
+% - "Laurentide", "Laurentide Ice Sheet", "LIS"
+% - "Greenland", "Greenland Ice Sheet", "GIS"
+% - "Patagonian", "Patagonian Ice Sheet", "PIS"
+% - "Antarctic", "Antarctic Ice Sheet", "AIS"
+% - "West Antarctic", "West Antarctic Ice Sheet", "WAIS"
+% - "East Antarctic", "East Antarctic Ice Sheet", "EAIS"
+%
+% Output arguments
+% icePoly - Polyshape object of the ice sheet
+% iceLonlat - Longitude-latitude vertices of the ice sheet polygon
+%
+% Author
+% 2026/04/07, En-Chi Lee (williameclee@arizona.edu)
+
+function [icePoly, iceLonlat] = icesheetPoly(icesheet)
+
+ arguments (Input)
+ icesheet {mustBeTextScalar}
+ end
+
+ arguments (Output)
+ icePoly (1, 1) polyshape
+ iceLonlat (:, 2) double
+ end
+
+ switch lower(icesheet)
+ case {"fennoscandian", "fennoscandian ice sheet", "fis"}
+ iceLonlat = [0, 45; 30, 45; 95, 70; 95, 90; 0, 90; NaN, NaN; 360, 68; 344, 58; 344, 47; 360, 47];
+ case {"laurentide", "laurentide ice sheet", "lis"}
+ iceLonlat = [240, 90; 185, 70; 185, 50; 215, 50; 230, 30; 290, 30; 315, 45; 315, 55; 300, 60; 300, 70; 286, 77; 300, 90];
+ case {"greenland", "greenland ice sheet", "gis"}
+ iceLonlat = [300, 60; 300, 70; 286, 77; 300, 90; 355, 90; 355, 75; 340, 70; 325, 60; 315, 55];
+ case {"patagonian", "patagonian ice sheet", "pis"}
+ iceLonlat = [275, -35; 275, -60; 300, -60; 300, -35];
+ case {"antarctic", "antarctic ice sheet", "ais"}
+ iceLonlat = [0, -90; 360, -90; 360, -60; 0, -60];
+ case {"west antarctic", "west antarctic ice sheet", "wais"}
+ iceLonlat = [185, -60; 166, -78; 153, -83; 162, -85; 200, -87; 240, -85; 290, -85; 306, -83; 318, -83; 340, -74; 340, -60];
+ case {"east antarctic", "east antarctic ice sheet", "eais"}
+ iceLonlat = [0, -90; 0, -60; 185, -60; 166, -78; 153, -83; 162, -85; 200, -87; 240, -85; 290, -85; 306, -83; 318, -83; 340, -74; 340, -60; 360, -60; 360, -90];
+ otherwise
+ error("Unrecognised ice sheet name: %s", icesheet)
+ end
+
+ icePoly = polyshape(iceLonlat(:, 1), iceLonlat(:, 2));
+end
diff --git a/oceandomains/indian.m b/domains/indian.m
similarity index 100%
rename from oceandomains/indian.m
rename to domains/indian.m
diff --git a/oceandomains/natlantic.m b/domains/natlantic.m
similarity index 100%
rename from oceandomains/natlantic.m
rename to domains/natlantic.m
diff --git a/oceandomains/npacific.m b/domains/npacific.m
similarity index 100%
rename from oceandomains/npacific.m
rename to domains/npacific.m
diff --git a/oceandomains/oceans.m b/domains/oceans.m
similarity index 100%
rename from oceandomains/oceans.m
rename to domains/oceans.m
diff --git a/oceandomains/pacific.m b/domains/pacific.m
similarity index 100%
rename from oceandomains/pacific.m
rename to domains/pacific.m
diff --git a/oceandomains/satlantic.m b/domains/satlantic.m
similarity index 100%
rename from oceandomains/satlantic.m
rename to domains/satlantic.m
diff --git a/oceandomains/spacific.m b/domains/spacific.m
similarity index 100%
rename from oceandomains/spacific.m
rename to domains/spacific.m
diff --git a/externalio/selen/convertSelenModel.m b/externalio/selen/convertSelenModel.m
new file mode 100644
index 0000000..679bdc6
--- /dev/null
+++ b/externalio/selen/convertSelenModel.m
@@ -0,0 +1,72 @@
+%% CONVERTSELENMODEL - Converts Selen outputs to the format used by SLEPIAN
+% Converts the geoid and vertical displacement outputs from the SELEN model
+% to the format used by SLEPIAN_DELTA. This is not intended to be called
+% directly by users.
+%
+% Syntax
+% convertSelenModel(modelFolder)
+%
+% Input arguments
+% modelFolder - folder containing the SELEN model outputs
+%
+% References
+% Spada, Giorgio & Melini, Daniele. (2019). SELEN4 (SELEN version 4.0): a
+% Fortran program for solving the gravitationally and topographically
+% self-consistent sea-level equation in glacial isostatic adjustment
+% modeling. Geosci. Model Dev.
+% https://doi.org/10.5194/gmd-12-5055-2019
+% https://github.com/geodynamics/selen
+% Tamisiea, Mark E. (2011). Ongoing glacial isostatic contributions to
+% observations of sea level change. Geophys. J. Int.
+% https://doi.org/10.1111/j.1365-246X.2011.05116.x
+%
+% Author
+% 2026/04/01, En-Chi Lee (williameclee@arizona.edu)
+
+function convertSelenModel(modelFolder)
+
+ arguments (Input)
+ modelFolder {mustBeFolder}
+ end
+
+ dataFolder = fullfile(modelFolder, "FPR");
+ % Detect the degree of the model from the model folder name
+ [~, modelFolderSelf] = fileparts(modelFolder);
+ modelName = regexp(modelFolderSelf, "RUN_([A-Za-z0-9-_]+)", "tokens", "once");
+ modelInfo = regexp(modelFolderSelf, "RUN_([A-Za-z0-9]+)-R(\d+)-L(\d+)-I(\d+)", "tokens", "once");
+
+ % Prepare grid for interpolation
+ L = str2double(modelInfo{3});
+ h = 360 / (2 * L);
+ lon = 0:h:360;
+ lat = -90:h:90;
+ [lonn, latt] = meshgrid(lon, lat);
+ %% Geoid and surface density
+ geoidPath = fullfile(dataFolder, "gdot.pix");
+ geoidXyz = readmatrix(geoidPath, "FileType", "text");
+ geoidXyz(:, 3) = geoidXyz(:, 3) / 1e3; % mm -> m
+ F = scatteredInterpolant(geoidXyz(:, 1), geoidXyz(:, 2), geoidXyz(:, 3), "nearest", "nearest");
+ geoidMesh = F(mod(lonn, 360), latt);
+ geoidPlm = xyz2plm_new(flip(geoidMesh), L);
+
+ geoidPlm(5, 3:4) = geoidPlm(5, 3:4) / 2.06; % See Appendix of Tamisiea (2011)
+
+ sdPlm = convertgravity(geoidPlm, "POT", "SD");
+
+ % Save the converted model
+ lmcosiM = sdPlm;
+ save(fullfile(modelFolder, sprintf("%s_SD.mat", modelName)), "lmcosiM");
+ lmcosiM = geoidPlm;
+ save(fullfile(modelFolder, sprintf("%s_POT.mat", modelName)), "lmcosiM");
+
+ %% Vertical displacement
+ vlmPath = fullfile(dataFolder, "udot.pix");
+ vlmXyz = readmatrix(vlmPath, "FileType", "text");
+ F = scatteredInterpolant(vlmXyz(:, 1), vlmXyz(:, 2), vlmXyz(:, 3), "nearest", "nearest");
+ vlmMesh = F(mod(lonn, 360), latt);
+ vlmPlm = xyz2plm_new(flip(vlmMesh), L);
+
+ % Save the converted model
+ lmcosiM = vlmPlm;
+ save(fullfile(modelFolder, sprintf("%s_VLM.mat", modelName)), "lmcosiM");
+end
diff --git a/externalio/selen/modSelenIceModel.m b/externalio/selen/modSelenIceModel.m
new file mode 100644
index 0000000..b53f87e
--- /dev/null
+++ b/externalio/selen/modSelenIceModel.m
@@ -0,0 +1,141 @@
+%% MODSELENICEMODEL - Modifies a ice model in SELEN format
+% This function modifies the ice model used in a SELEN run by changing the
+% volume change and time delay of a specified ice sheet. The modified model
+% is saved as a new file in the same format as the original model. This
+% function is intended to be used for sensitivity tests of the ice model in
+% SELEN.
+%
+% Syntax
+% modSelenIceModel(model, icesheet, newModel, volumeChangeFrac, timeDelay)
+%
+% Input arguments
+% model - Name of the SELEN model to modify (without the .pix extension)
+% The model must be located in the DATA folder of the SELEN
+% installation.
+% icesheet (optional) - Name of the ice sheet to modify
+% See ICESHEETPOLY for recognised ice sheet names.
+% The default ice sheet to modify is "fennoscandian".
+% newModel (optional) - Name of the new model to save
+% If empty, the new name is programmatically generated based on the
+% original model name and the modifications made.
+% volumeChangeFrac (optional) - Fractional change in the volume of the
+% specified ice sheet
+% For example, a value of 0.1 means a 10% increase in volume, while a
+% value of -0.2 means a 20% decrease in volume.
+% The default value is 0, meaning no change in volume.
+% timeDelay (optional) - Time delay in years for the specified ice sheet
+% A positive value means the ice sheet changes later than in the
+% original model, while a negative value means it changes earlier.
+% The default value is 0, meaning no time delay.
+%
+% See also
+% ICESHEETPOLY
+%
+% References
+% Spada, Giorgio & Melini, Daniele. (2019). SELEN4 (SELEN version 4.0): a
+% Fortran program for solving the gravitationally and topographically
+% self-consistent sea-level equation in glacial isostatic adjustment
+% modeling. Geosci. Model Dev.
+% https://doi.org/10.5194/gmd-12-5055-2019
+%
+% Author
+% 2026/04/07, En-Chi Lee (williameclee@arizona.edu)
+
+function modSelenIceModel(model, icesheet, newModel, volumeChangeFrac, timeDelay)
+
+ arguments (Input)
+ model {mustBeTextScalar} = "i6g-R44r-i"
+ icesheet {mustBeTextScalar} = "fennoscandian"
+ newModel {mustBeTextScalar} = ""
+ volumeChangeFrac (1, 1) double {mustBeReal, mustBeFinite} = 0
+ timeDelay (1, 1) double {mustBeInteger, mustBeReal, mustBeFinite} = 0
+ end
+
+ modelPath = fullfile(getenv("SELEN"), "DATA", sprintf("%s.pix", model));
+ icesheet = icesheetCode(icesheet);
+
+ ice = readmatrix(modelPath, "FileType", "text");
+ isIcesheet = labelIceSheet(ice(:, 2:3), icesheet);
+
+ iceDiff = ice(:, 5:end) - ice(:, end);
+
+ if volumeChangeFrac ~= 0
+ iceDiff(isIcesheet, :) = iceDiff(isIcesheet, :) * (1 + volumeChangeFrac);
+ end
+
+ if timeDelay > 0
+ iceDiff(isIcesheet, timeDelay + 1:end) = iceDiff(isIcesheet, 1:end - timeDelay);
+
+ iceDiff(isIcesheet, 1:timeDelay) = repmat(iceDiff(isIcesheet, 1), [1, timeDelay]);
+
+ elseif timeDelay < 0
+ iceDiff(isIcesheet, 1:end + timeDelay) = iceDiff(isIcesheet, -timeDelay + 1:end);
+ end
+
+ iceNew = [ice(:, 1:4), ice(:, end) + iceDiff];
+
+ if ~isempty(newModel)
+ newModel = sprintf("%s-%s_V%s_T%s", newModel, upper(icesheet), ...
+ replace(replace(sprintf("%+03.0f", volumeChangeFrac * 100), "+", "p"), "-", "n"), ...
+ replace(replace(sprintf("%+03d", timeDelay), "+", "p"), "-", "n"));
+ end
+
+ newModelPath = fullfile(getenv("SELEN"), "DATA", newModel);
+
+ fid = fopen(newModelPath, "w");
+ cleanUp = onCleanup(@() fclose(fid));
+
+ for i = 1:size(iceNew, 1)
+ fprintf(fid, "%12d %20f %20f %25f ", iceNew(i, 1:4));
+ fprintf(fid, " %11.0f", iceNew(i, 5:end));
+ fprintf(fid, "\n");
+ end
+
+end
+
+%% Subfunctions
+function code = icesheetCode(name)
+
+ arguments (Input)
+ name {mustBeTextScalar}
+ end
+
+ arguments (Output)
+ code (1, 1) string
+ end
+
+ switch lower(name)
+ case {"fennoscandian", "fennoscandian ice sheet", "fis"}
+ code = "FIS";
+ case {"laurentide", "laurentide ice sheet", "lis"}
+ code = "LIS";
+ case {"greenland", "greenland ice sheet", "gis"}
+ code = "GIS";
+ case {"patagonian", "patagonian ice sheet", "pis"}
+ code = "PIS";
+ case {"antarctic", "antarctic ice sheet", "ais"}
+ code = "AIS";
+ case {"west antarctic", "west antarctic ice sheet", "wais"}
+ code = "WAIS";
+ case {"east antarctic", "east antarctic ice sheet", "eais"}
+ code = "EAIS";
+ otherwise
+ error("Unrecognised ice sheet name: %s", name)
+ end
+
+end
+
+function isIcesheet = labelIceSheet(lonlat, icesheet)
+
+ arguments (Input)
+ lonlat (:, 2) {mustBeNumeric}
+ icesheet {mustBeTextScalar}
+ end
+
+ arguments (Output)
+ isIcesheet (:, 1) logical
+ end
+
+ isIcesheet = isinterior(icesheetPoly(icesheet), lonlat(:, 1), lonlat(:, 2));
+
+end
diff --git a/gia2plmt.m b/gia2plmt.m
index 0a62b59..14e1f6a 100644
--- a/gia2plmt.m
+++ b/gia2plmt.m
@@ -17,25 +17,26 @@
%
% Input arguments
% model - Name of the GIA model
-% - 'Steffen_ice6g_vm5a': A model computed by H. Steffen using the
+% - 'Steffen_ice6g_vm5a': A model computed by H. Steffen using the
% ice6g ice model and vm5a viscosity profile. Other models from
% this dataset are also available and use the original naming
% scheme. For example, 'Steffen_anu-ice_i72'.
% This family of models can also be specified as a cell array,
% e.g. {'Steffen', 'ice6g', 'vm5a'}.
-% - 'LM17.3': A model based on the data from the LM17.3 dataset.
-% - 'Paulson07': A model based on the ICE-5G ice load model of
+% - 'LM17.3': A model based on the data from the LM17.3 dataset.
+% - 'SELEN-modelname': A model from the SELEN software.
+% - 'Paulson07': A model based on the ICE-5G ice load model of
% Peltier (2004). Suitable for both Antarctica and Greenland. As
% corrected by Geruo A and J. Wahr.
% Please avoid using this model for oceans.
-% - 'Wangetal08': A model based on the older ICE-4G ice model, and
+% - 'Wangetal08': A model based on the older ICE-4G ice model, and
% viscosity which varies laterally. Suitable for Greenland.
-% - 'IJ05_R2': A model based on the Ivins et al. (2013) ice model.
+% - 'IJ05_R2': A model based on the Ivins et al. (2013) ice model.
% Suitable for Antarctica.
-% - 'IJ05': A model based on the Ivins and James (2005) ice model.
+% - 'IJ05': A model based on the Ivins and James (2005) ice model.
% Suitable for Antarctica.
-% - 'W12a_v1': A 'best' model from Whitehouse et al. (2012). Suitable
-% only for Antarctica.
+% - 'W12a_v1': A 'best' model from Whitehouse et al. (2012).
+% Suitable only for Antarctica.
% The input can also be the path to the model file.
% The default model is 'Steffen_ice6g_vm5a'.
% days - Number of days to calculate the GIA change for
@@ -55,8 +56,8 @@
% The default option is false.
%
% Output arguments
-% Plmt - GIA change in spherical harmonic format
-% PlmtU, PlmtL - Upper and lower bounds of the GIA change, if available
+% plmt - GIA change in spherical harmonic format
+% plmtU, plmtL - Upper and lower bounds of the GIA change, if available
%
% See also
% CORRECT4GIA
@@ -95,7 +96,9 @@
% PANGAEA, doi: 10.1594/PANGAEA.932462
%
% Last modified by
-% 2025/11/16, williameclee@arizona.edu (@williameclee)
+% 2026/04/01, En-Chi Lee (williameclee@arizona.edu)
+% - Added support for models from the SELEN software
+% 2025/11/16, En-Chi Lee (williameclee@arizona.edu)
function varargout = gia2plmt(varargin)
%% Initialisation
@@ -117,7 +120,7 @@
data = load(inputPath, 'lmcosiM', 'lmcosiU', 'lmcosiL');
if ~beQuiet
- fprintf('[ULMO>%s] Loaded %s GIA model\n', ...
+ fprintf('[Slepian>%s] Loaded %s GIA model\n', ...
callchaintext(callChain), inputPath, inputPath, upper(model));
end
@@ -134,13 +137,14 @@
if ~exist(altInputPath, 'file')
error('ULMO:LoadData:FileNotFound', ...
- 'GIA model %s not found at expected location:\n%s or %s', upper(model), inputPath, altInputPath);
+ 'GIA model %s not found at expected location:\n%s or %s', ...
+ upper(model), inputPath, altInputPath);
end
data = load(altInputPath, 'lmcosiM', 'lmcosiU', 'lmcosiL');
if ~beQuiet
- fprintf('[ULMO>%s] Loaded %s GIA model\n', ...
+ fprintf('[Slepian>%s] Loaded %s GIA model\n', ...
callchaintext(callChain), inputPath, inputPath, upper(model));
end
@@ -217,7 +221,7 @@
end
if ~beQuiet
- fprintf('[ULMO>%s] Saved %s GIA model\n', ...
+ fprintf('[Slepian>%s] Saved %s GIA model\n', ...
callchaintext(callChain), inputPath, inputPath, upper(model));
end
@@ -298,7 +302,7 @@
GIAtL = [];
if nargout > 1 && ~beQuiet
- warning(sprintf('ULMO:%s:NoBoundsToReturn', upper(mfilename)), ...
+ warning('ULMO:gia2plmt:NoBoundsToReturn', ...
'Upper and lower bounds are not available for model %s', upper(model));
end
@@ -385,8 +389,7 @@
end
if L ~= round(L)
- warning( ...
- sprintf('ULMO:%s:InvalidInput:NonIntegerDegree', upper(mfilename)), ...
+ warning('ULMO:gia2plmt:InvalidInput:NonIntegerDegree', ...
'The degree L must be an integer. Rounding to %d', round(L));
L = round(L);
end
@@ -417,10 +420,16 @@
error('GIA folder not found')
end
+ modelSource = "NaN";
+
if strncmp(model, 'Morrow', 6)
inputFolder = fullfile(inputFolder, model(1:6));
elseif strncmpi(model, 'Steffen', 7)
inputFolder = fullfile(inputFolder, 'Steffen21');
+ elseif strncmpi(model, 'selen-', 6) || strncmpi(model, 'selen_', 6)
+ modelSource = "selen";
+ model = regexprep(model, 'selen[-_]', '', 'ignorecase');
+ inputFolder = fullfile(inputFolder, 'Selen', sprintf("RUN_%s", model));
elseif strcmp(model, 'LM17.3')
inputFolder = fullfile(inputFolder, 'LM17.3');
@@ -434,6 +443,11 @@
% And the appropriate name
inputPath = fullfile(inputFolder, sprintf('%s_%s.mat', model, unit));
+
+ if ~exist(inputPath, 'file') && strcmpi(modelSource, "selen")
+ convertSelenModel(inputFolder);
+ end
+
end
function plmt = plm2plmt(plm, deltaYear)
diff --git a/giaz2plmt.m b/giaz2plmt.m
index f1a5ddf..06c71cf 100644
--- a/giaz2plmt.m
+++ b/giaz2plmt.m
@@ -14,11 +14,11 @@
%
% Input arguments
% model - Name of the GIA model
-% - Name of a model computed by H. Steffen. It can be specified as,
+% - Name of a model computed by H. Steffen. It can be specified as,
% e.g. 'Steffen_ice6g_vm5a' or {'Steffen', 'ice6g', 'vm5a'}.
-% - LM17.3 is also supported.
+% - LM17.3 is also supported.
% Other models are specified in the same way.
-% - The input can also be the path to the model file.
+% - The input can also be the path to the model file.
% The default model is 'Steffen_ice6g_vm5a'.
% Format: string or 1 x 3 cell array.
% years - Number of years to calculate the GIA change for
@@ -54,7 +54,9 @@
% PANGAEA, doi: 10.1594/PANGAEA.932462
%
% Last modified by
-% 2025/08/03, williameclee@arizona.edu (@williameclee)
+% 2026/04/01, En-Chi Lee (williameclee@arizona.edu)
+% - Added support for models from the SELEN software
+% 2025/08/03, En-Chi Lee (williameclee@arizona.edu)
function varargout = giaz2plmt(varargin)
%% Initialisation
@@ -86,8 +88,21 @@
wPlmt = wSph;
wUPlmt = [];
wLPlmt = [];
+ elseif strncmpi(model, 'selen-', 6) || strncmpi(model, 'selen_', 6)
+ model = regexprep(model, 'selen[-_]', '', 'ignorecase');
+ inputFolder = fullfile(getenv('IFILES'), 'GIA', 'Selen', sprintf("RUN_%s", model));
+ inputPath = fullfile(inputFolder, sprintf('%s_VLM.mat', model));
+
+ if ~exist(inputPath, 'file')
+ convertSelenModel(inputFolder);
+ end
+
+ load(inputPath, 'lmcosiM');
+ wPlmt = lmcosiM;
+ wUPlmt = [];
+ wLPlmt = [];
else
- error('ULMO:LoadData:FileNotFound', 'Unrecognised model name %s', upper(model));
+ error('Slepian:LoadData:FileNotFound', 'Unrecognised model name %s', upper(model));
end
if (size(wPlmt, 1) < addmup(L) || ...
@@ -238,7 +253,7 @@
% Make sure the file exists
if exist(inputPath, 'file') ~= 2
- error('ULMO:LoadData:FileNotFound', ...
+ error('Slepian:LoadData:FileNotFound', ...
'Model %s not found at %s', upper(model), inputPath);
end
diff --git a/plotting/formatlatticks.m b/plotting/formatlatticks.m
new file mode 100644
index 0000000..3a14ff1
--- /dev/null
+++ b/plotting/formatlatticks.m
@@ -0,0 +1,114 @@
+%% FORMATLATTICKS - Format latitude tick labels for a plot.
+% The tick labels will always have the same number of significant digits.
+%
+% Syntax
+% formatlatticks
+% Formats the y-axis tick labels of the current plot as latitude
+% ticks.
+% formatlatticks(lat)
+% Formats the y-axis tick labels of the current plot as latitude
+% ticks with the specified values.
+% formatlatticks(ax)
+% Formats the y-axis tick labels of the specified axes as latitude
+% ticks.
+% latL = formatlatticks(__)
+% Returns the formatted latitude tick labels as a cell array of
+% character vectors.
+%
+% Input arguments
+% lat - Latitude ticks in degrees
+% ax - Axes object to format
+%
+% Output arguments
+% latL - Formatted latitude tick labels
+%
+% Last modified by
+% 2024/10/11, williameclee@arizona.edu (@williameclee)
+
+function varargout = formatlatticks(varargin)
+
+ if nargin > 0 && isa(varargin{1}, 'matlab.graphics.axis.Axes')
+ ax = varargin{1};
+ varargin(1) = [];
+ else
+ ax = nan;
+ end
+
+ ip = inputParser;
+ addOptional(ip, 'lat', [], @(x) isnumeric(x) || iscell(x));
+ parse(ip, varargin{:});
+ lat = ip.Results.lat;
+
+ if iscell(lat)
+ lat = cell2mat(lat);
+ elseif isempty(lat)
+
+ if isnumeric(ax) && isnan(ax)
+ lat = yticks;
+ else
+ lat = yticks(ax);
+ end
+
+ end
+
+ if isnumeric(ax) && isnan(ax)
+ yticks(lat)
+ else
+ yticks(ax, lat)
+ end
+
+ if isscalar(lat)
+ latL = formatlattick(lat, nan);
+ elseif isnumeric(lat) || iscell(lat)
+
+ if iscell(lat)
+ lat = cell2mat(lat);
+ end
+
+ minDiff = min(diff(lat));
+
+ if minDiff >= 1
+ sigNum = 0;
+ else
+ sigNum = ceil(-log10(minDiff));
+ end
+
+ latL = arrayfun(@(x) formatlattick(x, sigNum), lat, 'UniformOutput', false);
+ else
+ error('Invalid input type %s', class(lat));
+ end
+
+ if nargout > 0
+ varargout = {latL};
+ return
+ end
+
+ clear varargout
+
+ if isnumeric(ax) && isnan(ax)
+ yticklabels(latL)
+ else
+ yticklabels(ax, latL)
+ end
+
+end
+
+%% Subfunctions
+function latL = formatlattick(lat, sigNum)
+ sgn = sign(lat);
+
+ if ~isnan(sigNum)
+ lat = sprintf(['%0.', num2str(sigNum), 'f'], abs(lat));
+ end
+
+ latL = num2str(lat);
+
+ if sgn > 0
+ latL = [latL, '°N'];
+ elseif sgn < 0
+ latL = [latL, '°S'];
+ else
+ latL = '0°';
+ end
+
+end
diff --git a/plotting/formatlonticks.m b/plotting/formatlonticks.m
new file mode 100644
index 0000000..10afa25
--- /dev/null
+++ b/plotting/formatlonticks.m
@@ -0,0 +1,114 @@
+%% FORMATLONTICKS - Format longitude tick labels for a plot
+%
+% Syntax
+% formatlonticks
+% Formats the x-axis tick labels of the current plot as longitude
+% ticks.
+% formatlonticks(lon)
+% Formats the x-axis tick labels of the current plot as longitude
+% ticks with the specified values.
+% formatlonticks(ax)
+% Formats the x-axis tick labels of the specified axes as longitude
+% ticks.
+% lonL = formatlonticks(__)
+% Returns the formatted longitude tick labels as a cell array of
+% character vectors.
+%
+% Input arguments
+% lon - Longitude ticks in degrees
+% ax - Axes object to format
+%
+% Output arguments
+% lonL - Formatted longitude tick labels
+%
+% Last modified by
+% 2024/10/25, williameclee@arizona.edu (@williameclee)
+
+function varargout = formatlonticks(varargin)
+
+ if nargin > 0 && isa(varargin{1}, 'matlab.graphics.axis.Axes')
+ ax = varargin{1};
+ varargin(1) = [];
+ else
+ ax = nan;
+ end
+
+ ip = inputParser;
+ addOptional(ip, 'lon', [], @(x) isnumeric(x) || iscell(x));
+ parse(ip, varargin{:});
+ lon = ip.Results.lon;
+
+ if iscell(lon)
+ lon = cell2mat(lon);
+ elseif isempty(lon)
+
+ if isnumeric(ax) && isnan(ax)
+ lon = xticks;
+ else
+ lon = xticks(ax);
+ end
+
+ end
+
+ if isnumeric(ax) && isnan(ax)
+ xticks(lon)
+ else
+ xticks(ax, lon)
+ end
+
+ if isscalar(lon)
+ lonL = formatlontick(lon, nan);
+ elseif isnumeric(lon) || iscell(lon)
+
+ if iscell(lon)
+ lon = cell2mat(lon);
+ end
+
+ minDiff = min(diff(lon));
+
+ if minDiff >= 1
+ sigNum = 0;
+ else
+ sigNum = ceil(-log10(minDiff));
+ end
+
+ lonL = arrayfun(@(x) formatlontick(x, sigNum), lon, 'UniformOutput', false);
+ else
+ error('Invalid input type %s', class(lon));
+ end
+
+ if nargout > 0
+ varargout = {lonL};
+ return
+ end
+
+ clear varargout
+
+ if isnumeric(ax) && isnan(ax)
+ xticklabels(lonL)
+ else
+ xticklabels(ax, lonL)
+ end
+
+end
+
+%% Subfunctions
+function lonL = formatlontick(lon, sigNum)
+
+ lon = mod(lon + 180, 360) - 180;
+
+ if ~isnan(sigNum)
+ lonL = sprintf(['%0.', num2str(sigNum), 'f'], abs(lon));
+ end
+
+ lonL = num2str(lonL);
+
+ if lon > 0
+ lonL = [lonL, '°E'];
+ elseif lon < 0
+ lonL = [lonL, '°W'];
+ else
+ lonL = '0°';
+ end
+
+end
diff --git a/plotting/parselonlatinputs.m b/plotting/parselonlatinputs.m
new file mode 100644
index 0000000..2c250a4
--- /dev/null
+++ b/plotting/parselonlatinputs.m
@@ -0,0 +1,85 @@
+%% PARSELONLATINPUTS - Parses inputs for any longitude and latitude inputs
+%
+% Syntax
+% [lon, lat, __, type] = parselonlatinputs(lon, lat, __)
+% [lon, lat, __, type] = parselonlatinputs(lonlat, __)
+% [lon, lat, __, type] = parselonlatinputs(p, __)
+% [lon, lat, __, type] = parselonlatinputs(geodomain, __)
+%
+% Input arguments
+% lon, lat - Numeric array of longitudes and latitudes
+% The two arrays must have the same size.
+% lonlat - N-by-2 (or 2-by-N) numeric array of longitude and latitude
+% coordinates
+% p - polyshape object
+% geodomain - GeoDomain object
+% Any other arguments are not processed and passed back to the caller.
+%
+% Output arguments
+% lon, lat - numeric array of longitudes and latitudes
+% type - which type of input was used
+% 0 - lon and lat
+% 1 - polyshape
+% 2 - GeoDomain
+%
+% Author
+% 2024/08/12, En-Chi Lee (williameclee@gmail.com)
+%
+% Last modified by
+% 2026/04/01, En-Chi Lee (williameclee@arizona.edu)
+% - Updated to documentation
+% 2024/10/18, En-Chi Lee (williameclee@gmail.com)
+
+function varargout = parselonlatinputs(varargin)
+ % Find longitude and latitude
+ if isnumeric(varargin{1})
+ typeflag = 0;
+
+ if length(varargin) >= 2 && isnumeric(varargin{2}) && ...
+ isequal(size(varargin{1}), size(varargin{2}))
+ % First input is lon, second is lat
+ lon = varargin{1};
+ lat = varargin{2};
+
+ % Remove lonlat from varargin
+ varargin(1:2) = [];
+
+ elseif any(size(varargin{1}) == 2)
+ % First input is lonlat
+ lonlat = varargin{1};
+
+ if size(lonlat, 1) == 2 && size(lonlat, 2) ~= 2
+ lonlat = lonlat';
+ end
+
+ lon = lonlat(:, 1);
+ lat = lonlat(:, 2);
+
+ varargin(1) = [];
+
+ else
+ error('Invalid input argument size')
+ end
+
+ elseif isa(varargin{1}, 'polyshape')
+ typeflag = 1;
+ % First input is a polyshape
+ [lon, lat] = boundary(varargin{1});
+ varargin(1) = [];
+ elseif isa(varargin{1}, 'GeoDomain')
+ % This class is found in the slepian_ulmo package
+ typeflag = 2;
+ % First input is a GeoDomain
+ domain = varargin{1};
+ lonlat = domain.Lonlat;
+ lon = lonlat(:, 1);
+ lat = lonlat(:, 2);
+ varargin(1) = [];
+ else
+ error('Invalid input argument type for the first argument')
+ end
+
+ typeflag = uint8(typeflag);
+
+ varargout = {lon, lat, varargin, typeflag};
+end
diff --git a/solvedegree1.m b/solvedegree1.m
index 628d110..cbeb474 100644
--- a/solvedegree1.m
+++ b/solvedegree1.m
@@ -96,11 +96,13 @@
% See also
% SOLVESLE, GRACE2PLMT
%
-% Authored by
-% 2025/03/18, williameclee@arizona.edu (@williameclee)
+% Author
+% 2025/03/18, En-Chi Lee (williameclee@arizona.edu)
%
% Last modified by
-% 2025/10/01, williameclee@arizona.edu (@williameclee)
+% 2026/04/01, En-Chi Lee (williameclee@arizona.edu)
+% - Made sure the waitbar is properly deleted
+% 2025/10/01, En-Chi Lee (williameclee@arizona.edu)
function varargout = solvedegree1(varargin)
%% Initialisation
@@ -163,6 +165,7 @@
%% Loading data
wbar = waitbar(0, 'Loading GRACE data', ...
"Name", upper(mfilename), "CreateCancelBtn", 'setappdata(gcbf,''canceling'',1)');
+ cleanup = onCleanup(@() delete(wbar));
[gracePlmt, graceStdPlmt, dates] = grace2plmt_new(pcenter, rlevel, Ldata, ...
"Unit", 'SD', "OutputFormat", 'timefirst', "TimeFormat", 'datetime', ...
diff --git a/solvesle.m b/solvesle.m
index 40f183a..f63c284 100644
--- a/solvesle.m
+++ b/solvesle.m
@@ -86,11 +86,13 @@
% by cryosphere and climate driven mass change. Geoscientific Model
% Development, 9(3):1087–1109. doi: 10.5194/gmd-9-1087-2016
%
-% Authored by
-% 2024/11/20, williameclee@arizona.edu (@williameclee)
+% Author
+% 2024/11/20, En-Chi Lee (williameclee@arizona.edu)
%
% Last modified by
-% 2026/02/12, williameclee@arizona.edu (@williameclee)
+% 2026/04/01, En-Chi Lee (williameclee@arizona.edu)
+% - Made sure the waitbar is properly deleted
+% 2026/02/12, En-Chi Lee (williameclee@arizona.edu)
% - Added support for ALMA3 love numbers
function varargout = solvesle(varargin)
@@ -107,6 +109,7 @@
if ~beQuiet
wbar = waitbar(0, 'Initialising', ...
"Name", upper(mfilename), "CreateCancelBtn", 'setappdata(gcbf,''canceling'',1)');
+ cleanup = onCleanup(@() delete(wbar));
end
% Load the ocean function if not provided