function [dyear] = day2year(day, dref)
%
%  Converts a number of days since dref into a decimal year
%  value. dref is the base year.
%
%  Written by Tom Sgouros.
%
% $Id: day2year.m,v 1.8 1999/05/03 17:40:59 tom Exp $


% The preceding empty line is important.

% $Log: day2year.m,v $
% Revision 1.8  1999/05/03 17:40:59  tom
% fixed spelling of my name.
%
% Revision 1.7  1998/11/30 00:49:09  dbyrne
%
%
% Fixed bug in unpack that deleted all returning names for datasets with
% no empty strings in the required fields.  Fixed vectorization of day2year
% and add workaround for a FLOP bug I found on ix86 machines.  Fixed calendar
% day to yearday conversion error in nscat30cat.m.  Added workaround for writeval
% bug in nscat30.m (2 server lines force 2 loaddods calls).  Fixed a bunch
% of bugs in new SSM/I code submitted by Rob Morris, JPL.  Added SSM/I, CZCS,
% Levitus (1982), World Ocean Atlas monthly, seasonal, and annual datasets and
% Mauna Loa CO2 record.  Fixed a logic error in urieccat.m. -- dbyrne, 98/11/29
%
% Revision 1.1  1998/05/17 14:10:44  dbyrne
% *** empty log message ***
%
% Revision 1.2  1997/11/26 19:01:49  jimg
% Changed comment about dref.
%
% Revision 1.1  1997/11/26 18:51:27  jimg
% Added.
%

if any(floor(day) == 1)
  % catch the first of the year and *force* it to be an integer
  % this catches a bug present on PC ix86 chips that I found for
  % converstion of Jan 1,1997.
  i = find(floor(day) == 1);
  dyear = nan*day;
  dyear(i) = dref(i)+(day(i)-floor(day(i)))/(365+isleap(dref));
end
k = find(floor(day) ~= 1);
dref = dref(k);
day = day(k);

ndays = [0   36524   73048  109572  146097  182621 ...
      219145  255669  292194  328718  365242  401766 ...
      438291  474815  511339  547863  584388  620912 ...
      657436  693960  730485];
for j = 1:length(k)
  if dref(j) ~= 1
    offset(j) = year2day(dref(j), 1);
    day(j) = day(j) + offset(j);
  end
  i = find(ndays < day(j));
  if ~isempty(i)
    i = max(i);
  else
    i = 21;
  end
  diff(j) = day(j) - ndays(i);
  nleaps(j) = floor(diff(j)/1461);
  dyear(k(j)) = floor(diff(j)/365);
  diff(j) = (diff(j)-nleaps(j)) - (dyear(k(j))*365);
  dyear(k(j)) = dyear(k(j)) + 1 + (100*(i-1));
  
  if isleap(dyear(k(j)))
    dyear(k(j)) = dyear(k(j)) + diff(j)/366.0;
  else
    dyear(k(j)) = dyear(k(j)) + diff(j)/365.0;
  end
end
return
