diff --git a/fvcom_prepro/get_AMM_tsobc.m b/fvcom_prepro/get_AMM_tsobc.m deleted file mode 100644 index 6d39544c148c6e3c5e4e3d1ac9f0c4c6c738d497..0000000000000000000000000000000000000000 --- a/fvcom_prepro/get_AMM_tsobc.m +++ /dev/null @@ -1,236 +0,0 @@ -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 diff --git a/fvcom_prepro/get_AMM_netCDF.m b/fvcom_prepro/get_POLCOMS_netCDF.m similarity index 93% rename from fvcom_prepro/get_AMM_netCDF.m rename to fvcom_prepro/get_POLCOMS_netCDF.m index 2d56e71c79c6571a3bb99b530cce4c9ae42ff076..53a13876d43e57f32071a48d5409723f7b623b96 100644 --- a/fvcom_prepro/get_AMM_netCDF.m +++ b/fvcom_prepro/get_POLCOMS_netCDF.m @@ -1,8 +1,8 @@ -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). diff --git a/fvcom_prepro/get_POLCOMS_tsobc.m b/fvcom_prepro/get_POLCOMS_tsobc.m index 4e996b38e104798e98a492c06251a19b2919150e..6c6298e56268969b319a1ad67b4bc39ede4bbe55 100644 --- a/fvcom_prepro/get_POLCOMS_tsobc.m +++ b/fvcom_prepro/get_POLCOMS_tsobc.m @@ -1,40 +1,35 @@ -function Mobj = get_POLCOMS_tsobc(Mobj, polcoms_ts, polcoms_z) -% Read temperature and salinity from POLCOMS NetCDF model output files and -% interpolate onto the open boundaries in Mobj. +function Mobj = get_POLCOMS_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_POLCOMS_tsobc(Mobj, polcoms_ts, polcoms_bathy, varlist) +% function Mobj = get_POLCOMS_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. -% polcoms_ts = Cell array of POLCOMS NetCDF file(s) in which 4D -% variables of temperature and salinity (called 'ETW' and -% 'x1X') exist. Their shape should be (y, x, sigma, time). -% polcoms_z = Cell array of POLCOMS NetCDF file(s) in which 4D -% variables of bathymetry and sigma layer thickness can be -% found. They should be called 'depth' and 'pdepth' -% respectively. +% 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 PML POLCOMS-ERSEM 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: -% - The POLCOMS NetCDF files should also contain 'time', 'lat' and 'lon' -% variables. % -% - If you supply multiple files in polcoms_ts, there are a few -% assumptions: +% - 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, @@ -44,15 +39,13 @@ function Mobj = get_POLCOMS_tsobc(Mobj, polcoms_ts, polcoms_z) % is just the input file times in Modified Julian Day. % % EXAMPLE USAGE -% Mobj = get_POLCOMS_forcing(Mobj, polcoms_ts, polcoms_z) +% Mobj = get_POLCOMS_tsobc(Mobj, ts, depth) % % Author(s): % Pierre Cazenave (Plymouth Marine Laboratory) % % Revision history -% 2013-01-09 First version based on the FVCOM shelf model -% get_POLCOMS_forcing.m script (i.e. not a function but a plain script). -% 2013-02-04 Add support for reading in multiple NetCDF files. +% 2013-02-07 First version based on get_POLCOMS_tsobc.m. % %========================================================================== @@ -64,147 +57,19 @@ if ftbverbose fprintf(['begin : ' subname '\n']) end -varlist = {'lon', 'lat', 'ETW', 'x1X', 'time'}; - -% Get the results. Check we have a cell array, and if we don't, assume it's -% a file name. - -if iscell(polcoms_ts) - todo = length(polcoms_ts); - todo2 = length(polcoms_z); - if todo2 ~= todo - error('Supply matching POLCOMS temperature/salinity and depth NetCDF files') - end - clear todo2 -else - todo = 1; -end - -for ii = 1:todo - - if ftbverbose - if iscell(polcoms_ts) - ftn = polcoms_ts{ii}; - fzn = polcoms_z{ii}; - else - ftn = polcoms_ts; - fzn = polcoms_z; - end - % Strip path from filename for the verbose output. - [~, basename, ext] = fileparts(ftn); - tmp_fn = [basename, ext]; - - if todo == 1 - fprintf('%s: extracting file %s... ', subname, tmp_fn) - else - fprintf('%s: extracting file %s (%i of %i)... ', subname, tmp_fn, ii, todo) - end - end - - nc = netcdf.open(ftn, 'NOWRITE'); - - for var = 1:numel(varlist) - - getVar = varlist{var}; - varid_pc = netcdf.inqVarID(nc, getVar); - - data = double(netcdf.getVar(nc, varid_pc, 'single')); - if ii == 1 - pc.(getVar).data = data; - else - if ndims(data) < 4 - if strcmpi(varlist{var}, 'time') - % If the dimension is time, we need to be a bit more - % clever since we'll need a concatenated time series - % (in which values are continuous and from which we - % can extract a sensible time). As such, we need to add - % the maximum of the existing time. On the first - % iteration, we should save ourselves the base time - % (from the units of the time variable). - pc.(getVar).data = [pc.(getVar).data; data + max(pc.(getVar).data)]; - else - % This should be a fixed set of values (probably lon or - % lat) in which case we don't need to append them, so - % just replace the existing values with those in the - % current NetCDF file. - pc.(getVar).data = data; - end - elseif ndims(data) == 4 - % Assume we're concatenating with time (so along the fourth - % dimesnion. - pc.(getVar).data = cat(4, pc.(getVar).data, data); - else - error('Unsupported number of dimensions in POLCOMS data') - end - end - % Try to get some units (important for the calculation of MJD). - try - if ii == 1 - units = netcdf.getAtt(nc, varid_pc, 'units'); - else - % Leave the units values alone so we always use the values - % from the first file. This is particularly important for - % the time calculation later on which is dependent on - % knowing the time origin of the first file. - continue - end - catch - units = []; - end - pc.(getVar).units = units; - end - - netcdf.close(nc) - - % Extract the bathymetry ('pdepth' is cell thickness, 'depth' is cell - % centre depth). We still need to append along the fourth dimension - % here too (I think depth and pdepth include the tidal component). - nc = netcdf.open(fzn, 'NOWRITE'); - varid_pc = netcdf.inqVarID(nc, 'depth'); - data = double(netcdf.getVar(nc, varid_pc, 'single')); - if ii == 1 - pc.depth.data = data; - else - pc.depth.data = cat(4, pc.depth.data, data); - end - try - pc.depth.units = netcdf.getAtt(nc, varid_pc, 'units'); - catch - pc.depth.units = []; - end - clear data - - varid_pc = netcdf.inqVarID(nc, 'pdepth'); - data = double(netcdf.getVar(nc, varid_pc, 'single')); - if ii == 1 - pc.pdepth.data = data; - else - pc.pdepth.data = cat(4, pc.pdepth.data, data); - end - try - pc.pdepth.units = netcdf.getAtt(nc, varid_pc, 'units'); - catch - pc.pdepth.units = []; - end - clear data - netcdf.close(nc) - - if ftbverbose - fprintf('done.\n') - end - -end +varlist = {'lon', 'lat', 'ETWD', 'x1XD', 'time', 'depth', 'pdepthD'}; % Data format: % -% pc.ETW.data and pc.x1X.data are y, x, sigma, time +% amm.ETWD.data and amm.x1XD.data are y, x, sigma, time % -[~, ~, nz, nt] = size(pc.ETW.data); +pc = get_POLCOMS_netCDF(ts, varlist); + +[~, ~, nz, nt] = size(pc.ETWD.data); % Make rectangular arrays for the nearest point lookup. [lon, lat] = meshgrid(pc.lon.data, pc.lat.data); -% Find the nearest POLCOMS point to each point in the FVCOM open boundaries fvlon = Mobj.lon(Mobj.obc_nodes(Mobj.obc_nodes ~= 0)); fvlat = Mobj.lat(Mobj.obc_nodes(Mobj.obc_nodes ~= 0)); @@ -220,9 +85,9 @@ if ftbverbose tic end for t = 1:nt - % Get the current 3D array of POLCOMS results. - pctemp3 = pc.ETW.data(:, :, :, t); - pcsal3 = pc.x1X.data(:, :, :, t); + % Get the current 3D array of PML POLCOMS-ERSEM results. + ammtemp3 = pc.ETWD.data(:, :, :, t); + ammsalt3 = pc.x1XD.data(:, :, :, t); % Preallocate the intermediate results arrays. itempz = nan(nf, nz); @@ -232,14 +97,14 @@ for t = 1:nt for j = 1:nz % Now extract the relevant layer from the 3D subsets. Transpose the % data to be (x, y) rather than (y, x). - pctemp2 = pctemp3(:, :, j)'; - pcsal2 = pcsal3(:, :, j)'; - pcdepth2 = squeeze(pc.depth.data(:, :, j, t))'; + ammtemp2 = ammtemp3(:, :, j)'; + ammsalt2 = ammsalt3(:, :, j)'; + ammdepth2 = squeeze(pc.depth.data(:, :, j, t))'; % Create new arrays which will be flattened when masking (below). - tpctemp2 = pctemp2; - tpcsal2 = pcsal2; - tpcdepth2 = pcdepth2; + tammtemp2 = ammtemp2; + tammsalt2 = ammsalt2; + tammdepth2 = ammdepth2; tlon = lon; tlat = lat; @@ -247,11 +112,11 @@ for t = 1:nt % 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 POLCOMS data is irregularly spaced). - mask = tpcdepth2 < -20000; - tpctemp2(mask) = []; - tpcsal2(mask) = []; - tpcdepth2(mask) = []; + % the PML POLCOMS-ERSEM 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 @@ -272,16 +137,17 @@ for t = 1:nt fy = fvlat(i); [~, ii] = sort(sqrt((tlon - fx).^2 + (tlat - fy).^2)); - % Get the n nearest nodes from POLCOMS (more? fewer?). + % Get the n nearest nodes from PML POLCOMS-ERSEM data (more? + % fewer?). ixy = ii(1:16); % Get the variables into static variables for the % parallelisation. plon = tlon(ixy); plat = tlat(ixy); - ptemp = tpctemp2(ixy); - psal = tpcsal2(ixy); - pdepth = tpcdepth2(ixy); + ptemp = tammtemp2(ixy); + psal = tammsalt2(ixy); + pdepth = tammdepth2(ixy); % Use a triangulation to do the horizontal interpolation. tritemp = TriScatteredInterp(plon', plat', ptemp', 'natural'); @@ -293,10 +159,10 @@ for t = 1:nt % 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 POLCOMS domain. Setting to the closest POLCOMS value.', fx, fy) - itempobc(i) = tpctemp2(ii(1)); - isalobc(i) = tpcsal2(ii(1)); - idepthobc(i) = tpcdepth2(ii(1)); + warning('FVCOM boundary node at %f, %f is outside the PML POLCOMS-ERSEM domain. Setting to the closest PML POLCOMS-ERSEM value.', fx, fy) + itempobc(i) = tammtemp2(ii(1)); + isalobc(i) = tammsalt2(ii(1)); + idepthobc(i) = tammdepth2(ii(1)); end end @@ -318,24 +184,24 @@ for t = 1:nt for pp = 1:nf % Get the FVCOM depths at this node tfz = Mobj.siglayz(oNodes(pp), :); - % Now get the interpolated POLCOMS depth at this node. + % Now get the interpolated PML POLCOMS-ERSEM depth at this node tpz = idepthz(pp, :); % Get the temperature and salinity values for this node and - % interpolate down the water column (from POLCOMS to FVCOM). - % TODO: Use csaps for the vertical interplation/subsampling at - % each location. Alternatively, the pchip interp1 method seems - % to do a decent job of the interpolation; it might be a more - % suitable candidate in the absence of csaps. In fact, the demo - % of csaps in the MATLAB documentation makes the interpolation - % look horrible (shaving off extremes). I think pchip is - % better. + % interpolate down the water column (from PML POLCOMS-ERSEM 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 POLCOMS domain. Skipping.', fvlon(pp), fvlat(pp)) + warning('FVCOM boundary node at %f, %f is outside the PML POLCOMS-ERSEM domain. Skipping.', fvlon(pp), fvlat(pp)) continue end end @@ -353,27 +219,20 @@ end Mobj.temperature = fvtemp; Mobj.salt = fvsal; -% Do we have to interpolate to the FVCOM time series? Looking at page 325 -% of the FVCOM manual, it looks like the temperature and salinity are on a -% different time sampling from the other example files (14 time steps vs. -% 3625 for _wnd.nc or 43922 for _julain_obc.nc (i.e. surface elevation at -% the boundary)). That's not to say those files were all used in the same -% model run... In the interim, just convert the current times to Modified -% Julian Day (this is a bit ugly). -pc.time.yyyymmdd = strtrim(regexp(pc.time.units, 'since', 'split')); -pc.time.yyyymmdd = str2double(regexp(pc.time.yyyymmdd{end}, '-', 'split')); -Mobj.ts_times = greg2mjulian(pc.time.yyyymmdd(1), pc.time.yyyymmdd(2), pc.time.yyyymmdd(3), 0, 0, 0) + (pc.time.data / 3600 / 24); +% Convert the current times to Modified Julian Day (this is a bit ugly). +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'))); -% % Convert the POLCOMS times to Modified Julian Day (this is a very ugly). -% pc.time.yyyymmdd = strtrim(regexp(pc.time.units, 'since', 'split')); -% pc.time.strtime = regexp(pc.time.yyyymmdd{end}, '-', 'split'); -% % This new version of the time has the year in a weird format (yr.#). We -% % thus need to split it again to get the decimal year (post-2000 only?). -% pc.time.strtimeyr = regexp(pc.time.strtime, '\.', 'split'); -% pc.time.yyyymmdd = str2double([pc.time.strtimeyr{1}(2), pc.time.strtime{2:3}]); -% pc.time.yyyymmdd(1) = pc.time.yyyymmdd(1) + 2000; % add full year. -% Mobj.ts_times = greg2mjulian(pc.time.yyyymmdd(1), pc.time.yyyymmdd(2), pc.time.yyyymmdd(3), 0, 0, 0) + (pc.time.data / 3600 / 24); +Mobj.ts_times = greg2mjulian(... + 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); if ftbverbose fprintf(['end : ' subname '\n']) -end \ No newline at end of file +end diff --git a/fvcom_prepro/interp_AMM2FVCOM.m b/fvcom_prepro/interp_POLCOMS2FVCOM.m similarity index 79% rename from fvcom_prepro/interp_AMM2FVCOM.m rename to fvcom_prepro/interp_POLCOMS2FVCOM.m index f0c62bf4b0af257cb3c69cfb3ab7531da58988c9..0b15cdef4706aa3312397d356b5fcab50d26c3f9 100644 --- a/fvcom_prepro/interp_AMM2FVCOM.m +++ b/fvcom_prepro/interp_POLCOMS2FVCOM.m @@ -1,9 +1,9 @@ -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);