      SUBROUTINE PMSDNUT2 (RMJD, PM)
*+
*  - - - - - - - - - - -
*   P M S D N U T 2
*  - - - - - - - - - - -
*
*  This routine is part of the International Earth Rotation and
*  Reference Systems Service (IERS) Conventions software collection.
*
*  This subroutine evaluates the model of polar motion for
*  a nonrigid Earth due to tidal gravitation. This polar motion
*  is equivalent to the so-called "subdiurnal nutation." The model
*  is a sum of a first order polynomial and 25 trigonometric terms
*  (15 long periodic and 10 quasi diurnal) with coefficients given
*  in Table 5.1a of the IERS Conventions (2010).
*
*     :------------------------------------------:
*     :                                          :
*     :                 IMPORTANT                :
*     :                                          :
*     : In the present version this subroutine   :
*     : neglects the linear trend and the long   :
*     : periodic terms of the expansion, for the :
*     : reasons explained in Section 5.5.1.1 of  :
*     : the IERS Conventions (2010).             :
*     : If the full expansion is needed, set the :
*     : parameter iband to 0 instead of 1,       :
*     : that is, replace the statement           :
*     :     PARAMETER ( iband = 1 )              :
*     : to  PARAMETER ( iband = 0 )              :            
*     :                                          :
*     :__________________________________________:
*
*  In general, Class 1, 2, and 3 models represent physical effects that
*  act on geodetic parameters while canonical models provide lower-level
*  representations or basic computations that are used by Class 1, 2, or
*  3 models.
* 
*  Status:  Class 1 model
*
*     Class 1 models are those recommended to be used a priori in the
*     reduction of raw space geodetic data in order to determine
*     geodetic parameter estimates.
*     Class 2 models are those that eliminate an observational
*     singularity and are purely conventional in nature.
*     Class 3 models are those that are not required as either Class
*     1 or 2.
*     Canonical models are accepted as is and cannot be classified as
*     a Class 1, 2, or 3 model.
*
*  Given:
*     rmjd        d      Time expressed as modified Julian date
*
*  Returned:
*     pm          d(2)      Vector of length 2 (Note 1)
*
*  Notes:
*
*  1) The polar motion coordinates (dx, dy) are expressed in
*     microarcseconds.
*
*  Called:
*     FUNDARG             Compute the angular fundamental arguments
*
*  Test case:
*     given input: rmjd = 54335D0 ( August 23, 2007 ) 
*
*     expected output: (dx) pm(1)  = 24.83144238273364834D0 microarcseconds
*                      (dy) pm(2) = -14.09240692041837661D0 microarcseconds
*
*  References:
*
*     Petit, G. and Luzum, B. (eds.), IERS Conventions (2010),
*     IERS Technical Note No. 36, BKG (2010)
*
*  Revisions:
*  2005 March       A.Brzezinski   Original code
*  2008 November 26 B.E.Stetzler   Initial changes to code
*  2008 December 01 B.E.Stetzler   Provided test case
*  2009 August   18 B.E.Stetzler   Capitalized all variables for FORTRAN
*                                  77 compatibility
*  2010 May      14 B.E.Stetzler   Replaced call to PMARGS to FUNDARG
*                                  for universal fundamental argument
*                                  subroutine
*  2010 May      17 B.E.Stetzler   Validated test case using internally
*                                  computed GMST and call to FUNDARG
*                                  matched previous external call to
*                                  PMARGS for all six parameters
*  2010 June     23 B.E.Stetzler   Modified coefficients of the long
*                                  and short period terms in polar 
*                                  motion and secular polar motion
*                                  rate to coincide with Table 5.1a
*  2011 October  13 B.E. Stetzler  Modification above yields difference
*                                  in output between the original PMsdnut
*                                  and PMSDNUT2.F at the 10^-2 muas level
*                                  Test case header corrected to match
*                                  output from PMSDNUT2.F
*-----------------------------------------------------------------------

      IMPLICIT NONE
      DOUBLE PRECISION RMJD, PM(2)

*         ----------------------------
*           D E F I N I T I O N S
*         ----------------------------
*  iband  - parameter defining the range of periods for the terms which
*           are included in computations; if equal to 1 only the quasi 
*           diurnal terms are computed, otherwise the full model
*  iarg   - array defining for each of the 25 trigonometric terms a set
*           of 6 integer multipliers of the fundamental angular arguments
*  arg    - vector of the following 6 fundamental arguments used to
*           compute the angular argument of the trigonometric functions
*           arg(1:6) = [ GMST+pi, el, elp, f, d, om ]; this vector is
*           evaluated by the subroutine FUNDARG which is called as an 
*           external subroutine.  Originally evaluated by the subroutine
*           PMARGS. 
*  period - array of periods of the trigonometric terms of expansion, in
*           mean solar days; only for a check - not used in computations
*  xs, xc - sine and cosine coefficients of the x coordinate of the pole,
*           in microarcseconds
*  ys, yc - sine and cosine coefficients of the y coordinate of the pole,
*           in microarcseconds
*  angle  - angular argument of the trigonometric functions
*           angle = Sum(i=1:6) iarg(i,j)*arg(i), for j=1,25

      INTEGER IBAND, I, J, JSTART
      PARAMETER ( IBAND = 1 )
      INTEGER IARG(6,25)
      DOUBLE PRECISION T, GMST, L, LP, F, D, OM
      DOUBLE PRECISION ARG(6)
      DOUBLE PRECISION PER(25), XS(25), XC(25), YS(25), YC(25)
      DOUBLE PRECISION ANGLE, XRATE, YRATE

