Due to a shift in policy, from 0900 GMT on Wednesday 14th July 2021, we will be disabling ssh access to the server for external users. External users who wish to continue to access code repositories on the server will need to switch to using https. This can be accomplished in the following way: 1) On the repo on gitlab, use the clone dialogue and select ‘Clone with HTTPS’ to get the address of the repo; 2) From within the checkout of your repo run: $ git remote set-url origin HTTPS_ADDRESS. Here, replace HTTPS_ADDRESS with the address you have just copied from GitLab. Pulls and pushes will now require you to enter a username and password rather than using a ssh key. If you would prefer not to enter a password each time, you might consider caching your login credentials.

write_FVCOM_forcing.m 21.1 KB
Newer Older
1
function write_FVCOM_forcing(Mobj, fileprefix, data, infos, fver)
2
% Write data out to FVCOM NetCDF forcing file.
3
%
4
% write_FVCOM_forcing(fvcom_forcing_file, data, infos, fver)
5
%
6 7 8
% DESCRIPTION:
%   Takes the given interpolated data (e.g. from grid2fvcom) and writes out
%   to a NetCDF file.
9
%
10 11 12 13 14 15 16 17 18 19
% INPUT:
%   Mobj - MATLAB mesh object
%   fileprefix - Output NetCDF file prefix (plus path) will be
%       fileprefix_{wnd,hfx,evap}.nc if fver is '3.1.0', otherwise output
%       files will be fileprefix_wnd.nc.
%   data - Struct of the data to be written out.
%   infos - Additional remarks to be written to the "infos" NetCDF variable
%   fver - Output for version 3.1.0 or 3.1.6. The latter means all the
%       forcing can go in a single file, the former needs separate files
%       for specific forcing data (wind, heating and precipitation).
20
%
21
% The fields in data may be called any of:
22
%     - 'u10', 'v10', 'uwnd', 'vwnd' - wind components
23
%     - 'slp'       - sea level pressure
24
%     - 'pevpr'     - evaporation
25
%     - 'prate'     - precipitation
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
%     - 'slp'       - sea level pressure
%     - 'nlwrs'     - net longwave radiation*,**
%     - 'nswrs'     - net shortwave radiation*,**
%     - 'shtfl'     - sensible heat net flux*,**
%     - 'lhtfl'     - latent heat net flux*,**
%     - 'lon'       - longitude (vector)
%     - 'lat'       - latitude (vector)
%     - 'x'         - eastings (vector)
%     - 'y'         - northings (vector)
%     - 'nshf'      - pre-computed net surface heat flux**
%
% Fields marked with an * are combined to form the "surface net heat flux"
% (nshf) as follows:
%
%   nshf = -nlwrs + -nswrs - lhtfl - shtfl;
%
% ** Alternatively, a new field 'nshf' (net surface heat flux) can be
% supplied, in which case shtfl and lhtfl are not necessary as their only
% function is in the calculation of net surface heat flux. This approach
% eliminates the need to interpolate so many variables onto the FVCOM grid,
% thus decreasing the time needed to generate the FVCOM input files.
47
%
48 49
% OUTPUT:
%   FVCOM wind speed forcing NetCDF file(s)
50
%
51 52
% Author(s):
%   Pierre Cazenave (Plymouth Marine Laboratory)
53 54 55
%   Karen Thurston (National Oceanography Centre, Liverpool)
%
% PWC Revision history:
56 57 58
%   2012-11-01 - First version based on the parts of grid2fvcom_U10V10.m
%   which dealt with writing the NetCDF file. This version now dynamically
%   deals with varying numbers of forcing data.
59 60 61 62 63
%   2012-11-09 - Add the correct calculation for the net surface heat flux.
%   2013-02-14 - Fix the surface net heat flux sign convention, include the
%   longwave radiation variable and add support for a new field in the
%   input struct ('nshf') which contains the pre-computed net surface heat
%   flux.
64 65 66 67
%
% KJT Revision history:
%   2013-01-16 - Added support for output of sea level pressure.
%
68 69 70
%==========================================================================

multi_out = false; % default to 3.1.6, single output file
71
if nargin < 4 || nargin > 5
72
    error('Incorrect number of arguments')
