c * * * periodic 1d electrostatic particle simulation kernel code * * *
c for IBM SP1

      program beps1k

      parameter( indx =   7, npx =    5120, npxb = 512)
cbenc1 parameter( indx =  11, npx =   409600, npxb = 40960)
cbenc2 parameter( indx =  14, npx =  4096000, npxb =   409600)
cbenc3 parameter( indx =  18, npx = 40960000, npxb =  4096000)

c Time
      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
cdm/sm mshare = (0,1) = (no,yes) architecture is shared memory
      parameter( indnvp= 2, idps =    2, idimp =   2, mshare =   0)
! 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+mshare*(nvp-1))
      parameter(npmax=(np/nvp)*1.01+250,nxpmx=(nx-1)/nvp+4)
! Grids/dim
      parameter(kxp=(nxh-1)/nvp+1,kblok=1+mshare*(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)
  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 for parallel processing
c ------------------------------------
      call PPINIT(idproc,nvp)
      kstrt = idproc + 1
      if (kstrt.eq.1) then
         open(unit=6,file='output1',form='formatted',status='unknown')
      endif
c initialize timer
      call TIMERA(-1,'total   ')
c initialize constants
      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
c        do 110 j = 1, nxpmx
c           if (j.le.nxp3) then
               q(j,k) = 0.
c           endif
  110    continue
  120 continue
      if (npx.gt.0) call PISTR1 (part,edges,npp,nps,vtx,zero,npx,nx,idim
     1p,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,id
     1imp,npmax,nblok,idps)
c ----------------------------------------
c deposit charge for initial distribution
c ----------------------------------------
      call PDOST1 (part,q,npp,noff,qme,idimp,npmax,nblok,nxpmx)
c     call pxost1 (part,q,tx,ty,npp,nxp,noff,nn,amx,dxu,dxl,qme,idimp,np
c    1max,nblok,nxpmx,npd)
c     call pyost1 (part,q,qt,npp,nxp,noff,qme,idimp,npmax,nblok,nxpmx,np
c    1d)
c add background ion density
      qi0 = -qme/affp
      do 180 k = 1, nblok
      do 170 j = 1, nxp(k)
c     do 170 j = 1, nxpmx
c     if (j.le.nxp(k)) then
         q(j+1,k) = q(j+1,k) + qi0
c     endif
  170 continue
  180 continue
c =======================================
c
c * * * start main iteration loop * * *
c
c =======================================
  500 if (nloop.le.itime) go to 2000
      if (kstrt.eq.1) 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 cppfp1 (q,qc,isign,scr,kstrt,nvp,nxpmx,nblok,kxp,kblok,idps)
      call pfft1r(qc,fc,isign,mixup,sct,indx,kstrt,kxp,kblok)
c ----------------------------------------
c calculate force/charge in fourier space
c ----------------------------------------
      call ppois1 (qc,fc,isign,ffc,ax,affp,we,nx,kstrt,kxp,kblok)
c -------------------------------------
c transform force/charge to real space
c -------------------------------------
      isign = 1
      call pfft1r(fc,qc,isign,mixup,sct,indx,kstrt,kxp,kblok)
c --------------------------------------------------------------------
c copy data from field to particle partition, and copy to guard cells
c --------------------------------------------------------------------
      call cppfp1 (fx,fc,isign,scr,kstrt,nvp,nxpmx,nblok,kxp,kblok,idps)
c ----------------------------------------
c particle push and charge density update
c ----------------------------------------
c     call timera(-1,'push    ')
c -----------------------------------
c initialize charge density to zero
c -----------------------------------
      do 1170 k = 1, nblok
      nxp3 = nxp(k) + 3
      do 1160 j = 1, nxp3
c     do 1160 j = 1, nxpmx
c     if (j.le.nxp3) then
      q(j,k) = 0.
c     endif
 1160 continue
 1170 continue
      wke = 0.
c ---------------
c push particles
c ---------------
      call ppush1(part,fx,npp,noff,qtme,dt,wke,idimp,npmax,nblok,nxpmx)
c ------------------------------------------------
c move particles into appropriate spatial regions
c ------------------------------------------------
      call pmove1 (part,edges,npp,sbufr,sbufl,rbufr,rbufl,ihole,jsr,jsl,
     1jss,nx,kstrt,nvp,idimp,npmax,nblok,idps,nbmax,ntmax,ierr)
c     call pxmov1 (part,edges,npp,sbufr,sbufl,rbufr,rbufl,ihole,jsr,jsl,
c    1jss,nx,kstrt,nvp,idimp,npmax,nblok,idps,nbmax,ntmax,maskp,ierr)
c ---------------
c deposit charge
c ---------------
      call pdost1 (part,q,npp,noff,qme,idimp,npmax,nblok,nxpmx)
c     call pxost1 (part,q,tx,ty,npp,nxp,noff,nn,amx,dxu,dxl,qme,idimp,np
c    1max,nblok,nxpmx,npd)
c     call pyost1 (part,q,qt,npp,nxp,noff,qme,idimp,npmax,nblok,nxpmx,np
c    1d)
c ----------------------
c add background charge
c ----------------------
      do 1190 k = 1, nblok
      do 1180 j = 1, nxp(k)
c     do 1180 j = 1, nxpmx
c     if (j.le.nxp(k)) then
      q(j+1,k) = q(j+1,k) + qi0
c     endif
 1180 continue
 1190 continue
c     call timera(1,'push    ')
c ------------------
c energy diagnostic
c ------------------
      wt = we + wke
      if (kstrt.eq.1) write (6,993) we, wke, wt
      itime = itime + 1
      go to 500
 2000 continue
c ====================================
c
c * * * end main iteration loop * * *
c
c ====================================
      if (kstrt.eq.1) write (6,992)
      call timera(1,'total   ')
      call ppexit
      stop
      end
C =========================================================================
c
c - - - this subroutine initializes parallel processing
c
C =========================================================================
      subroutine ppinit(idproc,nvp)
c idproc = processor id
c nvp = number of real or virtual processors requested
      character*8 jbname
c     integer PVMDEFAULT, PVMRAW, INTEGER4
c     integer CMMD_partition_size, CMMD_self_address
c     integer CMMD_set_io_mode, CMMD_independent
c     parameter(PVMDEFAULT=0,PVMRAW=1,INTEGER4=3,MAXTDS=1024)
c     parameter(CMMD_independent=1)
      common /pparms/ nproc, nprec, lgrp, jbname
c     common /pvtids/ idts
c     dimension idts(MAXTDS+1)
      dimension nbuf(4)
c     save pvtids
c nproc = number of real or virtual processors obtained
c nprec = (4,8) = number of bytes per word

csm this segment is used for shared memory computers
csm     nproc = nvp
csm     idproc = 0

cipsc/pvm/pvm3
      INCLUDE 'pvm/ppinit.pvm.f'

ccm5 this segment is used for distributed memory cm-5 computers
ccm5     nprec = 4
ccm5     nproc = CMMD_partition_size()
ccm5     idproc = CMMD_self_address()
ccm5     ierr = CMMD_set_io_mode(6,CMMD_independent)

c this segment is used for distributed memory sp1 computers
c lgrp = group identifier for all groups
      nprec = 4
      call mp_environ(nproc,idproc)
      call mp_task_query(nbuf,4,3)
      lgrp = nbuf(4)
      if (nproc.ne.nvp) then
         write (6,*) ' processor number error: nvp,nproc=',nvp,nproc
         call ppexit
         stop
      endif
      return
      end

cpvm/pvm3
      INCLUDE 'pvm/ppexit.pvm.f'

cpvm3
      INCLUDE 'pvm/gittid.pvm.f'
cpvm3 
      INCLUDE 'pvm/gitpid.pvm.f'

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   kstrt = starting data block number
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 - - - for 1d code, this subroutine calculates initial particle co-ordinate
c - - - and velocity, with uniform density and maxwellian velocity with drift,
c - - - for distributed data.
C =========================================================================
      subroutine pistr1 (part,edges,npp,nps,vtx,vdx,npx,nx,idimp,npmax,n
     1blok,idps)
c part(1,n,k) = position x of particle n in partition k
c part(2,n,k) = velocity vx of particle n in partition k
c edges(1,k) = left boundary of particle partition k
c edges(2,k) = right boundary of particle partition k
c npp(k) = number of particles in partition k
c nps(k) = starting address of particles in partition k
c vtx = thermal velocity of particles in x direction
c vdx = drift velocity of particles x direction
c npx = number of particles distributed in x direction
c idimp = size of phase space = 2
c nx = system length in x direction
c npmax = maximum number of particles in each partition
c nblok = number of particle partitions.
c idps = number of partition boundaries
c ranorm = gaussian random number with zero mean and unit variance
c with spatial decomposition
      character*8 jbname
c     integer CMMD_combiner_fadd
c     parameter(CMMD_combiner_fadd=9)
c     double precision CMMD_reduce_float
      double precision ranorm
      common /pparms/ nproc, nprec, lgrp, jbname
      external s_vadd
      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)
