c This function returns the julian date of easter at 12 noon for
c year in the gregorian calendar.
c (see also: http://aa.usno.navy.mil/AA/)
real*8 function easter (year)
integer year
ia = mod (year, 19)
ib = year / 100
ic = year - ib * 100
id = ib / 4
ie = ib - id * 4
if = (ib + 8) / 25
ig = (ib - if + 1) / 3
ih = mod (ia * 19 + ib - id - ig + 15 , 30)
ii = ic / 4
ik = ic - ii * 4
il = mod (32 + 2 * ie + 2 * ii - ih - ik , 7)
im = (ia + 11 * ih + 22 * il) / 451
iq = ih + il - 7 * im + 22
number_of_leap_years = year / 4 + year / 400 - year / 100
easter = 1721119. d0 + year * 365.d0
. + number_of_leap_years + iq
return
end