c
c (0) Initialize program (write in the code or read from a file, or from command line arguments.)
c
      character*10 weekday
c
      maxtries = 100 ! number of maximum years to be searched into the future.
      ntries = 0
      idirtime = 1 ! =1 to test future years, = -1 two test past ones.
c Find a day that accords to your query, in this case
c the 12 of June of year yyyy and save it in variables iyyy,  month, iday.
      iyyy = 2021
      month =  6
      iday =  26
      weekday = 'Saturday'
C
      print*,' '
      print*,' Finding matches for: ',weekday,iday,month,iyyy
      print*,' '
c(2) Convert this to Modified Julian Day (MJD), that is the number of days
c since some past epoch (actually the 0 hours of May 24, 1968 , 
c but it does not matter here what date it is.) Call this mjd0 and calculate it
c iwith the subroutine at the end of this comment and save as "mjd0" the returned
c integer argument "jday".
c
      call modjulianday(iday, iyyy, month, mjd0)
c
c(3) Convert increase the number of tries, counting this one,
c and increase  the year value:
    3 continue
      ntries = ntries+1

c     print*,' ntries ',ntries

      iyyy = iyyy+1*idirtime !  then convert to MJD with
c the subroutine below and save returned argument  "jday" as "mjd".
c
      call modjulianday(iday, iyyy, month, mjd)
c 
      mjdiff = mjd-mjd0
      dmjdiff = mjdiff
      dweeks = dmjdiff/7
      iweeks = dweeks
      resid = abs(dweeks-iweeks) ! Take the absolute value,
c so this works for searching the future and the past.
      isuccess = 0
      if(resid.lt.0.000001) isuccess = 1
      if(isuccess.eq.1) then
      print*,' There is a match on ',
     .weekday,iday,month,iyyy
      print*,' '
      endif
c or go to (3):
      if(ntries.lt.maxtries) go to 3
      stop
      end
C
      subroutine modjulianday(iday, iyear, month, mjd)
c
c     For calculating the julian day from day, month, year. 
c
      data offset /2440000/
c
      JDAY = iDAY-32075+1461*(iYEAR+4800+(MONTH-14)/12)/4+367*
     .(MONTH-2-(MONTH-14)/12*12)/12-3*((iYEAR+4900+(MONTH-14)/12)/100)/4
c
      mjd = jday-offset
c
      return
      end
c