c     do 30 j = 1, npmax
c     if ((j.ge.nps(k)).and.(j.le.npp(k))) then
      sum0 = sum0 + part(2,j,k)
c     endif
   30 continue
      sum1 = sum1 + sum0
   40 continue
cipsc this line is used for distributed memory intel computers
cipsc     call gssum(sum1,1,work1)
cpvm this line is used for distributed memory pvm/pvm3 clusters
cpvm     call psum(sum1,work1,1,1)
ccm5 this line is used for distributed memory cm-5 computers
ccm5     sum1 = CMMD_reduce_float(sum1,CMMD_combiner_fadd)

csp1 this line is used for distributed memory sp1 computers
      call mp_combine(sum1,work1,nprec,s_vadd,lgrp)
      sum1 = work1
      sum1 = sum1/float(npx) - vdx
      do 60 k = 1, nblok
      do 50 j = nps(k), npp(k)
c     do 50 j = 1, npmax
c     if ((j.ge.nps(k)).and.(j.le.npp(k))) then
      part(2,j,k) = part(2,j,k) - sum1
c     endif
   50 continue
   60 continue
      return
      end
C =========================================================================
c for 1d code, this subroutine calculates particle charge density
c using second-order spline interpolation, with periodic boundaries
c and distributed data.
c charge density is approximated by values at the nearest grid points
c q(n)=qm*(.75-dx**2), q(n+1)=.5*qm*(.5+dx)**2, q(n-1)=.5*qm*(.5-dx)**2
c where n = nearest grid point and dx = x-n
c qm = charge on particle, in units of e
c scalar version with spatial decomposition
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)
c     do 10 j = 1, npmax
c     if (j.le.npp(k)) then
      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
c     endif
   10 continue
   20 continue
      return
      end
C =========================================================================
c for 1d code, this subroutine calculates particle charge density
c using second-order spline interpolation, with periodic boundaries
c and distributed data.
c charge density is approximated by values at the nearest grid points
c q(n)=qm*(.75-dx**2), q(n+1)=.5*qm*(.5+dx)**2, q(n-1)=.5*qm*(.5-dx)**2
c where n = nearest grid point and dx = x-n
c qm = charge on particle, in units of e
c idimp = size of phase space = 2
c vector version with spatial decomposition
C =========================================================================
      subroutine pvost1 (part,q,npp,noff,nn,nu,nl,amx,dxu,dxl,qm,idimp,n
     1pmax,nblok,nxpmx,npd)
      dimension part(idimp,npmax,nblok), q(nxpmx,nblok)
      dimension npp(nblok), noff(nblok)
      dimension nn(npd,nblok), nu(npd,nblok), nl(npd,nblok)
      dimension amx(npd,nblok), dxu(npd,nblok), dxl(npd,nblok)
c ------------------------------
c find maximum number of blocks
c ------------------------------
      if (npmax.gt.npd) then
         ipmax = float(npmax - 1)/float(npd) + 1.
      else
         ipmax = 1
      endif
c --------------
c parallel loop
c --------------
      do 40 k = 1, nblok
      nnoff = noff(k) - 2
      npb = npd
      if (npp(k).gt.npd) then
         ipp = float(npp(k) - 1)/float(npd) + 1.
      else
         ipp = 1
      endif
c ------------------------------------
c outer loop over blocks of particles
c ------------------------------------
      do 30 j = 1, ipp
c do 30 j = 1, ipmax
      jb = (j - 1)*npd
      if (j.ge.ipp) npb = npp(k) - (ipp - 1)*npd
c ----------------------------
c find interpolation weights
c ----------------------------
      do 10 i = 1, npb
c     do 10 i = 1, npd
c     if ((j.le.ipp).and.(i.le.npb)) then
         nn(i,k) = part(1,i+jb,k) + .5
         dxl(i,k) = part(1,i+jb,k) - float(nn(i,k))
         nn(i,k) = nn(i,k) - nnoff
         amx(i,k) = qm*(.75 - dxl(i,k)*dxl(i,k))
         nu(i,k) = nn(i,k) + 1
         dxu(i,k) = .5*qm*(.5 + dxl(i,k))**2
         nl(i,k) = nn(i,k) - 1
         dxl(i,k) = .5*qm*(.5 - dxl(i,k))**2
c     endif
   10 continue
c deposit charge
      do 20 i = 1, npb
