write_FVCOM_forcing

 Write data out to FVCOM NetCDF forcing file.

 write_FVCOM_forcing(fvcom_forcing_file, data, infos, fver)

 DESCRIPTION:
   Takes the given interpolated data (e.g. from grid2fvcom) and writes out
   to a NetCDF file.

 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).

 The fields in data may be called any of:
     - 'u10', 'v10', 'uwnd', 'vwnd' - wind components
     - 'slp'       - sea level pressure
     - 'P_E'       - evaporation
     - 'prate'     - precipitation
     - 'nswrs'     - short wave radiation
     - 'shtfl'\
     - 'lhtfl' >   - combine to form "surface net heat flux"
     - 'nlwrs'/
     - 'time'
     - 'lon'
     - 'lat'
     - 'x'
     - 'y'

 OUTPUT:
   FVCOM wind speed forcing NetCDF file(s)


 Author(s):
   Pierre Cazenave (Plymouth Marine Laboratory)
   Karen Thurston (National Oceanography Centre, Liverpool)

 PWC Revision history:
   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.
   2012-11-09 - Add the correct calculation for the surface net heat flux.

 KJT Revision history:
   2013-01-16 - Added support for output of sea level pressure.

==========================================================================

0001 function write_FVCOM_forcing(Mobj, fileprefix, data, infos, fver)
0002 % Write data out to FVCOM NetCDF forcing file.
0003 %
0004 % write_FVCOM_forcing(fvcom_forcing_file, data, infos, fver)
0005 %
0006 % DESCRIPTION:
0007 %   Takes the given interpolated data (e.g. from grid2fvcom) and writes out
0008 %   to a NetCDF file.
0009 %
0010 % INPUT:
0011 %   Mobj - MATLAB mesh object
0012 %   fileprefix - Output NetCDF file prefix (plus path) will be
0013 %       fileprefix_{wnd,hfx,evap}.nc if fver is '3.1.0', otherwise output
0014 %       files will be fileprefix_wnd.nc.
0015 %   data - Struct of the data to be written out.
0016 %   infos - Additional remarks to be written to the "infos" NetCDF variable
0017 %   fver - Output for version 3.1.0 or 3.1.6. The latter means all the
0018 %       forcing can go in a single file, the former needs separate files
0019 %       for specific forcing data (wind, heating and precipitation).
0020 %
0021 % The fields in data may be called any of:
0022 %     - 'u10', 'v10', 'uwnd', 'vwnd' - wind components
0023 %     - 'slp'       - sea level pressure
0024 %     - 'P_E'       - evaporation
0025 %     - 'prate'     - precipitation
0026 %     - 'nswrs'     - short wave radiation
0027 %     - 'shtfl'\
0028 %     - 'lhtfl' >   - combine to form "surface net heat flux"
0029 %     - 'nlwrs'/
0030 %     - 'time'
0031 %     - 'lon'
0032 %     - 'lat'
0033 %     - 'x'
0034 %     - 'y'
0035 %
0036 % OUTPUT:
0037 %   FVCOM wind speed forcing NetCDF file(s)
0038 %
0039 %
0040 % Author(s):
0041 %   Pierre Cazenave (Plymouth Marine Laboratory)
0042 %   Karen Thurston (National Oceanography Centre, Liverpool)
0043 %
0044 % PWC Revision history:
0045 %   2012-11-01 - First version based on the parts of grid2fvcom_U10V10.m
0046 %   which dealt with writing the NetCDF file. This version now dynamically
0047 %   deals with varying numbers of forcing data.
0048 %   2012-11-09 - Add the correct calculation for the surface net heat flux.
0049 %
0050 % KJT Revision history:
0051 %   2013-01-16 - Added support for output of sea level pressure.
0052 %
0053 %==========================================================================
0054 
0055 multi_out = false; % default to 3.1.6, single output file
0056 if nargin < 4 || nargin > 5
0057     error('Incorrect number of arguments')
0058 elseif nargin == 5
0059     if strcmpi(fver, '3.1.0')
0060         multi_out = true;
0061     end
0062 end
0063 
0064 subname = 'write_FVCOM_forcing';
0065 
0066 global ftbverbose;
0067 if(ftbverbose)
0068   fprintf('\n')
0069   fprintf(['begin : ' subname '\n'])
0070 end
0071 
0072 tri = Mobj.tri;
0073 nNodes = Mobj.nVerts;
0074 nElems = Mobj.nElems;
0075 ntimes = numel(data.time);
0076 
0077 if strcmpi(Mobj.nativeCoords,'cartesian')
0078     x = Mobj.x;
0079     y = Mobj.y;
0080 else
0081     x = Mobj.lon;
0082     y = Mobj.lat;
0083 end
0084 
0085 % Create element vertices positions
0086 xc = nodes2elems(x, Mobj);
0087 yc = nodes2elems(y, Mobj);
0088 
0089 %--------------------------------------------------------------------------
0090 % Create the NetCDF header for the FVCOM forcing file
0091 %--------------------------------------------------------------------------
0092 
0093 if multi_out
0094     suffixes = {'_wnd', '_hfx', '_evap', '_air_press'};
0095 else
0096     suffixes = {'_wnd'};
0097 end
0098 
0099 for i=1:length(suffixes)
0100     nc = netcdf.create([fileprefix, suffixes{i}, '.nc'], 'clobber');
0101 
0102 %     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'type','FVCOM Forcing File')
0103     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'title','FVCOM Forcing File')
0104     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'institution','Plymouth Marine Laboratory')
0105     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'source','FVCOM grid (unstructured) surface forcing')
0106     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)'])
0107     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'references','http://fvcom.smast.umassd.edu, http://codfish.smast.umassd.edu')
0108     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.0')
0109 %     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'infos',infos)
0110     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'CoordinateSystem',Mobj.nativeCoords)
0111     netcdf.putAtt(nc,netcdf.getConstant('NC_GLOBAL'),'CoordinateProjection','init=epsg:4326') % WGS84?
0112 
0113     % Dimensions
0114     nele_dimid=netcdf.defDim(nc,'nele',nElems);
0115     node_dimid=netcdf.defDim(nc,'node',nNodes);
0116     three_dimid=netcdf.defDim(nc,'three',3);
0117     time_dimid=netcdf.defDim(nc,'time',netcdf.getConstant('NC_UNLIMITED'));
0118     datestrlen_dimid=netcdf.defDim(nc,'DateStrLen',26);
0119 
0120     % Space variables
0121     x_varid=netcdf.defVar(nc,'x','NC_FLOAT',node_dimid);
0122     netcdf.putAtt(nc,x_varid,'long_name','nodal x-coordinate');
0123     netcdf.putAtt(nc,x_varid,'units','meters');
0124 
0125     y_varid=netcdf.defVar(nc,'y','NC_FLOAT',node_dimid);
0126     netcdf.putAtt(nc,y_varid,'long_name','nodal y-coordinate');
0127     netcdf.putAtt(nc,y_varid,'units','meters');
0128 
0129     xc_varid=netcdf.defVar(nc,'xc','NC_FLOAT',nele_dimid);
0130     netcdf.putAtt(nc,xc_varid,'long_name','zonal x-coordinate');
0131     netcdf.putAtt(nc,xc_varid,'units','meters');
0132 
0133     yc_varid=netcdf.defVar(nc,'yc','NC_FLOAT',nele_dimid);
0134     netcdf.putAtt(nc,yc_varid,'long_name','zonal y-coordinate');
0135     netcdf.putAtt(nc,yc_varid,'units','meters');
0136 
0137     nv_varid=netcdf.defVar(nc,'nv','NC_INT',[nele_dimid, three_dimid]);
0138     netcdf.putAtt(nc,nv_varid,'long_name','nodes surrounding element');
0139 
0140     % Time variables
0141     time_varid=netcdf.defVar(nc,'time','NC_FLOAT',time_dimid);
0142     netcdf.putAtt(nc,time_varid,'long_name','time');
0143     netcdf.putAtt(nc,time_varid,'units','days since 1858-11-17 00:00:00');
0144     netcdf.putAtt(nc,time_varid,'format','modified julian day (MJD)');
0145     netcdf.putAtt(nc,time_varid,'time_zone','UTC');
0146 
0147     itime_varid=netcdf.defVar(nc,'Itime','NC_INT',time_dimid);
0148     netcdf.putAtt(nc,itime_varid,'units','days since 1858-11-17 00:00:00');
0149     netcdf.putAtt(nc,itime_varid,'format','modified julian day (MJD)');
0150     netcdf.putAtt(nc,itime_varid,'time_zone','UTC');
0151 
0152     itime2_varid=netcdf.defVar(nc,'Itime2','NC_INT',time_dimid);
0153     netcdf.putAtt(nc,itime2_varid,'units','msec since 00:00:00');
0154     netcdf.putAtt(nc,itime2_varid,'time_zone','UTC');
0155 
0156     % Since we have a dynamic number of variables in the struct, try to be a
0157     % bit clever about how to create the output variables.
0158     fnames = fieldnames(data);
0159     used_varids = cell(0);
0160     used_fnames = cell(0);
0161     used_dims = cell(0); % exclude time (assume all variables vary in time)
0162 
0163     for vv=1:length(fnames)
0164         % Need to check both whether we have the data but also whether
0165         % we're outputting to several NetCDF files. If so, we drop some
0166         % variables if we're in the wrong file loop.
0167         switch fnames{vv}
0168             case {'uwnd', 'u10'}
0169                 if strcmpi(suffixes{i}, '_wnd') || ~multi_out
0170                     % wind components (assume we have v if we have u)
0171 
0172                     % On the elements
0173                     u10_varid=netcdf.defVar(nc,'U10','NC_FLOAT',[nele_dimid, time_dimid]);
0174                     netcdf.putAtt(nc,u10_varid,'long_name','Eastward Wind Speed');
0175 %                     netcdf.putAtt(nc,u10_varid,'standard_name','Eastward Wind Speed');
0176                     netcdf.putAtt(nc,u10_varid,'units','m/s');
0177                     netcdf.putAtt(nc,u10_varid,'grid','fvcom_grid');
0178                     netcdf.putAtt(nc,u10_varid,'coordinates','');
0179                     netcdf.putAtt(nc,u10_varid,'type','data');
0180 
0181                     v10_varid=netcdf.defVar(nc,'V10','NC_FLOAT',[nele_dimid, time_dimid]);
0182                     netcdf.putAtt(nc,v10_varid,'long_name','Northward Wind Speed');
0183 %                     netcdf.putAtt(nc,v10_varid,'standard_name','Northward Wind Speed');
0184                     netcdf.putAtt(nc,v10_varid,'units','m/s');
0185                     netcdf.putAtt(nc,v10_varid,'grid','fvcom_grid');
0186                     netcdf.putAtt(nc,v10_varid,'coordinates','');
0187                     netcdf.putAtt(nc,v10_varid,'type','data');
0188 
0189                     uwind_varid=netcdf.defVar(nc,'uwind_speed','NC_FLOAT',[nele_dimid, time_dimid]);
0190                     netcdf.putAtt(nc,uwind_varid,'long_name','Eastward Wind Speed');
0191                     netcdf.putAtt(nc,uwind_varid,'standard_name','Wind Speed');
0192                     netcdf.putAtt(nc,uwind_varid,'units','m/s');
0193                     netcdf.putAtt(nc,uwind_varid,'grid','fvcom_grid');
0194                     netcdf.putAtt(nc,uwind_varid,'type','data');
0195 
0196                     vwind_varid=netcdf.defVar(nc,'vwind_speed','NC_FLOAT',[nele_dimid, time_dimid]);
0197                     netcdf.putAtt(nc,vwind_varid,'long_name','Northward Wind Speed');
0198                     netcdf.putAtt(nc,vwind_varid,'standard_name','Wind Speed');
0199                     netcdf.putAtt(nc,vwind_varid,'units','m/s');
0200                     netcdf.putAtt(nc,vwind_varid,'grid','fvcom_grid');
0201                     netcdf.putAtt(nc,vwind_varid,'type','data');
0202 
0203                     % On the nodes
0204 %                     u10_node_varid=netcdf.defVar(nc,'U10','NC_FLOAT',[node_dimid, time_dimid]);
0205 %                     netcdf.putAtt(nc,u10_node_varid,'long_name','Eastward 10-m Velocity');
0206 %                     netcdf.putAtt(nc,u10_node_varid,'standard_name','Eastward Wind Speed');
0207 %                     netcdf.putAtt(nc,u10_node_varid,'units','m/s');
0208 %                     netcdf.putAtt(nc,u10_node_varid,'grid','fvcom_grid');
0209 %                     netcdf.putAtt(nc,u10_node_varid,'type','data');
0210 %                     netcdf.putAtt(nc,u10_node_varid,'coordinates','');
0211 %
0212 %                     v10_node_varid=netcdf.defVar(nc,'V10','NC_FLOAT',[node_dimid, time_dimid]);
0213 %                     netcdf.putAtt(nc,v10_node_varid,'long_name','Northward 10-m Velocity');
0214 %                     netcdf.putAtt(nc,v10_node_varid,'standard_name','Northward Wind Speed');
0215 %                     netcdf.putAtt(nc,v10_node_varid,'units','m/s');
0216 %                     netcdf.putAtt(nc,v10_node_varid,'grid','fvcom_grid');
0217 %                     netcdf.putAtt(nc,v10_node_varid,'type','data');
0218 %                     netcdf.putAtt(nc,v10_node_varid,'coordinates','');
0219 
0220 %                     % Both node and element centred
0221 %                     used_varids = [used_varids, {'u10_varid', 'v10_varid', 'u10_node_varid', 'v10_node_varid'}];
0222 %                     used_fnames = [used_fnames, {'uwnd', 'vwnd', 'uwnd', 'vwnd'}];
0223 %                     used_dims = [used_dims, {'nElems', 'nElems', 'nNodes', 'nNodes'}];
0224 %                     % Only on the nodes
0225 %                     used_varids = [used_varids, {'u10_node_varid', 'v10_node_varid'}];
0226 %                     used_fnames = [used_fnames, {'uwnd', 'vwnd'}];
0227 %                     used_dims = [used_dims, {'nNodes', 'nNodes'}];
0228                     % Only on the elements
0229 %                     used_varids = [used_varids, {'u10_varid', 'v10_varid'}];
0230 %                     used_fnames = [used_fnames, {'uwnd', 'vwnd'}];
0231 %                     used_dims = [used_dims, {'nElems', 'nElems'}];
0232                     % Only on the elements (both U10/V10 and uwind_speed and
0233                     % vwind_speed).
0234                     used_varids = [used_varids, {'u10_varid', 'v10_varid', 'uwind_varid', 'vwind_varid'}];
0235                     used_fnames = [used_fnames, {'uwnd', 'vwnd', 'uwnd', 'vwnd'}];
0236                     used_dims = [used_dims, {'nElems', 'nElems', 'nElems', 'nElems'}];
0237                 end
0238 
0239             case 'slp'
0240                 if strcmpi(suffixes{i}, '_air_press') || ~multi_out
0241                     % Sea level pressure
0242                     slp_varid=netcdf.defVar(nc,'air_pressure','NC_FLOAT',[node_dimid, time_dimid]);
0243                     netcdf.putAtt(nc,slp_varid,'long_name','Surface air pressure');
0244                     netcdf.putAtt(nc,slp_varid,'units','Pa');
0245                     netcdf.putAtt(nc,slp_varid,'grid','fvcom_grid');
0246                     netcdf.putAtt(nc,slp_varid,'coordinates','');
0247                     netcdf.putAtt(nc,slp_varid,'type','data');
0248 
0249                     used_varids = [used_varids, 'slp_varid'];
0250                     used_fnames = [used_fnames, fnames{vv}];
0251                     used_dims = [used_dims, 'nNodes'];
0252                 end
0253 
0254             case 'P_E'
0255                 if strcmpi(suffixes{i}, '_evap') || ~multi_out
0256                     % Evaporation
0257                     pe_varid=netcdf.defVar(nc,'evap','NC_FLOAT',[node_dimid, time_dimid]);
0258                     netcdf.putAtt(nc,pe_varid,'long_name','Evaporation');
0259                     netcdf.putAtt(nc,pe_varid,'description','Evaporation, ocean lose water is negative');
0260                     netcdf.putAtt(nc,pe_varid,'units','m s-1');
0261                     netcdf.putAtt(nc,pe_varid,'grid','fvcom_grid');
0262                     netcdf.putAtt(nc,pe_varid,'coordinates','');
0263                     netcdf.putAtt(nc,pe_varid,'type','data');
0264 
0265                     used_varids = [used_varids, 'pe_varid'];
0266                     used_fnames = [used_fnames, fnames{vv}];
0267                     used_dims = [used_dims, 'nNodes'];
0268                 end
0269 
0270             case 'prate'
0271                 if strcmpi(suffixes{i}, '_evap') || ~multi_out
0272                     % Precipitation
0273                     prate_varid=netcdf.defVar(nc,'precip','NC_FLOAT',[node_dimid, time_dimid]);
0274                     netcdf.putAtt(nc,prate_varid,'long_name','Precipitation');
0275                     netcdf.putAtt(nc,prate_varid,'description','Precipitation, ocean lose water is negative');
0276                     netcdf.putAtt(nc,prate_varid,'units','m s-1');
0277                     netcdf.putAtt(nc,prate_varid,'grid','fvcom_grid');
0278                     netcdf.putAtt(nc,prate_varid,'coordinates','');
0279                     netcdf.putAtt(nc,prate_varid,'type','data');
0280 
0281                     used_varids = [used_varids, 'prate_varid'];
0282                     used_fnames = [used_fnames, fnames{vv}];
0283                     used_dims = [used_dims, 'nNodes'];
0284                 end
0285 
0286             case 'nswrs'
0287                 if strcmpi(suffixes{i}, '_hfx') || ~multi_out
0288                     % Shortwave radiation
0289                     nswrs_varid=netcdf.defVar(nc,'short_wave','NC_FLOAT',[node_dimid, time_dimid]);
0290                     netcdf.putAtt(nc,nswrs_varid,'long_name','Short Wave Radiation');
0291                     netcdf.putAtt(nc,nswrs_varid,'units','W m-2');
0292                     netcdf.putAtt(nc,nswrs_varid,'grid','fvcom_grid');
0293                     netcdf.putAtt(nc,nswrs_varid,'coordinates','');
0294                     netcdf.putAtt(nc,nswrs_varid,'type','data');
0295 
0296                     used_varids = [used_varids, 'nswrs_varid'];
0297                     used_fnames = [used_fnames, fnames{vv}];
0298                     used_dims = [used_dims, 'nNodes'];
0299                 end
0300 
0301             case {'shtfl', 'lhtfl', 'nlwrs'}
0302                 try
0303                     % We might have already made this attribute, so fail
0304                     % elegantly if we do. This is because we need to put
0305                     % all three of shtfl, lhtfl and nlwrs to make Surface
0306                     % Net Heat Flux.
0307                     if strcmpi(suffixes{i}, '_hfx') || ~multi_out
0308                         % Surface net heat flux
0309                         nhf_varid=netcdf.defVar(nc,'net_heat_flux','NC_FLOAT',[node_dimid, time_dimid]);
0310                         netcdf.putAtt(nc,nhf_varid,'long_name','Surface Net Heat Flux');
0311                         netcdf.putAtt(nc,nhf_varid,'units','W m-2');
0312                         netcdf.putAtt(nc,nhf_varid,'grid','fvcom_grid');
0313                         netcdf.putAtt(nc,nhf_varid,'coordinates','');
0314                         netcdf.putAtt(nc,nhf_varid,'type','data');
0315                     end
0316                 end
0317                 if strcmpi(suffixes{i}, '_hfx') || ~multi_out
0318                     % We need to save the current variable name even if we've
0319                     % already made its attribute.
0320                     used_varids = [used_varids, 'nhf_varid'];
0321                     used_fnames = [used_fnames, fnames{vv}];
0322                     used_dims = [used_dims, 'nNodes'];
0323                 end
0324 
0325             case {'time', 'lon', 'lat', 'x', 'y'}
0326                 continue
0327 
0328             otherwise
0329                 if(ftbverbose)
0330                     warning('Unknown or unused input data type: %s', fnames{vv})
0331                 end
0332         end
0333     end
0334 
0335     % End definitions
0336     netcdf.endDef(nc);
0337 
0338     % Put the easy ones in first.
0339     netcdf.putVar(nc,nv_varid, tri');
0340     netcdf.putVar(nc,time_varid,0,ntimes,data.time);
0341     netcdf.putVar(nc,itime_varid,0,ntimes,floor(data.time));
0342     netcdf.putVar(nc,itime2_varid,0,ntimes,mod(data.time,1)*24*3600*1000);
0343     netcdf.putVar(nc,x_varid,x);
0344     netcdf.putVar(nc,y_varid,y);
0345     netcdf.putVar(nc,xc_varid,xc);
0346     netcdf.putVar(nc,yc_varid,yc);
0347 
0348     % Now do the dynamic ones. Set the heat flux to not done (0) until we
0349     % hit one of the holy trinity (shtfl, lhtfl, nlwrs).
0350     hf_done = 0;
0351     for ff=1:length(used_fnames)
0352         if strcmpi(used_fnames{ff}, 'shtfl') || strcmpi(used_fnames{ff}, 'lhtfl') || strcmpi(used_fnames{ff}, 'nlwrs')
0353             hf_done = hf_done + 1;
0354             if hf_done == 3
0355                 % We've got all three heat parameters, so dump them into the file.
0356                 hf = data.shtfl.node + data.lhtfl.node + data.nlwrs.node;
0357                 netcdf.putVar(nc,nhf_varid,[0,0],[nNodes,ntimes],hf)
0358             end
0359         else
0360             % One of the other data sets for which we can simply dump the
0361             % existing array without waiting for other data
0362             if strcmpi(used_dims{ff}, 'nNodes')
0363                 eval(['netcdf.putVar(nc,',used_varids{ff},',[0,0],[',used_dims{ff},',ntimes],data.',used_fnames{ff},'.node);'])
0364             else
0365                 eval(['netcdf.putVar(nc,',used_varids{ff},',[0,0],[',used_dims{ff},',ntimes],data.',used_fnames{ff},'.data);'])
0366             end
0367         end
0368     end
0369     if hf_done ~= 3
0370         warning('Did not have all the required heat flux parameters. Need ''shtfl'', ''lhtfl'', ''nlwrs''')
0371     end
0372 
0373     % Close the NetCDF file(s)
0374     netcdf.close(nc);
0375 end