* Set constants

*  Arcseconds to radians
      DOUBLE PRECISION DAS2R
      PARAMETER ( DAS2R = 4.848136811095359935899141D-6 )

*  Arcseconds in a full circle
      DOUBLE PRECISION TURNAS
      PARAMETER ( TURNAS = 1296000D0 )

*  rmjd0   - modified Julian date of J2000
*  twopi   - 2*pi

      DOUBLE PRECISION RMJD0, PI, TWOPI
      PARAMETER ( RMJD0   = 51544.5D0                )
      PARAMETER ( PI      = 3.141592653589793238462643D0 )
      PARAMETER ( TWOPI   = 6.283185307179586476925287D0 )

*  Radians to seconds
      DOUBLE PRECISION RAD2SEC
      PARAMETER ( RAD2SEC = 86400D0/TWOPI            )

* Coefficients of the long periodic terms in polar motion
* Source: IERS Conventions (2010), Table 5.1a

      DATA 
     . (  (IARG(I,J),I=1,6),    PER(J),   XS(J),  XC(J),  YS(J),  YC(J),
     .                                                           J=1,15)
     ./ 0,  0, 0,  0,  0, -1, 6798.3837,    0.0,   0.6,   -0.1,   -0.1,
     .  0, -1, 0,  1,  0,  2, 6159.1355,    1.5,   0.0,   -0.2,    0.1,
     .  0, -1, 0,  1,  0,  1, 3231.4956,  -28.5,  -0.2,    3.4,   -3.9,
     .  0, -1, 0,  1,  0,  0, 2190.3501,   -4.7,  -0.1,    0.6,   -0.9,
     .  0,  1, 1, -1,  0,  0, 438.35990,   -0.7,   0.2,   -0.2,   -0.7,
     .  0,  1, 1, -1,  0, -1, 411.80661,    1.0,   0.3,   -0.3,    1.0,
     .  0,  0, 0,  1, -1,  1, 365.24219,    1.2,   0.2,   -0.2,    1.4,
     .  0,  1, 0,  1, -2,  1, 193.55971,    1.3,   0.4,   -0.2,    2.9,
     .  0,  0, 0,  1,  0,  2, 27.431826,   -0.1,  -0.2,    0.0,   -1.7,
     .  0,  0, 0,  1,  0,  1, 27.321582,    0.9,   4.0,   -0.1,   32.4,
     .  0,  0, 0,  1,  0,  0, 27.212221,    0.1,   0.6,    0.0,    5.1,
     .  0, -1, 0,  1,  2,  1, 14.698136,    0.0,   0.1,    0.0,    0.6,
     .  0,  1, 0,  1,  0,  1, 13.718786,   -0.1,   0.3,    0.0,    2.7,
     .  0,  0, 0,  3,  0,  3, 9.1071941,   -0.1,   0.1,    0.0,    0.9,
     .  0,  0, 0,  3,  0,  2, 9.0950103,   -0.1,   0.1,    0.0,    0.6/

* Coefficients of the quasi diurnal terms in polar motion
* Source: IERS Conventions (2010), Table 5.1a

      DATA 
     .(  (IARG(I,J),i=1,6),     PER(J),   XS(J),  XC(J),  YS(J),  YC(J),
     .                                                          J=16,25)
     ./ 1, -1, 0, -2,  0, -1, 1.1196992,   -0.4,   0.3,   -0.3,   -0.4,
     .  1, -1, 0, -2,  0, -2, 1.1195149,   -2.3,   1.3,   -1.3,   -2.3,
     .  1,  1, 0, -2, -2, -2, 1.1134606,   -0.4,   0.3,   -0.3,   -0.4,
     .  1,  0, 0, -2,  0, -1, 1.0759762,   -2.1,   1.2,   -1.2,   -2.1,
     .  1,  0, 0, -2,  0, -2, 1.0758059,  -11.4,   6.5,   -6.5,  -11.4,
     .  1, -1, 0,  0,  0,  0, 1.0347187,    0.8,  -0.5,    0.5,    0.8,
     .  1,  0, 0, -2,  2, -2, 1.0027454,   -4.8,   2.7,   -2.7,   -4.8,
     .  1,  0, 0,  0,  0,  0, 0.9972696,   14.3,  -8.2,    8.2,   14.3,
     .  1,  0, 0,  0,  0, -1, 0.9971233,    1.9,  -1.1,    1.1,    1.9,
     .  1,  1, 0,  0,  0,  0, 0.9624365,    0.8,  -0.4,    0.4,    0.8/

