!      Transforms a date expressed as "YEAR, DOY"  (Real*4) into 
!        - Modified Julian Day (MJD) at Oh UT (Real*8)
!        - Besselian years (Real*8)	
! 
!       f77 transdoy.f -o transdoy; rm transdoy.o 
!	
!                     Author: Bernard Guinot
!                     Contact: Martine Feissel-Vernier <mf.vernier@gmail.com>

	Implicit none

	Real*8 transdates
	Real*4 AN,DJ

	type*, 'YEAR  DOY  ( Ex : 1988 366)'
	accept*, an, dj
	type*, transdates(an,dj,1), transdates(an,dj,2)

	Stop
	End

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	DoublePrecision FUNCTION transdates(an,dj,kn)

	Implicit none

	CHARACTER*4 MMOIS
	Real*8 AMJD,AB
	Real*4 AN1,AJ,H1,P1,S1,AN2,DNJ,H2,P2,S2
	Real*4 AN,DJ
	Integer*4 kn

        COMMON/DATAGE/AMJD,AB,AN1,MMOIS,AJ,H1,P1,S1,AN2,DNJ,H2,P2,S2

	an1=0.
	amjd=0.
	ab=0.
	an2=an
	dnj=dj
	
	CALL DATEUR

	if(kn.eq.1) transdates=AMJD
	if(kn.eq.2) transdates=AB

	Return
	End

!       Appel du SOUS PROGRAMME DATEUR (Auteur: B. Guinot)

c Description
c -----------   DATE COMMUNE : AN1,MMOIS,AJ,H1,P1,S1
c                              ^^^
c               AN ET NUMERO DE JOUR : AN2,DNJ,H2,P2,S2
c                                      ^^^
c               DATE JULIENNE MODIFIEE : AMJD (DECIMAL,DOUBLE PRECISION)
c                                        ^^^^
c               ANNEES DE BESSEL : AB (DECIMAL,DOUBLE PRECISION)
c                                  ^^
c Emploi
c-------
c        Declarations :      CHARACTER*4 MMOIS
c                            DOUBLE PRECISION AMJD,AB
c      COMMON/DATAGE/AMJD,AB,AN1,MMOIS,AJ,H1,P1,S1,AN2,DNJ,H2,P2,S2
c
c        Definition : definir la date dans l'une des 4 notations
c                     mettre a zero les 3 autres valeurs soulignees ci-dessus
c                                                        ^^^^^^^^^^
c        Appel : CALL DATEUR     .==> les trois autres notations sont disponibles
C=======================================================================C
C            SOUS PROGRAMME DATEUR
C
C            DATE COMMUNE : AN1, MMOIS, AJ, H1, P1, S1
C            AN ET NUMERO DE JOUR : AN2, DNJ, H2, P2, S2
C            DATE JULIENNE MODIFIEE : AMJD(DECIMAL, DOUBLE PRECISION)
C            ANNEES DE BESSEL : AB(DECIMAL, DOUBLE PRECISION)
C=======================================================================C
	subroutine dateur
!                                        (Auteur: B. Guinot)
	
	CHARACTER*4 MOISA, MOISB, MOIS
	
	DOUBLE PRECISION AMJD0,ANNEE
	
	! zc avec pg. appelant
	DOUBLE PRECISION AMJD, AB
	CHARACTER*4 MMOIS
	COMMON/DATAGE/ AMJD, 
     @                 AB, 
     @                 AN1, 
     @                 MMOIS, 
     @                 AJ, 
     @                 H1, 
     @                 P1, 
     @                 S1, 
     @                 AN2, DNJ, H2, P2, S2
        
        DIMENSION MOISA(365),MOISB(366),MOIS(366)
	
        DATA MOISA/ 31*'JAN ',28*'FEB ',31*'MAR ',30*'APR ',
     @              31*'MAY ',30*'JUN ',31*'JUL ',31*'AUG ',
     @              30*'SEP ',31*'OCT ',30*'NOV ',31*'DEC '/,
     @       MOISB/ 31*'JAN ',29*'FEB ',31*'MAR ',30*'APR ',
     @              31*'MAY ',30*'JUN ',31*'JUL ',31*'AUG ',
     @              30*'SEP ',31*'OCT ',30*'NOV ',31*'DEC '/
C
C            CONVERSION BESSEL EN MJD ET VICE VERSA

C PROGRAMME DATEUR MODIFIE EN NOVEMBRE 1999 PAR D. GAMBIS