73
elseif nargin == 5
74 75 76 77 78 79 80 81 82 83 84 85 86 87
    if strcmpi(fver, '3.1.0')
        multi_out = true;
    end
end

subname = 'write_FVCOM_forcing';

global ftbverbose;
if(ftbverbose)
  fprintf('\n')
  fprintf(['begin : ' subname '\n'])
end

tri = Mobj.tri;
88
nNodes = Mobj.nVerts;
89 90 91 92 93 94 95 96 97 98
nElems = Mobj.nElems;
ntimes = numel(data.time);

if strcmpi(Mobj.nativeCoords,'cartesian')
    x = Mobj.x;
    y = Mobj.y;
else
    x = Mobj.lon;
    y = Mobj.lat;
end
99 100
% Create a string for each variable's coordinate attribute
coordString = sprintf('FVCOM %s coordinates', Mobj.nativeCoords);
101

102 103 104 105
% Create element vertices positions
xc = nodes2elems(x, Mobj);
yc = nodes2elems(y, Mobj);

106 107 108 109 110
%--------------------------------------------------------------------------
% Create the NetCDF header for the FVCOM forcing file
%--------------------------------------------------------------------------

if multi_out
111
    suffixes = {'_wnd', '_hfx', '_evap', '_air_press'};
112 113 114 115 116 117 118
else
    suffixes = {'_wnd'};
end

for i=1:length(suffixes)
    nc = netcdf.create([fileprefix, suffixes{i}, '.nc'], 'clobber');

119 120 121
%     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'type','FVCOM Forcing File')
    netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'title','FVCOM Forcing File')
    netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'institution','Plymouth Marine Laboratory')
122
    netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'source','FVCOM grid (unstructured) surface forcing')
123
    netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'history',['File created on ', datestr(now, 'yyyy-mm-dd HH:MM:SS'), ' with write_FVCOM_forcing.m from the fvcom-toolbox (https://github.com/pwcazenave/fvcom-toolbox)'])
124
    netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'references','http://fvcom.smast.umassd.edu, http://codfish.smast.umassd.edu')
125 126
    netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.0')
%     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'infos',infos)
127
    netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'CoordinateSystem',Mobj.nativeCoords)
128
    netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'CoordinateProjection','init=epsg:4326') % WGS84?
129 130 131

    % Dimensions
    nele_dimid=netcdf.defDim(nc,'nele',nElems);
132
    node_dimid=netcdf.defDim(nc,'node',nNodes);
133
    three_dimid=netcdf.defDim(nc,'three',3);
134
    time_dimid=netcdf.defDim(nc,'time',netcdf.getConstant('NC_UNLIMITED'));
135
    datestrlen_dimid=netcdf.defDim(nc,'DateStrLen',26);
136 137 138 139

    % Space variables
    x_varid=netcdf.defVar(nc,'x','NC_FLOAT',node_dimid);
    netcdf.putAtt(nc,x_varid,'long_name','nodal x-coordinate');
140
    netcdf.putAtt(nc,x_varid,'units','meters');
141 142 143

    y_varid=netcdf.defVar(nc,'y','NC_FLOAT',node_dimid);
    netcdf.putAtt(nc,y_varid,'long_name','nodal y-coordinate');
144
    netcdf.putAtt(nc,y_varid,'units','meters');
145

146
    xc_varid=netcdf.defVar(nc,'xc','NC_FLOAT',nele_dimid);
147
    netcdf.putAtt(nc,xc_varid,'long_name','zonal x-coordinate');
148
    netcdf.putAtt(nc,xc_varid,'units','meters');
149

150
    yc_varid=netcdf.defVar(nc,'yc','NC_FLOAT',nele_dimid);
151
    netcdf.putAtt(nc,yc_varid,'long_name','zonal y-coordinate');
152
    netcdf.putAtt(nc,yc_varid,'units','meters');
153

154
    nv_varid=netcdf.defVar(nc,'nv','NC_INT',[nele_dimid, three_dimid]);
155 156
    netcdf.putAtt(nc,nv_varid,'long_name','nodes surrounding element');

