module nbi_pllbuf

  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: zero
  implicit none
  private :: fp, zero
  save

  integer :: nbuf_dep = 0
  integer :: nbuf_orb = 0
  integer :: nbuf_ntr = 0
  integer :: nbuf_lst = 0
  integer :: nbuf_data = 0
  integer :: nmax_data = 0
  integer :: n_data = 0
  real(fp), dimension(:), allocatable :: r8buf,r8buf1

  private :: r8buf_expand

  interface r8put
     module procedure r8put_0d,r8put_1d,r8put_2d,r8put_3d,r8put_4d
  end interface

  interface r8get
     module procedure r8get_0d,r8get_1d,r8get_2d,r8get_3d,r8get_4d
  end interface

  interface r8nz
     module procedure r8nz_0d,r8nz_1d,r8nz_2d,r8nz_3d,r8nz_4d
  end interface

  interface r8zero
     module procedure r8zero_0d,r8zero_1d,r8zero_2d,r8zero_3d,r8zero_4d
  end interface

  contains
    subroutine r8buf_expand(ineed)
      integer, intent(in) :: ineed  ! minimum allocation

      real(fp), dimension(:), allocatable :: r8buf_tmp
      integer :: inew

      if(n_data.gt.0) then
         allocate(r8buf_tmp(n_data)); r8buf_tmp = r8buf(1:n_data)
      endif
      if(allocated(r8buf)) deallocate(r8buf)
      if(allocated(r8buf1)) deallocate(r8buf1)

      inew=max(nbuf_data,ineed,2*nmax_data)
      nmax_data=inew
      allocate(r8buf(nmax_data),r8buf1(nmax_data))

      if(n_data.gt.0) then
         r8buf(1:n_data)=r8buf_tmp
         deallocate(r8buf_tmp)
      endif
    end subroutine r8buf_expand

    subroutine r8put_0d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp) :: zdata

      integer :: ineed

      ineed=n_data+1
      if(ineed.gt.nmax_data) call r8buf_expand(ineed)

      iadr=iadr+1
      r8buf(iadr)=zdata
      n_data=iadr
    end subroutine r8put_0d

    subroutine r8get_0d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp) :: zdata

      iadr=iadr+1
      zdata=r8buf1(iadr)
    end subroutine r8get_0d

    subroutine r8put_1d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp), dimension(:) :: zdata

      integer :: isize,ineed,i1

      isize=size(zdata)
      ineed=n_data+isize
      if(ineed.gt.nmax_data) call r8buf_expand(ineed)

      do i1=1,isize
         iadr=iadr+1
         r8buf(iadr)=zdata(i1)
      enddo
      n_data=iadr
    end subroutine r8put_1d

    subroutine r8get_1d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp), dimension(:) :: zdata

      integer :: isize,i1

      isize=size(zdata)

      do i1=1,isize
         iadr=iadr+1
         zdata(i1)=r8buf1(iadr)
      enddo
    end subroutine r8get_1d

    subroutine r8put_2d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp), dimension(:,:) :: zdata

      integer :: isize,ineed,i1,i2,isizes(2)

      isize=1
      do i1=1,2
         isizes(i1)=size(zdata,i1)
         isize=isize*isizes(i1)
      enddo
      ineed=n_data+isize
      if(ineed.gt.nmax_data) call r8buf_expand(ineed)

      do i2=1,isizes(2)
         do i1=1,isizes(1)
            iadr=iadr+1
            r8buf(iadr)=zdata(i1,i2)
         enddo
      enddo
      n_data=iadr
    end subroutine r8put_2d

    subroutine r8get_2d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp), dimension(:,:) :: zdata

      integer :: i1,i2,isizes(2)

      do i1=1,2
         isizes(i1)=size(zdata,i1)
      enddo

      do i2=1,isizes(2)
         do i1=1,isizes(1)
            iadr=iadr+1
            zdata(i1,i2)=r8buf1(iadr)
         enddo
      enddo
    end subroutine r8get_2d

    subroutine r8put_3d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp), dimension(:,:,:) :: zdata

      integer :: isize,ineed,i1,i2,i3,isizes(3)

      isize=1
      do i1=1,3
         isizes(i1)=size(zdata,i1)
         isize=isize*isizes(i1)
      enddo
      ineed=n_data+isize
      if(ineed.gt.nmax_data) call r8buf_expand(ineed)

      do i3=1,isizes(3)
         do i2=1,isizes(2)
            do i1=1,isizes(1)
               iadr=iadr+1
               r8buf(iadr)=zdata(i1,i2,i3)
            enddo
         enddo
      enddo
      n_data=iadr
    end subroutine r8put_3d

    subroutine r8get_3d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp), dimension(:,:,:) :: zdata

      integer :: i1,i2,i3,isizes(3)

      do i1=1,3
         isizes(i1)=size(zdata,i1)
      enddo

      do i3=1,isizes(3)
         do i2=1,isizes(2)
            do i1=1,isizes(1)
               iadr=iadr+1
               zdata(i1,i2,i3)=r8buf1(iadr)
            enddo
         enddo
      enddo
    end subroutine r8get_3d

    subroutine r8put_4d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp), dimension(:,:,:,:) :: zdata

      integer :: isize,ineed,i1,i2,i3,i4,isizes(4)

      isize=1
      do i1=1,4
         isizes(i1)=size(zdata,i1)
         isize=isize*isizes(i1)
      enddo
      ineed=n_data+isize
      if(ineed.gt.nmax_data) call r8buf_expand(ineed)

      do i4=1,isizes(4)
         do i3=1,isizes(3)
            do i2=1,isizes(2)
               do i1=1,isizes(1)
                  iadr=iadr+1
                  r8buf(iadr)=zdata(i1,i2,i3,i4)
               enddo
            enddo
         enddo
      enddo
      n_data=iadr
    end subroutine r8put_4d

    subroutine r8get_4d(zdata,iadr)
      integer, intent(inout) :: iadr  ! address last used
      real(fp), dimension(:,:,:,:) :: zdata

      integer :: i1,i2,i3,i4,isizes(4)

      do i1=1,4
         isizes(i1)=size(zdata,i1)
      enddo

      do i4=1,isizes(4)
         do i3=1,isizes(3)
            do i2=1,isizes(2)
               do i1=1,isizes(1)
                  iadr=iadr+1
                  zdata(i1,i2,i3,i4)=r8buf1(iadr)
               enddo
            enddo
         enddo
      enddo
    end subroutine r8get_4d

    subroutine r8nz_0d(zdata,inum)
      real(fp), intent(in) :: zdata
      integer, intent(out) :: inum  ! number of non-zero values found

      nbuf_data=nbuf_data+1
      inum=0
      if(zdata.ne.zero) inum=1
    end subroutine r8nz_0d

    subroutine r8nz_1d(zdata,inum)
      real(fp), dimension(:), intent(in) :: zdata
      integer, intent(out) :: inum  ! number of non-zero values found

      integer :: isize,i1

      isize=size(zdata)
      nbuf_data=nbuf_data+isize

      inum=0
      do i1=1,isize
         if(zdata(i1).ne.zero) inum=inum+1
      enddo
    end subroutine r8nz_1d

    subroutine r8nz_2d(zdata,inum)
      real(fp), dimension(:,:), intent(in) :: zdata
      integer, intent(out) :: inum  ! number of non-zero values found

      integer :: i1,i2,isizes(2)

      do i1=1,2
         isizes(i1)=size(zdata,i1)
      enddo
      nbuf_data=nbuf_data+isizes(1)*isizes(2)

      inum=0
      do i2=1,isizes(2)
         do i1=1,isizes(1)
            if(zdata(i1,i2).ne.zero) inum=inum+1
         enddo
      enddo
    end subroutine r8nz_2d

    subroutine r8nz_3d(zdata,inum)
      real(fp), dimension(:,:,:), intent(in) :: zdata
      integer, intent(out) :: inum  ! number of non-zero values found

      integer :: i1,i2,i3,isizes(3)

      do i1=1,3
         isizes(i1)=size(zdata,i1)
      enddo
      nbuf_data=nbuf_data+isizes(1)*isizes(2)*isizes(3)

      inum=0
      do i3=1,isizes(3)
         do i2=1,isizes(2)
            do i1=1,isizes(1)
               if(zdata(i1,i2,i3).ne.zero) inum=inum+1
            enddo
         enddo
      enddo
    end subroutine r8nz_3d

    subroutine r8nz_4d(zdata,inum)
      real(fp), dimension(:,:,:,:), intent(in) :: zdata
      integer, intent(out) :: inum  ! number of non-zero values found

      integer :: i1,i2,i3,i4,isizes(4)

      do i1=1,4
         isizes(i1)=size(zdata,i1)
      enddo
      nbuf_data=nbuf_data+isizes(1)*isizes(2)*isizes(3)*isizes(4)

      inum=0
      do i4=1,isizes(4)
         do i3=1,isizes(3)
            do i2=1,isizes(2)
               do i1=1,isizes(1)
                  if(zdata(i1,i2,i3,i4).ne.zero) inum=inum+1
               enddo
            enddo
         enddo
      enddo
    end subroutine r8nz_4d

    subroutine r8zero_0d(zdata,inum)
      real(fp), intent(inout) :: zdata
      integer, intent(out) :: inum  ! dummy argument 

      inum=0
      zdata=zero

    end subroutine r8zero_0d

    subroutine r8zero_1d(zdata,inum)
      real(fp), dimension(:), intent(inout) :: zdata
      integer, intent(out) :: inum  ! dummy argument 

      integer :: isize,i1

      inum=0
      isize=size(zdata)

      do i1=1,isize
         zdata(i1)=zero
      enddo
    end subroutine r8zero_1d

    subroutine r8zero_2d(zdata,inum)
      real(fp), dimension(:,:), intent(inout) :: zdata
      integer, intent(out) :: inum  ! dummy argument 

      integer :: i1,i2,isizes(2)

      do i1=1,2
         isizes(i1)=size(zdata,i1)
      enddo

      inum=0
      do i2=1,isizes(2)
         do i1=1,isizes(1)
            zdata(i1,i2)=zero
         enddo
      enddo
    end subroutine r8zero_2d

    subroutine r8zero_3d(zdata,inum)
      real(fp), dimension(:,:,:), intent(inout) :: zdata
      integer, intent(out) :: inum  ! dummy argument 

      integer :: i1,i2,i3,isizes(3)

      do i1=1,3
         isizes(i1)=size(zdata,i1)
      enddo

      inum=0
      do i3=1,isizes(3)
         do i2=1,isizes(2)
            do i1=1,isizes(1)
               zdata(i1,i2,i3)=zero
            enddo
         enddo
      enddo
    end subroutine r8zero_3d

    subroutine r8zero_4d(zdata,inum)
      real(fp), dimension(:,:,:,:), intent(inout) :: zdata
      integer, intent(out) :: inum  ! dummy argument 

      integer :: i1,i2,i3,i4,isizes(4)

      do i1=1,4
         isizes(i1)=size(zdata,i1)
      enddo

      inum=0
      do i4=1,isizes(4)
         do i3=1,isizes(3)
            do i2=1,isizes(2)
               do i1=1,isizes(1)
                  zdata(i1,i2,i3,i4)=zero
               enddo
            enddo
         enddo
      enddo
    end subroutine r8zero_4d

end module nbi_pllbuf
