/*
*********************************************************************
* Convert Universal Time (in standard clock time) to calendar data;
* apply Daylight Saving Time to provide adjusted CLOCK time, when required.

*   Input
*     seconds       = universal time [s]
*     DST_hour_1    = hour in year, when daylight saving starts
*     DST_hour_2    = ditto ends
*   Output
*     year_no       = the number of the current year relative to 1900
*     day1_in_year  = day number within current year (0-364),
*                     possible leap day removed to prepare
*                     for comparison with rules.
*     hour1_in_year = hours from beginning of year, clock time
*     hour1_in_day  = clock hour in current day,
*                     corrected for daylight saving
*     month         = month in current year
*     day           = day in current month
*     weekday       = day in week, Monday = 1

  _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
  Revisions:

  PG 020301  Added func DST_Hour used for DST-correction

  LE 020220  Function StepF moved from schedule.nmf

*/

FUNCTION 

  VOID UT2Calendar (seconds, DST_hour_1, DST_hour_2,
                    year_no, day1_in_year, hour1_in_year,
                    hour1_in_day, month, day, week_day)
 
  LANGUAGE F77

  INPUT 
     FLOAT seconds, DST_hour_1, DST_hour_2;

  OUTPUT
     FLOAT year_no, day1_in_year, hour1_in_year, 
           hour1_in_day, month, day, week_day;

  CODE
      SUBROUTINE UT2Calendar (seconds, DST_hour_1, DST_hour_2,
     &   year_no, day1_in_year, hour1_in_year,
     &   hour1_in_day, month, day, week_day)

      REAL*8 seconds, DST_hour_1, DST_hour_2,
     &   year_no, day1_in_year, hour1_in_year,
     &   hour1_in_day, month, day, week_day

      INTEGER
     &   leapPeriods, yearInLeapPeriod, leapYear,
     &   dayFrom1901, dayInPeriod, dayFromJan1, dayFromMarch1,
     &   mm, dd

      REAL*8 hour_from_1900, hour_in_day, hour_in_year

      REAL*8 stepF

C45
* Wrong input (date outside 1901 -- 2099) is signalled by
* year_no < 0 (and all other output variables)

* Float quantities have underscore, integers are written together
* Corrected quantities are marked with a 1

* Quantities with possible leap day removed, used for calculation
* of day light saving
*   hour_in_year    hour in year, corrected to a 365 day year,
*                   standard time
*   hour1_in_year   ditto, clock time w possible DST
*   day1            day number, corrected to a 365 day year


* Get year no and day no (0-365) within year
      hour_from_1900 = seconds/3600
C45=  dayFrom1901 = Int (hour_from_1900/24) - 366
      dayFrom1901 = Int (hour_from_1900/24) - 365
C45   IF (dayFrom1901 .LT. 0) dayFrom1901 = 0

      leapPeriods = Int(dayFrom1901/1461)
      dayInPeriod = Mod(dayFrom1901, 1461)

      yearInLeapPeriod = Int(dayInPeriod/365) - Int(dayInPeriod/1460)
      dayFromJan1 = dayInPeriod - yearInLeapPeriod*365
      leapYear = Int(yearInLeapPeriod/3)

* Export year no (put back 1900)
      year_no = leapPeriods*4 + yearInLeapPeriod + 1
      IF (dayFrom1901 .LT. 0
     & .OR. year_no .LT. 1
     & .OR. year_no .GT. 199) THEN
        year_no = -1
        day1_in_year = -1
        hour1_in_year = -1
        hour1_in_day = -1
        month = -1
        day = -1
        week_day = -1
        RETURN
      ENDIF

* Hour in day and in year, std time
      hour_in_day = Mod(hour_from_1900, 24.d0)
      hour_in_year = dayFromJan1 * 24 + hour_in_day

