        ! Calcul du quaternion de rotation a partir des suborutines SOFA 
        ! Last modification: June 2023
	implicit none
	
	integer n,ndate,i,j,kdate,m,l, ii,nj
	parameter (n=20000)
	double precision mjd_v,xp_v,yp_v,UT1_UTC_v,lod_v,dX_v,dY_v,sxp_v,syp_v,sUT1_UTC_v,sLOD_v,sdX_v,sdY_v 		
        double precision mmjd(n),xp(n),yp(n),UT1_UTC(n),lod(n),dX(n),dY(n),UT1_TAI(n)
	double precision mjd_b(21),xp_b(21),yp_b(21),UT1_UTC_b(21),lod_b(21),dX_b(21),dY_b(21),UT1_TAI_b(21)
	double precision xpi,ypi,xpi_dt,ypi_dt,UT1_UTCi,lodi,dXi,dYi,UT1_TAIi

	double precision spline,tervon,interplin,mjd_to_TAIUTC,iau_S06, iau_SP00, iau_ERA00, iau_S00, tartuf, iau_obl06
        double precision pi,degrad,secrad,masrad
	double precision epsa,eps,psi,omega

	external iau_S06A,iau_SP00, iau_S00, iau_obl06

	double precision sCIO,sTIO,theta,theta1,X,Y,Z,X_mod,Y_mod
	double precision cor_x ,cor_y, cor_ut1, cor_lod, A, B, C, D 
	double precision MQ(3,3), XC(3), XC1(3), RP(3,3), RN(3,3), RPN(3,3)
	double precision QUA(4), QUAI(4), QUAS(4), DQUA(4), QUAC(4), QUAW(4), QPN(4), QS(4)
	double precision s,s1,s2,h,lat,lon,fac, x1, y1, z1

	double precision mjd_tai, d_mjd_tai
	double precision mjd_ut1, d_mjd_ut1, mjd_ut1_0h, jd, jd0, d_mjd
	double precision mjd_tt,  d_mjd_tt 
	double precision mjd_utc, d_mjd_utc
	
	double precision mjd1,    d_mjd1, mjd1_tai, mjd1_utc
	double precision mjd2,    d_mjd2, mjd2_tai, mjd2_utc
	
	double precision jd_ut1,  jd_ut1_0h, d_jd_ut1
	double precision jd_utc,  jd_utc_0h, d_jd_utc
	double precision jd_tai,  jd_tai_0h, d_jd_tai
	double precision jd_tt,   jd_tt_0h,  d_jd_tt 
	
	double precision dt, mjd_tai_dt, mjd_ut1_dt, mjd_tt_dt, jd_tt_0h_dt, d_jd_tt_dt, jd_ut1_0h_dt, d_jd_ut1_dt 
	double precision wx, wy, wz, wx1, wy1, wz1
	  
	integer*4 option,pm,ut,dP,hf, optq,optp
	integer an1,mois1,jour1,nh1,nm1,ns1,an2,mois2,jour2,nh2,nm2,ns2,an,mois,jour,nh,nm,ns,ncom 
        character*100 clat,clon,ch,coption,cpm,cut,cdP,chf, cjour1,cmois1,can1,ch1,cm1,cs1,cjour2,cmois2,can2,ch2,cm2,cs2,coptq,
     &    coptp,fileinput,kom,scale,sign
	integer IHMSF(4)
	COMMON/conversion/pi,secrad,degrad

	real :: start, finish
    	call cpu_time(start)   	

        ! Conversion factor
        pi     = 6.28318530717958648d0
        secrad = pi / (360.d0*3600.d0)
        degrad = pi / 360.d0	
	masrad = secrad / 1d3
	fileinput='/home/ftp/iers/series/opa/eopc04_IAU2000' !  serie 9000 depuis 1962 + prediction

        call getarg(1,can1)	    
	
	if (Len_trim(can1).eq.0) then
	write(6,*) '<<   Transf. de coordonnees du repere terrestre international Oxyz >>'
	write(6,*) '<<   au repere celeste international OXYZ et parametres derives    >>'
	write(6,*) 'Usage : rotation                     '
        write(6,*) '     [an1][mois1][jour1][h1][mn1][s1]'
        write(6,*) '     [an2][mois2][jour2][h2][mn2][s2]'
        write(6,*) '     [pas][option][pm][ut1][dP][hf] [lon][lat][optq][optp]'
	write(6,*) '[an1] .... 1ere date : annee  '
        write(6,*) '[mois1] .. 1ere date : mois   '
        write(6,*) '[jour1]... 1ere date : jour   '       
        write(6,*) '[h1]  .... 1ere date : heure  '	  
        write(6,*) '[mn1] .... 1ere date : minute '	  
        write(6,*) '[s1]  .... 1ere date : seconde'	  
        write(6,*) '[an2] .... date fin  : annee  '
        write(6,*) '[mois2] .. date fin  : mois   '
        write(6,*) '[jour2] .. date fin  : jour   '	 
        write(6,*) '[h2]  .... date fin  : heure  '    
        write(6,*) '[mn2] .... date fin  : minute '    
        write(6,*) '[s2]  .... date fin  : seconde'    
        write(6,*) '[pas]  ... pas temporel en seconde'       
	write(6,*) '[option] 1 matrice de rotation / quaternion de rotation'
	write(6,*) '	     2 axe Oz / axe (lat, long) projete sur le plan eq. J2000.0 OXY (arcsecond)'
	write(6,*) '	     3 vec. inst. rotation dans ITRF (nanorad/s)',
     &                        'mjd_tai,an,mois,jour,nh,nm,s,wx,wy,wz (TRF)/ wx1,wy1,wz1 (repere local)/Omega xp,Omega yp'
	write(6,*) '	     4 vec. inst. rotation dans ICRF (nanorad/s)'			    
	write(6,*) '	     5 precession axe de figure : Terre réelle & Terre Rigide Kinoshita-Souchay (1995)', 
     &	                      'mjd_tai,an,mois,jour,nh,nm,s, X, Y (Terre réelle), X1,Y1 (Terre rigide)' 				       
	write(6,*) '[pm] .... mouvement du pole   :  oui :1 /non :0'
	write(6,*) '[ut] .... UT1 - UTC 	  :  oui :1 /non :0'
	write(6,*) '[dP] .... ecarts pole celeste :  oui :1 /non :0'
	write(6,*) '[hf] .... corr. de marees diurnes et semi-diurnes :  oui :1 /non :0'
	write(6,*) '[lon] ... longitude 	 '
	write(6,*) '[lat] ... latitude  	 '
        write(6,*) '[optq] .. 0: ITRF-> ICRF / 1: CIF (Celestial Interm. Frame)-> ICRF / 2:TIF(Terretrial Interm. Frame)->CIF',
     &                       '/ 3: ICRF --> TIF / 4: ITRF --> CIF' 
	write(6,*) '[optp] .. 1: quaternion approximé / 2: quaternion exact'
	write(6,*) 'Sortie dans trans.res' 
	stop     
	endif

        call getarg(2,cmois1)
        call getarg(3,cjour1)   
        call getarg(4,ch1)
        call getarg(5,cm1)
        call getarg(6,cs1)        
        call getarg(7,can2)
        call getarg(8,cmois2)
        call getarg(9,cjour2)
        call getarg(10,ch2)
        call getarg(11,cm2)
        call getarg(12,cs2)
        call getarg(13,ch)
        call getarg(14,coption)
        call getarg(15,cpm)
        call getarg(16,cut)        
	call getarg(17,cdP)
	call getarg(18,chf)
	call getarg(19,clon)
	call getarg(20,clat)
	call getarg(21,coptq)
	call getarg(22,coptp)
	
	read(cjour1,*)  jour1
	read(cmois1,*)  mois1
	read(can1,*)    an1
	read(ch1,*)     nh1
	read(cm1,*)     nm1
	read(cs1,*)     s1 		
	read(cjour2,*)  jour2
	read(cmois2,*)  mois2
	read(can2,*)    an2
	read(ch2,*)     nh2
	read(cm2,*)     nm2
	read(cs2,*)     s2
	read(ch,*) h
	read(coption,*) option
	read(cpm,*) pm
	read(cut,*) ut
	read(cdP,*) dP
	read(chf,*) hf
	read(clon,*) lon
	read(clat,*) lat
	read(coptq,*) optq
	read(coptp,*) optp
	
	lon = lon * degrad
	lat = lat * degrad

	
	
! conversion en jour julien 
	
	ns1 = dint(s1)
	ns2 = dint(s2)
    		
        ! call DATEDJ (jour1,mois1,an1,jd,mjd1) 
	
	call iau_CAL2JD( an1, mois1, jour1, jd0, mjd1, j )
	
	d_mjd1   = dble(nh1)/24.D0 + dble(nm1)/60.D0/24.D0+ s1/86400.D0   ! fraction de jour de la date de debut
	mjd1_utc = mjd1 + d_mjd1
	mjd1_tai = mjd1_utc + mjd_to_TAIUTC(mjd1_utc)/86400.d0            ! conversion en echelle de temps TAI
	
	call iau_CAL2JD( an2, mois2, jour2, jd0, mjd2, j )
	
	d_mjd2    = dble(nh2)/24.D0 + dble(nm2)/60.D0/24.D0 + s2/86400.D0 
	mjd2_utc = mjd2 + d_mjd2
	mjd2_tai = mjd2_utc + mjd_to_TAIUTC(mjd2_utc)/86400.d0            ! conversion en echelle de temps TAI	
	
	i = 0     		
	ncom = 13 + abs(an1 - 1962) * 365  ! le debut de la serie 9000 est 1962
      
! -------------------------------------
! lecture et selection des EOP observes
! -------------------------------------

	open(10,file=fileinput, status='old')
	
	do j =1,ncom
	read(10,81,end=9) kom
	enddo	
81 	format(a100)
9	continue
	
	do j=1,100000	
		read(10,*,end=10) an,mois,jour,mjd_v,xp_v,yp_v,UT1_UTC_v,lod_v,dX_v,dY_v,sxp_v,syp_v,sUT1_UTC_v,sLOD_v,sdX_v,sdY_v
	      
		if ((mjd_v.ge.mjd1-10.d0).and.(mjd_v.le.mjd2+10.d0)) then
			i = i + 1
			mmjd(i)    = mjd_v 
			xp(i)      = xp_v
			yp(i)      = yp_v
			UT1_TAI(i) = UT1_UTC_v - mjd_to_TAIUTC(mjd_v)! calcul de UT1-TAI pour faire  l'interpolation sur une echelle de temps
			!UT1_TAI(i) = UT1_UTC_v - tartuf(mjd_v)! calcul de UT1-TAI pour faire  l'interpolation sur une echelle de temps
			lod(i)     = lod_v
			dX(i)      = dX_v
			dY(i)      = dY_v
			
		endif	
	enddo
		
10 	continue
	close(10)
	
	kdate = i  ! nombre de dates selectionnees		

 	! boucle temporelle - interpolation des EOP avec le pas h  

	ndate = dint((mjd2_tai-mjd1_tai)*86400.d0/h)+1
		
	do i=1,ndate	
	       
		d_mjd   = dble(i-1) * h / 86400.d0               ! fraction de jour
		mjd_tai = mjd1_tai + d_mjd                       ! en TAI				
		mjd_utc = mjd1_utc + d_mjd                       ! en UTC	

		d_mjd_tai = mjd_tai - dint(mjd_tai)
		d_mjd_utc = mjd_utc - dint(mjd_utc)
		
		jd_tai_0h = dint(mjd_tai) + 2400000.5 
		
		jd_utc_0h = dint(mjd_utc) + 2400000.5
			
		if (kdate.ne.0) then ! EOP dans la periode selectionnee

			j = 0
	        	do l = 1,kdate	 
				if ( (mmjd(l).ge.mjd_utc - 10.0).and.(mmjd(l).le.mjd_utc + 10.0) ) then  ! selection de 10 points avant et apres la date centrale
					j = j+1	
					mjd_b(j)     = mmjd(l) 
					xp_b(j)      = xp(l)
					yp_b(j)      = yp(l) 
					UT1_TAI_b(j) = UT1_TAI(l)
					lod_b(j)     = lod(l) 
					dX_b(j)      = dX(l) 
					dY_b(j)      = dY(l) 
				endif
					 	  
			enddo	
			nj =j 
	     		if (nj.lt.3) then 
				xpi      = 0.d0
				ypi      = 0.d0 
				UT1_UTCi = 0.d0
				lodi     = 0.d0 
				dXi      = 0.d0
				dYi      = 0.d0		
			else
				! -----------------------  
				! interpolation a mjd_tai 
				! ----------------------- 
	        		xpi	 = spline(mjd_b,xp_b,nj,mjd_tai)
	        		ypi	 = spline(mjd_b,yp_b,nj,mjd_tai)
	       	 		UT1_TAIi = spline(mjd_b,UT1_TAI_b,nj,mjd_tai) 		
				UT1_UTCi = UT1_TAIi + mjd_to_TAIUTC(mjd_utc)  ! on reforme UT1 - UTC	
				!UT1_UTCi = UT1_TAIi + tartuf(mjd_utc)        ! on reforme UT1 - UTC	
				lodi	 = spline(mjd_b,lod_b,nj,mjd_tai)
	        		dXi	 = spline(mjd_b,dX_b,nj,mjd_tai)
	        		dYi	 = spline(mjd_b,dY_b,nj,mjd_tai)				       
			
			      
			endif
		 	      
        	else ! aucun EOP dans la periode selectionnee
	 		xpi      = 0.d0
			ypi      = 0.d0 
			UT1_UTCi = 0.d0
			lodi     = 0.d0 
			dXi      = 0.d0
			dYi      = 0.d0	
         	endif 						      
                  
	      
		if (hf.eq.1) then ! corrections marees oceaniques		        
			call PMUT1_OCEANS(mjd_tai,cor_x,cor_y,cor_ut1,cor_lod)
			xpi	 = xpi      + cor_x 	 
			ypi      = ypi	    + cor_y 
			UT1_UTCi = UT1_UTCi + cor_ut1		
			lodi	 = lodi     + cor_lod 
		endif
	
		if (pm.eq.0) then
			xpi      = 0.d0
			ypi      = 0.d0 
		endif
	
		if (ut.eq.0) then
			UT1_UTCi = 0.d0
			lodi     = 0.d0
	
		endif
			
		if (dP.eq.0) then
			dXi      = 0.d0
			dYi      = 0.d0
		endif 
		
		! write(777,*) mjd_tai, dXi,dYi, dP, ut 