c do 20 i = 1, npd
c if ((j.le.ipp).and.(i.le.npb)) then
         q(nn(i,k),k) = q(nn(i,k),k) + amx(i,k)
         q(nu(i,k),k) = q(nu(i,k),k) + dxu(i,k)
         q(nl(i,k),k) = q(nl(i,k),k) + dxl(i,k)
c     endif
   20 continue
   30 continue
   40 continue
      return
      end
C =========================================================================
c for 1d code, this subroutine calculates particle charge density
c using second-order spline interpolation, with periodic boundaries
c and distributed data.
c charge density is approximated by values at the nearest grid points
c q(n)=qm*(.75-dx**2), q(n+1)=.5*qm*(.5+dx)**2, q(n-1)=.5*qm*(.5-dx)**2
c where n = nearest grid point and dx = x-n
c qu, ql = scratch arrays
c qm = charge on particle, in units of e
c idimp = size of phase space = 2
c npd = size of scratch buffers for vectorized charge deposition
c optimized version with spatial decomposition
C =========================================================================
      subroutine pxost1 (part,q,qu,ql,npp,nxp,noff,nn,amx,dxu,dxl,qm,idi
     1mp,npmax,nblok,nxpmx,npd)
      dimension part(idimp,npmax,nblok), q(nxpmx,nblok)
      dimension qu(nxpmx,nblok), ql(nxpmx,nblok)
      dimension npp(nblok), nxp(nblok), noff(nblok)
      dimension nn(npd,nblok)
      dimension amx(npd,nblok), dxu(npd,nblok), dxl(npd,nblok)
c find maximum number of blocks
      if (npmax.gt.npd) then
         ipmax = float(npmax - 1)/float(npd) + 1.
      else
         ipmax = 1
      endif
c --------------
c parallel loop
c --------------
      do 60 k = 1, nblok
c zero out scratch arrays
         nxp3 = nxp(k) + 3
         do 10 j = 1, nxp3
c do 10 j = 1, nxpmx
c if (j.le.nxp3) then
            qu(j,k) = 0.
            ql(j,k) = 0.
c endif
   10    continue
         nnoff = noff(k) - 2
         npb = npd
         if (npp(k).gt.npd) then
            ipp = float(npp(k) - 1)/float(npd) + 1.
         else
            ipp = 1
         endif
c ------------------------------------
c outer loop over blocks of particles
c ------------------------------------
         do 40 j = 1, ipp
c do 40 j = 1, ipmax
            jb = (j - 1)*npd
            if (j.ge.ipp) npb = npp(k) - (ipp - 1)*npd
c ---------------------------
c find interpolation weights
c ---------------------------
            do 20 i = 1, npb
c do 20 i = 1, npd
c if ((j.le.ipp).and.(i.le.npb)) then
               nn(i,k) = part(1,i+jb,k) + .5
              dxl(i,k) = part(1,i+jb,k) - float(nn(i,k))
              nn(i,k) = nn(i,k) - nnoff
              amx(i,k) = qm*(.75 - dxl(i,k)*dxl(i,k))
              dxu(i,k) = .5*qm*(.5 + dxl(i,k))**2
              dxl(i,k) = .5*qm*(.5 - dxl(i,k))**2
c endif
   20       continue
c deposit charge
              do 30 i = 1, npb
c do 30 i = 1, npd
c if ((j.le.ipp).and.(i.le.npb)) then
                 q(nn(i,k),k) = q(nn(i,k),k) + amx(i,k)
                 qu(nn(i,k),k) = qu(nn(i,k),k) + dxu(i,k)
                 ql(nn(i,k),k) = ql(nn(i,k),k) + dxl(i,k)
c endif
   30         continue
   40    continue
c -----------------
c add partial sums
c -----------------
      nxp3 = nxp(k) + 3
      nxp2 = nxp(k) + 2
      do 50 j = 2, nxp2
c do 50 j = 2, nxpmx
c if (j.le.nxp2) then
         q(j,k) = q(j,k) + qu(j-1,k) + ql(j+1,k)
c endif
   50 continue
      q(1,k) = ql(2,k)
      q(nxp3,k) = qu(nxp2,k)
   60 continue
      return
      end
C =========================================================================
c for 1d code, this subroutine calculates particle charge density
c using second-order spline interpolation, with periodic boundaries
c and distributed data, using t. yabe's memory duplication scheme, as in
c a. nishiguchi, s. orii, and t. yabe, j. comput. phys. 61, 519 (1985).
c charge density is approximated by values at the nearest grid points
c q(n)=qm*(.75-dx**2), q(n+1)=.5*qm*(.5+dx)**2, q(n-1)=.5*qm*(.5-dx)**2
c where n = nearest grid point and dx = x-n
c qm = charge on particle, in units of e
c idimp = size of phase space = 2
c vector version with spatial decomposition
C =========================================================================
      subroutine pyost1 (part,q,qt,npp,nxp,noff,qm,idimp,npmax,nblok,nxp
     1mx,npd)
      dimension part(idimp,npmax,nblok), q(nxpmx,nblok)
      dimension npp(nblok), nxp(nblok), noff(nblok)
      dimension qt(npd,nxpmx,nblok)
c -------------------------------
c find maximum number of blocks
c -------------------------------
      if (npmax.gt.npd) then
         ipmax = float(npmax - 1)/float(npd) + 1.
      else
         ipmax = 1
      endif
c -------------------------------
c parallel loop
c -------------------------------
      do 70 k = 1, nblok
      nnoff = noff(k) - 2
      npb = npd
      if (npp(k).gt.npd) then
         ipp = float(npp(k) - 1)/float(npd) + 1.
      else
         ipp = 1
      endif
c -------------------------------
c zero out duplicate charge arrays
c -------------------------------
      nxp3 = nxp(k) + 3
      do 20 l = 1, nxp3
c do 20 l = 1, nxpmx
      do 10 j = 1, npd
c if (l.le.nxp3) then
         qt(j,l,k) = 0.
c endif
   10 continue
   20 continue
c -------------------------------
c outer loop over blocks of particles
c -------------------------------
      do 40 j = 1, ipp
c do 40 j = 1, ipmax
         jb = (j - 1)*npd
         if (j.ge.ipp) npb = npp(k) - (ipp - 1)*npd
c -----------------------------------------------
c find interpolation weights and deposit charge
c -----------------------------------------------
      do 30 i = 1, npb
c do 30 i = 1, npd
c if ((j.le.ipp).and.(i.le.npb)) then
         nn = part(1,i+jb,k) + .5
         dx = part(1,i+jb,k) - float(nn)
         nn = nn - nnoff
         qt(i,nn,k) = qt(i,nn,k) + qm*(.75 - dx*dx)
         nn = nn + 1
         qt(i,nn,k) = qt(i,nn,k) + .5*qm*(.5 + dx)**2
         nn = nn - 2
         qt(i,nn,k) = qt(i,nn,k) + .5*qm*(.5 - dx)**2
