	implicit none
	
	real*8 mjd,mjd0,mjd1,tdj,s,tartuf
	integer jour,mois,an,nh,nm,i,j
	
c	write(6,*) 'mjd='
c
c	read(5,*) mjd

c	mjd0= 49533.d0 
c	mjd0 = mjd0 + tartuf(mjd0)/86400.d0 ! 29/6/1994
	
c	do i=1,10
c	mjd = mjd0 + dble(i-1) * 1.d0/2.d0 ! TAI 

	mjd0= 53736.d0 - 2.d0 / 86400.d0 
	mjd0 = mjd0 + tartuf(mjd0)/86400.d0 ! 29/6/1994
	
	do i=1,10
	mjd = mjd0 + dble(i-1) * 1.d0/86400.d0 ! TAI 

c	mjd0= 49534.d0 - 2.d0 / 86400.d0 
c	mjd0 = mjd0 + tartuf(mjd0)/86400.d0 ! 29/6/1994
	
c	do i=1,10
c	mjd = mjd0 + dble(i-1) * 1.d0/86400.d0 ! TAI 
		
	call MJDDATE(mjd,jour,mois,an,nh,nm,s)  ! entree TAI / sortie UTC
	
c	write(6,1000) tdj,jour,mois,an,nh,nm,s
	write(6,1100) mjd,jour,mois,an,nh,nm,s
	
	enddo
	
1000	format('tdj=',f16.8,/,'jour=',i3,/,'mois=',i3,/,'an=',i5,/,'nh=',i3,/,'nm=',i3,/,'s=',f5.2)			
1100	format('tdj=',f16.8,' jour=',i3,' mois=',i3,' an=',i5,' nh=',i3,' nm=',i3,' s=',f7.3)			

	return
	end
		
C=======================================================================C

      subroutine MJDDATE (mjd,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
      real*8 tdj,s,a,b,c,d,e,z,z1,t,t1,f,alpha,mjd,tartuf,mjd_tai,mjd_utc,tdj_tai,tdj_utc
      real*8 mjd_intercal_tai,mjd_intercal_utc,sintercal 
      integer*4 jour,mois,an,nh,nm,nsint,i
* 
* --- Initialisation --------------------------------------------------- 
* 
       jour = 0 
       mois = 0 
       an   = 0 
       nh   = 0 
       nm   = 0       
       s    = 0.d0  

* prise en compte des secondes intercalaires 
 	 
       mjd_tai = mjd 
    
* --- Conversion ------------------------------------------------------- 
* 
       t  = mjd_tai + 0.5d0 + 2400000.5
       z  = dint(t)
       t1 = mjd_tai 
       z1  = dint(t1) 
       f  = t1 - z1 
* 
       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 = int(f) 
       f  = (f-dble(nh))*60.d0 
       nm = int(f) 
       f  = (f-dble(nm))*60.d0 
       s  = f 

        nsint  = tartuf(mjd) 
	
	do i = 1,1000
	mjd_utc = mjd -  dble(nsint) /86400.d0 
	nsint = tartuf(mjd_utc) 
	enddo
	write(6,1020) nsint/86400.d0
1020	format('nsint (jours)=',f19.8)	
      ! calcul de la date civile (necessairement en UTC) 
        s = s - nsint
      	if (s.le.-1d-4) then
           
	if (tartuf(mjd_utc).eq.tartuf(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

  
* --- Sortie ----------------------------------------------------------- 
* 
      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 REAL*8(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

