function datjd (year,month,day)
! Compute Julian date from citizen year, month and day.
! Work for all years, inspired by Numerical Recipes.
real(double), parameter :: break = 15 + 31*(10 + 12*1582)
real(double), intent( in ) :: year,month,day
real(double) :: datjd,y,m,d,a
if( abs(year) < epsilon(year) ) stop "datjd: there is no year zero."
y = year
if( y < 0.0 ) y = y + 1
if( month > 2.0 )then
m = month + 1.0
d = day
else
y = y - 1
m = month + 13.0
d = day
endif
datjd = int(365.25*y) + int(30.6001*m) + d + 1720994.5
if( d + 31*(m + 12*y) ≥ break )then
a = int(y/100.0)
datjd = datjd + 2.0 - a + int(a/4.0)
endif
end function datjd