C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/datesub.for
C
C======DATESUB.FOR with Sample Drivers.
C      COLLECTED AND PUT TOGETHER JANUARY 1972, H. D. KNOBLE .
C      ORIGINAL REFERENCES ARE CITED IN EACH ROUTINE.
C
      INTEGER YYYY,MM,DD,JD,WD,DDD,MMA,DDA,NDIFF,I
      INTEGER*2 YYYY2,MM2,DD2
C
C====IDAY IS A COMPANION TO CALEND; GIVEN A CALENDAR DATE, YYYY, MM,
C           DD, IDAY IS RETURNED AS THE DAY OF THE YEAR.
C           EXAMPLE: IDAY(1984,4,22)=113
      IDAY(YYYY,MM,DD)=3055*(MM+2)/100-(MM+10)/13*2-91
     ,                 +(1-(MOD(YYYY,4)+3)/4+(MOD(YYYY,100)+99)/100
     ,                 -(MOD(YYYY,400)+399)/400)*(MM+10)/13+DD
C
C====IZLR(YYYY,MM,DD) GIVES THE WEEKDAY NUMBER 0=SUNDAY, 1=MONDAY,
C      ... 6=SATURDAY.  EXAMPLE: IZLR(1970,1,1)=4=THURSDAY
      IZLR(YYYY,MM,DD)=MOD((13*(MM+10-(MM+10)/13*12)-1)/5+DD+77
     ,            +5*(YYYY+(MM-14)/12-(YYYY+(MM-14)/12)/100*100)/4
     ,            + (YYYY+(MM-14)/12)/400-(YYYY+(MM-14)/12)/100*2,7)
C
C  Compute date this year for changing clocks back to EST.
C  See program: estdst.f90
C
C  Is this a leap year? I.e. is 12/31/yyyy the 366th day of the year?
      CALL GETDAT(YYYY2,MM2,DD2)
C---GETDAT is builtin using some Compilers.
      YYYY=YYYY2
      IF(IDAY(YYYY,12,31).EQ.366) THEN
         WRITE(*,*) YYYY,' is a Leap Year'
      ELSE
         WRITE(*,*) YYYY,' is not a Leap Year'
      ENDIF
C
C  DAYSUB SHOULD RETURN: 1970, 1, 1, 4, 1
      CALL DAYSUB(JD(1970,1,1),YYYY,MM,DD,WD,DDD)
      IF(YYYY.NE.1970.OR.MM.NE.1.OR.DD.NE.1.OR.WD.NE.4.OR.DDD.NE.1)
     *  THEN
        WRITE(*,*)'DAYSUB Failed; YYYY,MM,DD,WD,DDD=',YYYY,MM,DD,WD,DDD
        STOP 1
      ENDIF
C  DIFFERENCE BETWEEN TWO SAME MONTHS AND DAYS OVER 1 LEAP YEAR IS 366.
      NDIFF=NDAYS(5,22,1984,5,22,1983)
      IF(NDIFF.NE.366) THEN
        WRITE(*,*) 'NDAYS FAILED; NDIFF=',NDIFF
      ELSE
C  RECOVER MONTH AND DAY FROM YEAR AND DAY NUMBER.
        CALL CALEND(YYYY,DDD,MMA,DDA)
        IF(MMA.NE.1.AND.DDA.NE.1) THEN
          WRITE(*,*) 'CALEND FAILED; MMA,DDA=',MMA,DDA
                                  ELSE
          WRITE(*,*) '** DATE MANIPULATION SUBROUTINES SIMPLE TEST OK.'
        END IF
      END IF
      STOP
      END
      SUBROUTINE CALEND(YYYY,DDD,MM,DD)
C=====CALEND WHEN GIVEN A VALID YEAR, YYYY, AND DAY OF THE
C             YEAR, DDD, RETURNS THE MONTH, MM, AND DAY OF THE
C             MONTH, DD.
C             SEE ACM ALGORITHM 398, TABLELESS DATE CONVERSION, BY
C             DICK STONE, CACM 13(10):621.
      INTEGER YYYY,DDD,MM,DD,T
      T=0
      IF(MOD(YYYY,4).EQ.0) T=1
C-------THE FOLLOWING STATEMENT IS NECESSARY IF YYYY IS LESS TNAN
C           1900 OR GREATER THAN 2100.
      IF(MOD(YYYY,400).NE.0.AND.MOD(YYYY,100).EQ.0) T=0
      DD=DDD
      IF(DDD.GT.59+T) DD=DD+2-T
      MM=((DD+91)*100)/3055
      DD=(DD+91)-(MM*3055)/100
      MM=MM-2
C-------MM WILL BE CORRECT IFF DDD IS CORRECT FOR YYYY.
      IF(MM.GE.1 .AND. MM.LE.12) RETURN
      WRITE(*,1) DDD
1     FORMAT('0$$$CALEND: DAY OF THE YEAR INPUT =',I11,
     ,       ' IS OUT OF RANGE.')
      STOP 8
      END
      SUBROUTINE CDATE(JD,YYYY,MM,DD)
