C* -------------------------------------------------------------------
C
C* The following code is Dr. Decyk's 1-D PIC code for single-processor
C* computers. [KD]
C*
C* -------------------------------------------------------------------

c * * * periodic 1d electrostatic particle simulation kernel code * * *
c 1D sm version  for public  
c This is a simple 1d skeleton particle-in-cell code designed for
c exploring new computer architectures.  it contains the critical pieces
c needed for depositing charge, advancing particles, and solving the
c field.  the code moves only electrons, with periodic electrostatic
c forces obtained by solving poisson's equation with fast fourier
c transforms.  the only diagnostic is particle and field energy.
cc portable gcpic kernel code, using algorithm described in:
c p. c. liewer and v. k. decyk, j. computational phys. 85, 302 (1989).
c written by viktor k. decyk, ucla
c
      program beps1k
c       parameter( indx =  14, npx =   4096000, npxb = 409600)
       parameter( indx =   7, npx =    5120, npxb = 512)

c       parameter( indx =   7, npx =   64768, npxb = 512)
c        parameter( indx =  18, npx =   409600, npxb = 40960)
c       parameter( indx =  18, npx =  4096000, npxb =   409600)
c r6
c        parameter( indx =  13, npx =  6096000, npxb =   409600)
c      parameter( indx =  18, npx = 40960000, npxb =  4096000)

c Time
      parameter( tend =  5.000, dt = 0.2000000e+00)
c      parameter( tend =  45.000, dt = 0.2000000e+00)
c Velocity
      parameter( vtx =   1.000,  vdx= 5.000, vtdx =   0.500)
      parameter( avdx=   5.000)
c Arch , vprocs, boundaries
      parameter( indnvp= 2, idps =    2, idimp =   2)
c we changed indnvp from 2 to 3
! Electrons
      parameter(np=npx+npxb)
      parameter(nx=2**indx,nxh=nx/2)
! Time steps
      parameter(nloop=tend/dt+.0001)
! Vprocs, partitions, particle/partition, 
      parameter(nvp=2**indnvp,nblok=1+1*(nvp-1))
      parameter(npmax=(np/nvp)*1.01+250,nxpmx=(nx-1)/nvp+4)
! Grids/dim
      parameter(kxp=(nxh-1)/nvp+1,kblok=1+1*(nxh/kxp-1))
      parameter(nbmax=1+(2*(npx*vtx+npxb*vtdx)+1.4*npxb*avdx)*dt/nx)
      parameter(ntmax=2*nbmax)

      parameter(npd=128)
      complex ffc, fc, qc, sct
      common /large/ part

      dimension part(idimp,npmax,nblok)
      dimension maskp(npmax,nblok)
      dimension q(nxpmx,nblok), fx(nxpmx,nblok)
      dimension qc(kxp,kblok), fc(kxp,kblok)
      dimension ffc(kxp,kblok)
      dimension mixup(kxp,kblok), sct(kxp,kblok)
      dimension edges(idps,nblok)
      dimension nxp(nblok), noff(nblok)
      dimension npp(nblok), nps(nblok)

      dimension sbufl(idimp,nbmax,nblok), sbufr(idimp,nbmax,nblok)
      dimension rbufl(idimp,nbmax,nblok), rbufr(idimp,nbmax,nblok)

      dimension ihole(ntmax,nblok)
      dimension jsl(idps,nblok), jsr(idps,nblok)
      dimension jss(idps,nblok), scr(idps,nblok)

      dimension nn(npd,nblok), nu(npd,nblok), nl(npd,nblok)
      dimension amx(npd,nblok), dxu(npd,nblok), dxl(npd,nblok)

      dimension tx(nxpmx,nblok), ty(nxpmx,nblok)
      double precision timebegin,timestop,t1,t2
      
  991 format (5h t = ,i7)
  992 format (19h * * * q.e.d. * * *)
  993 format (34h field, kinetic, total energies = ,3e14.7)

      data qme,ax /-1.,.8666667/
c ------------------------------------
c initialize 
c ------------------------------------
      kstrt = 1
c      open(unit=6,file='output1',form='formatted',status='unknown')
      
      itime = 0
      anx = float(nx)
      qtme = qme*dt
      affp = anx/float(np)
      zero = 0.
c ------------------------------
c calculate partition variables
c ------------------------------
      call dcomp1(edges,nxp,noff,nx,kstrt,nvp,idps,nblok)
c prepare fft tables

      isign = 0
      call pfft1r(qc,fc,isign,mixup,sct,indx,kstrt,kxp,kblok)
c calculate form factors
      call ppois1(qc,fc,isign,ffc,ax,affp,we,nx,kstrt,kxp,kblok)

c -------------------------------------------------------------
c initial density profile and maxwellian velocity distribution
c background electrons
c -------------------------------------------------------------
      do 120 k = 1, nblok
         nps(k) = 1
         npp(k) = 0
         nxp3 = nxp(k) + 3
         do 110 j = 1, nxp3
               q(j,k) = 0.
  110    continue
  120 continue
      if (npx.gt.0) call pistr1(part,edges,npp,nps,vtx,zero,npx,nx,
     1   idimp, npmax,nblok,idps)
c -----------------
c beam electrons
c -----------------
      do 140 k = 1, nblok
      nps(k) = npp(k) + 1
  140 continue
      if (npxb.gt.0) call pistr1 (part,edges,npp,nps,vtdx,vdx,npxb,nx,
     1  idimp, npmax,nblok,idps)

