!
!  $Author: tomita $
!  $Date: 2007/08/01 20:09:58 $
!  $Revision: 1.1.1.1 $
!
MODULE InputParameters

   IMPLICIT NONE

   PRIVATE

   INTEGER, PARAMETER, PUBLIC :: &
            r4 = SELECTED_REAL_KIND(6)  ! Kind for 32-bits Real Numbers
   INTEGER, PARAMETER, PUBLIC :: &
            r8 = SELECTED_REAL_KIND(15) ! Kind for 64-bits Real Numbers
   INTEGER, PARAMETER,PUBLIC  :: i4 = SELECTED_INT_KIND(9)   ! Kind for 32-bits Integer Numbers

   REAL (KIND=r8), PUBLIC :: Undef, Lon0, Lat0

   INTEGER, PUBLIC :: Imax, Jmax, Idim_sst, Jdim_sst,kdim,Idim_xice, Jdim_xice

   LOGICAL, PUBLIC :: PolarMean, GrADS, Linear

   LOGICAL, PUBLIC :: FlagInput(5), FlagOutput(5)

   INTEGER, DIMENSION (:,:), ALLOCATABLE, PUBLIC :: MaskInput

   CHARACTER (LEN=7), PUBLIC :: nLats='.G     '

   CHARACTER (LEN=11),PUBLIC :: DirPreIn='pre/datain/'

   CHARACTER (LEN=12), PUBLIC :: DirPreOut='pre/dataout/'

   CHARACTER (LEN=13), PUBLIC :: DirModelIn='model/datain/'

   CHARACTER (LEN=18), PUBLIC :: VarNameOut='AOI_SSTNOAA'

   CHARACTER (LEN=32), PUBLIC :: VarName

   CHARACTER (LEN=528), PUBLIC :: DirMainInput, DirMainOutputmon,DirMainOutputwkl

   CHARACTER (LEN=10), PUBLIC :: LabelI='YYYYMMDDHH'
   CHARACTER (LEN=10), PUBLIC :: LabelF='YYYYMMDDHH'

   CHARACTER (LEN=64) :: NameNML='AOI_SSTNOAA.nml'

   INTEGER, PUBLIC :: nferr=0    ! Standard Error Print Out
   INTEGER, PUBLIC :: nfinp=5    ! Standard Read In
   INTEGER, PUBLIC :: nfprt=6    ! Standard Print Out
   INTEGER, PUBLIC :: nfclm=10   ! To Read Climatological Roughness Length Data
   INTEGER, PUBLIC :: nfout=20   ! To Write Intepolated Climatological Roughness Length Data
   INTEGER, PUBLIC :: nfctl=30   ! To Write Output Data Description
   INTEGER, PUBLIC :: nfvgm=40   ! To Write Unformatted Gaussian Grid Vegetation Mask

   PUBLIC :: InitInputParameters


CONTAINS


SUBROUTINE InitInputParameters ()

   IMPLICIT NONE

   INTEGER :: ios

   NAMELIST /InputDim/ Imax, Jmax,kdim,LabelI,LabelF, &
                       GrADS, Linear, VarName, DirMainInput, DirMainOutputmon,DirMainOutputwkl

   Imax=192
   Jmax=96
   Idim_sst=0
   Jdim_sst=0
   Idim_xice=0
   Jdim_xice=0
   kdim=8
   GrADS=.TRUE.
   Linear=.TRUE.
   VarName='AOI_SSTNOAA'
   LabelI='          '
   LabelF='          '
   DirMainInput='./ '
   DirMainOutputmon='./ '
   DirMainOutputwkl='./ '

   OPEN (UNIT=nfinp, FILE='./'//TRIM(NameNML), &
         FORM='FORMATTED', ACCESS='SEQUENTIAL', &
         ACTION='READ', STATUS='OLD', IOSTAT=ios)
   IF (ios /= 0) THEN
      WRITE (UNIT=nferr, FMT='(3A,I4)') &
            ' ** (Error) ** Open file ', &
              './'//TRIM(NameNML), &
            ' returned IOStat = ', ios
      STOP  ' ** (Error) **'
   END IF
   READ  (UNIT=nfinp, NML=InputDim)
   CLOSE (UNIT=nfinp)

   WRITE (UNIT=nfprt, FMT='(/,A)')  ' &InputDim'
   WRITE (UNIT=nfprt, FMT='(A,I6)') '     Imax = ', Imax
   WRITE (UNIT=nfprt, FMT='(A,I6)') '     Jmax = ', Jmax
   WRITE (UNIT=nfprt, FMT='(A,I6)') '     kdim = ', kdim
   WRITE (UNIT=nfprt, FMT='(A,L6)') '    GrADS = ', GrADS
   WRITE (UNIT=nfprt, FMT='(A,L6)') '   Linear = ', Linear
   WRITE (UNIT=nfprt, FMT='(A)')    '  VarName = '//TRIM(VarName)
   WRITE (UNIT=nfprt, FMT='(A)')    '  LabelI = '//TRIM(LabelI)
   WRITE (UNIT=nfprt, FMT='(A)')    '  LabelF = '//TRIM(LabelF)
   WRITE (UNIT=nfprt, FMT='(A)')    '  DirMainInput = '//TRIM(DirMainInput)
   WRITE (UNIT=nfprt, FMT='(A)')    '  DirMainOutputwkl = '//TRIM(DirMainOutputwkl)
   WRITE (UNIT=nfprt, FMT='(A)')    '  DirMainOutputmon = '//TRIM(DirMainOutputmon)

   WRITE (UNIT=nfprt, FMT='(A,/)')  ' /'

   Undef=-999.0_r8

   ! For Linear Interpolation
   Lon0=0.0_r8  ! Start at Prime Meridian
   Lat0=90.0_r8 ! Start at North Pole

   ! For Area Weighted Interpolation
   ALLOCATE (MaskInput(Idim_sst,Jdim_sst))
   MaskInput=1
   PolarMean=.FALSE.
   FlagInput(1)=.TRUE.   ! Start at North Pole
   FlagInput(2)=.TRUE.   ! Start at Prime Meridian
   FlagInput(3)=.TRUE.   ! Latitudes Are at Center of Box
   FlagInput(4)=.TRUE.   ! Longitudes Are at Center of Box
   FlagInput(5)=.FALSE.  ! Regular Grid
   FlagOutput(1)=.TRUE.  ! Start at North Pole
   FlagOutput(2)=.TRUE.  ! Start at Prime Meridian
   FlagOutput(3)=.FALSE. ! Latitudes Are at North Edge of Box
   FlagOutput(4)=.TRUE.  ! Longitudes Are at Center of Box
   FlagOutput(5)=.TRUE.  ! Gaussian Grid

   WRITE (nLats(3:7), '(I5.5)') Jmax

END SUBROUTINE InitInputParameters

END MODULE InputParameters

