      SUBROUTINE NUTATE( T, DPSI, DEPS, DDPSI, DDEPS )
C
C     GIVEN DYNAMIC TIME T (JULIAN MILLENIA SINCE J2000), THIS ROUTINE
C     CALCULATES THE NUTATION IN LONGITUDE, DPSI, AND IN THE OBLIQUITY
C     OF THE ECLIPTIC, DEPS, AS WELL AS THE ASSOCIATED VELOCITY CORREC-
C     TIONS, DDPSI AND DDEPS.  DPSI AND DEPS ARE RETURNED IN RADIANS,
C     AND DDPSI AND DDEPS IN RADIANS/JULIAN MILLENIA.
C
C     REFERENCE: P. K. SEIDELMANN, 1980 IAU THEORY OF NUTATION,
C     CELESTIAL MECHANICS V. 27 (1982), PP. 79-106
C
C     B. KNAPP, 1999-06-11, 2000-09-01, 2001-04-24 (i*4)
C
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      IMPLICIT NONE
C
C     INPUT
      REAL*8 T
C
C     OUTPUT
      REAL*8 DPSI, DEPS, DDPSI, DDEPS
C
C     CONSTANTS
      REAL*8 D2R, S2R
      PARAMETER (D2R = 3.14159265358979324D0/180.D0)
      PARAMETER (S2R = 1.D-4/3600.D0*D2R)
C        (CONVERT UNITS OF 10000 ARC SECONDS TO RADIANS)
C
C     LOCAL VARIABLES
      REAL*8 M, E, F, D, O, DM, DE, DF, DD, DO, L, DL,
     &   SL, CL, PSI_K, EPS_K
      INTEGER*4 K
C
C     NUTATION COEFFICIENTS (CONSTANTS)
      INTEGER*4 NTERMS
      PARAMETER (NTERMS=106)