c ----------------------------------------
c deposit charge for initial distribution
c ----------------------------------------

      call pdost1 (part,q,npp,noff,qme,idimp,npmax,nblok,nxpmx)
c add background ion density
      qi0 = -qme/affp
      do 180 k = 1, nblok
      do 170 j = 1, nxp(k)
         q(j+1,k) = q(j+1,k) + qi0
  170 continue
  180 continue

      call timer(timebegin)
c =======================================
c
c * * * start main iteration loop * * *
c
c =======================================
  500 if (nloop.le.itime) go to 2000
      write (6,991) itime
c ----------------------------------
c transform charge to fourier space
c ----------------------------------

 
      isign = -1

c -------------------------------------------------------------------
c copy data from particle to field partition, and add up guard cells
c -------------------------------------------------------------------
      call timer(t1)
      call cppfp1 (q,qc,isign,scr,kstrt,nvp,nxpmx,nblok,kxp,kblok,idps)
      call timer(t2)
      print*,'after cppfp',t2-t1

      call timer(t1)
      call pfft1r(qc,fc,isign,mixup,sct,indx,kstrt,kxp,kblok)
      call timer(t2)
      print*,'after ifft',t2-t1

c ----------------------------------------
c calculate force/charge in fourier space
c ----------------------------------------
      call timer(t1)
      call ppois1 (qc,fc,isign,ffc,ax,affp,we,nx,kstrt,kxp,kblok)
      call timer(t2)
      print*,'after ppois',t2-t1
c -------------------------------------
c transform force/charge to real space
c -------------------------------------
      isign = 1
      call timer(t1)
      call pfft1r(fc,qc,isign,mixup,sct,indx,kstrt,kxp,kblok)
      call timer(t2)
      print*,'after dfft',t2-t1
c --------------------------------------------------------------------
c copy data from field to particle partition, and copy to guard cells
c --------------------------------------------------------------------
       call timer(t1)
      call cppfp1 (fx,fc,isign,scr,kstrt,nvp,nxpmx,nblok,kxp,kblok,idps)
      call timer(t2)
      print*,'after cppfp',t2-t1
c ----------------------------------------
c particle push and charge density update
c ----------------------------------------
c -----------------------------------
c initialize charge density to zero
c -----------------------------------
      do 1170 k = 1, nblok
      nxp3 = nxp(k) + 3
      do 1160 j = 1, nxp3
      q(j,k) = 0.
 1160 continue
 1170 continue
      wke = 0.
c ---------------
c push particles
c ---------------
      call timer(t1)
      call ppush1(part,fx,npp,noff,qtme,dt,wke,idimp,npmax,nblok,nxpmx)
      call timer(t2)
      print*, 'after ppush',t2-t1
c ------------------------------------------------
c move particles into appropriate spatial regions
c ------------------------------------------------' 
      call timer(t1)
      call pmove1 (part,edges,npp,sbufr,sbufl,rbufr,rbufl,ihole,jsr,jsl,
     1jss,nx,kstrt,nvp,idimp,npmax,nblok,idps,nbmax,ntmax,ierr)
      call timer(t2)
      print*,'after move',t2-t1
c ---------------
c deposit charge
c ---------------
      call timer(t1)
      call pdost1 (part,q,npp,noff,qme,idimp,npmax,nblok,nxpmx)
      call timer(t2)
      print*,'after pdost',t2-t1
c ----------------------
c add background charge
c ----------------------
      do 1190 k = 1, nblok
      do 1180 j = 1, nxp(k)
      q(j+1,k) = q(j+1,k) + qi0
 1180 continue
 1190 continue
c ------------------
c energy diagnostic
c ------------------
      wt = we + wke
      write (6,993) we, wke, wt
      itime = itime + 1
      go to 500
 2000 continue
c ====================================
c
c * * * end main iteration loop * * *
c
c ====================================
      write (6,992)
      call timer(timestop)
      print*, 'Total Time = ', timestop-timebegin,'seconds.' 
      stop
      end

      subroutine timer(t)

        integer val(8)
        double precision t

        call date_and_time(values=val)
        t = dble(val(8))*1d-3 + dble(val(7)) +          
     *   dble(val(6))*60d0 + dble(val(5))*3600d0

      end 





C =========================================================================
c
c  This subroutine determines spatial boundaries for particle
c  decomposition, calculates number of grid points in each spatial
c  region, and the offset of these grid points from the global address
c
C =========================================================================
      subroutine dcomp1(edges,nxp,noff,nx,kstrt,nvp,idps,nblok)
      dimension edges(idps,nblok)
      dimension nxp(nblok), noff(nblok)
      ks = kstrt - 2
      at1 = float(nx)/float(nvp)
      do 10 k = 1, nblok
         kb = k + ks
         edges(1,k) = at1*float(kb)
         noff(k) = edges(1,k) + .5
         edges(2,k) = at1*float(kb + 1)
         kr = edges(2,k) + .5
         nxp(k) = kr - noff(k)
   10 continue
      return
      end
C =========================================================================
c
c  This subroutine calculates initial particle co-ordinate
c  and velocity, with uniform density and maxwellian velocity with drift.
C =========================================================================
      subroutine pistr1 (part,edges,npp,nps,vtx,vdx,npx,nx,idimp,npmax,
     1nblok,idps)
c ranorm = gaussian random number with zero mean and unit variance
c with spatial decomposition
      double precision ranorm
      dimension part(idimp,npmax,nblok)
      dimension edges(idps,nblok), npp(nblok), nps(nblok)