! vitesse de rotation axiale rad/s
		omega = 7.2921151467064d-5*(1.d0 - lodi/86400.d0)
	
			
	      	!write(6,*) mjd_tai,lodi,omega		
		mjd_ut1    = mjd_utc + UT1_UTCi / 86400.d0 ! date julienne modifiee UT1 	  
		mjd_ut1_0h = dint(mjd_UT1)                 ! date julienne modifiee UT1 a 0h		  

		jd_ut1     = mjd_ut1    + 2400000.5d0      ! date julienne UT1  	  
		jd_ut1_0h  = mjd_ut1_0h + 2400000.5d0      ! date julienne UT1   0hUT1
 		d_jd_ut1   = mjd_ut1 - mjd_ut1_0h          ! fraction de jour UT1 a partir de Oh UT1 = mjd_UT1 - mjd_UT1_0h
		
		mjd_tt     = mjd_tai      + 32.184/86400.d0  		
		jd_tt_0h   = dint(mjd_tt) + 2400000.5d0  ! date julienne TT  a OhTT 	
		d_jd_tt    = mjd_tt - dint(mjd_tt) 
			
		call  iau_D2DTF ('UTC', 9, jd_utc_0h, d_mjd_utc, an, mois, jour, IHMSF, j  )
		 
	                                     
! ###########################
! coordonnees célestes du CIP
! ###########################
      		call iau_XY06 ( jd_tt_0h, d_jd_tt, X_mod, Y_mod)   
		
		X = X_mod + dXi * secrad  
		Y = Y_mod + dYi * secrad 
		Z = dsqrt( 1.d0 - X**2 - Y**2) 
		
! ################################
! coordonnees du CIP dans la Terre
! ################################
		xpi = xpi * secrad  
		ypi = ypi * secrad  

!		write(77,777) mjd_utc, mjd_tai, mjd_ut1, xpi, ypi, ut1_utci, dXi, dYi, X_mod, Y_mod
!777             format(3f20.9,7f20.15)	
  		
! #########################
! origine non tournante CIO
! #########################
		!write(6,*) 'jd_UT1',jd_UT1
		sCIO = iau_S00(jd_tt_0h, d_jd_tt, X,Y )     ! rad 	
					
! #########################
! origine non tournante TIO
! #########################
		sTIO = iau_SP00(jd_tt_0h, d_jd_tt )          ! rad
		
! #############################
! angle de rotation de la Terre
! #############################	        	   
      		theta = iau_ERA00(jd_UT1_0h, d_jd_UT1)  ! rad		 	
      			
! ######################################
! quaternion d'orientation et sa derivee
! ######################################
	        			
		if (option.le.4) then
			 
			! write(6,*) xpi,ypi,dXi,dYi
			call quaternion(X,Y,Z,sCIO,theta,sTIO,xpi,ypi,QUA,MQ,optq,optp)
							
			if (option.eq.1) write(6,1000) mjd_tai,an,mois,jour, IHMSF(1),IHMSF(2),IHMSF(3),IHMSF(4), 
     &                  MQ(1,1), MQ(1,2), MQ(1,3),MQ(2,1),MQ(2,2),MQ(2,3),MQ(3,1),MQ(3,2), MQ(3,3),QUA(1),QUA(2),QUA(3),QUA(4)		
1000			format(f22.13,i6,2(i3),3(i3),'.',I9.9, 9f18.14, 4f18.14)
 		               
	       		if (option.eq.2) then
				! coordonnees de l'axe Oz dans le repere celeste	
				x1 = MQ(1,3) / secrad   ! conversion en seconde
				y1 = MQ(2,3) / secrad   ! conversion en seconde
				z1 = MQ(3,3) / secrad   ! conversion en seconde
							
				write(6,1010) mjd_tai,an,mois,jour, IHMSF(1),IHMSF(2),IHMSF(3),IHMSF(4), x1, y1, z1
1010				format(f22.13,i6,2(i3),3(i3),'.',I9.9, 3f18.6)

			endif

			if ((option.eq.3).or.(option.eq.4)) then

				dt = 5.d0  ! pas de la derivation temporelle en seconde 
				
				do ii=1, 2
					if (ii.eq.1) dt = 0.0 
					if (ii.eq.2) dt = 5.d0 
					
					mjd_tai_dt   = mjd_tai + dt / 86400.d0
					
					mjd_ut1_dt   = mjd_ut1 + dt / 86400.d0 
					mjd_tt_dt    = mjd_tt  + dt / 86400.d0 
				
					jd_tt_0h_dt  = dint(mjd_tt_dt)  + 2400000.5d0 ! date julienne TT  a OhTT 		  
					d_jd_tt_dt   = mjd_tt_dt  - dint(mjd_tt_dt)   ! fraction date julienne TT  a partir de 0hTT  
				  						
					jd_ut1_0h_dt = dint(mjd_ut1_dt) + 2400000.5d0 ! date julienne UT1 a OhUT1		
					d_jd_ut1_dt  = mjd_ut1_dt - dint(mjd_ut1_dt)  ! fraction date julienne UT1 a partir de 0hUT1  
			
       	   				if ((pm.ne.0).and.(kdate.ne.0)) then
						xpi_dt = spline(mjd_b,xp_b,nj,mjd_tai_dt) 
	        				ypi_dt = spline(mjd_b,yp_b,nj,mjd_tai_dt)  
						if (hf.eq.1) then ! corrections marees oceaniques		        
							call PMUT1_OCEANS(mjd_tai,cor_x,cor_y,cor_ut1,cor_lod)
								xpi_dt = xpi_dt + cor_x 	 
								ypi_dt = ypi_dt + cor_y 
							endif                             			       	       
					else
						xpi_dt=0.d0
						ypi_dt=0.d0
					endif
					
					xpi_dt = xpi_dt * secrad
					ypi_dt = ypi_dt * secrad
					 
      					!write(6,*) 'couc', nj,xpi,ypi
					
					theta = iau_ERA00(jd_UT1_0h_dt, d_jd_UT1_dt)  ! rad	
					
					if (ii.eq.1)                                  theta1 = theta	
					if ( (ii.eq.2).and.(theta - theta1.le.0.d0) ) theta  = theta + pi   
						
					call iau_XY06 ( jd_tt_0h_dt, d_jd_tt_dt, X_mod, Y_mod)   
		
					X = X_mod + dXi * secrad  
					Y = Y_mod + dYi * secrad 
					Z = dsqrt( 1.d0 - X**2 - Y**2) 
					
					if (ii.eq.1) call quaternion(X,Y,Z,sCIO,theta,sTIO,xpi_dt,ypi_dt,QUAI,MQ,optq,optp)
					if (ii.eq.2) call quaternion(X,Y,Z,sCIO,theta,sTIO,xpi_dt,ypi_dt,QUAS,MQ,optq,optp)				  
					
				
				enddo	
			
		        	do ii=1,4                           
						DQUA(ii) = ( QUAS(ii) - QUAI(ii) ) /  dt 											
				enddo
				
	
				call qua_con(QUA,QUAC)
				
				 	
				if (option.eq.3) then ! vecteur instantane de rotation dans le repere terrestre 

! ############################################################
! Quaternion de rotation instantannée dans le repere terrestre 
! q_{\omega} = 2 \overline{q} \dot{q}
! ############################################################
					
					call qua_mul(QUAC,DQUA,QUAW)
				
					call qua_mul_sca(QUAW,2.d0)
					
					wx = QUAW(2) * 1d9  ! nanorad/s
					wy = QUAW(3) * 1d9
					wz = omega   * 1d9
					
					!write(6,*) mjd_tai, wx, wy, wz 
					wx1 =  dsin(lat) * dcos(lon) * wx + dsin(lat) * dsin(lon) * wy - dcos(lat) * wz
					wy1 = -dsin(lon) * wx + dcos(lon) * wy 
					wz1 =  dcos(lat) * dcos(lon) * wx + dcos(lat) * dsin(lon) * wy + dsin(lat) * wz
	
					x1 =  xpi * 7.292115d-5 * 1d9 ! composantes associees aux coordonnées du pôle observees nanorad/s
					y1 = -ypi * 7.292115d-5 * 1d9
		
					write(6,1020) mjd_tai,an,mois,jour,IHMSF(1),IHMSF(2),IHMSF(3),IHMSF(4),
     &                                  wx, wy, wz, wx1, wy1, wz1, x1, y1 
1020				        format(f22.13,i6,2(i3),3(i3),'.',i9.9, 8f18.7)

			
				endif
	
	 
				if (option.eq.4) then ! vecteur instantane de rotation dans OXYZ
		
! ############################################################
! Quaternion de rotation instantannée dans le repere celeste 
! q_{\Omega} = 2 \dot{q} \overline{q} 
! ############################################################
				
					call qua_mul(DQUA, QUAC, QUAW)
				
					call qua_mul_sca(QUAW,2.d0)
					
					wx = QUAW(2) * 1d9
					wy = QUAW(3) * 1d9
					wz = QUAW(4) * 1d9	
	
c composantes du vecteur instantane de rotation dans le repere terrestre 
					write(6,1030)  mjd_tai,an, mois, jour,IHMSF(1),IHMSF(2),IHMSF(3), IHMSF(4), wx, wy, wz
1030				        format(f22.13,i6,2(i3),3(i3),'.',I9.9,3f18.7)


				endif
	
			endif  ! fin de la condition : option 3 ou 4

		endif  ! fin de la condition : option <=4 

! ------------------------------------------------------------------------------------------------------------------------	
! option= 5: precession-nutation de l'axe de figure d'une Terre réelle et d'une Terre rigide (modele de Kinoshita Souchay)
! ------------------------------------------------------------------------------------------------------------------------
		if (option.eq.5) then	
			
			
			! matrice de precession-nutation KS (Terre rigide) du CRF vers le TRF
			
			! IAU 1976 precession matrix.
       			call iau_pmat06 ( jd_tt_0h, d_jd_tt, RP )

			! Kinoshita Souchay Earth Rigid Model of nutation components.
       			call iau_nut_ks ( jd_tt_0h, d_jd_tt, psi, eps )
       
			! Obliquity (IAU 1980).
       			epsa = iau_obl06 ( jd_tt_0h, d_jd_tt )

			! Nutation matrix
      		 	CALL iau_numat ( epsa, psi, eps, RN )

			! Form observed precession-nutation matrix.
       			CALL iau_rxr ( RN, RP, RPN )
						

			! les coordonnes dans le CRF de l'axe de figure sont [X,Y,Z] = transposee(PN) * [0,0,1] 
			X1 = RPN(3,1) / secrad ! conversion en seconde
			Y1 = RPN(3,2) / secrad ! conversion en seconde
			write(6,1040) 	mjd_tai,an, mois, jour, IHMSF(1),IHMSF(2),IHMSF(3),IHMSF(4), X/secrad, Y/secrad, X1,Y1 
1040			format(f22.13,i6,2(i3),3(i3),'.',I9.9, 4f20.6)
				 
		endif
	         
	enddo 
	

	call cpu_time(finish)
        !print '("quaternion CPU Time = ",f12.6," seconds.")',finish-start
	
	end 

!#########################################################################
	subroutine quaternion(X,Y,Z,sCIO,theta,sTIO,xpi,ypi,QUA,MQ,optq,optp)
