module n3date
  implicit none

  private
  public:: nwp_gettime, nwp_systime
  public:: nwp_ymdhm2seq, nwp_ymdh2seq, nwp_ymd2seq
  public:: nwp_seq2ymdhm, nwp_seq2ymdh, nwp_seq2ymd
  public:: isotime

  integer, parameter, private:: DAY = 60 * 24
 
contains

  character(16) function isotime(seq) result(result)
    use n3fmt, only: fmt_decimal
    integer, intent(in):: seq
    integer:: y, m, d, h, mn, stat
    character(16):: buf
    call nwp_seq2ymdhm(y, m, d, h, mn, seq)
    buf = '####-##-##T##:##'
    call fmt_decimal(y, buf(1:4))
    call fmt_decimal(m, buf(6:7))
    call fmt_decimal(d, buf(9:10))
    call fmt_decimal(h, buf(12:13))
    call fmt_decimal(mn, buf(15:16))
    result = buf
  end function

  subroutine nwp_gettime(tcname, y, m, d, h, mn, iseq)
    use n3io1, only: iounit_new
    character(*), intent(in):: tcname
    integer, intent(out):: y, m, d, h, mn, iseq
    integer:: iostat, iunit
    call iounit_new(iunit)
    if (iunit < 0) goto 900
    open(unit=iunit, file=tcname, iostat=iostat)
    if (iostat /= 0) goto 900
    read(iunit, '(I4,I3,I3,I3,I3)', iostat=iostat) y, m, d, h, mn
    if (iostat /= 0) goto 900
    close(iunit, iostat=iostat)
    call nwp_ymdhm2seq(y, m, d, h, mn, iseq)
    return
    900 continue
    call nwp_systime(y, m, d, h, mn, iseq)
  end subroutine

  subroutine nwp_systime(y, m, d, h, mn, iseq)
    integer, intent(out):: y, m, d, h, mn, iseq
    integer:: b(8)
    call date_and_time(values=b)
    call nwp_ymdhm2seq(b(1), b(2), b(3), b(5), b(6), iseq)
    iseq = iseq + b(4)
    call nwp_seq2ymdhm(y, m, d, h, mn, iseq)
  end subroutine

  subroutine nwp_ymdhm2seq(y, m, d, h, mn, iseq)
    integer, intent(in):: y, m, d, h, mn
    integer, intent(out):: iseq
    call nwp_ym2seq(y, m, iseq)
    iseq = iseq + mn + 60 * (h + 24 * d)
  end subroutine

  subroutine nwp_ymdh2seq(y, m, d, h, iseq)
    integer, intent(in):: y, m, d, h
    integer, intent(out):: iseq
    call nwp_ym2seq(y, m, iseq)
    iseq = iseq + 60 * (h + 24 * d)
  end subroutine

  subroutine nwp_ymd2seq(y, m, d, iseq)
    integer, intent(in):: y, m, d
    integer, intent(out):: iseq
    call nwp_ym2seq(y, m, iseq)
    iseq = iseq + DAY * d
  end subroutine

  ! CURRENTLY PRIVATE
  ! WORKS YEARS FROM -2282 TO 5883.
  subroutine nwp_ym2seq(y, m, iseq)
    integer, intent(in):: y, m
    integer, intent(out):: iseq
    integer:: yr, mo, dseq
    integer, parameter:: DAYS_PER_400YEARS = 36525 * 4 - 3
    if (m < 3) then
      mo = m + 12 + 1
      yr = y - 1
    else
      yr = y
      mo = m + 1
    endif
    dseq = (306 * mo) / 10 + 72620
    do, while (yr < 2000)
      yr = yr + 400
      dseq = dseq - DAYS_PER_400YEARS
    enddo
    do, while (yr >= 2400)
      yr = yr - 400
      dseq = dseq + DAYS_PER_400YEARS
    enddo
    dseq = dseq + 36525 * (yr - 2000) / 100 - (yr - 2000) / 100
    iseq = dseq * DAY
  end subroutine

  subroutine nwp_seq2ymdhm(y, m, d, h, mn, iseq)
    integer, intent(out):: y, m, d, h, mn
    integer, intent(in):: iseq
    integer:: residual
    call nwp_seq2ym(y, m, iseq, residual)
    d = residual / DAY
    residual = mod(residual, DAY)
    h = residual / 60
    residual = mod(residual, 60)
    mn = residual
  end subroutine

  subroutine nwp_seq2ymdh(y, m, d, h, iseq)
    integer, intent(out):: y, m, d, h
    integer, intent(in):: iseq
    integer:: residual
    call nwp_seq2ym(y, m, iseq, residual)
    d = residual / DAY
    residual = mod(residual, 60)
    h = residual / 60
  end subroutine

  subroutine nwp_seq2ymd(y, m, d, iseq)
    integer, intent(out):: y, m, d
    integer, intent(in):: iseq
    integer:: residual
    call nwp_seq2ym(y, m, iseq, residual)
    d = residual / DAY
  end subroutine

  ! CURRENTLY PRIVATE
  ! WORKS YEARS FROM -2282 TO 5883.
  subroutine nwp_seq2ym(y, m, iseq, residual)
    integer, intent(out):: y, m, residual
    integer, intent(in):: iseq
    integer, parameter:: MINS_PER_400YEARS = (36524 * 4 + 1) * DAY
    integer, parameter:: DAYS_PER_CENTURY = 36524
    integer, parameter:: DAYS_PER_OLYMPIAD = 365 * 4 + 1
    integer, parameter:: ORIGIN2K = 104749920 ! 2000-03-01
    integer:: d_res, ic, iy, r2
    !
    ! ADJUST ALL DATES INTO 2000-03-01/2400-02-29
    !
    if (iseq < ORIGIN2K) then
      residual = iseq + MINS_PER_400YEARS - ORIGIN2K
      r2 = modulo(residual, MINS_PER_400YEARS)
      y = 1600 - 400 * ((r2 - residual) / MINS_PER_400YEARS)
      residual = r2
    else
      ! y = 2000
      residual = iseq - ORIGIN2K
      y = 2000 + 400 * (residual / MINS_PER_400YEARS)
      residual = modulo(residual, MINS_PER_400YEARS)
    endif
    !
    ! THEN THE OPERATION UNIT BECOMES DAY
    !
    d_res = residual / DAY
    residual = mod(residual, DAY)
    !
    ! DECIDE CENTURY ... 3rd CENTURY IS LONGER BY 1 DAY
    !
    ic = (d_res / DAYS_PER_CENTURY) * 100
    d_res = mod(d_res, DAYS_PER_CENTURY)
    if (ic == 400) then
      ic = ic - 100
      d_res = d_res + DAYS_PER_CENTURY
    endif
    y = y + ic
    !
    ! DECIDE OLYMPIAD
    !
    y = y + (d_res / DAYS_PER_OLYMPIAD) * 4
    d_res = mod(d_res, DAYS_PER_OLYMPIAD)
    !
    ! DECIDE YEAR ... YEAR MAY BE LONGER THAN 365 DAYS
    !
    iy = d_res / 365
    d_res = mod(d_res, 365)
    if (iy == 4) then
      iy = 3
      d_res = d_res + 365
    endif
    y = y + iy
    !
    ! DECIDE MONTH
    !
    m = (d_res * 10 + 4) / 306 + 3
    d_res = d_res - (306 * m + 6) / 10 + 93
    residual = residual + d_res * DAY
    if (m > 12) then
      m = m - 12
      y = y + 1
    endif
  end subroutine

end module