SUBROUTINE ZENITH_ANGLE( JD, GEOLAT, GEOLON, P, ZETA, STATUS )
C
C     GIVEN JULIAN DATE JD (UT), THE GEOGRAPHIC LATITUDE AND LONGITUDE
C     (DEG) OF A POINT ON THE EARTH, AND A SOLAR SYSTEM BODY NUMBER, P,
C     THIS ROUTINE RETURNS THE ZENITH ANGLE OF THE BODY, ZETA (DEG).
C
C     B. KNAPP, 1995-03-06
C               1996-02-05 (USE FUNCTION ZENITH_DIST.FOR)
C               2001-04-20 (NEW EPHEM INTERFACE)
C
C     BODY NUMBERS: 0=SUN, 1=MERCURY, 2=VENUS, 4=MARS, 5=JUPITER,
C                   6=SATURN, 7=URANUS, 8=NEPTUNE, 9=PLUTO, 10=MOON
C
C     THE OUTPUT VARIABLE STATUS IS RETURNED NON-ZERO IF AN ERROR
C     OCCURRED.
C
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      IMPLICIT NONE
C
C     INPUT
      REAL*8 JD, GEOLAT, GEOLON
      INTEGER*4 P
C
C     OUTPUT
      REAL*8 ZETA
      INTEGER*4 STATUS
C
C     LOCAL
      REAL*8 PV(6)
C
C     CONSTANTS
      INTEGER*4 APPARENT, EQUATORIAL, SPHERICAL
      PARAMETER (APPARENT=2, EQUATORIAL=1, SPHERICAL=1)
C
C     FUNCTION
      REAL*8 ZENITH_DIST
C
C     EXTERNALS
      EXTERNAL EPHEM, ZENITH_DIST
C
C
C     IN CASE ANYTHING GOES WRONG, RETURN ZETA = -180.
      ZETA = -180.
C
C     GET THE APPARENT RA (PV(1)) AND DEC (PV(2)) OF THE BODY
      CALL EPHEM( JD, P, APPARENT, EQUATORIAL, SPHERICAL, PV, STATUS )
      IF ( STATUS .NE. 0 ) RETURN
C
C     COMPUTE THE ZENITH DISTANCE
      ZETA = ZENITH_DIST( JD, GEOLAT, GEOLON, PV(1), PV(2) )
C
      STATUS = 0
      RETURN
      END