module integral_data

  use typedef, only: &
       real64, int64

  implicit none

  save

  integer(kind=int64), parameter :: &
       nsavedefault = 9223372036854775807_int64, &
       nsizedefault = 1024_int64

  integer(kind=int64) :: &
       ndata = -1, &
       ndim = -1, &
       nsize = -1, &
       nsave = -1


  real(kind=real64), dimension(:), allocatable, target :: &
       coord
  real(kind=real64), dimension(:,:), allocatable, target :: &
       data

contains

  subroutine clear()

    if (allocated(coord)) deallocate(coord)
    if (allocated(data)) deallocate(data)
    ndata = -1
    ndim = -1
    nsize = -1
    nsave = -1

  end subroutine clear

  subroutine get_ndata(ndata_, ndim_, mstep_)

    use typedef, only: &
         int64

    implicit none

    integer(kind=int64), intent(out) :: &
         ndata_, ndim_, mstep_

    mstep_ = ndata
    ndata_ = min(ndata, nsize)
    ndim_ = ndim

  end subroutine get_ndata

  subroutine get_data(coord_, data_)

    use typedef, only: &
         int64, real64

    implicit none

    real(kind=real64), intent(out), dimension(:) :: &
         coord_
    real(kind=real64), intent(out), dimension(:,:) :: &
         data_

    ! local variables

    integer(kind=int64) :: &
         m1, m2, n1, k

    m1 = size(coord_, 1)
    m2 = size(data_, 2)
    n1 = size(data_, 1)

    k = min(nsave, ndata)

    if (m1 < k) error stop 'coord too small'
    if (m2 < k) error stop 'data dim 2 too small'
    if (n1 < ndim) error stop 'data dim 1 too small'

    if (ndata > nsize) then
       k = mod(ndata, nsize)
       coord_(1:nsize) = cshift(coord(1:nsize),shift=k)
       data_(1:ndim,1:nsize) = cshift(data(1:ndim,1:nsize),shift=k,dim=2)
    else
       coord_(1:ndata) = coord(1:ndata)
       data_(1:ndim,1:ndata) = data(1:ndim,1:ndata)
    endif

  end subroutine get_data


  subroutine initialize(ndim_, nsave_, nsize_)

    use typedef, only: &
         int64

    implicit none

    integer(kind=int64), intent(in) :: &
         ndim_

    integer(kind=int64), intent(in), optional :: &
         nsize_, nsave_

    call clear()

    ndim = ndim_
    if (present(nsize_)) then
       nsize = nsize_
    else
       nsize = nsizedefault
    endif
    if (present(nsave_)) then
       nsave = nsave_
    else
       nsave = nsavedefault
    endif
    nsize = min(nsize, nsave)
    ndata = 0

    allocate(coord(nsize), data(ndim, nsize))

  end subroutine initialize


  subroutine grow(increase_)

    use typedef, only: &
         int64, real64

    implicit none

    real(kind=real64), parameter :: &
         GOLDEN = 0.5d0 * (sqrt(5.d0) - 1.d0)

    integer(kind=int64), intent(in), optional :: &
         increase_

    integer(kind=int64) :: &
         increase, newsize

    real(kind=real64), dimension(:), allocatable :: &
         coord_temp
    real(kind=real64), dimension(:,:), allocatable :: &
         data_temp

    if (nsize == 0) error stop '[grow] ERROR: not initialized'

    if (present(increase_)) then
       increase = increase_
    else
       increase = max(1024, int(nsize * GOLDEN))
    endif

    newsize = min(nsize + increase, nsave)

    ALLOCATE(coord_temp(newsize), data_temp(ndim, newsize))
    coord_temp(1:nsize) = coord(1:nsize)
    data_temp(:,1:nsize) = data(:,1:nsize)
    call MOVE_ALLOC(coord_temp, coord)
    call MOVE_ALLOC(data_temp, data)

    nsize = newsize

  end subroutine grow


  subroutine append(t, y)

    use typedef, only: &
         int64, real64

    implicit none

    real(kind=real64), intent(in) :: &
         t
    real(kind=real64), intent(in), dimension(:) :: &
         y
    integer(kind=int64) :: &
         k

    if (nsize <= 0) error stop '[append] ERROR: not initialized'

    if (size(y, 1) /= ndim) error stop '[append] ERROR: dimension mismatch'

    if ((ndata < nsave) .and. (ndata == nsize)) call grow

    k = mod(ndata, nsave) + 1
    ndata = ndata + 1
    coord(k) = t
    data(:,k) = y(:)

  end subroutine append


end module integral_data