157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
    % 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,'format','modified julian day (MJD)');
    netcdf.putAtt(nc,time_varid,'time_zone','UTC');

    itime_varid=netcdf.defVar(nc,'Itime','NC_INT',time_dimid);
    netcdf.putAtt(nc,itime_varid,'units','days since 1858-11-17 00:00:00');
    netcdf.putAtt(nc,itime_varid,'format','modified julian day (MJD)');
    netcdf.putAtt(nc,itime_varid,'time_zone','UTC');

    itime2_varid=netcdf.defVar(nc,'Itime2','NC_INT',time_dimid);
    netcdf.putAtt(nc,itime2_varid,'units','msec since 00:00:00');
    netcdf.putAtt(nc,itime2_varid,'time_zone','UTC');

173
    % Since we have a dynamic number of variables in the struct, try to be a
174
    % bit clever about how to create the output variables.
175 176 177 178 179 180 181 182 183 184 185 186 187 188
    fnames = fieldnames(data);
    used_varids = cell(0);
    used_fnames = cell(0);
    used_dims = cell(0); % exclude time (assume all variables vary in time)

    for vv=1:length(fnames)
        % Need to check both whether we have the data but also whether
        % we're outputting to several NetCDF files. If so, we drop some
        % variables if we're in the wrong file loop.
        switch fnames{vv}
            case {'uwnd', 'u10'}
                if strcmpi(suffixes{i}, '_wnd') || ~multi_out
                    % wind components (assume we have v if we have u)

189 190 191
                    % On the elements
                    u10_varid=netcdf.defVar(nc,'U10','NC_FLOAT',[nele_dimid, time_dimid]);
                    netcdf.putAtt(nc,u10_varid,'long_name','Eastward Wind Speed');
192
%                     netcdf.putAtt(nc,u10_varid,'standard_name','Eastward Wind Speed');
193 194
                    netcdf.putAtt(nc,u10_varid,'units','m/s');
                    netcdf.putAtt(nc,u10_varid,'grid','fvcom_grid');
195
                    netcdf.putAtt(nc,u10_varid,'coordinates',coordString);
196
                    netcdf.putAtt(nc,u10_varid,'type','data');
197 198 199

                    v10_varid=netcdf.defVar(nc,'V10','NC_FLOAT',[nele_dimid, time_dimid]);
                    netcdf.putAtt(nc,v10_varid,'long_name','Northward Wind Speed');
200
%                     netcdf.putAtt(nc,v10_varid,'standard_name','Northward Wind Speed');
201 202
                    netcdf.putAtt(nc,v10_varid,'units','m/s');
                    netcdf.putAtt(nc,v10_varid,'grid','fvcom_grid');
203
                    netcdf.putAtt(nc,v10_varid,'coordinates',coordString);
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
                    netcdf.putAtt(nc,v10_varid,'type','data');

                    uwind_varid=netcdf.defVar(nc,'uwind_speed','NC_FLOAT',[nele_dimid, time_dimid]);
                    netcdf.putAtt(nc,uwind_varid,'long_name','Eastward Wind Speed');
                    netcdf.putAtt(nc,uwind_varid,'standard_name','Wind Speed');
                    netcdf.putAtt(nc,uwind_varid,'units','m/s');
                    netcdf.putAtt(nc,uwind_varid,'grid','fvcom_grid');
                    netcdf.putAtt(nc,uwind_varid,'type','data');

                    vwind_varid=netcdf.defVar(nc,'vwind_speed','NC_FLOAT',[nele_dimid, time_dimid]);
                    netcdf.putAtt(nc,vwind_varid,'long_name','Northward Wind Speed');
                    netcdf.putAtt(nc,vwind_varid,'standard_name','Wind Speed');
                    netcdf.putAtt(nc,vwind_varid,'units','m/s');
                    netcdf.putAtt(nc,vwind_varid,'grid','fvcom_grid');
                    netcdf.putAtt(nc,vwind_varid,'type','data');
219 220 221 222 223 224 225 226

                    % On the nodes
