Commit 728b457f authored by Pierre Cazenave's avatar Pierre Cazenave
Browse files

Propagate the changes from get_POLCOMS_tsobc.m to this file to support reading...

Propagate the changes from get_POLCOMS_tsobc.m to this file to support reading in multiple NetCDF files and concatenating their data along the time axis. This also includes a section to replace the times in the restart file with an arbitrary time period, but this is commented out. This is useful if the restart file does not belong to the same time period as the data and you essentially want the restart file as a template only
parent b432f085
......@@ -22,12 +22,12 @@ function replace_FVCOM_restart_vars(Mobj, polcoms_ts, polcoms_z, start_date, fv_
% - Mobj.siglayz - sigma layer depths for all model
% nodes.
% - Mobj.lon, Mobj.lat - node coordinates (long/lat).
% 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.
% 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.
% start_date = Gregorian start date array (YYYY, MM, DD, hh, mm, ss).
% fv_restart = full path to the FVCOM restart file.
% varlist = cell array of variables to extract from the POLCOMS
......@@ -49,6 +49,7 @@ function replace_FVCOM_restart_vars(Mobj, polcoms_ts, polcoms_z, start_date, fv_
%
% Revision history
% 2013-01-28 First version (based on initial_tsobc.m).
% 2013-02-05 Add support for reading in multiple NetCDF files.
%
%==========================================================================
......@@ -64,49 +65,131 @@ end
% Extract the POLCOMS data specified in varlist
%--------------------------------------------------------------------------
if ftbverbose
fprintf('%s : read POLCOMS data... ', subname)
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
nc = netcdf.open(polcoms_ts, 'NOWRITE');
for ii = 1:todo
for var=1:numel(varlist)
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\n', subname, tmp_fn)
else
fprintf('%s: extracting file %s (%i of %i)... ', subname, tmp_fn, ii, todo)
end
end
getVar = varlist{var};
varid_pc = netcdf.inqVarID(nc, getVar);
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)
data = netcdf.getVar(nc, varid_pc, 'single');
pc.(getVar).data = double(data);
% Try to get some units (important for the calculation of MJD).
% 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
units = netcdf.getAtt(nc,varid_pc,'units');
catch err
units = [];
pc.depth.units = netcdf.getAtt(nc, varid_pc, 'units');
catch
pc.depth.units = [];
end
pc.(getVar).units = units;
end
clear data getVar varid_pc var
clear data
netcdf.close(nc)
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
% Extract the bathymetry ('pdepth' is cell thickness, 'depth' is cell
% centre depth).
nc = netcdf.open(polcoms_z, 'NOWRITE');
varid_pc = netcdf.inqVarID(nc, 'depth');
pc.depth.data = double(netcdf.getVar(nc, varid_pc, 'single'));
try
pc.depth.units = netcdf.getAtt(nc, varid_pc, 'units');
catch err
pc.depth.units = [];
end
varid_pc = netcdf.inqVarID(nc, 'pdepth');
pc.pdepth.data = double(netcdf.getVar(nc, varid_pc, 'single'));
try
pc.pdepth.units = netcdf.getAtt(nc, varid_pc, 'units');
catch err
pc.pdepth.units = [];
end
netcdf.close(nc)
% Data format:
%
......@@ -122,13 +205,16 @@ netcdf.close(nc)
% 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.
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);
% 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);
% Given our intput time (in start_date), find the nearest time
% index for the POLCOMS data.
......@@ -155,10 +241,9 @@ pctempz = nan(nx, ny, fz);
pcsalz = nan(nx, ny, fz);
for xi = 1:nx
% For the parallel loop, get all the y and z dimension data for the
% current x position (temperature, salinity and depth).
% N.B. The POLCOMS data is stored y, x, z, t (first two dimensions the
% wrong way around).
% Get all the y and z dimension data for the current x position
% (temperature, salinity and depth). N.B. The POLCOMS data is stored y,
% x, z, t (first two dimensions the wrong way around).
xtemp = squeeze(pc.ETW.data(:, xi, :, tidx));
xsalt = squeeze(pc.x1X.data(:, xi, :, tidx));
xdepth = squeeze(pc.depth.data(:, xi, :, tidx));
......@@ -379,9 +464,24 @@ for ii = 1:numvars
sfvdata(:, :, tt) = startdata + (ss(tt) .* td);
end
end
% Replace the values with the scaled interpolated values.
netcdf.putVar(ncout, varid, sfvdata)
% % We might also want to replace the time. If so, uncomment these
% % lines to replace with an arbitrary time period.
% elseif strcmpi(varname, 'time')
% tmp_start_time = greg2mjulian(start_date(1), start_date(2), start_date(3) - 7, start_date(4), start_date(5), start_date(6));
% tmp_time = tmp_start_time:(tmp_start_time + nt - 1);
% netcdf.putVar(ncout, varid, tmp_time)
% elseif strcmpi(varname, 'Times')
% tmp_time = [];
% for i = 1:nt;
% tmp_time = [tmp_time, sprintf('%-026s', datestr(datenum(start_date) - 7 + (i - 1), 'yyyy-mm-dd HH:MM:SS.FFF'))];
% end
% netcdf.putVar(ncout, varid, tmp_time)
% elseif strcmpi(varname, 'Itime')
% tmp_start_time = greg2mjulian(start_date(1), start_date(2), start_date(3) - 7, start_date(4), start_date(5), start_date(6));
% tmp_time = tmp_start_time:(tmp_start_time + nt - 1);
% netcdf.putVar(ncout, varid, floor(tmp_time))
else
% We need to check if the dimension is unlimited, and use a start
% and end with netcdf.putVar if it is. This is largely because
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment