snippet.matlab
% 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);
  • configuration/make_restart.m.txt
  • Last modified: 2017/04/18 09:11
  • by pica