!*****************************************************************************************
!> author: Jacob Williams
!
!  Conversion factors.
!
!# See also
!
!  1. A. Thompson and B. N. Taylor, "NIST Special Publication 811:
!     Guide for the use of the International System of Units".
!     http://www.nist.gov/pml/pubs/sp811/

    module conversion_module

    implicit none

    INTEGER, PARAMETER, PRIVATE :: &
            r4 = SELECTED_REAL_KIND(6)  ! Kind for 32-bits Real Numbers
    INTEGER, PARAMETER, PRIVATE :: &
            r8 = SELECTED_REAL_KIND(15) ! Kind for 64-bits Real Numbers
    INTEGER, PARAMETER,PRIVATE  :: i4 = SELECTED_INT_KIND(9)   ! Kind for 32-bits Integer Numbers

    real(r8),parameter,PRIVATE :: one        = 1.0_r8
    real(r8),parameter,PRIVATE :: pi         = acos(-one)
    public

    !metric/imperial:
    real(r8),parameter :: lbm2kg  = 0.45359237_r8         !! exact
    real(r8),parameter :: lbf2N   = 4.4482216152605_r8    !! exact
    real(r8),parameter :: ft2m    = 0.3048_r8             !! exact
    real(r8),parameter :: mile2km = 1.609344_r8           !! exact
    real(r8),parameter :: nmi2km  = 1.852_r8              !! exact
    real(r8),parameter :: slug2kg = lbf2N/ft2m     !! approximately 14.593902937206362
    real(r8),parameter :: kg2lbm  = one/lbm2kg     !! approximately 2.2046226218487757
    real(r8),parameter :: N2lbf   = one/lbf2N      !! approximately 0.2248089430997105
    real(r8),parameter :: m2ft    = one/ft2m       !! approximately 3.280839895013123
    real(r8),parameter :: km2mile = one/mile2km    !! approximately 0.621371192237334
    real(r8),parameter :: km2nmi  = one/nmi2km     !! approximately 0.5399568034557235
    real(r8),parameter :: kg2slug = ft2m/lbf2N     !! approximately 0.06852176585679176

    !angles:
    real(r8),parameter :: deg2rad = pi/180.0_r8
    real(r8),parameter :: rad2deg = 180.0_r8/pi

    !metric:
    real(r8),parameter :: km2m = 1000.0_r8
    real(r8),parameter :: m2km = one/km2m
    real(r8),parameter :: au2m = 149597870700.0_r8 !! IAU 2012 defined value

    !time:
    real(r8),parameter :: min2sec      = 60.0_r8
    real(r8),parameter :: hr2min       = 60.0_r8
    real(r8),parameter :: day2hr       = 24.0_r8
    real(r8),parameter :: year2day     = 365.25_r8         !! julian year
    real(r8),parameter :: century2day  = year2day*100.0_r8 !! julian century
    real(r8),parameter :: deg2arcmin   = 60.0_r8
    real(r8),parameter :: deg2arcsec   = 3600.0_r8
    real(r8),parameter :: hr2sec       = hr2min*min2sec
    real(r8),parameter :: day2min      = day2hr*hr2min
    real(r8),parameter :: day2sec      = day2min*min2sec
    real(r8),parameter :: century2sec  = century2day*day2sec
    real(r8),parameter :: day2year     = one/year2day
    real(r8),parameter :: day2century  = one/century2day
    real(r8),parameter :: hr2day       = one/day2hr
    real(r8),parameter :: sec2hr       = one/hr2sec
    real(r8),parameter :: sec2day      = one/day2sec
    real(r8),parameter :: sec2century  = one/century2sec
    real(r8),parameter :: arcmin2deg   = one/deg2arcmin
    real(r8),parameter :: arcsec2deg   = one/deg2arcsec
CONTAINS 

  Subroutine Init_conversion_module()
  END Subroutine Init_conversion_module

END module conversion_module
    

!*****************************************************************************************
!> author: Jacob Williams
!
!  Time conversion routines.

    module time_module

    implicit none

    private
   INTEGER, PARAMETER, PUBLIC :: &
            r4 = SELECTED_REAL_KIND(6)  ! Kind for 32-bits Real Numbers
   INTEGER, PARAMETER, PUBLIC :: &
            r8 = SELECTED_REAL_KIND(15) ! Kind for 64-bits Real Numbers
   INTEGER, PARAMETER,PUBLIC  :: i4 = SELECTED_INT_KIND(9)   ! Kind for 32-bits Integer Numbers

    !parameters:
    real(r8),parameter :: jd_j2000 = 2451545.0_r8  !! julian date of J2000 epoch

    interface julian_date
        !! calendar date to julian date
        module procedure :: julian_date_realsec, &
                            julian_date_intsec
    end interface

    interface julian_date_to_calendar_date
        module procedure :: calendar_date_realsec
    end interface

    !public routines:
    public :: julian_day
    public :: weekday_nm
    public :: julian_date
    public :: et_to_jd
    public :: jd_to_et
    public :: jd_to_mjd
    public :: mjd_to_jd
    public :: julian_date_to_calendar_date

    !test routine:
    public :: time_module_test

    contains

    pure integer function weekday_nm(y,m,d)
    ! Returns the day of the week calculated using Zeller's congruence.
    ! Returned value is an integer scalar in the range [0-6], such that:
    !
    ! 0: Sunday
    ! 1: Monday
    ! 2: Tuesday
    ! 3: Wednesday
    ! 4: Thursday
    ! 5: Friday
    ! 6: Saturday
    integer,intent(in) :: y   !! year (YYYY)
    integer,intent(in) :: m   !! month (MM)
    integer,intent(in) :: d   !! day (DD)

!    class(datetime), intent(in) :: self
    integer :: year, month,day, j, k

    year  = y
    month = m
    day   = d
    if (month <= 2) then
      month = month + 12
      year  = year - 1
    end if

    j = year / 100
    k = mod(year, 100)

    weekday_nm = mod( day + ((month + 1) * 26) / 10 + k + k / 4 + j / 4 + 5 * j, 7) -1

    if (weekday_nm < 0) weekday_nm = 6

  end function weekday_nm