C
      REAL*8 CM(NTERMS)   /
     &    0,  0,  0,  0,  0,  1,  0,  0,  1,  0,  1,  0, -1,  1, -1, -1,
     &    0,  1, -2,  2,  0,  2,  1,  2,  0, -1,  0, -1,  0,  0,  0,  1,
     &    0, -1,  2,  1,  0,  0,  0,  1, -2,  2,  0,  1,  1,  2,  0,  0,
     &    1,  2,  0,  0,  1,  0,  1,  3,  1, -1,  0, -2,  1,  1,  1,  0,
     &    1, -2, -1,  2,  0, -1,  1,  3, -2,  1,  2, -1,  1, -2,  0,  0,
     &    0,  2, -1,  0,  1,  2,  1,  0,  3,  1,  0, -1,  0,  0,  0,  1,
     &    0,  1,  2, -1,  2,  1,  1,  2,  0,  0
     &                    /
      REAL*8 CE(NTERMS)   /
     &    0,  0,  0,  0,  1,  0,  1,  0,  0, -1,  0,  0,  0,  0,  0,  0,
     &    0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  2,  2,  0,
     &   -1,  0,  0,  0,  1, -1,  0,  1,  0,  0,  0,  0,  0,  0,  0, -1,
     &   -1,  0,  1,  1,  0,  0,  0,  0, -1, -1, -1,  0, -1,  0,  1,  0,
     &    0,  0,  0,  0, -2,  0,  1,  0,  0,  0,  0,  0,  1,  0,  1,  1,
     &    0,  1,  0,  1,  0,  0, -1,  0,  0,  0,  1, -1,  0,  0,  1,  0,
     &   -1,  1,  0,  0,  0,  0,  0,  0,  0,  1
     &                    /
      REAL*8 CF(NTERMS)   /
     &    0,  2,  2,  0,  0,  0,  2,  2,  2,  2,  0,  2,  2,  0,  0,  2,
     &    0,  2,  2,  0,  2,  2,  2,  0,  2,  2,  2,  0,  0,  2,  0,  0,
     &    0,  2, -2,  2,  2,  2,  2,  0,  0,  2,  0,  2,  0,  2,  0,  2,
     &    0,  0,  2,  0, -2,  0,  0,  2,  2,  2,  2,  2,  0,  2,  0,  2,
     &    0,  0,  2,  0,  2,  2,  2,  0,  2,  2,  2,  0,  2,  2,  0, -2,
     &   -2,  0,  0,  2,  0,  2,  0,  4,  2,  2,  2,  0, -2,  2,  0, -2,
     &    2,  0, -2,  4,  0,  0, -2,  0,  2,  0
     &                    /
      REAL*8 CD(NTERMS)   /
     &    0, -2,  0,  0,  0,  0, -2,  0,  0, -2, -2, -2,  0,  0,  0,  2,
     &    2,  0,  0, -2,  2,  0, -2,  0,  0,  0, -2,  2,  0, -2,  0, -2,
     &    0,  2,  0,  2,  0,  0,  2, -2,  2, -2,  2, -2,  2,  0, -2, -2,
     &    0, -2, -2, -2,  0,  1, -1,  0,  0,  2,  2,  0, -1,  0,  0,  1,
     &    0,  0,  4,  0, -2, -2,  0,  0,  2,  2, -2,  0, -2,  4,  0,  2,
     &    2, -2,  1, -2,  2,  2, -2, -2, -2, -2,  0,  2,  0, -1,  2, -2,
     &    0, -2,  0,  0, -4, -4,  2,  2,  4,  1
     &                    /
      REAL*8 CO(NTERMS)   /
     &    1,  2,  2,  2,  0,  0,  2,  1,  2,  2,  0,  1,  2,  1,  1,  2,
     &    0,  1,  1,  0,  2,  2,  2,  0,  0,  1,  0,  1,  1,  2,  0,  1,
     &    1,  1,  0,  2,  2,  2,  1,  0,  1,  2,  1,  1,  0,  1,  1,  1,
     &    0,  1,  1,  0,  0,  0,  0,  2,  2,  2,  2,  2,  0,  0,  0,  2,
     &    2,  1,  2,  1,  1,  1,  2,  0,  2,  1,  1,  2,  2,  2,  2,  0,
     &    1,  0,  1,  0,  1,  2,  0,  2,  2,  0,  1,  1,  1,  2,  0,  0,
     &    1,  1,  1,  2,  0,  0,  0,  0,  2,  0
     &                    /
      REAL*8 PSI0(NTERMS) /
     &  -171996, -13187,  -2274,   2062,   1426,    712,   -517,   -386,
     &     -301,    217,   -158,    129,    123,     63,    -58,    -59,
     &       63,    -51,     46,     48,    -38,    -31,     29,     29,
     &       26,     21,    -22,     16,    -15,    -16,     17,    -13,
     &      -12,    -10,     11,     -8,      7,     -7,     -7,     -7,
     &       -6,      6,     -6,      6,      6,     -5,     -5,     -5,
     &        5,      4,      4,     -4,      4,     -4,     -4,     -3,
     &       -3,     -3,     -3,     -3,     -3,      3,     -3,      2,
     &       -2,     -2,     -2,      2,     -2,     -2,      2,      2,
     &        1,     -1,      1,      1,      1,     -1,      1,     -1,
     &        1,      1,      1,     -1,     -1,     -1,      1,      1,
     &        1,     -1,      1,      1,     -1,     -1,     -1,     -1,
     &       -1,     -1,      1,      1,     -1,     -1,     -1,      1,
     &       -1,      1
     &                    /
      REAL*8 PSI1(NTERMS) /
     &    -1742,    -16,     -2,      2,    -34,      1,     12,     -4,
     &        0,     -5,      0,      1,      0,      1,      1,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      1,     -1,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0
     &                    /
      REAL*8 EPS0(NTERMS) /
     &    92025,   5736,    977,   -895,     54,     -7,    224,    200,
     &      129,    -95,     -1,    -70,    -53,    -33,     32,     26,
     &       -2,     27,    -24,      1,     16,     13,    -12,     -1,
     &       -1,    -10,      0,     -8,      9,      7,      0,      7,
     &        6,      5,      0,      3,     -3,      3,      3,      0,
     &        3,     -3,      3,     -3,      0,      3,      3,      3,
     &        0,     -2,     -2,      0,      0,      0,      0,      1,
     &        1,      1,      1,      1,      0,      0,      0,     -1,
     &        1,      1,      1,     -1,      1,      1,     -1,      0,
     &       -1,      1,     -1,     -1,     -1,      1,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0
     &                    /
      REAL*8 EPS1(NTERMS) /
     &       89,    -31,     -5,      5,     -1,      0,     -6,      0,
     &       -1,      3,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0,      0,      0,      0,      0,      0,      0,
     &        0,      0
     &                    /