* clock time = f(std time)
*
* CT takes an 1 hour jump when DST goes on. When DST has to go off,
* we let it grow at half pace for two hours, rather then setting it
* back 1 hour.
*
*  clock  ^
*  time   |                   /
*         |                 /
*         |               /
*         |            .
*         |         /
*         |       /
*         |     /
*         |    |
*         |   /
*         | /
*         | ---|-----|---|-------->
*            on        off        std time [hr]
*                 off-2

      hour1_in_year = hour_in_year +
     &   stepF (hour_in_year - DST_hour_1) -
     &   0.5*stepF (hour_in_year - (DST_hour_2-2)) *
     &      (hour_in_year-(DST_hour_2-2)) +
     &   0.5*stepF (hour_in_year - DST_hour_2) *
     &      (hour_in_year-DST_hour_2)

* Hour in day and day in year, clock time
      hour1_in_day = Mod(hour1_in_year, 24.d0)
      day1_in_year = Int(hour1_in_year/24)

* Calendar data
* Leap year = 31+30+31+30+31 + 31+30+31+30+31 + 31+28(29)
*             m  a  m  j  j    a  s  o  n  d    j  f
      dayFromMarch1 = Mod(day1_in_year + 365-59, 365.d0+leapYear)
      mm = Int(dayFromMarch1/153) * 5
      dd = Mod(dayFromMarch1, 153)

      mm = mm + Int(dd/61) * 2
      dd = Mod(Int(dd), 61)

*             j  f  m  a  m  j    j  a  s  o  n  d
      mm = 2 + mm + Int(dd/31)
      dd = Mod(dd, 31)

* Deliver in ranges (1..12), (1..31), (1..7)
*     month = 1 + mm - Int(mm/12) * 12
      month = 1 + Mod(mm, 12)
      day = 1 + dd
      week_day = 1
     &  + Mod (Int ( (hour_from_1900
     &               +  stepF (hour_in_year - DST_hour_1)
     &               -  stepF (hour_in_year - DST_hour_2)
     &               ) / 24
     &             ) + 0
     &        , 7
     &        )


ctrc51      write(7,97) leapPeriods, yearInLeapPeriod, leapYear,
ctrc51     &   dayFrom1901, dayInPeriod, dayFromJan1, dayFromMarch1,
ctrc51     &   hour_from_1900, hour_in_year, hour1_in_year, hour1_in_day,
ctrc51     &   dst_hour_1, dst_hour_2
ctrc51   97 format (/, '#     lp     ylp      ly     d01',
ctrc51     &   '     dip,     dj1     dm1',
ctrc51     &   /, 1x, 7i8,
ctrc51     &   /, 1x, 'hr_fr_1900    hr_in_year    hr1_in_year   ',
ctrc51     &          'hr1_in_day    dst_hr_1      dst_hr_2',
ctrc51     &   /, 1x, 6g14.6)

      RETURN
      END

  END_CODE

END_FUNCTION

/*
*********************************************************************
  Step function = 1 IF X >= 0
*/

FUNCTION

  FLOAT StepF (X)

  LANGUAGE F77

  INPUT
     FLOAT X;

CODE

      FUNCTION StepF (X)

      DOUBLE PRECISION StepF, X

      IF (X .GE. 0.) THEN
         StepF = 1.
      ELSE
         StepF = 0.
      ENDIF

      RETURN
      END

END_CODE

END_FUNCTION


/*
*********************************************************************
  DST_Correction
*/

FUNCTION

  FLOAT DST_Hour (hour, DST_start, DST_end)

  LANGUAGE F77

  INPUT
     FLOAT hour, DST_start, DST_end;

  CODE


      FUNCTION DST_Hour(hour, DST_start, DST_end)
      DOUBLE PRECISION hour, DST_start, DST_end, DST_hour, StepF

      DST_Hour = hour -
     &   StepF (hour - DST_start) +
     &   StepF (hour - (DST_end-1) ) *
     &      (hour - (DST_end-1) ) -
     &   StepF (hour - DST_end) *
     &      (hour - DST_end)

      RETURN
      END

END_CODE

END_FUNCTION