!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 2/3/2015
!
!  Convert ephemeris time (seconds from J2000 epoch) to Julian date.

    pure function et_to_jd(et) result(jd)

    use conversion_module, only: sec2day

    implicit none

    real(r8),intent(in) :: et   !! ephemeris time [sec from J2000 epoch]
    real(r8)            :: jd   !! Julian date [days]

    jd = jd_j2000 + et*sec2day

    end function et_to_jd
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 3/19/2016
!
!  Convert Julian date to ephemeris time (seconds from J2000 epoch).

    pure function jd_to_et(jd) result(et)

    use conversion_module, only: day2sec

    implicit none

    real(r8),intent(in) :: jd   !! Julian date [days]
    real(r8)            :: et   !! ephemeris time [sec from J2000 epoch]

    et = (jd - jd_j2000) * day2sec

    end function jd_to_et
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 3/15/2015
!
!  Converts Julian date to Modified Julian date.
!
!### Reference
!   * [USNO](http://tycho.usno.navy.mil/mjd.html)

    pure function jd_to_mjd(jd) result(mjd)

    implicit none

    real(r8)            :: mjd  !! modified julian date
    real(r8),intent(in) :: jd   !! julian date

    mjd = jd - 2400000.5_r8

    end function jd_to_mjd
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 3/15/2015
!
!  Converts Modified Julian date to Julian date.
!
!### Reference
!   * [USNO](http://tycho.usno.navy.mil/mjd.html)

    pure function mjd_to_jd(mjd) result(jd)

    implicit none

    real(r8)            :: jd   !! julian date
    real(r8),intent(in) :: mjd  !! modified julian date

    jd = mjd + 2400000.5_r8

    end function mjd_to_jd
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!
!  Returns the Julian day number (i.e., the Julian date at Greenwich noon)
!  on the specified YEAR, MONTH, and DAY.
!
!  Valid for any Gregorian calendar date producing a
!  Julian date greater than zero.
!
!### Reference
!   * [USNO](http://aa.usno.navy.mil/faq/docs/JD_Formula.php)

    pure integer function julian_day(y,m,d)

    implicit none

    integer,intent(in) :: y   !! year (YYYY)
    integer,intent(in) :: m   !! month (MM)
    integer,intent(in) :: d   !! day (DD)
    integer    :: julian_year   !! day (DD)
    integer    :: julian_full   !! day (DD)

    julian_year= 1-32075+1461*(y+4800+(1-14)/12)/4+367*&
                 (1-2-(1-14)/12*12)/12-3*((y+4900+(1-14)/12)/100)/4

    julian_full = d-32075+1461*(y+4800+(m-14)/12)/4+367*&
                 (m-2-(m-14)/12*12)/12-3*((y+4900+(m-14)/12)/100)/4

     julian_day=julian_full-julian_year + 1
    end function julian_day
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 1/21/2015
!
!  Returns the Julian date for the specified YEAR, MONTH, DAY, HR, MIN, SEC.
!
!  Valid for any Gregorian calendar date producing a
!  Julian date greater than zero.
!
!### History
!  * JW : 10/4/2017 : moved main code to [[julian_date_realsec]] routine.

    pure function julian_date_intsec(y,m,d,hour,minute,second) result(julian_date)

    implicit none

    real(r8)           :: julian_date
    integer,intent(in) :: y
    integer,intent(in) :: m
    integer,intent(in) :: d
    integer,intent(in) :: hour
    integer,intent(in) :: minute
    integer,intent(in) :: second

    ! call the other routine:
    julian_date = julian_date_realsec(y,m,d,hour,minute,real(second,r8))

    end function julian_date_intsec
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 1/21/2015
!
!  Returns the Julian date for the specified YEAR, MONTH, DAY, HR, MIN, SEC.
!
!  Valid for any Gregorian calendar date producing a
!  Julian date greater than zero.
!
!### History
!  * JW : 10/4/2017 : made `second` a real value & renamed routine.

    pure function julian_date_realsec(y,m,d,hour,minute,second) result(julian_date)

    implicit none

    real(r8)            :: julian_date
    integer,intent(in)  :: y
    integer,intent(in)  :: m
    integer,intent(in)  :: d
    integer,intent(in)  :: hour
    integer,intent(in)  :: minute
    real(r8),intent(in) :: second

    integer :: julian_day_number

    julian_day_number = julian_day(y,m,d)

    julian_date = real(julian_day_number,r8) + &
                    (hour-12.0_r8)/24.0_r8 + &
                    minute/1440.0_r8 + &
                    second/86400.0_r8

    end function julian_date_realsec
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!
!  Returns the year, month, day, hr, min, sec for the specified Julian date.
!
!### See also
!  * https://aa.usno.navy.mil/faq/JD_formula.html
!  * http://www.davidgsimpson.com/software/jd2greg_f90.txt

    pure subroutine calendar_date_realsec(julian_date,year,month,day,hrs,min,sec)

    implicit none

    real(r8),intent(in)  :: julian_date !! julian date
    integer,intent(out)  :: year
    integer,intent(out)  :: month
    integer,intent(out)  :: day
    integer,intent(out)  :: hrs
    integer,intent(out)  :: min
    real(r8),intent(out) :: sec

    integer :: i,j,k,l,n,jd
    real(r8) :: frac_day

    jd = int(julian_date)

    l = jd+68569
    n = 4*l/146097
    l = l-(146097*n+3)/4
    i = 4000*(l+1)/1461001
    l = l-1461*i/4+31
    j = 80*l/2447
    k = l-2447*j/80
    l = j/11
    j = j+2-12*l
    i = 100*(n-49)+i+l

    year  = i
    month = j
    day   = k

    frac_day = julian_date - real(jd,r8) + 0.5_r8

    hrs = int(frac_day*24.0_r8)
    min = int((frac_day - hrs/24.0_r8) * 1440.0_r8)
    sec = (frac_day - hrs/24.0_r8 - min/1440.0_r8) * 86400.0_r8

    if (sec == 60.0_r8) then
        sec = 0.0_r8
        min = min + 1
     end if

     if (min == 60) then
        min = 0
        hrs = hrs + 1
     end if

    end subroutine calendar_date_realsec
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 1/21/2015
!
!  Test routine for the Julian date routines.

    subroutine time_module_test()

    implicit none

    real(r8) :: jd, sec
    integer :: year,month,day,hrs,min

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' time_module_test'
    write(*,*) '---------------'
    write(*,*) ''

    ! JD = 2451545.0
    jd = julian_date(2000,1,1,12,0,0)

    call calendar_date_realsec(jd,year,month,day,hrs,min,sec)

    write(*,*) 'jd    ', jd
    write(*,*) 'year  ', year
    write(*,*) 'month ', month
    write(*,*) 'day   ', day
    write(*,*) 'hrs   ', hrs
    write(*,*) 'min   ', min
    write(*,*) 'sec   ', sec

    if (year/=2000)  error stop 'error: incorrect year'
    if (month/=1)    error stop 'error: incorrect month'
    if (day/=1)      error stop 'error: incorrect day'
    if (hrs/=12)     error stop 'error: incorrect hrs'
    if (min/=0)      error stop 'error: incorrect min'
    if (sec/=0.0_r8) error stop 'error: incorrect sec'

    write(*,*) 'PASSED'

    end subroutine time_module_test
!*****************************************************************************************
    end module time_module
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************
!*****************************************************************************************


module datetime_module

  use, intrinsic :: iso_fortran_env, only: int64, real32, real64, &
                                           stderr => error_unit
  use, intrinsic :: iso_c_binding, only: c_char, c_int, c_null_char

  implicit none

  private

  public :: datetime, timedelta, clock
  public :: date2num
  public :: datetimeRange
  public :: daysInMonth
  public :: daysInYear
  public :: isLeapYear
  public :: num2date
  public :: strptime
  public :: tm2date
  public :: tm_struct
  public :: c_strftime
  public :: c_strptime

  real(real64), parameter :: zero = 0._real64, one = 1._real64

  ! Constant multipliers to transform a number of some time unit to another
  real(real64), parameter :: d2h = 24._real64 ! day -> hour
  real(real64), parameter :: h2d = one / d2h ! hour -> day
  real(real64), parameter :: d2m = d2h * 60._real64 ! day -> minute
  real(real64), parameter :: m2d = one / d2m ! minute -> day
  real(real64), parameter :: m2h = one / 60 ! minute -> hour
  real(real64), parameter :: s2d = m2d / 60 ! second -> day
  real(real64), parameter :: d2s = 86400._real64 ! day -> second
  real(real64), parameter :: h2s = 3600._real64 ! hour -> second
  real(real64), parameter :: s2h = one / h2s ! second -> hour
  real(real64), parameter :: m2s = 60._real64 ! minute -> second
  real(real64), parameter :: s2m = one / m2s ! second -> minute

  integer, parameter :: MAXSTRLEN = 99 ! maximum string length for strftime

  type :: datetime

    private

    integer :: year        = 1 ! year [1-HUGE(year)]
    integer :: month       = 1 ! month in year [1-12]
    integer :: day         = 1 ! day in month [1-31]
    integer :: hour        = 0 ! hour in day [0-23]
    integer :: minute      = 0 ! minute in hour [0-59]
    integer :: second      = 0 ! second in minute [0-59]
    integer :: millisecond = 0 ! milliseconds in second [0-999]
    real(real64) :: tz = 0 ! timezone offset from UTC [hours]

  contains

    ! getter functions
    procedure, pass(self), public :: getYear
    procedure, pass(self), public :: getMonth
    procedure, pass(self), public :: getDay
    procedure, pass(self), public :: getHour
    procedure, pass(self), public :: getMinute
    procedure, pass(self), public :: getSecond
    procedure, pass(self), public :: getMillisecond
    procedure, pass(self), public :: getTz

    ! public methods
    procedure, pass(self), public :: isocalendar
    procedure, pass(self), public :: isoformat
    procedure, pass(self), public :: isValid
    procedure, nopass,     public :: now
    procedure, pass(self), public :: secondsSinceEpoch
    procedure, pass(self), public :: strftime
    procedure, pass(self), public :: tm
    procedure, pass(self), public :: tzOffset
    procedure, pass(self), public :: isoweekday
    procedure, pass(self), public :: isoweekdayLong
    procedure, pass(self), public :: isoweekdayShort
    procedure, pass(self), public :: utc
    procedure, pass(self), public :: weekday
    procedure, pass(self), public :: weekdayLong
    procedure, pass(self), public :: weekdayShort
    procedure, pass(self), public :: yearday

    ! private methods
    procedure, pass(self), private :: addMilliseconds
    procedure, pass(self), private :: addSeconds
    procedure, pass(self), private :: addMinutes
    procedure, pass(self), private :: addHours
    procedure, pass(self), private :: addDays

    ! operator overloading procedures
    procedure, pass(d0), private :: datetime_plus_timedelta
    procedure, pass(d0), private :: timedelta_plus_datetime
    procedure, pass(d0), private :: datetime_minus_datetime
    procedure, pass(d0), private :: datetime_minus_timedelta
    procedure, pass(d0), private :: datetime_eq
    procedure, pass(d0), private :: datetime_neq
    procedure, pass(d0), private :: datetime_gt
    procedure, pass(d0), private :: datetime_ge
    procedure, pass(d0), private :: datetime_lt
    procedure, pass(d0), private :: datetime_le

    generic :: operator(+)  => datetime_plus_timedelta, &
                               timedelta_plus_datetime
    generic :: operator(-)  => datetime_minus_datetime, &
                               datetime_minus_timedelta
    generic :: operator(==) => datetime_eq
    generic :: operator(/=) => datetime_neq
    generic :: operator(>)  => datetime_gt
    generic :: operator(>=) => datetime_ge
    generic :: operator(<)  => datetime_lt
    generic :: operator(<=) => datetime_le

  end type datetime

  interface datetime
    module procedure :: datetime_constructor
  endinterface datetime

  type :: timedelta
    private

    integer :: days = 0
    integer :: hours = 0
    integer :: minutes = 0
    integer :: seconds = 0
    integer :: milliseconds = 0

  contains

    procedure, pass(self), public :: getDays
    procedure, pass(self), public :: getHours
    procedure, pass(self), public :: getMinutes
    procedure, pass(self), public :: getSeconds
    procedure, pass(self), public :: getMilliseconds

    procedure, public :: total_seconds

    procedure, private :: timedelta_plus_timedelta
    procedure, private :: timedelta_minus_timedelta
    procedure, private :: unary_minus_timedelta
    procedure, private :: timedelta_eq
    procedure, private :: timedelta_neq
    procedure, private :: timedelta_gt
    procedure, private :: timedelta_ge
    procedure, private :: timedelta_lt
    procedure, private :: timedelta_le

    generic :: operator(+)  => timedelta_plus_timedelta
    generic :: operator(-)  => timedelta_minus_timedelta, unary_minus_timedelta
    generic :: operator(==) => timedelta_eq
    generic :: operator(/=) => timedelta_neq
    generic :: operator(>)  => timedelta_gt
    generic :: operator(>=) => timedelta_ge
    generic :: operator(<)  => timedelta_lt
    generic :: operator(<=) => timedelta_le

  end type timedelta

  interface timedelta
    module procedure :: timedelta_constructor
  endinterface timedelta

  type,bind(c) :: tm_struct
    ! Derived type for compatibility with C and C++ struct tm.
    ! Enables calling strftime and strptime using iso_c_binding.
    ! See http://www.cplusplus.com/reference/ctime/tm for reference.
    integer(c_int) :: tm_sec = 0 ! Seconds [0-60] (1 leap second)
    integer(c_int) :: tm_min = 0 ! Minutes [0-59]
    integer(c_int) :: tm_hour = 0 ! Hours [0-23]
    integer(c_int) :: tm_mday = 0 ! Day [1-31]
    integer(c_int) :: tm_mon = 0 ! Month [0-11]
    integer(c_int) :: tm_year = 0 ! Year - 1900
    integer(c_int) :: tm_wday = 0 ! Day of week [0-6]
    integer(c_int) :: tm_yday = 0 ! Days in year [0-365]
    integer(c_int) :: tm_isdst = 0 ! DST [-1/0/1]
  end type tm_struct

  interface

    function c_strftime(str, slen, format, tm) bind(c, name='strftime') result(rc)
      ! Returns a formatted time string, given input time struct and format.
      ! See https://www.cplusplus.com/reference/ctime/strftime for reference.
      import :: c_char, c_int, tm_struct
      character(kind=c_char), intent(out) :: str(*) ! result string
      integer(c_int), value, intent(in) :: slen ! string length
      character(kind=c_char), intent(in) :: format(*) ! time format
      type(tm_struct), intent(in) :: tm ! tm_struct instance
      integer(c_int) :: rc ! return code
    end function c_strftime

    function c_strptime(str,format,tm) bind(c,name='strptime') result(rc)
      ! Interface to POSIX strptime.
      ! Returns a time struct object based on the input time string str and format.
      ! See http://man7.org/linux/man-pages/man3/strptime.3.html for reference.
      import :: c_char, c_int, tm_struct
      character(kind=c_char), intent(in) :: str(*) ! input string
      character(kind=c_char), intent(in) :: format(*) ! time format
      type(tm_struct), intent(out) :: tm ! result tm_struct
      integer(c_int) :: rc ! return code
    end function c_strptime

  end interface


  type :: clock
    type(datetime) :: startTime
    type(datetime) :: stopTime
    type(datetime) :: currentTime
    type(timedelta) :: tickInterval
    logical :: alarm = .false.
    logical :: started = .false.
    logical :: stopped = .false.
  contains
    procedure :: reset
    procedure :: tick
  end type clock

contains


  pure elemental subroutine reset(self)
    ! Resets the clock to its start time.
    class(clock), intent(in out) :: self
    self % currentTime = self % startTime
    self % started = .false.
    self % stopped = .false.
  end subroutine reset


  pure elemental subroutine tick(self)
    ! Increments the currentTime of the clock instance by one tickInterval.
    class(clock), intent(in out) :: self
    if (self % stopped) return
    if (.not. self % started) then
      self % started = .true.
      self % currentTime = self % startTime
    end if
    self % currentTime = self % currentTime + self % tickInterval
    if (self % currentTime >= self % stopTime) self % stopped = .true.
  end subroutine tick


  pure elemental type(datetime) function datetime_constructor( &
    year, month, day, hour, minute, second, millisecond, tz)
    ! Constructor function for the `datetime` class.
    integer, intent(in), optional :: year, month, day, hour, minute, second, millisecond
    real(real64), intent(in), optional :: tz ! timezone offset in hours

    datetime_constructor % year = 1
    if (present(year)) datetime_constructor % year = year

    datetime_constructor % month = 1
    if (present(month)) datetime_constructor % month = month

    datetime_constructor % day = 1
    if (present(day)) datetime_constructor % day = day

    datetime_constructor % hour = 0
    if (present(hour)) datetime_constructor % hour = hour

    datetime_constructor % minute = 0
    if (present(minute)) datetime_constructor % minute = minute

    datetime_constructor % second = 0
    if (present(second)) datetime_constructor % second = second

    datetime_constructor % millisecond = 0
    if (present(millisecond)) datetime_constructor % millisecond = millisecond

    datetime_constructor % tz = 0
    if (present(tz)) datetime_constructor % tz = tz

  end function datetime_constructor


  pure elemental integer function getYear(self)
    ! Returns the year component
    class(datetime), intent(in) :: self
    getYear = self % year
  end function getYear


  pure elemental integer function getMonth(self)
    ! Returns the year component
    class(datetime), intent(in) :: self
    getMonth = self % month
  end function getMonth


  pure elemental integer function getDay(self)
    ! Returns the year component
    class(datetime), intent(in) :: self
    getDay = self % day
  end function getDay


  pure elemental integer function getHour(self)
    ! Returns the year component
    class(datetime), intent(in) :: self
    getHour = self % hour
  end function getHour


  pure elemental integer function getMinute(self)
    ! Returns the year component
    class(datetime), intent(in) :: self
    getMinute = self % minute
  end function getMinute


  pure elemental integer function getSecond(self)
    ! Returns the year component
    class(datetime), intent(in) :: self
    getSecond = self % second
  end function getSecond


  pure elemental integer function getMillisecond(self)
    ! Returns the year component
    class(datetime), intent(in) :: self
    getMillisecond = self % millisecond
  end function getMillisecond


  pure elemental real(real64) function getTz(self)
    ! Returns the timezone offset component
    class(datetime), intent(in) :: self
    getTz = self % tz
  end function getTz


  pure elemental subroutine addMilliseconds(self, ms)
    ! Adds an integer number of milliseconds to self. Called by `datetime`
    ! addition (`+`) and subtraction (`-`) operators.
    class(datetime), intent(in out) :: self
    integer, intent(in) :: ms
    self % millisecond = self % millisecond + ms
    do
      if (self % millisecond >= 1000) then
        call self % addSeconds(self % millisecond / 1000)
        self % millisecond = mod(self % millisecond, 1000)
      else if (self % millisecond < 0) then
        call self % addSeconds(self % millisecond / 1000 - 1)
        self % millisecond = mod(self % millisecond, 1000) + 1000
      else
        exit
      end if
    end do
  end subroutine addMilliseconds


  pure elemental subroutine addSeconds(self, s)
    ! Adds an integer number of seconds to self. Called by `datetime`
    ! addition (`+`) and subtraction (`-`) operators.
    class(datetime), intent(in out) :: self
    integer, intent(in) :: s
    self % second = self % second + s
    do
      if (self % second >= 60) then
        call self % addMinutes(self % second / 60)
        self % second = mod(self % second, 60)
      else if (self % second < 0) then
        call self % addMinutes(self % second / 60 - 1)
        self % second = mod(self % second, 60) + 60
      else
        exit
      end if
    end do
  end subroutine addSeconds


  pure elemental subroutine addMinutes(self,m)
    ! Adds an integer number of minutes to self. Called by `datetime`
    ! addition (`+`) and subtraction (`-`) operators.
    class(datetime), intent(in out) :: self
    integer, intent(in) :: m
    self % minute = self % minute + m
    do
      if (self % minute >= 60) then
        call self % addHours(self % minute / 60)
        self % minute = mod(self % minute, 60)
      else if (self % minute < 0) then
        call self % addHours(self % minute / 60 - 1)
        self % minute = mod(self % minute, 60) + 60
      else
        exit
      end if
    end do
  end subroutine addMinutes


  pure elemental subroutine addHours(self,h)
    ! Adds an integer number of hours to self. Called by `datetime`
    ! addition (`+`) and subtraction (`-`) operators.
    class(datetime), intent(in out) :: self
    integer, intent(in) :: h
    self % hour = self % hour + h
    do
      if (self % hour >= 24) then
        call self % addDays(self % hour / 24)
        self % hour = mod(self % hour, 24)
      else if (self % hour < 0) then
        call self % addDays(self % hour / 24 - 1)
        self % hour = mod(self % hour, 24) + 24
      else
        exit
      end if
    end do
  end subroutine addHours


  pure elemental subroutine addDays(self, d)
    ! Adds an integer number of dayss to self. Called by `datetime`
    ! addition (`+`) and subtraction (`-`) operators.
    class(datetime), intent(in out) :: self
    integer, intent(in) :: d
    integer :: daysInCurrentMonth
    self % day = self % day + d
    do
      daysInCurrentMonth = daysInMonth(self % month, self % year)
      if (self % day > daysInCurrentMonth) then
        self % day = self % day - daysInCurrentMonth
        self % month = self % month+1
        if (self % month > 12) then
          self % year = self % year + self % month/12
          self % month = mod(self % month, 12)
        end if
      else if (self % day < 1) then
        self % month = self % month-1
        if (self % month < 1) then
          self % year = self % year + self % month / 12 - 1
          self % month = 12 + mod(self % month, 12)
        end if
        self % day = self % day + daysInMonth(self % month, self % year)
      else
        exit
      end if
    end do
  end subroutine addDays


  pure elemental character(23) function isoformat(self,sep)
    ! Returns character string with time in ISO 8601 format.
    class(datetime), intent(in) :: self
    character, intent(in), optional :: sep
    character :: separator

    separator = 'T'
    if (present(sep)) separator = sep

    ! TODO below is a bit cumbersome and was implemented
    ! at a time before the interface to strftime. Now we
    ! could do something like:
    !
    ! isoformat = self % strftime('%Y-%m-%d'//separator//'%H:%M:%S')
    !
    isoformat = int2str(self % year,       4)//'-'//      &
                int2str(self % month,      2)//'-'//      &
                int2str(self % day,        2)//separator//&
                int2str(self % hour,       2)//':'//      &
                int2str(self % minute,     2)//':'//      &
                int2str(self % second,     2)//'.'//      &
                int2str(self % millisecond,3)

  end function isoformat


  pure elemental logical function isValid(self)
    ! Checks whether the `datetime` instance has valid component values.
    ! Returns `.true.` if the `datetime` instance is valid, and `.false.`
    ! otherwise.
    class(datetime), intent(in) :: self

    ! assume valid
    isValid = .true.

    if (self % year < 1) then
      isValid = .false.
      return
    end if

    if (self % month < 1 .or. self % month > 12) then
      isValid = .false.
      return
    end if

    if (self % day < 1 .or. &
       self % day > daysInMonth(self % month,self % year)) then
      isValid = .false.
      return
    end if

    if (self % hour < 0 .or. self % hour > 23) then
      isValid = .false.
      return
    end if

    if (self % minute < 0 .or. self % minute > 59) then
      isValid = .false.
      return
    end if

    if (self % second < 0 .or. self % second > 59) then
      isValid = .false.
      return
    end if

    if (self % millisecond < 0 .or. self % millisecond > 999) then
      isValid = .false.
      return
    end if

  end function isValid


  type(datetime) function now()
    ! Returns a `datetime` instance with current time.
    character(5) :: zone
    integer :: values(8)
    integer :: hour, minute

    ! Obtain local machine time zone information
    call date_and_time(zone=zone, values=values)

    read(zone(1:3), '(i3)') hour
    read(zone(4:5), '(i2)') minute

    now = datetime(year = values(1), month = values(2), day = values(3), &
                   hour = values(5), minute = values(6), second = values(7), &
                   millisecond = values(8))

    now % tz = hour + minute * m2h

  end function now


  pure elemental integer function weekday(self)
    ! Returns the day of the week calculated using Zeller's congruence.
    ! Returned value is an integer scalar in the range [0-6], such that:
    !
    ! 0: Sunday
    ! 1: Monday
    ! 2: Tuesday
    ! 3: Wednesday
    ! 4: Thursday
    ! 5: Friday
    ! 6: Saturday
    class(datetime), intent(in) :: self
    integer :: year, month, j, k

    year  = self % year
    month = self % month

    if (month <= 2) then
      month = month + 12
      year  = year - 1
    end if

    j = year / 100
    k = mod(year, 100)

    weekday = mod(self % day + ((month + 1) * 26) / 10 + k + k / 4 + j / 4 + 5 * j, 7) -1

    if (weekday < 0) weekday = 6

  end function weekday


  pure elemental integer function isoweekday(self)
    ! Returns the day of the week per ISO 8601 returned from weekday().
    ! Returned value is an integer scalar in the range [1-7].
    class(datetime), intent(in) :: self
    isoweekday = self % weekday()
    if (isoweekday == 0) isoweekday = 7
  end function isoweekday


  pure elemental character(9) function weekdayLong(self)
    ! Returns the full name of the day of the week.
    class(datetime), intent(in) :: self
    character(9), parameter :: &
      days(*) = ['Sunday   ', 'Monday   ', 'Tuesday  ','Wednesday', &
                 'Thursday ', 'Friday   ', 'Saturday ']
    weekdayLong = days(self % weekday() + 1)
  end function weekdayLong


  pure elemental character(9) function isoweekdayLong(self)
    ! Returns the full name of the day of the week for ISO 8601
    ! ordered weekdays.
    class(datetime), intent(in) :: self
    character(9), parameter :: &
      days(7) = ['Monday   ','Tuesday  ','Wednesday','Thursday ', &
                 'Friday   ','Saturday ','Sunday   ']
    isoweekdayLong = days(self % isoweekday())
  end function isoweekdayLong


  pure elemental character(3) function weekdayShort(self)
    ! Returns the short (3-letter) name of the day of the week.
    class(datetime), intent(in) :: self
    character(3), parameter :: days(7) = ['Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat']
    weekdayShort = days(self % weekday() + 1)
  end function weekdayShort


  pure elemental character(3) function isoweekdayShort(self)
    ! Returns the short (3-letter) name of the day of the week
    ! based on ISO 8601 ordering.
    class(datetime), intent(in) :: self
    character(3), parameter :: days(7) = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
    isoweekdayShort = days(self % isoweekday())
  end function isoweekdayShort


  function isocalendar(self)
    ! Returns an array of 3 integers, year, week number, and week day,
    ! as defined by ISO 8601 week date. Essentially a wrapper around C
    ! `strftime` function.
    class(datetime), intent(in) :: self
    integer :: isocalendar(3)
    integer :: year, week, wday, rc
    character(20) :: string

    rc = c_strftime(string, len(string), '%G %V %u' // c_null_char, self % tm())

    read(string(1:4), '(i4)') year
    read(string(6:7), '(i2)') week
    read(string(9:9), '(i1)') wday

    isocalendar = [year, week, wday]

  end function isocalendar


  integer(int64) function secondsSinceEpoch(self)
    ! Returns an integer number of seconds since the UNIX Epoch (1 Jan 1970).
    ! Since Windows does not have strftime('%s'), we implement this using
    ! datetime itself.
    class(datetime), intent(in) :: self
    type(timedelta) :: delta
    type(datetime) :: this_time, unix_time

    this_time = datetime(self%year, self%month, self%day, &
                         self%hour, self%minute, self%second)
    unix_time = datetime(1970, 1, 1, 0, 0, 0)
    delta = this_time - unix_time
    secondsSinceEpoch = delta%total_seconds()

  end function secondsSinceEpoch


  function strftime(self, format)
    ! Wrapper around C and C++ `strftime` function.
    class(datetime), intent(in) :: self
    character(*), intent(in)  :: format
    character(:), allocatable :: strftime
    integer :: rc
    character(MAXSTRLEN) :: resultString
    resultString = ""
    rc = c_strftime(resultString, len(resultString), trim(format) // c_null_char, &
                    self % tm())
    if (rc == 0) write(stderr,*) "datetime:strftime failure, format: ", trim(format)
    strftime = resultString(1:len_trim(resultString)-1)  !< strip null
  end function strftime


  pure elemental type(tm_struct) function tm(self)
    ! Returns a `tm_struct` instance of the current `datetime`.
    class(datetime), intent(in) :: self
    tm % tm_sec = self % second
    tm % tm_min = self % minute
    tm % tm_hour = self % hour
    tm % tm_mday = self % day
    tm % tm_mon = self % month - 1
    tm % tm_year = self % year - 1900
    tm % tm_wday = self % weekday()
    tm % tm_yday = self % yearday() - 1
    tm % tm_isdst = -1
  end function tm


  pure elemental character(5) function tzOffset(self)
    ! Returns a character string with timezone offset in hours from UTC,
    ! in format +/-[hh][mm].
    class(datetime), intent(in) :: self
    integer :: hours,minutes

    if (self % tz < 0) then
      tzOffset(1:1) = '-'
    else
      tzOffset(1:1) = '+'
    end if

    hours = int(abs(self % tz))
    minutes = nint((abs(self % tz) - hours) * 60)

    if (minutes == 60) then
      minutes = 0
      hours = hours + 1
    end if

    write(tzOffset(2:5), '(2i2.2)') hours, minutes

  end function tzOffset


  pure elemental type(datetime) function utc(self)
    ! Returns the `datetime` instance at Coordinated Universal Time (UTC).
    class(datetime), intent(in) :: self
    integer :: hours, minutes, sgn
    hours = int(abs(self % tz))
    minutes = nint((abs(self % tz) - hours) * 60)
    sgn = int(sign(one, self % tz))
    utc = self - timedelta(hours=sgn * hours, minutes=sgn * minutes)
    utc % tz = 0
  end function utc


  pure elemental integer function yearday(self)
    ! Returns the integer day of the year (ordinal date).
    class(datetime), intent(in) :: self
    integer :: month
    yearday = 0
    do month = 1, self % month-1
      yearday = yearday + daysInMonth(month, self % year)
    end do
    yearday = yearday + self % day
  end function yearday


  pure elemental function datetime_plus_timedelta(d0,t) result(d)
    ! Adds a `timedelta` instance to a `datetime` instance, and returns a
    ! new `datetime` instance. Overloads the operator `+`.
    class(datetime), intent(in) :: d0
    class(timedelta), intent(in) :: t
    type(datetime) :: d

    integer :: milliseconds, seconds, minutes, hours, days

    d = datetime(year = d0 % getYear(),               &
                 month = d0 % getMonth(),             &
                 day = d0 % getDay(),                 &
                 hour = d0 % getHour(),               &
                 minute = d0 % getMinute(),           &
                 second = d0 % getSecond(),           &
                 millisecond = d0 % getMillisecond(), &
                 tz = d0 % getTz())

    milliseconds = t % getMilliseconds()
    seconds = t % getSeconds()
    minutes = t % getMinutes()
    hours = t % getHours()
    days = t % getDays()

    if (milliseconds /= 0) call d % addMilliseconds(milliseconds)
    if (seconds /= 0) call d % addSeconds(seconds)
    if (minutes /= 0) call d % addMinutes(minutes)
    if (hours /= 0) call d % addHours(hours)
    if (days /= 0) call d % addDays(days)

  end function datetime_plus_timedelta


  pure elemental function timedelta_plus_datetime(t,d0) result(d)
    ! Adds a `timedelta` instance to a `datetime` instance, and returns a
    ! new `datetime` instance. Overloads the operator `+`.
    class(timedelta), intent(in) :: t
    class(datetime), intent(in) :: d0
    type(datetime) :: d
    d = d0 + t
  end function timedelta_plus_datetime


  pure elemental function datetime_minus_timedelta(d0,t) result(d)
    ! Subtracts a `timedelta` instance from a `datetime` instance and
    ! returns a new `datetime` instance. Overloads the operator `-`.
    class(datetime), intent(in) :: d0
    class(timedelta), intent(in) :: t
    type(datetime) :: d
    d = d0 + (-t)
  end function datetime_minus_timedelta


  pure elemental function datetime_minus_datetime(d0,d1) result(t)
    ! Subtracts a `datetime` instance from another `datetime` instance,
    ! and returns a `timedelta` instance. Overloads the operator `-`.
    class(datetime), intent(in) :: d0, d1
    type(timedelta) :: t
    real(real64) :: daysDiff
    integer :: days,hours,minutes,seconds,milliseconds
    integer :: sign_

    daysDiff = date2num(d0)-date2num(d1)

    if (daysDiff < 0) then
      sign_ = -1
      daysDiff = ABS(daysDiff)
    else
      sign_ = 1
    end if

    days         = int(daysDiff)
    hours        = int((daysDiff-days)*d2h)
    minutes      = int((daysDiff-days-hours*h2d)*d2m)
    seconds      = int((daysDiff-days-hours*h2d-minutes*m2d)*d2s)
    milliseconds = nint((daysDiff-days-hours*h2d-minutes*m2d&
                                 -seconds*s2d)*d2s*1e3_real64)

    t = timedelta(sign_*days,sign_*hours,sign_*minutes,sign_*seconds,&
                  sign_*milliseconds)

  end function datetime_minus_datetime


  pure elemental logical function datetime_gt(d0,d1) result(res)
    ! `datetime` comparison operator that returns `.true.` if `d0` is
    ! greater than `d1` and `.false.` otherwise. Overloads the
    ! operator `>`.
    class(datetime), intent(in) :: d0, d1
    type(datetime) :: d0_utc, d1_utc

    ! Convert to UTC before making comparison
    d0_utc = d0 % utc()
    d1_utc = d1 % utc()

    ! Compare years
    if (d0_utc % year > d1_utc % year) then
      res = .true.
    else if (d0_utc % year < d1_utc % year) then
      res = .false.
    else

      ! Compare months
      if (d0_utc % month > d1_utc % month) then
        res = .true.
      else if (d0_utc % month < d1_utc % month) then
        res = .false.
      else

        ! Compare days
        if (d0_utc % day > d1_utc % day) then
          res = .true.
        else if (d0_utc % day < d1_utc % day) then
          res = .false.
        else

          ! Compare hours
          if (d0_utc % hour > d1_utc % hour) then
            res = .true.
          else if (d0_utc % hour < d1_utc % hour) then
            res = .false.
          else

            ! Compare minutes
            if (d0_utc % minute > d1_utc % minute) then
              res = .true.
            else if (d0_utc % minute < d1_utc % minute) then
              res = .false.
            else

              ! Compare seconds
              if (d0_utc % second > d1_utc % second) then
                res = .true.
              else if (d0_utc % second < d1_utc % second) then
                res = .false.
              else

                ! Compare milliseconds
                if (d0_utc % millisecond > d1_utc % millisecond) then
                  res = .true.
                else
                  res = .false.
                end if

              end if
            end if
          end if
        end if
      end if
    end if

  end function datetime_gt


  pure elemental logical function datetime_lt(d0,d1) result(res)
    ! `datetime` comparison operator that returns `.true.` if `d0` is
    ! less than `d1` and `.false.` otherwise. Overloads the operator `<`.
    class(datetime), intent(in) :: d0, d1
    res = d1 > d0
  end function datetime_lt


  pure elemental logical function datetime_eq(d0,d1) result(res)
    ! `datetime` comparison operator that returns `.true.` if `d0` is
    ! equal to `d1` and `.false.` otherwise. Overloads the operator `==`.
    class(datetime), intent(in) :: d0, d1
    type(datetime) :: d0_utc, d1_utc

    ! Convert to UTC before making comparison
    d0_utc = d0 % utc()
    d1_utc = d1 % utc()

    res = d0_utc % year        == d1_utc % year   .and. &
          d0_utc % month       == d1_utc % month  .and. &
          d0_utc % day         == d1_utc % day    .and. &
          d0_utc % hour        == d1_utc % hour   .and. &
          d0_utc % minute      == d1_utc % minute .and. &
          d0_utc % second      == d1_utc % second .and. &
          d0_utc % millisecond == d1_utc % millisecond

  end function datetime_eq


  pure elemental logical function datetime_neq(d0,d1) result(res)
    ! `datetime` comparison operator that eturns `.true.` if `d0` is
    ! not equal to `d1` and `.false.` otherwise. Overloads the operator `/=`.
    class(datetime), intent(in) :: d0, d1
    res = .not. d0 == d1
  end function datetime_neq


  pure elemental logical function datetime_ge(d0,d1) result(res)
    ! `datetime` comparison operator. Returns `.true.` if `d0` is greater
    ! than or equal to `d1` and `.false.` otherwise. Overloads the
    ! operator `>=`.
    class(datetime), intent(in) :: d0, d1
    res = d0 > d1 .or. d0 == d1
  end function datetime_ge


  pure elemental logical function datetime_le(d0,d1) result(res)
    ! `datetime` comparison operator. Returns `.true.` if `d0` is less
    ! than or equal to `d1`, and `.false.` otherwise. Overloads the
    ! operator `<=`.
    class(datetime), intent(in) :: d0, d1
    res = d1 > d0 .or. d0 == d1
  end function datetime_le


  pure elemental logical function isLeapYear(year)
    ! Returns `.true.` if year is leap year and `.false.` otherwise.
    integer, intent(in) :: year
    isLeapYear = (mod(year,4) == 0 .and. .not. mod(year,100) == 0)&
            .or. (mod(year,400) == 0)
  end function isLeapYear


  pure function datetimeRange(d0, d1, t)
    ! Given start and end `datetime` instances `d0` and `d1` and time
    ! increment as `timedelta` instance `t`, returns an array of
    ! `datetime` instances. The number of elements is the number of whole
    ! time increments contained between datetimes `d0` and `d1`.
    type(datetime), intent(in) :: d0, d1
    type(timedelta), intent(in) :: t
    real(real64) :: datenum0, datenum1, eps, increment
    type(datetime), allocatable :: datetimeRange(:)
    integer :: n, nm
    eps = 1e-10_real64
    datenum0 = date2num(d0)
    datenum1 = date2num(d1)
    increment = t % total_seconds() * s2d
    nm = floor((datenum1 - datenum0 + eps) / increment) + 1
    allocate(datetimeRange(nm))
    do n = 1, nm
      datetimeRange(n) = num2date(datenum0 + (n - 1) * increment)
    end do
  end function datetimeRange


  pure elemental integer function daysInMonth(month,year)
    ! Given integer month and year, returns an integer number
    ! of days in that particular month.
    integer, intent(in) :: month, year

    integer, parameter :: days(*) = [31, 28, 31, 30, 31, 30, &
                                     31, 31, 30, 31, 30, 31]

    if (month < 1 .or. month > 12) then
      ! Should raise an error and abort here, however we want to keep
      ! the pure and elemental attributes. Make sure this function is
      ! called with the month argument in range.
      daysInMonth = 0
      return
    end if

    if (month == 2 .and. isLeapYear(year)) then
      daysInMonth = 29
    else
      daysInMonth = days(month)
    end if

  end function daysInMonth


  pure elemental integer function daysInYear(year)
    ! Returns the number of days in year.
    integer, intent(in) :: year
    if (isLeapYear(year)) then
      daysInYear = 366
    else
      daysInYear = 365
    end if
  end function daysInYear


  pure elemental real(real64) function date2num(d)
    ! Given a datetime instance d, returns number of days since
    ! `0001-01-01 00:00:00`, taking into account the timezone offset.
    type(datetime), intent(in) :: d
    type(datetime) :: d_utc
    integer :: year

    ! Convert to UTC first
    d_utc = d % utc()

    ! d_utc % year must be positive:
    if (d_utc % year < 1) then
      date2num = 0
      return
    end if

    date2num = 0
    do year = 1,d_utc % year-1
      date2num = date2num + daysInYear(year)
    end do

    date2num = date2num          &
             + d_utc % yearday() &
             + d_utc % hour*h2d  &
             + d_utc % minute*m2d&
             + (d_utc % second+1e-3_real64*d_utc % millisecond)*s2d

  end function date2num


  pure elemental type(datetime) function num2date(num)
    ! Given number of days since `0001-01-01 00:00:00`, returns a
    ! correspoding `datetime` instance.
    real(real64), intent(in) :: num
    integer :: year, month, day, hour, minute, second, millisecond
    real(real64) :: days, totseconds

    ! num must be positive
    if (num < 0) then
      num2date = datetime(1)
      return
    end if

    days = num

    year = 1
    do
      if (int(days) <= daysInYear(year))exit
      days = days-daysInYear(year)
      year = year+1
    end do

    month = 1
    do
      if (inT(days) <= daysInMonth(month,year))exit
      days = days-daysInMonth(month,year)
      month = month+1
    end do

    day         = int(days)
    totseconds  = (days-day)*d2s
    hour        = int(totseconds*s2h)
    minute      = int((totseconds-hour*h2s)*s2m)
    second      = int(totseconds-hour*h2s-minute*m2s)
    millisecond = nint((totseconds-int(totseconds))*1e3_real64)

    num2date = datetime(year,month,day,hour,minute,second,millisecond,tz=zero)

    ! Handle a special case caused by floating-point arithmethic:
    if (num2date % millisecond == 1000) then
      num2date % millisecond = 0
      call num2date % addSeconds(1)
    end if

    if (num2date % second == 60) then
      num2date % second = 0
      call num2date % addMinutes(1)
    end if
    if (num2date % minute == 60) then
      num2date % minute = 0
      call num2date % addHours(1)
    end if
    if (num2date % hour == 24) then
      num2date % hour = 0
      call num2date % addDays(1)
    end if

  end function num2date


  type(datetime) function strptime(str,format)
    ! A wrapper function around C/C++ strptime function.
    ! Returns a `datetime` instance.
    character(*), intent(in) :: str, format
    integer :: rc
    type(tm_struct) :: tm
    rc = c_strptime(trim(str) // c_null_char, trim(format) // c_null_char, tm)
    strptime = tm2date(tm)
  end function strptime


  pure elemental type(datetime) function tm2date(ctime)
    ! Given a `tm_struct` instance, returns a corresponding `datetime`
    ! instance.
    type(tm_struct), intent(in) :: ctime

    tm2date % millisecond = 0
    tm2date % second      = ctime % tm_sec
    tm2date % minute      = ctime % tm_min
    tm2date % hour        = ctime % tm_hour
    tm2date % day         = ctime % tm_mday
    tm2date % month       = ctime % tm_mon+1
    tm2date % year        = ctime % tm_year+1900
    tm2date % tz          = 0

  end function tm2date


  pure function int2str(i, length)
    ! Converts an integer `i` into a character string of requested length,
    ! pre-pending zeros if necessary.
    integer, intent(in) :: i, length
    character(length) :: int2str
    character(2) :: string
    write(string, '(i2)') length
    write(int2str, '(i' // string // '.' // string //')') i
  end function int2str


  pure elemental type(timedelta) function timedelta_constructor( &
    days, hours, minutes, seconds, milliseconds)
    ! Constructor function for the `timedelta` class.
    integer, intent(in), optional :: days, hours, minutes, seconds, milliseconds

    timedelta_constructor % days = 0
    if (present(days)) timedelta_constructor % days = days

    timedelta_constructor % hours = 0
    if (present(hours)) timedelta_constructor % hours = hours

    timedelta_constructor % minutes = 0
    if (present(minutes)) timedelta_constructor % minutes = minutes

    timedelta_constructor % seconds = 0
    if (present(seconds)) timedelta_constructor % seconds = seconds

    timedelta_constructor % milliseconds = 0
    if (present(milliseconds)) timedelta_constructor % milliseconds = milliseconds

  end function timedelta_constructor


  pure elemental integer function getDays(self)
    ! Returns the number of days.
    class(timedelta), intent(in) :: self
    getDays = self % days
  end function getDays


  pure elemental integer function getHours(self)
    ! Returns the number of hours.
    class(timedelta), intent(in) :: self
    getHours = self % hours
  end function getHours


  pure elemental integer function getMinutes(self)
    ! Returns the number of minutes.
    class(timedelta), intent(in) :: self
    getMinutes = self % minutes
  end function getMinutes


  pure elemental integer function getSeconds(self)
    ! Returns the number of seconds.
    class(timedelta), intent(in) :: self
    getSeconds = self % seconds
  end function getSeconds


  pure elemental integer function getMilliseconds(self)
    ! Returns the number of milliseconds.
    class(timedelta), intent(in) :: self
    getMilliseconds = self % milliseconds
  end function getMilliseconds


  pure elemental real(real64) function total_seconds(self)
    ! Returns a total number of seconds contained in a `timedelta`
    ! instance.
    class(timedelta), intent(in) :: self
    total_seconds = self % days*86400._real64 &
                  + self % hours*3600._real64 &
                  + self % minutes*60._real64 &
                  + self % seconds            &
                  + self % milliseconds*1e-3_real64
  end function total_seconds


  pure elemental type(timedelta) function timedelta_plus_timedelta(t0,t1) result(t)
    ! Adds two `timedelta` instances together and returns a `timedelta`
    ! instance. Overloads the operator `+`.
    class(timedelta), intent(in) :: t0, t1
    t = timedelta(days         = t0 % days         + t1 % days,    &
                  hours        = t0 % hours        + t1 % hours,   &
                  minutes      = t0 % minutes      + t1 % minutes, &
                  seconds      = t0 % seconds      + t1 % seconds, &
                  milliseconds = t0 % milliseconds + t1 % milliseconds)
  end function timedelta_plus_timedelta


  pure elemental type(timedelta) function timedelta_minus_timedelta(t0,t1) result(t)
    ! Subtracts a `timedelta` instance from another. Returns a
    ! `timedelta` instance. Overloads the operator `-`.
    class(timedelta), intent(in) :: t0, t1
    t = t0 + (-t1)
  end function timedelta_minus_timedelta


  pure elemental type(timedelta) function unary_minus_timedelta(t0) result(t)
    ! Takes a negative of a `timedelta` instance. Overloads the operator `-`.
    class(timedelta), intent(in) :: t0
    t % days         = -t0 % days
    t % hours        = -t0 % hours
    t % minutes      = -t0 % minutes
    t % seconds      = -t0 % seconds
    t % milliseconds = -t0 % milliseconds
  end function unary_minus_timedelta


  pure elemental logical function timedelta_eq(td0,td1) result(res)
    ! `timedelta` object comparison operator. Returns `.true.` if `td0`
    ! is equal to `td1` and `.false.` otherwise. Overloads the operator
    ! `==`.
    class(timedelta), intent(in) :: td0, td1
    res = td0 % total_seconds() == td1 % total_seconds()
  end function timedelta_eq


  pure elemental logical function timedelta_neq(td0,td1) result(res)
    ! `timedelta` object comparison operator. Returns `.true.` if `td0`
    ! is not equal to `td1` and `.false.` otherwise. Overloads the
    ! operator `/=`.
    class(timedelta), intent(in) :: td0, td1
    res = td0 % total_seconds() /= td1 % total_seconds()
  end function timedelta_neq


  pure elemental logical function timedelta_gt(td0,td1) result(res)
    ! `timedelta` object comparison operator. Returns `.true.` if
    ! `td0` is greater than `td1` and `.false.` otherwise. Overloads the
    ! operator `>`.
    class(timedelta), intent(in) :: td0, td1
    res = td0 % total_seconds() > td1 % total_seconds()
  end function timedelta_gt


  pure elemental logical function timedelta_ge(td0,td1) result(res)
    ! `timedelta` object comparison operator. Returns `.true.` if `td0`
    ! is greater than or equal to `td1` and `.false.` otherwise.
    ! Overloads the operator >=.
    class(timedelta), intent(in) :: td0, td1
    res = td0 % total_seconds() >= td1 % total_seconds()
  end function timedelta_ge


  pure elemental logical function timedelta_lt(td0,td1) result(res)
    ! `timedelta` object comparison operator. Returns `.true.` if `td0`
    ! is less than `td1` and `.false.` otherwise. Overloads the operator
    ! `<`.
    class(timedelta), intent(in) :: td0, td1
    res = td0 % total_seconds() < td1 % total_seconds()
  end function timedelta_lt


  pure elemental logical function timedelta_le(td0,td1) result(res)
    ! `timedelta` object comparison operator. Returns `.true.` if `td0`
    ! is less than or equal to `td1` and `.false.` otherwise. Overloads
    ! the operator `<=`.
    class(timedelta), intent(in) :: td0, td1
    res = td0 % total_seconds() <= td1 % total_seconds()
  end function timedelta_le

end module datetime_module