Due to a shift in policy, from 0900 GMT on Wednesday 14th July 2021, we will be disabling ssh access to the server for external users. External users who wish to continue to access code repositories on the server will need to switch to using https. This can be accomplished in the following way: 1) On the repo on gitlab, use the clone dialogue and select ‘Clone with HTTPS’ to get the address of the repo; 2) From within the checkout of your repo run: $ git remote set-url origin HTTPS_ADDRESS. Here, replace HTTPS_ADDRESS with the address you have just copied from GitLab. Pulls and pushes will now require you to enter a username and password rather than using a ssh key. If you would prefer not to enter a password each time, you might consider caching your login credentials.

get_POLCOMS_meanflow.m 6.57 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
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_nodes - list of open boundary node inidices.
%               - Mobj.nObcNodes - number of nodes 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 two new fields (called meanflow,
%           and mf_times). meanflow has a size of (sum(Mobj.nObcNodes),
%           time). The time dimension is determined based on the input
%           NetCDF file. The mf_time 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.
48 49
%    2013-02-26 Add interpolation of the u and v vectors separately and
%    then calculate the interpolated velocity at the end.
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
%
%==========================================================================

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);

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);

75 76 77
% fvmf = nan(nf, nt); % FVCOM interpolated mean flow
fvmfu = nan(nf, nt); % FVCOM interpolated mean flow vectors
fvmfv = nan(nf, nt); % FVCOM interpolated mean flow vectors
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93

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);

    % Calculate the velocity and direction (?) from the u and v vectors.
    % This approach means we don't have to interpolate u and v separately
    % (which is expensive) but just do velocity on its own.
94
%     pcvel3 = sqrt(pcu3.^2 + pcv3.^2);
95 96 97 98 99 100 101
       
    % For the velocity, it appears we don't need to have a depth varying
    % value (instead mean flow is scaled according to the values in MFDIST,
    % although I'm not yet sure exactly what those values are...). So,
    % we'll calculate a depth average velocity here for the time being.
    % Transpose the array here so it's (x, y) rather than the original (y,
    % x).
102 103 104
%     pcvel3mean = mean(pcvel3, 3)';
    pcu3mean = mean(pcu3, 3)';
    pcv3mean = mean(pcv3, 3)';
105 106 107 108 109

    % 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;
    
110 111 112
%     tpcvel3mean = pcvel3mean;
    tpcu3mean = pcu3mean;
    tpcv3mean = pcv3mean;
113 114 115
    tlon = lon;
    tlat = lat;
    % Clear out the land values and flatten the arrays.
116 117 118
%     tpcvel3mean(mask) = [];
    tpcu3mean(mask) = [];
    tpcv3mean(mask) = [];
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
    tlon(mask) = [];
    tlat(mask) = [];
    
    % Speed up the tightest loop with a parallelized loop.
    parfor i = 1:nf
        % Now we can do each position within the depth averaged array.

        fx = fvlon(i);
        fy = fvlat(i);

        [~, ii] = sort(sqrt((tlon - fx).^2 + (tlat - fy).^2));
        % 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);
138 139 140
%         pvel = tpcvel3mean(ixy);
        pu = tpcu3mean(ixy);
        pv = tpcv3mean(ixy);
141 142

        % Use a triangulation to do the horizontal interpolation.
143 144 145 146 147 148
%         trivel = TriScatteredInterp(plon', plat', pvel', 'natural');
        triu = TriScatteredInterp(plon', plat', pu', 'natural');
        triv = TriScatteredInterp(plon', plat', pv', 'natural');
%         ivelobc(i) = trivel(fx, fy);
        iuobc(i) = triu(fx, fy);
        ivobc(i) = triv(fx, fy);
149 150 151

        % Check if we have NaNs (mostly if the position is outside the
        % model domain).
152
        if isnan(ivelobc(i)) || isnan(iuobc(i)) || isnan(ivobc(i))
153
            warning('FVCOM boundary node at %f, %f is outside the PML POLCOMS-ERSEM domain. Setting to the closest PML POLCOMS-ERSEM value.', fx, fy)
154 155 156
%             ivelobc(i) = tpcvel2(ii(1));
            iuobc(i) = tpcu2(ii(1));
            ivobc(i) = tpcv2(ii(1));
157 158 159 160 161
        end
    end

    % The horizontally-interpolated values in the final FVCOM results
    % array.
162 163 164
%     fvmf(:, t) = ivelobc;
    fvmfu(:, t) = iuobc;
    fvmfv(:, t) = ivobc;
165 166 167 168 169 170 171 172 173 174

    if ftbverbose
        fprintf('done.\n')
    end
end
if ftbverbose
    toc
end

% Stick the values in the mesh structure.
175 176 177 178 179
Mobj.meanflow_u = fvmfu;
Mobj.meanflow_v = fvmfv;
% Mobj.velocity = fvmf;
Mobj.velocity = sqrt(Mobj.meanflow_u.^2 + Mobj.meanflow_v.^2);

180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199

% 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