function varargout = writeNet(ncfile,X,Y,varargin)
%writeNet   Write network nodes for an unstructured net from a curvilinear mesh
% or unstructured mesh generated by D-Flow FM or RGF-grid (incl. netlinks)
%
%  dflowfm.writeNet(ncfile,X,Y,<Links>,<Z>) 
%
%  If <Links> is provided, the X, Y & Links (arrays) are directly written
%  to a netCDF file with correct CF and DflowFM attributes. If <Links> is
%  not provided, meshgrid is called on x and y (matrices), and 
%  <Links> is calculated (this takes some time and can only be used for
%  curvylinear-type grids). Note that x, y & Links can easily be read with
%  dflowfm.readNet (see Example 2 below). If bathymetry is to be added as 
%  well, also add the <z> array, which can also be read with dflowfm.readNet
%  (again, see Example 2 below). Note that this is not possible without the
%  <Links> array.
%
% This file can be used as input for dflowfm.fillDep.
%
% Example 1:
%
%  ncfile = 'writeNet_test_net.nc';
%  [X,Y]=meshgrid(1:3,1:4);
%  X([4 9])=nan;Y([4 9])=nan;
%  deletefile(ncfile)
%  dflowfm.writeNet(ncfile,X,Y)
%  dflowfm.plotNet (ncfile)
%  D =    nc2struct(ncfile)
%  hold on
%  poly_bi_plot(D.NetLink,D.NetNode_x,D.NetNode_y,'r','linewidth',2);
%
% Example 2:
%
%  delft3d_grid = wlgrid('read','wlgrid.grd');
%  dflowfm.writeNet('fm_grid_from_d3d_grid_net.nc',delft3d_grid.X,delft3d_grid.Y);
%
% Example 3:
%
%  FM_mesh = dflowfm.readNet('D_Flow_FM_net.nc');
%  FM_mesh.node.x = FM_mesh.node.x + 50; (*)
%  FM_mesh.node.z = FM_mesh.node.x - 1; (**)
%  dflowfm.writeNet('New_D_Flow_FM_net.nc',FM_mesh.node.x,FM_mesh.node.y,FM_mesh.edge.NetLink,FM_mesh.node.z);
%
%  (*)  Manual adjustment, e.g. shift of 50 meters, for projected (UTM) zones only!
%  (**) Manual adjustment, e.g. shift of bed = z-level (1 meter lower)
%
% See also: dflowfm, delft3d, fillDep, quat2net, dflowfm.readNet,
% dflowfm.mike2dflowfm

%   --------------------------------------------------------------------
%   Copyright (C) 2012 Deltares
%       Gerben de Boer
%
%       <g.j.deboer@deltares.nl>
%
%       Deltares
%       P.O. Box 177
%       2600 MH Delft
%       The Netherlands
%
%   Revision:   Freek Scheel 2013
%
%     - Compatibility with readNet and writeNet (Addition of netLinks input)
%     - Z-levels can be added as well
%     - When overwriting a file, it is first deleted (prevents errors due
%       to net-cdf standard reading and writing approach when files are
%       different in size)
%     - Some minor tweaks and few issues resolved
%
%   This library is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this library.  If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $

   OPT.institution = '';
   OPT.xname       = 'NetNode_x';
   OPT.yname       = 'NetNode_y';
   OPT.zname       = 'NetNode_z';

   if nargin==0
      varargout = {OPT};
      return
   end
   
   if (nargin)>3
       NetLink = varargin{1};
       varargin{1} = [];
   else
       NetLink     = [];
   end
   
   if (nargin)>4
       Z = varargin{2};
       varargin{2} = [];
   else
       Z     = [];
   end
   
   if (nargin)>5
       NetLinkType = varargin{3};
       varargin{2} = [];
   else
       NetLinkType     = [];
   end
   
   if verLessThan('matlab','7.12.0.635')
      error('At least Matlab release R2011a is required for writing netCDF files due tue NCWRITESCHEMA.')
   end

%% determine links between active nodes

      if isempty(NetLink)
      NetNode_mask = ~isnan(X) & ~isnan(Y);
      nNetNode = sum(NetNode_mask(:));
      if length(size(X)==2)
          [ContourLink,NetLink]=quat2net(X,Y,'sub2ind',0); % subsind is used because we skip nan coordinates here
          vals = find(isnan(X));
          
% ORDERS OF MAGNITUDE TOO SLOW, use sub2ind option in quat2net
%           %% adjust pointer for removal of nans
%           % this step takes much time !!!
%           adjust = 0.*NetLink;
%           for i=1:length(vals)
%              adjust = adjust + (NetLink >=vals(i));
%           end
%           NetLink = NetLink - adjust;
          %%
          
      else
          warning('NetLink calculation only implemented for curvi-linear meshes')
      end
      X = X(NetNode_mask(:));
      Y = Y(NetNode_mask(:));
      else
      nNetNode = length(X);
      NetLink = NetLink';
      end      

