Skip to content

GitLab

fvcom-toolbox
fvcom-toolbox
Commit d59b9c5a authored by Pierre Cazenave's avatar Pierre Cazenave
Browse files
Options

Remove the old POLCOMS versions of the scripts and replace with those designed... 

Remove the old POLCOMS versions of the scripts and replace with those designed against the daily mean files from the AMM domain. Adjust the comments and variable names accordingly
parent 9ef1a54e
Expand all Show whitespace changes
Inline Side-by-side
Showing with 99 additions and 476 deletions
fvcom_prepro/get_AMM_tsobc.m deleted 100644 → 0
View file @ 9ef1a54e
function Mobj = get_AMM_tsobc(Mobj, ts)
% Read temperature and salinity from the PML POLCOMS-ERSEM NetCDF model
% output files and interpolate onto the open boundaries in Mobj.
%
% function Mobj = get_AMM_tsobc(Mobj, ts, polcoms_bathy, varlist)
%
% DESCRIPTION:
% Interpolate temperature and salinity values onto the FVCOM open
% boundaries at all sigma levels.
%
% INPUT:
% Mobj = MATLAB mesh structure which must contain:
% - Mobj.siglayz - sigma layer depths for all model nodes.
% - Mobj.lon, Mobj.lat - node coordinates (lat/long).
% - Mobj.obc_nodes - list of open boundary node inidices.
% - Mobj.nObcNodes - number of nodes in each open boundary.
% ts = Cell array of AMM NetCDF file(s) in which 4D
% variables of temperature and salinity (called 'ETWD' and
% 'x1XD') exist. Their shape should be (y, x, sigma, time).
%
% NOTES:
%
% - If you supply multiple files in ts, there are a few assumptions:
%
% - Variables are only appended if there are 4 dimensions; fewer than
% that, and the values are assumed to be static across all the given
% files (e.g. longitude, latitude etc.). The fourth dimension is
% time.
% - The order of the files given should be chronological.
%
% - The NetCDF files used here are those from the PML POLCOMS-ERSEM model
% output.
%
% OUTPUT:
% Mobj = MATLAB structure in which three new fields (called temperature,
% salinity and ts_time). temperature and salinity have sizes
% (sum(Mobj.nObcNodes), sigma, time). The time dimension is
% determined based on the input NetCDF file. The ts_time variable
% is just the input file times in Modified Julian Day.
%
% EXAMPLE USAGE
% Mobj = get_AMM_forcing(Mobj, ts, depth)
%
% Author(s):
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
% 2013-02-07 First version based on get_POLCOMS_tsobc.m.
%
%==========================================================================
subname = 'get_AMM_tsobc';
global ftbverbose;
if ftbverbose
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end
varlist = {'lon', 'lat', 'ETWD', 'x1XD', 'time', 'depth', 'pdepthD'};
% Data format:
%
% amm.ETWD.data and amm.x1XD.data are y, x, sigma, time
%
amm = get_AMM_netCDF(ts, varlist);
[~, ~, nz, nt] = size(amm.ETWD.data);
% Make rectangular arrays for the nearest point lookup.
[lon, lat] = meshgrid(amm.lon.data, amm.lat.data);
fvlon = Mobj.lon(Mobj.obc_nodes(Mobj.obc_nodes ~= 0));
fvlat = Mobj.lat(Mobj.obc_nodes(Mobj.obc_nodes ~= 0));
% Number of boundary nodes
nf = sum(Mobj.nObcNodes);
% Number of sigma layers.
fz = size(Mobj.siglayz, 2);
fvtemp = nan(nf, fz, nt); % FVCOM interpolated temperatures
fvsal = nan(nf, fz, nt); % FVCOM interpolated salinities
if ftbverbose
tic
end
for t = 1:nt
% Get the current 3D array of AMM results.
ammtemp3 = amm.ETWD.data(:, :, :, t);
ammsalt3 = amm.x1XD.data(:, :, :, t);
% Preallocate the intermediate results arrays.
itempz = nan(nf, nz);
isalz = nan(nf, nz);
idepthz = nan(nf, nz);
for j = 1:nz
% Now extract the relevant layer from the 3D subsets. Transpose the
% data to be (x, y) rather than (y, x).
ammtemp2 = ammtemp3(:, :, j)';
ammsalt2 = ammsalt3(:, :, j)';
ammdepth2 = squeeze(amm.depth.data(:, :, j, t))';
% Create new arrays which will be flattened when masking (below).
tammtemp2 = ammtemp2;
tammsalt2 = ammsalt2;
tammdepth2 = ammdepth2;
tlon = lon;
tlat = lat;
% Create and apply a mask to remove values outside the domain. This
% inevitably flattens the arrays, but it shouldn't be a problem
% since we'll be searching for the closest values in such a manner
% as is appropriate for an unstructured grid (i.e. we're assuming
% the AMM data is irregularly spaced).
mask = tammdepth2 > 20000;
tammtemp2(mask) = [];
tammsalt2(mask) = [];
tammdepth2(mask) = [];
% Also apply the masks to the position arrays so we can't even find
% positions outside the domain, effectively meaning if a value is
% outside the domain, the nearest value to the boundary node will
% be used.
tlon(mask) = [];
tlat(mask) = [];
% Preallocate the intermediate results arrays.
itempobc = nan(nf, 1);
isalobc = nan(nf, 1);
idepthobc = nan(nf, 1);
% Speed up the tightest loop with a parallelized loop.
parfor i = 1:nf
% Now we can do each position within the 2D layer.
fx = fvlon(i);
fy = fvlat(i);
[~, ii] = sort(sqrt((tlon - fx).^2 + (tlat - fy).^2));
% Get the n nearest nodes from AMM (more? fewer?).
ixy = ii(1:16);
% Get the variables into static variables for the
% parallelisation.
plon = tlon(ixy);
plat = tlat(ixy);
ptemp = tammtemp2(ixy);
psal = tammsalt2(ixy);
pdepth = tammdepth2(ixy);
% Use a triangulation to do the horizontal interpolation.
tritemp = TriScatteredInterp(plon', plat', ptemp', 'natural');
trisal = TriScatteredInterp(plon', plat', psal', 'natural');
triz = TriScatteredInterp(plon', plat', pdepth', 'natural');
itempobc(i) = tritemp(fx, fy);
isalobc(i) = trisal(fx, fy);
idepthobc(i) = triz(fx, fy);
% Check all three, though if one is NaN, they all will be.
if isnan(itempobc(i)) || isnan(isalobc(i)) || isnan(idepthobc(i))
warning('FVCOM boundary node at %f, %f is outside the AMM domain. Setting to the closest AMM value.', fx, fy)
itempobc(i) = tammtemp2(ii(1));
isalobc(i) = tammsalt2(ii(1));
idepthobc(i) = tammdepth2(ii(1));
end
end
% Put the results in this intermediate array.
itempz(:, j) = itempobc;
isalz(:, j) = isalobc;
idepthz(:, j) = idepthobc;
end
% Now we've interpolated in space, we can interpolate the z-values
% to the sigma depths.
oNodes = Mobj.obc_nodes(Mobj.obc_nodes ~= 0);
for zi = 1:fz
% Preallocate the output arrays
fvtempz = nan(nf, fz);
fvsalz = nan(nf, fz);
for pp = 1:nf
% Get the FVCOM depths at this node
tfz = Mobj.siglayz(oNodes(pp), :);
% Now get the interpolated AMM depth at this node
tpz = idepthz(pp, :);
% Get the temperature and salinity values for this node and
% interpolate down the water column (from AMM to FVCOM). I had
% originally planned to use csaps for the vertical
% interplation/subsampling at each location. However, the demo
% of csaps in the MATLAB documentation makes the interpolation
% look horrible (shaving off extremes). interp1 provides a
% range of interpolation schemes, of which pchip seems to do a
% decent job of the interpolation (at least qualitatively).
if ~isnan(tpz)
fvtempz(pp, :) = interp1(tpz, itempz(pp, :), tfz, 'linear', 'extrap');
fvsalz(pp, :) = interp1(tpz, isalz(pp, :), tfz, 'linear', 'extrap');
else
warning('Should never see this... ') % because we test for NaNs when fetching the values.
warning('FVCOM boundary node at %f, %f is outside the AMM domain. Skipping.', fvlon(pp), fvlat(pp))
continue
end
end
end
% The horizontally- and vertically-interpolated values in the final
% FVCOM results array.
fvtemp(:, :, t) = fvtempz;
fvsal(:, :, t) = fvsalz;
end
if ftbverbose
toc
end
Mobj.temperature = fvtemp;
Mobj.salt = fvsal;
% Convert the current times to Modified Julian Day (this is a bit ugly).
amm.time.all = strtrim(regexp(amm.time.units, 'since', 'split'));
amm.time.datetime = strtrim(regexp(amm.time.all{end}, ' ', 'split'));
amm.time.ymd = str2double(strtrim(regexp(amm.time.datetime{1}, '-', 'split')));
amm.time.hms = str2double(strtrim(regexp(amm.time.datetime{2}, ':', 'split')));
Mobj.ts_times = greg2mjulian(...
amm.time.ymd(1), ...
amm.time.ymd(2), ...
amm.time.ymd(3), ...
amm.time.hms(1), ...
amm.time.hms(2), ...
amm.time.hms(3)) + (amm.time.data / 3600 / 24);
if ftbverbose
fprintf(['end : ' subname '\n'])
end
fvcom_prepro/get_AMM_netCDF.m → fvcom_prepro/get_POLCOMS_netCDF.m
View file @ d59b9c5a
function ncdata = get_AMM_netCDF(files, varlist)
function ncdata = get_POLCOMS_netCDF(files, varlist)
% Read temperature and salinity from NetCDF model output files and
% interpolate onto the open boundaries in Mobj.
%
% function struct = get_AMM_netCDF(Mobj, files, varlist)
% function struct = get_POLCOMS_netCDF(Mobj, files, varlist)
%
% DESCRIPTION:
% Extract variables in varlist to a struct from the files given in
......@@ -19,7 +19,7 @@ function ncdata = get_AMM_netCDF(files, varlist)
% blank for that variable.
%
% EXAMPLE USAGE
% S = get_AMM_netCDF({'/tmp/2000.nc', '/tmp/2001.nc', {'temp', 'salt'})
% S = get_POLCOMS_netCDF({'/tmp/2000.nc', '/tmp/2001.nc', {'temp', 'salt'})
%
% NOTES:
%
......@@ -43,7 +43,7 @@ function ncdata = get_AMM_netCDF(files, varlist)
%
%==========================================================================
subname = 'get_AMM_netCDF';
subname = 'get_POLCOMS_netCDF';
global ftbverbose;
if ftbverbose
......@@ -112,7 +112,7 @@ for ii = 1:todo
% a time dimension.
ncdata.(getVar).data = cat(ndims(data), ncdata.(getVar).data, data);
else
error('Unsupported number of dimensions in AMM data')
error('Unsupported number of dimensions in PML POLCOMS-ERSEM data')
end
end
% Try to get some units (important for the calculation of MJD).
......
fvcom_prepro/get_POLCOMS_tsobc.m
View file @ d59b9c5a
This diff is collapsed.
fvcom_prepro/interp_AMM2FVCOM.m → fvcom_prepro/interp_POLCOMS2FVCOM.m
View file @ d59b9c5a
function Mobj = interp_AMM2FVCOM(Mobj, ts, start_date, varlist)
function Mobj = interp_POLCOMS2FVCOM(Mobj, ts, start_date, varlist)
% Use an FVCOM restart file to seed a model run with spatially varying
% versions of otherwise constant variables (temperature and salinity only
% for the time being).
%
% function interp_AMM2FVCOM(Mobj, ts, start_date, fv_restart, varlist)
% function interp_POLCOMS2FVCOM(Mobj, ts, start_date, fv_restart, varlist)
%
% DESCRIPTION:
% FVCOM does not yet support spatially varying temperature and salinity
......@@ -32,7 +32,7 @@ function Mobj = interp_AMM2FVCOM(Mobj, ts, start_date, varlist)
% temperature).
%
% EXAMPLE USAGE
% interp_AMM2FVCOM(Mobj, '/tmp/ts.nc', '2006-01-01 00:00:00', ...
% interp_POLCOMS2FVCOM(Mobj, '/tmp/ts.nc', '2006-01-01 00:00:00', ...
% {'lon', 'lat', 'ETWD', 'x1XD', 'time'})
%
% Author(s):
......@@ -43,7 +43,7 @@ function Mobj = interp_AMM2FVCOM(Mobj, ts, start_date, varlist)
%
%==========================================================================
subname = 'interp_AMM2FVCOM';
subname = 'interp_POLCOMS2FVCOM';
global ftbverbose;
if ftbverbose
......@@ -57,29 +57,29 @@ end
% Data format:
%
% amm.ETWD.data and amm.x1XD.data are y, x, sigma, time
% pc.ETWD.data and pc.x1XD.data are y, x, sigma, time
%
amm = get_AMM_netCDF(ts, varlist);
pc = get_POLCOMS_netCDF(ts, varlist);
% Number of sigma layers.
[fn, fz] = size(Mobj.siglayz);
% Make rectangular arrays for the nearest point lookup.
[lon, lat] = meshgrid(amm.lon.data, amm.lat.data);
[lon, lat] = meshgrid(pc.lon.data, pc.lat.data);
% Convert the current times to Modified Julian Day (this is a bit ugly).
amm.time.all = strtrim(regexp(amm.time.units, 'since', 'split'));
amm.time.datetime = strtrim(regexp(amm.time.all{end}, ' ', 'split'));
amm.time.ymd = str2double(strtrim(regexp(amm.time.datetime{1}, '-', 'split')));
amm.time.hms = str2double(strtrim(regexp(amm.time.datetime{2}, ':', 'split')));
pc.time.all = strtrim(regexp(pc.time.units, 'since', 'split'));
pc.time.datetime = strtrim(regexp(pc.time.all{end}, ' ', 'split'));
pc.time.ymd = str2double(strtrim(regexp(pc.time.datetime{1}, '-', 'split')));
pc.time.hms = str2double(strtrim(regexp(pc.time.datetime{2}, ':', 'split')));
Mobj.ts_times = greg2mjulian(...
amm.time.ymd(1), ...
amm.time.ymd(2), ...
amm.time.ymd(3), ...
amm.time.hms(1), ...
amm.time.hms(2), ...
amm.time.hms(3)) + (amm.time.data / 3600 / 24);
pc.time.ymd(1), ...
pc.time.ymd(2), ...
pc.time.ymd(3), ...
pc.time.hms(1), ...
pc.time.hms(2), ...
pc.time.hms(3)) + (pc.time.data / 3600 / 24);
% Given our intput time (in start_date), find the nearest time
% index for the regularly gridded data.
......@@ -98,13 +98,13 @@ if ftbverbose
end
% Permute the arrays to be x by y rather than y by x.
temperature = permute(squeeze(amm.ETWD.data(:, :, :, tidx)), [2, 1, 3]);
salinity = permute(squeeze(amm.x1XD.data(:, :, :, tidx)), [2, 1, 3]);
depth = permute(squeeze(amm.depth.data(:, :, :, tidx)), [2, 1, 3]);
temperature = permute(squeeze(pc.ETWD.data(:, :, :, tidx)), [2, 1, 3]);
salinity = permute(squeeze(pc.x1XD.data(:, :, :, tidx)), [2, 1, 3]);
depth = permute(squeeze(pc.depth.data(:, :, :, tidx)), [2, 1, 3]);
mask = depth(:, :, end) >= 0; % land is positive.
amm.tempz = grid_vert_interp(Mobj, lon, lat, temperature, depth, mask);
amm.salz = grid_vert_interp(Mobj, lon, lat, salinity, depth, mask);
pc.tempz = grid_vert_interp(Mobj, lon, lat, temperature, depth, mask);
pc.salz = grid_vert_interp(Mobj, lon, lat, salinity, depth, mask);
if ftbverbose
fprintf('done.\n')
......@@ -125,8 +125,8 @@ fvtemp = nan(fn, fz);
fvsalt = nan(fn, fz);
for zi = 1:fz
% Set up the interpolation object.
ft = TriScatteredInterp(lon(:), lat(:), reshape(amm.tempz(:, :, zi), [], 1), 'natural');
fs = TriScatteredInterp(lon(:), lat(:), reshape(amm.salz(:, :, zi), [], 1), 'natural');
ft = TriScatteredInterp(lon(:), lat(:), reshape(pc.tempz(:, :, zi), [], 1), 'natural');
fs = TriScatteredInterp(lon(:), lat(:), reshape(pc.salz(:, :, zi), [], 1), 'natural');
% Interpolate temperature and salinity onto the unstructured grid.
fvtemp(:, zi) = ft(Mobj.lon, Mobj.lat);
fvsalt(:, zi) = fs(Mobj.lon, Mobj.lat);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment