      SUBROUTINE clas_transf(tps,tuin,tut1,xp,yp,dpsi,deps,iclef1,xq)
C--------------------------------------------------------------------------C
C    Matrix transformation from the terrestrial reference system (TRS) to  C
C    the celestial reference system (CRS) at the date of observation       C
C                [CRS]= XQ [TRS]= prec * nut * uxy [TRS]                   C
C    Option (1) corresponding to the classical procedure,                  C
C      see IERS Conventions (1996), Chapter 5, p.20-25.                    C
C                                                                          C
C    (a-h,o-z) REAL*8                                                      C
C                                                                          C
C      INPUT :                                                             C
C       tps= Julian TDB date  (in Julian days)                             C
C       tuin= Julian UT1 date  (in Julian days)                            C
C       tut1= Universal time UT1, ((hx60)+m)x60)+sec +UT1-UTC, (in second) C 
C       xp, yp= "polar coordinates" of the CEP in the TRS  (in arcsecond)  C
C                see IERS Conventions (1996), p.21                         C 
C       dpsi, deps= celestial pole offsets (in arcsecond)                  C
C                   see IERS Conventions (1996), p. 22                     C
C       iclef1= 0  ==> complete equation of equinoxes                      C
C             <>0  ==> incomplete equation of equinoxes                    C
C                      (used in the analyses performed prior to 1997)      C
C      (see "IERS Conventions (1996)", IERS Technical Note 21, pp. 20-22)  C
C                                                                          C
C     OUTPUT :                                                             C
C       xq(3,3)= transformation matrix from the [TRS] to the [CRS]         C 
C                dimension(3,3), (in radian)                               C
C                                                                          C
C     COMMON / conversion /                                                C
C       tpi= 2*pi                                                          C
C       secrad= arc second in radian                                       C
C       hourad= second of time in radian                                   C
C                                                                          C
C     SUBROUTINE INTERFACE :                                               C
C       called subroutines : precession, wahr, prduxy, multm               C
C--------------------------------------------------------------------------C

       IMPLICIT NONE
       INTEGER iclef1
       REAL*8 tps,tuin,tut1,xp,yp,dpsi,deps,xq(3,3)
       REAL*8 tpi,secrad,hourad
       REAL*8 prec(3,3),nut(3,3),uxy(3,3)
       REAL*8 xx(3,3)

       COMMON/conversion/tpi,secrad,hourad


            ! Conversion factor
        tpi= 6.28318530717958648d0
        secrad= tpi / (360.d0*3600.d0)
        hourad= tpi / (24.d0*3600.d0)

        CALL precession(tps,prec)
        CALL wahr(tps,dpsi,deps,nut)
        CALL prduxy(tps,tuin,tut1,nut,xp,yp,uxy,iclef1)

        CALL multm(nut,uxy,xx,3,3)
        CALL multm(prec,xx,xq,3,3)
 
        RETURN
      END
C===========================================================================C


      SUBROUTINE prduxy(tps,tuin,tut1,nut,xp,yp,uxy,iclef1)
C-------------------------------------------------------------------------C
C      INPUT :                                                            C
C       tps= Julian TDB date  (in Julian days)                            C
C       tuin= Julian UT1 date  (in Julian days)                           C
C       tut1= Universal time UT1, ((h*60)+m)*60+sec+UT1-UTC, (in second)  C
C       nut(3,3)= nutation matrix  (in radian)                            C
C       xp, yp= "polar coordinates" of the CEP in the TRS  (in arcsecond) C
C                see IERS Conventions (1996), p.21                        C
C       iclef1= 0  ==> complete equation of equinoxes                     C
C             <>0  ==> incomplete equation of equinoxes                   C
C                see IERS Conventions (1996), p.21-22                     C
C                                                                         C
C     OUTPUT :                                                            C
C       uxy(3,3)= R3(-gst)*R1(yp)*R2(xp) matrix  (in radian)              C
C                                                                         C
C     COMMON / conversion /                                               C
C       tpi= 2*pi                                                         C
C       secrad= arc second in radian                                      C
C       hourad= second of time in radian                                  C
C                                                                         C
C     SUBROUTINE INTERFACE :                                              C
C        caller subroutine : clas_transf                                  C
C-------------------------------------------------------------------------C

       IMPLICIT NONE
       INTEGER iclef1
       REAL*8 tps,tuin,tut1,nut(3,3),xp,yp,uxy(3,3)
       REAL*8 tpi,secrad,hourad
       REAL*8 t,tu,h,omega,eq_equi
       REAL*8 ch,sh,cxp,sxp,cyp,syp

       COMMON/conversion/tpi,secrad,hourad


        t= (tps-2451545.d0) / 36525.d0
        tu= (tuin-2451545.d0) / 36525.d0

        IF (iclef1.EQ.0) THEN
          omega= (450160.28d0 - (1296000.d0*5.d0 + 482890.539)*t
     .         + 7.455d0*t**2 + 0.008d0*t**3)*secrad
          eq_equi= (0.00265d0*DSIN(omega) - 0.00009d0*DSIN(2.d0*omega))*secrad 
        ELSE
          eq_equi= 0.d0
        ENDIF

        h= (tut1 + (6.d0*60.d0 + 41.d0)*60.d0 + 50.54841d0
     .     + 8640184.812866d0*tu + 0.093104d0*tu**2 - 6.2D-6*tu**3)*hourad
     .     + DATAN(nut(1,2)/nut(1,1)) + eq_equi

        ch= DCOS(h)
        sh= DSIN(h)
        cxp= DCOS(xp*secrad)
        sxp= DSIN(xp*secrad)
        cyp= DCOS(yp*secrad)
        syp= DSIN(yp*secrad)

        uxy(1,1)= ch * cxp
        uxy(1,2)= ch * sxp * syp - sh * cyp
        uxy(1,3)= - ch * sxp * cyp - sh * syp
        uxy(2,1)= sh * cxp
        uxy(2,2)= sh * sxp * syp + ch * cyp
        uxy(2,3)= - sh * sxp * cyp + ch * syp
        uxy(3,1)= sxp
        uxy(3,2)= - cxp * syp
        uxy(3,3)= cxp * cyp

        RETURN
      END
C===========================================================================C


      SUBROUTINE precession(tps,prec)
C--------------------------------------------------------------C
C     INPUT :                                                  C
C      tps= Julian TDB date  (in Julian days)                  C
C                                                              C
C    OUTPUT :                                                  C
C      prec(3,3)= precession matrix  (in radian)               C
C                                                              C
C     COMMON / conversion /                                    C
C       tpi= 2*pi                                              C
C       secrad= arc second in radian                           C
C       hourad= second of time in radian                       C
C                                                              C
C     SUBROUTINE INTERFACE :                                   C
C        caller subroutine : clas_transf                       C
C--------------------------------------------------------------C

       IMPLICIT NONE
       REAL*8 tps,prec(3,3)
       REAL*8 tpi,secrad,hourad
       REAL*8 t,tt,ttt
       REAL*8 xia,theta,za
       REAL*8 cxia,sxia,ctheta,stheta,cza,sza

       COMMON/conversion/tpi,secrad,hourad


C     ************************************************************************
C      Numerical expression for the precession quantities xia,theta,za
C      IERS Conventions (1996), p. 22
C     ************************************************************************

        t= (tps-2451545.d0) / 36525.d0
        tt= t**2
        ttt= tt*t

        xia= (2306.2181d0*t + 0.30188d0*tt + 0.017998d0*ttt) * secrad
        theta= (2004.3109d0*t - 0.42665d0*tt - 0.041833d0*ttt) * secrad
        za= (2306.2181d0*t + 1.09468d0*tt + 0.018203d0*ttt) * secrad

        cxia= DCOS(xia)
        sxia= DSIN(xia)
        ctheta= DCOS(theta)
        stheta= DSIN(theta)
        cza= DCOS(za)
        sza= DSIN(za)

        prec(1,1)= cxia * ctheta * cza - sxia * sza 
        prec(1,2)= cxia * ctheta * sza + sxia * cza
        prec(1,3)= cxia * stheta
        prec(2,1)= - sxia * ctheta * cza - cxia * sza
        prec(2,2)= - sxia * ctheta * sza + cxia * cza
        prec(2,3)= - sxia * stheta
        prec(3,1)= - stheta * cza
        prec(3,2)= - stheta * sza
        prec(3,3)= ctheta

        RETURN
      END
C===========================================================================C


      SUBROUTINE wahr(tps,dpsi,deps,nut)
