      program julian
      integer day,month,year
      real*8 date,gmst,jdate,bepoch,lyear,jd1900
      character*80 text,answer
      lyear=365.242198781D0
      jd1900=2415020.31352D0
      one=1.0
      zero=0.0
  10  call questn('JD->D -1-, D->JD -2-, BY->D -3-, 0=stop: ',answer)
      read(answer,*)icom
      goto (91,92,93)icom
      if(icom.eq.0)stop
      goto 10
  91  call questn('Enter Julian date! ',answer)
      read(answer,*,err=91,end=91)jdate     
  911 year=1950
      month=1
      day=1
      call julda(year,month,day,date,gmst)
      if(jdate.ge.date)then
       isign=1
      else
       isign=-1
      endif
      do i=1,100
       year=1950+(i-1)*isign
       call julda(year,month,day,date,gmst)
       ndays=int(jdate-date)
       if(abs(ndays).lt.365)goto 100
      enddo
  100 if(isign.lt.0)then
       year=year-1
       call julda(year,month,day,date,gmst)
       ndays=int(jdate-date)
      endif
      do i=1,ndays
       call nextday(day,month,year)
      enddo
      fracod=jdate-int(jdate)-0.5
      if(fracod.lt.0)fracod=fracod+1.0
      ihour=int(fracod*24.0)
      imin=int((fracod*24.0-float(ihour))*60.0)
      write(text,1)year,month,day,ihour,imin
    1 format('Year= ',i4,' month= ',i2,' day= ',i2,', 'i2,' UT ',i2)
      call putout(text)
      goto 10
   92 call questn('Enter year month day! ',answer)
      read(answer,*,end=92,err=92)year,month,day
      call julda(year,month,day,date,gmst)
      write(text,2)date
    2 format('Julian date (0 UT) is = ',f11.3)
      call putout(text)
      goto 10
   93 call questn('Enter Besselian epoch! ',answer)
      if(answer(1:1).eq.'n')goto 111
      read(answer,*,end=93,err=93)bepoch
      jdate=((bepoch-1900.D0)*lyear)+jd1900
      goto 911
  111 end 
      
      subroutine questn(prompt,answer)
      CHARACTER*(*) answer
      CHARACTER*(*) PROMPT
      IF (PROMPT.NE.' ') WRITE (*,10) PROMPT
 10   FORMAT(1x,A,$)
      READ (5, '(A)', END=1000) answer
 1000 CONTINUE
      end

        SUBROUTINE JULDA(YEAR,MONTH,DAY,DATE,GMST)
C
C               Returns DATE = Julian Date and GMST = Greenwich Mean
C               sidereal time at 0 hrs U.T. on day DAY/MONTH/YEAR.
C               TJP 1977 August 5
C               NB Accuracy not tested; GMST is correct at least to
C               nearest second
C
        INTEGER DAY,MONTH,YEAR
        INTEGER MOFF(12)
        REAL*8 DATE,GMST,DU,SMD,T,DFLOAT,DMOD
        DATA MOFF/0,31,59,90,120,151,181,212,243,273,304,334/
C
C               JD number at 12 hrs UT on Jan 0 of year YEAR:
C               (Gregorian Calendar)
C
        NYRM1=YEAR-1
        IC=NYRM1/100
        JD=1721425 + 365*NYRM1 + NYRM1/4 - IC + IC/4
C
C               Number of days from Standard Epoch 1900 Jan 0.5
C               (JD 2415020.0) to Jan 0.0 of YEAR:
C
        DU=DFLOAT(JD-2415020) - 0.5D0
C
C               Day number; is it a leap year?
C
        NDAYS=MOFF(MONTH) + DAY
        IF(MONTH.GT.2 .AND. MOD(YEAR,4).EQ.0 .AND. MOD(YEAR,100).NE.0)
     1     NDAYS=NDAYS+1
C
C               UT from Epoch to today (days), (centuries)
C
        SMD=DU+DFLOAT(NDAYS)
        T=SMD/36525.D0
C
C               Greenwich Mean Sidereal Time
C
        GMST= DFLOAT(6*3600 + 38*60) +45.836D0
     1     + 8 640 184.542D0*T  +  0.0929D0*T**2
        GMST=DMOD(GMST,86400D0)
C
C
        DATE=2415020.D0+SMD
C
        RETURN
        END


C*PUTOUT -- write line on standard output
C+
      SUBROUTINE PUTOUT(TEXT)
      CHARACTER*(*) TEXT
C
C This subroutine writes one line on the standard output; the text to be
C written is supplied as a character-string argument.
C
C Argument:
C  TEXT   (input)  : character string for output.
C
C Subroutines required:
C  Fortran formatted I/O.
C
C History:
C  1987 Nov 11 - TJP
C-----------------------------------------------------------------------
      WRITE (*,'(1x,A)') TEXT
      END

      SUBROUTINE NEXTDAY ( IDAY, IMON, IYEAR )
C
C   Increments the date by one
C
      IMPLICIT UNDEFINED(A-Z)
      INTEGER *4 IDAY, IMON, IYEAR, DAYS(12,0:1), LEAP
      DATA DAYS / 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31,
     $            31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /
      LEAP = MIN ( 1, MOD( IYEAR, 4 ) )
      IDAY = IDAY + 1
      IF ( IDAY .GT. DAYS(IMON,LEAP) ) THEN
          IMON = IMON + 1
          IDAY = 1
          IF ( IMON .EQ. 13 ) THEN
              IMON = 1
              IYEAR = IYEAR + 1
          END IF
      END IF
      RETURN
      END