!
!	input   PRT X,Y,Z, sCIO,theta,sTIO,xpi,ypi
!       	optq = 0 : quaternion total                 : ITRF --> ICRF
!       	optq = 1 : quaternion PN                    : CIF  --> ICRF
!       	optq = 2 : quaternion diurnal rotation      : TIF  --> CIF
!       	optq = 3 : quaternion PM                    : ITRF --> TIF
!       	optq = 4 : quaternion PM + diurnal rotation : ITRF --> CIF
!               optp     :  1 expression approchee / 2 expression exacte 
!        output QUA      : quaternion de rotation 
!               MQ       : matrice de rotation 
	
	implicit none
	double precision X, Y, Z, sCIO, theta, theta1, sTIO, xpi, ypi, A, B, C, D, fac, QUA(4), MQ(3,3)
	integer optq,optp
		
	
	if (optq.eq.0) theta1 = (theta + sTIO - sCIO)/2.d0
	
	if (optq.eq.1) then  ! quaternion de precession-nutation
		theta1 = -sCIO/2.d0
		xpi   = 0.d0
		ypi   = 0.d0		
	endif
	
	if (optq.eq.2) then  ! quaternion de rotation 
		X     = 0.d0
		Y     = 0.d0 
		Z     = 1.d0
		theta1 = theta/2.d0
		xpi   = 0.d0
		ypi   = 0.d0		
	endif
	
	if (optq.eq.3) then  ! quaternion mouvement du pôle du repere terrestre intermediaire à l'ITRF
		theta1 = sTIO/2.d0
		X  = 0.d0
		Y  = 0.d0
		Z  = 1.d0		
	endif

	if (optq.eq.4) then  ! quaternion ITRF --> CIF
	        theta1 = (theta + sTIO)/2.d0
		X  = 0.d0
		Y  = 0.d0
		Z  = 1.d0		
	endif
		
	fac   = dsqrt(1.d0/2.d0/(1+Z))
	
	if (optp.eq.1) then	! expression approchée
		QUA(1) = fac * ( dcos(theta1) * (  1.d0 + Z + X * xpi/2.d0 - Y * ypi/2.d0  ) + dsin(theta1) * ( X * ypi/2.d0  + Y * xpi/2.d0  ) )
		QUA(2) = fac * ( dcos(theta1) * (-(1.d0 + Z) * ypi/2.d0       - Y )          + dsin(theta1) * ( (1.d0 + Z) * xpi/2.d0 + X ) )
		QUA(3) = fac * ( dcos(theta1) * (-(1.d0 + Z) * xpi/2.d0 + X )                + dsin(theta1) * (-(1.d0 + Z) * ypi/2.d0 + Y) )
		QUA(4) = fac * ( dcos(theta1) * ( X * ypi/2.d0 + Y * xpi/2.d0 )    + dsin(theta1) * ( 1.d0 + Z - X * xpi/2.d0   + Y * ypi/2.d0) )
	endif
				
	if (optp.eq.2) then ! expression exacte		
		A =  dcos(xpi/2d0) * dcos(ypi/2d0)
		B =  dcos(xpi/2d0) * dsin(ypi/2d0)
		C =  dsin(xpi/2d0) * dcos(ypi/2d0)
		D = -dsin(xpi/2d0) * dsin(ypi/2d0)
		
		QUA(1) = fac * ( dcos(theta1) * (  (1.d0 + Z) * A + X * C - Y * B  ) + dsin(theta1) * ( -(1.d0 + Z) * D + X * B + Y * C ) )  
		QUA(2) = fac * ( dcos(theta1) * ( -(1.d0 + Z) * B + X * D - Y * A  ) + dsin(theta1) * (  (1.d0 + Z) * C + X * A + Y * D ) )
		QUA(3) = fac * ( dcos(theta1) * ( -(1.d0 + Z) * C + X * A + Y * D  ) + dsin(theta1) * ( -(1.d0 + Z) * B - X * D + Y * A ) )
		QUA(4) = fac * ( dcos(theta1) * (  (1.d0 + Z) * D + X * B + Y * C  ) + dsin(theta1) * (  (1.d0 + Z) * A - X * C + Y * B ) )  
	endif
	
	! ###################
	! matrice equivalente
	! ###################
	MQ(1,1) = QUA(1)**2 + QUA(2)**2 - QUA(3)**2 - QUA(4)**2 
	MQ(1,2) = 2.D0 * ( QUA(2) * QUA(3) - QUA(4) * QUA(1) )	
	MQ(1,3) = 2.D0 * ( QUA(2) * QUA(4) + QUA(3) * QUA(1) )	
	MQ(2,1) = 2.D0 * ( QUA(2) * QUA(3) + QUA(4) * QUA(1) )	
	MQ(2,2) = QUA(1)**2 - QUA(2)**2 + QUA(3)**2 - QUA(4)**2   
	MQ(2,3) = 2.D0 * ( QUA(3) * QUA(4) - QUA(2) * QUA(1) )				   
	MQ(3,1) = 2.D0 * ( QUA(2) * QUA(4) - QUA(3) * QUA(1) )	
	MQ(3,2) = 2.D0 * ( QUA(3) * QUA(4) + QUA(2) * QUA(1) )	
	MQ(3,3) = QUA(1)**2 - QUA(2)**2 - QUA(3)**2  + QUA(4)**2   
	
	
	end

! ########################################################################################################
! Function qui calcule TAI - UTC en fonction du mjd
! Entrée : mjd (double)
! Sortie : TAI_moins_UTC (double)
! Ce programme remplace l'ancien programme 
! Pour les sauts de seconde, il le fichier /obs/ftp/iers/bul/bulc/Leap_Second.dat
 
	function mjd_to_TAIUTC(mjd)
			
	implicit none
	double precision  mjd, mjd_to_TAIUTC
	character*256 ligne
	character*1 caractere_debut_ligne
	character*6 MJD_FICHIER
			
	double precision TABLEAU_SAUTS_SECONDE(1000)  ! on prévoit 1000 saut de seconde dans le fichier
	integer NB_SAUTS_SECONDE
			
	integer i
			
	double precision TAI_moins_UTC
			
	NB_SAUTS_SECONDE=0
			
	open(4, status='old', file='/obs/ftp/iers/bul/bulc/Leap_Second.dat',err=500)
			
20	read (4,*, end=30) MJD_FICHIER
				
	caractere_debut_ligne=MJD_FICHIER(1:1)
			
	if((caractere_debut_ligne.NE."#").AND.(MJD_FICHIER.NE."")) then
		NB_SAUTS_SECONDE = NB_SAUTS_SECONDE + 1
		read(MJD_FICHIER,*) TABLEAU_SAUTS_SECONDE(NB_SAUTS_SECONDE)	
	endif
				
	goto 20
			

30	continue

	close(4)
		
	if(mjd >= 41317.0)  then
 												
		 TAI_moins_UTC=9.0
				
		do 40 i=1,NB_SAUTS_SECONDE	
			if(mjd>= TABLEAU_SAUTS_SECONDE(i)) TAI_moins_UTC = TAI_moins_UTC + 1.d0						
40		continue
			
	else
							
		if(mjd < 41317)  TAI_moins_UTC = 4.213170d0 + (mjd - 39126.d0)*0.002592d0  ! jusqu'au 1 jan  1972
	   	if(mjd < 39887)  TAI_moins_UTC = 4.313170d0 + (mjd - 39126.d0)*0.002592d0  ! jusqu'au 1 fev  1968
		if(mjd < 39126)  TAI_moins_UTC = 3.840130d0 + (mjd - 38761.d0)*0.001296d0  ! jusqu'au 1 jan  1966
		if(mjd < 39004)  TAI_moins_UTC = 3.740130d0 + (mjd - 38761.d0)*0.001296d0  ! jusqu'au 1 sept 1965
	   	if(mjd < 38942)  TAI_moins_UTC = 3.640130d0 + (mjd - 38761.d0)*0.001296d0  ! jusqu'au 1 juil 1965
		if(mjd < 38820)  TAI_moins_UTC = 3.540130d0 + (mjd - 38761.d0)*0.001296d0  ! jusqu'au 1 mars 1965
	   	if(mjd < 38761)  TAI_moins_UTC = 3.440130d0 + (mjd - 38761.d0)*0.001296d0  ! jusqu'au 1 jan  1965
	   	if(mjd < 38639)  TAI_moins_UTC = 3.340130d0 + (mjd - 38761.d0)*0.001296d0  ! jusqu'au 1 sept 1964
	   	if(mjd < 38486)  TAI_moins_UTC = 3.240130d0 + (mjd - 38761.d0)*0.001296d0  ! jusqu'au 1 avr  1964
	   	if(mjd < 38395)  TAI_moins_UTC = 1.945858d0 + (mjd - 37665.d0)*0.0011232d0 ! jusqu'au 1 jan  1964
	   	if(mjd < 38334)  TAI_moins_UTC = 1.845858d0 + (mjd - 37665.d0)*0.0011232d0 ! jusqu'au 1 nov  1963
	   	if(mjd < 37665)  TAI_moins_UTC = 1.372818d0 + (mjd - 37300.d0)*0.001296d0  ! jusqu'au 1 jan  1962
	   	if(mjd < 37512)  TAI_moins_UTC = 1.422818d0 + (mjd - 37300.d0)*0.001296d0  ! jusqu'au 1 aout 1961
	   	if(mjd < 37300)  TAI_moins_UTC = 0.d0
				
						
	endif
	   								
	mjd_to_TAIUTC = TAI_moins_UTC
			
400	continue
							
	RETURN
			
500  	write(*,*) "PROBLEME OUVERTURE FICHIER SAUT DE SECONDE !??!"
	goto 400	
	
	return 			
	END 
	
       
! #######################################################################
! Programme spline:interpolation par une fonction spline cubique complete

	function SPLINE(X,Y,n,xs)

! entree :   X     : tableau des dates
!            xs    : date a interpolee
!            Y     : tableau des valeurs a interpoler
!            n     : dimensions des tableaux x et y
!            ns    : dimensions des tableaux xs et ys
!	     alpha : valeur de la derivee du signal en t0
!            beta  : valeur de la derivee du signal en tn 
! sortie     spline: valeurs interpolees a la date dates xs
! #######################################################################

        implicit none	
	integer i,n,n1,j
	parameter (n1=10000)	
	double precision X(n),Y(n)
	double precision S(n1),T(n1),d(n1),a(n1),b(n1),c(n1)
	double precision M(n1),h(n1)
	double precision alpha,beta,xs,spline
	

! La determination des coefficients des polynomes qui constituent la fonction
! spline consiste a resoudre un systeme tridiagonal par la methode du double
! balayage de Cholesky


	do i=1,n-1
	h(i) = x(i+1)-x(i)
	enddo
	
	alpha =(y(2)-y(1))/h(1)
	beta  =(y(n)-y(n-1))/h(n-1)
	
	b(1) = h(1)/3.d0
	c(1) = h(1)/6.d0
	d(1)=(y(2)-y(1))/h(1)-alpha
	
	S(1) =  d(1)/b(1)
	T(1) = -c(1)/b(1)

	do i=2,n-1	
	
	a(i) = h(i-1)/6.d0
	b(i) = (h(i-1)+h(i))/3.d0
	c(i) = h(i)/6.d0
	d(i) = (y(i+1)-y(i))/h(i) - (y(i)-y(i-1))/h(i-1)
	
	S(i) = (d(i)-a(i)*S(i-1)) / (b(i)+a(i)*T(i-1))
	T(i) =-c(i)/(b(i)+a(i)*T(i-1))	
	enddo
	

	a(n) = h(n-1)/6.d0
	b(n) = h(n-1)/3.d0
	
	d(n) = beta - (y(n)-y(n-1))/h(n-1)
	
! calcul des coeeficients M(i)
		
	M(n)=(d(n)-a(n)*S(n-1))/(a(n)*T(n-1)+b(n))

	do i=n-1,1,-1
	M(i)=S(i)+T(i)*M(i+1)	 
	enddo


	if  (xs.lt.x(1)) then
	j = 1
	     spline = 1.d0 / h(j)
     & *(  M(j)  /6.d0 * (x(j+1)-xs)**3  
     &    -M(j+1)/6.d0 * (x(j  )-xs)**3 
     &    +( y(j  )-M(j)   *h(j)**2 /6.d0 ) * (x(j+1)-xs) 
     &    -( y(j+1)-M(j+1) *h(j)**2 /6.d0 ) * (x(j  )-xs)  )
 
	endif

	if  ( (xs.ge.x(1)).and.(xs.lt.x(n-1)) ) then
	
	do j=1,n-2
	if  ( (xs.ge.x(j)).and.(xs.lt.x(j+1)) ) then

	     spline = 1.d0 / h(j)
     &  *(  M(j)  /6.d0 * (x(j+1)-xs)**3  
     &    -M(j+1)/6.d0 * (x(j  )-xs)**3 
     &    +( y(j  )-M(j)   *h(j)**2 /6.d0 ) * (x(j+1)-xs) 
     &    -( y(j+1)-M(j+1) *h(j)**2 /6.d0 ) * (x(j  )-xs)  )
          
	endif
	enddo
	endif
	
	if (xs.ge.x(n-1)) then
        j = n-1 
	     spline = 1.d0 / h(j)
     &  * (  M(j)  /6.d0 * (x(j+1)-xs)**3  
     &    -M(j+1)/6.d0 * (x(j  )-xs)**3 
     &    +( y(j  )-M(j)   *h(j)**2 /6.d0 ) * (x(j+1)-xs) 
     &    -( y(j+1)-M(j+1) *h(j)**2 /6.d0 ) * (x(j  )-xs)  )

	endif
	
	return

	end
	
!##################################
! produit de deux quaternions
	subroutine qua_mul(Q1,Q2,Q)
!##################################
	
	implicit none
	double precision Q1(4), Q2(4), Q(4)
	
	Q(1) = Q1(1) * Q2(1) - Q1(2) * Q2(2) - Q1(3) * Q2(3) - Q1(4) * Q2(4) 
	Q(2) = Q1(1) * Q2(2) + Q1(2) * Q2(1) + Q1(3) * Q2(4) - Q1(4) * Q2(3) 
	Q(3) = Q1(1) * Q2(3) - Q1(2) * Q2(4) + Q1(3) * Q2(1) + Q1(4) * Q2(2) 
	Q(4) = Q1(1) * Q2(4) + Q1(2) * Q2(3) - Q1(3) * Q2(2) + Q1(4) * Q2(1) 
	
	return 
	end
	
