subroutine mjdn_to_yd(mjdn,iy,id)
c
c given modified julian day number (mjdn) calculates year and day of
c year, on the gregorian (modern) calendar.  Works for all dates AD.
c could be made to work for dates from and including any year divisible by 400
c (eg -400 = 401 BC) by subtracting the mjdn of the last day of the
c previous year, instead of the last day of 1 BC as used here, and then
c subtracting the corresponding number of years from iy at the end.
c
      implicit none
      integer mjdn,iy,id,m2,m3,m4,m5,m6,ndiv
c
      iy=1
      m2=mjdn+678575
c
c -678575 is the mjdn of the last day of 1 BC (ie day 365 year 0)
c
      ndiv=m2/146097
c
c there are 146097 days in 400 years (4*36524 + 1)
c
      iy=iy+400*ndiv
      m3=m2-ndiv*146097
      ndiv=m3/36524
c
c there are 36524 days in 100 years if not including a year divisible by 400
c (76*365+24*366)
c
      iy=iy+100*ndiv
      m4=m3-ndiv*36524
      ndiv=m4/1461
c
c there are 1461 days in 4 years including a leap year (3*365 + 366)
c leap years are years divisible by 4 but not 100 unless also by 400.
c
      iy=iy+4*ndiv
      m5=m4-ndiv*1461
      ndiv=m5/365
      iy=iy+ndiv
      m6=m5-ndiv*365
      id=m6
c      id=m6+1
      if(m5.eq.0 .and. (m4.ne.0 .or. m3.eq.0))then
        iy=iy-1
        id=366
      endif
      if(id.eq.0)then
        iy=iy-1
        id=365
      endif
      return
      end