C--------------------------------------------------------------C
C     INPUT :                                                  C
C      tps= Julian TDB date  (in Julian days)                  C 
C      dpsi, deps= celestial pole offsets (in arcsecond)       C
C                   see IERS Conventions (1996) p. 22-25       C
C                                                              C
C    OUTPUT :                                                  C
C      nut(3,3)= nutation matrix  (in radian)                  C
C                                                              C
C     COMMON / conversion /                                    C
C       tpi= 2*pi                                              C
C       secrad= arc second in radian                           C
C       hourad= second of time in radian                       C
C                                                              C
C     SUBROUTINE INTERFACE :                                   C
C        caller subroutine : clas_transf                       C
C--------------------------------------------------------------C 

       IMPLICIT NONE
       INTEGER marg(9,106),matarg(9,106)
       INTEGER i,j
       REAL*8 tps,dpsi,deps,nut(3,3)
       REAL*8 tpi,secrad,hourad
       REAL*8 t,tt,ttt,r
       REAL*8 dangle(5),angle,deltapsi,deltaeps,epsa,eps
       REAL*8 cpsi,spsi,ceps,seps,cepsa,sepsa

       COMMON/conversion/tpi,secrad,hourad

       EQUIVALENCE (marg(1,1),matarg(1,1))

C      *********************************************
C       Table 5.1, IERS Conventions (1996), p. 23-25
C      *********************************************

        DATA marg /
     @     0,  0,  0,  0,  1, -171996, -1742, 92025,   89,
     @     0,  0,  2, -2,  2,  -13187,   -16,  5736,  -31,
     @     0,  0,  2,  0,  2,   -2274,    -2,   977,   -5,
     @     0,  0,  0,  0,  2,    2062,     2,  -895,    5,
     @     0,  1,  0,  0,  0,    1426,   -34,    54,   -1,
     @     1,  0,  0,  0,  0,     712,     1,    -7,    0,
     @     0,  1,  2, -2,  2,    -517,    12,   224,   -6,
     @     0,  0,  2,  0,  1,    -386,    -4,   200,    0,
     @     1,  0,  2,  0,  2,    -301,     0,   129,   -1,
     @     0, -1,  2, -2,  2,     217,    -5,   -95,    3,
     @     1,  0,  0, -2,  0,    -158,     0,    -1,    0,
     @     0,  0,  2, -2,  1,     129,     1,   -70,    0,
     @    -1,  0,  2,  0,  2,     123,     0,   -53,    0,
     @     1,  0,  0,  0,  1,      63,     1,   -33,    0,
     @     0,  0,  0,  2,  0,      63,     0,    -2,    0,
     @    -1,  0,  2,  2,  2,     -59,     0,    26,    0,
     @    -1,  0,  0,  0,  1,     -58,    -1,    32,    0,
     @     1,  0,  2,  0,  1,     -51,     0,    27,    0,
     @     2,  0,  0, -2,  0,      48,     0,     1,    0,
     @    -2,  0,  2,  0,  1,      46,     0,   -24,    0,
     @     0,  0,  2,  2,  2,     -38,     0,    16,    0,
     @     2,  0,  2,  0,  2,     -31,     0,    13,    0,
     @     2,  0,  0,  0,  0,      29,     0,    -1,    0,
     @     1,  0,  2, -2,  2,      29,     0,   -12,    0,
     @     0,  0,  2,  0,  0,      26,     0,    -1,    0,
     @     0,  0,  2, -2,  0,     -22,     0,     0,    0,
     @    -1,  0,  2,  0,  1,      21,     0,   -10,    0,
     @     0,  2,  0,  0,  0,      17,    -1,     0,    0,
     @     0,  2,  2, -2,  2,     -16,     1,     7,    0,
     @    -1,  0,  0,  2,  1,      16,     0,    -8,    0,
     @     0,  1,  0,  0,  1,     -15,     0,     9,    0,
     @     1,  0,  0, -2,  1,     -13,     0,     7,    0,
     @     0, -1,  0,  0,  1,     -12,     0,     6,    0,
     @     2,  0, -2,  0,  0,      11,     0,     0,    0,
     @    -1,  0,  2,  2,  1,     -10,     0,     5,    0,
     @     1,  0,  2,  2,  2,      -8,     0,     3,    0,
     @     0, -1,  2,  0,  2,      -7,     0,     3,    0,
     @     0,  0,  2,  2,  1,      -7,     0,     3,    0,
     @     1,  1,  0, -2,  0,      -7,     0,     0,    0,
     @     0,  1,  2,  0,  2,       7,     0,    -3,    0,
     @    -2,  0,  0,  2,  1,      -6,     0,     3,    0,
     @     0,  0,  0,  2,  1,      -6,     0,     3,    0,
     @     2,  0,  2, -2,  2,       6,     0,    -3,    0,
     @     1,  0,  0,  2,  0,       6,     0,     0,    0,
     @     1,  0,  2, -2,  1,       6,     0,    -3,    0,
     @     0,  0,  0, -2,  1,      -5,     0,     3,    0,
     @     0, -1,  2, -2,  1,      -5,     0,     3,    0,
     @     2,  0,  2,  0,  1,      -5,     0,     3,    0,
     @     1, -1,  0,  0,  0,       5,     0,     0,    0,
     @     1,  0,  0, -1,  0,      -4,     0,     0,    0,
     @     0,  0,  0,  1,  0,      -4,     0,     0,    0,
     @     0,  1,  0, -2,  0,      -4,     0,     0,    0,
     @     1,  0, -2,  0,  0,       4,     0,     0,    0,
     @     2,  0,  0, -2,  1,       4,     0,    -2,    0,
     @     0,  1,  2, -2,  1,       4,     0,    -2,    0,
     @     1,  1,  0,  0,  0,      -3,     0,     0,    0,
     @     1, -1,  0, -1,  0,      -3,     0,     0,    0,
     @    -1, -1,  2,  2,  2,      -3,     0,     1,    0,
     @     0, -1,  2,  2,  2,      -3,     0,     1,    0,
     @     1, -1,  2,  0,  2,      -3,     0,     1,    0,
     @     3,  0,  2,  0,  2,      -3,     0,     1,    0,
     @    -2,  0,  2,  0,  2,      -3,     0,     1,    0,
     @     1,  0,  2,  0,  0,       3,     0,     0,    0,
     @    -1,  0,  2,  4,  2,      -2,     0,     1,    0,
     @     1,  0,  0,  0,  2,      -2,     0,     1,    0,
     @    -1,  0,  2, -2,  1,      -2,     0,     1,    0,
     @     0, -2,  2, -2,  1,      -2,     0,     1,    0,
     @    -2,  0,  0,  0,  1,      -2,     0,     1,    0,
     @     2,  0,  0,  0,  1,       2,     0,    -1,    0,
     @     3,  0,  0,  0,  0,       2,     0,     0,    0,
     @     1,  1,  2,  0,  2,       2,     0,    -1,    0,
     @     0,  0,  2,  1,  2,       2,     0,    -1,    0,
     @     1,  0,  0,  2,  1,      -1,     0,     0,    0,
     @     1,  0,  2,  2,  1,      -1,     0,     1,    0,
     @     1,  1,  0, -2,  1,      -1,     0,     0,    0,
     @     0,  1,  0,  2,  0,      -1,     0,     0,    0,
     @     0,  1,  2, -2,  0,      -1,     0,     0,    0,
     @     0,  1, -2,  2,  0,      -1,     0,     0,    0,
     @     1,  0, -2,  2,  0,      -1,     0,     0,    0,
     @     1,  0, -2, -2,  0,      -1,     0,     0,    0,
     @     1,  0,  2, -2,  0,      -1,     0,     0,    0,
     @     1,  0,  0, -4,  0,      -1,     0,     0,    0,
     @     2,  0,  0, -4,  0,      -1,     0,     0,    0,
     @     0,  0,  2,  4,  2,      -1,     0,     0,    0,
     @     0,  0,  2, -1,  2,      -1,     0,     0,    0,
     @    -2,  0,  2,  4,  2,      -1,     0,     1,    0,
     @     2,  0,  2,  2,  2,      -1,     0,     0,    0,
     @     0, -1,  2,  0,  1,      -1,     0,     0,    0,
     @     0,  0, -2,  0,  1,      -1,     0,     0,    0,
     @     0,  0,  4, -2,  2,       1,     0,     0,    0,
     @     0,  1,  0,  0,  2,       1,     0,     0,    0,
     @     1,  1,  2, -2,  2,       1,     0,    -1,    0,
     @     3,  0,  2, -2,  2,       1,     0,     0,    0,
     @    -2,  0,  2,  2,  2,       1,     0,    -1,    0,
     @    -1,  0,  0,  0,  2,       1,     0,    -1,    0,
     @     0,  0, -2,  2,  1,       1,     0,     0,    0,
     @     0,  1,  2,  0,  1,       1,     0,     0,    0,
     @    -1,  0,  4,  0,  2,       1,     0,     0,    0,
     @     2,  1,  0, -2,  0,       1,     0,     0,    0,
     @     2,  0,  0,  2,  0,       1,     0,     0,    0,
     @     2,  0,  2, -2,  1,       1,     0,    -1,    0,
     @     2,  0, -2,  0,  1,       1,     0,     0,    0,
     @     1, -1,  0, -2,  0,       1,     0,     0,    0,
     @    -1,  0,  0,  1,  1,       1,     0,     0,    0,
     @    -1, -1,  0,  2,  1,       1,     0,     0,    0,
     @     0,  1,  0,  1,  0,       1,     0,     0,    0/


