! part2d class for QuickPIC Open Source 1.0
! update: 04/18/2016

      module part2d_class

      use perrors_class
      use parallel_pipe_class
      use spect2d_class
      use fdist2d_class
      use ufield2d_class
      use part2d_lib
      use hdf5io_class
      use mpi
               
      implicit none

      private

      public :: part2d

      type part2d

         private

         class(spect2d), pointer, public :: sp => null()
         class(perrors), pointer, public :: err => null()
         class(parallel_pipe), pointer, public :: p => null()
!
! qbm = particle charge/mass ratio
! dt = time interval between successive calculations
! ci = reciprical of velocity of light
! xdim = dimension of the particle coordinates
! nbmax = size of buffer for passing particles between processors
! np = total number of particles
! npp = number of particles in current partition
! npmax = maximum number of particles in each partition
! part(:,:) = initial particle coordinates
! ppart(:,:,:) = particle coordinates for OpenMP
! nppmx, nppmx0, nbmaxp, ntmaxp, npbmx, irc, ncl, ihole, kpic = parameters for OpenMP
!         
         real :: qbm, dt, ci
         integer :: npmax, nbmax, np, xdim, npp = 0
         real, dimension(:,:), pointer :: part => null()
         real, dimension(:,:,:), pointer :: ppart => null()
         integer :: nppmx, nppmx0, nbmaxp, ntmaxp, npbmx, irc = 0
         integer, dimension(:,:), pointer :: ncl => null()
         integer, dimension(:,:,:), pointer :: ihole => null()
         integer, dimension(:), pointer :: kpic => null()
         
         contains
         
         generic :: new => init_part2d
         generic :: renew => renew_part2d
         generic :: del => end_part2d
         generic :: qdp => qdeposit
         generic :: amjdp => amjdeposit
         generic :: push => partpush
         generic :: pmv => pmove
         generic :: extpsi => extractpsi
         generic :: pcp => partcopy
         generic :: pcb => partcopyback
         generic :: psend => pipesend_part2d
         generic :: precv => piperecv_part2d
         generic :: wr => writehdf5_part2d
         procedure, private :: init_part2d, renew_part2d
         procedure, private :: end_part2d
         procedure, private :: qdeposit
         procedure, private :: amjdeposit
         procedure, private :: partpush
         procedure, private :: pmove
         procedure, private :: extractpsi
         procedure, private :: partcopy
         procedure, private :: partcopyback
         procedure, private :: pipesend_part2d
         procedure, private :: piperecv_part2d, writehdf5_part2d
                           
      end type 

      save      

      character(len=10) :: class = 'part2d:'
      character(len=128) :: erstr
! parameters for OpenMP
      integer :: mx = 16, my = 16
      integer :: mx1, myp1, mxyp1
! ppbuff = buffer array for reordering tiled particle array
      real, dimension(:,:,:), allocatable :: ppbuff
      integer :: szpbuf = 0
! sbufl/sbufr = particle buffers sent to nearby processors
! rbufl/rbufr = particle buffers received from nearby processors
      real, dimension(:,:), allocatable :: sbufl, sbufr, rbufl, rbufr
      integer :: szbufs = 0
! ncll/nclr/mcll/mclr = number offsets send/received from processors
      integer, dimension(:,:), allocatable :: ncll, nclr, mcll, mclr
      integer :: sznbufs = 0
      
      contains
!
      subroutine init_part2d(this,pp,perr,psp,pf,fd,qbm,dt,ci,xdim,s)
      
         implicit none
         
         class(part2d), intent(inout) :: this
         class(spect2d), intent(in), pointer :: psp
         class(perrors), intent(in), pointer :: perr
         class(parallel_pipe), intent(in), pointer :: pp
         class(fdist2d), intent(inout) :: pf
         class(ufield2d), intent(in), pointer :: fd
         real, intent(in) :: qbm, dt, ci, s
         integer, intent(in) :: xdim

