Read temperature and salinity from POLCOMS NetCDF model output files and
interpolate onto the open boundaries in Mobj.
function Mobj = get_POLCOMS_tsobc(Mobj, polcoms_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 = POLCOMS NetCDF file in which 4D variables of temperature
and salinity (called 'ETW' and 'x1X') exist. Their shape
should be (y, x, sigma, time).
polcoms_z = POLCOMS NetCDF file in which 4D variables of bathymetry
and sigma layer thickness can be found. They should be
called 'depth' and 'pdepth' respectively.
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_POLCOMS_forcing(Mobj, polcoms_ts, polcoms_z)
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).
==========================================================================0001 function Mobj = get_POLCOMS_tsobc(Mobj, polcoms_ts, polcoms_z) 0002 % Read temperature and salinity from POLCOMS NetCDF model output files and 0003 % interpolate onto the open boundaries in Mobj. 0004 % 0005 % function Mobj = get_POLCOMS_tsobc(Mobj, polcoms_ts, polcoms_bathy, varlist) 0006 % 0007 % DESCRIPTION: 0008 % Interpolate temperature and salinity values onto the FVCOM open 0009 % boundaries at all sigma levels. 0010 % 0011 % INPUT: 0012 % Mobj = MATLAB mesh structure which must contain: 0013 % - Mobj.siglayz - sigma layer depths for all model 0014 % nodes. 0015 % - Mobj.lon, Mobj.lat - node coordinates (lat/long). 0016 % - Mobj.obc_nodes - list of open boundary node inidices. 0017 % - Mobj.nObcNodes - number of nodes in each open 0018 % boundary. 0019 % polcoms_ts = POLCOMS NetCDF file in which 4D variables of temperature 0020 % and salinity (called 'ETW' and 'x1X') exist. Their shape 0021 % should be (y, x, sigma, time). 0022 % polcoms_z = POLCOMS NetCDF file in which 4D variables of bathymetry 0023 % and sigma layer thickness can be found. They should be 0024 % called 'depth' and 'pdepth' respectively. 0025 % 0026 % OUTPUT: 0027 % Mobj = MATLAB structure in which three new fields (called temperature, 0028 % salinity and ts_time). temperature and salinity have sizes 0029 % (sum(Mobj.nObcNodes), sigma, time). The time dimension is 0030 % determined based on the input NetCDF file. The ts_time variable 0031 % is just the input file times in Modified Julian Day. 0032 % 0033 % EXAMPLE USAGE 0034 % Mobj = get_POLCOMS_forcing(Mobj, polcoms_ts, polcoms_z) 0035 % 0036 % Author(s): 0037 % Pierre Cazenave (Plymouth Marine Laboratory) 0038 % 0039 % Revision history 0040 % 2013-01-09 First version based on the FVCOM shelf model 0041 % get_POLCOMS_forcing.m script (i.e. not a function but a plain script). 0042 % 0043 %========================================================================== 0044 0045 subname = 'get_POLCOMS_forcing'; 0046 0047 global ftbverbose; 0048 if ftbverbose 0049 fprintf('\n') 0050 fprintf(['begin : ' subname '\n']) 0051 end 0052 0053 varlist = {'lon', 'lat', 'ETW', 'x1X', 'time'}; 0054 0055 % Get the results 0056 nc = netcdf.open(polcoms_ts, 'NOWRITE'); 0057 0058 for var=1:numel(varlist) 0059 0060 getVar = varlist{var}; 0061 varid_pc = netcdf.inqVarID(nc, getVar); 0062 0063 data = netcdf.getVar(nc, varid_pc, 'single'); 0064 pc.(getVar).data = double(data); 0065 % Try to get some units (important for the calculation of MJD). 0066 try 0067 units = netcdf.getAtt(nc,varid_pc,'units'); 0068 catch 0069 units = []; 0070 end 0071 pc.(getVar).units = units; 0072 end 0073 0074 netcdf.close(nc) 0075 0076 % Extract the bathymetry ('pdepth' is cell thickness, 'depth' is cell 0077 % centre depth). 0078 nc = netcdf.open(polcoms_z, 'NOWRITE'); 0079 varid_pc = netcdf.inqVarID(nc, 'depth'); 0080 pc.depth.data = double(netcdf.getVar(nc, varid_pc, 'single')); 0081 try 0082 pc.depth.units = netcdf.getAtt(nc, varid_pc, 'units'); 0083 catch 0084 pc.depth.units = []; 0085 end 0086 varid_pc = netcdf.inqVarID(nc, 'pdepth'); 0087 pc.pdepth.data = double(netcdf.getVar(nc, varid_pc, 'single')); 0088 try 0089 pc.pdepth.units = netcdf.getAtt(nc, varid_pc, 'units'); 0090 catch 0091 pc.pdepth.units = []; 0092 end 0093 netcdf.close(nc) 0094 0095 % Data format: 0096 % 0097 % pc.ETW.data and pc.x1X.data are y, x, sigma, time 0098 % 0099 [~, ~, nz, nt] = size(pc.ETW.data); 0100 0101 % Make rectangular arrays for the nearest point lookup. 0102 [lon, lat] = meshgrid(pc.lon.data, pc.lat.data); 0103 0104 % Find the nearest POLCOMS point to each point in the FVCOM open boundaries 0105 fvlon = Mobj.lon(Mobj.obc_nodes(Mobj.obc_nodes ~= 0)); 0106 fvlat = Mobj.lat(Mobj.obc_nodes(Mobj.obc_nodes ~= 0)); 0107 0108 % Number of boundary nodes 0109 nf = sum(Mobj.nObcNodes); 0110 % Number of sigma layers. 0111 fz = size(Mobj.siglayz, 2); 0112 0113 % itemp = nan(nf, nz, nt); % POLCOMS interpolated temperatures 0114 % isal = nan(nf, nz, nt); % POLCOMS interpolated salinities 0115 fvtemp = nan(nf, fz, nt); % FVCOM interpolated temperatures 0116 fvsal = nan(nf, fz, nt); % FVCOM interpolated salinities 0117 0118 if ftbverbose 0119 tic 0120 end 0121 for t = 1:nt 0122 % Get the current 3D array of POLCOMS results. 0123 pctemp3 = pc.ETW.data(:, :, :, t); 0124 pcsal3 = pc.x1X.data(:, :, :, t); 0125 0126 % Preallocate the intermediate results arrays. 0127 itempz = nan(nf, nz); 0128 isalz = nan(nf, nz); 0129 idepthz = nan(nf, nz); 0130 0131 for j = 1:nz 0132 % Now extract the relevant layer from the 3D subsets. Transpose the 0133 % data to be (x, y) rather than (y, x). 0134 pctemp2 = pctemp3(:, :, j)'; 0135 pcsal2 = pcsal3(:, :, j)'; 0136 pcdepth2 = squeeze(pc.depth.data(:, :, j, t))'; 0137 0138 % Create new arrays which will be flattened when masking (below). 0139 tpctemp2 = pctemp2; 0140 tpcsal2 = pcsal2; 0141 tpcdepth2 = pcdepth2; 0142 tlon = lon; 0143 tlat = lat; 0144 0145 % Create and apply a mask to remove values outside the domain. This 0146 % inevitably flattens the arrays, but it shouldn't be a problem 0147 % since we'll be searching for the closest values in such a manner 0148 % as is appropriate for an unstructured grid (i.e. we're assuming 0149 % the POLCOMS data is irregularly spaced). 0150 mask = tpcdepth2 < -20000; 0151 tpctemp2(mask) = []; 0152 tpcsal2(mask) = []; 0153 tpcdepth2(mask) = []; 0154 % Also apply the masks to the position arrays so we can't even find 0155 % positions outside the domain, effectively meaning if a value is 0156 % outside the domain, the nearest value to the boundary node will 0157 % be used. 0158 tlon(mask) = []; 0159 tlat(mask) = []; 0160 0161 % Preallocate the intermediate results arrays. 0162 itempobc = nan(nf, 1); 0163 isalobc = nan(nf, 1); 0164 idepthobc = nan(nf, 1); 0165 0166 % Speed up the tightest loop with a parallelized loop. 0167 parfor i = 1:nf 0168 % Now we can do each position within the 2D layer. 0169 0170 fx = fvlon(i); 0171 fy = fvlat(i); 0172 0173 [~, ii] = sort(sqrt((tlon - fx).^2 + (tlat - fy).^2)); 0174 % Get the n nearest nodes from POLCOMS (more? fewer?). 0175 ixy = ii(1:16); 0176 0177 % Get the variables into static variables for the 0178 % parallelisation. 0179 plon = tlon(ixy); 0180 plat = tlat(ixy); 0181 ptemp = tpctemp2(ixy); 0182 psal = tpcsal2(ixy); 0183 pdepth = tpcdepth2(ixy); 0184 0185 % Use a triangulation to do the horizontal interpolation. 0186 tritemp = TriScatteredInterp(plon', plat', ptemp', 'natural'); 0187 trisal = TriScatteredInterp(plon', plat', psal', 'natural'); 0188 triz = TriScatteredInterp(plon', plat', pdepth', 'natural'); 0189 itempobc(i) = tritemp(fx, fy); 0190 isalobc(i) = trisal(fx, fy); 0191 idepthobc(i) = triz(fx, fy); 0192 0193 % Check all three, though if one is NaN, they all will be. 0194 if isnan(itempobc(i)) || isnan(isalobc(i)) || isnan(idepthobc(i)) 0195 warning('FVCOM boundary node at %f, %f is outside the POLCOMS domain. Setting to the closest POLCOMS value.', fx, fy) 0196 itempobc(i) = tpctemp2(ii(1)); 0197 isalobc(i) = tpcsal2(ii(1)); 0198 idepthobc(i) = tpcdepth2(ii(1)); 0199 end 0200 end 0201 0202 % Put the results in this intermediate array. 0203 itempz(:, j) = itempobc; 0204 isalz(:, j) = isalobc; 0205 idepthz(:, j) = idepthobc; 0206 end 0207 0208 % Now we've interpolated in space, we can interpolate the z-values 0209 % to the sigma depths. 0210 oNodes = Mobj.obc_nodes(Mobj.obc_nodes ~= 0); 0211 for zi = 1:fz 0212 0213 % Preallocate the output arrays 0214 fvtempz = nan(nf, fz); 0215 fvsalz = nan(nf, fz); 0216 0217 for pp = 1:nf 0218 % Get the FVCOM depths at this node 0219 tfz = Mobj.siglayz(oNodes(pp), :); 0220 % Now get the interpolated POLCOMS depth at this node 0221 tpz = idepthz(pp, :); 0222 0223 % Get the temperature and salinity values for this node and 0224 % interpolate down the water column (from POLCOMS to FVCOM). 0225 % TODO: Use csaps for the vertical interplation/subsampling at 0226 % each location. Alternatively, the pchip interp1 method seems 0227 % to do a decent job of the interpolation; it might be a more 0228 % suitable candidate in the absence of csaps. In fact, the demo 0229 % of csaps in the MATLAB documentation makes the interpolation 0230 % look horrible (shaving off extremes). I think pchip is 0231 % better. 0232 if ~isnan(tpz) 0233 fvtempz(pp, :) = interp1(tpz, itempz(pp, :), tfz, 'linear', 'extrap'); 0234 fvsalz(pp, :) = interp1(tpz, isalz(pp, :), tfz, 'linear', 'extrap'); 0235 else 0236 warning('Should never see this... ') % because we test for NaNs when fetching the values. 0237 warning('FVCOM boundary node at %f, %f is outside the POLCOMS domain. Skipping.', fvlon(pp), fvlat(pp)) 0238 continue 0239 end 0240 end 0241 end 0242 0243 % The horizontally-interpolated values in the final results array. 0244 % itemp(:, :, t) = itempz; 0245 % isal(:, :, t) = isalz; 0246 % The horizontally- and vertically-interpolated values in the final 0247 % FVCOM results array. 0248 fvtemp(:, :, t) = fvtempz; 0249 fvsal(:, :, t) = fvsalz; 0250 end 0251 if ftbverbose 0252 toc 0253 end 0254 0255 Mobj.temperature = fvtemp; 0256 Mobj.salt = fvsal; 0257 0258 % Do we have to interpolate to the FVCOM time series? Looking at page 325 0259 % of the FVCOM manual, it looks like the temperature and salinity are on a 0260 % different time sampling from the other example files (14 time steps vs. 0261 % 3625 for _wnd.nc or 43922 for _julain_obc.nc (i.e. surface elevation at 0262 % the boundary)). That's not to say those files were all used in the same 0263 % model run... In the interim, just convert the current times to Modified 0264 % Julian Day (this is a bit ugly). 0265 % pc.time.yyyymmdd = strtrim(regexp(pc.time.units, 'since', 'split')); 0266 % pc.time.yyyymmdd = str2double(regexp(pc.time.yyyymmdd{end}, '-', 'split')); 0267 % Mobj.ts_times = greg2mjulian(pc.time.yyyymmdd(1), pc.time.yyyymmdd(2), pc.time.yyyymmdd(3), 0, 0, 0) + (pc.time.data / 3600 / 24); 0268 0269 % Convert the POLCOMS times to Modified Julian Day (this is a very ugly). 0270 pc.time.yyyymmdd = strtrim(regexp(pc.time.units, 'since', 'split')); 0271 pc.time.strtime = regexp(pc.time.yyyymmdd{end}, '-', 'split'); 0272 % This new version of the time has the year in a weird format (yr.#). We 0273 % thus need to split it again to get the decimal year (post-2000 only?). 0274 pc.time.strtimeyr = regexp(pc.time.strtime, '\.', 'split'); 0275 pc.time.yyyymmdd = str2double([pc.time.strtimeyr{1}(2), pc.time.strtime{2:3}]); 0276 pc.time.yyyymmdd(1) = pc.time.yyyymmdd(1) + 2000; % add full year. 0277 Mobj.ts_times = greg2mjulian(pc.time.yyyymmdd(1), pc.time.yyyymmdd(2), pc.time.yyyymmdd(3), 0, 0, 0) + (pc.time.data / 3600 / 24); 0278 0279 if ftbverbose 0280 fprintf(['end : ' subname '\n']) 0281 end