c endif
   30 continue
   40 continue
c -------------------------------
c add up duplicate charge arrays
c -------------------------------
      do 60 l = 1, nxp3
c do 60 l = 1, nxpmx
         sum1 = 0.
         do 50 j = 1, npd
c if (l.le.nxp3) then
            sum1 = sum1 + qt(j,l,k)
c endif
   50    continue
         q(l,k) = q(l,k) + sum1
   60 continue
   70 continue
      return
      end
C =========================================================================
c for 1d code, this subroutine calculates particle charge density
c using second-order spline interpolation, with periodic boundaries
c and distributed data.
c charge density is approximated by values at the nearest grid points
c q(n)=qm*(.75-dx**2), q(n+1)=.5*qm*(.5+dx)**2, q(n-1)=.5*qm*(.5-dx)**2
c where n = nearest grid point and dx = x-n
c qu, ql = scratch arrays
c nn = scratch address array for vectorized charge deposition
c qm = charge on particle, in units of e
c idimp = size of phase space = 2
c special version with spatial decomposition
C =========================================================================
      subroutine pzost1 (part,q,qu,ql,npp,nxp,noff,nn,amx,dxu,dxl,qm,idi
     1mp,npmax,nblok,nxpmx,npd)
      dimension part(idimp,npmax,nblok), q(nxpmx,nblok)
      dimension qu(nxpmx,nblok), ql(nxpmx,nblok)
      dimension npp(nblok), nxp(nblok), noff(nblok)
      dimension nn(npd,nblok)
      dimension amx(npd,nblok), dxu(npd,nblok), dxl(npd,nblok)
      nxp1 = nxpmx - 1
c find maximum number of blocks
      if (npmax.gt.npd) then
         ipmax = float(npmax - 1)/float(npd) + 1.
      else
         ipmax = 1
      endif
c ----------------
c parallel loop
c ----------------
      do 70 k = 1, nblok
c ------------------------
c zero out scratch arrays
c ------------------------
      nxp3 = nxp(k) + 3
      do 10 j = 1, nxp3
c     do 10 j = 1, nxpmx
c     if (j.le.nxp3) then
      qu(j,k) = 0.
      ql(j,k) = 0.
c     endif
   10 continue
      nnoff = noff(k) - 2
      npb = npd
      if (npp(k).gt.npd) then
         ipp = float(npp(k) - 1)/float(npd) + 1.
      else
         ipp = 1
      endif
c ------------------------------------
c outer loop over blocks of particles
c ------------------------------------
      do 50 j = 1, ipp
c     do 50 j = 1, ipmax
      jb = (j - 1)*npd
      if (j.ge.ipp) npb = npp(k) - (ipp - 1)*npd
c ------------------------------------
c find interpolation weights
c ------------------------------------
      do 20 i = 1, npb
c do 20 i = 1, npd
c if ((j.le.ipp).and.(i.le.npb)) then
      nn(i,k) = part(1,i+jb,k) + .5
      dxl(i,k) = part(1,i+jb,k) - float(nn(i,k))
      nn(i,k) = nn(i,k) - nnoff
      amx(i,k) = qm*(.75 - dxl(i,k)*dxl(i,k))
      dxu(i,k) = .5*qm*(.5 + dxl(i,k))**2
      dxl(i,k) = .5*qm*(.5 - dxl(i,k))**2
c endif
   20 continue
c ------------------------------------
c deposit charge
c ------------------------------------
      do 40 l = 2, nxp1
      do 30 i = 1, npb
c do 30 i = 1, npd
c if ((j.le.ipp).and.(i.le.npb)) then
      if (nn(i,k).eq.l) then
         q(l,k) = q(l,k) + amx(i,k)
         qu(l,k) = qu(l,k) + dxu(i,k)
         ql(l,k) = ql(l,k) + dxl(i,k)
      endif
c endif
   30 continue
   40 continue
   50 continue
c ------------------------------------
c add partial sums
c ------------------------------------
      nxp3 = nxp(k) + 3
      nxp2 = nxp(k) + 2
      do 60 j = 2, nxp2
c do 60 j = 2, nxpmx
c if (j.le.nxp2) then
         q(j,k) = q(j,k) + qu(j-1,k) + ql(j+1,k)
c endif
   60 continue
      q(1,k) = ql(2,k)
      q(nxp3,k) = qu(nxp2,k)
   70 continue
      return
      end
C =========================================================================
c thus subroutine copies data between particle and field partitions,
c where the field and particle partitions are assumed to be the same.
c in going from particle to field partitions, data from guard cells
c is added to the field partition.  in going from field to particle
c partitions, data is copied to the guard cells.
c f(j,k) = real data for grid j in particle partition k.  the number of
c grids per partition is uniform and includes three extra guard cells.
c fc(j,k) = complex data for grid j in field partition k.  the number of
c grids per partition is uniform and maximum number of partitions is nxh
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 scr(idps,k) = scratch array for particle partition k
c kstrt = starting data block number
c kxp = number of complex grids in each field partition.
C =========================================================================
      subroutine cppfp1 (f,fc,isign,scr,kstrt,nvp,nxpmx,nblok,kxp,kblok,
     1idps)
      character*8 jbname
c     integer PVMRAW, REAL4, COMPLEX8
c     integer CMMD_send_async, CMMD_receive_block
c     parameter(PVMRAW=1,REAL4=4,COMPLEX8=5)
      common /pparms/ nproc, nprec, lgrp, jbname
      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
c ---------------------------------------
c add guard cells
c ---------------------------------------
csm this loop is used for shared memory computers
csm     do 50 k = 1, kblok
csm     kl = k + ks
csm     if (kl.lt.1) then
csm        scr(1,k) = f(2*kxp+2,kl+nvp)
csm        scr(2,k) = f(2*kxp+3,kl+nvp)
csm     else
csm        scr(1,k) = f(2*kxp+2,kl)
csm        scr(2,k) = f(2*kxp+3,kl)
csm     endif
csm     fc(1,k) = fc(1,k) + cmplx(scr(1,k),scr(2,k))
csm     kr = k + ks + 2
csm     if (kr.gt.nvp) then
csm        scr(1,k) = f(1,kr-nvp)
csm     else
csm        scr(1,k) = f(1,kr)
csm     endif
csm     fc(kxp,k) = fc(kxp,k) + cmplx(0.,scr(1,k))
csm  50 continue