c particle distribution
      at1 = float(nx)/float(npx)
      do 20 j = 1, npx
c uniform density profile
         xt = at1*(float(j) - .5)
c maxwellian velocity distribution
         vxt = vtx*ranorm(d)
         do 10 k = 1, nblok
         if ((xt.ge.edges(1,k)).and.(xt.lt.edges(2,k))) then
            npp(k) = npp(k) + 1
            part(1,npp(k),k) = xt
            part(2,npp(k),k) = vxt
         endif
   10    continue
   20 continue
c add correct drift
      sum1 = 0.
      do 40 k = 1, nblok
      sum0 = 0.
      do 30 j = nps(k), npp(k)
         sum0 = sum0 + part(2,j,k)
   30 continue
      sum1 = sum1 + sum0
   40 continue

      sum1 = sum1/float(npx) - vdx
      do 60 k = 1, nblok
         do 50 j = nps(k), npp(k)
            part(2,j,k) = part(2,j,k) - sum1
   50    continue
   60 continue
      return
      end
C =========================================================================
c
c  This subroutine calculates particle charge density
c  using second-order spline interpolation, with periodic boundaries.
c
C =========================================================================
      subroutine pdost1 (part,q,npp,noff,qm,idimp,npmax,nblok,nxpmx)
      dimension part(idimp,npmax,nblok), q(nxpmx,nblok)
      dimension npp(nblok), noff(nblok)
      do 20 k = 1, nblok
         nnoff = noff(k) - 2
         do 10 j = 1, npp(k)
            nn = part(1,j,k) + .5
            dx = part(1,j,k) - float(nn)
            nn = nn - nnoff
            q(nn,k) = q(nn,k) + qm*(.75 - dx*dx)
            nn = nn + 1
            q(nn,k) = q(nn,k) + .5*qm*(.5 + dx)**2
            nn = nn - 2
            q(nn,k) = q(nn,k) + .5*qm*(.5 - dx)**2
   10    continue
   20 continue
      return
      end

C =========================================================================
c
c thus subroutine copies data between particle and field partitions.
c if isign = -1, data is copied from particle to field partition.
c if isign = 1, data is copied from field to particle partition.
c
C =========================================================================
      subroutine cppfp1 (f,fc,isign,scr,kstrt,nvp,nxpmx,nblok,kxp,kblok,
     1idps)
      complex fc
      dimension f(nxpmx,nblok), fc(kxp,kblok)
      dimension scr(idps,nblok)
      ks = kstrt - 2
      moff = nxpmx*nvp
      if (isign) 10, 60, 70
c ---------------------------------------
c copy from particle to field partition
c ---------------------------------------
   10 do 30 k = 1, kblok
         do 20 j = 1, kxp
            fc(j,k) = cmplx(f(2*j,k),f(2*j+1,k))
   20    continue
   30 continue


      do 50 k = 1, kblok
        kl = k + ks
        if (kl.lt.1) then
           scr(1,k) = f(2*kxp+2,kl+nvp)
           scr(2,k) = f(2*kxp+3,kl+nvp)
        else
           scr(1,k) = f(2*kxp+2,kl)
           scr(2,k) = f(2*kxp+3,kl)
        endif
        fc(1,k) = fc(1,k) + cmplx(scr(1,k),scr(2,k))

        tterol=scr(1,k)


        kr = k + ks + 2
        if (kr.gt.nvp) then
           scr(1,k) = f(1,kr-nvp)
        else
           scr(1,k) = f(1,kr)
        endif
        fc(kxp,k) = fc(kxp,k) + cmplx(0.,scr(1,k))
         
  50  continue

   60  return
 
c --------------------------------------
c copy from field to particle partition
c --------------------------------------
   70 do 90 k = 1, kblok
         do 80 j = 1, kxp
            f(2*j,k) = real(fc(j,k))
            f(2*j+1,k) = aimag(fc(j,k))
   80    continue
   90 continue

      do 100 k = 1, kblok
         kl = k + ks
         if (kl.lt.1) then
            f(1,k) = aimag(fc(kxp,kl+nvp))
         else
            f(1,k) = aimag(fc(kxp,kl))
         endif
         kr = k + ks + 2
         if (kr.gt.nvp) then
            f(2*kxp+2,k) = real(fc(1,kr-nvp))
            f(2*kxp+3,k) = aimag(fc(1,kr-nvp))
         else
            f(2*kxp+2,k) = real(fc(1,kr))
            f(2*kxp+3,k) = aimag(fc(1,kr))
         endif
 100  continue

      return
      end

C =========================================================================
c  This subroutine solves 1d poisson's equation in fourier space for
c  force/charge (or convolution of electric field over particle shape)
c  with periodic boundary conditions.
c  The equation used is:
c    fx(k) = -sqrt(-1)*k*g(k)*s(k), where k = 2pi*j/nx, j = fourier mode,
c    g(k) = (affp/k**2)*s(k), and s(k) = exp(-(k*ax)**2/2), except for
c    fx(k=0) = fx(k=pi) = 0.
c
c  If isign is not equal to 0, force/charge is calculated.
c   aimag(ffc(j,k)) = particle shape factor s for fourier mode jj - 1
c   real(ffc(j,k)) = potential green's function g for fourier mode jj - 1
c
C =========================================================================
      subroutine ppois1 (q,fx,isign,ffc,ax,affp,wp,nx,kstrt,kxp,kblok)
      complex q, fx, ffc
      dimension q(kxp,kblok), fx(kxp,kblok)
      dimension ffc(kxp,kblok)
      nxh = nx/2
      ks = kstrt - 2
      dnx = 6.28318530717959/float(nx)
      if (isign.ne.0) go to 30
      if (kstrt.gt.nxh) return
