C** time_manip.f 
C** 
C** Subroutines and functions for manipulating time. 
C**
C** Time formats: 
C      Gregorian Day : days since some start date.
C                      KD (integer) day counter; 
C                      SECKD (real) elapsed seconds in day KD
C
C      DMY : Day Month Year seconds 
C            IDD (integer)  current day (= 1 for first day of the month)
C            IMM (integer)  current month (1-12) 
C            IYear (integer) current year (e.g. 1999)
C            SECKD (real) elapsed seconds in day IDD. 
C            the seconds counter is identical in DMY and Gregorian Day.
C
C      UTC0: Elapsed decimal days in the specified year. 
C            IYear (integer) current year (e.g. 1999)
C            YDAYUTC0 - year day in utc0 format.
C            part of a day are represented as fractional days. 
C            e.g. Noon Jan 1 is day 0.5. 
C
C Format Conversion Routines
C    DMY to Gregorian Day :  call GDAY(IDD,IMM,IYear,KD)
C    Gregorian Day to DMY :  call  DMY(IDD,IMM,IYear,KD)
C    Gregorian Day to UTC0:  call GDAY_UTC0(KD,SECKD,IYEAR,YDAYUTC0)
C    UTC0 to Gregorian Day:  call UTC0_GDAY(IYEAR,YDAYUTC0,KD,SECKD)
C
C Gregorian Time manipulation 
C    Compute difference between two Gregorian times: 
C                    call gday_diff(kd1,sec1,kd2,sec2,kd_diff,sec_diff)
C    Compare two Gregorian times: igday_compare(kd1,sec1,kd2,sec2)
C    Update Gregorian time where seckd >=0 :     call up_date(kd,seckd)
C    Update Gregorian time for arbitrary seckd : call up_date2(kd,seckd)
C
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE GDAY(IDD,IMM,IYear,KD)
C!
C!  GIVEN DAY,MONTH,(EACH 2 DIGITS) and year (four digits), GDAY RETURNS
C!  THE DAY#, KD BASED ON THE GREGORIAN CALENDAR.
C!  THE GREGORIAN CALENDAR, CURRENTLY 'UNIVERSALLY' IN USE WAS
C!  INITIATED IN EUROPE IN THE SIXTEENTH CENTURY. NOTE THAT GDAY
C!  IS VALID ONLY FOR GREGORIAN CALENDAR DATES.
C
C   KD=1 CORRESPONDS TO JANUARY 1, 0000
C	
c 	Note that the Gregorian reform of the Julian calendar 
c	omitted 10 days in 1582 in order to restore the date
c	of the vernal equinox to March 21 (the day after
c	Oct 4, 1582 became Oct 15, 1582), and revised the leap 
c	year rule so that centurial years not divisible by 400
c	were not leap years.
c
C   THIS ROUTINE WAS WRITTEN BY EUGENE NEUFELD, AT IOS, IN JUNE 1990.
C
      INTEGER NDP(13)
      INTEGER NDM(12)
	
      DATA NDP/0,31,59,90,120,151,181,212,243,273,304,334,365/
      DATA NDM/31,28,31,30,31,30,31,31,30,31,30,31/
C!
	lp=6
c make iyy and icc variables
        icc=iyear/100
        iyy=iyear-icc*100

C!  TEST FOR INVALID INPUT:
      IF(ICC.LT.0)THEN
	 WRITE(LP,5000)ICC
	 STOP
      ENDIF
      IF(IYY.LT.0.OR.IYY.GT.99)THEN
	 WRITE(LP,5010)IYY
	 STOP
      ENDIF
      IF(IMM.LE.0.OR.IMM.GT.12)THEN
	 WRITE(LP,5020)IMM
	 STOP
      ENDIF
      IF(IDD.LE.0)THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.NE.2.AND.IDD.GT.NDM(IMM))THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.EQ.2.AND.IDD.GT.29)THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.EQ.2.AND.IDD.GT.28.AND.((IYY/4)*4-IYY.NE.0.OR.(IYY.EQ.0.AND
     .    .(ICC/4)*4-ICC.NE.0)))THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
5000  FORMAT(' INPUT ERROR. ICC = ',I7)
5010  FORMAT(' INPUT ERROR. IYY = ',I7)
5020  FORMAT(' INPUT ERROR. IMM = ',I7)
5030  FORMAT(' INPUT ERROR. IDD = ',I7)
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST CENTURY:
      KD = ICC*36524 + (ICC+3)/4
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST YEAR:
      KD = KD + IYY*365 + (IYY+3)/4
C!
C!  ADJUST FOR CENTURY RULE:
C!  (VIZ. NO LEAP-YEARS ON CENTURYS EXCEPT WHEN THE 2-DIGIT
C!  CENTURY IS DIVISIBLE BY 4.)
      IF(IYY.GT.0.AND.(ICC-(ICC/4)*4).NE.0) KD=KD-1
C!  KD NOW TRULY REPRESENTS THE DAY# OF THE LAST DAY OF LAST YEAR.
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST MONTH:
      KD = KD + NDP(IMM)
