function ncwind2delft3d(srcpath,tgtpath,tgtbase)
% NCWIND2DELFT3D convert wind NetCDF files to Delft3D meteo files
%    NCWIND2DELFT3D(SRCPATH,TGTPATH,TGTBASE) reads files wind_speed.nc and
%    wind_direction.nc in SRCPATH, and writes TGTBASE.amu, TGTBASE.amv, and
%    optionally TGTBASE.grd files in TGTPATH folder.
%
%    This function requires that WLGRID (included in Delft3D-MATLAB
%    toolbox) is on the search path as well as the MEXNC and SNCTOOLS
%    toolboxes needed to read the NetCDF file.
%
%    Example:
%    OSS_TRUNK = 'd:/delft3d_trunk';
%    % Delft3D-MATLAB toolbox
%    addpath([OSS_TRUNK '/src/tools_lgpl/matlab/quickplot/progsrc'])
%    % NetCDF toolboxes
%    addpath([OSS_TRUNK '/src/third_party_open/netcdf/matlab/mexnc'])
%    addpath([OSS_TRUNK '/src/third_party_open/netcdf/matlab/snctools'])
%    javaaddpath([OSS_TRUNK '/src/third_party_open/netcdf/matlab/netcdfAll-4.1.jar'])
%    ncwind2delft3d('d:/','d:/','wind')

ncsp = fullfile(srcpath,'wind_speed.nc');
ncdr = fullfile(srcpath,'wind_direction.nc');
grd  = fullfile(tgtpath,[tgtbase '.grd']);
uwnd = fullfile(tgtpath,[tgtbase '.amu']);
vwnd = fullfile(tgtpath,[tgtbase '.amv']);
nodata = 999.999;

% get lon/lat
lon = nc_varget(ncdr,'LON791_819');
lat = nc_varget(ncdr,'LAT996_1023');

% get data
wspeed   = nc_varget(ncsp,'WS');
wfromdir = nc_varget(ncdr,'WD');
% assume wfromdir is defined in degrees from the north, positive in
% clockwise direction
%
%             0
%             |
% ==>  -90 ---*--- 90  <==
%             |
%            180
%
% tranform to wtodir in radians from east, positive in anti-clockwise
% direction
%
%           pi/2
%            |
% <==  pi ---*--- 0  ==>
%            |
%          -pi/2
%
wtodir   = (90-wfromdir)*pi/180-pi;
% compute wind velocity components
u     = wspeed.*cos(wtodir);
v     = wspeed.*sin(wtodir);

% check if lon/lat are one dimensional
curvi = sum(size(lon)>1)>1;
if curvi
    % NO: curvilinear grid
    wlgrid('write','filename',grd,'coordinatesystem','spherical','x',lon,'y',lat)
else
    % YES: rectilinear grid
    dlon = diff(lon);
    dlat = diff(lat);
    if all(abs(dlon-dlon(1))<1e-7) && all(abs(dlat-dlat(1))<1e-7)
        % simple equidistant rectilinear grid
    else
        % non-equidistant rectilinear grid -> should be treated as curvi
        % grid fof Delft3D
        curvi = true;
        wlgrid('write','filename',grd,'coordinatesystem','spherical','x',lon,'y',lat)
    end
end

dataline_format = [repmat('%7.3f ',1,size(u,1)) '\n'];
u(isnan(u)) = nodata;
v(isnan(v)) = nodata;

% open files
ufid = fopen(uwnd,'w');
vfid = fopen(vwnd,'w');
uvfd = [ufid vfid];

% write file headers
mfprintf(uvfd,'FileVersion      = 1.03\n');
if curvi
    mfprintf(uvfd,'Filetype         = meteo_on_curvilinear_grid\n');
    mfprintf(uvfd,'grid_file        = wind.grd\n');
    mfprintf(uvfd,'first_data_value = grid_ulcorner\n');
    mfprintf(uvfd,'data_row         = grid_row\n');
else
    mfprintf(uvfd,'Filetype         = meteo_on_equidistant_grid\n');
    mfprintf(uvfd,'n_cols           = %i\n',size(u,1));
    mfprintf(uvfd,'n_rows           = %i\n',size(u,2));
    mfprintf(uvfd,'grid_unit        = degree\n');
    % code currently assumes lon and lat are increasing
    mfprintf(uvfd,'x_llcenter       = %g\n',lon(1));
    mfprintf(uvfd,'dx               = %g\n',dlon(1));
    mfprintf(uvfd,'y_llcenter       = %g\n',lat(1));
    mfprintf(uvfd,'dy               = %g\n',dlat(1));
end
mfprintf(uvfd,'NODATA_value     = %7.3f\n',nodata);
mfprintf(uvfd,'n_quantity       = 1\n');
fprintf (ufid,'quantity1        = x_wind\n');
fprintf (vfid,'quantity1        = y_wind\n');
mfprintf(uvfd,'unit1            = m s-1\n');

% write time blocks
mfprintf(uvfd,'TIME = 0 minutes since 1900-01-01 00:00:00 +00:00\n');
fprintf (ufid,dataline_format,fliplr(u));
fprintf (vfid,dataline_format,fliplr(v));
mfprintf(uvfd,'TIME = 0 minutes since 2100-01-01 00:00:00 +00:00\n');
fprintf (ufid,dataline_format,fliplr(u));
fprintf (vfid,dataline_format,fliplr(v));

% close files
fclose(ufid);
fclose(vfid);

% support function to write two (or more) similar files
function mfprintf(fids,varargin)
for i=1:length(fids)
    fprintf(fids(i),varargin{:});
end