c ----------------------------
c prepare form factor array
c ----------------------------
      do 20 k = 1, kblok
         koff = kxp*(k + ks) - 1
         do 10 j = 1, kxp
            dkx = dnx*float(j + koff)
            at1 = exp(-.5*(dkx*ax)**2)
            if (dkx.eq.0.) then
               ffc(j,k) = cmplx(0.,1.)
            else
               ffc(j,k) = cmplx(affp*at1/(dkx*dkx),at1)
            endif
   10    continue
   20 continue
      return
c ---------------------------------------------
c calculate force/charge and sum field energy
c ---------------------------------------------
   30 wp = 0.
      if (kstrt.gt.nxh) go to 60
      do 50 k = 1, kblok
         koff = kxp*(k + ks) - 1
         sum1 = 0.
         do 40 j = 1, kxp
            at1 = real(ffc(j,k))*aimag(ffc(j,k))
            at2 = dnx*float(j + koff)*at1
            fx(j,k) = at2*cmplx(aimag(q(j,k)),-real(q(j,k)))
            sum1 = sum1 + at1*real(q(j,k)*conjg(q(j,k)))
   40    continue
         wp = wp + sum1
   50 continue
   60 continue

      wp = float(nx)*wp
      return
      end

C =========================================================================
c
c  This subroutine updates particle co-ordinate and velocity
c  using leap-frog scheme in time and second-order spline interpolation
c  in space, with periodic boundary conditions.
c  equations used are:
c    v(t+dt/2) = v(t-dt/2) + (q/m)*fx(x(t))*dt, where q/m is charge/mass,
c    x(t+dt) = x(t) + v(t+dt/2)*dt
c  fx(x(t)) is approximated by interpolation from the nearest grid points
c
C =========================================================================
      subroutine ppush1(part,fx,npp,noff,qtm,dt,ek,idimp,npmax,nblok,
     1nxpmx)
      dimension part(idimp,npmax,nblok), fx(nxpmx,nblok)
      dimension npp(nblok), noff(nblok)
      sum1 = 0.
      do 20 k = 1, nblok
      nnoff = noff(k) - 2
      do 10 j = 1, npp(k)
c ----------------------------
c find interpolation weights
c ----------------------------
         nn = part(1,j,k) + .5
         dx = part(1,j,k) - float(nn)
         nn = nn - nnoff
c ----------------------------
c find acceleration
c ----------------------------
         dx = (.75 - dx*dx)*fx(nn,k) + .5*(fx(nn+1,k)*(.5 + dx)**2 
     1        + fx(nn-1,k)*(.5 - dx)**2)
c ----------------------------
c new velocity
c ----------------------------
         dx = part(2,j,k) + qtm*dx
c ----------------------------
c average kinetic energy
c ----------------------------
         sum1 = sum1 + (part(2,j,k) + dx)**2
         part(2,j,k) = dx
c ----------------------------
c new position
c ----------------------------
         part(1,j,k) = part(1,j,k) + dx*dt
   10 continue
   20 continue

c normalize kinetic energy
      ek = ek + .125*sum1
      return
      end

C =========================================================================
c
c This subroutine moves particles into appropriate spatial regions
c
C =========================================================================
      subroutine pmove1 (part,edges,npp,sbufr,sbufl,rbufr,rbufl,ihole,js
     1r,jsl,jss,nx,kstrt,nvp,idimp,npmax,nblok,idps,nbmax,ntmax,ierr)
      dimension part(idimp,npmax,nblok)
      dimension edges(idps,nblok), npp(nblok)
      dimension sbufl(idimp,nbmax,nblok), sbufr(idimp,nbmax,nblok)
      dimension rbufl(idimp,nbmax,nblok), rbufr(idimp,nbmax,nblok)
      dimension jsl(idps,nblok), jsr(idps,nblok), jss(idps,nblok)
      dimension ihole(ntmax,nblok)
      anx = float(nx)
      ks = kstrt - 2
      iter = 4
c ----------------------------------------------------------------
c debugging section: count total number of particles before move
c ----------------------------------------------------------------
      npr = 0
      do 5 k = 1, nblok
         npr = npr + npp(k)
    5 continue

c --------------------------
c buffer outgoing particles
c --------------------------
   10 do 30 k = 1, nblok
         kb = k + ks
         jsl(1,k) = 0
         jsr(1,k) = 0
         jss(2,k) = 0
         do 20 j = 1, npp(k)
            xt = part(1,j,k)

            if (xt.ge.edges(2,k)) then
               if (jsr(1,k).lt.nbmax) then
                  jsr(1,k) = jsr(1,k) + 1
                  if ((kb+1).eq.nvp) xt = xt - anx
                  sbufr(1,jsr(1,k),k) = xt
                  sbufr(2,jsr(1,k),k) = part(2,j,k)
                  ihole(jsl(1,k)+jsr(1,k),k) = j
               else
                  jss(2,k) = 1
               endif

            elseif (xt.lt.edges(1,k)) then
               if (jsl(1,k).lt.nbmax) then
                  jsl(1,k) = jsl(1,k) + 1
                  if (kb.eq.0) xt = xt + anx
                  sbufl(1,jsl(1,k),k) = xt
                  sbufl(2,jsl(1,k),k) = part(2,j,k)
                  ihole(jsl(1,k)+jsr(1,k),k) = j
               else
                  jss(2,k) = 1
            endif
         endif
   20    continue
         jss(1,k) = jsl(1,k) + jsr(1,k)
   30 continue