C=====GIVEN A JULIAN DAY NUMBER, NNNNNNNN, YYYY,MM,DD ARE RETURNED AS
C              AS THE CALENDAR DATE. JD=NNNNNNNN IS THE JULIAN DATE
C              FROM AN EPOCK IN THE VERY DISTANT PAST.  SEE CACM
C              1968 11(10):657, LETTER TO THE EDITOR BY FLIEGEL AND
C              VAN FLANDERN.
C    EXAMPLE CALL CDATE(2440588,YYYY,MM,DD) RETURNS 1970 1 1 .
C
      INTEGER JD,YYYY,MM,DD,L,N
      L=JD+68569
      N=4*L/146097
      L=L-(146097*N + 3)/4
      YYYY=4000*(L+1)/1461001
      L=L-1461*YYYY/4+31
      MM=80*L/2447
      DD=L-2447*MM/80
      L=MM/11
      MM=MM + 2 - 12*L
      YYYY=100*(N-49) + YYYY + L
      RETURN
      END
      SUBROUTINE DAYSUB(JD,YYYY,MM,DD,WD,DDD)
C=====GIVEN JD, A JULIAN DAY # (SEE ASF JD), THIS ROUTINE
C        CALCULATES DD, THE DAY NUMBER OF THE MONTH; MM, THE MONTH
C        NUMBER; YYYY THE YEAR; WD THE WEEKDAY NUMBER, AND DDD
C        THE DAY NUMBER OF THE YEAR.
C        ARITHMETIC STATEMENT FUNCTIONS 'IZLR' AND 'IDAY' ARE TAKEN
C        FROM REMARK ON ALGORITHM 398, BY J. DOUGLAS ROBERTSON,
C        CACM 15(10):918.
C
C   EXAMPLE: CALL DAYSUB(2440588,YYYY,MM,DD,WD,DDD) YIELDS 1970 1 1 4 1.
C
      INTEGER JD,YYYY,MM,DD,WD,DDD
C
C====IZLR(YYYY,MM,DD) GIVES THE WEEKDAY NUMBER 0=SUNDAY, 1=MONDAY,
C      ... 6=SATURDAY.  EXAMPLE: IZLR(1970,1,1)=4=THURSDAY
C
      IZLR(YYYY,MM,DD)=MOD((13*(MM+10-(MM+10)/13*12)-1)/5+DD+77
     ,            +5*(YYYY+(MM-14)/12-(YYYY+(MM-14)/12)/100*100)/4
     ,            + (YYYY+(MM-14)/12)/400-(YYYY+(MM-14)/12)/100*2,7)
C
C====IDAY IS A COMPANION TO CALEND; GIVEN A CALENDAR DATE, YYYY, MM,
C           DD, IDAY IS RETURNED AS THE DAY OF THE YEAR.
C           EXAMPLE: IDAY(1984,4,22)=113
C
      IDAY(YYYY,MM,DD)=3055*(MM+2)/100-(MM+10)/13*2-91
     ,                 +(1-(MOD(YYYY,4)+3)/4+(MOD(YYYY,100)+99)/100
     ,                 -(MOD(YYYY,400)+399)/400)*(MM+10)/13+DD
C
      CALL CDATE(JD,YYYY,MM,DD)
      WD=IZLR(YYYY,MM,DD)
      DDD=IDAY(YYYY,MM,DD)
      RETURN
      END
      FUNCTION JD(YYYY,MM,DD)
      INTEGER YYYY,MM,DD
C              DATE ROUTINE JD(YYYY,MM,DD) CONVERTS CALENDER DATE TO
C              JULIAN DATE.  SEE CACM 1968 11(10):657, LETTER TO THE
C              EDITOR BY HENRY F. FLIEGEL AND THOMAS C. VAN FLANDERN.
C    EXAMPLE JD(1970,1,1)=2440588
      JD=DD-32075+1461*(YYYY+4800+(MM-14)/12)/4
     ,         +367*(MM-2-((MM-14)/12)*12)/12-3*
     ,         ((YYYY+4900+(MM-14)/12)/100)/4
      RETURN
      END
      FUNCTION NDAYS(MM1,DD1,YYYY1, MM2,DD2,YYYY2)
      INTEGER YYYY1,MM1,DD1,YYYY2,MM2,DD2
C=====NDAYS IS RETURNED AS THE NUMBER OF DAYS BETWEEN TWO
C              DATES; THAT IS  MM1/DD1/YYYY1 MINUS MM2/DD2/YYYY2,
C              WHERE DATEI AND DATEJ HAVE ELEMENTS MM, DD, YYYY.
C-------NDAYS WILL BE POSITIVE IFF DATE1 IS MORE RECENT THAN DATE2.
      NDAYS=JD(YYYY1,MM1,DD1)-JD(YYYY2,MM2,DD2)
      RETURN
      END