function [tw,lonw,latw,wx,wy,taux,tauy]=loaddailywinds(year);

% function [tw,lons,lats,wx,wy,taux,tauy]=loaddailywinds(year);
%
% load daily NCEP wind speed and stress from netCDF files.

direc='/phod007/lumpkin/ncep/dailywinds/';
wfilename=['wnd.sig995.',num2str(year),'.nc'];
tfilename=['flx.sfc.gauss.',num2str(year),'.nc'];

fprintf('Loading daily NCEP winds for %d.\n',year);

presentdir=pwd;
eval(['cd ',direc]);

% gunzip the files
eval(['!gunzip u',wfilename,'.gz']);
eval(['!gunzip v',wfilename,'.gz']);
eval(['!gunzip u',tfilename,'.gz']);
eval(['!gunzip v',tfilename,'.gz']);

% load time, lon, lat 
eval(['t=getnc(''u',tfilename,''',''time'');']);
eval(['t2=getnc(''u',wfilename,''',''time'');']);
eval(['lons=getnc(''u',tfilename,''',''lon'');']);
eval(['lats=getnc(''u',tfilename,''',''lat'');']);
eval(['lonw=getnc(''u',wfilename,''',''lon'');']);
eval(['latw=getnc(''u',wfilename,''',''lat'');']);

eval(['windx=getnc(''u',wfilename, ...
        ''',''uwnd'',[-1 -1 -1],[-1 -1 -1],', ...
        '[-1 -1 -1],[2 3 1],1);']);
eval(['windy=getnc(''v',wfilename, ...
        ''',''vwnd'',[-1 -1 -1],[-1 -1 -1],', ...
        '[-1 -1 -1],[2 3 1],1);']);
eval(['taux=getnc(''u',tfilename, ...
        ''',''uflx'',[-1 -1 -1],[-1 -1 -1],', ...
        '[-1 -1 -1],[2 3 1],1);']);
eval(['tauy=getnc(''v',tfilename, ...
        ''',''vflx'',[-1 -1 -1],[-1 -1 -1],', ...
        '[-1 -1 -1],[2 3 1],1);']);

% change the sign of taux, tauy: I want the stress
% on the ocean, not on the atmosphere.
taux=-taux;tauy=-tauy;

% gzip the files
eval(['!gzip u',wfilename]);
eval(['!gzip v',wfilename]);
eval(['!gzip u',tfilename]);
eval(['!gzip v',tfilename]);

eval(['cd ',presentdir]);

% most recent year may not have same end
% date for wind speed and wind stress.
% truncate to guarantee same t grid
q=find(t<min(t2) | t>max(t2));
if (length(q))
    t(q)=[]; taux(:,:,q)=[]; tauy(:,:,q)=[];
end;
q=find(t2<min(t) | t2>max(t));
if (length(q));
    t2(q)=[]; windx(:,:,q)=[]; windy(:,:,q)=[];
end;
clear t2 q;

% convert t from hours to julian days
t=(t-17067072)/24+julian([1948,1,1,12,0,0]);

% interpolate wind speeds onto wind stress grid
wx=NaN*ones(size(taux));
wy=wx;

fprintf('interpolating speed onto stress grid ... ');
for i=1:length(t)
    wx(:,:,i)=interp2(lonw,latw',windx(:,:,i),lons,lats');
    wy(:,:,i)=interp2(lonw,latw',windy(:,:,i),lons,lats');
end;
fprintf('done.\n');
windx=wx;windy=wy;
lonw=lons;latw=lats;tw=t;
% fill NaNs from interpolation at edge
wx(:,length(lonw),:)=.5*(wx(:,length(lonw)-1,:)+ ...
    wx(:,1,:));
wy(:,length(lonw),:)=.5*(wy(:,length(lonw)-1,:)+ ...
    wy(:,1,:));

% finally, convert from lon=0:360 to lon=-180:180
q=find(lonw>180);
lonw(q)=lonw(q)-360;
[lonw,q]=sort(lonw);
wx=wx(:,q,:);
wy=wy(:,q,:);
taux=taux(:,q,:);
tauy=tauy(:,q,:);