%                     u10_node_varid=netcdf.defVar(nc,'U10','NC_FLOAT',[node_dimid, time_dimid]);
%                     netcdf.putAtt(nc,u10_node_varid,'long_name','Eastward 10-m Velocity');
%                     netcdf.putAtt(nc,u10_node_varid,'standard_name','Eastward Wind Speed');
%                     netcdf.putAtt(nc,u10_node_varid,'units','m/s');
%                     netcdf.putAtt(nc,u10_node_varid,'grid','fvcom_grid');
%                     netcdf.putAtt(nc,u10_node_varid,'type','data');
227
%                     netcdf.putAtt(nc,u10_node_varid,'coordinates',coordString);
228
%
229 230 231 232 233 234
%                     v10_node_varid=netcdf.defVar(nc,'V10','NC_FLOAT',[node_dimid, time_dimid]);
%                     netcdf.putAtt(nc,v10_node_varid,'long_name','Northward 10-m Velocity');
%                     netcdf.putAtt(nc,v10_node_varid,'standard_name','Northward Wind Speed');
%                     netcdf.putAtt(nc,v10_node_varid,'units','m/s');
%                     netcdf.putAtt(nc,v10_node_varid,'grid','fvcom_grid');
%                     netcdf.putAtt(nc,v10_node_varid,'type','data');
235
%                     netcdf.putAtt(nc,v10_node_varid,'coordinates',coordString);
236 237

%                     % Both node and element centred
238 239
%                     used_varids = [used_varids, {'u10_varid', 'v10_varid', 'u10_node_varid', 'v10_node_varid'}];
%                     used_fnames = [used_fnames, {'uwnd', 'vwnd', 'uwnd', 'vwnd'}];
240 241 242 243 244 245
%                     used_dims = [used_dims, {'nElems', 'nElems', 'nNodes', 'nNodes'}];
%                     % Only on the nodes
%                     used_varids = [used_varids, {'u10_node_varid', 'v10_node_varid'}];
%                     used_fnames = [used_fnames, {'uwnd', 'vwnd'}];
%                     used_dims = [used_dims, {'nNodes', 'nNodes'}];
                    % Only on the elements
246 247 248 249 250 251 252 253
%                     used_varids = [used_varids, {'u10_varid', 'v10_varid'}];
%                     used_fnames = [used_fnames, {'uwnd', 'vwnd'}];
%                     used_dims = [used_dims, {'nElems', 'nElems'}];
                    % Only on the elements (both U10/V10 and uwind_speed and
                    % vwind_speed).
                    used_varids = [used_varids, {'u10_varid', 'v10_varid', 'uwind_varid', 'vwind_varid'}];
                    used_fnames = [used_fnames, {'uwnd', 'vwnd', 'uwnd', 'vwnd'}];
                    used_dims = [used_dims, {'nElems', 'nElems', 'nElems', 'nElems'}];
254 255
                end

256 257 258 259 260 261 262
            case 'slp'
                if strcmpi(suffixes{i}, '_air_press') || ~multi_out
                    % Sea level pressure
                    slp_varid=netcdf.defVar(nc,'air_pressure','NC_FLOAT',[node_dimid, time_dimid]);
                    netcdf.putAtt(nc,slp_varid,'long_name','Surface air pressure');
                    netcdf.putAtt(nc,slp_varid,'units','Pa');
                    netcdf.putAtt(nc,slp_varid,'grid','fvcom_grid');
263
                    netcdf.putAtt(nc,slp_varid,'coordinates',coordString);
264 265 266 267 268 269 270
                    netcdf.putAtt(nc,slp_varid,'type','data');

                    used_varids = [used_varids, 'slp_varid'];
                    used_fnames = [used_fnames, fnames{vv}];
                    used_dims = [used_dims, 'nNodes'];
                end

271 272 273
            case 'P_E'
                if strcmpi(suffixes{i}, '_evap') || ~multi_out
                    % Evaporation
274
                    pe_varid=netcdf.defVar(nc,'evap','NC_FLOAT',[node_dimid, time_dimid]);
275
                    netcdf.putAtt(nc,pe_varid,'long_name','Evaporation');
276
                    netcdf.putAtt(nc,pe_varid,'description','Evaporation, ocean lose water is negative');
277 278
                    netcdf.putAtt(nc,pe_varid,'units','m s-1');
                    netcdf.putAtt(nc,pe_varid,'grid','fvcom_grid');
279
                    netcdf.putAtt(nc,pe_varid,'coordinates',coordString);
280 281 282 283
                    netcdf.putAtt(nc,pe_varid,'type','data');

                    used_varids = [used_varids, 'pe_varid'];
                    used_fnames = [used_fnames, fnames{vv}];