c --------------------------------
c check for full buffer condition
c --------------------------------
      ibflg = 0
      do 40 k = 1, nblok
         ibflg = ibflg + jss(2,k)
   40 continue

c -------------------------
c copy particle buffers
c -------------------------
   50 iter = iter + 4
      do 80 k = 1, nblok

      kl = k + ks
      if (kl.lt.1) then
        jsl(2,k) = jsr(1,kl+nvp)
      else
        jsl(2,k) = jsr(1,kl)
      endif

      do 60 j = 1, jsl(2,k)
         if (kl.lt.1) then
            rbufl(1,j,k) = sbufr(1,j,kl+nvp)
            rbufl(2,j,k) = sbufr(2,j,kl+nvp)
         else
            rbufl(1,j,k) = sbufr(1,j,kl)
            rbufl(2,j,k) = sbufr(2,j,kl)
         endif
   60 continue

       kr = k + ks + 2
       if (kr.gt.nvp) then
         jsr(2,k) = jsl(1,kr-nvp)
       else
         jsr(2,k) = jsl(1,kr)
       endif

      do 70 j = 1, jsr(2,k)
      if (kr.gt.nvp) then
        rbufr(1,j,k) = sbufl(1,j,kr-nvp)
        rbufr(2,j,k) = sbufl(2,j,kr-nvp)
      else
        rbufr(1,j,k) = sbufl(1,j,kr)
        rbufr(2,j,k) = sbufl(2,j,kr)
      endif
  70  continue

   80 continue
c ------------------------------------------- 
c check if particles must be passed further
c ------------------------------------------- 
      nps = 0
      do 110 k = 1, nblok
      jsl(1,k) = 0
      do 90 j = 1, jsr(2,k)
         if (rbufr(1,j,k).lt.edges(1,k)) jsl(1,k) = 1
   90 continue
      jsr(1,k) = 0
      do 100 j = 1, jsl(2,k)
      if (rbufl(1,j,k).ge.edges(2,k)) jsr(1,k) = 1
  100 continue
      nps = nps + (jsl(1,k) + jsr(1,k))
  110 continue

      if (nps.eq.0) go to 150
c -------------------------------------------
c remove particles which do not belong here
c -------------------------------------------
      do 140 k = 1, nblok
      kb = k + ks
      jsl(1,k) = 0
      jss(2,k) = 0
      do 120 j = 1, jsr(2,k)
      xt = rbufr(1,j,k)
      if (xt.lt.edges(1,k)) then
         jsl(1,k) = jsl(1,k) + 1
         if (kb.eq.0) xt = xt + anx
         sbufl(1,jsl(1,k),k) = xt
         sbufl(2,jsl(1,k),k) = rbufr(2,j,k)
      else
         jss(2,k) = jss(2,k) + 1
         rbufr(1,jss(2,k),k) = xt
         rbufr(2,jss(2,k),k) = rbufr(2,j,k)
      endif
  120 continue
      jsr(2,k) = jss(2,k)

      jsr(1,k) = 0
      jss(2,k) = 0
      do 130 j = 1, jsl(2,k)
      xt = rbufl(1,j,k)
      if (xt.ge.edges(2,k)) then
         jsr(1,k) = jsr(1,k) + 1
         if ((kb+1).eq.nvp) xt = xt - anx
         sbufr(1,jsr(1,k),k) = xt
         sbufr(2,jsr(1,k),k) = rbufl(2,j,k)
      else
         jss(2,k) = jss(2,k) + 1
         rbufl(1,jss(2,k),k) = xt
         rbufl(2,jss(2,k),k) = rbufl(2,j,k)
      endif
  130 continue
      jsl(2,k) = jss(2,k)
  140 continue
c -------------------------------------------
c check if move would overflow particle array
c -------------------------------------------
  150 ierr = 0
      do 160 k = 1, nblok
      if ((npp(k)+jsl(2,k)+jsr(2,k)-jss(1,k)).gt.npmax) then
         jss(2,k) = 1
      else
         jss(2,k) = 0
      endif
      ierr = ierr + jss(2,k)
  160 continue

      if (ierr.gt.0) then
         write (6,*) 'particle overflow error, ierr =', ierr
         return
      endif
      do 200 k = 1, nblok
c ------------------------------------------------------
c distribute incoming particles from buffers
c ------------------------------------------------------
      jss(2,k) = min0(jss(1,k),jsl(2,k))
      do 170 j = 1, jss(2,k)
      part(1,ihole(j,k),k) = rbufl(1,j,k)
      part(2,ihole(j,k),k) = rbufl(2,j,k)
  170 continue
      if (jss(1,k).gt.jsl(2,k)) then
         jss(2,k) = min0(jss(1,k)-jsl(2,k),jsr(2,k))
      else
         jss(2,k) = jsl(2,k) - jss(1,k)
      endif
      do 180 j = 1, jss(2,k)

      if (jss(1,k).gt.jsl(2,k)) then
         part(1,ihole(j+jsl(2,k),k),k) = rbufr(1,j,k)
         part(2,ihole(j+jsl(2,k),k),k) = rbufr(2,j,k)
      else
         part(1,j+npp(k),k) = rbufl(1,j+jss(1,k),k)
         part(2,j+npp(k),k) = rbufl(2,j+jss(1,k),k)
      endif
  180 continue
      if (jss(1,k).le.jsl(2,k)) then
         npp(k) = npp(k) + (jsl(2,k) - jss(1,k))
         jss(1,k) = jsl(2,k)
      endif
      jss(2,k) = jss(1,k) - (jsl(2,k) + jsr(2,k))
      if (jss(2,k).gt.0) then
         jss(1,k) = (jsl(2,k) + jsr(2,k))
         jsr(2,k) = jss(2,k)
      else
         jss(1,k) = jss(1,k) - jsl(2,k)
         jsr(2,k) = -jss(2,k)
      endif
      do 190 j = 1, jsr(2,k)
