subroutine julian(yr,day)
! returns present yr and day of year
integer i,day,yr,vec(8),molen(12)
data molen /31,28,31,30,31,30,31,31,30,31,30,31/
call date_and_time(values=vec)
molen(2)=28   ;yr=vec(1)  ; day=vec(3)  !day of mo
if(mod(yr,4).eq.0) then
 molen(2)=29          !normal leap
 if(mod(yr,100).eq.0) molen(2)=28
 if(mod(yr,400).eq.0) molen(2)=29
endif
do i=1,vec(2)-1
 day=day+molen(i)   !sum preceeding month lengths
enddo
return
end