284
                    used_dims = [used_dims, 'nNodes'];
285
                end
286

287 288 289
            case 'prate'
                if strcmpi(suffixes{i}, '_evap') || ~multi_out
                    % Precipitation
290
                    prate_varid=netcdf.defVar(nc,'precip','NC_FLOAT',[node_dimid, time_dimid]);
291
                    netcdf.putAtt(nc,prate_varid,'long_name','Precipitation');
292
                    netcdf.putAtt(nc,prate_varid,'description','Precipitation, ocean lose water is negative');
293 294
                    netcdf.putAtt(nc,prate_varid,'units','m s-1');
                    netcdf.putAtt(nc,prate_varid,'grid','fvcom_grid');
295
                    netcdf.putAtt(nc,prate_varid,'coordinates',coordString);
296 297 298 299
                    netcdf.putAtt(nc,prate_varid,'type','data');

                    used_varids = [used_varids, 'prate_varid'];
                    used_fnames = [used_fnames, fnames{vv}];
300
                    used_dims = [used_dims, 'nNodes'];
301 302 303 304 305
                end

            case 'nswrs'
                if strcmpi(suffixes{i}, '_hfx') || ~multi_out
                    % Shortwave radiation
306
                    nswrs_varid=netcdf.defVar(nc,'short_wave','NC_FLOAT',[node_dimid, time_dimid]);
307 308 309
                    netcdf.putAtt(nc,nswrs_varid,'long_name','Short Wave Radiation');
                    netcdf.putAtt(nc,nswrs_varid,'units','W m-2');
                    netcdf.putAtt(nc,nswrs_varid,'grid','fvcom_grid');
310
                    netcdf.putAtt(nc,nswrs_varid,'coordinates',coordString);
311 312 313 314
                    netcdf.putAtt(nc,nswrs_varid,'type','data');

                    used_varids = [used_varids, 'nswrs_varid'];
                    used_fnames = [used_fnames, fnames{vv}];
315
                    used_dims = [used_dims, 'nNodes'];
316 317
                end

318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343
            case 'nlwrs'
                if strcmpi(suffixes{i}, '_hfx') || ~multi_out
                    % Longwave radiation
                    nswrs_varid=netcdf.defVar(nc,'long_wave','NC_FLOAT',[node_dimid, time_dimid]);
                    netcdf.putAtt(nc,nswrs_varid,'long_name','Long Wave Radiation');
                    netcdf.putAtt(nc,nswrs_varid,'units','W m-2');
                    netcdf.putAtt(nc,nswrs_varid,'grid','fvcom_grid');
                    netcdf.putAtt(nc,nswrs_varid,'coordinates',coordString);
                    netcdf.putAtt(nc,nswrs_varid,'type','data');

                    used_varids = [used_varids, 'nlwrs_varid'];
                    used_fnames = [used_fnames, fnames{vv}];
                    used_dims = [used_dims, 'nNodes'];
                end

            case {'shtfl', 'lhtfl', 'nshf'} % , 'nlwrs', 'nswrs'}
                % We can't trigger on nlwrs and nswrs here because they're
                % the triggers for the net longwave and shortwave variables
                % above. Instead, we need to use the latent heat flux
                % variables as the triggers for the net heat flux.
                % We also need to check for the existence of the "net
                % surface heat flux ('nshf')" field which can be created before
                % calling grid2fvcom. This approach means there's fewer
                % calls to the (expensive) interpolation as the net surface
                % heat flux is calculated before being interpolated onto
                % the FVCOM grid.
344 345 346
                try
                    % We might have already made this attribute, so fail
                    % elegantly if we do. This is because we need to put
347 348
                    % all four of shtfl, lhtfl, nlwrs and nswrs to make
                    % Surface Net Heat Flux.
349 350
                    if strcmpi(suffixes{i}, '_hfx') || ~multi_out
                        % Surface net heat flux
351
                        nhf_varid=netcdf.defVar(nc,'net_heat_flux','NC_FLOAT',[node_dimid, time_dimid]);
352 353 354
                        netcdf.putAtt(nc,nhf_varid,'long_name','Surface Net Heat Flux');
                        netcdf.putAtt(nc,nhf_varid,'units','W m-2');
                        netcdf.putAtt(nc,nhf_varid,'grid','fvcom_grid');