c ------------------------------------------------------
c holes left over
c fill up remaining holes in particle array with particles from bottom
c ------------------------------------------------------
      if (jss(2,k).gt.0) then
         j1 = npp(k) - j + 1
         j2 = jss(1,k) + jss(2,k) - j + 1
         if (j1.gt.ihole(j2,k)) then
c ------------------------------------------------------
c move particle only if it is below current hole
c ------------------------------------------------------
            part(1,ihole(j2,k),k) = part(1,j1,k)
            part(2,ihole(j2,k),k) = part(2,j1,k)
         endif
      else
c ------------------------------------------------------
c no more holes
c distribute remaining particles from right into bottom
c ------------------------------------------------------
         part(1,j+npp(k),k) = rbufr(1,j+jss(1,k),k)
         part(2,j+npp(k),k) = rbufr(2,j+jss(1,k),k)
      endif
  190 continue
      if (jss(2,k).gt.0) then
         npp(k) = npp(k) - jsr(2,k)
      else
         npp(k) = npp(k) + jsr(2,k)
      endif
      jss(1,k) = 0
  200 continue
c ------------------------------------------------------
c check if any particles have to be passed further
c ------------------------------------------------------
      if (nps.gt.0) go to 50
c check if buffer overflowed and more particles remain to be checked
      if (ibflg.gt.0) go to 10
c debugging section: count total number of particles after move
      nps = 0
      do 205 k = 1, nblok
      nps = nps + npp(k)
  205 continue

      if (npr.ne.nps) write (6,*) 'number error, np-old,np-new=',npr,nps
      return
      end

C =========================================================================
c
c  This subroutine performs a one dimensional real to complex fast
c  fourier transform and its inverse, using complex arithmetic,
c  for data which is distributed in blocks
c    f = input and output data
c    g = scratch array
c    indx = power of 2 which determines length of transform, nx = 2**indx
c
c    if isign = 0, the fft tables are prepared
c    if isign = -1, an inverse fourier transform is performed
c
c    f(n,m) = (1/nx)*sum(f(j,k)*exp(-sqrt(-1)*2pi*nn*jj/nx))
c    if isign = 1, a forward fourier transform is performed
c    f(j,k) = sum(f(n,m)*exp(sqrt(-1)*2pi*nn*jj/nx))
c      where nn = n + kxp*(m - 1) and jj = j + kxp*(k - 1)
c written by viktor k. decyk, ucla
C =========================================================================
      subroutine pfft1r(f,g,isign,mixup,sct,indx,kstrt,kxp,kblok)
      complex f, g, sct, t1
      dimension f(kxp,kblok), g(kxp,kblok)
      dimension mixup(kxp,kblok), sct(kxp,kblok)
      nx = 2**indx
      nxh = nx/2
      if (kstrt.gt.nxh) return
      indx1 = indx - 1
      kbs = nxh/kxp
      kxp2 = kxp + 2
      ks = kstrt - 2
      dnx = 6.28318530717959/float(nx)
      if (isign) 60, 10, 220
c ---------------------------
c prepare fft tables
c bit-reverse index table:
c mixup(j,k) = 1 + reversed bits of (j + kxp*(k + kstrt - 2) - 1)
c ---------------------------
   10 call bitrv1(f,g,mixup,isign,indx1,kstrt,kxp,kblok)
c ---------------------------
c local sine/cosine table for the angles 2*n*pi/kxp
c ---------------------------
      do 50 k = 1, kblok
      do 40 j = 1, kxp
         arg = dnx*float((nxh*(j - 1))/kxp)
         sct(j,k) = cmplx(cos(arg),-sin(arg))
   40 continue
   50 continue
      return
c ---------------------------
c inverse fourier transform
c bit-reverse array elements to temporary
c ---------------------------
   60 call bitrv1(f,g,mixup,isign,indx1,kstrt,kxp,kblok)
 
      do 160 l = 1, indx1
      nxs = 2**(l - 1)
      kxs = nxs/kxp
c ---------------------------
c local calculation
c ---------------------------
      if (kxs.eq.0) then
         do 110 k = 1, kblok
         do 90 j = 1, kxp
            lb = (j - 1)/nxs
            if (lb.eq.(2*(lb/2))) then
               f(j,k) = g(j+nxs,k)
            else
               f(j,k) = g(j-nxs,k)
            endif
   90    continue
c ---------------------------
c perform reduction
c ---------------------------
         km = kxp/nxs
         do 100 j = 1, kxp
            kb = j - 1
            lb = kb/nxs
            t1 = sct(1+km*(kb-nxs*lb),k)
            if (lb.eq.(2*(lb/2))) then
               g(j,k) = g(j,k) + t1*f(j,k)
            else
               g(j,k) = f(j,k) - t1*g(j,k)
            endif
  100    continue
  110    continue
      else
