subroutine lpteop( time_mjd, dut, dlod )

*  Routine for computing variations in earth rotation caused by
*  long-period (zonal) body tides and ocean tides.
*  For a given time, the routine returns the tidally induced 
*  variation in UT1 and in length-of-day (LOD).
*
*  Arguments:
*
*    time_mjd - IN -  desired time, in decimal modified Julian day.
*    dut      - OUT - anomaly in UT1-TAI, in microseconds.
*    dlod     - OUT - excess LOD, in microseconds (per day).
*
*  Notes:
*   1. The routine that computes required astronomical mean longitudes
*     here ignores the difference between UTC and TT, which is acceptable
*     in this particular application.
*   2. Currently reads an external table of 80 long-period spectral lines.
*     User must modify OPEN statement to point to this dataset.
*
*  R. Ray  NASA/GSFC  12 May 2012
*

      implicit none
      double precision, intent(in)  :: time_mjd
      double precision, intent(out) :: dut, dlod

*  local variables...
      integer    mt
      parameter (mt=80)
      double precision a(2,mt),b(2,mt),shpnp(5),beta(7),akk(7),s1,s2
      double precision a1,a2,b1,b2,prd,cteamp,arg,tmjd
      integer          k(7,mt), kk(7),lu,i,j
      double precision pi,rad
      parameter (pi=3.14159265359d0)
      parameter (rad=pi/180.d0)
      save init,k,a,b
      logical init
      data init/.true./

      tmjd = time_mjd

      if (init) then
         init = .false.
***********************************************************************
*
*  MODIFY THIS OPEN STATEMENT TO FIND  table.d  ON YOUR MACHINE
*
         lu = 2
         open( lu,  file='./table.d' ,  status='old' )

***********************************************************************
***********************************************************************
         read(lu,*)
         do i=1,mt
            read(lu,7) j,kk,prd,cteamp,a1,a2,b1,b2
 7          format(i3,1x,7i2,f12.0,f11.0,2f11.0,2f9.0)
            k(:,i) = kk
            a(1,i) = a1
            a(2,i) = a2
            b(1,i) = b1
            b(2,i) = b2
         end do
         close(lu)
      endif

      call astro5z( tmjd, shpnp )
      beta(1) = 0.d0
      beta(7) = 90.d0
      beta(2:6) = shpnp
      beta(5) = -beta(5)  ! The 5th Doodson number is for N' = -N.
      s1=0.d0; s2=0.d0
      do i=1,mt
         akk = k(:,i)
         arg = DOT_PRODUCT( beta, akk )*rad
         s1 = s1 + a(1,i)*cos(arg) + a(2,i)*sin(arg)
         s2 = s2 + b(1,i)*cos(arg) + b(2,i)*sin(arg)
      end do
      dut = s1
      dlod = s2
      return
      end

*****************************************************************


      SUBROUTINE ASTRO5Z( TIME, SHPNP )
*
*---------------------------------------------------------------------
*  Computes the 5 basic astronomical mean longitudes  s, h, p, N, p'.
*
*  Note N is not N', i.e. N is decreasing with time.
*
*  TIME is UTC in decimal Modified Julian Day (MJD).
*  All longitudes returned in degrees.
*
*  R. D. Ray, NASA/GSFC   August 2003
*
*  Most of the formulae for mean longitudes are extracted from 
*  Jean Meeus, Astronomical Algorithms, 2nd ed., 1998.  
*  Page numbers below refer to this book.
*
*  Warning: This version of ASTRO5 ignores ET - UTC difference!
*           (ok for long-period tides)
*
*---------------------------------------------------------------------
*
      IMPLICIT NONE
      DOUBLE PRECISION TIME, SHPNP(5)

      DOUBLE PRECISION TJD,T,CIRCLE,DEL,TJLAST,D
      PARAMETER       (CIRCLE=360.0D0)

*     Compute time argument in centuries relative to J2000
*     ----------------------------------------------------
      TJD = TIME + 2400000.5D0
      T = ( TJD - 2451545.d0 )/36525.d0
*
*     mean longitude of moon (p.338)
*     ------------------------------
      SHPNP(1) = (((-1.53388d-8*T + 1.855835d-6)*T - 1.5786d-3)*T +
     *            481267.88123421d0)*T + 218.3164477d0
*
*     mean elongation of moon (p.338)
*     -------------------------------
      D = (((-8.8445d-9*T + 1.83195d-6)*T - 1.8819d-3)*T +
     *          445267.1114034d0)*T + 297.8501921d0
*
*     mean longitude of sun
*     ---------------------
      SHPNP(2) = SHPNP(1) - D
*
*     mean longitude of lunar perigee (p.343)
*     ---------------------------------------
      SHPNP(3) = ((-1.249172d-5*T - 1.032d-2)*T + 4069.0137287d0)*T +
     *            83.3532465d0
*
*     mean longitude of ascending lunar node (p.144)
*     ----------------------------------------------
      SHPNP(4) = ((2.22222d-6*T + 2.0708d-3)*T - 1934.136261d0)*T +
     *            125.04452d0
*
*     mean longitude of solar perigee (Simon et al., 1994)
*     ----------------------------------------------------
      SHPNP(5) = 282.94d0 + 1.7192d0 * T

      SHPNP = MODULO( SHPNP, CIRCLE )
      RETURN
      END