C!
C!  ADJUST FOR LEAP YEARS:
      IF(IMM.GT.2.AND.((IYY/4)*4-IYY).EQ.0.AND.((IYY.NE.0).OR.
     .   (((ICC/4)*4-ICC).EQ.0)))   KD=KD+1
C!  KD NOW TRULY REPRESENTS THE DAY# OF THE LAST DAY OF THE LAST
C!  MONTH.
C!
C!  CALCULATE THE CURRENT DAY#:
      KD = KD + IDD
      RETURN
C!
C!
      ENTRY DMY(IDD,IMM,IYear,KD)
C!
C!  GIVEN THE (GREGORIAN) DAY#, KD, AS CALCULATED ABOVE IN THIS ROUTINE,
C!  ENTRY DMY RETURNS THE (GREGORIAN) DAY, MONTH, YEAR AND CENTURY.
C!
C!  TEST FOR VALID INPUT:
      IF(KD.LE.0) WRITE(LP,5040)KD
5040  FORMAT(' KD = ',I7,'  INVALID INPUT. DMY STOP.')
C!
C!  SAVE KD
      KKD=KD
C!  CALCULATE ICC AND SUBTRACT THE NUMBER OF DAYS REPRESENTED BY ICC
C!  FROM KKD
C!  JFH IS THE NUMBER OF 400 YEAR INTERVALS UP TO KKD
C!  JCC IS THE NUMBER OF ADDITIONAL CENTURIES UP TO KKD
      JFH = KKD/146097
      KKD = KKD - JFH*146097
      IF(KKD.LT.36525)THEN
	 JCC = 0
      ELSE
	 KKD = KKD - 36525
	 JCC = 1 + KKD/36524
	 KKD = KKD - (JCC-1)*36524
      END IF
      ICC = 4*JFH + JCC
      IF(KKD.EQ.0)THEN
	 ICC = ICC-1
	 IYY = 99
	 IMM = 12
	 IDD = 31
c	 RETURN
	go to 110
      ENDIF
C!
C!  CALCULATE IYY. JFY IS THE NUMBER OF FOUR YEAR INTERVALS IN THE
C!  CURRENT CENTURY. THE FIRST FOUR YEAR INTERVAL IS SHORT (1460 DAYS
C!  RATHER THAN 1461)IF THE CURRENT CENTURY IS NOT DIVISIBLE BY 4, AND
C!  IN THIS CASE JCC.NE.0 AS CALCULATED ABOVE.
C!
C!  CALCULATE JFY:
      JFY = 0
      IF(JCC.EQ.0)GOTO 10
      IF(KKD.LT.1460)GOTO 10
      JFY = 1
      KKD = KKD - 1460
10    KK = KKD/1461
      JFY = JFY + KK
      KKD = KKD - KK*1461
C!
C!  CALCULATE JYY, THE REMAINING YEARS OF THE CURRENT CENTURY UP TO THE
C!  CURRENT DAY:
      JYY = 0
C!  THE NEXT YEAR IS NOT A LEAP YEAR IF JFY=0 AND JCC.NE.0.
      IF(JFY.EQ.0.AND.JCC.NE.0)GOTO 20
      IF(KKD.LT.366)GOTO 30
      JYY = 1
      KKD = KKD - 366
20    JYYY = KKD/365
      JYY = JYY + JYYY
      KKD = KKD - JYYY*365
30    IYY = 4*JFY + JYY
      IF(KKD.EQ.0) THEN
	 IYY=IYY-1
	 IMM=12
	 IDD=31
c	 RETURN
	go to 110
      END IF
C!
C!  SET L=1 IF WE HAVE A LEAP YEAR.
      L=0
      IF(IYY-(IYY/4)*4.NE.0)GOTO 40
      IF(IYY.EQ.0.AND.(ICC-(ICC/4)*4).NE.0)GOTO 40
      L=1
C!
C!  CALCULATE IMM AND IDD
40    IF(KKD.GT.31) GOTO 50
      IMM=1
      IDD=KKD
c      RETURN
       go to 110
C!
50    IF(KKD.GT.59)GOTO 60
      IMM = 2
      IDD = KKD-31
c      RETURN
	go to 110
C!
60    IF(KKD.GT.60)GOTO 70
      IF(L.EQ.0)GOTO 70
      IMM = 2
      IDD = 29
c      RETURN
	go to 110
C!
70    IF(L.EQ.1) KKD=KKD-1
      DO 80 I=4,13
	 IF(KKD.GT.NDP(I))GOTO 80
	 IMM = I-1
	 IDD = KKD - NDP(I-1)
c	 RETURN
	go to 110
C!
80    CONTINUE

 110 	iyear=icc*100+iyy
	return

90    WRITE(LP,5050)
5050  FORMAT(' ERROR IN DMY.')
      STOP
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE GDAY_UTC0(KD,SECKD,IYEAR,YDAYUTC0)
C-----------------------------------------------------------------------
      CALL DMY(ID,IM,IYEAR,KD)
      CALL GDAY(01,01,IYEAR,KDNEWYEAR)
      YDAYUTC0=KD-KDNEWYEAR+SECKD/3600.0/24.0
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE UTC0_GDAY(IYEAR,YDAYUTC0,KD,SECKD)
C-----------------------------------------------------------------------
      CALL GDAY(01,01,IYEAR,KDNEWYEAR)
      KD=KDNEWYEAR+INT(AINT(YDAYUTC0))
      SECKD=24.0*3600.0*(YDAYUTC0-AINT(YDAYUTC0))
      RETURN
      END