c ---------------------------
c copy data
c ---------------------------
         do 130 k = 1, kblok
            kb = k + ks
            lb = kb/kxs
            kb = kb + 1
        do 120 j = 1, kxp
        if (lb.eq.(2*(lb/2))) then
           f(j,k) = g(j,kb+kxs)
        else
           f(j,k) = g(j,kb-kxs)
        endif
 120    continue

  130    continue

c ------------------
c perform reduction
c ------------------
         km = nxh/nxs
         dns = dnx*float(km)
         do 150 k = 1, kblok
            kb = k + ks
            lb = kb/kxs
            kb = kxp*(kb - kxs*lb) - 1
            lb = lb - 2*(lb/2)
            do 140 j = 1, kxp
               arg = dns*float(j + kb)
               t1 = cmplx(cos(arg),-sin(arg))
               if (lb.eq.0) then
                  g(j,k) = g(j,k) + t1*f(j,k)
               else
                  g(j,k) = f(j,k) - t1*g(j,k)
               endif
  140       continue
  150    continue
      endif
  
  160 continue



c unscramble coefficients and normalize result
      do 180 k = 1, kblok
         kb = k + ks
         kk = kbs - kb
      do 170 j = 2, kxp
         f(j,k) = g(j,kk)
 170  continue
      if (kb.gt.0) f(1,k) = g(1,kk+1)

  180  continue

c ------------------
c perform reduction
c ------------------
         ani = .5/float(nx)
         do 210 k = 1, kblok
            kb = k + ks
            koff = kxp*kb - 1
            do 190 j = 2, kxp
               jj = kxp2 - j
               arg = dnx*float(j + koff)
               g(j,k) = ani*((g(j,k) + conjg(f(jj,k))) 
     1            - (g(j,k) 
     2            - conjg(f(jj,k)))*cmplx(sin(arg),cos(arg)))
  190       continue
            if (kb.eq.0) then
               g(1,k) = 2.*ani*cmplx(real(g(1,k)) + 
     1                   aimag(g(1,k)),real(g(1,k))
     2                   - aimag(g(1,k)))
            else
               arg = dnx*float(1 + koff)
               g(1,k) = ani*((g(1,k) + conjg(f(1,k))) - 
     1                   (g(1,k) - conjg(f(1,k)))*
     2                   cmplx(sin(arg),cos(arg)))
            endif
c copy result to destination
            do 200 j = 1, kxp
               f(j,k) = g(j,k)
  200       continue
  210 continue

      return
c ------------------
c forward fourier transform
c scramble coefficients
c ------------------
  220 do 240 k = 1, kblok
         kb = k + ks
         kk = kbs - kb
      do 230 j = 2, kxp
         g(j,k) = f(j,kk)
  230 continue
      if (kb.gt.0) g(1,k) = f(1,kk+1)

  240 continue
c ------------------
c perform reduction
c ------------------
      do 270 k = 1, kblok
         kb = k + ks
         koff = kxp*kb - 1
         do 250 j = 2, kxp
            jj = kxp2 - j
            arg = dnx*float(j + koff)
            f(j,k) = (f(j,k) + conjg(g(jj,k))) - 
     1               (f(j,k) - conjg(g(jj,k)))*cmplx(sin(arg),-cos(arg))
  250    continue
         if (kb.eq.0) then
            f(1,k) = cmplx(real(f(1,k)) + aimag(f(1,k)),real(f(1,k)) - 
     1                aimag(f(1,k)))
         else
            arg = dnx*float(1 + koff)
            f(1,k) = (f(1,k) + conjg(g(1,k))) - 
     1               (f(1,k) - conjg(g(1,k)))*cmplx(sin(arg),-cos(arg))
         endif
c copy result to temporary
         do 260 j = 1, kxp
            g(j,k) = f(j,k)
  260    continue
  270 continue
c ----------------------------
c bit-reverse array elements
c ----------------------------
      call bitrv1(g,f,mixup,isign,indx1,kstrt,kxp,kblok)
      do 370 l = 1, indx1
         nxs = 2**(l - 1)
         kxs = nxs/kxp
c -------------------------
c local calculation
c -------------------------
         if (kxs.eq.0) then
c copy local data
            do 320 k = 1, kblok
               do 300 j = 1, kxp
                  lb = (j - 1)/nxs
                  if (lb.eq.(2*(lb/2))) then
                     g(j,k) = f(j+nxs,k)
                  else
                     g(j,k) = f(j-nxs,k)
                  endif
  300          continue
c perform reduction
               km = kxp/nxs
               do 310 j = 1, kxp
                 kb = j - 1
                 lb = kb/nxs
                 t1 = conjg(sct(1+km*(kb-nxs*lb),k))
                 if (lb.eq.(2*(lb/2))) then
                    f(j,k) = f(j,k) + t1*g(j,k)
                 else
                    f(j,k) = g(j,k) - t1*f(j,k)
                 endif
  310          continue
  320       continue
         else
            do 340 k = 1, kblok
               kb = k + ks
               lb = kb/kxs
               kb = kb + 1
               lb = lb - 2*(lb/2)
        do 330 j = 1, kxp
           if (lb.eq.0) then
              g(j,k) = f(j,kb+kxs)
           else
              g(j,k) = f(j,kb-kxs)
           endif
 330    continue

  340    continue
