C   PROGRAM BARTELDATES
C
C.. A program to calculate the start day of each Bartel period,
C   beginning at Jan. 16, 1972 for the number of periods entered.
C   N. Schofield   3/17/1992
c   3/18/99 - new YDMD subroutine added to correct calculation of month
c   and day for leap years. Y2K compliance verified. 
c   Changed output record format.     J. F. Cooper  3/17/1999
C..Initialize
      IYR=72
      NTOTDAYS = -11
      ICENT=19
      WRITE(6,*), ' Enter number of Bartel periods for which dates ',
     *   'are desired.'
      READ (5,'(I6)') LIMIT
C
C..Start loop
      DO 200 I=0,LIMIT
      LYR=0
      NTOTDAYS=NTOTDAYS+27
      NBARTEL=1894+I
C
C..Check for end of year
      IF (MOD(IYR,4).EQ.0) LYR=1
      IF (NTOTDAYS.LE.(365+LYR)) GO TO 100
      IYR=IYR+1
      NTOTDAYS=NTOTDAYS-(365+LYR)
      IF (IYR.LE.99) GO TO 100
      IYR=IYR-100
      ICENT=ICENT+1
C
C..Calculate day & month
  100 IDD=1000*IYR+NTOTDAYS
      CALL YDMD(IDD,IMD)
      IYEAR=IYR+ICENT*100
      IM=(IMD-IYR*10000)/100
      IDAY=MOD(IMD,100)
      WRITE (6,150) NBARTEL,IYEAR,NTOTDAYS,IM,IDAY
      WRITE (1,150) NBARTEL,IYEAR,NTOTDAYS,IM,IDAY
  150 FORMAT (' Bartel cycle = ',I6,' year =',I4,' d.o.y. =',I3,
     .   ' month =',I2,' d.o.m. =',I2)
  200 CONTINUE
      END
C      SUBROUTINE YMDD(IMD,IDD)
C      DIMENSION MON(13)
C      DATA MON/1,32,60,91,121,152,182,213,244,274,305,335,366/
C     Change date from YYMMDD to YYDDD
C      IY = IMD/10000
C      IM = (IMD-IY*10000)/100
C      ID = MOD(IMD,100)
C      IDD = IY*1000+MON(IM)+ID-1
C      LYR = MOD(IY+3,4)/3
C      IF(IM.GT.2) IDD=IDD+LYR
C     RETURN
C
      SUBROUTINE YDMD(IDD,IMD)
C..Change date from YYDDD to YYMMDD
      DIMENSION MON(13),MONLYR(13)
      DATA MON/1,32,60,91,121,152,182,213,244,274,305,335,366/
      DATA MONLYR/1,32,61,92,122,153,183,214,245,275,306,336,367/
      IY = IDD/1000
      ID = MOD(IDD,1000)
      LYR = MOD(IY+3,4)/3
      DO 20 M=2,13
        M1=M
      IF (ID.LT.MON(M).AND.LYR.EQ.0) GO TO 30
   20 IF (ID.LT.MONLYR(M).AND.LYR.EQ.1) GO TO 30
   30 IM = M1-1
      IF(LYR.EQ.0) IDY = ID-MON(IM)+1
      IF(LYR.EQ.1) IDY = ID-MONLYR(IM)+1
      IMD = IY*10000+IM*100+IDY
      RETURN
      END