subroutine day_code(month,day,year,code)
c
c=======================================================================
c                                                                    ===
c  This subroutine computes a code for the day of the week, given    ===
c  the date. This code is good for date after:                       ===
c                                                                    ===
c                              January 1, 1752 AD                    ===
c                                                                    ===
c  the year the Gregorian calander was adopted in Britian and the    ===
c  American colonies.                                                ===
c                                                                    ===
c  On Input:                                                         ===
c                                                                    ===
c     DAY       The day of the month (integer).                      ===
c     MONTH     The month, 1=January, 2=February, ... (integer).     ===
c     YEAR      The year, including the century (integer).           ===
c                                                                    ===
c  On Output:                                                        ===
c                                                                    ===
c     CODE      A code for the corresponding day of the week         ===
c               (integer):                                           ===
c               CODE = 0  =>  Sunday                                 ===
c               CODE = 1  =>  Monday                                 ===
c               CODE = 2  =>  Tuesday                                ===
c               CODE = 3  =>  Wednesday                              ===
c               CODE = 4  =>  Thursday                               ===
c               CODE = 5  =>  Friday                                 ===
c               CODE = 6  =>  Saturday                               ===
c                                                                    ===
c     Calls:  none                                                   ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      logical leap_flag
      integer base_cen,base_qcen,base_qyear,base_year,bym1_dec31,code,
     *        day,feb_end,leap,month,n,no_day,no_yr,nqy,nyc,nyqc,year
      integer month_day(12)
      parameter (base_cen=1700,base_qcen=1600,base_qyear=1748,
     *           base_year=1752,bym1_dec31=5,feb_end=59) 
      data month_day /31,28,31,30,31,30,31,31,30,31,30,31/
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Compute the number of years since the base year, the number of
c  years since the beginning of the base century and the number of
c  years since the beginning of the base 400 year.
c-----------------------------------------------------------------------
c
      no_yr=year-base_year
      nqy=year-base_qyear
      nyc=year-base_cen
      nyqc=year-base_qcen
c
c-----------------------------------------------------------------------
c  Compute the number of leapdays in that time.  Determine if this
C  is a leap year.
c-----------------------------------------------------------------------
c
      leap=nqy/4-nyc/100+nyqc/400
      leap_flag=((mod(nqy,4).eq.0).and.(mod(nyc,100).ne.0)).or.
     *           (mod(nyqc,400).eq.0)
c
c-----------------------------------------------------------------------
c  Compute the number of days this year.  The leap year corrections
c  are:
c        Jan. 1 - Feb. 28   Have not had the leap day counted above.
c        Feb.29             Counting leap day twice.
c-----------------------------------------------------------------------
c
      no_day=day
      do 10 n=1,month-1
        no_day=no_day+month_day(n)
  10  continue
      if (leap_flag.and.(no_day.le.feb_end))  no_day=no_day-1
      if (leap_flag.and.(month.eq.2).and.(day.eq.29)) no_day=no_day-1
c
c-----------------------------------------------------------------------
c  Compute the total number of days since Jan. 1 of the base year,
c  exclusive of the 364 day per year which represent an even 52
c  weeks.  Actually, only need to do the addition mod 7.
c-----------------------------------------------------------------------
c
      no_day=mod(no_day,7)+mod(leap,7)+mod(no_yr,7)+bym1_dec31
c
c-----------------------------------------------------------------------
c  Get the day of the week code.
c-----------------------------------------------------------------------
c
      code=mod(no_day,7)
      return
      end