%% 1a Create file (add all NEFIS 'map-version' group info)

      %% Add overall meta info
      %  http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#description-of-file-contents

      nc = struct('Name','/','Format','classic'); % hard-coded to allow Dflow-FM library to read it
      nc.Attributes(    1) = struct('Name','title'              ,'Value',  'D-Flow-FM grid file');
      nc.Attributes(end+1) = struct('Name','institution'        ,'Value',  OPT.institution);
      nc.Attributes(end+1) = struct('Name','references'         ,'Value',  'http://svn.oss.deltares.nl');
      nc.Attributes(end+1) = struct('Name','source'             ,'Value',  '$HeadURL$ $Id$');
      nc.Attributes(end+1) = struct('Name','history'            ,'Value', ['Created on ',datestr(now,'yyyy-mm-ddTHH:MM:SS')]);
      nc.Attributes(end+1) = struct('Name','Conventions'        ,'Value',  'CF-1.4:Deltares-0.1');

      nc.Attributes(end+1) = struct('Name','terms_for_use'      ,'Value', ['These data can be used freely for research purposes provided that the following source is acknowledged: ',OPT.institution]);
      nc.Attributes(end+1) = struct('Name','disclaimer'         ,'Value',  'This data is made available in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.');

%% 2 Create dimensions

      nc.Dimensions(1) = struct('Name','nNetNode'       ,'Length',nNetNode);
      nc.Dimensions(2) = struct('Name','nNetLink'       ,'Length',size(NetLink,1));
      nc.Dimensions(3) = struct('Name','nNetLinkPts'    ,'Length',2);
      nc.Dimensions(4) = struct('Name','nBndLink'       ,'Length',1);
      nc.Dimensions(5) = struct('Name','nNetElem'       ,'Length',1);
      nc.Dimensions(6) = struct('Name','nNetElemMaxNode','Length',7);

%% 3 Create variables

      ifld     = 0;
      
      ifld     = ifld + 1;clear attr dims
      attr(    1)  = struct('Name', 'units'         , 'Value', 'm');
      attr(end+1)  = struct('Name', 'standard_name' , 'Value', 'projection_x_coordinate');
      attr(end+1)  = struct('Name', 'long_name'     , 'Value', 'netnodal x-coordinate');
      attr(end+1)  = struct('Name', 'grid_mapping'  , 'Value', 'projected_coordinate_system');
      attr(end+1)  = struct('Name', 'actual_range'  , 'Value', [nan nan]);
      nc.Variables(ifld) = struct('Name'       , OPT.xname, ...
                                  'Datatype'   , 'double', ...
                                  'Dimensions' , nc.Dimensions(1), ...
                                  'Attributes' , attr,...
                                  'FillValue'  , []); % this doesn't do anything

      ifld     = ifld + 1;clear attr dims
      attr(    1)  = struct('Name', 'units'         , 'Value', 'm');
      attr(end+1)  = struct('Name', 'standard_name' , 'Value', 'projection_y_coordinate');
      attr(end+1)  = struct('Name', 'long_name'     , 'Value', 'netnodal x-coordinate');
      attr(end+1)  = struct('Name', 'grid_mapping'  , 'Value', 'projected_coordinate_system');
      attr(end+1)  = struct('Name', 'actual_range'  , 'Value', [nan nan]);
      nc.Variables(ifld) = struct('Name'       , OPT.yname, ...
                                  'Datatype'   , 'double', ...
                                  'Dimensions' , nc.Dimensions(1), ...
                                  'Attributes' , attr,...
                                  'FillValue'  , []); % this doesn't do anything

      ifld     = ifld + 1;clear attr dims
      attr(    1)  = struct('Name', 'units'         , 'Value', 'm');
      attr(end+1)  = struct('Name', 'positive'      , 'Value', 'up');
      attr(end+1)  = struct('Name', 'standard_name' , 'Value', 'sea_floor_depth');
      attr(end+1)  = struct('Name', 'long_name'     , 'Value', 'Bottom level at net nodes (flow element''s corners)');
      attr(end+1)  = struct('Name', 'coordinates'   , 'Value', [OPT.xname ' ' OPT.yname]);
      attr(end+1)  = struct('Name', 'grid_mapping'  , 'Value', 'projected_coordinate_system');
      attr(end+1)  = struct('Name', 'actual_range'  , 'Value', [nan nan]);