cipsc/pvm/pvm3
c     INCLUDE 'pvm/cppfp1.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
      do 50 k = 1, kblok
      kr = k + ks + 2
      kl = k + ks
      if (kl.lt.1) then
         call mp_recv(scr,2*nprec,kl+nvp-1,moff+1,msid)
      else
         call mp_recv(scr,2*nprec,kl-1,moff+1,msid)
      endif
      if (kr.gt.nvp) then
         call mp_bsend(f(2*kxp+2,k),2*nprec,kr-nvp-1,moff+1)
      else
         call mp_bsend(f(2*kxp+2,k),2*nprec,kr-1,moff+1)
      endif
      call mp_wait(msid,nbytes)
      fc(1,k) = fc(1,k) + cmplx(scr(1,k),scr(2,k))
      if (kr.gt.nvp) then
         call mp_recv(scr,nprec,kr-nvp-1,moff+2,msid)
      else
         call mp_recv(scr,nprec,kr-1,moff+2,msid)
      endif
      if (kl.lt.1) then
         call mp_bsend(f(1,k),nprec,kl+nvp-1,moff+2)
      else
         call mp_bsend(f(1,k),nprec,kl-1,moff+2)
      endif
      call mp_wait(msid,nbytes)
      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
c --------------------
c copy to guard cells
c -------------------- 

csm this loop is used for shared memory computers
csm     do 100 k = 1, kblok
csm     kl = k + ks
csm     if (kl.lt.1) then
csm        f(1,k) = aimag(fc(kxp,kl+nvp))
csm     else
csm        f(1,k) = aimag(fc(kxp,kl))
csm     endif
csm     kr = k + ks + 2
csm     if (kr.gt.nvp) then
csm        f(2*kxp+2,k) = real(fc(1,kr-nvp))
csm        f(2*kxp+3,k) = aimag(fc(1,kr-nvp))
csm     else
csm        f(2*kxp+2,k) = real(fc(1,kr))
csm        f(2*kxp+3,k) = aimag(fc(1,kr))
csm     endif
csm  100 continue

cipsc/pvm/pvm3
c     INCLUDE 'pvm/cppfp1_2.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
      do 100 k = 1, kblok
      kr = k + ks + 2
      kl = k + ks
      scr(1,k) = aimag(fc(kxp,k))
      if (kl.lt.1) then
         call mp_recv(f(1,k),nprec,kl+nvp-1,moff+3,msid)
      else
         call mp_recv(f(1,k),nprec,kl-1,moff+3,msid)
      endif
      if (kr.gt.nvp) then
         call mp_bsend(scr,nprec,kr-nvp-1,moff+3)
      else
         call mp_bsend(scr,nprec,kr-1,moff+3)
      endif
      call mp_wait(msid,nbytes)
      if (kr.gt.nvp) then
         call mp_recv(f(2*kxp+2,k),2*nprec,kr-nvp-1,moff+4,msid)
      else
         call mp_recv(f(2*kxp+2,k),2*nprec,kr-1,moff+4,msid)
      endif
      if (kl.lt.1) then
        call mp_bsend(fc(1,k),2*nprec,kl+nvp-1,moff+4)
      else
         call mp_bsend(fc(1,k),2*nprec,kl-1,moff+4)
      endif
      call mp_wait(msid,nbytes)
  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, for distributed data.
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 kstrt = starting data block number
c if isign = 0, form factor array is prepared
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 ax = half-width of particle in x direction
c affp = normalization constant = nx/np, where np = number of particles
c electric field energy is also calculated, using
c wp = nx*sum((affp/k**2)*|q(k)*s(k)|**2)
c nx = system length in x direction
C =========================================================================
      subroutine ppois1 (q,fx,isign,ffc,ax,affp,wp,nx,kstrt,kxp,kblok)
      character*8 jbname
ccm5     integer CMMD_combiner_fadd
ccm5     parameter(CMMD_combiner_fadd=9)
ccm5     double precision CMMD_reduce_float
      common /pparms/ nproc, nprec, lgrp, jbname
      complex q, fx, ffc
      external s_vadd
      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

cipsc this line is used for distributed memory intel computers
cipsc     call gssum(wp,1,work1)
cpvm this line is used for distributed memory pvm/pvm3 clusters
cpvm     call psum(wp,work1,1,1)
ccm5 this line is used for distributed memory cm-5 computers
ccm5     wp = CMMD_reduce_float(wp,CMMD_combiner_fadd)

c this line is used for distributed memory sp1 computers
      call mp_combine(wp,work1,nprec,s_vadd,lgrp)
      wp = work1
      wp = float(nx)*wp
      return
      end
C =========================================================================
c for 1d code, 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, for distributed data.
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 and x(t+dt) = x(t) + v(t+dt/2)*dt
c fx(x(t)) is approximated by interpolation from the nearest grid points
c fx(x) = (.75-dx**2)*fx(n)+.5*(fx(n+1)*(.5+dx)**2+fx(n-1)*(.5-dx)**2)
c where n = nearest grid point and dx = x-n
c qtm = particle charge/mass ratio times dt
c kinetic energy/mass at time t is also calculated, using
c ek = .125*sum((v(t+dt/2)+v(t-dt/2))**2)
c idimp = size of phase space = 2
c scalar version with spatial decomposition
C =========================================================================
      subroutine ppush1(part,fx,npp,noff,qtm,dt,ek,idimp,npmax,nblok,nxp
     1mx)
      character*8 jbname
ccm5     integer CMMD_combiner_fadd
ccm5     parameter(CMMD_combiner_fadd=9)
ccm5     double precision CMMD_reduce_float
      common /pparms/ nproc, nprec, lgrp, jbname
      external s_vadd
      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 do 10 j = 1, npmax
c if (j.le.npp(k)) then
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 + fx(nn-
     11,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
c     endif
   10 continue
   20 continue
cipsc this line is used for distributed memory intel computers
cipsc     call gssum(sum1,1,work1)
cpvm this line is used for distributed memory pvm/pvm3 clusters
cpvm     call psum(sum1,work1,1,1)
ccm5 this line is used for distributed memory cm-5 computers
ccm5     sum1 = CMMD_reduce_float(sum1,CMMD_combiner_fadd)

csp1 this line is used for distributed memory sp1 computers
      call mp_combine(sum1,work1,nprec,s_vadd,lgrp)
      sum1 = work1

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


c vector version with spatial decomposition :: push routine
c      INCLUDE 'vector/pvpsh1.f'

c optimized version with spatial decomposition :: push routine
c      INCLUDE 'vector/pxpsh1.f'

C =========================================================================
c this subroutine moves particles into appropriate spatial regions
c sbufl = buffer for particles being sent to left processor
c sbufr = buffer for particles being sent to right processor
c nx = system length in x direction
c kstrt = starting data block number
c idimp = size of phase space = 2
c idps = number of partition boundaries
c nbmax =  size of buffers for passing particles between processors
c ntmax =  size of hole array for particles leaving processors
c ierr = (0,1) = (no,yes) error condition exists
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)
      character*8 jbname