355
                        netcdf.putAtt(nc,nhf_varid,'coordinates',coordString);
356 357 358 359 360 361
                        netcdf.putAtt(nc,nhf_varid,'type','data');
                    end
                end
                if strcmpi(suffixes{i}, '_hfx') || ~multi_out
                    % We need to save the current variable name even if we've
                    % already made its attribute.
362 363
                    used_varids = [used_varids, 'nhf_varid'];
                    used_fnames = [used_fnames, fnames{vv}];
364
                    used_dims = [used_dims, 'nNodes'];
365
                end
366

367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386
            case {'time', 'lon', 'lat', 'x', 'y'}
                continue

            otherwise
                if(ftbverbose)
                    warning('Unknown or unused input data type: %s', fnames{vv})
                end
        end
    end

    % End definitions
    netcdf.endDef(nc);

    % Put the easy ones in first.
    netcdf.putVar(nc,nv_varid, tri');
    netcdf.putVar(nc,time_varid,0,ntimes,data.time);
    netcdf.putVar(nc,itime_varid,0,ntimes,floor(data.time));
    netcdf.putVar(nc,itime2_varid,0,ntimes,mod(data.time,1)*24*3600*1000);
    netcdf.putVar(nc,x_varid,x);
    netcdf.putVar(nc,y_varid,y);
387 388
    netcdf.putVar(nc,xc_varid,xc);
    netcdf.putVar(nc,yc_varid,yc);
389

390
    % Now do the dynamic ones. Set the heat flux to not done (0) until we
391
    % hit one of the holy quad (shtfl, lhtfl, nlwrs and nswrs).
392
    hf_done = 0;
393
    nshf = 0;
394
    for ff=1:length(used_fnames)
395
        if strcmpi(used_fnames{ff}, 'shtfl') || strcmpi(used_fnames{ff}, 'lhtfl') || strcmpi(used_fnames{ff}, 'nlwrs') || strcmpi(used_fnames{ff}, 'nswrs')
396
            hf_done = hf_done + 1;
397 398 399 400 401 402 403 404 405 406
            if hf_done == 4 && nshf == 0
                % We've got all four heat parameters, so dump them into the
                % file. We have to flip the signs of the net fluxes and
                % subtract the latent and sensible heat fluxes because
                % NCEP's convention is positive upwards (out of the ocean)
                % and negative downwards (into the ocean), whereas FVCOM is
                % positive downwards (into the ocean) and negative upwards
                % (out of the ocean).
                hf = -data.nlwrs.node + -data.nswrs.node - data.lhtfl.node - data.shtfl.node;
                %hf = data.shtfl.node + data.lhtfl.node + data.nlwrs.node + data.nswrs.node;
407
                netcdf.putVar(nc,nhf_varid,[0,0],[nNodes,ntimes],hf)
408
            end
409 410 411 412 413 414 415 416
        elseif strcmpi(used_fnames{ff}, 'nshf')
            % We have pre-computed net surface heat flux, in which case set
            % hf_done to 4 and put the data into the netCDF. Also set the
            % nshf variable 1 to stop the net surface heat flux variable
            % being overwritten above.
            hf_done = 4;
            nshf = 1;
            netcdf.putVar(nc, nhf_varid, [0, 0], [nNodes, ntimes], data.nshf.node)
417
        else
418 419
            % One of the other data sets for which we can simply dump the
            % existing array without waiting for other data
420
        	if strcmpi(used_dims{ff}, 'nNodes')
421
                eval(['netcdf.putVar(nc,',used_varids{ff},',[0,0],[',used_dims{ff},',ntimes],data.',used_fnames{ff},'.node);'])
422
            else
423
                eval(['netcdf.putVar(nc,',used_varids{ff},',[0,0],[',used_dims{ff},',ntimes],data.',used_fnames{ff},'.data);'])
424
            end
425 426
        end
    end
427 428
    if hf_done ~= 4
        warning('Did not have all the required heat flux parameters. Need ''shtfl'', ''lhtfl'', ''nlwrs'' and ''nwsrs''.')
429
    end
430

431
    % Close the NetCDF file(s)
432
    netcdf.close(nc);
433
end
434 435 436 437

if ftbverbose
    fprintf(['end   : ' subname '\n'])
end