subroutine gregorian (julian,iday,month,iyear,iyday)
c
c=======================================================================
c                                                                    ===
c  This routine converts Julian day number to calendar (Gregorian)   ===
c  date.                                                             ===
c                                                                    ===
c  On Input:                                                         ===
c                                                                    ===
c     JULIAN   Julian day (integer)                                  ===
c                                                                    ===
c  On Ouput:                                                         ===
c                                                                    ===
c     IDAY     day of the month (integer)                            ===
c     IYDAY    year day (integer)                                    ===
c     IYEAR    year (integer)                                        ===
c     MONTH    month of the year (integer)                           ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
       integer iday,igreg,ileap,iyd(13),iyday,iydl(13),iyear,ja,jalpha,
     &         jb,jc,jd,je,julian,month
       parameter (igreg=2299161)
       data iyd  /1,32,60,91,121,152,182,213,244,274,305,335,366/
       data iydl /1,32,61,92,122,153,183,214,245,275,306,336,367/
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
       if (julian.ge.igreg) then
         jalpha=int((float(julian-1867216)-0.25)/36524.25)
         ja=julian+1+jalpha-int(0.25*float(jalpha))
       else
         ja=julian
       endif
       jb=ja+1524
       jc=int(6680.+(float(jb-2439870)-122.1)/365.25)
       jd=365*jc+int(0.25*float(jc))
       je=int(float(jb-jd)/30.6001)
       iday=jb-jd-int(30.6001*float(je))
       month=je-1
       if (month.gt.12) month=month-12
       iyear=jc-4715
       if (month.gt.2) iyear=iyear-1
       if (iyear.le.0)iyear=iyear-1
       ileap=mod(iyear,4)
       if (ileap.eq.0) then
           iyday=iydl(month)+iday-1
       else
           iyday=iyd(month)+iday-1
       endif
       return
       end