!##################################
!  quaternion conjuge
	subroutine qua_con(Q,QC)
!##################################
	
	implicit none
	double precision Q(4), QC(4)
	
	QC(1) =  Q(1) 
	QC(2) = -Q(2)
	QC(3) = -Q(3) 
	QC(4) = -Q(4)  
	
	return 
	end	

!##################################
!  quaternion cmutliplication par un scalaire
	subroutine qua_mul_sca(Q,sca)
!##################################
	
	implicit none
	double precision Q(4),sca
	
	Q(1) = Q(1) * sca 
	Q(2) = Q(2) * sca
	Q(3) = Q(3) * sca 
	Q(4) = Q(4) * sca  
	
	return 
	end


c##################################
c  quaternion cmutliplication par un scalaire
	subroutine qua_matequiv(Q,MQ)
c
c entree : QUA(4) quaternion de rotation 
c sortie : Q(3,3) matrice de transformation de coordonnees associée
c##################################
	
	implicit none
	double precision Q(4), MQ(3,3)

	MQ(1,1) = Q(1)**2 + Q(2)**2 - Q(3)**2 - Q(4)**2 
	MQ(1,2) = 2.D0 * ( Q(2) * Q(3) - Q(4) * Q(1) )	
	MQ(1,3) = 2.D0 * ( Q(2) * Q(4) + Q(3) * Q(1) )	
	MQ(2,1) = 2.D0 * ( Q(2) * Q(3) + Q(4) * Q(1) )	
	MQ(2,2) = Q(1)**2 - Q(2)**2 + Q(3)**2 - Q(4)**2   
	MQ(2,3) = 2.D0 * ( Q(3) * Q(4) - Q(2) * Q(1) )				   
	MQ(3,1) = 2.D0 * ( Q(2) * Q(4) - Q(3) * Q(1) )	
	MQ(3,2) = 2.D0 * ( Q(3) * Q(4) + Q(2) * Q(1) )	
	MQ(3,3) = Q(1)**2 - Q(2)**2 - Q(3)**2  + Q(4)**2    

	end

	
C----------------------------------------------------------------
      SUBROUTINE PMUT1_OCEANS (rjd,cor_x,cor_y,cor_ut1,cor_lod)
