
!=========================================================================================================================================
SUBROUTINE UNPACK_UT(packed_time, time_ut)
!=========================================================================================================================================
!WARNING:  Character variable 'time_ut' must be at least 19 bytes long in calling program unit
!WARNING:  VALID ONLY FOR START DATE = 1970 (current for Windows OS as at Apr 2013)
!WARNING:  VALID ONLY FOR YEARS TO ~2038  ('packed_time' requires more than 4 bytes for dates beyond ~2038, unless start date changed)
!NB: Coding accounts for 'centruy/400' leap year rule for year 2100 only.  But not invoked since currently limited to years before ~2038. 
!D Spurr, 8 April 2013
!=========================================================================================================================================
IMPLICIT NONE
    integer, intent(IN) :: packed_time  
    character(*), intent(OUT) :: time_ut
    integer :: century, year, days  
    integer :: quadrennials_since_1970                !number of 4 year periods since start of 1970
    integer :: quadrennials_since_2102                !number of 4 year periods since start of 2102
    integer :: residual 
    integer :: normal_year_secs =  31536000           !seconds per normal year (365 days) 
    integer :: leap_year_secs   =  31622400           !seconds per leap year (366 days)
    integer :: quadrennial_secs = 126230400           !seconds per 4 years (365 *3 + 366 = 1461 days)
    integer :: nonleap_quadrennial_secs = 126144000   !seconds per 4 non leap years (365 *4 = 1460 days)
    integer :: mth_start_normal(12) = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]   !Days in normal year before start of month
    integer :: mth_start_leap(12)   = [0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335]   !Days in leap year before start of month
!=========================================================================================================================================

    quadrennials_since_1970 = INT (packed_time / quadrennial_secs)
    if(quadrennials_since_1970 < 32) then
        residual = packed_time - quadrennials_since_1970 * quadrennial_secs
    elseif(quadrennials_since_1970 == 32) then  !NB: ### non-leap 'century years' from 2200 on not modelled ### 
        residual = packed_time - 32 * quadrennial_secs
        if(residual > nonleap_quadrennial_secs) then
            quadrennials_since_1970 = 33
            residual = residual - nonleap_quadrennial_secs
        endif
    else !quadrennials_since_1970 > 33
        residual = packed_time - 32 * quadrennial_secs - nonleap_quadrennial_secs
        quadrennials_since_2102 = INT (residual / quadrennial_secs)
        quadrennials_since_1970 = quadrennials_since_2102 + 33
        residual = packed_time - quadrennials_since_2102 * quadrennial_secs
    endif
    
    year = 1970 + quadrennials_since_1970 * 4
    
    if(residual < normal_year_secs) then !yr1 in 4 yr cycle
        call DATETIME_FOR_YEAR (year, residual, mth_start_normal, time_ut)
    else
        year = year + 1
        residual = residual - normal_year_secs
        if(residual < normal_year_secs) then  !yr2 in 4 yr cycle 
            call  DATETIME_FOR_YEAR (year, residual, mth_start_normal, time_ut)
        else
            year = year + 1
            residual = residual - normal_year_secs
            if(residual < leap_year_secs) then  !yr3 in 4 yr cycle (a leap year)
                century = INT(year/100) * 100
                if(year == century) then
                    if(MOD(century,400) > 0) then
                        call  DATETIME_FOR_YEAR (year, residual, mth_start_normal, time_ut)
                    else
                        call  DATETIME_FOR_YEAR (year, residual, mth_start_leap, time_ut)
                    endif
                else
                    call  DATETIME_FOR_YEAR (year, residual, mth_start_leap, time_ut)
                endif
            else
                year = year + 1
                residual = residual - leap_year_secs
                if(residual < normal_year_secs) then  !yr4 in 4 yr cycle
                    call  DATETIME_FOR_YEAR (year, residual, mth_start_normal, time_ut)
                else
                    !ERROR - should never reach here (code error if do)
                endif !yr4
            endif !yr3
        endif !yr2
    endif !yr1
        
RETURN

END SUBROUTINE UNPACK_UT


    
!==============================================================================
! CALCULATE UT FOR 'year' AND 'secs_since_prev_year'
! NB: Time '00:00:00' is treated as within the following day
!==============================================================================
SUBROUTINE DATETIME_FOR_YEAR (year, secs_since_prev_year, mth_start, time_ut)
IMPLICIT NONE
    integer, intent(IN) :: year  
    integer, intent(IN) :: secs_since_prev_year  
    integer, intent(IN) :: mth_start(12)  !Days in year before start of month
    character(*), intent(OUT) :: time_ut
    integer :: residual  
    integer :: days, mth, day, hr, mins, secs
    integer :: day_secs = 86400   !seconds per day
    
    days = INT (secs_since_prev_year/day_secs)  !days in year before current date
    residual = secs_since_prev_year - days * day_secs

    mth = SUM(MIN(MAX(days + 1 - mth_start,0),1))
    day  = days + 1 - mth_start(mth)
    
    hr  =  INT (residual/3600)
    residual = residual - hr*3600
    mins = INT (residual/60)
    residual = residual - mins*60
    secs = INT (residual)
    write(time_ut,'(i4,5(a,i2.2))') year,'_', mth,'_', day, ' ', hr,':', mins,':', secs
    
END SUBROUTINE DATETIME_FOR_YEAR