C   CORRESPONDANCE DES DATES INITIALES
C
C     AB = 1750        AMJD = 360 233.4836
C     AB = 1800        AMJD = 378 495.5935
C     AB = 1850        AMJD = 396 757.7035
C     AB = 1900        AMJD = 415 019.8134
C     AB = 1950        AMJD = 433 281.9234
C     AB = 2000        AMJD = 451 544.0333


	ANNEE=365.242199D0
	AB0=1950.
	AMJD0=33281.9234
	
	 
c	print*,'AN1,AMJD,ab,an2',AN1,AMJD,ab,an2
      IF (AN1.GT.1.) GO TO 100
      IF (AN2.GT.1.) GO TO 202
      IF (AMJD.GT.1.D0) GO TO 300
C
      AMJD = AMJD0+ANNEE*(AB-AB0)
      GO TO 450
  300 continue
  

  	AB = AB0+(AMJD-AMJD0)/ANNEE
C
C            CONVERSION MJD EN AN ET NUMERO DE JOUR
C
  450 MJD = AMJD
      MJDZ = 33282
C  .........................    DO 1 I=1,50
      DO 1 I=1,500
      INCR = 365
      J = I+1
      IF(MOD(J,4).EQ.0 ) INCR = 366
      MJDZ = MJDZ + INCR
      IF ((MJD-MJDZ).LT.0) GO TO 2
    1 CONTINUE
2	CONTINUE
    
	NAN = IFIX(AB0)+I-1
      AN2 = FLOAT(NAN)
      NDNJ = MJD-MJDZ+INCR+1


      DNJ = FLOAT(NDNJ)
      FH = (AMJD - FLOAT(MJD))*24.
      H2 = FLOAT (IFIX(FH))
      FM = (FH - H2)*60.
      P2 = FLOAT (IFIX(FM))
      S2 = (FM -P2)*60.
c	TYPE*,'NAN,AN2,NJDN,amjd', NAN,AN2,NDNJ,amjd

C
C            CONVERSION AN ET NUMERO DE JOUR EN DATE COMMUNE
C
  202 NAN = AN2
      NDNJ = IFIX(DNJ)
  200 IF (MOD(NAN,4).EQ.0 .AND.NAN.NE.1800.AND.NAN.NE.1900) GO TO 20
      DO 21 I = 1,365
   21 MOIS(I) = MOISA(I)
      GO TO 22
   20 DO 23 I = 1,366
   23 MOIS(I) = MOISB(I)
   22 I = 0
      J = 0
   30 I = I+1
      J = J+1
      IF (J.GE.28) GO TO 31
      GO TO 32
   31 IF(MOIS(I).NE.MOIS(I-1)) GO TO 33
      GO TO 32
   33 J = 1
   32 AJ = FLOAT(J)
      IF (I.EQ.NDNJ) GO TO 35
      GO TO 30
   35 AN1 = AN2
      MMOIS = MOIS(I)
      H1 = H2
      P1 = P2
      S1 = S2
      IF (AB.NE.0.) GO TO 1000
C
C            CONVERSION AN ET NUMERO DE JOUR EN MJD ET
C            ANNEE DE BESSEL
C
  150 I = IFIX(AN2 - 1950)
      MJD = I*365+(I+1)/4+33281
      AMJD = DFLOAT(MJD)+DBLE(DNJ)+H2/24.D0+P2/1440.D0+S2/86400.D0
c      TYPE*,'AN2,AB0,I,MJD,DNJ,AMJD,ANNEE',AN2,AB0,I,MJD,DNJ,AMJD,ANNEE
      AB = AB0+(AMJD-AMJD0)/ANNEE   
      GO TO 1000
C
C            CALCUL DU NUMERO DE JOUR D'APRES LA DATE COMMUNE
C
  100 NAN = AN1

      IF (MOD(NAN,4).EQ.0 .AND.NAN.NE.1800.AND.NAN.NE.1900) GO TO 60
      DO 50 I = 1,365
      IF (MOISA(I).EQ.MMOIS) GO TO 52
   50 CONTINUE
   60 DO 51 I = 1,366
      IF (MOISB(I).EQ.MMOIS) GO TO 52
   51 CONTINUE
   52 DNJ = FLOAT (I) - 1.+AJ
      AN2 = AN1
      H2 = H1
      P2 = P1
      S2 = S1
      GO TO 150
 1000 RETURN
      END