warning([mfilename,'agree with developer on _FillValue'])
      attr(end+1)  = struct('Name', '_FillValue'    , 'Value', -999); % this initializes at NaN rather than 9.9692e36
      nc.Variables(ifld) = struct('Name'       , OPT.zname, ...
                                  'Datatype'   , 'double', ...
                                  'Dimensions' , nc.Dimensions(1), ...
                                  'Attributes' , attr,...
                                  'FillValue'  , []); % this doesn't do anything

      ifld     = ifld + 1;clear attr dims
      attr(    1)  = struct('Name', 'standard_name' , 'Value', 'netlink');
      attr(end+1)  = struct('Name', 'long_name'     , 'Value', 'link between two netnodes');
      nc.Variables(ifld) = struct('Name'       , 'NetLink', ...
                                  'Datatype'   , 'int32', ...
                                  'Dimensions' , nc.Dimensions([3 2]), ...
                                  'Attributes' , attr,...
                                  'FillValue'  , []); % this doesn't do anything

      ifld     = ifld + 1;clear attr dims
      attr(    1)  = struct('Name', 'long_name'     , 'Value', 'type of netlink');
      attr(end+1)  = struct('Name', 'valid_range'   , 'Value', int32([0 2]));
      attr(end+1)  = struct('Name', 'flag_values'   , 'Value', int32([0 1 2]));
      attr(end+1)  = struct('Name', 'flag_meanings' , 'Value', 'closed_link_between_2D_nodes link_between_1D_nodes link_between_2D_nodes');
      nc.Variables(ifld) = struct('Name'       , 'NetLinkType', ...
                                  'Datatype'   , 'int32', ...
                                  'Dimensions' , nc.Dimensions([2]), ...
                                  'Attributes' , attr,...
                                  'FillValue'  , []); % this doesn't do anything

%% 4 Create netCDF file

% Check if file already exists, if so, it must be deleted first, due to 
% net-ncf file and save approach

      if exist(ncfile)==2
          warning(' ');
          disp('File already exists, will be overwritten!');
          delete(ncfile);
      end
      
      if length(ncfile)<8
          warning('DfLow-FM grid file should end with "_net.nc"');
      elseif ~strcmpi(ncfile(end-6:end),'_net.nc')
          warning('DfLow-FM grid file should end with "_net.nc"');
      end

      ncwriteschema(ncfile, nc);

%% 5 Fill variables (always)
%    Data is initialized as NaN due to attribute '_FillValue' in ncwriteschema.
      
      ncwrite(ncfile,OPT.xname,X(:));ncwriteatt(ncfile,OPT.xname,'actual_range',[min(X(:)) max(X(:))]);
      ncwrite(ncfile,OPT.yname,Y(:));ncwriteatt(ncfile,OPT.yname,'actual_range',[min(Y(:)) max(Y(:))]);
      if ~isempty(Z)
          ncwrite(ncfile,OPT.zname,Z(:));ncwriteatt(ncfile,OPT.zname,'actual_range',[min(Z(:)) max(Z(:))]);
      end
      ncwrite(ncfile,'NetLink',NetLink');
      
      if ~isempty(NetLinkType)
        ncwrite(ncfile,'NetLinkType',NetLinkType);
      else
        ncwrite(ncfile,'NetLinkType',repmat(int32(2),[1 size(NetLink,1)]));
      end

      disp(['D-Flow FM grid file saved as ' ncfile])
      
% dimensions:
% 	nNetNode = 42984 ;                                                // node.n
% 	nNetLink = 84880 ;
% 	nNetLinkPts = 2 ;
% 	nBndLink = 2208 ;
% 	nNetElem = 41888 ;                                                // face.FlowElemSize
% 	nNetElemMaxNode = 7 ;
% 
% 
% variables:
% 	// Preference 'PRESERVE_FVD':  false,
% 	// dimensions consistent with ncBrowse, not with native MATLAB netcdf package.
% 	double NetNode_x(nNetNode), shape = [42984]
% 		NetNode_x:units = "m" 
% 		NetNode_x:standard_name = "projection_x_coordinate" 
% 		NetNode_x:long_name = "netnodal x-coordinate" 
% 	double NetNode_y(nNetNode), shape = [42984]
% 		NetNode_y:units = "m" 
% 		NetNode_y:standard_name = "projection_y_coordinate" 
% 		NetNode_y:long_name = "netnodal y-coordinate" 
% 	double NetNode_z(nNetNode), shape = [42984]
%
% 	int32 NetLink(nNetLink,nNetLinkPts), shape = [84880 2]
% 		NetLink:standard_name = "netlink" 
% 		NetLink:long_name = "link between two netnodes" 
% 	int32 NetLinkType(nNetLink), shape = [84880]
% 		NetLinkType:long_name = "type of netlink" 
% 		NetLinkType:valid_range = 0 2 d
% 		NetLinkType:flag_values = 0 1 2 d
% 		NetLinkType:flag_meanings = "closed_link_between_2D_nodes link_between_1D_nodes link_between_2D_nodes" 
% 	int32 NetElemNode(nNetElem,nNetElemMaxNode), shape = [41888 7]
% 		NetElemNode:long_name = "Mapping from net cell to net nodes." 
% 	int32 BndLink(nBndLink), shape = [2208]
% 		BndLink:long_name = "Netlinks that compose the net boundary." 
% 
% 
% //global attributes:
% 		:institution = "Deltares" 
% 		:references = "http://www.deltares.nl" 
% 		:source = "UNSTRUC v1.0.11.11857:1191, model" 
% 		:history = "Created on 2010-07-12T23:48:59+0200, UNSTRUC" 
% 		:Conventions = "CF-1.4:Deltares-0.1" 
% }