! local data
         character(len=18), save :: sname = 'init_part2d:'
         integer :: xtras, noff, nxyp, nx, npmax, nbmax
                  
         this%sp => psp
         this%err => perr
         this%p => pp

         call this%err%werrfl2(class//sname//' started')

         this%qbm = qbm
         this%dt = dt
         this%ci = ci
         this%xdim = xdim
         npmax = pf%getnpmax()
         this%npmax = npmax
         nbmax = max(int(0.01*npmax),100)
         this%nbmax = nbmax
         noff = fd%getnoff()
         nxyp = fd%getnd2p()
         nx = fd%getnd1p()
         
         allocate(this%part(xdim,npmax))
         mx1 = (nx - 1)/mx + 1
         myp1 = (nxyp - 1)/my + 1; mxyp1 = mx1*myp1
         allocate(this%kpic(mxyp1))
         
         call pf%dist(this%part,this%npp,fd,s)

! find number of particles in each of mx, my tiles: updates kpic, nppmx
         call PPDBLKP2L(this%part,this%kpic,this%npp,noff,this%nppmx,&
         &this%xdim,npmax,mx,my,mx1,mxyp1,this%irc)
! check for errors
         if (this%irc /= 0) then
            write (erstr,*) 'PPDBLKP2L error, irc=', this%irc
            call this%err%equit(class//sname//erstr); return
         endif
         
         xtras = 10
!    
! allocate vector particle data
         this%nppmx0 = (1.0 + xtras)*this%nppmx
         this%ntmaxp = xtras*this%nppmx
         this%npbmx = xtras*this%nppmx
         this%nbmaxp = 0.25*mx1*this%npbmx
         allocate(this%ppart(xdim,this%nppmx0,mxyp1))
         allocate(this%ncl(8,mxyp1))
         allocate(this%ihole(2,this%ntmaxp+1,mxyp1))
!
! copy ordered particle data for OpenMP
         call PPPMOVIN2L(this%part,this%ppart,this%kpic,this%npp,noff,&
         &this%nppmx0,this%xdim,npmax,mx,my,mx1,mxyp1,this%irc)
! check for errors
         if (this%irc /= 0) then
            write (erstr,*) 'PPPMOVIN2L overflow error, irc=', this%irc
            call this%err%equit(class//sname//erstr); return
         endif
! 
! sanity check
!         call PPPCHECK2L(this%ppart,this%kpic,noff,nxyp,this%xdim,&
!         &this%nppmx0,nx,mx,my,mx1,myp1,this%irc)
! check error
!         if (this%irc /= 0) then
!            write (erstr,*) 'PPPCHECK2L error: irc=', this%irc
!            call this%err%equit(class//sname//erstr); return
!         endif

         call this%err%werrfl2(class//sname//' ended')

      end subroutine init_part2d
!
      subroutine end_part2d(this)
          
         implicit none
         
         class(part2d), intent(inout) :: this
         character(len=18), save :: sname = 'end_part2d:'

         call this%err%werrfl2(class//sname//' started')
         deallocate(this%part,this%ppart,this%ncl,this%ihole,this%kpic)
         call this%err%werrfl2(class//sname//' ended')
         
         return
         
      end subroutine end_part2d
!
      subroutine renew_part2d(this,pf,fd,s)
      
         implicit none
         
         class(part2d), intent(inout) :: this
         class(fdist2d), intent(inout) :: pf
         class(ufield2d), pointer, intent(in) :: fd
         real, intent(in) :: s

! local data
         character(len=18), save :: sname = 'renew_part2d:'
         integer :: noff, prof
                  
         call this%err%werrfl2(class//sname//' started')
         
         noff = fd%getnoff()
         prof = pf%getnpf()         
         
         call pf%dist(this%part,this%npp,fd,s)
         
         call PPDBLKP2L(this%part,this%kpic,this%npp,noff,this%nppmx,&
         &this%xdim,this%npmax,mx,my,mx1,mxyp1,this%irc)
! check for errors
         if (this%irc /= 0) then
            write (erstr,*) 'PPDBLKP2L error, irc=', this%irc
            call this%err%equit(class//sname//erstr); return
         endif
         
! copy ordered particle data for OpenMP
         call PPPMOVIN2L(this%part,this%ppart,this%kpic,this%npp,noff,&
         &this%nppmx0,this%xdim,this%npmax,mx,my,mx1,mxyp1,this%irc)
! check for errors
         if (this%irc /= 0) then
            write (erstr,*) 'PPPMOVIN2L overflow error, irc=', this%irc
            call this%err%equit(class//sname//erstr); return
         endif


         call this%err%werrfl2(class//sname//' ended')

      end subroutine renew_part2d
!      
      subroutine qdeposit(this,q)
! deposit the charge density      
      
         implicit none
         
         class(part2d), intent(in) :: this
         class(ufield2d), pointer, intent(in) :: q
! local data
         character(len=18), save :: sname = 'qdeposit:'
         real, dimension(:,:,:), pointer :: pq => null()
         integer :: noff, nxv, nypmx
                  
         call this%err%werrfl2(class//sname//' started')
         
         pq => q%getrf()
         noff = q%getnoff()
         nxv = size(pq,2)
         nypmx = size(pq,3)
         
         call PPGPPOST2L(this%ppart,pq(1,:,:),this%kpic,noff,&
         &this%xdim,this%nppmx0,mx,my,nxv,nypmx,mx1,mxyp1)         
         
         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine qdeposit
!      
      subroutine amjdeposit(this,ef,bf,psit,cu,amu,dcu,dex)
! deposit the current, acceleration and momentum flux      
      
         implicit none
         
         class(part2d), intent(inout) :: this
         class(ufield2d), pointer, intent(in) :: cu, amu, dcu
         class(ufield2d), pointer, intent(in) :: ef, bf, psit
         real, intent(in) :: dex
         character(len=18), save :: sname = 'amjdeposit'
! local data
         real, dimension(:,:,:), pointer :: pef => null(), pbf => null()
         real, dimension(:,:,:), pointer :: ppsit => null(), pcu => null()
         real, dimension(:,:,:), pointer :: pamu => null(), pdcu => null()
         integer :: noff, nyp, nx, nxv, nypmx

         call this%err%werrfl2(class//sname//' started')
         
         pef => ef%getrf(); pbf => bf%getrf()
         ppsit => psit%getrf(); pcu => cu%getrf()
         pamu => amu%getrf(); pdcu => dcu%getrf()
         noff = ef%getnoff()
         nxv = size(pef,2); nypmx = size(pef,3)
         nx = ef%getnd1(); nyp = ef%getnd2p()
         
         call PPGRDCJPPOST2L_QP(this%ppart,pef,pbf,ppsit(1,:,:),pcu,pdcu,&
         &pamu,this%kpic,noff,nyp,this%qbm, this%dt,this%ci,this%xdim,&
         &this%nppmx0,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex)
         
         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine amjdeposit
!      
      subroutine partpush(this,ef,bf,psit,dex)
      
         implicit none
         
         class(part2d), intent(inout) :: this
         class(ufield2d), pointer, intent(in) :: ef, bf, psit
         real, intent(in) :: dex
         character(len=18), save :: sname = 'partpush'
! local data
         real, dimension(:,:,:), pointer :: pef => null(), pbf => null()
         real, dimension(:,:,:), pointer :: ppsit => null()
         integer :: noff, nyp, nx, ny, nxv, nypmx
         real :: ek
         
         call this%err%werrfl2(class//sname//' started')

         pef => ef%getrf(); pbf => bf%getrf()
         ppsit => psit%getrf()
         noff = ef%getnoff(); ny = ef%getnd2()
         nxv = size(pef,2); nypmx = size(pef,3)
         nx = ef%getnd1(); nyp = ef%getnd2p()

         call PPGRBPPUSHF23L_QP(this%ppart,pef,pbf,ppsit(1,:,:),this%kpic,&
         &this%ncl,this%ihole,noff,nyp,this%qbm,this%dt,this%dt,this%ci,ek,&
         &this%xdim,this%nppmx0,nx,ny,mx,my,nxv,nypmx,mx1,mxyp1,this%ntmaxp,&
         &this%irc,dex)
         if (this%irc /= 0) then
            write (erstr,*) 'PPGRBPPUSHF23L_QP error, irc=', this%irc
            call this%err%equit(class//sname//erstr); return
         endif

         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine partpush
!
      subroutine pmove(this,fd)
      
         implicit none
         
         class(part2d), intent(inout) :: this
         class(ufield2d), pointer, intent(in) :: fd
         character(len=18), save :: sname = 'pmove:'
! local data
! list = (true,false) = list of particles leaving tiles found in push
         integer :: noff, nyp, nx, ny, nxv, nypmx, kstrt, nvp
         integer :: npbmx, nbmax, idimp, nppmx, ntmax, irc
         real, dimension(:,:,:), pointer :: ppart => null()
         integer, dimension(:,:), pointer :: ncl => null()
         integer, dimension(:,:,:), pointer :: ihole => null()
         integer, dimension(:), pointer :: kpic => null()
         logical :: list = .true.
         
         call this%err%werrfl2(class//sname//' started')
         
         noff = fd%getnoff(); ny = fd%getnd2()
         nx = fd%getnd1(); nyp = fd%getnd2p()
         npbmx = this%npbmx; nbmax = this%nbmaxp
         idimp = this%xdim; nppmx = this%nppmx0
         ntmax = this%ntmaxp; ppart => this%ppart
         kstrt = this%p%getlkstrt(); nvp = this%p%getlnvp()
         ncl => this%ncl; ihole => this%ihole; kpic => this%kpic
         irc = this%irc
! check if required size of buffer has increased
         if (szpbuf < idimp*npbmx*mxyp1) then
            if (szpbuf /= 0) deallocate(ppbuff)
! allocate new buffer
            allocate(ppbuff(idimp,npbmx,mxyp1))
            szpbuf = idimp*npbmx*mxyp1
         endif
! check if required size of buffers has increased
         if (szbufs < idimp*nbmax) then
            if (szbufs /= 0) deallocate(sbufl,sbufr,rbufl,rbufr)
! allocate new buffers
            allocate(sbufl(idimp,nbmax),sbufr(idimp,nbmax))
            allocate(rbufl(idimp,nbmax),rbufr(idimp,nbmax))
            szbufs = idimp*nbmax
         endif
! check if required size of buffers has increased
         if (sznbufs < 3*mx1) then
            if (sznbufs /= 0) deallocate(ncll,nclr,mcll,mclr)
! allocate new buffers
            allocate(ncll(3,mx1),nclr(3,mx1),mcll(3,mx1),mclr(3,mx1))
            sznbufs = 3*mx1
         endif
!
! first part of particle reorder on x and y cell with mx, my tiles:
! list of particles leaving tile already calculated by push
         if (list) then
! updates: ppart, ppbuff, sbufl, sbufr, ncl, ncll, nclr, irc
            call PPPORDERF2LA(ppart,ppbuff,sbufl,sbufr,ncl,ihole,ncll,nclr,&
            &idimp,nppmx,mx1,myp1,npbmx,ntmax,nbmax,irc)
            if (irc /= 0) then
               write (erstr,*) kstrt,'mporderf2a error: ntmax, irc=',ntmax,irc            
            endif
! calculate list of particles leaving tile
         else
! updates ppart, ppbuff, sbufl, sbufr, ncl, ihole, ncll, nclr, irc
            call PPPORDER2LA(ppart,ppbuff,sbufl,sbufr,kpic,ncl,ihole,ncll,nclr&
            &,noff,nyp,idimp,nppmx,nx,ny,mx,my,mx1,myp1,npbmx,ntmax,nbmax,irc)
            if (irc /= 0) then
               write (erstr,*) kstrt,'mporder2a error: ntmax, irc=',ntmax,irc
            endif
         endif
         if (irc /= 0) then
               call this%err%equit(class//sname//erstr)
               return
         endif
!
! move particles into appropriate spatial regions with MPI:
! updates rbufr, rbufl, mcll, mclr
         call PPPMOVE2(sbufr,sbufl,rbufr,rbufl,ncll,nclr,mcll,mclr,kstrt,&
         &nvp,idimp,nbmax,mx1)
!
! second part of particle reorder on x and y cell with mx, my tiles:
! updates ppart, kpic
         call PPPORDER2LB(ppart,ppbuff,rbufl,rbufr,kpic,ncl,ihole,mcll,mclr,&
         &idimp,nppmx,mx1,myp1,npbmx,ntmax,nbmax,irc)
         
         if (irc /= 0) then
            write (erstr,*) kstrt,'mporder2b error: nppmx, irc=',nppmx,irc
            call this%err%equit(class//sname//erstr)
            stop
         endif

         this%irc = irc

         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine pmove
!      
      subroutine extractpsi(this,psi,dex)
      
         implicit none
         
         class(part2d), intent(inout) :: this
         class(ufield2d), pointer, intent(in) :: psi
         real, intent(in) :: dex
         character(len=18), save :: sname = 'extractpsi'
! local data
         real, dimension(:,:,:), pointer :: ppsi
         integer :: noff, nyp, nx, nxv, nypmx
         
         call this%err%werrfl2(class//sname//' started')

         ppsi => psi%getrf(); noff = psi%getnoff()
         nyp = psi%getnd2p(); nx = psi%getnd1()
         nxv = size(ppsi,2); nypmx = size(ppsi,3)
         
         call WPGPSIPOST2L_QP(this%ppart,ppsi(1,:,:),this%kpic,this%qbm,noff,&
         &nyp,this%xdim,this%nppmx0,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex)

         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine extractpsi
!
      subroutine partcopy(this,fd)
      
         implicit none
         
         class(part2d), intent(inout) :: this
         class(ufield2d), pointer, intent(in) :: fd
! local data
         character(len=18), save :: sname = 'partcopy:'
         integer :: noff
         
         call this%err%werrfl2(class//sname//' started')

         noff = fd%getnoff()         
         
! copy ordered particle data for OpenMP
         call PPPMOVIN2L(this%part,this%ppart,this%kpic,this%npp,noff,&
         &this%nppmx0,this%xdim,this%npmax,mx,my,mx1,mxyp1,this%irc)
! check for errors
         if (this%irc /= 0) then
            write (erstr,*) 'PPPMOVIN2L overflow error, irc=', this%irc
            call this%err%equit(class//sname//erstr); return
         endif
                  
         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine partcopy
!
      subroutine partcopyback(this)
      
         implicit none
         
         class(part2d), intent(inout) :: this
! local data
         character(len=18), save :: sname = 'partcopyback:'
         
         call this%err%werrfl2(class//sname//' started')

         this%irc = 0

         call PPPCOPYOUT2(this%part,this%ppart,this%kpic,this%npp,&
         &this%npmax,this%nppmx0,this%xdim,mxyp1,this%irc)

         if (this%irc /= 0) then
            write (erstr,*) 'PPPCOPYOUT2 overflow error, irc=', this%irc
            call this%err%equit(class//sname//erstr); return
         endif

         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine partcopyback  
!
      subroutine pipesend_part2d(this,tag,id)
      
         implicit none
         
         class(part2d), intent(inout) :: this
         integer, intent(in) :: tag
         integer, intent(inout) :: id
! local data
         character(len=18), save :: sname = 'pipesend_part2d:'
         integer :: des, ierr
         
         
         call this%err%werrfl2(class//sname//' started')
         
         des = this%p%getidproc()+this%p%getlnvp()
         
         if (des >= this%p%getnvp()) then
            id = MPI_REQUEST_NULL         
            call this%err%werrfl2(class//sname//' ended')
            return
         endif
         
         call this%pcb()
                  
         call MPI_ISEND(this%part,this%npp*this%xdim,this%p%getmreal(),&
         &des,tag,this%p%getlworld(),id,ierr)

! check for errors
         if (ierr /= 0) then
            write (erstr,*) 'MPI_ISEND failed'
            call this%err%equit(class//sname//erstr); return
         endif

         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine pipesend_part2d
!      
      subroutine piperecv_part2d(this,fd,tag)
      
         implicit none
         
         class(part2d), intent(inout) :: this
         class(ufield2d), pointer, intent(in) :: fd
         integer, intent(in) :: tag
! local data
         character(len=18), save :: sname = 'piperecv_part2d:'
         integer, dimension(10) :: istat
         integer :: nps, id, des, ierr
         
         
         call this%err%werrfl2(class//sname//' started')

         des = this%p%getidproc()-this%p%getlnvp()
         
         if (des < 0) then
            call this%err%werrfl2(class//sname//' ended')
            return
         endif

         call MPI_IRECV(this%part,this%npmax*this%xdim,this%p%getmreal(),&
         &des,tag,this%p%getlworld(),id,ierr)

         call MPI_WAIT(id,istat,ierr)
         
         call MPI_GET_COUNT(istat,this%p%getmreal(),nps,ierr)

         this%npp = nps/this%xdim
         
         call this%pcp(fd)
         
! check for errors
         if (ierr /= 0) then
            write (erstr,*) 'MPI failed'
            call this%err%equit(class//sname//erstr); return
         endif
         
         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine piperecv_part2d
!      
      subroutine writehdf5_part2d(this,file,delta)

         implicit none
         
         class(part2d), intent(inout) :: this
         class(hdf5file), intent(in) :: file
         real, dimension(2), intent(in) :: delta
! local data
         character(len=18), save :: sname = 'writehdf5_part2d:'
         integer :: ierr

         call this%err%werrfl2(class//sname//' started')                  
         call this%pcb()
         call pwpart(this%p,this%err,file,this%part,this%npp,1,delta,ierr)
         call this%err%werrfl2(class//sname//' ended')
      
      end subroutine writehdf5_part2d
!      
      end module part2d_class