function [doy]=day_of_year(year, month, day)
% Calculates doy from year, month day
mon_length=[31,28,31,30,31,30,31,31,30,31,30,31]';
day_0=[0; cumsum(mon_length(1:11))];
day_0v=[day_0(1:2);day_0(3:12)+1];
f_v=find( year==400*floor(year/400) | (year==4*floor(year/4) & year~=100*floor(year/100)));
month_v=month(f_v);
day_v=day(f_v);
doy=day_0(month)+day;
doy(f_v)=day_0v(month_v)+day_v;