c ------------------
c perform reduction
c ------------------
         km = nxh/nxs
         dns = dnx*float(km)
         do 360 k = 1, kblok
         kb = k + ks
         lb = kb/kxs
         kb = kxp*(kb - kxs*lb) - 1
         lb = lb - 2*(lb/2)
         do 350 j = 1, kxp
         arg = dns*float(j + kb)
         t1 = cmplx(cos(arg),sin(arg))
         if (lb.eq.0) then
            f(j,k) = f(j,k) + t1*g(j,k)
         else
            f(j,k) = g(j,k) - t1*f(j,k)
         endif
  350    continue
  360    continue
      endif
  370 continue
      return
      end
C =========================================================================
c this subroutine performs bit-reversal of distributed complex data,
c that is, it calculates  
c    g(j,k) = f(jj,kk), where
c    kk = (mixup(j,k)-1)/kxp+1 and jj = mixup(j,k)-kxp*(kk-1), and
c    mixup(j,k) = 1 + reversed bits of (j + kxp*(k + kstrt - 2) - 1)
c    f = complex input array
c    g = complex output array
c    mixup = array of bit reversed addresses
c    if isign = 0, bit-reverse array is prepared
c    if isign = 1, bit-reversal is performed
c written by viktor k. decyk, ucla
C =========================================================================
      subroutine bitrv1(f,g,mixup,isign,indx,kstrt,kxp,kblok)
      complex f, g
      dimension f(kxp,kblok), g(kxp,kblok), mixup(kxp,kblok)
      nx = 2**indx
      if (kstrt.gt.nx) return
      ks = kstrt - 2
      if (isign.ne.0) go to 40
c ---------------------------------
c prepare bit-reverse index table
c ---------------------------------
      do 30 k = 1, kblok
         koff = kxp*(k + ks) - 1
         do 20 j = 1, kxp
            lb = j + koff
            ll = 0
            do 10 l = 1, indx
               jb = lb/2
               it = lb - 2*jb
               lb = jb
               ll = 2*ll + it
   10       continue
            mixup(j,k) = ll + 1
   20    continue
   30 continue
      return
c ---------------------------------
c bit-reverse array elements
c ---------------------------------
  40  do 60 k = 1, kblok
      do 50 j = 1, kxp
        ll = mixup(j,k)
        kk = (ll - 1)/kxp + 1
        jj = ll - kxp*(kk - 1)
        g(j,k) = f(jj,kk)
  50  continue
  60  continue

      return
      end
C =========================================================================
c this program calculates a random number y from a gaussian distribution
c with zero mean and unit variance, according to the method of
c mueller and box:
c    y(k) = (-2*ln(x(k)))**1/2*sin(2*pi*x(k+1))
c    y(k+1) = (-2*ln(x(k)))**1/2*cos(2*pi*x(k+1)),
c where x is a random number uniformly distributed on (0,1).
c written for the ibm by viktor k. decyk, ucla
C =========================================================================
      function ranorm(d)
      integer r1,r2,r4,r5
      double precision ranorm,h1l,h1u,h2l,r0,r3,asc,bsc,temp
      save iflg,r1,r2,r4,r5,h1l,h1u,h2l,r0
      data r1,r2,r4,r5 /885098780,1824280461,1396483093,55318673/
      data h1l,h1u,h2l /65531.0d0,32767.0d0,65525.0d0/
      data iflg,r0 /0,0.0d0/
      if (iflg.eq.0) go to 10
      ranorm = r0
      r0 = 0.0d0
      iflg = 0
      return
   10 isc = 65536
      asc = dfloat(isc)
      bsc = asc*asc
      i1 = r1 - (r1/isc)*isc
      r3 = h1l*dfloat(r1) + asc*h1u*dfloat(i1)
      i1 = r3/bsc
      r3 = r3 - dfloat(i1)*bsc
      bsc = 0.5d0*bsc
      i1 = r2/isc
      isc = r2 - i1*isc
      r0 = h1l*dfloat(r2) + asc*h1u*dfloat(isc)
      asc = 1.0d0/bsc
      isc = r0*asc
      r2 = r0 - dfloat(isc)*bsc
      r3 = r3 + (dfloat(isc) + 2.0d0*h1u*dfloat(i1))
      isc = r3*asc
      r1 = r3 - dfloat(isc)*bsc
      temp = dsqrt(-2.0d0*dlog((dfloat(r1) + dfloat(r2)*asc)*asc))
      isc = 65536
      asc = dfloat(isc)
      bsc = asc*asc
      i1 = r4 - (r4/isc)*isc
      r3 = h2l*dfloat(r4) + asc*h1u*dfloat(i1)
      i1 = r3/bsc
      r3 = r3 - dfloat(i1)*bsc
      bsc = 0.5d0*bsc
      i1 = r5/isc
      isc = r5 - i1*isc
      r0 = h2l*dfloat(r5) + asc*h1u*dfloat(isc)
      asc = 1.0d0/bsc
      isc = r0*asc
      r5 = r0 - dfloat(isc)*bsc
      r3 = r3 + (dfloat(isc) + 2.0d0*h1u*dfloat(i1))
      isc = r3*asc
      r4 = r3 - dfloat(isc)*bsc
      r0 = 6.28318530717959d0*((dfloat(r4) + dfloat(r5)*asc)*asc)
      ranorm = temp*dsin(r0)
      r0 = temp*dcos(r0)
      iflg = 1
      return
      end


