SUBROUTINE MOON (DAY, MOOLON,MOOLAT,MOODIS)
C****
C**** Input: DAY = days measured since 2000 January 1, hour 0
C****
C**** Output: MOOLON = Earth's LONgitude directly beneath the MOOn
C****         MOOLAT = Earth's LATitude directly beneath the MOOn
C****         MOODIS = DIStance from MOOn to Earth in Earth radii units
C****
C**** Constants: OBLIQ = dihedral angle between Earth's equatorial
C****                    plane and orbital plane of year 2000
C****           TAofVE = true anomaly of vernal equinox of year 2000
C****           MAofVE = mean anomaly of vernal equinox of year 2000
C****           EDAYzY = tropical year = Earth days per year = 365.2425
C****           VE200D = days from 2000 January 1, hour 0 till vernal
C****                    equinox of year 2000 = 31 + 29 + 19 + 7.5/24
C****
C**** Intermediate quantities:
C****       RA = right ascension of the Moon =
C****          = longitude of Moon in geocentric equatorial coordinates
C****   VEQLON = longitude of Greenwich Meridion in geocentric
C****            equatorial coordinates at vernal equinox of year 2000
C****   ROTATE = change in longitude in geocentric nonrotating
C****            equatorial coordinates from a point's location on
C****            vernal equinox to is current location where the point
C****            is fixed on rotating Earth
C****
C**** Angles, longitude and latitude are measured in radians unless
C**** the variable is terminated by a capital D in which case it is
C**** measured in degrees
C****
      IMPLICIT  REAL*8 (A-H,M,O-Z)
      PARAMETER (TWOPI=6.283185307179586477d0, PIz180=TWOPI/360.,
     *           OBLIQ=.409101d0, TAofVE=1.345728d0, MAofVE=1.313255d0,
     *           EDAYzY=365.2425d0, VE200D=79.3125)
      INTEGER*4 RADEG,RAMIN, MLDEG,MLMIN
      REAL*8    ECLONN(3,0:6),ECLATN(3,4),HOPARN(3,0:4)
C****
      DATA ECLONN /218.32d0,   0.   ,  481267.883d0,
     *               6.29d0, 134.9d0,  477198.85d0,
     *              -1.27d0, 259.2d0, -413335.38d0,
     *                .66d0, 235.7d0,  890534.23d0,
     *                .21d0, 269.9d0,  954397.70d0,
     *               -.19d0, 357.5d0,   35999.05d0,
     *               -.11d0, 186.6d0,  966404.05d0/
      DATA ECLATN /  5.13d0,  93.3d0,  483202.03d0,
     *                .28d0, 228.2d0,  960400.87d0,
     *               -.28d0, 318.3d0,    6003.18d0,
     *               -.17d0, 217.6d0, -407332.20d0/
      DATA HOPARN / .9508d0,   0.   ,       0.    ,
     *              .0518d0, 134.9d0,  477198.85d0,
     *              .0095d0, 259.2d0, -413335.38d0,
     *              .0078d0, 235.7d0,  890534.23d0,
     *              .0028d0, 269.9d0,  954397.70d0/
C**** Calculate time in Julian centuries
      TJC = (DAY-.5) / (100.*365.25)
C****
C**** Geocentric Ecliptic coordinates:
C**** EcLat = angle between an object and its perpendicular projection
C****         onto the Ecliptic plane with Earth center as angle vertex
C**** EcLon = angle from vernal equinox vector (ray from Earth to Sun
C****         at vernal equinox) to projection onto the ecliptic plane
C****         of vector of current Earth to object
C****
C**** Calculate ecliptic longitude ECLOND in degrees and radians ECLON
      ECLOND = ECLONN(1,0) + ECLONN(3,0)*TJC
      DO 10 I=1,6
   10 ECLOND = ECLOND +
     +         ECLONN(1,I)*SIN(PIz180*(ECLONN(2,I)+ECLONN(3,I)*TJC))
      ECLOND = MOD (ECLOND, 360.d0)
      ECLON  = ECLOND*PIz180
C**** Calculate ecliptic latitude ECLATD in degrees and radians ECLAT
      ECLATD = 0.
      DO 20 I=1,4
   20 ECLATD = ECLATD +
     +         ECLATN(1,I)*SIN(PIz180*(ECLATN(2,I)+ECLATN(3,I)*TJC))
      ECLAT  = ECLATD*PIz180
C**** Calculate horizontal parallax in degrees HOPARD and radians HOPAR
      HOPARD = HOPARN(1,0)
      DO 30 I=1,4
   30 HOPARD = HOPARD +
     +         HOPARN(1,I)*COS(PIz180*(HOPARN(2,I)+HOPARN(3,I)*TJC))
      HOPAR  = HOPARD*PIz180
C**** Calculate MOODIS in Earth radii
      MOODIS = 1. / SIN(HOPAR)
C**** Calculate geocentric ecliptic rectangular coordinates
      XEC = COS(ECLAT)*COS(ECLON)
      YEC = COS(ECLAT)*SIN(ECLON)
      ZEC = SIN(ECLAT)
C****
C**** Geocentric nonrotating Equatorial coordinates:
C**** Declination = angle between an object and its perpendicular
C****               projection onto the Equatorial plane with Earth
C****               center as angle vertex
C**** RightAscention = angle from the vernal equinox vector (ray from
C****                  Earth to Sun at vernal equinox) to projection
C****                  onto equatorial plane of vector of current Earth
C****                  to object
C****
C**** Calculate geocentric equatorial rectangular coordinates
      XEQ = XEC
      YEQ = YEC*COS(OBLIQ) - ZEC*SIN(OBLIQ)
      ZEQ = ZEC*COS(OBLIQ) + YEC*SIN(OBLIQ)
C**** Calculate right ascension and MOOLON
      RA     = ATAN2 (YEQ,XEQ)
      VEQLON = TWOPI*VE200D - TWOPI/2. + MAofVE - TAofVE  !  modulo 2*Pi
      ROTATE = TWOPI*(DAY-VE200D)*(EDAYzY+1.)/EDAYzY
      MOOLON = MODULO (RA-ROTATE-VEQLON, TWOPI)
               IF(MOOLON.gt.TWOPI/2.)  MOOLON = MOOLON - TWOPI
C**** Calculate declination = MLAT in radians
      MOOLAT = ASIN (ZEQ)
      RETURN
      END