c     integer PVMRAW, INTEGER4, REAL4
c     parameter(PVMRAW=1,INTEGER4=3,REAL4=4)
ccm5     integer CMMD_send_async, CMMD_receive_block, CMMD_reduce_int
ccm5     integer CMMD_combiner_add
ccm5     parameter(CMMD_combiner_add=0)
      common /pparms/ nproc, nprec, lgrp, jbname
      external i_vadd
      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

cipsc this line is used for distributed memory intel computers
cipsc     call gisum(npr,1,iwork1)
cpvm this line is used for distributed memory pvm/pvm3 clusters
cpvm     call pisum(npr,iwork1,1,1)
ccm5 this line is used for distributed memory cm-5 computers
ccm5     npr = CMMD_reduce_int(npr,CMMD_combiner_add)

csp1 this segment is used for distributed memory sp1 computers
      call mp_combine(npr,iwork1,4,i_vadd,lgrp)
      npr = iwork1
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)
c     do 20 j = 1, npmax
c     if (j.le.npp(k)) then
      xt = part(1,j,k)
c ------------------------------
c particles going to the right
c ------------------------------
      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
c ----------------------------
c particles going to the left
c ----------------------------
      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
c     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

cipsc this line is used for distributed memory intel computers
cipsc     call gisum(ibflg,1,iwork1)
cpvm this line is used for distributed memory pvm/pvm3 clusters
cpvm     call pisum(ibflg,iwork1,1,1)
ccm5 this line is used for distributed memory cm-5 computers
ccm5     ibflg = CMMD_reduce_int(ibflg,CMMD_combiner_add)

csp1 this segment is used for distributed memory sp1 computers
      call mp_combine(ibflg,iwork1,4,i_vadd,lgrp)
      ibflg = iwork1

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

csm this segment is used for shared memory computers
csm     kl = k + ks
csm     if (kl.lt.1) then
csm        jsl(2,k) = jsr(1,kl+nvp)
csm     else
csm        jsl(2,k) = jsr(1,kl)
csm     endif

cipsc/pvm/pvm3/cm5

c      INCLUDE 'pvm/pmove1.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
      kr = k + ks + 2
      kl = k + ks
      if (kl.lt.1) then
         call mp_recv(jsl(2,k),4,kl+nvp-1,iter-3,msid)
      else
         call mp_recv(jsl(2,k),4,kl-1,iter-3,msid)
      endif
      if (kr.gt.nvp) then
         call mp_bsend(jsr(1,k),4,kr-nvp-1,iter-3)
      else
         call mp_bsend(jsr(1,k),4,kr-1,iter-3)
      endif
      call mp_wait(msid,nbytes)

csm this loop is used for shared memory computers
csm     do 60 j = 1, jsl(2,k)
csm     do 60 j = 1, nbmax
csm     if (j.le.jsl(2,k)) then
csm     if (kl.lt.1) then
csm        rbufl(1,j,k) = sbufr(1,j,kl+nvp)
csm        rbufl(2,j,k) = sbufr(2,j,kl+nvp)
csm     else
csm        rbufl(1,j,k) = sbufr(1,j,kl)
csm        rbufl(2,j,k) = sbufr(2,j,kl)
csm     endif
csm     endif
csm  60 continue

cipsc/pvm/pvm3/cm5
c    INCLUDE 'pvm/pmove1_2.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
      if (kl.lt.1) then
         call mp_recv(rbufl,idimp*nprec*jsl(2,k),kl+nvp-1,iter-2,msid)
      else
         call mp_recv(rbufl,idimp*nprec*jsl(2,k),kl-1,iter-2,msid)
      endif
      if (kr.gt.nvp) then
         call mp_bsend(sbufr,idimp*nprec*jsr(1,k),kr-nvp-1,iter-2)
      else
         call mp_bsend(sbufr,idimp*nprec*jsr(1,k),kr-1,iter-2)
      endif
      call mp_wait(msid,nbytes)
c --------------------------
c get particles from right
c --------------------------
csm this segment is used for shared memory computers
csm     kr = k + ks + 2
csm     if (kr.gt.nvp) then
csm        jsr(2,k) = jsl(1,kr-nvp)
csm     else
csm        jsr(2,k) = jsl(1,kr)
csm     endif

cipsc/pvm/pvm3/cm5
c   INCLUDE 'pvm/pmove1_3.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
      if (kr.gt.nvp) then
         call mp_recv(jsr(2,k),4,kr-nvp-1,iter-1,msid)
      else
         call mp_recv(jsr(2,k),4,kr-1,iter-1,msid)
      endif
      if (kl.lt.1) then
         call mp_bsend(jsl(1,k),4,kl+nvp-1,iter-1)
      else
         call mp_bsend(jsl(1,k),4,kl-1,iter-1)
      endif
      call mp_wait(msid,nbytes)
c this loop is used for shared memory computers
c     do 70 j = 1, jsr(2,k)
c     do 70 j = 1, nbmax
c     if (j.le.jsr(2,k)) then
c     if (kr.gt.nvp) then
c        rbufr(1,j,k) = sbufl(1,j,kr-nvp)
c        rbufr(2,j,k) = sbufl(2,j,kr-nvp)
c     else
c        rbufr(1,j,k) = sbufl(1,j,kr)
c        rbufr(2,j,k) = sbufl(2,j,kr)
c     endif
c     endif
c  70 continue

cipsc/pvm/pvm3/cm5
c      INCLUDE 'pvm/pmove1_4.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
      if ((k+ks+2).gt.nvp) then
         call mp_recv(rbufr,idimp*nprec*jsr(2,k),kr-nvp-1,iter,msid)
      else
         call mp_recv(rbufr,idimp*nprec*jsr(2,k),kr-1,iter,msid)
      endif
      if ((k+ks).lt.1) then
         call mp_bsend(sbufl,idimp*nprec*jsl(1,k),kl+nvp-1,iter)
      else
         call mp_bsend(sbufl,idimp*nprec*jsl(1,k),kl-1,iter)
      endif
      call mp_wait(msid,nbytes)
   80 continue
c ------------------------------------------- 
c check if particles must be passed further
c ------------------------------------------- 
      nps = 0
      do 110 k = 1, nblok
c ------------------------------------------------------------------ 
c check if any particles coming from right must go further to left
c ------------------------------------------------------------------ 
      jsl(1,k) = 0
      do 90 j = 1, jsr(2,k)
c     do 90 j = 1, nbmax
c     if (j.le.jsr(2,k)) then
      if (rbufr(1,j,k).lt.edges(1,k)) jsl(1,k) = 1
c     endif
   90 continue
c ------------------------------------------------------------------ 
c check if any particles coming from left must go further to right
c ------------------------------------------------------------------ 
      jsr(1,k) = 0
      do 100 j = 1, jsl(2,k)