C
C     POLYNOMIALS (IN DEGREES) FOR THE FUNDAMENTAL ARGUMENTS (CONSTANTS)
      INTEGER*4 PDEG
      PARAMETER (PDEG=3)
C
C     MEAN ANOMALY OF THE MOON
      REAL*8 PM(0:PDEG) /  1.349629813888889D+2,  4.771988673980556D+6,
     &                     8.697222222222222D-1,  1.777777777777778D-2 /
C
C     MEAN ANOMALY OF THE SUN (EARTH)
      REAL*8 PE(0:PDEG) /  3.575277233333333D+2,  3.599905034000000D+5,
     &                    -1.602777777777778D-2, -3.333333333333333D-3 /
C
C     MOON'S ARGUMENT OF LATITUDE
      REAL*8 PF(0:PDEG) /  9.327191027777778D+1,  4.832020175380556D+6,
     &                    -3.682500000000000D-1,  3.055555555555556D-3 /
C
C     MEAN ELONGATION OF THE MOON FROM THE SUN
      REAL*8 PD(0:PDEG) /  2.978503630555556D+2,  4.452671114800000D+6,
     &                    -1.914166666666667D-1,  5.277777777777778D-3 /
C
C     LONGITUDE OF THE ASCENDING NODE OF THE MOON'S MEAN ORBIT ON THE
C     ECLIPTIC, MEASURED FROM THE MEAN EQUINOX OF THE DATE
      REAL*8 PO(0:PDEG) /  1.250445222222222D+2, -1.934136260833333D+4,
     &                     2.070833333333333D-1,  2.222222222222222D-3 /
C
      SAVE CM, CE, CF, CD, CO, PSI0, PSI1, EPS0, EPS1,
     &     PM, PE, PF, PD, PO
C
C     EVALUATE THE FUNDAMENTAL ARGUMENTS
      M  = PM(0)+T*(PM(1)+T*(     PM(2)+T*     PM(3)))
      DM =          PM(1)+T*(2.D0*PM(2)+T*3.D0*PM(3))
      E  = PE(0)+T*(PE(1)+T*(     PE(2)+T*     PE(3)))
      DE =          PE(1)+T*(2.D0*PE(2)+T*3.D0*PE(3))
      F  = PF(0)+T*(PF(1)+T*(     PF(2)+T*     PF(3)))
      DF =          PF(1)+T*(2.D0*PF(2)+T*3.D0*PF(3))
      D  = PD(0)+T*(PD(1)+T*(     PD(2)+T*     PD(3)))
      DD =          PD(1)+T*(2.D0*PD(2)+T*3.D0*PD(3))
      O  = PO(0)+T*(PO(1)+T*(     PO(2)+T*     PO(3)))
      DO =          PO(1)+T*(2.D0*PO(2)+T*3.D0*PO(3))
C
C     COMPUTE THE SUM OF ALL TERMS
      DPSI = 0.D0
      DEPS = 0.D0
      DDPSI = 0.D0
      DDEPS = 0.D0
      DO 1 K=NTERMS,1,-1
         L  = MOD( CM(K)* M + CE(K)* E + CF(K)* F + CD(K)* D + CO(K)* O,
     &      360.D0 )*D2R
         DL =    ( CM(K)*DM + CE(K)*DE + CF(K)*DF + CD(K)*DD + CO(K)*DO
     &             )*D2R
         SL = SIN( L )
         CL = COS( L )
         PSI_K = PSI0(K) + PSI1(K)*T
         EPS_K = EPS0(K) + EPS1(K)*T
         DPSI = DPSI + PSI_K*SL
         DEPS = DEPS + EPS_K*CL
         DDPSI = DDPSI + PSI_K*CL*DL + SL*PSI1(K)
         DDEPS = DDEPS - EPS_K*SL*DL + CL*EPS1(K)
    1 CONTINUE
C
      DPSI = DPSI*S2R
      DEPS = DEPS*S2R
      DDPSI = DDPSI*S2R
      DDEPS = DDEPS*S2R
C
      RETURN
      END
