program g2y
c
c  This is a simple driver program
c
      implicit none
      integer*4 gps_week,year,month,mday,hour,minute
      real*8    sec_of_week, second
      
      gps_week = 842
      sec_of_week = 475200.0d0
      
      call gps_to_ymdhms(gps_week, sec_of_week,
     $ year,month,mday,hour,minute,second)
     
      write(6,100)gps_week, sec_of_week,
     $ year,month,mday,hour,minute,second
100   format('GPS, YMDHMS:',i9,2x,f18.10,3x,5i5,f18.12)

      gps_week = 0
      sec_of_week = 0.0d0
      call ymdhms_to_gps(year,month,mday,hour,minute,second,
     $ gps_week, sec_of_week)
      write(6,100)gps_week, sec_of_week,
     $ year,month,mday,hour,minute,second
     
      stop
      end
      
c-------------------------------------------------------------------------     
      subroutine gps_to_ymdhms(gps_week, sec_of_week,
     $ year,month,mday,hour,minute,second)
      implicit none
      integer*4 gps_week,year,month,mday,hour,minute
      real*8    sec_of_week, second
      integer*4 jan61980, jan11901
      real*8    sec_per_day
      integer*4 guess, more, yday, mjd, days_fr_jan11901
      integer*4 leap_month_day, delta_yrs, num_four_yrs
      integer*4 regu_month_day, days_left, years_so_far
      real*8    fmjd
      
      dimension regu_month_day(13)
      dimension leap_month_day(13)
      data regu_month_day / 0, 31, 59, 90, 120, 151, 181, 212,
     $ 243, 273, 304, 334, 365 /
      data leap_month_day / 0, 31, 60, 91, 121, 152, 182, 213,
     $ 244, 274, 305, 335, 366 /
      
      jan61980 = 44244
      jan11901 = 15385
      sec_per_day = 86400.0d0
      
      mjd = gps_week*7 + sec_of_week/sec_per_day + jan61980
      fmjd = dmod(sec_of_week, sec_per_day)/sec_per_day
c      write(6,100) mjd, fmjd
c100   format('mjd, fmjd:', i9,2x,f16.14)
      
      days_fr_jan11901 = mjd - jan11901
      num_four_yrs = days_fr_jan11901/1461
      years_so_far = 1901 + 4*num_four_yrs
      days_left = days_fr_jan11901 - 1461*num_four_yrs
      delta_yrs = days_left/365 - days_left/1460
      
      year = years_so_far + delta_yrs
      yday = days_left - 365*delta_yrs + 1
      hour = fmjd*24.0d0
      minute = fmjd*1440.0d0 - hour*60.0d0
      second = fmjd*86400.0d0 - hour*3600.0d0 - minute*60.0d0
      
c      write(6,'('' year is:''i9)') year
c      write(6,'('' day-of-year is:''i9)') yday
      
      guess = yday*0.0320d0 + 1
      more = 0
      if( mod(year,4) .eq. 0 ) then
        if( (yday - leap_month_day(guess+1)) .gt. 0 )more = 1
        month = guess + more
        mday = yday - leap_month_day(guess+more)
      else
        if( (yday - regu_month_day(guess+1)) .gt. 0 )more = 1
        month = guess + more
        mday = yday - regu_month_day(guess+more)
      endif
      
c      write(6,'('' month is:''i9)') month
c      write(6,'('' day-of-month is:''i9)') mday
      
      return
      end
c-----------------------------------------------------------------------     
      subroutine ymdhms_to_gps( year,month,mday,hour,minute,second,
     $                          gps_week, sec_of_week)
      implicit none
      integer*4 gps_week,year,month,mday,hour,minute
      real*8    sec_of_week, second
      integer*4 jan61980, jan11901
      real*8    sec_per_day
      integer*4 yday, mjd, leap_month_day, regu_month_day
      real*8    fmjd
      
      dimension regu_month_day(12)
      dimension leap_month_day(12)
      data regu_month_day / 0, 31, 59, 90, 120, 151, 181, 212,
     $ 243, 273, 304, 334 /
      data leap_month_day / 0, 31, 60, 91, 121, 152, 182, 213,
     $ 244, 274, 305, 335 /
      
      jan61980 = 44244
      jan11901 = 15385
      sec_per_day = 86400.0d0
      
      if( mod(year,4) .eq. 0 ) then
        yday = leap_month_day(month) + mday
      else
        yday = regu_month_day(month) + mday
      endif
	
      mjd = ((year - 1901)/4)*1461 +  mod((year - 1901),4)*365
     $ + yday - 1 + jan11901
      fmjd = ((second/60.0d0 + minute)/60.0d0 + hour)/24.0d0
c      write(6,100) mjd, fmjd
c100   format(' ymdhms_to_gps  mjd,fmjd:', i9,2x,f16.14)
      
      gps_week = (mjd - JAN61980)/7
      sec_of_week = ( (mjd - JAN61980) - gps_week*7 + fmjd )*sec_per_day
      
      return
      end