function Mobj = get_POLCOMS_meanflow(Mobj, files) % Read mean flow from the PML POLCOMS-ERSEM NetCDF AMM model output files % and interpolate onto the open boundaries in Mobj. % % function Mobj = get_POLCOMS_meanflow(Mobj, ts, polcoms_bathy, varlist) % % DESCRIPTION: % Interpolate u and v flow vectors to calculate daily mean flow at the % FVCOM open boundaries at all sigma levels. % % INPUT: % Mobj = MATLAB mesh structure which must contain: % - Mobj.lon, Mobj.lat - node coordinates (lat/long). % - Mobj.obc_elements - list of open boundary element % inidices. % - Mobj.nObcElements - number of elements in each open boundary. % files = Cell array of PML POLCOMS-ERSEM NetCDF file(s) in which 4D % variables of u and v velocity components (called 'ucurD' and % 'vcurD') exist. Their shape should be (y, x, sigma, time). % % NOTES: % % - If you supply multiple files in files, 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 six new fields (mf_time, meanflow_u, % meanflow_v, meanflow_ubar, meanflow_ubar and velocity). % meanflow_ubar, meanflow_vbar and velocity have sizes of % (sum(Mobj.nObcNodes), time); meanflow_u and meanflow_v have % sizes (sum(Mobj.nObcNodes), siglay, time). The time dimension % is determined based on the input NetCDF file. The mf_times % variable is just the input file times in Modified Julian Day. % % EXAMPLE USAGE % Mobj = get_POLCOMS_meanflow(Mobj, files) % % Author(s): % Pierre Cazenave (Plymouth Marine Laboratory) % % Revision history % 2013-02-20 First version. % 2013-02-26 Add interpolation of the u and v vectors separately and % then calculate the interpolated velocity at the end. % 2013-02-27 Change the vertical interpolation to be scaled within the % POLCOMS-ERSEM depth range for the current node. The net result is that % the vertical profiles are squashed or stretched to fit within the % FVCOM depths. This means the full profile structure is maintained in % the resulting FVCOM boundary input despite the differing depths at the % FVCOM boundary node. % 2013-02-28 Change the interpolation to occur on the open boundary % elements rather than on the open boundary nodes. This requires a % couple of extra steps in the pre-processing (notably running % find_boundary_elements) as well as transferring the sigma depths to % the element centres. % %========================================================================== subname = 'get_POLCOMS_meanflow'; global ftbverbose; if ftbverbose fprintf('\n') fprintf(['begin : ' subname '\n']) end varlist = {'lon', 'lat', 'ucurD', 'vcurD', 'time', 'depth', 'pdepthD'}; pc = get_POLCOMS_netCDF(files, varlist); [~, ~, nz, nt] = size(pc.ucurD.data); [lon, lat] = meshgrid(pc.lon.data, pc.lat.data); % We need positions of the open boundary element centroids (i.e. the MCE % centroid) at which to calculate the mean flow velocity. As such, we need % to first calculate a number of parameters at the element centres from the % nodal values we typically collect (e.g. positions, sigma layers etc.). Do % all that here. lonc = nodes2elems(Mobj.lon, Mobj); latc = nodes2elems(Mobj.lat, Mobj); % For the sigma levels, we need to wrap the conversion in a loop. siglayzc = nan([Mobj.nElems, size(Mobj.siglayz, 2)]); for zz = 1:size(Mobj.siglayz, 2) siglayzc(:, zz) = nodes2elems(Mobj.siglayz(:, zz), Mobj); end oElements = [Mobj.read_obc_elements{:}]; % Find the FVCOM positions for all the open boundary elements. fvlon = lonc(oElements); fvlat = latc(oElements); % Number of open boundary elements ne = sum(Mobj.nObcElements); % Number of sigma layers. fz = size(siglayzc, 2); fvmfu = nan(ne, fz, nt); % FVCOM interpolated flow vector components fvmfv = nan(ne, fz, nt); % FVCOM interpolated flow vector components if ftbverbose tic end for t = 1:nt if ftbverbose fprintf('%s : %i of %i timesteps... ', subname, t, nt) end % Get the current 3D array of PML POLCOMS-ERSEM results. pcu3 = pc.ucurD.data(:, :, :, t); pcv3 = pc.vcurD.data(:, :, :, t); pcz3 = pc.depth.data(:, :, :, t); % Preallocate the intermediate results arrays. iuz = nan(ne, nz); ivz = nan(ne, nz); izz = nan(ne, nz); % We need to create a mask to eliminate land values and apply it to the % depth averaged values. mask = squeeze(pc.depth.data(:, :, end, t))' >= 0; % Mask the longitude and latitude here (since it's depth independent, % there's no point doing it inside the next loop). tlon = lon; tlat = lat; tlon(mask) = []; tlat(mask) = []; for z = 1:nz % Turns out we do need vertical velocity for the mean flow, so % iterate through each vertical layer before interpolating % horizontally. The vertical interpolation happens after the % horizontal has been completed. Transpose the arrays to be (x, y) % rather than (y, x). pcu2 = pcu3(:, :, z)'; pcv2 = pcv3(:, :, z)'; pcz2 = pcz3(:, :, z)'; % Flatten the arrays through the masking. pcu2(mask) = []; pcv2(mask) = []; pcz2(mask) = []; % Preallocate the intermediate results arrays. iuobc = nan(ne, 1); ivobc = nan(ne, 1); izobc = nan(ne, 1); % Speed up the tightest loops with a parallelized loop. parfor i = 1:ne % Now we can do each position within the current depth layer. fx = fvlon(i); fy = fvlat(i); [~, ii] = sort(sqrt((tlon - fx).^2 + (tlat - fy).^2)); % Get the n nearest elements 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); pu = pcu2(ixy); pv = pcv2(ixy); pz = pcz2(ixy); % Use a triangulation to do the horizontal interpolation. triu = TriScatteredInterp(plon', plat', pu', 'natural'); triv = TriScatteredInterp(plon', plat', pv', 'natural'); triz = TriScatteredInterp(plon', plat', pz', 'natural'); iuobc(i) = triu(fx, fy); ivobc(i) = triv(fx, fy); izobc(i) = triz(fx, fy); % Check if we have NaNs (mostly if the position is outside the % model domain). if isnan(iuobc(i)) || isnan(ivobc(i)) warning('FVCOM boundary element at %f, %f is outside the PML POLCOMS-ERSEM domain. Setting to the closest PML POLCOMS-ERSEM value.', fx, fy) iuobc(i) = pcu2(ii(1)); ivobc(i) = pcv2(ii(1)); izobc(i) = pcz2(ii(1)); end end % Put the results in this intermediate array. iuz(:, z) = iuobc; ivz(:, z) = ivobc; izz(:, z) = izobc; end % Now we've interpolated in space, we can interpolate the z-values % to the sigma depths. % Preallocate the output arrays fvuz = nan(ne, fz); fvvz = nan(ne, fz); for pp = 1:ne % Get the FVCOM depths at this element tfz = siglayzc(oElements(pp), :); % Now get the interpolated PML POLCOMS-ERSEM depth at this element tpz = izz(pp, :); % To ensure we get the full vertical expression of the vertical % profiles, we need to normalise the POLCOMS-ERSEM depths to the % FVCOM range for the current element. This is because in instances % where FVCOM depths are shallower (e.g. in coastal regions), if we % don't normalise the depths, we end up truncating the vertical % profile. This approach ensures we always use the full vertical % profile, but we're potentially squeezing it into a smaller depth. A = max(tpz); B = min(tpz); C = max(tfz); D = min(tfz); norm_tpz = (((D - C) * (tpz - A)) / (B - A)) + C; % Get the u and v velocity values for this elment and 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) fvuz(pp, :) = interp1(norm_tpz, iuz(pp, :), tfz, 'linear', 'extrap'); fvvz(pp, :) = interp1(norm_tpz, ivz(pp, :), tfz, 'linear', 'extrap'); else warning('Should never see this... ') % because we test for NaNs when fetching the values. warning('FVCOM boundary element at %f, %f is outside the PML POLCOMS-ERSEM domain. Skipping.', fvlon(pp), fvlat(pp)) continue end end % The horizontally- and vertically-interpolated values in the final % FVCOM results array. fvmfu(:, :, t) = fvuz; fvmfv(:, :, t) = fvvz; if ftbverbose fprintf('done.\n') end end if ftbverbose toc end % Stick the values in the mesh structure. Mobj.meanflow_u = fvmfu; Mobj.meanflow_v = fvmfv; % Now we have the 3D arrays, create depth averaged velocities too Mobj.meanflow_ubar = squeeze(mean(Mobj.meanflow_u, 2)); Mobj.meanflow_vbar = squeeze(mean(Mobj.meanflow_v, 2)); % Depth averaged velocity Mobj.velocity = mean(sqrt(Mobj.meanflow_u(:, :, end).^2 + Mobj.meanflow_v(:, :, end).^2), 2); % 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'))); Mobj.mf_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 % Check the interpolation along the boundary at the last time step. N.B. % Since the depths are scaled from the POLCOMS range to the FVCOM range, % the profiles won't match perfectly in these plots. This is because the % interpolated FVCOM profiles use the full water column from the POLCOMS % data rather than truncating it at the FVCOM depth. % % Map the open boundaries with the depth averaged velocity as a background % figure(1) % clf % imagesc(pc.lon.data, pc.lat.data, mean(sqrt(pcu3.^2 + pcv3.^2), 3)') % shading flat % axis('equal', 'tight') % set(gca,'YDir','normal') % colorbar % caxis([0, 0.05]) % hold on % plot(fvlon, fvlat, 'wo') % axis([min(fvlon) - 0.1, max(fvlon) + 0.1, min(fvlat) - 0.1, max(fvlat) + 0.1]) % % for i = 1:ne % % % Add the current element's position and value as a scatter point to % % the map plot. % scatter(fvlon(i), fvlat(i), 50, Mobj.velocity(i), 'filled') % % % Do vertical profiles of u, v and velocity. % figure(2) % clf % % subplot(3,1,1) % % FVCOM vertical profile. Has to be (i, :, end) because the % % corresponding POLCOMS data isn't stored as a function of time (i.e. % % iuz, ivz and izz are all for the last time step only). % fvuz = Mobj.meanflow_u(i, :, end); % fvz = siglayzc(oElements(i), :); % plot(fvuz, fvz) % hold on % % The interpolated POLCOMS vertical profile (last time step only) % plot(iuz(i, :), izz(i, :), 'g') % % The depth-averaged velocity (again, last time step only) % fvubar = Mobj.meanflow_ubar(i, end); % plot([fvubar, fvubar], [min(fvz), max(fvz)], 'k') % xlim([-0.1, 0.2]) % ylim([-100, 0]) % xlabel('Mean flow u-velocity (m^{3}s^{-1})') % ylabel('Depth (m)') % title('u-component') % legend('FVCOM', 'POLCOMS', 'Mean FVCOM', 'Location', 'SouthEast') % legend('boxoff') % % subplot(3,1,2) % % FVCOM vertical profile. Has to be (i, :, end) because the % % corresponding POLCOMS data isn't stored as a function of time (i.e. % % iuz, ivz and izz are all for the last time step only). % fvvz = Mobj.meanflow_v(i, :, end); % fvz = siglayzc(oElements(i), :); % plot(fvvz, fvz) % hold on % % The interpolated POLCOMS vertical profile (last time step only) % plot(ivz(i, :), izz(i, :), 'g') % % The depth-averaged velocity (again, last time step only) % fvvbar = Mobj.meanflow_vbar(i, end); % plot([fvvbar, fvvbar], [min(fvz), max(fvz)], 'k') % xlim([-0.1, 0.2]) % ylim([-100, 0]) % xlabel('Mean flow v-velocity (m^{3}s^{-1})') % ylabel('Depth (m)') % title('v-component') % legend('FVCOM', 'POLCOMS', 'Mean FVCOM', 'Location', 'SouthEast') % legend('boxoff') % % subplot(3,1,3) % % FVCOM vertical profile. Has to be (i, :, end) because the % % corresponding POLCOMS data isn't stored as a function of time (i.e. % % iuz, ivz and izz are all for the last time step only). % fvvelz = sqrt(Mobj.meanflow_u(i, :, end).^2 + Mobj.meanflow_v(i, :, end).^2); % fvz = siglayzc(oElements(i), :); % plot(fvvelz, fvz) % hold on % % The interpolated POLCOMS vertical profile (last time step only) % plot(sqrt(iuz(i, :).^2 + ivz(i, :).^2), izz(i, :), 'g') % % The depth-averaged velocity (again, last time step only) % fvvelbar = Mobj.velocity(i); % plot([fvvelbar, fvvelbar], [min(fvz), max(fvz)], 'k') % xlim([-0.1, 0.2]) % ylim([-100, 0]) % xlabel('Mean flow velocity (m^{3}s^{-1})') % ylabel('Depth (m)') % title('velocity') % legend('FVCOM', 'POLCOMS', 'Mean FVCOM', 'Location', 'SouthEast') % legend('boxoff') % pause(0.1) % end