% Create a new restart file from a donor template file. The donor
% restart file has a single time step only and is generated with
% an ncks cut of a single time step from an FVCOM restart file:
% ncks -dtime,0 source_restart.nc output_restart.nc
%
% Pierre Cazenave (Plymouth Marine Laboratory)
matlabrc
close all
clc
global ftbverbose
ftbverbose = 1; % be noisy
% Set up the MATLAB fvcom-toolbox.
addpath('/home/user/fvcom-toolbox/utilities')
addpath('/home/user/fvcom-toolbox/fvcom_prepro/')
conf.base = '/home/user/run/';
% Case name for the model inputs and outputs
conf.casename = 'casename_v01';
% Model time ([Y,M,D,h,m,s])
conf.year = 2005;
conf.month = 01;
% Set start time for the single time step restart file.
conf.startDate = [conf.year, conf.month, 01, 00, 00, 00];
conf.donor = fullfile(conf.base, 'output', ...
conf.casename, ...
sprintf('%04d', conf.year), ...
sprintf('%02d', conf.month), ...
sprintf('%s_restart_donor_0001.nc', conf.casename));
%% Read the donor restart file for time and spatial data.
vars = {'lon', 'lat', 'x', 'y', 'nv', ...
'siglay', 'siglev', 'h', ...
'time', 'Times', ...
'salinity', 'temp'};
for v = 1:length(vars)
if strcmpi(vars{v}, 'Times')
Mobj.(vars{v}) = ncread(conf.donor, vars{v});
else
Mobj.(vars{v}) = double(ncread(conf.donor, vars{v}));
end
end
clear v
Mobj.have_lonlat = true;
% Fix the Times
Mobj.Times = Mobj.Times';
Mobj.ts_times = Mobj.time;
% Make depth resolved sigma depths.
Mobj.siglayz = Mobj.siglay .* ...
repmat(Mobj.h, [1, size(Mobj.siglay, 2)]);
[conf.year, conf.month, conf.day, ~, ~, ~] = ...
mjulian2greg(Mobj.time);
conf.dates = [conf.year, conf.month, conf.day, 0, 0, 0];
% Round to the nearest month.
month = conf.month + sign(conf.day - 15 + double(~(conf.month - 2)))...
+ double(~(conf.day - 15 + double(~(conf.month - 2))));
conf.coarse.year = conf.year + double((month - 12) == 1) - ...
double((1 - month) == 1);
conf.coarse.month = mod(month, 12) + 12 * double(~mod(month, 12));
%% Prepare the data
% Use HYCOM model output to interpolate the restart file.
coarse = get_HYCOM_forcing(Mobj, modelTime, ...
{'temperature', 'salinity'});
coarse.lon = double(coarse.lon);
coarse.lat = double(coarse.lat);
% Interpolate the coarse data on the FVCOM grid. Do temperature/salinity
% only since velocities are non-tidal in most coarse models.
varlist = {'lon', 'lat', 'temperature', 'salinity'};
Mobj = interp_HYCOM2FVCOM(Mobj, coarse, conf.dates, varlist);
% Sometimes we end up with ridiculous values in the sea which
% are land values. We need to replace them with something much more
% sensible here. Use the same approach as from interp_HYCOM2FVCOM where we
% just use the nearest sensible value.
for vv = 1:length(varlist)
currvar = varlist{vv};
switch currvar
case {'lon', 'lat'}
continue
otherwise
% We've got data of some sort. Assume we don't want
% anything above 50 units (works for temperature and
% salinity as well as velocity; won't work for things
% like density).
threshold = 50; % a threshold with which you're comfortable
if max(Mobj.restart.(currvar)(:)) > 50
warning(['We''ve got stupidly high ', ...
'''Mobj.restart.%s'' values. Fixing ', ...
'those above %f here.'], currvar)
fvtmp = Mobj.restart.(currvar);
badmask = Mobj.restart.(currvar) > 100;
fvtmp(badmask) = nan;
[nx, nz] = size(fvtmp);
fvidx = 1:nx;
flon = Mobj.lon;
flat = Mobj.lat;
for zi = 1:nz
fvnanidx = fvidx(isnan(fvtmp(:, zi)));
fvfinidx = fvidx(~isnan(fvtmp(:, zi)));
for ni = fvnanidx
% Find the nearest non-nan data.
xx = flon(ni);
yy = flat(ni);
[~, di] = min(sqrt((flon(fvfinidx) - xx).^2 + ...
(flat(fvfinidx) - yy).^2));
fvtmp(ni, zi) = fvtmp(fvfinidx(di), zi);
end
end
Mobj.restart.(currvar) = fvtmp;
clear fvtmp badmask nx nz fvidx flon flag zi ...
fvnanidx fvfinidx ni xx yy di
end
end
end
clear currvar vv
% Rename the restart field names with FVCOM variable names. This
% approach is a bit overkill since most are the same, but it's more
% easily expanded if other variables need to be interpolated later.
replacements.temperature = 'temp';
replacements.salinity = 'salinity';
fields = fieldnames(Mobj.restart);
for ff = 1:length(fields)
% Don't change the fields if they're the same name.
if isfield(replacements, fields{ff}) && ...
~strcmpi(fields{ff}, replacements.(fields{ff}))
Mobj.restart.(replacements.(fields{ff})) = ...
Mobj.restart.(fields{ff});
Mobj.restart = rmfield(Mobj.restart, fields{ff});
end
end
% Have a look see at the interpolated data.
% figure(100)
% clf
% patch('Vertices', [Mobj.lon, Mobj.lat], 'Faces', Mobj.nv,...
% 'Cdata', Mobj.restart.salinity(:, 1), 'edgecolor', 'none', ...
% 'facecolor', 'interp');
% colorbar
% axis('square', 'tight')
% title('Surface salinity')
% caxis([33, 36])
%% Write out all the required files.
conf.outbase = fullfile(conf.base, 'input', conf.casename, ...
sprintf('%04d', conf.year), ...
sprintf('%02d', conf.month));
write_FVCOM_restart(conf.donor, ...
fullfile(conf.outbase, ...
sprintf('%s_%04d-%02d_%s_restart.nc', ...
conf.casename, conf.year, conf.month, conf.source)), ...
Mobj.restart);