* Rate of secular polar motion, in microarcseconds per year
* Source: IERS Conventions (2010), Table 5.1a

      DATA XRATE, YRATE / -3.8, -4.3/

* Compute the periodical part of the model
* Coordinates of the pole are set to zero first
      PM(1) = 0D0
      PM(2) = 0D0

* Evaluate the vector of the fundamental arguments
* arg(1:6) = [ GMST+pi, el, elp, f, d, om ] at t = rmjd

*  Convert the input epoch to Julian centuries of TDB since J2000
      T = (RMJD-RMJD0)/36525D0

*  Compute GMST + pi
      GMST = MOD (   67310.54841D0 +
     .               T*( (8640184.812866D0 + 3155760000D0) +
     .               T*( 0.093104D0 +
     .               T*( -0.0000062 ))), 86400D0 )

      CALL FUNDARG ( T, L, LP, F, D, OM )

      ARG(1) = GMST / RAD2SEC + PI
      ARG(1) = DMOD( ARG(1), TWOPI )
      ARG(2) = L
      ARG(3) = LP
      ARG(4) = F
      ARG(5) = D
      ARG(6) = OM 

      IF (IBAND.EQ.1) THEN
        JSTART = 16
      ELSE
        JSTART = 1
      ENDIF
      DO 20 J=JSTART,25

* For the j-th term of the trigonometric expansion, compute the angular
* argument angle of sine and cosine functions as a linear integer
* combination of the 6 fundamental arguments
        ANGLE = 0D0
        DO 10 I=1,6
          ANGLE = ANGLE + IARG(I,J) * ARG(I)
   10   CONTINUE
        ANGLE = DMOD( ANGLE, TWOPI )


* Compute contribution from the j-th term to the polar motion coordinates
        PM(1) = PM(1) + XS(J)*DSIN(ANGLE) + XC(J)*DCOS(ANGLE)
        PM(2) = PM(2) + YS(J)*DSIN(ANGLE) + YC(J)*DCOS(ANGLE)
   20 CONTINUE
      IF (IBAND.EQ.1) RETURN

* Add the secular term of the model
      PM(1) = PM(1) + XRATE * (RMJD-RMJD0) / 365.25D0
      PM(2) = PM(2) + YRATE * (RMJD-RMJD0) / 365.25D0

      RETURN

*  Finished.

*+----------------------------------------------------------------------
*
*  Copyright (C) 2008
*  IERS Conventions Center
*
*  ==================================
*  IERS Conventions Software License
*  ==================================
*
*  NOTICE TO USER:
*
*  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING TERMS AND CONDITIONS
*  WHICH APPLY TO ITS USE.
*
*  1. The Software is provided by the IERS Conventions Center ("the
*     Center").
*
*  2. Permission is granted to anyone to use the Software for any
*     purpose, including commercial applications, free of charge,
*     subject to the conditions and restrictions listed below.
*
*  3. You (the user) may adapt the Software and its algorithms for your
*     own purposes and you may distribute the resulting "derived work"
*     to others, provided that the derived work complies with the
*     following requirements:
*
*     a) Your work shall be clearly identified so that it cannot be
*        mistaken for IERS Conventions software and that it has been
*        neither distributed by nor endorsed by the Center.
*
*     b) Your work (including source code) must contain descriptions of
*        how the derived work is based upon and/or differs from the
*        original Software.
*
*     c) The name(s) of all modified routine(s) that you distribute
*        shall be changed.
* 
*     d) The origin of the IERS Conventions components of your derived
*        work must not be misrepresented; you must not claim that you
*        wrote the original Software.
*
*     e) The source code must be included for all routine(s) that you
*        distribute.  This notice must be reproduced intact in any
*        source distribution. 
*
*  4. In any published work produced by the user and which includes
*     results achieved by using the Software, you shall acknowledge
*     that the Software was used in obtaining those results.
*
*  5. The Software is provided to the user "as is" and the Center makes
*     no warranty as to its use or performance.   The Center does not
*     and cannot warrant the performance or results which the user may
*     obtain by using the Software.  The Center makes no warranties,
*     express or implied, as to non-infringement of third party rights,
*     merchantability, or fitness for any particular purpose.  In no
*     event will the Center be liable to the user for any consequential,
*     incidental, or special damages, including any lost profits or lost
*     savings, even if a Center representative has been advised of such
*     damages, or for any claim by any third party.
*
*  Correspondence concerning IERS Conventions software should be
*  addressed as follows:
*
*                     Gerard Petit
*     Internet email: gpetit[at]bipm.org
*     Postal address: IERS Conventions Center
*                     Time, frequency and gravimetry section, BIPM
*                     Pavillon de Breteuil
*                     92312 Sevres  FRANCE
*
*     or
*
*                     Brian Luzum
*     Internet email: brian.luzum[at]usno.navy.mil
*     Postal address: IERS Conventions Center
*                     Earth Orientation Department
*                     3450 Massachusetts Ave, NW
*                     Washington, DC 20392
*
*
*-----------------------------------------------------------------------
      END