c     do 100 j = 1, nbmax
c     if (j.le.jsl(2,k)) then
      if (rbufl(1,j,k).ge.edges(2,k)) jsr(1,k) = 1
c     endif
  100 continue
      nps = nps + (jsl(1,k) + jsr(1,k))
  110 continue

cipsc this line is used for distributed memory intel computers
cipsc     call gisum(nps,1,iwork1)
cpvm this line is used for distributed memory pvm/pvm3 clusters
cpvm     call pisum(nps,iwork1,1,1)
ccm5 this line is used for distributed memory cm-5 computers
ccm5     nps = CMMD_reduce_int(nps,CMMD_combiner_add)

csp1 this segment is used for distributed memory sp1 computers
      call mp_combine(nps,iwork1,4,i_vadd,lgrp)
      nps = iwork1
      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
c -------------------------------------------
c first check particles coming from right
c -------------------------------------------
      jsl(1,k) = 0
      jss(2,k) = 0
      do 120 j = 1, jsr(2,k)
c     do 120 j = 1, nbmax
c     if (j.le.jsr(2,k)) then
      xt = rbufr(1,j,k)
c -------------------------------------------
c particles going to the left
c -------------------------------------------
      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)
c -------------------------------------------
c particles staying here
c -------------------------------------------
      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
c     endif
  120 continue
      jsr(2,k) = jss(2,k)
c -------------------------------------------
c next check particles coming from left
c -------------------------------------------
      jsr(1,k) = 0
      jss(2,k) = 0
      do 130 j = 1, jsl(2,k)
c     do 130 j = 1, nbmax
c     if (j.le.jsl(2,k)) then
      xt = rbufl(1,j,k)
c -------------------------------------------
c particles going to the right
c -------------------------------------------
      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)
c -------------------------------------------
c particles staying here
c -------------------------------------------
      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
c     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

cipsc this line is used for distributed memory intel computers
cipsc     call gisum(ierr,1,iwork1)
cpvm this line is used for distributed memory pvm/pvm3 clusters
cpvm     call pisum(ierr,iwork1,1,1)
ccm5 this line is used for distributed memory cm-5 computers
ccm5     ierr = CMMD_reduce_int(ierr,CMMD_combiner_add)

csp1 this segment is used for distributed memory sp1 computers
      call mp_combine(ierr,iwork1,4,i_vadd,lgrp)
      ierr = iwork1
      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 distribute particles coming from left into holes
c ------------------------------------------------------
      jss(2,k) = min0(jss(1,k),jsl(2,k))
      do 170 j = 1, jss(2,k)
c     do 170 j = 1, nbmax
c     if (j.le.jss(2,k)) then
      part(1,ihole(j,k),k) = rbufl(1,j,k)
      part(2,ihole(j,k),k) = rbufl(2,j,k)
c     endif
  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)
c     do 180 j = 1, nbmax
c     if (j.le.jss(2,k)) then
c ------------------------------------------------------
c no more particles coming from left
c distribute particles coming from right into holes
c ------------------------------------------------------
      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
c ------------------------------------------------------
c no more holes
c distribute remaining particles from left into bottom
c ------------------------------------------------------
         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
c     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     do 190 j = 1, ntmax
c     if (j.le.jsr(2,k)) then
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
c     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

cipsc this line is used for distributed memory intel computers
cipsc     call gisum(nps,1,iwork1)
cpvm this line is used for distributed memory pvm/pvm3 clusters
cpvm     call pisum(nps,iwork1,1,1)
ccm5 this line is used for distributed memory cm-5 computers
ccm5     nps = CMMD_reduce_int(nps,CMMD_combiner_add)

csp1 this segment is used for distributed memory sp1 computers
      call mp_combine(nps,iwork1,4,i_vadd,lgrp)
      nps = iwork1
      if (npr.ne.nps) write (6,*) 'number error, np-old,np-new=',npr,nps
      return
      end

c optimized version for vector processor :: move
c     INCLUDE 'vector/pxmov1.f'

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 if isign = 0, the fft tables are prepared
c if isign = -1, an inverse fourier transform is performed
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 kxp = number of data values per block
c kstrt = starting data block number
c mixup = array of bit reversed addresses
c sct = sine/cosine table
c the real data is stored in a complex array of length nx/2,
c with the odd/even points stored in the real/imaginary parts.
c the real part of fourier mode nx/2 is stored in the imaginary part
c of mode zero.
c written by viktor k. decyk, ucla
c parallel version
C =========================================================================
      subroutine pfft1r(f,g,isign,mixup,sct,indx,kstrt,kxp,kblok)
      character*8 jbname
c     integer PVMRAW, COMPLEX8
ccm5     integer CMMD_send_async, CMMD_receive_block
c     parameter(PVMRAW=1,COMPLEX8=5)
      common /pparms/ nproc, nprec, lgrp, jbname
      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
c ---------------------------
c copy local data
c ---------------------------
         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
c ---------------------------
c distributed calculation
c ---------------------------
      else
c ---------------------------
c copy distributed data
c ---------------------------
         do 130 k = 1, kblok
         kb = k + ks
         lb = kb/kxs
         kb = kb + 1
csm this loop is used for shared memory computers
csm        do 120 j = 1, kxp
csm        if (lb.eq.(2*(lb/2))) then
csm           f(j,k) = g(j,kb+kxs)
csm        else
csm           f(j,k) = g(j,kb-kxs)
csm        endif
csm 120    continue

cipsc/pvm/pvm3/cm5
c   INCLUDE 'pvm/pfft1r.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
         if (lb.eq.(2*(lb/2))) then
            call mp_recv(f,2*nprec*kxp,kb+kxs-1,l+kxp,msid)
            call mp_bsend(g,2*nprec*kxp,kb+kxs-1,l+kxp)
         else
            call mp_recv(f,2*nprec*kxp,kb-kxs-1,l+kxp,msid)
            call mp_bsend(g,2*nprec*kxp,kb-kxs-1,l+kxp)
         endif
         call mp_wait(msid,nbytes)
  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 this segment is used for distributed memory sp1 computers
c unscramble coefficients and normalize result
c copy distributed data
      do 180 k = 1, kblok
      kb = k + ks
      kk = kbs - kb
csm this loop is used for shared memory computers
csm     do 170 j = 2, kxp
csm     f(j,k) = g(j,kk)
csm 170 continue
c     if (kb.gt.0) f(1,k) = g(1,kk+1)

cipsc/pvm/pvm3/cm5 
c       INCLUDE 'pvm/pfft1r_2.pvm.f'

