Algorithm for the computation of the Easter date

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