get_POLCOMS_tsobc.m 8.67 KB
Newer Older
1 2 3
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.
4
%
5
% function Mobj = get_POLCOMS_tsobc(Mobj, ts, polcoms_bathy, varlist)
6 7 8 9 10 11
%
% DESCRIPTION:
%    Interpolate temperature and salinity values onto the FVCOM open
%    boundaries at all sigma levels.
%
% INPUT:
12 13 14 15 16 17 18 19 20
%   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).
%
21 22
% NOTES:
%
23
%   - If you supply multiple files in ts, there are a few assumptions:
24 25 26 27 28 29
%
%       - 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.
30 31 32
% 
%   - The NetCDF files used here are those from the PML POLCOMS-ERSEM model
%   output.
33
%
34
% OUTPUT:
35 36
%    Mobj = MATLAB structure in which three new fields (called temperature,
%           salinity and ts_time). temperature and salinity have sizes
37
%           (sum(Mobj.nObcNodes), sigma, time). The time dimension is
38 39
%           determined based on the input NetCDF file. The ts_time variable
%           is just the input file times in Modified Julian Day.
40 41
%
% EXAMPLE USAGE
42
%    Mobj = get_POLCOMS_tsobc(Mobj, ts, depth)
43 44 45 46 47
%
% Author(s):
%    Pierre Cazenave (Plymouth Marine Laboratory)
%
% Revision history
48
%    2013-02-07 First version based on get_POLCOMS_tsobc.m.
49 50 51
%
%==========================================================================

Pierre Cazenave's avatar
Pierre Cazenave committed
52
subname = 'get_POLCOMS_tsobc';
53 54 55 56 57 58 59

global ftbverbose;
if ftbverbose
    fprintf('\n')
    fprintf(['begin : ' subname '\n'])
end

60
varlist = {'lon', 'lat', 'ETWD', 'x1XD', 'time', 'depth', 'pdepthD'};
61 62 63

% Data format:
% 
64
%   amm.ETWD.data and amm.x1XD.data are y, x, sigma, time
65
% 
66 67 68
pc = get_POLCOMS_netCDF(ts, varlist);

[~, ~, nz, nt] = size(pc.ETWD.data);
69 70

% Make rectangular arrays for the nearest point lookup.
71
[lon, lat] = meshgrid(pc.lon.data, pc.lat.data);
72 73 74 75 76 77

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);
78 79
% Number of sigma layers.
fz = size(Mobj.siglayz, 2);
80 81 82 83

fvtemp = nan(nf, fz, nt); % FVCOM interpolated temperatures
fvsal = nan(nf, fz, nt); % FVCOM interpolated salinities

84 85 86
if ftbverbose
    tic
end
87
for t = 1:nt
88 89 90
    % Get the current 3D array of PML POLCOMS-ERSEM results.
    ammtemp3 = pc.ETWD.data(:, :, :, t);
    ammsalt3 = pc.x1XD.data(:, :, :, t);
91 92 93 94 95 96 97 98 99
    
    % 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).
100 101 102
        ammtemp2 = ammtemp3(:, :, j)';
        ammsalt2 = ammsalt3(:, :, j)';
        ammdepth2 = squeeze(pc.depth.data(:, :, j, t))';
103 104
       
        % Create new arrays which will be flattened when masking (below).
105 106 107
        tammtemp2 = ammtemp2;
        tammsalt2 = ammsalt2;
        tammdepth2 = ammdepth2;
108 109 110 111 112 113 114
        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
115 116 117 118 119
        % the PML POLCOMS-ERSEM data is irregularly spaced).
        mask = tammdepth2 > 20000;
        tammtemp2(mask) = [];
        tammsalt2(mask) = [];
        tammdepth2(mask) = [];
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
        % 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));
140 141
            % Get the n nearest nodes from PML POLCOMS-ERSEM data (more?
            % fewer?).
142 143 144 145 146 147
            ixy = ii(1:16);

            % Get the variables into static variables for the
            % parallelisation.
            plon = tlon(ixy);
            plat = tlat(ixy);
148 149 150
            ptemp = tammtemp2(ixy);
            psal = tammsalt2(ixy);
            pdepth = tammdepth2(ixy);
151 152 153 154 155 156 157 158 159 160 161
            
            % 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))
162 163 164 165
                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));
166 167 168
            end
        end
        
169
        % Put the results in this intermediate array.
170 171 172 173 174 175 176 177 178
        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
179 180 181 182 183

        % Preallocate the output arrays
        fvtempz = nan(nf, fz);
        fvsalz = nan(nf, fz);

184 185
        for pp = 1:nf
            % Get the FVCOM depths at this node
186
            tfz = Mobj.siglayz(oNodes(pp), :);
187
            % Now get the interpolated PML POLCOMS-ERSEM depth at this node
188
            tpz = idepthz(pp, :);
189

190
            % Get the temperature and salinity values for this node and
191 192 193 194 195 196 197 198
            % 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).
199
            if ~isnan(tpz)
200 201
                fvtempz(pp, :) = interp1(tpz, itempz(pp, :), tfz, 'linear', 'extrap');
                fvsalz(pp, :) = interp1(tpz, isalz(pp, :), tfz, 'linear', 'extrap');
202
            else
203
                warning('Should never see this... ') % because we test for NaNs when fetching the values.
204
                warning('FVCOM boundary node at %f, %f is outside the PML POLCOMS-ERSEM domain. Skipping.', fvlon(pp), fvlat(pp))
205 206 207 208 209 210 211 212 213 214
                continue
            end
        end
    end
    
    % The horizontally- and vertically-interpolated values in the final
    % FVCOM results array.
    fvtemp(:, :, t) = fvtempz;
    fvsal(:, :, t) = fvsalz;
end
215 216 217
if ftbverbose
    toc
end
218 219 220 221

Mobj.temperature = fvtemp;
Mobj.salt = fvsal;

222 223 224 225 226
% 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')));
227

228 229 230 231 232 233 234
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);
235

236 237
if ftbverbose
    fprintf(['end   : ' subname '\n'])
238
end