
      SUBROUTINE INDEP_LOOP_PPUSH1(part,fx,npp,noff,
     *   qtme,dt,sum_0,idimp,npmax,nblok,nxpmx) 
	
    	REAL PART(idimp,npmax,nblok)
	INTEGER NOFF(nblok),NPP(nblok)
	REAL FX(nxpmx,nblok)
	REAL QTME,DT
	REAL SUM_0(nblok)
        
	integer lb,ub,nn,j,k,nnoff

	real dx
        include 'pglocal.f'
        k = 1
   	sum_0(k) = 0.
        nnoff = noff(k) - 2
        do 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)+ fx(nn-1,k)*((.5 - dx)**2))
C ----------------------------
C new velocity
C ----------------------------
   	   dx = part(2,j,k) + qtme*dx
C ----------------------------
C average kinetic energy
C ----------------------------
           sum_0(k) = sum_0(k) + (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
      	enddo
        return

        end 



        SUBROUTINE permute(F,g,mixup,indx,kxp,kblok)

        complex f,g
	dimension F(kxp,kblok),G(kxp,kblok)
        dimension mixup(kxp,kblok)  

	integer myprocid
	integer dest,source,nproc,dest_id
	complex global_f(kxp,kblok)
        include 'pglocal.f'
	
	nproc = kblok
	myprocid = pghpf_myprocnum()

C ---------------------------
C inverse fourier transform
C bit-reverse array elements to temporary
C ---------------------------
       
         do j=1,kxp
           global_f(j,myprocid+1) = f(j,1)
         enddo
       

         kcnt = myprocid
         dest   = mod(myprocid+1,nproc)
         source = mod(myprocid-1+nproc,nproc) 

         do i=1,nproc-1
             if (mod (myprocid,2) .eq. 0) then
                call pghpf_csend(dest,global_f(1,kcnt+1),kxp,1,
     *                         PGLCL_CPLX8)
             else
                ikcnt = mod(kcnt-1+nproc,nproc)
                call pghpf_crecv(source,global_f(1,ikcnt+1),kxp,1,
     *                         PGLCL_CPLX8)
             endif

             if (mod (myprocid,2) .eq. 1) then
                call pghpf_csend(dest,global_f(1,kcnt+1),kxp,1,
     *                         PGLCL_CPLX8)
             else
                ikcnt = mod(kcnt-1+nproc,nproc)
                call pghpf_crecv(source,global_f(1,ikcnt+1),kxp,1,
     *                         PGLCL_CPLX8)   
             endif
             kcnt = ikcnt
         enddo
         
         k=1
         do j = 1, kxp
            ll = mixup(j,k)
            kk = (ll - 1)/kxp + 1
            jj = ll - kxp*(kk - 1)
            g(j,k) = global_f(jj,kk)
         enddo

         return
         end





	SUBROUTINE FFT(F,g,isign,mixup,sct,indx,kxp,kblok)



        complex f,g,sct
	dimension F(kxp,kblok),G(kxp,kblok)
        dimension sct(kxp,kblok),mixup(kxp,kblok)  
     

  	INTEGER LB,UB,sign,info
	integer myprocid
	integer dest,source,nproc,dest_id,lind
	Real m1,m2 
	complex t1,global_f(kxp,kblok)
	parameter(TAG=1)
        include 'pglocal.f'
	
	nproc = kblok
	nx = 2**indx
	nxh = nx/2
	indx1 = indx-1
	kbs = nxh/kxp
	ks = -1
	dnx = 6.28318530717959/float(nx)
	sign = isign
	temp= log10(float(kxp))/log10(2.0)
        lind = temp
C important note: log10(2048.0)/log10(2.0) is not equal to 11
C with this xlf compiler if target variable is integer
C Therefore we assign the result to real var.    

	myprocid = pghpf_myprocnum()
        lb = 1
        ub = 1
C ---------------------------
C inverse fourier transform
C bit-reverse array elements to temporary
C ---------------------------
           
         do  l = 1, indx1
            nxs = 2**(l - 1)
            kxs = nxs/kxp
C ---------------------------
C local calculation
C ---------------------------

            if (kxs.eq.0) then
               k=1
               do 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
               enddo              
C ---------------------------
C perform reduction
C ---------------------------
               km = kxp/nxs
               do j = 1, kxp
                   kb = j - 1
                   lb = kb/nxs
                   if (sign .eq. -1) then
                      t1 = sct(1+km*(kb-nxs*lb),k)
                   else if (sign .eq.1) then
                      t1 = conjg(sct(1+km*(kb-nxs*lb),k) )
                   endif
                   
                   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
                enddo
            else
C ---------------------------
C copy data
C ---------------------------
               k = 1
               if (btest(myprocid,l-lind-1) .eqv. .false.) then
                 dest_id = ibset(myprocid,l-lind-1)
                 dest = dest_id
                 call pghpf_csend(dest,g(1,k),kxp,1,PGLCL_CPLX8)
               else
                 dest_id = ibclr(myprocid,l-lind-1)
                 source = dest_id
                 call pghpf_crecv(source,f(1,k),kxp,1,PGLCL_CPLX8)
               endif

               if (btest(myprocid,l-lind-1) .eqv.  .true.) then
                 dest_id = ibclr(myprocid,l-lind-1)
                 dest = dest_id
                 call pghpf_csend(dest,g(1,k),kxp,1,PGLCL_CPLX8)
               else
                 dest_id = ibset(myprocid,l-lind-1)                       
                 source = dest_id
                 call pghpf_crecv(source,f(1,k),kxp,1,PGLCL_CPLX8)
               endif

C ------------------
C perform reduction
C ------------------
               km = nxh/nxs
               dns = dnx*float(km)
               do k = lb,ub
                  kb = myprocid+1 + ks 
                  llb = kb/kxs
                  kb = kxp*(kb - kxs*llb) - 1
                  llb = llb - 2*(llb/2)
                  do j = 1, kxp
                     arg = dns*float(j + kb)
                     t1 = cmplx(cos(arg),sign*sin(arg))
                     if (llb.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
                  enddo
                enddo
            endif

        enddo
C do l=1,indx1
  
        return
        end 



        SUBROUTINE SCRAMBLE(F,G,isign,indx,kxp,kblok,SHIFTED)

        COMPLEX F(kxp,kblok),G(kxp,kblok)
        COMPLEX SHIFTED(kblok)

        INTEGER LB,UB,sign,mytid,info
        integer myinst,tids(kblok),bufid,myprocid
        integer dest,source,nproc
        Real m1,m2
        complex first
        parameter(TAG=1)

        include 'pglocal.f'

        LB = 1
        UB = 1

        nproc = kblok
        nx = 2**indx
        nxh = nx/2
        kbs = nxh/kxp
        ks = -1
        dnx = 6.28318530717959/float(nx)
        kxp2 = kxp +2
        sign = -isign

        myprocid = pghpf_myprocnum()

C calculate multipliers used in scarambling
        if (isign .eq. -1) then
           m1 = 0.5/float(nx)
           m2 = 2.0*m1
        else if (isign .eq. 1) then
           m1 = 1.0
           m2 = 1.0
        endif

        do k=lb,ub
           kb = myprocid+1+ks  
           kk = kbs-kb

           if (myprocid .lt. nproc/2) then
              source = kk-1
              call pghpf_crecv(source,f(2,k),kxp-1,1,
     *                         PGLCL_CPLX8)
              call pghpf_crecv(source,first,1,1,
     *                         PGLCL_CPLX8)
              if (myprocid .gt. 0) f(1,k) = first
           else
              dest = kk-1
              call pghpf_csend(dest,g(2,k),kxp-1,1,
     *                         PGLCL_CPLX8)
              call pghpf_csend(dest,shifted(k),1,1,
     *                         PGLCL_CPLX8)
           endif

           if (myprocid .ge. nproc/2) then
             source = kk-1
             call pghpf_crecv(source,f(2,k),kxp-1,1,
     *                         PGLCL_CPLX8)
             call pghpf_crecv(source,first,1,1,
     *                         PGLCL_CPLX8)
             f(1,k) = first
           else
              dest = kk-1
              call pghpf_csend(dest,g(2,k),kxp-1,1,
     *                         PGLCL_CPLX8)
              call pghpf_csend(dest,shifted(k),1,1,
     *                         PGLCL_CPLX8)
           endif
        enddo


C ------------------
C perform reduction
C ------------------
          do k = lb,ub
             kb = myprocid+1+ks   
             koff = kxp*kb - 1
             do j = 2, kxp
                jj = kxp2 - j
                arg = dnx*float(j + koff)
                g(j,k) = m1*((g(j,k) + conjg(f(jj,k))) 
     *              - (g(j,k) 
     *              - conjg(f(jj,k)))*cmplx(sin(arg),sign*cos(arg)))
             enddo
             if (kb .eq. 0) then
                g(1,k) = m2*cmplx(real(g(1,k)) + 
     *                     aimag(g(1,k)),real(g(1,k)) 
     *                     - aimag(g(1,k)))
             else
                arg = dnx*float(1 + koff)
                g(1,k) = m1*((g(1,k) + conjg(f(1,k))) - 
     *                     (g(1,k) - conjg(f(1,k)))* 
     *                     cmplx(sin(arg),sign*cos(arg)))
             endif
C copy result to destination
             do j = 1, kxp
                f(j,k) = g(j,k)
             enddo
        enddo
        return
        end 



	SUBROUTINE INDEP_LOOP_PDOST1(PART,Q,NOFF,NPP,QME,idimp,
     *     npmax,kblok,nblok,nxpmx)

   	REAL  PART(idimp,npmax,nblok)
    	INTEGER NOFF(nblok),NPP(nblok)
        REAL  Q(nxpmx,nblok)
        REAL QME 
        include 'pglocal.f'

        INTEGER LB,UB,p
   
        LB = 1
        UB = 1

        p = pghpf_myprocnum()

        k=1
        nnoff = noff(k) - 2
        do 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) + qme*(.75 - dx*dx)
            nn = nn + 1
            q(nn,k) = q(nn,k) + .5*qme*((.5 + dx)**2)
            nn = nn - 2
            q(nn,k) = q(nn,k) + .5*qme*((.5 - dx)**2)
       enddo

       RETURN
       END 

       subroutine do_forall(f,fc,kxp,kblok)
       real f(kxp,kblok)
       complex fc(kxp,kblok)
       
       k = 1
       do j=1,kxp
          f(2*j,k) = real(fc(j,k))
          f(2*j+1,k) = aimag(fc(j,k))
       enddo
       return
       end
       
C Each processor in our application will get one blok of data for any
C variable data. I mean that processor k will get all the k'th partiton
C of each data

	SUBROUTINE INDEP_LOOP_PIST1_1(NPP,NPS,PART,SUM_0,idimp,nblok,
     *   npmax)
        INTEGER NPP(nblok),NPS(nblok)
        REAL    SUM_0(nblok)
        REAL PART(idimp,npmax,nblok)

        INTEGER J
        integer lb,ub,pnum
        include 'pglocal.f'
        lb = 1
        ub = 1

        pnum = pghpf_myprocnum()

C each processor will get only one chunk of data, k=1 lbound = ubound
C for each partitioned data

        DO k=lb,ub

            SUM_0(k)= 0
            DO J=NPS(k),NPP(k)
               SUM_0(k) = SUM_0(k) + PART(2,J,k)
            ENDDO
        ENDDO
        RETURN
        END 


        SUBROUTINE INDEP_LOOP_PIST1_2(NPP,NPS,PART,SUM_1,idimp,
     *    nblok,npmax)

        INTEGER NPP(nblok),NPS(nblok)
        REAL PART(idimp,npmax,nblok)
        REAL    SUM_1

        INTEGER J
        INTEGER UB,LB
        include 'pglocal.f'
        lb = 1
        ub = 1
        DO K=LB, UB
           DO J=NPS(K),NPP(K)
              PART(2,J,K) = PART(2,J,K) - SUM_1
           ENDDO
        ENDDO
    
        RETURN
        END 


        SUBROUTINE INDEP_LOOP_PPOIS1_1(FFC,NX,KXP,AFFP,AX,kblok)
        REAL  AFFP,AX
        COMPLEX FFC(kxp,kblok)

        integer pnum
        INTEGER LB,UB
        integer lindex(2),g_index(2)
        data lindex/2*1/
        include 'pglocal.f'
        dnx = 6.28318530717959/float(nx)
C ----------------------------
C prepare form factor array
C ----------------------------
        lb = 1
        ub = 1
        pnum = pghpf_myprocnum() + 1
        DO K=LB,UB 
            koff = kxp*(pnum -1) - 1
            do  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
            enddo

       ENDDO
      
       return
       END 

      SUBROUTINE INDEP_LOOP_PPOIS1_2(FFC,QC,FC,NX,KXP,
     *     AFFP,AX,sum_0,kblok,nblok)

       REAL  AFFP,AX
       COMPLEX FFC(kxp,kblok),QC(kxp,kblok)
       COMPLEX FC(kxp,kblok)
       REAL sum_0(nblok)        
       include 'pglocal.f'
         integer pnum
         integer lb,ub
         integer lindex(2),g_index(2)
         data lindex/2*1/

C ---------------------------------------------
C calculate force/charge and sum field energy
C ---------------------------------------------
         dnx = 6.28318530717959/float(nx)    

         k = 1
         pnum = pghpf_myprocnum()+1
         koff = kxp*(pnum-1) - 1
         sum_0(k) = 0.
         do  j = 1, kxp
            at1 = real(ffc(j,k))*aimag(ffc(j,k))
            at2 = dnx*float(j + koff)*at1
            fc(j,k) = at2*cmplx(aimag(qc(j,k)),-real(qc(j,k)))
            sum_0(k) = sum_0(k) + at1*real(qc(j,k)*conjg(qc(j,k)))
         enddo


        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(mixup,isign,indx,kxp,kblok,sct)

        complex  sct(kxp,kblok)
        integer  mixup(kxp,kblok)
        integer isign,indx,kxp,kblok

        integer llb,ub,myprocnum
        real arg,dnx
        include 'pglocal.f'
       
        llb = 1
        ub = 1
        myprocnum = pghpf_myprocnum()

 
      nx = 2**indx
      nxh = nx/2
      dnx = 6.28318530717959/float(nx)
      ks = -1

C ---------------------------------
C prepare bit-reverse index table
C ---------------------------------
      do  k = llb,ub
         koff = kxp*myprocnum - 1
         do  j = 1, kxp
            lb = j + koff
            ll = 0
            do  l = 1, indx
               jb = lb/2
               it = lb - 2*jb
               lb = jb
               ll = 2*ll + it
            enddo
            mixup(j,k) = ll + 1
         enddo
      enddo


       
      do k = llb,ub
        do j = 1, kxp
           arg = dnx*float((nxh*(j - 1))/kxp)
           sct(j,k) = cmplx(cos(arg),-sin(arg))
        enddo
      enddo

      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 dist_particle (part,edges,npp,npx,vtx,
     *    nblok,idimp,npmax,nx)


       integer  npx,nblok,idimp,npmax,nx
       real   vtx
       real part(idimp,npmax,nblok)
       real edges(idimp,nblok)
       integer npp(nblok)

       double precision ranorm
       integer myinst,bufid,myprocid,mytid,tids(nblok)
       integer dest,source,nproc,lb,ub,offset,p,num_of_part
       include 'pglocal.f'
       parameter (TAG=1)
      
       nproc = nblok
       lb =  1
       ub = 1
        
       myprocid = pghpf_myprocnum()
       item_cnt = npx/nblok

C particle distribution
       if (myprocid .eq.nproc-1) then
          at1 = float(nx)/float(npx)
          do p=0,nproc-1
              num_of_part = 0
              offset = npp(1)
              do j = 1,item_cnt 
                 xt = at1*(float(j+p*item_cnt) - 0.5)
                 vxt = vtx*ranorm(d)
                 num_of_part = num_of_part +1
                 part(1,num_of_part+offset,1) = xt
                 part(2,num_of_part+offset,1) = vxt
              enddo
              if (p .lt. nproc-1 ) then
                call pghpf_csend(p,item_cnt,1,1,PGLCL_INT4)
                call pghpf_csend(p,part(1,offset+1,1),2*item_cnt,1,
     *              PGLCL_REAL4)
              else
                 npp(1) = npp(1) + item_cnt
              endif
           enddo
       else
            call pghpf_crecv(nproc-1,num_of_part,1,1,PGLCL_INT4)
            offset = npp(lb)+1
            call pghpf_crecv(nproc-1,part(1,offset,lb),2*num_of_part,1,
     *             PGLCL_REAL4)   
            npp(lb)= npp(lb) + num_of_part
       endif
     
       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 


      SUBROUTINE INDEP_LOOP_PMOVE1_1(PART,EDGES,NPP,  
     *   SBUFR,SBUFL,IHOLE,JSR,JSL,JSS,NX,NVP,NBMAX,idimp,npmax,
     *   nblok,ntmax)

      REAL PART(idimp,npmax,nblok)
      REAL SBUFR(idimp,nbmax,nblok),SBUFL(idimp,nbmax,nblok)
      REAL  EDGES(idimp,nblok)
      INTEGER JSR(idimp,nblok),JSL(idimp,nblok),JSS(idimp,nblok)
      INTEGER IHOLE(ntmax,nblok)
      INTEGER NX,NVP,NBMAX
      INTEGER NPP(nblok)
      include 'pglocal.f'

C --------------------------
C buffer outgoing particles
C --------------------------

      integer lb,ub,myprocid

      
      lb = 1
      ub = 1


      myprocid = pghpf_myprocnum()
      anx = float(nx)
      ks = -1

      k = 1

      kb = myprocid+1+ks 
      jsl(1,k) = 0
      jsr(1,k) = 0
      jss(2,k) = 0
      do  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
      enddo
      jss(1,k) = jsl(1,k) + jsr(1,k)

      return 

      END  

      SUBROUTINE INDEP_LOOP_PMOVE1_2(EDGES,rbufr,
     *     rbufl,JSR,JSL,ISUM_0,idimp,nbmax,nblok)

        REAL RBUFR(idimp,nbmax,nblok),RBUFL(idimp,nbmax,nblok)
        REAL  EDGES(idimp,nblok)
        INTEGER  JSR(idimp,nblok),JSL(idimp,nblok)
        integer isum_0(nblok)
        include 'pglocal.f'

       integer lb,ub

       lb = 1
       ub = 1 

C -------------------------------------------
C check if particles must be passed further
C -------------------------------------------


      DO k=lb,ub

      jsl(1,k) = 0
      do  j = 1, jsr(2,k)
         if (rbufr(1,j,k).lt.edges(1,k)) jsl(1,k) = 1
      enddo

      jsr(1,k) = 0
      do  j = 1, jsl(2,k)
         if (rbufl(1,j,k).ge.edges(2,k)) jsr(1,k) = 1
      enddo
      ISUM_0(k) = (jsl(1,k) + jsr(1,k))
      ENDDO

      return
      END 

	  SUBROUTINE INDEP_LOOP_PMOVE1_3(EDGES,SBUFR,SBUFL,
     *    RBUFR,RBUFL,JSR,JSL,JSS,NVP,NX,idimp,nbmax,nblok)


      REAL SBUFR(idimp,nbmax,nblok),SBUFL(idimp,nbmax,nblok)
      REAL RBUFR(idimp,nbmax,nblok),RBUFL(idimp,nbmax,nblok)
      REAL  EDGES(idimp,nblok)
      INTEGER  JSR(idimp,nblok),JSL(idimp,nblok),JSS(idimp,nblok)
      INTEGER  NX,NVP

      integer lb,ub,myprocid
      include 'pglocal.f'
      lb = 1
      ub = 1
      myprocid = pghpf_myprocnum()

C -------------------------------------------
C remove particles which do not belong here
C -------------------------------------------
      anx = float(nx)
      ks = -1

      DO k=lb,ub

      kb = myprocid+1+ks   
      jsl(1,k) = 0
      jss(2,k) = 0
      do 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
      enddo

      jsr(2,k) = jss(2,k)

      jsr(1,k) = 0
      jss(2,k) = 0
      do  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
      enddo
      jsl(2,k) = jss(2,k)
      ENDDO

      return
      END 

      SUBROUTINE INDEP_LOOP_PMOVE1_4(JSR,JSL,JSS,NPMAX, 
     *           ISUM_0,NPP,idimp,nblok)

      INTEGER  JSR(idimp,nblok),JSL(idimp,nblok),JSS(idimp,nblok)
      integer isum_0(nblok),npp(nblok)
      INTEGER NPMAX
      include 'pglocal.f'

      
      integer lb,ub

      lb = 1
      ub = 1

      DO k=lb,ub
    
      it = (npp(k)+jsl(2,k)+jsr(2,k)-jss(1,k))

      if (it.gt.npmax) then
          jss(2,k) = 1
      else
          jss(2,k) = 0
      endif

      ISUM_0(k) = jss(2,k)
      ENDDO

      return
      END 

      SUBROUTINE INDEP_LOOP_PMOVE1_5(PART,EDGES,NPP,
     *  RBUFR,RBUFL,IHOLE,JSR,JSL,JSS,idimp,npmax,nblok,
     *  nbmax,ntmax)

      REAL PART(idimp,npmax,nblok)
      REAL RBUFR(idimp,nbmax,nblok),RBUFL(idimp,nbmax,nblok)
      REAL  EDGES(idimp,nblok)
      INTEGER  JSR(idimp,nblok),JSL(idimp,nblok),JSS(idimp,nblok)
      integer ihole(ntmax,nblok)
      integer npp(nblok)
      include 'pglocal.f'

      integer lb,ub

      lb = 1
      ub = 1

C ------------------------------------------------------
C distribute incoming particles from buffers
C ------------------------------------------------------
      DO k=lb,ub

      jss(2,k) = min0(jss(1,k),jsl(2,k))
      do  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)
      enddo

      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 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
      enddo

      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 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
      enddo

      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

      ENDDO

      return
      END  