C
C    This subroutine provides, in time domain, the diurnal/subdiurnal
C    tidal effets on polar motion ("), UT1 (s) and LOD (s). The tidal terms,
C    listed in the program above, have been extracted from the procedure   
C    ortho_eop.f coed by Eanes in 1997.
C    
C    N.B.:  The fundamental lunisolar arguments are those of Simon et al.  
C
C    These corrections should be added to "average"
C    EOP values to get estimates of the instantaneous values.
C
C     PARAMETERS ARE :
C     rjd      - epoch of interest given in mjd
C     cor_x    - tidal correction in x (sec. of arc)
C     cor_y    - tidal correction in y (sec. of arc)
C     cor_ut1  - tidal correction in UT1-UTC (sec. of time)
C     cor_lod  - tidal correction in length of day (sec. of time)
C
C     coded by Ch. Bizouard (2002), initially coded by McCarthy and 
C     D.Gambis(1997) for the 8 prominent tidal waves.  
      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      
      INTEGER nlines
      PARAMETER(nlines=71)
      DOUBLE PRECISION ARG(6),    ! Array of the tidal arguments   
     .                 DARG(6)    ! Array of their time derivative 
      REAL*4 XCOS(nlines),XSIN(nlines),
     .YCOS(nlines),YSIN(nlines),UTCOS(nlines),UTSIN(nlines)
      INTEGER NARG(nlines,6)
      
      halfpi = 1.5707963267948966d0
      secrad=2.d0*halfpi/(180.d0*3600.d0)	

c  Oceanic tidal terms present in x (microas),y(microas),ut1(microseconds)       
c  NARG(j,6) : Multipliers of GMST+pi and Delaunay arguments. 
	
       data( 
     & NARG(j,1),NARG(j,2),NARG(j,3),NARG(j,4),NARG(j,5),NARG(j,6),
     & XSIN(j),XCOS(j),YSIN(j),YCOS(j),UTSIN(j),UTCOS(j),j=1,nlines)/
     &1,-1, 0,-2,-2,-2,  -0.05,   0.94,  -0.94,  -0.05,  0.396, -0.078,
     &1,-2, 0,-2, 0,-1,   0.06,   0.64,  -0.64,   0.06,  0.195, -0.059,
     &1,-2, 0,-2, 0,-2,   0.30,   3.42,  -3.42,   0.30,  1.034, -0.314,
     &1, 0, 0,-2,-2,-1,   0.08,   0.78,  -0.78,   0.08,  0.224, -0.073,
     &1, 0, 0,-2,-2,-2,   0.46,   4.15,  -4.15,   0.45,  1.187, -0.387,
     &1,-1, 0,-2, 0,-1,   1.19,   4.96,  -4.96,   1.19,  0.966, -0.474,
     &1,-1, 0,-2, 0,-2,   6.24,  26.31, -26.31,   6.23,  5.118, -2.499,
     &1, 1, 0,-2,-2,-1,   0.24,   0.94,  -0.94,   0.24,  0.172, -0.090,
     &1, 1, 0,-2,-2,-2,   1.28,   4.99,  -4.99,   1.28,  0.911, -0.475,
     &1, 0, 0,-2, 0, 0,  -0.28,  -0.77,   0.77,  -0.28, -0.093,  0.070,
     &1, 0, 0,-2, 0,-1,   9.22,  25.06, -25.06,   9.22,  3.025, -2.280,
     &1, 0, 0,-2, 0,-2,  48.82, 132.91,-132.90,  48.82, 16.020,-12.069,
     &1,-2, 0, 0, 0, 0,  -0.32,  -0.86,   0.86,  -0.32, -0.103,  0.078,
     &1, 0, 0, 0,-2, 0,  -0.66,  -1.72,   1.72,  -0.66, -0.194,  0.154,
     &1,-1, 0,-2, 2,-2,  -0.42,  -0.92,   0.92,  -0.42, -0.083,  0.074,
     &1, 1, 0,-2, 0,-1,  -0.30,  -0.64,   0.64,  -0.30, -0.057,  0.050,
     &1, 1, 0,-2, 0,-2,  -1.61,  -3.46,   3.46,  -1.61, -0.308,  0.271,
     &1,-1, 0, 0, 0, 0,  -4.48,  -9.61,   9.61,  -4.48, -0.856,  0.751,
     &1,-1, 0, 0, 0,-1,  -0.90,  -1.93,   1.93,  -0.90, -0.172,  0.151,
     &1, 1, 0, 0,-2, 0,  -0.86,  -1.81,   1.81,  -0.86, -0.161,  0.137,
     &1, 0,-1,-2, 2,-2,   1.54,   3.03,  -3.03,   1.54,  0.315, -0.189,
     &1, 0, 0,-2, 2,-1,  -0.29,  -0.58,   0.58,  -0.29, -0.062,  0.035,
     &1, 0, 0,-2, 2,-2,  26.13,  51.25, -51.25,  26.13,  5.512, -3.095,
     &1, 0, 1,-2, 2,-2,  -0.22,  -0.42,   0.42,  -0.22, -0.047,  0.025,
     &1, 0,-1, 0, 0, 0,  -0.61,  -1.20,   1.20,  -0.61, -0.134,  0.070,
     &1, 0, 0, 0, 0, 1,   1.54,   3.00,  -3.00,   1.54,  0.348, -0.171,
     &1, 0, 0, 0, 0, 0, -77.48,-151.74, 151.74, -77.48,-17.620,  8.548,
     &1, 0, 0, 0, 0,-1, -10.52, -20.56,  20.56, -10.52, -2.392,  1.159,
     &1, 0, 0, 0, 0,-2,   0.23,   0.44,  -0.44,   0.23,  0.052, -0.025,
     &1, 0, 1, 0, 0, 0,  -0.61,  -1.19,   1.19,  -0.61, -0.144,  0.065,
     &1, 0, 0, 2,-2, 2,  -1.09,  -2.11,   2.11,  -1.09, -0.267,  0.111,
     &1,-1, 0, 0, 2, 0,  -0.69,  -1.43,   1.43,  -0.69, -0.288,  0.043,
     &1, 1, 0, 0, 0, 0,  -3.46,  -7.28,   7.28,  -3.46, -1.610,  0.187,
     &1, 1, 0, 0, 0,-1,  -0.69,  -1.44,   1.44,  -0.69, -0.320,  0.037,
     &1, 0, 0, 0, 2, 0,  -0.37,  -1.06,   1.06,  -0.37, -0.407, -0.005,
     &1, 2, 0, 0, 0, 0,  -0.17,  -0.51,   0.51,  -0.17, -0.213, -0.005,
     &1, 0, 0, 2, 0, 2,  -1.10,  -3.42,   3.42,  -1.09, -1.436, -0.037,
     &1, 0, 0, 2, 0, 1,  -0.70,  -2.19,   2.19,  -0.70, -0.921, -0.023,
     &1, 0, 0, 2, 0, 0,  -0.15,  -0.46,   0.46,  -0.15, -0.193, -0.005,
     &1, 1, 0, 2, 0, 2,  -0.03,  -0.59,   0.59,  -0.03, -0.396, -0.024,
     &1, 1, 0, 2, 0, 1,  -0.02,  -0.38,   0.38,  -0.02, -0.253, -0.015,
     &2,-3, 0,-2, 0,-2,  -0.49,  -0.04,   0.63,   0.24, -0.089, -0.011,
     &2,-1, 0,-2,-2,-2,  -1.33,  -0.17,   1.53,   0.68, -0.224, -0.032,
     &2,-2, 0,-2, 0,-2,  -6.08,  -1.61,   3.13,   3.35, -0.637, -0.177,
     &2, 0, 0,-2,-2,-2,  -7.59,  -2.05,   3.44,   4.23, -0.745, -0.222,
     &2, 0, 1,-2,-2,-2,  -0.52,  -0.14,   0.22,   0.29, -0.049, -0.015,
     &2,-1,-1,-2, 0,-2,   0.47,   0.11,  -0.10,  -0.27,  0.033,  0.013,
     &2,-1, 0,-2, 0,-1,   2.12,   0.49,  -0.41,  -1.23,  0.141,  0.058,
     &2,-1, 0,-2, 0,-2, -56.87, -12.93,  11.15,  32.88, -3.795, -1.556,
     &2,-1, 1,-2, 0,-2,  -0.54,  -0.12,   0.10,   0.31, -0.035, -0.015,
     &2, 1, 0,-2,-2,-2, -11.01,  -2.40,   1.89,   6.41, -0.698, -0.298,
     &2, 1, 1,-2,-2,-2,  -0.51,  -0.11,   0.08,   0.30, -0.032, -0.014,
     &2,-2, 0,-2, 2,-2,   0.98,   0.11,  -0.11,  -0.58,  0.050,  0.022,
     &2, 0,-1,-2, 0,-2,   1.13,   0.11,  -0.13,  -0.67,  0.056,  0.025,
     &2, 0, 0,-2, 0,-1,  12.32,   1.00,  -1.41,  -7.31,  0.605,  0.266,
     &2, 0, 0,-2, 0,-2,-330.15, -26.96,  37.58, 195.92,-16.195, -7.140,
     &2, 0, 1,-2, 0,-2,  -1.01,  -0.07,   0.11,   0.60, -0.049, -0.021,
     &2,-1, 0,-2, 2,-2,   2.47,  -0.28,  -0.44,  -1.48,  0.111,  0.034,
     &2, 1, 0,-2, 0,-2,   9.40,  -1.44,  -1.88,  -5.65,  0.425,  0.117,
     &2,-1, 0, 0, 0, 0,  -2.35,   0.37,   0.47,   1.41, -0.106, -0.029,
     &2,-1, 0, 0, 0,-1,  -1.04,   0.17,   0.21,   0.62, -0.047, -0.013,
     &2, 0,-1,-2, 2,-2,  -8.51,   3.50,   3.29,   5.11, -0.437, -0.019,
     &2, 0, 0,-2, 2,-2,-144.13,  63.56,  59.23,  86.56, -7.547, -0.159,
     &2, 0, 1,-2, 2,-2,   1.19,  -0.56,  -0.52,  -0.72,  0.064,  0.000,
     &2, 0, 0, 0, 0, 1,   0.49,  -0.25,  -0.23,  -0.29,  0.027, -0.001,
     &2, 0, 0, 0, 0, 0, -38.48,  19.14,  17.72,  23.11, -2.104,  0.041,
     &2, 0, 0, 0, 0,-1, -11.44,   5.75,   5.32,   6.87, -0.627,  0.015,
     &2, 0, 0, 0, 0,-2,  -1.24,   0.63,   0.58,   0.75, -0.068,  0.002,
     &2, 1, 0, 0, 0, 0,  -1.77,   1.79,   1.71,   1.04, -0.146,  0.037,
     &2, 1, 0, 0, 0,-1,  -0.77,   0.78,   0.75,   0.45, -0.064,  0.017,
     &2, 0, 0, 2, 0, 2,  -0.33,   0.62,   0.65,   0.19, -0.049,  0.018/

      T = (rjd - 51544.5D0)/36525.0D0  ! julian century

C Arguments in the following order : chi=GMST+pi,l,lp,F,D,Omega
C et leur derivee temporelle 

      ARG(1) = (67310.54841d0 +
     .        (876600d0*3600d0 + 8640184.812866d0)*T +
     .         0.093104d0*T**2 -
     .         6.2d-6*T**3)*15.0d0 + 648000.0d0
      ARG(1)=dmod(ARG(1),1296000d0)*secrad 
   
      DARG(1) = (876600d0*3600d0 + 8640184.812866d0 
     .         + 2.d0 * 0.093104d0 * T - 3.d0 * 6.2d-6*T**2)*15.d0
      DARG(1) = DARG(1)* secrad / 36525.0D0   ! rad/day


      ARG(2) = -0.00024470d0*T**4 + 0.051635d0*T**3 + 31.8792d0*T**2
     .  + 1717915923.2178d0*T + 485868.249036d0
      ARG(2) = DMOD(ARG(2),1296000d0)*secrad
      
      DARG(2) = -4.d0*0.00024470d0*T**3 + 3.d0*0.051635d0*T**2 
     .  + 2.d0*31.8792d0*T + 1717915923.2178d0 
      DARG(2) = DARG(2)* secrad / 36525.0D0   ! rad/day

      ARG(3) = -0.00001149d0*T**4 - 0.000136d0*T**3
     .  -  0.5532d0*T**2
     .  + 129596581.0481d0*T + 1287104.79305d0
      ARG(3) = DMOD(ARG(3),1296000d0)*secrad

      DARG(3) = -4.D0*0.00001149d0*T**3 - 3.d0*0.000136d0*T**2
     .  -  2.D0*0.5532d0*T + 129596581.0481d0
      DARG(3) = DARG(3)* secrad / 36525.0D0   ! rad/day
          
      ARG(4) = 0.00000417d0*T**4 - 0.001037d0*T**3 - 12.7512d0*T**2
     .  + 1739527262.8478d0*T + 335779.526232d0
      ARG(4) = DMOD(ARG(4),1296000d0)*secrad

      DARG(4) = 4.d0*0.00000417d0*T**3 - 3.d0*0.001037d0*T**2 
     .- 2.d0 * 12.7512d0*T + 1739527262.8478d0 
      DARG(4) = DARG(4)* secrad / 36525.0D0   ! rad/day
    
      ARG(5) = -0.00003169d0*T**4 + 0.006593d0*T**3 - 6.3706d0*T**2
     .  + 1602961601.2090d0*T + 1072260.70369d0
      ARG(5) = DMOD(ARG(5),1296000d0)*secrad

      DARG(5) = -4.d0*0.00003169d0*T**3 + 3.d0*0.006593d0*T**2
     . - 2.d0 * 6.3706d0*T + 1602961601.2090d0
      DARG(5) = DARG(5)* secrad / 36525.0D0   ! rad/day

      ARG(6) = -0.00005939d0*T**4 + 0.007702d0*T**3
     .  + 7.4722d0*T**2
     .  - 6962890.2665d0*T + 450160.398036d0
      ARG(6) = DMOD(ARG(6),1296000d0)*secrad

      DARG(6) = -4.d0*0.00005939d0*T**3 + 3.d0 * 0.007702d0*T**2
     .  + 2.d0 * 7.4722d0*T - 6962890.2665d0
      DARG(6) = DARG(6)* secrad / 36525.0D0   ! rad/day

C CORRECTIONS

	cor_x  = 0.d0
	cor_y  = 0.d0
	cor_ut1= 0.d0
	cor_lod= 0.d0

 	do j=1,nlines
 	
	ag  = 0.d0
 	dag = 0.d0
		do i=1,6
 		ag  = ag  + dble(narg(j,i))*ARG(i)
 		dag = dag + dble(narg(j,i))*DARG(i)
		enddo
	ag=dmod(ag,4.d0*halfpi)

        cor_x  =cor_x   +dble(XCOS(j))*dcos(ag) +dble(XSIN(j))*dsin(ag)
        cor_y  =cor_y   +dble(YCOS(j))*dcos(ag) +dble(YSIN(j))*dsin(ag)
        cor_ut1=cor_ut1 +dble(UTCOS(j))*dcos(ag)+dble(UTSIN(j))*dsin(ag)
        cor_lod=cor_lod -(-dble(UTCOS(j))*dsin(ag) 
     &                    +dble(UTSIN(j))*dcos(ag) ) * dag   	 

        enddo
  
      cor_x = cor_x * 1.0d-6   ! arcsecond (")
      cor_y = cor_y * 1.0d-6   ! arcsecond (")
      cor_ut1 = cor_ut1 * 1.0d-6   ! second (s)
      cor_lod = cor_lod * 1.0d-6   ! second (s)
 
      RETURN
      END
      	
C----------------------------------------------------------------
      SUBROUTINE PM_GRAVI (rjd,cor_x,cor_y)
C
C    This subroutine provides, in time domain, the diurnal
C    lunisolar effet on polar motion (")
C    
C    N.B.:  The fundamental lunisolar arguments are those of Simon et al.  
C
C    These corrections should be added to "average"
C    EOP values to get estimates of the instantaneous values.
C
C     PARAMETERS ARE :
C     rjd      - epoch of interest given in mjd
C     cor_x    - tidal correction in x (sec. of arc)
C     cor_y    - tidal correction in y (sec. of arc)
C
C     coded by Ch. Bizouard (2002)
      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      
      INTEGER nlines
      PARAMETER(nlines=10)
      DOUBLE PRECISION ARG(6)    ! Array of the tidal arguments   
      REAL*4 XCOS(nlines),XSIN(nlines),YCOS(nlines),YSIN(nlines)
      INTEGER NARG(nlines,6)
      
      halfpi = 1.5707963267948966d0
      secrad=2.d0*halfpi/(180.d0*3600.d0)	

c  Diurnal lunisolar tidal terms present in x (microas),y(microas)      
c  NARG(j,6) : Multipliers of GMST+pi and Delaunay arguments. 
	
       data( 
     & NARG(j,1),NARG(j,2),NARG(j,3),NARG(j,4),NARG(j,5),NARG(j,6),
     & XSIN(j),XCOS(j),YSIN(j),YCOS(j),j=1,nlines)/    
     & 1,-1, 0,-2, 0,-1,    -.44,   .25,   -.25,  -.44,
     & 1,-1, 0,-2, 0,-2,   -2.31,  1.32,  -1.32, -2.31,
     & 1, 1, 0,-2,-2,-2,    -.44,   .25,   -.25,  -.44,
     & 1, 0, 0,-2, 0,-1,   -2.14,  1.23,  -1.23, -2.14,
     & 1, 0, 0,-2, 0,-2,  -11.36,  6.52,  -6.52,-11.36,
     & 1,-1, 0, 0, 0, 0,     .84,  -.48,    .48,   .84,
     & 1, 0, 0,-2, 2,-2,   -4.76,  2.73,  -2.73, -4.76,
     & 1, 0, 0, 0, 0, 0,   14.27, -8.19,   8.19, 14.27,
     & 1, 0, 0, 0, 0,-1,    1.93, -1.11,   1.11,  1.93,
     & 1, 1, 0, 0, 0, 0,     .76,  -.43,    .43,   .76/
 
      T = (rjd - 51544.5D0)/36525.0D0  ! julian century

C Arguments in the following order : chi=GMST+pi,l,lp,F,D,Omega
C et leur derivee temporelle 

      ARG(1) = (67310.54841d0 +
     .        (876600d0*3600d0 + 8640184.812866d0)*T +
     .         0.093104d0*T**2 -
     .         6.2d-6*T**3)*15.0d0 + 648000.0d0
      ARG(1)=dmod(ARG(1),1296000d0)*secrad 
   

      ARG(2) = -0.00024470d0*T**4 + 0.051635d0*T**3 + 31.8792d0*T**2
     .  + 1717915923.2178d0*T + 485868.249036d0
      ARG(2) = DMOD(ARG(2),1296000d0)*secrad
      

      ARG(3) = -0.00001149d0*T**4 - 0.000136d0*T**3
     .  -  0.5532d0*T**2
     .  + 129596581.0481d0*T + 1287104.79305d0
      ARG(3) = DMOD(ARG(3),1296000d0)*secrad

          
      ARG(4) = 0.00000417d0*T**4 - 0.001037d0*T**3 - 12.7512d0*T**2
     .  + 1739527262.8478d0*T + 335779.526232d0
      ARG(4) = DMOD(ARG(4),1296000d0)*secrad

    
      ARG(5) = -0.00003169d0*T**4 + 0.006593d0*T**3 - 6.3706d0*T**2
     .  + 1602961601.2090d0*T + 1072260.70369d0
      ARG(5) = DMOD(ARG(5),1296000d0)*secrad

  
      ARG(6) = -0.00005939d0*T**4 + 0.007702d0*T**3
     .  + 7.4722d0*T**2
     .  - 6962890.2665d0*T + 450160.398036d0
      ARG(6) = DMOD(ARG(6),1296000d0)*secrad


C CORRECTIONS

	cor_x  = 0.d0
	cor_y  = 0.d0

 	do j=1,nlines
 	
	ag  = 0.d0
		do i=1,6
 		ag  = ag  + dble(narg(j,i))*ARG(i)
		enddo
	ag=dmod(ag,4.d0*halfpi)

        cor_x  =cor_x   +dble(XCOS(j))*dcos(ag) +dble(XSIN(j))*dsin(ag)
        cor_y  =cor_y   +dble(YCOS(j))*dcos(ag) +dble(YSIN(j))*dsin(ag) 

        enddo
  
      cor_x = cor_x * 1.0d-6   ! arcsecond (")
      cor_y = cor_y * 1.0d-6   ! arcsecond (")
 
      RETURN

      END


C=======================================================================C

      subroutine MJDDATE(mjd_tai,jour,mois,an,nh,nm,s) 
*----------------------------------------------------------------------- 
* 
*     Ref : GF9803 
* 
* --- Objet ------------------------------------------------------------- 
* 
*     Changement de calendrier : 
*     Date julienne modifiee TAI -> Date civile grégorienne UTC
*     prenant en compte les secondes intercalaires. 
* 
* --- Entree ------------------------------------------------------------ 
* 
*     mjd         Date julienne modifiee (reel dp) en TAI
*
* --- Sortie : Date civile gregorienne a la fraction de seconde pres ----
* 
*     jour        Quantieme du mois (entier). 
*     mois        Numero du mois (entier). 
*     an          Millesime astronomique de l'annee (entier). 
* 
*     nh          Nombre d'heures (entier). 
*     nm          Nombre de minutes (entier). 
*     s           Nombre de secondes (reel). 
* 
* --- Remarque --------------------------------------------------------- 
* 
*     Date du calendier (Gregorien/Julien) : 
*     Le calendier Gregorien commence le 15/10/1582 a 0h (DJ2299160.5). 
*     Avant cette date on utilise le calendrier Julien. 
* 
*     On utilise le formulaire de J. Meeus. 
*     La date julienne tdj doit etre superieure ou egale a 0. 
*     Si tdj<0, alors jour=mois=an=nh=nm=s=0. 
* 
*     L'annee kan est comptee en millesimes astronomiques. 
*     L'annee 1 avant J.C. des historiens correspond a l'an 0 des 
*     astronomes et ainsi de suite. 
*     Ex : l'an 80 avant J.C. correspond a l'an astronomique -79. 
* 
* --- Declarations ----------------------------------------------------- 
* 

      implicit none
      double precision tdj,s,a,b,c,d,e,z,t,f,alpha,mjd_to_TAIUTC,mjd_tai,mjd_utc,tdj_tai,tdj_utc
      double precision mjd_intercal_tai,mjd_intercal_utc,sintercal 
      integer jour,mois,an,nh,nm,nsint,i
* 
* --- Initialisation --------------------------------------------------- 
* 
       jour = 0 
       mois = 0 
       an   = 0 
       nh   = 0 
       nm   = 0       
       s    = 0.d0  
    
* --- Conversion ------------------------------------------------------- 
* 
       t  = mjd_tai + 0.5d0 + 2400000.5d0
       z  = dint(t)
    
       f  = mjd_tai - dint(mjd_tai) 
* 
       if (z.lt.2299161.d0) then 
       		a=z 
       else 
       		alpha = dint((z-1867216.25d0)/36524.25d0) 
         	a     = z + 1.d0 + alpha - dint(alpha/4.d0) 
       endif 
* 
       b = a + 1524.d0 
       c = dint((b-122.1d0)/365.25d0) 
       d = dint(365.25d0*c) 
       e = dint((b-d)/30.6001d0) 
* 
       jour = b - d - dint(30.6001d0*e) 
* 
       if (e.lt.13.5d0) then 
       		mois = e - 1.d0 
       else 
         	mois = e - 13.d0 
       endif 
* 
       if (mois.lt.2.5d0) then 
         	an = c - 4715.d0 
       else 
         	an = c - 4716.d0 
       endif 
* 
        f  = f * 24.d0 
        nh = dint(f) 
        f  = (f-dble(nh))*60.d0 
        nm = int(f) 
        f  = (f-dble(nm))*60.d0 
        s  = f 
 	
	
        nsint  = mjd_to_TAIUTC(mjd_tai) ! TAI - UTC
	
	do i = 1,10
		mjd_utc = mjd_tai -  dble(nsint) /86400.d0 
		nsint   = mjd_to_TAIUTC(mjd_utc)
	enddo
       
        ! calcul de la date civile (necessairement en UTC)         
	s = s - nsint
      
	if (s.le.-1d-4) then           
		if (mjd_to_TAIUTC(mjd_utc).eq.mjd_to_TAIUTC(mjd_utc-2.d0/86400.d0)) then  ! pas de saut de seconde autour de mjd_utc
			s  = s + 60.d0 ! cas normal 
		else ! pas de saut de seconde	
			s  = s + 61.d0 ! pour un saut de seconde
		endif
	
		nm = nm - 1 
		if (nm.eq.-1) then
			nm = 59
			nh = nh - 1
			if (nh.eq.-1) then   
                       		nh = 23 
                        	jour = jour - 1
				if (jour.eq.0) then
				
				mois = mois - 1			 
				if ((mois.eq.1).or.(mois.eq.3).or.(mois.eq.5).or.(mois.eq.7).or.(mois.eq.8).or.(mois.eq.10).or.(mois.eq.12)) then
					jour = 31				
				endif 

				if ((mois.eq.4).or.(mois.eq.6).or.(mois.eq.9).or.(mois.eq.11)) then
					jour = 30
				endif 

				if ((mois.eq.2).and.(mod(an,4).ne.0)) then ! annee non divisible par 4  - non bissextile
					jour = 28 
				endif  

				if ( (mois.eq.2).and.(mod(an,4).eq.0).and.(mod(an,100).eq.0).and.(mod(an,400).ne.0)) then ! non bissextile
					jour = 28
				endif  

				if ((mois.eq.2).and.(mod(an,4).eq.0).and.(mod(an,100).eq.0).and.(mod(an,400).eq.0)) then  ! bissextile
					jour = 29
				endif  

				if ((mois.eq.2).and.(mod(an,4).eq.0).and.(mod(an,100).ne.0).and.(mod(an,400).ne.0)) then  ! bissextile
					jour = 29
				endif  
				
				if (mois.eq.0) then
					mois = 12
					jour = 31
					an = an - 1
				endif
			        
				endif	
			endif	
	        endif
	   endif
	return 
	end

c #####################################################################
c This function performs the computation of the time shift TAI-UTC, 
c in seconds, to be removed from UT1-UTC = UT1 - TAI + TAI - UTC 
c in order to get UT1-TAI

c input  : D = Date in Modified Julian Day
c output : TARTUF = TAI-UTC 
c
c C'est la copie de TARTUF se trouvant dans la livrairie Lib90 de eoppc2
c ######################################################################
  	
       function TARTUF(D)
       implicit double precision (A-H,O-Z)
       double precision D,TARTUF
       double precision A(50),B(50),C(50),T(50)
   	
       A(1)=1.422818d0
       A(2)=1.372818d0
       A(3)=1.845858d0
       A(4)=1.945858d0
       A(5)=3.240130d0
      
       DO 10 I=6,11
  10   A(I)=A(I-1)+0.1d0
   
       A(12)=4.313170d0
       A(13)=4.213170d0
      
       DO 11 I=14,50
  11   A(I)=I-4
   
       B(1)=37300.d0
       B(2)=37300.d0
       B(3)=37665.d0
       B(4)=37665.d0
      
       DO 12 I=5,11
  12   B(I)=38761.d0
   
       B(12)=39126.d0
       B(13)=39126.d0
   
       DO 13 I=14,50
  13   B(I)=0.d0
   
       ! coefficients lineaires  
       DO 14 I=1,11
 14    C(I)=0.001296d0
   
       C(3)=0.0011232d0
       C(4)=0.0011232d0
       C(12)=0.002592d0
       C(13)=0.002592d0
 
       ! A partir du 1er juillet 1972 il n'y a plus de termes lineaires	
       DO 15 I=14,50
 15    C(I)=0.d0
             
       TA=37300.d0		  ! 010161
       T(1)=37512.d0		  ! 010861
       T(2)=37665.d0		  ! 010162
       T(3)=38334.d0		  ! 011163
       T(4)=38395.d0		  ! 010164
       T(5)=38486.d0		  ! 010464
       T(6)=38639.d0		  ! 010964
       T(7)=38761.d0		  ! 010165
       T(8)=38820.d0		  ! 010365
       T(9)=38942.d0		  ! 010765
       T(10)=39004.d0		  ! 010965
       T(11)=39126.d0		  ! 010166
       T(12)=39887.d0		  ! 010268
       T(13)=41317.d0		  ! 010172
       T(14)=41499.d0		  ! 010772
       T(15)=41683.d0		  ! 010173
C
C   LA BOUCLE 16 EST A JOUR POUR LE 1.7.82
C
       DO 16 i=16,22
       		T(I)=T(I-1)+365.d0
       		IF(MOD((I-15),4).EQ.0) T(I)=T(I)+1.d0  
       		! prise en compte des annees bissextiles
 16    CONTINUE
     
       T(23) = 44786.d0   ! 010781
       T(24) = 45151.d0   ! 010782
       T(25) = 45516.d0   ! 010783
       T(26) = 46247.d0	  ! 010785
       T(27) = 47161.d0	  ! 010188
       T(28) = 47892.d0   ! 010190
       T(29) = 48257.d0   ! 010191
       T(30) = 48804.d0   ! 010792
       T(31) = 49169.d0   ! 010793
       T(32) = 49534.d0   ! 010794
       T(33) = 50083.d0   ! 010196
       T(34) = 50630.d0   ! 010797
       T(35) = 51179.d0   ! 010199
       T(36) = 53736.d0   ! 010106
      
       DO  i = 37,50
       		T(i) = 53736.d0 + 20.D0 * 365.25d0  
       enddo
       
       TARTUF = 0.d0
       IF(D.LT.TA) THEN 
       WRITE(6,*) 'TAI is not defined before this date'
       GOTO 18
       TARTUF = 999.d0
       ENDIF
       
       DO 17 I=1,50     
	      IF(D.GE.T(I)) GOTO 17  
	      TARTUF = A(I) + (D-B(I))*C(I) 
 	      GOTO 18
   17  CONTINUE

       WRITE(6,*) ' TAI - UTC cannot be given after the date', T(i)
       TARTUF = 999.d0
   18  RETURN
       END         



c ########################################################
      SUBROUTINE iau_nut_ks ( EPOCH1, EPOCH2, DPSI, DEPS )
*+
*  - - - - - - - - - -
*   i a u _ N U T  _ K S 
*  - - - - - - - - - -
*
* Complement to SOFA by Ch. Bizouard, August 2008 
* 
*  Nutation, Rigid Earth Nutation of Kinoshita-Souchay model.
*
*  Status:  canonical model.
*
*  Given:
*     EPOCH1,EPOCH2   d      TDB starting epoch (Note 1)
*
*  Returned:
*     DPSI            d      nutation in longitude (radians)
*     DEPS            d      nutation in obliquity (radians)
*
*  Notes:
*
*  1) The epoch EPOCH1+EPOCH2 is a Julian Date, apportioned in any
*     convenient way between the two arguments.  For example,
*     JD(TDB)=2450123.7 could be expressed in any of these ways,
*     among others:
*
*            EPOCH1        EPOCH2
*
*         2450123.7D0        0D0        (JD method)
*          2451545D0      -1421.3D0     (J2000 method)
*         2400000.5D0     50123.2D0     (MJD method)
*         2450123.5D0       0.2D0       (date & time method)
*
*     The JD method is the most natural and convenient to use in
*     cases where the loss of several decimal digits of resolution
*     is acceptable.  The J2000 method is best matched to the way
*     the argument is handled internally and will deliver the
*     optimum resolution.  The MJD method and the date & time methods
*     are both good compromises between resolution and convenience.
*
*  2) The nutation components are with respect to the ecliptic of
*     date.
*
*  Called:
*     iau_ANP     normalize radians to range -pi to +pi
*
*
*
*-----------------------------------------------------------------------

      IMPLICIT NONE

      DOUBLE PRECISION EPOCH1, EPOCH2, DPSI, DEPS

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

*  2Pi
      DOUBLE PRECISION D2PI
      PARAMETER ( D2PI = 6.283185307179586476925287D0 )

*  Units of 0.1 milliarcsecond to radians
      DOUBLE PRECISION U2R
      PARAMETER ( U2R = DAS2R/1D4 )

*  Reference epoch (J2000), JD
      DOUBLE PRECISION DJ0
      PARAMETER ( DJ0 = 2451545D0 )

*  Days per Julian century
      DOUBLE PRECISION DJC
      PARAMETER ( DJC = 36525D0 )

      DOUBLE PRECISION T, EL, ELP, F, D, OM, DP, DE, ARG, S, C, PSIC, EPSS
      INTEGER I, J

      DOUBLE PRECISION iau_ANPM


*  ------------------------------------------------
*  Table of multiples of arguments and coefficients
*  ------------------------------------------------
*
*  The coefficient values are in 0.1 mas units and the rates of change
*  are in mas per Julian millennium.

      REAL X(11,389)

*                Multiple of            Longitude        Obliquity
*           L    L'   F    D  Omega   coeff. of sin    coeff. of cos

      DATA ((X(I,J),I=1,11),J=1,10)/
     :  -2.,  0.,  0.,  2.,  1.,      0.3,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  0.,  0.,      0.2,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0.,  1., -1.,  1.,      1.2,    0.0,    -0.6,  0.0,     0.6,     0.3,
     :  -1.,  1.,  0.,  1.,  2.,     -0.5,    0.0,     0.2,  0.0,     0.0,     0.0,
     :  -1.,  0.,  1.,  0.,  3.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0., -1.,  1.,  0.,      0.3,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  1.,     -0.5,    0.0,     0.2,  0.0,    -0.4,    -0.2,
     :  -1.,  0.,  0.,  1.,  0.,      0.0,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0., -2.,  2.,  0.,      0.0,    0.0,     0.0,  0.0,     0.2,     0.0,
     :  -2.,  0.,  0.,  2.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J= 11, 20)/
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  0., -2.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  0.,  2.,  0.,     -0.4,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :  -1.,  0.,  0.,  1.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0., -1.,  1.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  1.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  1., -1.,  1.,      0.2,    0.0,    -0.1,  0.0,     0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  1., -1.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J= 21, 30)/
     :   0.,  0.,  0.,  0.,  1.,      0.1,    0.0,     0.0,  0.0,     0.1,     0.0,
     :  -1.,  0.,  0.,  0.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  1.,      0.1,    0.0,     0.0,  0.0,     0.3,     0.2,
     :   0.,  0.,  0.,  0.,  1.,-172806.0,-1742.8, 92279.1, 90.8,     1.4,    -0.3,
     :   0., -2.,  2., -2.,  1.,     -1.3,    0.0,     0.7,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.3,     0.2,
     :   0.,  0.,  0.,  0.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2.,  0.,  0.,  2.,  0.,     -0.4,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  0.,  0.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J= 31, 40)/
     :   0.,  0.,  0.,  0.,  1.,      0.0,    0.0,     0.0,  0.0,     0.1,     0.0,
     :  -1.,  1.,  0.,  1.,  1.,      1.3,    0.0,    -0.7,  0.0,     0.0,     0.0,
     :  -1.,  0.,  1.,  0.,  2.,      0.0,    0.0,     0.0,  0.0,     0.4,     0.2,
     :   0.,  0.,  1., -1.,  1.,      0.3,    0.0,    -0.1,  0.0,     0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.5,    0.0,     0.0,  0.0,     0.4,     0.0,
     :   0.,  0.,  1., -1.,  1.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.1,     0.1,
     :   0.,  0.,  1., -1.,  1.,     -0.1,    0.0,     0.0,  0.0,    -0.1,    -0.1,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -2.,  1.,      0.0,    0.0,     0.0,  0.0,    -0.4,    -0.2/
      DATA ((X(I,J),I=1,11),J= 41, 50)/
     :  -2.,  0.,  0.,  2.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  1., -1.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.2,     0.1,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.2,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0.,  1., -1.,  1.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :  -2.,  0.,  0.,  2.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  2.,   2090.1,    2.0,  -903.4,  4.7,     0.1,     0.0,
     :   0.,  2., -2.,  2.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  0., -2.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J= 51, 60)/
     :  -1.,  1.,  0.,  1.,  0.,      3.3,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  1.,  0.,  1.,      0.0,    0.0,     0.0,  0.0,    -1.1,    -0.9,
     :   0.,  0.,  0.,  0.,  0.,     -0.2,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      2.1,    0.0,     1.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.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.3,     0.1,
     :   0.,  0.,  1., -1.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  3.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J= 61, 70)/
     :  -1.,  0.,  1.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,     0.3,     0.0,
     :  -1., -1.,  2., -1.,  1.,      0.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.3,     0.1,
     :   0.,  0.,  1., -1.,  1.,      0.0,    0.0,     0.0,  0.0,    -0.1,    -0.1,
     :   0.,  0.,  0.,  0.,  0.,     -1.2,    0.0,     0.5,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -2.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -2.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2.,  2.,  0.,  2.,  0.,     -0.2,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2.,  0.,  2.,  0.,  2.,     -3.1,    0.0,     1.3,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,     0.6,     0.0/
      DATA ((X(I,J),I=1,11),J= 71, 80)/
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  2.,  0.,  1.,     45.1,    0.5,   -24.1, -0.1,     0.0,     0.0,
     :  -2.,  0.,  2.,  0.,  0.,    -11.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.1,     0.0,
     :   2.,  0., -2.,  0.,  1.,      0.8,    0.0,    -0.4,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.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.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.5,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,     0.1,    -0.1/
      DATA ((X(I,J),I=1,11),J= 81, 90)/
     :   0.,  0.,  0.,  0.,  0.,      1.5,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.1,    -0.1,
     :  -2., -1.,  2.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2., -1.,  2.,  0.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.5,    0.0,     0.2,  0.0,     0.0,     0.0,
     :  -2.,  1.,  0.,  2.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.2,    0.0,     0.0,  0.0,     0.3,     0.0,
     :   1.,  0.,  0., -1.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J= 91,100)/
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,     0.3,     0.0,
     :  -1.,  0.,  0.,  1.,  0.,      2.8,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  1.,  0.,     -0.2,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  1.,  0.,  0.,  2.,      0.7,    0.0,    -0.3,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -1.2,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.4,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  1.,  1.,      0.6,    0.0,    -0.3,  0.0,     0.0,     0.0,
     :   0.,  1.,  0.,  0.,  1.,    -14.8,   -0.2,     7.9,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=101,110)/
     :   0.,  1.,  0.,  0.,  0.,    -17.2,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  1.,  0.,  0.,  0.,   1272.3,  -30.9,    -1.4,  0.0,     0.0,     0.0,
     :   0., -1.,  2., -2.,  2.,      0.5,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :   0., -1.,  2., -2.,  2.,    212.6,   -5.1,   -92.1,  2.8,     0.0,     0.0,
     :   0., -1.,  0.,  0.,  1.,     -9.8,    0.1,     5.2,  0.1,     0.0,     0.0,
     :   0., -1.,  2., -2.,  1.,     -4.8,   -0.1,     2.6, -0.1,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0., -1.,  0.,  0.,  2.,      0.5,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :   0.,  1., -2.,  2.,  0.,     -0.8,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=111,120)/
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2., -1.,  1.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   1.,  0., -2.,  1.,  0.,     -0.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.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.6,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2.,  1.,  2.,  0.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.1,     0.0/
      DATA ((X(I,J),I=1,11),J=121,130)/
     :   2.,  0.,  0., -2.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  0., -2.,  1.,      4.0,    0.0,    -2.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.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.1,     0.1,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  0.,  2.,  0.,    -44.9,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  0.,  2.,  1.,     -5.2,   -0.1,     2.8,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.4,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=131,140)/
     :  -2.,  0.,  0.,  2.,  2.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  2.,  0.,  0.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -2.,  3.,      1.2,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,    -0.1,  0.0,    -0.4,    -0.2,
     :   0.,  2.,  0.,  0.,  0.,     -0.2,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  2.,  0.,  0.,  0.,     15.9,   -0.8,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -2.,  2.,    -21.9,    0.0,     9.5, -0.1,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=141,150)/
     :   0.,  0.,  2., -2.,  2., -12751.9,  -12.9,  5523.8,-28.6,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.1,  0.0,    -0.4,     0.2,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.1,  0.0,    -0.1,     0.0,
     :   0., -2.,  0.,  0.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -2.,  1.,    125.1,    1.3,   -66.3, -0.1,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0., -2.,  2.,  0.,     20.6,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=151,160)/
     :   0.,  0., -2.,  2.,  1.,      0.8,    0.0,    -0.4,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.2,    -0.1,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.3,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   2.,  1.,  0., -2.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2., -1.,  0.,  2.,  0.,     -1.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2., -1.,  0.,  2.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,    -0.1,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,    -0.1,    -0.1/
      DATA ((X(I,J),I=1,11),J=161,170)/
     :   1.,  0.,  2., -3.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.5,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  1.,  2., -2.,  2.,     -0.6,    0.0,     0.3,  0.0,     0.0,     0.0,
     :   0.,  1.,  2., -2.,  2.,   -498.8,   11.9,   215.9, -6.5,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  1.,  2., -2.,  1.,      3.5,    0.0,    -1.9,  0.0,     0.0,     0.0,
     :   0., -1., -2.,  2.,  0.,      0.6,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.4,    0.0,    -0.2,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=171,180)/
     :   2.,  0.,  2., -4.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  2.,  2., -2.,  2.,    -15.2,    0.7,     6.6, -0.4,     0.0,     0.0,
     :   0.,  2.,  2., -2.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  3.,  2., -2.,  2.,     -0.4,    0.0,     0.2,  0.0,     0.0,     0.0,
     :  -1.,  0., -2.,  4.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  2., -4.,  1.,     -0.4,    0.0,     0.2,  0.0,     0.0,     0.0,
     :  -1., -2.,  0.,  2.,  0.,      0.3,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  2., -4.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -3.,  0.,  0.,  4.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=181,190)/
     :  -1.,  1.,  2., -2.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -3.,  2.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :  -1., -1.,  0.,  2.,  1.,      0.7,    0.0,    -0.4,  0.0,     0.0,     0.0,
     :  -1., -1.,  0.,  2.,  0.,      7.0,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   1.,  1.,  0., -2.,  1.,     -0.6,    0.0,     0.3,  0.0,     0.0,     0.0,
     :  -2.,  0.,  0.,  3.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0., -2.,  2.,  0.,     -0.6,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2., -2.,  1.,     -1.9,    0.0,     1.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2., -2.,  2.,      0.3,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   1., -2.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=191,200)/
     :  -1.,  0.,  0.,  2.,  2.,     -0.4,    0.0,     0.2,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  2.,  1.,     14.7,    0.1,    -7.7,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  2.,  0.,    149.6,    0.1,    -1.9,  0.0,     0.0,     0.0,
     :   1.,  0.,  0., -2.,  1.,    -12.1,   -0.1,     6.6,  0.0,     0.0,     0.0,
     :   1.,  0.,  0., -2.,  2.,      0.3,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   1., -1.,  0.,  0.,  1.,      0.5,    0.0,    -0.3,  0.0,     0.0,     0.0,
     :   1., -1.,  0.,  0.,  0.,      4.5,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -1.,  1.,  0.,  0.,  1.,     -0.3,    0.0,     0.2,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  1.,  1.,     -0.4,    0.0,     0.2,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  1.,  0.,     -3.8,    0.0,     0.1,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=201,210)/
     :   0.,  0.,  0.,  1.,  0.,     -0.2,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0., -1.,  1.,      0.3,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :  -1.,  1.,  0.,  2.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :  -1.,  1.,  0.,  2.,  0.,     -0.9,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1., -1.,  2.,  0.,  2.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   1., -1.,  0., -2.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1., -1.,  2.,  0.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -3.,  0.,  2.,  0.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  0.,  0.,  2.,     -1.9,    0.0,     0.8,  0.0,     0.0,     0.0,
     :   1.,  0.,  0.,  0.,  1.,     61.1,    0.6,   -31.9,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=211,220)/
     :   1.,  0.,  0.,  0.,  0.,    677.7,    0.7,    -9.7,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  0.,  1.,    -54.3,   -0.6,    29.7,  0.0,     0.0,     0.0,
     :   0.,  1.,  0.,  1.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  0.,  2.,      1.3,    0.0,    -0.6,  0.0,     0.0,     0.0,
     :   0.,  1.,  0.,  1.,  0.,      0.5,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0., -1.,  2., -1.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  1.,  0.,  1.,      0.0,    0.0,     0.0,  0.0,    -0.2,    -0.1,
     :   0.,  0.,  1.,  0.,  0.,      0.0,    0.0,     0.0,  0.0,     0.1,     0.0,
     :  -1.,  0.,  2.,  0.,  2.,    119.2,    0.1,   -51.4,  0.3,     0.0,     0.0,
     :  -1.,  2.,  0.,  2.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=221,230)/
     :  -1.,  0.,  2.,  0.,  1.,     19.8,    0.2,   -10.3,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2.,  0.,  0.,     -3.9,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   1.,  0., -2.,  0.,  1.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   1.,  1.,  0.,  0.,  1.,     -0.3,    0.0,     0.2,  0.0,     0.0,     0.0,
     :   1.,  1.,  0.,  0.,  0.,     -3.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   1., -1.,  2., -2.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1., -1.,  0.,  0.,  1.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -1.,  2.,     -0.7,    0.0,     0.3,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -1.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :  -1.,  1.,  2.,  0.,  2.,      0.4,    0.0,    -0.2,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=231,240)/
     :  -1.,  1.,  2.,  0.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   1.,  0.,  2., -2.,  2.,     27.7,    0.0,   -11.9,  0.1,     0.0,     0.0,
     :   1.,  0.,  2., -2.,  1.,      5.6,    0.1,    -2.9,  0.0,     0.0,     0.0,
     :  -1.,  0., -2.,  2.,  0.,      0.7,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  4., -2.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  1.,  2., -2.,  2.,      1.2,    0.0,    -0.5,  0.0,     0.0,     0.0,
     :   1.,  1.,  2., -2.,  1.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -2., -1.,  0.,  4.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -4.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0., -2.,  0.,  2.,  0.,      0.2,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=241,250)/
     :  -2.,  0.,  0.,  4.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  0.,  4.,  0.,      1.2,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  0., -4.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0., -1.,  0.,  2.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0., -1.,  0.,  2.,  0.,      4.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  1.,  0., -2.,  1.,     -0.3,    0.0,     0.2,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  3.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2., -1.,  2.,  2.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -2.,  0.,  2., -2.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  2.,  2.,     -0.4,    0.0,     0.2,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=251,260)/
     :   0.,  0.,  0.,  2.,  1.,     -6.2,   -0.1,     3.2,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  2.,  0.,     60.4,    0.1,    -1.6,  0.0,     0.0,     0.0,
     :   0., -2.,  2.,  0.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0., -2.,  1.,     -4.6,   -0.1,     2.6,  0.0,     0.0,     0.0,
     :   0.,  0.,  0., -2.,  2.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  2.,  2.,  2.,      1.4,    0.0,    -0.6,  0.0,     0.0,     0.0,
     :  -2.,  0.,  2.,  2.,  1.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  2.,  2.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2., -1.,  0.,  0.,  0.,      0.4,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  0.,  1.,  0.,     -0.3,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=261,270)/
     :   0.,  1.,  0.,  2.,  1.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0., -1.,  2.,  0.,  2.,     -6.9,    0.2,     3.0,  0.2,     0.0,     0.0,
     :   0.,  1.,  0.,  2.,  0.,     -0.6,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0., -1.,  2.,  0.,  1.,     -0.6,    0.0,     0.3,  0.0,     0.0,     0.0,
     :   0., -1.,  0., -2.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2.,  1.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  0.,  0.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  0.,  0.,  1.,      2.1,    0.0,    -1.1,  0.0,     0.0,     0.0,
     :   2.,  0.,  0.,  0.,  0.,     27.9,    0.0,    -0.8,  0.0,     0.0,     0.0,
     :  -2.,  0.,  0.,  0.,  1.,     -2.1,    0.0,     1.2,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=271,280)/
     :  -2.,  0.,  0.,  0.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  1.,  0.,  1.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  2.,  0.,  3.,      0.2,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  2.,  0.,  2.,  -2215.3,   -2.0,   949.3, -4.6,     0.0,     0.0,
     :  -1.,  0.,  2.,  0.,  2.,     -0.1,    0.0,     0.1,  0.0,    -0.1,     0.0,
     :   1.,  0.,  2.,  0.,  2.,      0.1,    0.0,    -0.1,  0.0,    -0.1,     0.0,
     :   0.,  0.,  2.,  0.,  1.,   -378.6,   -3.5,   194.2,  0.2,     0.0,     0.0,
     :   0.,  0.,  2.,  0.,  0.,     24.7,    0.0,    -0.7,  0.0,     0.0,     0.0,
     :   0.,  0., -2.,  0.,  1.,     -0.6,    0.0,     0.3,  0.0,     0.0,     0.0,
     :  -2.,  0.,  4.,  0.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=281,290)/
     :   2.,  1.,  0.,  0.,  0.,     -0.3,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2., -1.,  2., -2.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  2., -1.,  2.,     -0.3,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  1.,  2.,  0.,  2.,      7.4,   -0.2,    -3.2, -0.2,     0.0,     0.0,
     :   0.,  1.,  2.,  0.,  1.,      0.8,    0.0,    -0.4,  0.0,     0.0,     0.0,
     :   2.,  0.,  2., -2.,  2.,      6.3,    0.0,    -2.7,  0.0,     0.0,     0.0,
     :   2.,  0.,  2., -2.,  1.,      1.0,    0.0,    -0.5,  0.0,     0.0,     0.0,
     :  -2.,  0., -2.,  2.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  2.,  2.,  0.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  4., -2.,  2.,      0.9,    0.0,    -0.4,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=291,300)/
     :   0.,  0.,  4., -2.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   2.,  1.,  2., -2.,  2.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -1., -1.,  0.,  4.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  4.,  1.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0.,  4.,  0.,      1.3,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -1., -2.,  2.,  2.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  0., -4.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   1., -1.,  0.,  2.,  0.,      0.5,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1., -1.,  2.,  2.,  2.,     -2.8,    0.0,     1.2,  0.0,     0.0,     0.0,
     :  -1., -1.,  2.,  2.,  1.,     -0.4,    0.0,     0.2,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=301,310)/
     :   1.,  0.,  0.,  2.,  1.,     -0.9,    0.0,     0.5,  0.0,     0.0,     0.0,
     :   1.,  0.,  0.,  2.,  0.,      6.3,    0.0,    -0.3,  0.0,     0.0,     0.0,
     :  -1.,  0.,  0., -2.,  1.,     -0.5,    0.0,     0.2,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2.,  2.,  2.,    -58.5,   -0.1,    25.0, -0.1,     0.0,     0.0,
     :  -1.,  0.,  2.,  2.,  1.,    -10.1,   -0.1,     5.1,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2.,  2.,  0.,      0.6,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1., -1.,  2.,  0.,  2.,     -2.8,    0.0,     1.2,  0.0,     0.0,     0.0,
     :   1.,  1.,  0.,  2.,  0.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1., -1.,  2.,  0.,  1.,     -0.4,    0.0,     0.2,  0.0,     0.0,     0.0,
     :   0.,  0.,  2.,  1.,  2.,      1.6,    0.0,    -0.7,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=311,320)/
     :   0.,  0.,  2.,  1.,  1.,      0.3,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -1.,  1.,  2.,  2.,  2.,      0.6,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :  -1.,  1.,  2.,  2.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   3.,  0.,  0.,  0.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   3.,  0.,  0.,  0.,  0.,      1.5,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :  -3.,  0.,  0.,  0.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   1.,  0.,  2.,  0.,  2.,   -295.8,   -0.3,   126.1, -0.6,     0.0,     0.0,
     :   1.,  0.,  2.,  0.,  1.,    -51.0,   -0.4,    25.7,  0.0,     0.0,     0.0,
     :   1.,  0.,  2.,  0.,  0.,      3.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  1.,  2.,  1.,  2.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=321,330)/
     :  -1.,  0., -2.,  0.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :  -1.,  0.,  4.,  0.,  2.,      1.1,    0.0,    -0.5,  0.0,     0.0,     0.0,
     :  -1.,  0.,  4.,  0.,  1.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   2.,  0.,  2., -1.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  1.,  2.,  0.,  2.,      2.4,    0.0,    -1.0,  0.0,     0.0,     0.0,
     :   1.,  1.,  2.,  0.,  1.,      0.4,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :   3.,  0.,  2., -2.,  2.,      0.9,    0.0,    -0.4,  0.0,     0.0,     0.0,
     :   3.,  0.,  2., -2.,  1.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   1.,  0.,  4., -2.,  2.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0., -1.,  0.,  4.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=331,340)/
     :  -2., -1.,  2.,  4.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  4.,  1.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  4.,  0.,      0.5,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0., -2.,  2.,  2.,  2.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  2.,  4.,  2.,     -1.2,    0.0,     0.5,  0.0,     0.0,     0.0,
     :  -2.,  0.,  2.,  4.,  1.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0., -1.,  2.,  2.,  2.,     -2.6,    0.0,     1.1,  0.0,     0.0,     0.0,
     :   0., -1.,  2.,  2.,  1.,     -0.4,    0.0,     0.2,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2.,  3.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  0.,  2.,  0.,      0.6,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=341,350)/
     :   2.,  0.,  0.,  2.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :  -2.,  0.,  0., -2.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  2.,  2.,  2.,    -38.1,    0.0,    16.2,  0.0,     0.0,     0.0,
     :   0.,  0.,  2.,  2.,  1.,     -6.6,   -0.1,     3.3, -0.1,     0.0,     0.0,
     :   0.,  0.,  2.,  2.,  0.,      0.4,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2., -1.,  2.,  0.,  2.,     -0.5,    0.0,     0.2,  0.0,     0.0,     0.0,
     :   1.,  0.,  2.,  1.,  2.,      0.3,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   2., -1.,  2.,  0.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  2.,  1.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  1.,  2.,  2.,  2.,      0.5,    0.0,    -0.2,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=351,360)/
     :   0.,  1.,  2.,  2.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   4.,  0.,  0.,  0.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  2.,  0.,  2.,    -30.7,    0.0,    13.0, -0.1,     0.0,     0.0,
     :   2.,  0.,  2.,  0.,  1.,     -5.3,    0.0,     2.6,  0.0,     0.0,     0.0,
     :   2.,  0.,  2.,  0.,  0.,      0.3,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  4.,  0.,  2.,      0.2,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   2.,  1.,  2.,  0.,  2.,      0.4,    0.0,    -0.2,  0.0,     0.0,     0.0,
     :   2.,  1.,  2.,  0.,  1.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   4.,  0.,  2., -2.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1., -1.,  2.,  4.,  2.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=361,370)/
     :   1.,  0.,  0.,  4.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2.,  4.,  2.,     -1.5,    0.0,     0.7,  0.0,     0.0,     0.0,
     :  -1.,  0.,  2.,  4.,  1.,     -0.3,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   1., -1.,  2.,  2.,  2.,     -0.6,    0.0,     0.2,  0.0,     0.0,     0.0,
     :   1., -1.,  2.,  2.,  1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   3.,  0.,  0.,  2.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   1.,  0.,  2.,  2.,  2.,     -7.7,    0.0,     3.2,  0.0,     0.0,     0.0,
     :   1.,  0.,  2.,  2.,  1.,     -1.3,    0.0,     0.7,  0.0,     0.0,     0.0,
     :   1.,  0.,  2.,  2.,  0.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   3., -1.,  2.,  0.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=371,380)/
     :   1.,  1.,  2.,  2.,  2.,      0.1,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   3.,  0.,  2.,  0.,  2.,     -2.9,    0.0,     1.2,  0.0,     0.0,     0.0,
     :   3.,  0.,  2.,  0.,  1.,     -0.5,    0.0,     0.2,  0.0,     0.0,     0.0,
     :   0., -1.,  2.,  4.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  2.,  4.,  2.,     -0.7,    0.0,     0.3,  0.0,     0.0,     0.0,
     :   0.,  0.,  2.,  4.,  1.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   2., -1.,  2.,  2.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   2.,  0.,  2.,  2.,  2.,     -1.1,    0.0,     0.5,  0.0,     0.0,     0.0,
     :   2.,  0.,  2.,  2.,  1.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   4.,  0.,  2.,  0.,  2.,     -0.3,    0.0,     0.1,  0.0,     0.0,     0.0/
      DATA ((X(I,J),I=1,11),J=381,389)/
     :   1.,  0.,  2.,  4.,  2.,     -0.2,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   3.,  0.,  2.,  2.,  2.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0., -3.,  0., -3.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0., -2.,  0., -2.,     -0.3,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0., -2.,  0., -1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0.,  0.,      0.4,    0.0,    -0.1,  0.0,     0.0,     0.0,
     :   0.,  0., -2.,  2., -2.,     -0.1,    0.0,     0.1,  0.0,     0.0,     0.0,
     :   0.,  0.,  0.,  0., -1.,     -0.1,    0.0,     0.0,  0.0,     0.0,     0.0,
     :   0.,  0.,  2., -2.,  2.,      0.1,    0.0,     0.0,  0.0,     0.0,     0.0/





* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*  Interval between fundamental epoch J2000.0 and given epoch (JC).
      T = ( ( EPOCH1-DJ0 ) + EPOCH2 ) / DJC

*
*  FUNDAMENTAL ARGUMENTS in the FK5 reference system
*

*  Mean longitude of the Moon minus mean longitude of the Moon's
*  perigee.
       EL = iau_ANPM ( ( 485866.733D0 + ( 715922.633D0 +
     :                  ( 31.310D0 + 0.064D0 * T ) * T ) * T ) * DAS2R
     :                + MOD(1325D0*T, 1D0) * D2PI )

*  Mean longitude of the Sun minus mean longitude of the Sun's perigee.
       ELP = iau_ANPM ( ( 1287099.804D0 + ( 1292581.224D0 +
     :                   ( -0.577D0 -0.012D0 * T ) * T ) * T ) * DAS2R
     :                 + MOD(99D0*T, 1D0) * D2PI )

*  Mean longitude of the Moon minus mean longitude of the Moon's node.
       F = iau_ANPM ( ( 335778.877D0 + ( 295263.137D0 +
     :                 ( -13.257D0 + 0.011D0 * T ) * T ) * T ) * DAS2R
     :               + MOD(1342D0*T, 1D0) * D2PI )

*  Mea n elongation of the Moon from the Sun.
        D = iau_ANPM ( ( 1072261.307D0 + ( 1105601.328D0 +
     :                 ( -6.891D0 + 0.019D0 * T ) * T ) * T ) * DAS2R
     :               + MOD(1236D0*T, 1D0) * D2PI )

