subroutine mjdn_to_ymd(mjdn,iy,imo,id)
c
c given modified julian day number (mjdn) calculates year, month, and day.
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,imo,id,id2,m2,m3,m4,m5,m6,ndiv,ly,
     +        molens(12),mosum,i
      data molens/31,28,31,30,31,30,31,31,30,31,30,31/
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
      ly=0
      if(mod(iy,4).eq.0 .and.
     +     (mod(iy,100).ne.0 .or. mod(iy,400).eq.0) )ly=1
      molens(2)=28+ly
      imo=1
      mosum=0
      id2=id
      do i=1,11
        mosum=mosum+molens(i)
        if(id.gt.mosum)then
          id2=id-mosum
          imo=i+1
        endif
      enddo
      id=id2
      return
      end