Commit 64766a18 authored by Ricardo Torres's avatar Ricardo Torres 💬

Merge branch 'dev' of gitlab.ecosystem-modelling.pml.ac.uk:fvcom/fvcom-toolbox into dev

parents 1e87475b 66d4a2d4
function [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType,plotFig)
function Mobj = add_obc_nodes_list(Mobj, Nodelist, ObcName, ObcType, plotFig)
% Add a set of obc nodes comprising a single obc boundary to Mesh structure
% Using a list of nodes
% using a list of nodes.
%
% [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType)
% Mobj = add_obc_nodes_list(Mobj, Nodelist, ObcName, ObcType)
%
% DESCRIPTION:
% Select using ginput the set of nodes comprising an obc
% Add a set of open boundary nodes for a given boundary using a list of
% nodes.
%
% INPUT
% Mobj = Matlab mesh object
% Nlist = List of nodes
% Mobj = Matlab mesh object with fields:
% - nVerts - number of nodes in the domain
% - nativeCoords - coordinates in which the model runs (only for
% plotting the figure).
% - x, y, lon, lat - mesh node coordinates (either cartesian or
% spherical) (only for plotting the figure).
% - tri - model grid triangulation (only for plotting the figure)
% - h - model grid depths (only for plotting the figure)
% Nodelist = List of nodes
% ObcName = Name of the Open Boundary
% ObcType = FVCOM Flag for OBC Type
% plotFig = [optional] show a figure of the mesh (1 = yes)
......@@ -19,7 +26,8 @@ function [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType,plotFig)
% Mobj = Matlab mesh object with an additional obc nodelist
%
% EXAMPLE USAGE
% Mobj = add_obc_nodes_list(Mobj,Nlist,'OpenOcean')
% Nodelist = 1:100;
% Mobj = add_obc_nodes_list(Mobj, Nodelist, 'OpenOcean')
%
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
......@@ -32,10 +40,10 @@ function [Mobj] = add_obc_nodes_list(Mobj,Nlist,ObcName,ObcType,plotFig)
% prevent it from sorting the values it returns. Amended by Pierre to
% support pre-2012 versions of MATLAB whilst giving the same result.
% 2015-02-23 Output number of nodes if the verbose flag is set.
% 2017-08-31 Update the help to clarify what's needed.
%
%==========================================================================
subname = 'add_obc_nodes_list';
[~, subname] = fileparts(mfilename('fullpath'));
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s\n', subname)
......@@ -51,13 +59,10 @@ end
%--------------------------------------------------------------------------
% Make this works in versions of MATLAB older than 2012a (newer versions
% can just use unique(A, 'stable'), but checking versions is a pain).
[~, Nidx] = unique(Nlist);
Nlist = Nlist(sort(Nidx));
[~, Nidx] = unique(Nodelist);
Nodelist = Nodelist(sort(Nidx));
if max(Nlist) > Mobj.nVerts
fprintf('Your open boundary node number exceed the total number of nodes in the domain\n');
error('stopping...')
end
assert(max(Nodelist) < Mobj.nVerts, 'Your open boundary node number exceed the total number of nodes in the domain\n')
%--------------------------------------------------------------------------
% Plot the mesh
......@@ -76,16 +81,16 @@ if plotFig == 1
'Cdata', Mobj.h, 'edgecolor', 'k', 'facecolor', 'interp')
hold on
whos Nlist
plot(x(Nlist), y(Nlist), 'ro');
plot(x(Nodelist), y(Nodelist), 'ro');
axis('equal', 'tight')
title('open boundary nodes');
end
% add to mesh object
npts = numel(Nlist);
npts = numel(Nodelist);
Mobj.nObs = Mobj.nObs + 1;
Mobj.nObcNodes(Mobj.nObs) = npts;
Mobj.obc_nodes(Mobj.nObs,1:npts) = Nlist;
Mobj.obc_nodes(Mobj.nObs,1:npts) = Nodelist;
Mobj.obc_name{Mobj.nObs} = ObcName;
Mobj.obc_type(Mobj.nObs) = ObcType;
......
......@@ -27,27 +27,25 @@ function [Mobj] = estimate_ts(Mobj,u,zeta)
% latitude and longitudes. Also add arguments to the function to define
% current velocity and tidal amplitudes.
% 2016-02018 Update the help to reflect the required variables.
% 2017-09-01 Tidy the code a bit.
%
%==============================================================================
subname = 'estimate_ts';
[~, subname] = fileparts(mfilename('fullpath'));
global ftbverbose
if(ftbverbose)
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end;
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
%------------------------------------------------------------------------------
% Set constants
%------------------------------------------------------------------------------
g = 9.81; %gravitational acceleration
%u = 3.0; %u-velocity
%zeta = 11.0; %tide amp
g = 9.81; % gravitational acceleration
if(~Mobj.have_bath)
error('can''t estimate the time step without bathymetry')
end;
end
%------------------------------------------------------------------------------
% Compute the time step estimate
......@@ -83,15 +81,17 @@ for i=1:nElems
dpth = max(dpth,1);
ts(nds) = min(ts(nds),lside(i)/(sqrt(g*dpth) + u));
end
if(ftbverbose); fprintf('minimum time step: %f seconds\n',min(ts)); end;
if ftbverbose
fprintf('minimum time step: %f seconds\n',min(ts));
end
Mobj.ts = ts;
Mobj.have_ts = true;
Mobj.lside=lside;
if(ftbverbose)
fprintf(['end : ' subname '\n'])
if ftbverbose
fprintf('end : %s\n', subname)
end
function [km]=haversine(lat1,lon1,lat2,lon2)
function km=haversine(lat1,lon1,lat2,lon2)
% Haversine function to calculate first order distance measurement. Assumes
% spherical Earth surface. Lifted from:
%
......
......@@ -120,7 +120,6 @@ hdx = max([max(diff(lon(:,1))),max(diff(lat(1,:)))]);
% Number of sigma layers.
fz = size(Mobj.siglayz, 2);
% Make a 3D array of the HYCOM depths and mask where we don't have data.
% This can then be used in the interpolation instead of trying to deal with
% this on the fly.
......@@ -483,7 +482,7 @@ for v = 1:length(fields)
% Find and remove NaNs.
parfor pp=1:fz
test = fvtempz(:,pp);
test = fvtempz(:, pp);
if any(isnan(test))
igood = ~isnan(test);
ftri = scatteredInterpolant(fvlon(igood), fvlat(igood), test(igood), 'nearest', 'nearest');
......
......@@ -344,7 +344,7 @@ for ff = 1:fv_nr
% the original river location (in fvcom_xy).
[~, n_idx] = sort(sqrt( ...
(nemo_xy(ff, 1) - Mobj.lon(n_tri)).^2 ...
+ (nemo_xy(ff, 2) - Mobj.lon(n_tri)).^2));
+ (nemo_xy(ff, 2) - Mobj.lat(n_tri)).^2));
[row_2, ~] = find(Mobj.tri == n_tri(n_idx(1)));
if length(n_idx) > 1
......
function Mobj = interp_sst_assimilation(Mobj, conf, output_file)
% Interpolate SST values to a given FVCOM mesh object.
%
% Mobj = interp_sst_assimilation(Mobj, year, sst_dir, file_pattern)
%
% DESCRIPTION:
% Interpolate SST data from remote sensing data onto the supplied model
% grid.
%
% INPUT:
% Mobj - MATLAB mesh object containing fields:
% lon, lat - node coordinates (spherical)
% lonc, latc - element coordinates (spherical)
% tri - element triangulation table
% conf - struct with fields:
% sst_dir - directory containing the SST data
% sst_pattern - file name pattern for the SST data
% year - year for which to generate SST data
% output_file - path to which to output the netCDF
%
% OUTPUT:
% FVCOM data assimilation SST netCDF file.
% Mobj - input MATLAB mesh object with added 'assim.sst' field, with
% fields:
% data - the SST data
% time - the SST time series
%
% EXAMPLE USAGE:
% Mobj = read_fvcom_mesh('casename.grd');
% conf.sst_dir = '/home/user/GHRSST/';
% conf.sst_pattern = '-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc';
% interp_sst_assimilation(Mobj, conf, 'casename_sstgrd.nc');
%
% Author(s)
% Ricardo Torres (Plymouth Marine Laboratory)
% Pierre Cazenave (Plymouth Marine Laboratory)
%
% History:
% 2017-08-02 Turned the script into a function.
[~, subname] = fileparts(mfilename('fullpath'));
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s \n', subname)
end
year = conf.year;
% SST files.
sst_dir = conf.sst_dir;
sst_pattern = conf.sst_pattern;
time_span = datenum(year-1,12,31):datenum(year+1,01,01);
file_list = cell(length(time_span), 1);
% Add last day from the previous year and the first day of the following
% year.
for aa = 1:length(time_span)
file_list{aa} = fullfile(sst_dir, ...
datestr(time_span(aa),'yyyy'), ...
sprintf(sst_pattern, ...
str2num(datestr(time_span(aa), 'yyyy')), ...
str2num(datestr(time_span(aa), 'mm')), ...
str2num(datestr(time_span(aa), 'dd'))));
if ~exist(file_list{aa}, 'file')
error('We are missing a file (%s) from this year (%04d)', ...
file_list{aa}, datestr(time_span(aa),'yyyy'))
end
end
% Read SST data files and interpolate each to the FVCOM mesh
lon = ncread(file_list{1},'lon');
lat = ncread(file_list{1},'lat');
mask = ncread(file_list{1},'mask');
[lonm,latm]=meshgrid(lon,lat);
lonm=lonm';
latm=latm';
lonm = lonm(mask==1);
latm = latm(mask==1);
time = zeros(length(file_list),1);
sst = zeros(Mobj.nVerts,length(file_list),1,'single');
fvcomlon = Mobj.lon;
fvcomlat = Mobj.lat;
if license('test', 'Distrib_Computing_Toolbox')
if isempty(gcp('nocreate'))
% Force pool to be local in case we have remote pools available.
parpool('local');
end
end
if ftbverbose
fprintf('Progress:\n');
fprintf([repmat('.', 1, length(file_list)), '\n\n']);
end
parfor ff = 1:length(file_list)
if ftbverbose
fprintf('\b|\n');
end
sst_eo = ncread(file_list{ff}, 'analysed_sst') - 273.15;
mask = ncread(file_list{ff}, 'mask');
lon = ncread(file_list{ff}, 'lon');
lat = ncread(file_list{ff}, 'lat');
[lonm, latm] = meshgrid(lon, lat);
lonm = lonm';
latm = latm';
lonm = lonm(mask == 1);
latm = latm(mask == 1);
time_eo = ncread(file_list{ff}, 'time');
time_eo_units = ncreadatt(file_list{ff}, 'time', 'units');
t0str = textscan(time_eo_units, 'seconds since %s%s');
t0 = datenum([strtrim(t0str{1}{1}), strtrim(t0str{2}{1})], 'yyyy-mm-ddHH:MM:SS');
time_out = (t0 + double(time_eo/(60*60*24)));
sst_eo = sst_eo(mask == 1);
% Build interpolant
ft = scatteredInterpolant(double(lonm), double(latm), sst_eo, 'nearest', 'linear');
sst(:,ff) = ft(fvcomlon, fvcomlat);
time(ff) = time_out + 0.5; % fvcom expects these to be at mid-day
end
ntimes = length(time);
% Do the times.
[sYr, sMon, sDay, sHr, sMin, sSec] = datevec(time);
MJDtime = greg2mjulian(sYr, sMon, sDay, sHr, sMin, sSec);
% Create netCDF file
nc = netcdf.create(output_file, 'clobber');
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'year', num2str(year))
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'title','FVCOM SST 1km merged product File')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'institution','Plymouth Marine Laboratory')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'source','FVCOM grid (unstructured) surface forcing')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'history', sprintf('File created with %s from the MATLAB fvcom-toolbox', subname))
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'references','http://fvcom.smast.umassd.edu, http://codfish.smast.umassd.edu')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.0')
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'CoordinateSystem',Mobj.nativeCoords)
netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'CoordinateProjection','init=WGS84') % WGS84?
% Dimensions
nele_dimid=netcdf.defDim(nc,'nele',Mobj.nElems);
node_dimid=netcdf.defDim(nc,'node',Mobj.nVerts);
three_dimid=netcdf.defDim(nc,'three',3);
time_dimid=netcdf.defDim(nc,'time',netcdf.getConstant('NC_UNLIMITED'));
datestrlen_dimid=netcdf.defDim(nc,'DateStrLen',26);
% Space variables
lon_varid=netcdf.defVar(nc,'lon','NC_FLOAT',node_dimid);
netcdf.putAtt(nc,lon_varid,'long_name','nodal longitude');
netcdf.putAtt(nc,lon_varid,'units','degrees_easdt');
lat_varid=netcdf.defVar(nc,'lat','NC_FLOAT',node_dimid);
netcdf.putAtt(nc,lat_varid,'long_name','nodal latitude');
netcdf.putAtt(nc,lat_varid,'units','degrees_north');
lonc_varid=netcdf.defVar(nc,'lonc','NC_FLOAT',nele_dimid);
netcdf.putAtt(nc,lonc_varid,'long_name','zonal longitude');
netcdf.putAtt(nc,lonc_varid,'units','meters');
latc_varid=netcdf.defVar(nc,'latc','NC_FLOAT',nele_dimid);
netcdf.putAtt(nc,latc_varid,'long_name','zonal latitude');
netcdf.putAtt(nc,latc_varid,'units','meters');
nv_varid=netcdf.defVar(nc,'nv','NC_INT',[nele_dimid, three_dimid]);
netcdf.putAtt(nc,nv_varid,'long_name','nodes surrounding element');
% Time variables
time_varid=netcdf.defVar(nc,'time','NC_FLOAT',time_dimid);
netcdf.putAtt(nc,time_varid,'long_name','time');
netcdf.putAtt(nc,time_varid,'units','days since 1858-11-17 00:00:00');
netcdf.putAtt(nc,time_varid,'delta_t','0000-00-00 01:00:00')
netcdf.putAtt(nc,time_varid,'format','modified julian day (MJD)');
netcdf.putAtt(nc,time_varid,'time_zone','UTC');
times_varid=netcdf.defVar(nc,'Times','NC_CHAR',[datestrlen_dimid,time_dimid]);
netcdf.putAtt(nc,times_varid,'long_name','Calendar Date');
netcdf.putAtt(nc,times_varid,'format','String: Calendar Time');
netcdf.putAtt(nc,times_varid,'time_zone','UTC');
sst_varid = netcdf.defVar(nc, 'sst', 'NC_FLOAT', [node_dimid, time_dimid]);
netcdf.putAtt(nc, sst_varid, 'long_name', 'sea surface Temperature');
netcdf.putAtt(nc, sst_varid, 'units', 'Celsius Degree');
netcdf.putAtt(nc, sst_varid, 'grid', 'fvcom_grid');
netcdf.putAtt(nc, sst_varid, 'coordinates', Mobj.nativeCoords);
netcdf.putAtt(nc, sst_varid, 'type', 'data');
% End definitions
netcdf.endDef(nc);
% Put the easy ones in first.
netcdf.putVar(nc, nv_varid, Mobj.tri);
netcdf.putVar(nc,lon_varid,Mobj.lon);
netcdf.putVar(nc,lat_varid,Mobj.lat);
netcdf.putVar(nc,lonc_varid,Mobj.lonc);
netcdf.putVar(nc,latc_varid,Mobj.latc);
netcdf.putVar(nc,time_varid,0,ntimes,MJDtime);
nStringOut = char();
[nYr, nMon, nDay, nHour, nMin, nSec] = mjulian2greg(MJDtime);
for i=1:ntimes
nDate = [nYr(i), nMon(i), nDay(i), nHour(i), nMin(i), nSec(i)];
nStringOut = [nStringOut, sprintf('%04i/%02i/%02i %02i:%02i:%09.6f', nDate)];
end
netcdf.putVar(nc,times_varid,[0, 0], [26, ntimes], nStringOut);
netcdf.putVar(nc, sst_varid, [0, 0], [Mobj.nVerts, ntimes], sst)
fprintf('done.\n')
% Close the netCDF file(s)
netcdf.close(nc);
Mobj.assim.sst.data = sst;
Mobj.assim.sst.time = time;
% Plot sst if needed.
% for aa = 1:size(sst,2)
% plot_field(Mobj,sst(:,aa))
% pause
% end
if ftbverbose
fprintf('end : %s \n', subname)
end
......@@ -420,6 +420,7 @@ end
%%
nml.NML_NESTING = [];
nml.NML_NESTING.NESTING_ON = 'F';
nml.NML_NESTING.FABM_NESTING_ON = 'F';
nml.NML_NESTING.NESTING_BLOCKSIZE= 10;
fmt.NML_NESTING.NESTING_BLOCKSIZE.format= '%u';
nml.NML_NESTING.NESTING_TYPE= 3;
......
......@@ -92,7 +92,11 @@ ObcNodes = tmpObcNodes(tmpObcNodes~=0)';
%--------------------------------------------------------------------------
% Sanity check on input and dimensions
%--------------------------------------------------------------------------
nTimes = numel(MJD);
if ischar(MJD(1))
nTimes = size(MJD, 1);
else
nTimes = numel(MJD);
end
if ftbverbose; fprintf('Number of time steps %d\n',nTimes); end
nObcs = numel(ObcNodes);
......@@ -166,6 +170,28 @@ netcdf.endDef(nc);
% write data
netcdf.putVar(nc,nobc_varid,ObcNodes);
netcdf.putVar(nc,iint_varid,0,nTimes,1:nTimes);
if strtime
% If out MJD data is characters, assume we've already got a suitable
% array of Time strings. Use those to create an MJD array to write to
% netCDF. This is sometimes preferable to having MJD as an array of
% floats in the case where we've read in a 'time' variable from a
% netCDF file and its precision is insufficient to actually store the
% times properly. netCDF, otherwise, create one assuming we've actually
% got Modified Julian Days. If we've been given an array of floats,
% then just dump those to netCDF as before.
if ischar(MJD(1))
nStringOut = MJD';
MJD = datenum(nStringOut', 'YYYY-mm-dd HH:MM:SS.FFF') - 678942;
else
nStringOut = char();
[nYr, nMon, nDay, nHour, nMin, nSec] = mjulian2greg(MJD);
for i=1:nTimes
nDate = [nYr(i), nMon(i), nDay(i), nHour(i), nMin(i), nSec(i)];
nStringOut = [nStringOut, sprintf('%04i/%02i/%02i %02i:%02i:%09.6f', nDate)];
end
end
netcdf.putVar(nc,Times_varid,nStringOut);
end
if floattime
netcdf.putVar(nc,time_varid,0,nTimes,MJD);
end
......@@ -173,15 +199,6 @@ if inttime
netcdf.putVar(nc,itime_varid,floor(MJD));
netcdf.putVar(nc,itime2_varid,0,nTimes,round(mod(MJD,1) * 24 * 60 * 60 * 1000));
end
if strtime
nStringOut = char();
[nYr, nMon, nDay, nHour, nMin, nSec] = mjulian2greg(MJD);
for i=1:nTimes
nDate = [nYr(i), nMon(i), nDay(i), nHour(i), nMin(i), nSec(i)];
nStringOut = [nStringOut, sprintf('%04i/%02i/%02i %02i:%02i:%09.6f', nDate)];
end
netcdf.putVar(nc,Times_varid,nStringOut);
end
netcdf.putVar(nc,elevation_varid,Mobj.surfaceElevation);
% close file
......
......@@ -323,8 +323,8 @@ siglay_varid = netcdf.defVar(nc, 'siglay', 'NC_FLOAT', ...
netcdf.putAtt(nc, siglay_varid, 'long_name', 'Sigma Layers');
netcdf.putAtt(nc, siglay_varid, 'standard_name', 'ocean_sigma/general_coordinate');
netcdf.putAtt(nc, siglay_varid, 'positive', 'up');
netcdf.putAtt(nc, siglay_varid, 'valid_min', '-1.f');
netcdf.putAtt(nc, siglay_varid, 'valid_max', '0.f');
netcdf.putAtt(nc, siglay_varid, 'valid_min', -1);
netcdf.putAtt(nc, siglay_varid, 'valid_max', 0);
netcdf.putAtt(nc, siglay_varid, 'formula_terms', 'sigma: siglay eta: zeta depth: h');
siglayc_varid = netcdf.defVar(nc, 'siglay_center', 'NC_FLOAT', ...
......@@ -332,8 +332,8 @@ siglayc_varid = netcdf.defVar(nc, 'siglay_center', 'NC_FLOAT', ...
netcdf.putAtt(nc, siglayc_varid, 'long_name', 'Sigma Layers');
netcdf.putAtt(nc, siglayc_varid, 'standard_name', 'ocean_sigma/general_coordinate');
netcdf.putAtt(nc, siglayc_varid, 'positive', 'up');
netcdf.putAtt(nc, siglayc_varid, 'valid_min', '-1.f');
netcdf.putAtt(nc, siglayc_varid, 'valid_max', '0.f');
netcdf.putAtt(nc, siglayc_varid, 'valid_min', -1);
netcdf.putAtt(nc, siglayc_varid, 'valid_max', 0);
netcdf.putAtt(nc, siglayc_varid, 'formula_terms', 'sigma: siglay_center eta: zeta_center depth: h_center');
siglev_varid = netcdf.defVar(nc, 'siglev', 'NC_FLOAT', ...
......@@ -341,8 +341,8 @@ siglev_varid = netcdf.defVar(nc, 'siglev', 'NC_FLOAT', ...
netcdf.putAtt(nc, siglev_varid, 'long_name', 'Sigma Levels');
netcdf.putAtt(nc, siglev_varid, 'standard_name', 'ocean_sigma/general_coordinate');
netcdf.putAtt(nc, siglev_varid, 'positive', 'up');
netcdf.putAtt(nc, siglev_varid, 'valid_min', '-1.f');
netcdf.putAtt(nc, siglev_varid, 'valid_max', '0.f');
netcdf.putAtt(nc, siglev_varid, 'valid_min', -1);
netcdf.putAtt(nc, siglev_varid, 'valid_max', 0);
netcdf.putAtt(nc, siglev_varid, 'formula_terms', 'sigma:siglev eta: zeta depth: h');
siglevc_varid = netcdf.defVar(nc, 'siglev_center', 'NC_FLOAT', ...
......@@ -350,8 +350,8 @@ siglevc_varid = netcdf.defVar(nc, 'siglev_center', 'NC_FLOAT', ...
netcdf.putAtt(nc, siglevc_varid, 'long_name', 'Sigma Layers');
netcdf.putAtt(nc, siglevc_varid, 'standard_name', 'ocean_sigma/general_coordinate');
netcdf.putAtt(nc, siglevc_varid, 'positive', 'up');
netcdf.putAtt(nc, siglevc_varid, 'valid_min', '-1.f');
netcdf.putAtt(nc, siglevc_varid, 'valid_max', '0.f');
netcdf.putAtt(nc, siglevc_varid, 'valid_min', -1);
netcdf.putAtt(nc, siglevc_varid, 'valid_max', 0);
netcdf.putAtt(nc, siglevc_varid, 'formula_terms', 'sigma: siglev_center eta: zeta_center depth: h_center');
h_varid = netcdf.defVar(nc, 'h', 'NC_FLOAT', ...
......
......@@ -65,13 +65,16 @@ function write_FVCOM_restart(fv_restart, out_restart, indata, varargin)
% 2013-02-15 Fix bug wherein only the last field in the new data would
% only be added to the output netCDF file.
% 2013-03-13 Make the time rewriting optional and not just commented out.
% 2014-02-04 Incorporate Karen's functionality (see revision history
% below), but with the ability to retain the existing behaviour (where a
% new start time is still optional). User-specified constants are also
% supported but instead of specifying a new input argument, if a single
% scalar value is given in the input struct but the output is non-scalar
% (i.e. an array), then that scalar is tiled to the size of the expected
% output array.
% 2014-01-23 Add functionality to specify length of time series in output
% file (Karen Amoudry).
% 2014-01-31 Add functionality to replace user-specified variables in the
% output file with user-specified constant values (Karen Amoudry).
% 2014-02-04 Incorporate Karen's functionality (see above) but with the
% ability to retain the existing behaviour (where a new start time is
% still optional). User-specified constants are also supported but
% instead of specifying a new input argument, if a single scalar value is
% given in the input struct but the output is non-scalar (i.e. an array),
% then that scalar is tiled to the size of the expected output array.
% 2014-02-13 Fix output when only a single time step is in the input
% template restart file (fv_restart). Single time step files have the
% indata put in the time step (previously the ramping from model to
......@@ -80,16 +83,11 @@ function write_FVCOM_restart(fv_restart, out_restart, indata, varargin)
% 2014-07-08 Fix the creation of the Itime2 variable (milliseconds since
% midnight). Also update the help to better describe how the optional
% 'out_date' argument works.
%
% KJA Revision history:
% 2014-01-23 Add functionality to specify length of time series in output
% file.
% 2014-01-31 Add functionality to replace user-specified variables in the
% output file with user-specified constant values.
% 2017-08-01 Add netCDF4 compression to save on file size.
%
%==========================================================================
subname = 'write_FVCOM_restart';
[~, subname] = fileparts(mfilename('fullpath'));
global ftbverbose
if ftbverbose
......@@ -140,8 +138,11 @@ end
fnames = fieldnames(indata);
nf = length(fnames);
if exist(out_restart, 'file')
delete(out_restart)
end
nc = netcdf.open(fv_restart, 'NOWRITE');
ncout = netcdf.create(out_restart, 'clobber');
ncout = netcdf.create(out_restart, 'NETCDF4');
[numdims, numvars, numglobatts, unlimdimID] = netcdf.inq(nc);
......@@ -182,6 +183,8 @@ for ii = 1:numvars
% Create the variables.
varid = netcdf.defVar(ncout, varname, xtype, varDimIDs);
% Enable compression.
netcdf.defVarDeflate(ncout, varid, true, true, 7);
% Get each attribute and add it to the current variable.
for j = 1:varAtts
......
......@@ -81,9 +81,9 @@ netcdf.putAtt(nc, tau_cd_varid, 'units', 'N/m^2');
extra_fields = setdiff(fieldnames(sediment), {'tau_ce', 'tau_cd', 'erate'});
for f = 1:length(extra_fields)
sediment_name = strrep(extra_fields{f}, '_bedfrac', '');
sediment_name = extra_fields{f};
if ftbverbose
fprintf('Adding extra variable: %s\n', sediment_name)
fprintf('Defining extra variable: %s\n', sediment_name)
end
eval(sprintf('%s_varid = netcdf.defVar(nc, ''%s'', ''NC_DOUBLE'', node_dimid);\n', sediment_name, sediment_name))
eval(sprintf('netcdf.putAtt(nc, %s_varid, ''long_name'', ''Fraction of %s in surface layer'');\n', sediment_name, sediment_name))
......@@ -99,6 +99,15 @@ netcdf.putVar(nc, erate_varid, sediment.erate);
netcdf.putVar(nc, tau_ce_varid, sediment.tau_ce);
netcdf.putVar(nc, tau_cd_varid, sediment.tau_cd);
% Add the extra data.
for f = 1:length(extra_fields)
sediment_name = extra_fields{f};
if ftbverbose
fprintf('Adding extra variable: %s\n', sediment_name)
end
eval(sprintf('netcdf.putVar(nc, %s_varid, sediment.(sediment_name));', sediment_name))
end
% close netCDF
netcdf.close(nc)
......
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