SUBROUTINE yyyydoy_inc(yyyydoy)

c PURPOSE
c   Accepts 7-character yyyydoy and returns the next day in sequence
c   This algorithm uses all the usual fortran tricks with bits & bytes,
c   equivalences, etc. etc, and is y2k-friendly. 
c   Note: 2000 is a leap year, because of the '400' rule, while
c   1900 and 2100 are not.

c CALLING PARAMETER
c   yydoy is 5-byte character, input and output.

c HISTORY
c   writen 12/97 Edward A. Celarier, Software Corp of America
c

      CHARACTER*7 yyyydoy

      CHARACTER*4 doy
      INTEGER*4   idoy /0/
      BYTE        bdoy(4)
      EQUIVALENCE(doy, idoy, bdoy)

      CHARACTER*4 yyyy
      INTEGER*4   iyyyy
      BYTE        byyyy(4)
      EQUIVALENCE(yyyy, iyyyy, byyyy)

      INTEGER*4 d364 /'00030604'x/, mask /'0000000f'x/
      INTEGER*4 maskddd /'000f0f0f'x/
      INTEGER*2 m3 /'0003'x/
      INTEGER*4 w, ww

      doy= yyyydoy(4:7)

      w= AND(idoy , d364)
      w= XOR(w, d364)           !!! w is 0 if doy is 364, 365, or 366
      ww= AND(idoy, mask)
      
      IF(w .NE. 0) THEN
        IF(ww .NE. 9) THEN
          bdoy(4)=bdoy(4)+1
          yyyydoy(7:7)=doy(4:4)
          RETURN
        END IF
        idoy= AND(idoy,maskddd)
        ww= 100*bdoy(2) + 10*bdoy(3) + 10
        WRITE(yyyydoy(5:7),'(i3.3)') ww
        RETURN
      END IF

      IF(ww .EQ. 4) THEN
        bdoy(4)= bdoy(4)+1
        yyyydoy(7:7)=doy(4:4)
        RETURN
      END IF

      read(yyyydoy(1:4), '(i4)') iyyyy
      IF(ww .EQ. 5) THEN
       
        IF((AND(iyyyy, m3) .EQ. 0)) THEN
          yyyydoy(7:7)='6'
          RETURN
        END IF

      END IF

      WRITE(yyyydoy(1:4), '(i4.4)') iyyyy+1
      yyyydoy(5:7)='001'
      RETURN
      END