!
! Date (year-day) manipulation utilities:
!
      integer function ichiyd(iyd)
      implicit none
      integer,intent(in) :: iyd
      integer :: ndmon(12) =
c                 J  F  M  A  M  J  J  A  S  O  N  D
     |          (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: iyr,iday,leap
c
c Check validity of yyddd iyd:
! 4/01: change to yyyyddd iyd (e.g., 2001080 for day 80 of 2001):
c
      ichiyd = 0
      iyr = iyd/1000
!     if (iyr.lt.1.or.iyr.gt.99) then
      if (iyr.lt.1900.or.iyr.gt.2100) then
        write(6,"('>>> ichiyd: bad year=',i4)") iyr
        ichiyd = 1
        return
      endif
      iday = iyd-iyr*1000
      leap = 0
      ndmon(2) = 28
      if (mod(iyr,4).eq.0) then
        leap = 1
        ndmon(2) = 29
      endif
      if (iday.lt.1.or.(leap.le.0.and.iday.gt.365).or.
     +                 (leap.gt.0.and.iday.gt.366)) then
        write(6,"('>>> ichiyd: bad iday=',i5,' (iyr=',i2,')')")
     +    iday,iyr
        ichiyd = 1
        return
      endif
      return
      end
!-----------------------------------------------------------------------
      integer function inciyd(iyd,ndays)
      implicit none
!
! Args:
      integer,intent(in) :: iyd,ndays
!
! Local:
      integer :: i,leap,iyr,iday,inc
      integer,external :: ichiyd
c
c Increment iyd (yyddd) by ndays (ndays may be +/-):
c
      if (ndays.eq.0) then
        inciyd = iyd
        return
      endif
      if (ichiyd(iyd).ne.0) then
        write(6,"('>>> inciyd: bad iyd=',i6)") iyd
        return
      endif
      leap = 0
      iyr = iyd/1000
      if (mod(iyr,4).eq.0) leap = 1
      iday = iyd-iyr*1000
      inc = 1
      if (ndays.lt.0) inc = -1
      do i=1,abs(ndays)
        iday = iday+inc
        if (leap.le.0.and.iday.gt.365) then
          iday = 1
          iyr = iyr+1
        endif
        if (leap.gt.0.and.iday.gt.366) then
          iday = 1
          iyr = iyr+1
        endif
        if (iday.lt.1) then
          iyr = iyr-1
          if (mod(iyr,4).eq.0) then
            iday = 366
          else
            iday = 365
          endif
        endif
      enddo
      inciyd = iyr*1000+iday
      return
      end
!-----------------------------------------------------------------------
      subroutine iyd2date(iyd,imo,ida,iyr)
      implicit none
      integer,intent(in) :: iyd
      integer,intent(out) :: imo,ida,iyr 
c
c Convert given iyd (yyddd) to date, returning month, day, and
c   2-digit year in imo,ida,iyr:
c
      integer :: ndmon(12) =
c                 J  F  M  A  M  J  J  A  S  O  N  D
     |          (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: iday,id,m
      integer,external :: ichiyd
c
      imo = 0
      ida = 0
      iyr = 0
      if (ichiyd(iyd).ne.0) then
        write(6,"('iyd2date: bad iyd=',i5)") iyd
        return
      endif
      iyr = iyd/1000
      ndmon(2) = 28
      if (mod(iyr,4).eq.0) ndmon(2) = 29
      iday = iyd-iyr*1000
      id = 0
      do m=1,12
        id = id+ndmon(m)
        if (id.eq.iday) then
          id = ndmon(m)
          goto 100
        endif
        if (id.gt.iday) then
          id = ndmon(m)-id+iday
          goto 100
        endif
      enddo
      write(6,"('>>> iyd2date: could not find date for iyd=',i5)") iyd
      imo = 0
      ida = 0
      iyr = 0
      return
 100  continue
      imo = m
      ida = id
      return
      end
!-----------------------------------------------------------------------
      integer function iyddif(iyd0,iyd1)
      implicit none
      integer,intent(in) :: iyd0,iyd1
      integer :: ndmon(12) =
c                 J  F  M  A  M  J  J  A  S  O  N  D
     |          (/31,28,31,30,31,30,31,31,30,31,30,31/)
      integer :: iyr0,iyr1,i,ida0,ida1
      integer,external :: ichiyd
c
c Return number of days from iyd0 to iyd1:
c
      if (iyd0.eq.iyd1) then
        iyddif = 0
        return
      endif
      if (iyd0.gt.iyd1) then
        write(6,"('iyddif: iyd0 > iyd1 (iyd0,1=',2i6,')')")
     +    iyd0,iyd1
        iyddif = -1
        return
      endif
      if (ichiyd(iyd0).ne.0) then
        write(6,"('>>> iyddif: bad iyd0=',i6)") iyd0
        iyddif = -1
        return
      endif
      if (ichiyd(iyd1).ne.0) then
        write(6,"('>>> iyddif: bad iyd1=',i6)") iyd1
        iyddif = -1
        return
      endif
      iyr0 = iyd0/1000
      iyr1 = iyd1/1000
      iyddif = 0
      if (iyr1.gt.iyr0) then
        do i=iyr0,iyr1-1
          if (mod(i,4).ne.0) then
            iyddif = iyddif+365
          else
            iyddif = iyddif+366
          endif
        enddo
      endif
      ida0 = iyd0-iyr0*1000
      ida1 = iyd1-iyr1*1000
      iyddif = iyddif+(ida1-ida0)
      return
      end