C     ************************************************************************
C      Fundamental arguments of the IAU 1980 Theory of Nutation,
C      IERS Conventions (1996), p. 23 
C        dangle(1)= l        mean anomaly of the Moon
C        dangle(2)= l'       mean anomaly of the Sun
C        dangle(3)= F        mean longitude of the Moon - omega
C        dangle(4)= D        mean elongation of the Moon from the sun
C        dangle(5)= omega    mean longitude of the ascending node of the Moon
C        r= 360 degrees = 1 296 000.0"
C     ************************************************************************

        t= (tps-2451545.d0) / 36525.d0
        tt= t**2
        ttt= tt*t
        r= 1296000.d0

        dangle(1)= 485866.733d0 + (1325.d0*r + 715922.633d0)*t +
     .            31.310d0*tt + 0.064d0*ttt
        dangle(2)= 1287099.804d0 + (99.d0*r + 1292581.224d0)*t -
     .            0.577d0*tt - 0.012d0*ttt
        dangle(3)= 335778.877d0 + (1342.d0*r + 295263.137d0)*t -
     .            13.257d0*tt + 0.011d0*ttt
        dangle(4)= 1072261.307d0 + (1236.d0*r + 1105601.328d0)*t -
     .            6.891d0*tt + 0.019d0*ttt
        dangle(5)= 450160.280d0 - (5.d0*r + 482890.539d0)*t +
     .            7.455d0*tt + 0.008d0*ttt
 
           ! conversion " in radian
        DO i=1,5
          dangle(i)= dangle(i)*secrad
        ENDDO

        deltapsi=0.d0
        deltaeps=0.d0

        DO i=1,106

          angle=0.d0
          DO j=1,5
            angle= angle + DFLOAT(matarg(j,i))*dangle(j)
          ENDDO
          angle= dmod(angle,tpi)

          deltapsi= deltapsi + ( DFLOAT(matarg(6,i)) + 
     .             DFLOAT(matarg(7,i))/10.d0*t ) * DSIN(angle) / 1.d4
          deltaeps= deltaeps + ( DFLOAT(matarg(8,i)) +
     .             DFLOAT(matarg(9,i))/10.d0*t ) * DCOS(angle) / 1.d4
        ENDDO

        deltapsi= (deltapsi + dpsi) * secrad
        deltaeps= (deltaeps + deps) * secrad

        epsa= (84381.448d0 - 46.815d0*t - 0.00059d0*tt +
     .       0.001813d0*ttt) * secrad
        eps= epsa + deltaeps

        cpsi= DCOS(deltapsi)
        spsi= DSIN(deltapsi)
        ceps= DCOS(eps)
        seps= DSIN(eps)
        cepsa= DCOS(epsa)
        sepsa= DSIN(epsa)

        nut(1,1)= cpsi
        nut(1,2)= ceps * spsi
        nut(1,3)= seps * spsi
        nut(2,1)= - cepsa * spsi
        nut(2,2)= cepsa * ceps * cpsi + sepsa * seps 
        nut(2,3)= cepsa * seps * cpsi - sepsa * ceps
        nut(3,1)= - sepsa * spsi
        nut(3,2)= sepsa * ceps * cpsi - cepsa * seps
        nut(3,3)= sepsa * seps * cpsi + cepsa * ceps

        RETURN
      END
C===========================================================================C


      SUBROUTINE multm(x1,x2,xr,ni,nj)
C---------------------------------------------------------------C
C  Multiplication of matrices                                   C
C    xr(ni,nj) = x1(ni,ni) * x2(ni,nj)                          C
C                                                               C
C     (a-h,o-z) REAL*8                                          C
C                                                               C
C     SUBROUTINE INTERFACE :                                    C
C        caller subroutine : clas_transf                        C
C---------------------------------------------------------------C

       IMPLICIT NONE
       INTEGER ni,nj
       INTEGER i,j,k
       REAL*8 x1(ni,ni),x2(ni,nj),xr(ni,nj)
       REAL*8 cal


        DO i=1,ni
          DO j=1,nj
            cal= 0.d0
            DO k=1,ni
              cal= cal + x1(i,k) * x2(k,j)
            ENDDO
            xr(i,j)= cal
          ENDDO
        ENDDO

        RETURN
      END
C===========================================================================C