*  Longitude of the mean ascending node of the lunar orbit on the
*  ecliptic, measured from the mean equinox of date.
       OM = iau_ANPM( ( 450160.280D0 + ( -482890.539D0 +
     :                 ( 7.455D0 + 0.008D0 * T ) * T ) * T ) * DAS2R
     :               + MOD( -5D0*T, 1D0) * D2PI )

*  ---------------
*  Nutation series
*  ---------------

*  Change time argument from centuries to millennia.
       T = T / 10D0
 
*  Initialize nutation components.
       DP = 0D0
       DE = 0D0

*  Sum the nutation terms, ending with the biggest.
       DO 1 J=389,1,-1

*     Form argument for current term.
         ARG = DBLE(X(1,J)) * EL
     :       + DBLE(X(2,J)) * ELP
     :       + DBLE(X(3,J)) * F
     :       + DBLE(X(4,J)) * D
     :       + DBLE(X(5,J)) * OM

*     Accumulate current nutation term.
        S = DBLE(X(6,J)) + DBLE(X(7,J)) * T
        C = DBLE(X(8,J)) + DBLE(X(9,J)) * T
       
* termes dephases  
	PSIC = DBLE(X(10,J)) 
        EPSS = DBLE(X(11,J)) 
        IF ( S .NE. 0D0 ) DP = DP + S * DSIN(ARG) 
        IF ( C .NE. 0D0 ) DE = DE + C * DCOS(ARG) 

        IF ( PSIC .NE. 0D0 ) DP = DP + PSIC * DCOS(ARG) 
        IF ( EPSS .NE. 0D0 ) DE = DE + EPSS * DSIN(ARG) 
	 