C
C***********************************************************************
C***********************************************************************
C gday_diff.f
C
C  compute difference of two times: kd1,sec1 and kd2,sec2
C  Time1 minus Time2
C  return values:  difference expressed in days and seconds. 
C
C  Convention: seconds are always positive. 
C              if kd1 = kd2 and sec1-sec2 = -1 
c              then kd_d = 0 and sec_d = 86399.0 
C
C
C  what do I do if kd2 is less than kd1?
C  what is the sign if the seconds if kd1 - kd2 < 0 
C
C  I am content to have sec_d < 0 if kd_d = 0 
C
C  Careful testing is required if kd_d <0 and up_date2 is called.
C 
      subroutine gday_diff(kd1,sec1,kd2,sec2,kd_d,sec_d)
	 integer kd1, kd2,kd_d
	 real    sec1, sec2, sec_d

	 integer kd1_l , kd2_l
	 real    sec1_l, sec2_l 

C create local copies
	 kd1_l = kd1
	 sec1_l = sec1
	 kd2_l = kd2
	 sec2_l = sec2

C insist that sec1_l and sec2_l fall in the range 0, ..., 86400 seconds
	 call up_date2(kd1_l, sec1_l)
	 call up_date2(kd2_l, sec2_l)

c write(*,*)kd1_l, sec1_l,kd2_l, sec2_l

	 kd_d = kd1_l - kd2_l
	 sec_d = sec1_l - sec2_l
	 call up_date2(kd_d,sec_d)

      return 
      end
C
C***********************************************************************
C***********************************************************************
C igday_compare.f
C
C  compare two times: kd1,sec1 and kd2,sec2
C  return values:  1 kd1,sec1 > kd2,sec2
C                  0  kd1,sec1 = kd2,sec2
C                  -1  kd1,sec1 < kd2,sec
C 
      integer function igday_compare(kd1,sec1,kd2,sec2)
	 integer kd1, kd2
	 real    sec1, sec2

	 integer kd1_l , kd2_l
	 real    sec1_l, sec2_l 

C create local copies
	 kd1_l = kd1
	 sec1_l = sec1
	 kd2_l = kd2
	 sec2_l = sec2

C insist that sec1_l and sec2_l fall in the range 0, ..., 86400 seconds
	 call up_date2(kd1_l, sec1_l)
	 call up_date2(kd2_l, sec2_l)

c write(*,*)kd1_l, sec1_l,kd2_l, sec2_l

	 if(kd1_l.lt.kd2_l) then 
	     igday_compare = -1
	 elseif(kd1_l.gt.kd2_l) then 
	     igday_compare = 1
	 else 
	    if(sec1_l.lt.sec2_l) then 
	       igday_compare = -1
	    elseif(sec1_l.gt.sec2_l) then 
	       igday_compare = 1
	    else
	       igday_compare = 0
	    endif
	 endif
      return 
      end
C
C***********************************************************************
C***********************************************************************
	Subroutine up_date(kd,ssec)
C 
C this subroutine updates the current gregorian day number kd
C and resets the time, ssec, accordingly
C it inherently assumes that the variable is updated every day!

      integer kd
      real ssec
	
      kd=kd+aint(ssec/8.64E4)
      ssec=amod(ssec,8.64E4)
 
      return
      end
C
C***********************************************************************
C***********************************************************************
	subroutine up_date2(kd,ssec)
C 
C Updates the current gregorian day number kd
C and resets ssec (the second counter), accordingly. 
C This routine deals with the case where time increases *and*
C the case where time decreases. 
C Oct 29, 1998 C. Hannah 
C 
C Based on Up_date from QUODDY4.1_1 
C
C Nov 18, 1998  : CGH
C  checked ibm,sgi, sun, dec-alpha in each case kd and ssec
C    were updated correctly, i.e. after call to up_date2
C    0 <= ssec < 86400 
C However during the development of this subroutine I left a note
C that:
C   kd does not change until the clock tic after midnight 
C           i.e.  ssec > 86400. second  
C Need to watch out. 
C   
C
C Note that in single precision round off starts to be a problem
C at the 0.01 second level.  Not clear how big the problem can be.
C Simple solution is to make the second counter (ssec) double precision. 
C 
C Comment from the original (that I do not understand)
C it inherently assumes that the variable is updated every day!

      integer kd
      real ssec

      real day2sec
      parameter (day2sec = 86400.0)
	
      if(ssec.gt.0.0) then 
         kd = kd + int(ssec/day2sec) 
         ssec = mod(ssec,day2sec)
      elseif(ssec.lt.0.0) then 
         kd = kd + int(ssec/day2sec) - 1 
         ssec = mod(ssec,day2sec) + day2sec
      else
      endif
 
      return
      end
C
C***********************************************************************
C***********************************************************************