c this segment is used for distributed memory sp1 computers
      if (kxp.gt.1) then
         call mp_recv(f(2,1),2*nprec*(kxp-1),kk-1,indx+kxp,msid)
         call mp_bsend(g(2,1),2*nprec*(kxp-1),kk-1,indx+kxp)
         call mp_wait(msid,nbytes)
      endif
      if (kb.gt.0) then
         call mp_recv(f,2*nprec,kk,indx+kxp+1,msid)
         call mp_bsend(g,2*nprec,kk,indx+kxp+1)
         call mp_wait(msid,nbytes)
      endif
  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))) - (g(j,k) - conjg(f(jj,k))
     1)*cmplx(sin(arg),cos(arg)))
  190 continue
      if (kb.eq.0) then
         g(1,k) = 2.*ani*cmplx(real(g(1,k)) + aimag(g(1,k)),real(g(1,k))
     1 - aimag(g(1,k)))
      else
         arg = dnx*float(1 + koff)
         g(1,k) = ani*((g(1,k) + conjg(f(1,k))) - (g(1,k) - conjg(f(1,k)
     1))*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 copy distributed data
c ------------------
  220 do 240 k = 1, kblok
      kb = k + ks
      kk = kbs - kb
csm this loop is used for shared memory computers
csm     do 230 j = 2, kxp
csm     g(j,k) = f(j,kk)
csm 230 continue
csm     if (kb.gt.0) g(1,k) = f(1,kk+1)

cipsc/pvm/pvm3/cm5
c     INCLUDE 'pvm/pfft1r_3.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
      if (kxp.gt.1) then
         call mp_recv(g(2,1),2*nprec*(kxp-1),kk-1,indx+kxp,msid)
         call mp_bsend(f(2,1),2*nprec*(kxp-1),kk-1,indx+kxp)
         call mp_wait(msid,nbytes)
      endif
      if (kb.gt.0) then
         call mp_recv(g,2*nprec,kk,indx+kxp+1,msid)
         call mp_bsend(f,2*nprec,kk,indx+kxp+1)
         call mp_wait(msid,nbytes)
      endif
  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
c -------------------------
c distributed calculation
c -------------------------
         else
c copy distributed data
            do 340 k = 1, kblok
               kb = k + ks
               lb = kb/kxs
               kb = kb + 1
               lb = lb - 2*(lb/2)
csm this loop is used for shared memory computers
csm        do 330 j = 1, kxp
csm        if (lb.eq.0) then
csm           g(j,k) = f(j,kb+kxs)
csm        else
csm           g(j,k) = f(j,kb-kxs)
csm        endif
csm 330    continue

cipsc/pvm/pvm3/cm5
c        INCLUDE 'pvm/pfft1r_4.pvm.f'

csp1 this segment is used for distributed memory sp1 computers
         if (lb.eq.0) then
            call mp_recv(g,2*nprec*kxp,kb+kxs-1,l+kxp,msid)
            call mp_bsend(f,2*nprec*kxp,kb+kxs-1,l+kxp)
         else
            call mp_recv(g,2*nprec*kxp,kb-kxs-1,l+kxp,msid)
            call mp_bsend(f,2*nprec*kxp,kb-kxs-1,l+kxp)
         endif
         call mp_wait(msid,nbytes)
  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 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 indx = power of 2 which determines length of data, nx = 2**indx
c kstrt = starting data block number
c kblok = number of data blocks
c kxp = number of data values per block
c written by viktor k. decyk, ucla
C =========================================================================
      subroutine bitrv1(f,g,mixup,isign,indx,kstrt,kxp,kblok)
      character*8 jbname
c     integer PVMRAW, COMPLEX8
c     parameter(PVMRAW=1,COMPLEX8=5)
ccm5     integer CMMD_send_async, CMMD_receive_block
      common /pparms/ nproc, nprec, lgrp, jbname
      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 ---------------------------------
csm this segment is used for shared memory computers
csm  40 do 60 k = 1, kblok
csm     do 50 j = 1, kxp
csm     ll = mixup(j,k)
csm     kk = (ll - 1)/kxp + 1
csm     jj = ll - kxp*(kk - 1)
csm     g(j,k) = f(jj,kk)
csm  50 continue
csm  60 continue

cdm this segment is used for distributed memory computers
cdm nvp = number of processors involved
   40 nvp = nx/kxp
cdm kxpp = length of data per message
      kxpp = (kxp - 1)/nvp + 1
cdm kxb = number of messages sent
      kxb = kxp/kxpp
      nvpp = (nvp - 1)/kxp + 1
cdm pack buffers
      do 60 k = 1, kblok
      do 50 j = 1, kxp
      g(j,k) = f((mixup(j,k)-1)/nvp+1,k)
   50 continue
   60 continue
cdm send and receive buffers
      do 80 k = 1, kblok
      do 70 i = 1, kxb
      ir = ieor((mixup(i,k) - 1)/kxp,((k+ks)/nvpp)*nvpp) + 1
      is = ir
      kk = (mixup((ir-1)/nvpp+1,k) - 1)/kxp

csp1 this segment is used for distributed memory sp1 computers
      call mp_recv(f(1+(kxp*kk)/nvp,k),2*kxpp*nprec,ir-1,kxp-kxb+i,msid)
      call mp_bsend(g(1+(kxp*(is-1))/nvp,k),2*kxpp*nprec,is-1,kxp-kxb+i)
      call mp_wait(msid,nbytes)
   70 continue
   80 continue
c unpack buffers
      do 110 k = 1, kblok
      do 100 i = 1, kxb
      do 90 j = 1, kxpp
      g(i+kxb*(j-1),k) = f(j+kxpp*(i-1),k)
   90 continue
  100 continue
  110 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
C =========================================================================
c
c - - - this subroutine performs timing
c       icntrl = (-1,0,1) = (initialize,ignore,read) clock
c                clock should be initialized before it is read!
c       chr = character variable for labeling timings
c
C =========================================================================
      subroutine timera(icntrl,chr)
c written for ibm sp1
      character*8 chr, jbname
      integer jclock, nclock
      double precision clock, work1
      external d_vmax
      common /pparms/ nproc, nprec, lgrp, jbname
      save jclock
   91 format (1x,a8,1x,'time = ',e14.7,1x,'sec')
      data jclock /0/
      if (icntrl.eq.0) return
      if (icntrl.eq.1) go to 10
c initialize clock
      call mp_sync(lgrp)
      jclock = mclock(d)
      return
c read clock and write time difference from last clock initialization
   10 nclock = mclock(d)
      clock = float(nclock - jclock)/100.
      call mp_combine(clock,work1,8,d_vmax,lgrp)
      clock = work1
      call mp_environ(iscr,idproc)
      if (idproc.eq.0) write (6,91) chr, clock
      return
      end

c parallel sum of a vector for pvm
c     INCLUDE 'vector/psum.f'

c parallel sum of a vector for pvm
c  INCLUDE 'vector/pisum.f'