*       Next term.
 1      CONTINUE

*  Convert results from 0.1 mas units to radians.
        DPSI = DP * U2R
        DEPS = DE * U2R

        return
        end  

    
	include '../Sofa/20210512/f77/src/anpm.for'
	include '../Sofa/20210512/f77/src/sp00.for' 
	include '../Sofa/20210512/f77/src/s00.for'
	include '../Sofa/20210512/f77/src/s06.for'
	include '../Sofa/20210512/f77/src/s06a.for'
	include '../Sofa/20210512/f77/src/era00.for'

	include '../Sofa/20210512/f77/src/fad03.for'  
	include '../Sofa/20210512/f77/src/fal03.for'   
	include '../Sofa/20210512/f77/src/faf03.for'   	
	include '../Sofa/20210512/f77/src/fama03.for'  
	include '../Sofa/20210512/f77/src/fane03.for' 
	include '../Sofa/20210512/f77/src/fapa03.for'  
	include '../Sofa/20210512/f77/src/faur03.for'
	include '../Sofa/20210512/f77/src/fae03.for'  
	include '../Sofa/20210512/f77/src/faju03.for'  
	include '../Sofa/20210512/f77/src/falp03.for'  
	include '../Sofa/20210512/f77/src/fame03.for'  
	include '../Sofa/20210512/f77/src/faom03.for'  
	include '../Sofa/20210512/f77/src/fasa03.for'  
	include '../Sofa/20210512/f77/src/fave03.for'
	include '../Sofa/20210512/f77/src/anp.for'
	include '../Sofa/20210512/f77/src/pnm06a.for'
	include '../Sofa/20210512/f77/src/bp06.for'
	include '../Sofa/20210512/f77/src/bpn2xy.for'
	include '../Sofa/20210512/f77/src/pfw06.for'
	include '../Sofa/20210512/f77/src/nut06a.for'
	include '../Sofa/20210512/f77/src/fw2m.for'
	include '../Sofa/20210512/f77/src/pmat06.for'
	include '../Sofa/20210512/f77/src/numat.for'
	include '../Sofa/20210512/f77/src/tr.for'
	include '../Sofa/20210512/f77/src/rxr.for'
	include '../Sofa/20210512/f77/src/cr.for'
	include '../Sofa/20210512/f77/src/obl06.for'
	include '../Sofa/20210512/f77/src/ir.for'
	include '../Sofa/20210512/f77/src/rz.for'
	include '../Sofa/20210512/f77/src/rx.for'
	include '../Sofa/20210512/f77/src/rxp.for'	
	include '../Sofa/20210512/f77/src/cp.for'
	include '../Sofa/20210512/f77/src/nut00a.for'
	include '../Sofa/20210512/f77/src/cal2jd.for'
	include '../Sofa/20210512/f77/src/d2tf.for'
	include '../Sofa/20210512/f77/src/d2dtf.for'
	include '../Sofa/20210512/f77/src/jd2cal.for'
	include '../Sofa/20210512/f77/src/dat.for'
	include '../Sofa/20210512/f77/src/xy06.for'
	include '../Sofa/20210512/f77/src/xys06a.for'

	     
	
