      Program beps2k
      include 'mpif.h'
      double precision dseed 
      parameter(indx=7,indy=8,npx=10*128,npy=10*256)
      parameter(npxb=3*128,npyb=3*256)
      parameter(tend=65.0000,dt=0.20)
      parameter(vty=1.0,vdy=5.0,vtdy=0.5,avdy=5.0)
      parameter(indnvp=3,idps=2,idimp=4,nvp=2**indnvp)
      parameter(npxy=npx*npy,npxyb=npxb*npyb,np=npxy+npxyb)
      parameter(nx=2**indx,ny=2**indy,nxh=nx/2,nyh=ny/2)
      parameter(nxv=nx+2,nyv=ny+2,nxvh=nxv/2)
      parameter(nloop=tend/dt+.0001)
      parameter(npmax=(np/nvp)*1.20+7000,nypmx=(ny-1)/nvp+4)
      parameter(kyp=(ny-1)/nvp+1,kxp=(nxh-1)/nvp+1)
      parameter(nmx=nx*(ny/nx)+ny*(nx/ny),nxy=nmx/(2-nx/nmx-ny/nmx))
      parameter(nxyh=nxy/2,nmxh=nxh*(ny/nxh)+ny*(nxh/ny))
      parameter(nxhy=nmxh/(2-nxh/nmxh-ny/nmxh))
      parameter(nbmax=1+(2*(npxy*vty+npxyb*vtdy)+2.*npxyb*avdy)*dt/ny)
      complex qc,fxc,fyc,qt,fxt,fyt,ffc,bs,br,sct,sc
      dimension tarray(2)
      dimension x(npmax),y(npmax),u(npmax),v(npmax)
      dimension q(nxv,nypmx),fx(nxv,nypmx),fy(nxv,nypmx)
      dimension qc(nxvh,kyp),fxc(nxvh,kyp),fyc(nxvh,kyp)
      dimension qt(nyv,kxp),fxt(nyv,kxp),fyt(nyv,kxp)
      dimension ffc(nyh,kxp),mixup(nxhy),sct(nxyh)
      dimension bs(kxp,kyp),br(kxp,kyp),edges(idps)
      dimension sbufl(idimp,nbmax),sbufr(idimp,nbmax)
      dimension rbufl(idimp,nbmax),rbufr(idimp,nbmax)
      dimension sc(nxvh,idps),sb(nxv,idps)
      data vtx,qme,vdx,vtdx /1.0,-1.0,0.0,.5/
      data ax,ay,dseed/.866667,.866667,123457.D0/
      data ntpose /1/

      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ierr)

      iout=myid+8
      kstrt = myid + 1

      itime = 0
      qtme = qme*dt
      affp = float(nx*ny)/float(np)
      zero = 0.0
      call dcomp2(edges,nyp,noff,ny,kstrt,nvp,idps)

      isign = 0
      call pfft2r(qc,qt,bs,br,isign,ntpose,mixup,sct,indx,indy,kstrt,
     &            nxvh,nyv,kxp,kyp,nxhy,nxyh)
      call ppois2(qt,fxt,fyt,isign,ffc,ax,ay,affp,we,nx,ny,kstrt,nyv,
     &            kxp,nyh)

      nyp3 = nyp + 3
      do 120 k = 1, nyp3
      do 110 j = 1, nx
      q(j,k) = 0.0
  110 continue
  120 continue

      npp=0
      if (npxy.gt.0) call pistr2(x,y,u,v,edges,npp,vtx,vty,zero,zero,
     &     npx,npy,nx,ny,idimp,npmax,idps,dseed)

      if (npxyb.gt.0) call pistr2(x,y,u,v,edges,npp,vtdx,vtdy,vdx,vdy,
     &     npxb,npyb,nx,ny,idimp,npmax,idps,dseed)
      call opdost2(x,y,q,npp,noff,qme,nx,idimp,npmax,nxv,nypmx)
 
      qi0 = -qme/affp
      do 230 k = 1, nyp
      do 220 j = 1, nx
      q(j,k+1)=q(j,k+1)+qi0
  220 continue
  230 continue     
      tpu=0.
      tmi=0.
      tde=0.
      toh=0.

c
c * * * start main iteration loop * * *
c

  500 if (itime.ge.nloop) go to 2000
      isign = -1
      tt1=dtime(tarray)
      call cppfp2(q,qc,isign,sb,sc,kstrt,nvp,nx,nxv,nypmx,nxvh,
     &            kyp,idps)
      call pfft2r(qc,qt,bs,br,isign,ntpose,mixup,sct,indx,indy,kstrt,
     &            nxvh,nyv,kxp,kyp,nxhy,nxyh)
      call ppois2(qt,fxt,fyt,isign,ffc,ax,ay,affp,we,nx,ny,kstrt,nyv,
     &            kxp,nyh)

      isign=1
      call pfft2r(fxc,fxt,bs,br,isign,ntpose,mixup,sct,indx,indy,kstrt,
     &            nxvh,nyv,kxp,kyp,nxhy,nxyh)
      call pfft2r(fyc,fyt,bs,br,isign,ntpose,mixup,sct,indx,indy,kstrt,
     &            nxvh,nyv,kxp,kyp,nxhy,nxyh)
      call cppfp2(fx,fxc,isign,sb,sc,kstrt,nvp,nx,nxv,nypmx,nxvh,
     &            kyp,idps)
      call cppfp2(fy,fyc,isign,sb,sc,kstrt,nvp,nx,nxv,nypmx,nxvh,
     &            kyp,idps)
      nyp3 = nyp + 3
      do 1240 k = 1, nyp3
      do 1230 j = 1, nx
      q(j,k) = 0.0
 1230 continue
 1240 continue
      tt2=dtime(tarray)
      toh=toh+tt2
      wke = 0.0
      tt1=dtime(tarray)
      call  oppush2(x,y,u,v,fx,fy,npp,noff,qtme,dt,wke,nx,ny,idimp,
     &  npmax,nxv,nypmx,sbufl,sbufr,nbmax,edges,il,ir,ii,idps,
     &  kstrt,nvp)
      tt2=dtime(tarray)
      tpu=tpu+tt2 
      tt1=dtime(tarray)
      call pmove2(x,y,u,v,npp,idimp,npmax,sbufl,sbufr,
     &  rbufl,rbufr,nbmax,kstrt,il,ir,ii,nvp)
      tt2=dtime(tarray)
      tmi=tmi+tt2
      call opdost2(x,y,q,npp,noff,qme,nx,idimp,npmax,nxv,nypmx)
      tt2=dtime(tarray)
      tde=tde+tt2
      tt1=dtime(tarray)
      do 1270 k = 1, nyp
      do 1260 j = 1, nx
      q(j,k+1) = q(j,k+1) + qi0
 1260 continue
 1270 continue
      tt2=dtime(tarray)
      toh=toh+tt2

      wt = we + wke
      itime = itime + 1
      write(iout,*) 'itime=',itime,we,wke,wt

      go to 500
 2000 continue

c
c * * * end main iteration loop * * *
c
      tto=tde+tpu+toh+tmi
      write(iout,*) 'The total time =',tto
      write(iout,*) 'The push time =',tpu
      write(iout,*) 'The migrate time =',tmi
      write(iout,*) 'The despost time=',tde
      write(iout,*) 'The other time is=',toh
      call MPI_FINALIZE(ierr)
      stop
      end



      subroutine dcomp2(edges,nyp,noff,ny,kstrt,nvp,idps)
      dimension edges(idps)
      at1 = float(ny)/float(nvp)
      kb = kstrt - 1
      edges(1) = at1*float(kb)
      noff = edges(1) + .5
      edges(2) = at1*float(kb + 1)
      kr = edges(2) + .5
      nyp = kr - noff
      return
      end



      subroutine pistr2(x,y,u,v,edges,npp,vtx,vty,vdx,vdy,npx,npy,nx,
     &   ny,idimp,npmax,idps,dseed)
      double precision dseed
      dimension x(npmax),y(npmax),u(npmax),v(npmax)
      dimension edges(idps)
      at1 = float(nx)/float(npx)
      at2 = float(ny)/float(npy)
      do 30 k = 1, npy
      yt = at2*(float(k) - .5)
      do 20 j = 1, npx
      xt = at1*(float(j) - .5)
      vxt = vdx+vtx*(randdb(dseed)-0.5)
      vyt = vdy+vty*(randdb(dseed)-0.5)
      if ((yt.ge.edges(1)).and.(yt.lt.edges(2))) then
      npt = npp + 1
      x(npt) = xt
      y(npt) = yt
      u(npt) = vxt
      v(npt) = vyt
      npp = npt
      endif
   20 continue
   30 continue
      return
      end



      subroutine opdost2(x,y,q,npp,noff,qm,nx,idimp,npmax,nxv,nypmx)
      dimension x(npmax),y(npmax),q(nxv,nypmx)
      mnoff = noff - 2
      do 10 j = 1, npp
      nn = x(j) + .5
      dxl = x(j) - float(nn)
      nn = nn + 1
      if (nn.gt.nx) nn = nn - nx
      amx = qm*(.75 - dxl*dxl)
      np = nn + 1
      if (np.gt.nx) np = np - nx
      dxp = .5*qm*(.5 + dxl)**2
      nl = nn - 1
      if (nl.lt.1) nl = nl + nx
      dxl = .5*qm*(.5 - dxl)**2
      mm = y(j) + .5
      dyl = y(j) - float(mm)
      mm = mm - mnoff
      amy = .75 - dyl*dyl
      mp = mm + 1
      dyp = .5*(.5 + dyl)**2
      ml = mm - 1
      dyl = .5*(.5 - dyl)**2
      q(nn,mm) = q(nn,mm) + amx*amy
      q(np,mm) = q(np,mm) + dxp*amy
      q(nl,mm) = q(nl,mm) + dxl*amy
      q(nn,mp) = q(nn,mp) + amx*dyp
      q(np,mp) = q(np,mp) + dxp*dyp
      q(nl,mp) = q(nl,mp) + dxl*dyp
      q(nn,ml) = q(nn,ml) + amx*dyl
      q(np,ml) = q(np,ml) + dxp*dyl
      q(nl,ml) = q(nl,ml) + dxl*dyl
   10 continue
      return
      end




      subroutine cppfp2(f,fc,isign,sb,sc,kstrt,nvp,nx,nxv,nypmx,nxvh,
     &     kyp,idps)
      include 'mpif.h'
      complex fc,sc
      dimension f(nxv,nypmx),fc(nxvh,kyp)
      dimension sb(nxv,idps),sc(nxvh,idps)
      dimension istatus(MPI_STATUS_SIZE)
      nxh = nx/2
      moff = nypmx*nvp
      if (isign.eq.-1) then
      do 30 k = 1, kyp
      do 20 j = 1, nxh
      fc(j,k) = cmplx(f(2*j-1,k+1),f(2*j,k+1))
   20 continue
   30 continue
      kr = kstrt + 1
      if (kr.gt.nvp) kr = kr - nvp
      kl = kstrt - 1
      if (kl.lt.1) kl = kl + nvp
      call MPI_SEND(f(1,kyp+2),2*nxv,MPI_DOUBLE_PRECISION,kr-1,moff+1,
     &              MPI_COMM_WORLD,ierr)
      call MPI_RECV(sb,2*nxv,MPI_DOUBLE_PRECISION,kl-1,moff+1,
     &              MPI_COMM_WORLD,istatus,ierr)
      do 60 j = 1, nxh
      fc(j,1) = fc(j,1) + cmplx(sb(2*j-1,1),sb(2*j,1))
      fc(j,2) = fc(j,2) + cmplx(sb(2*j-1,2),sb(2*j,2))
   60 continue
      call MPI_SEND(f,nxv,MPI_DOUBLE_PRECISION,kl-1,moff+2,
     &              MPI_COMM_WORLD,ierr)
      call MPI_RECV(sb,nxv,MPI_DOUBLE_PRECISION,kr-1,moff+2,
     &              MPI_COMM_WORLD,istatus,ierr)
      do 80 j = 1, nxh
      fc(j,kyp) = fc(j,kyp) + cmplx(sb(2*j-1,1),sb(2*j,1))
   80 continue
      elseif(isign.eq.1) then
      do 130 k = 1, kyp
      do 120 j = 1, nxh
      f(2*j-1,k+1) = real(fc(j,k))
      f(2*j,k+1) = aimag(fc(j,k))
  120 continue
  130 continue
      kr = kstrt + 1
      if (kr.gt.nvp) kr = kr - nvp
      kl = kstrt - 1
      if (kl.lt.1) kl = kl + nvp
      call MPI_SEND(fc(1,kyp),nxvh,MPI_DOUBLE_COMPLEX,kr-1,moff+3,
     &               MPI_COMM_WORLD,ierr)
      call MPI_RECV(sc,nxvh,MPI_DOUBLE_COMPLEX,kl-1,moff+3,
     &               MPI_COMM_WORLD,istatus,ierr)
      do 233 i=1,nxvh
      f(2*i-1,1)=real(sc(i,1))
      f(2*i,1)=aimag(sc(i,1))
 233  continue
      call MPI_SEND(fc,2*nxvh,MPI_DOUBLE_COMPLEX,kl-1,moff+4,
     &               MPI_COMM_WORLD,ierr)
      call MPI_RECV(sc,2*nxvh,MPI_DOUBLE_COMPLEX,kr-1,moff+4,
     &               MPI_COMM_WORLD,istatus,ierr)
      do 333 i=1,nxvh
      f(2*i-1,kyp+2)=real(sc(i,1))
      f(2*i,kyp+2)=aimag(sc(i,1))
      f(2*i-1,kyp+3)=real(sc(i,2))
      f(2*i,kyp+3)=aimag(sc(i,2))
 333  continue
      endif
      return
      end



      subroutine ppois2(q,fx,fy,isign,ffc,ax,ay,affp,wp,nx,ny,kstrt,
     &      nyv,kxp,nyhd)
      include 'mpif.h'
      complex q,fx,fy,ffc,zero
      dimension q(nyv,kxp),fx(nyv,kxp),fy(nyv,kxp)
      dimension ffc(nyhd,kxp)
      nxh = nx/2
      nyh = ny/2
      ny2 = ny + 2
      dnx = 6.28318530717959/float(nx)
      dny = 6.28318530717959/float(ny)
      zero = cmplx(0.0,0.0)
      if (isign.eq.0) then
      if (kstrt.gt.nxh) return
      joff = kxp*(kstrt - 1) - 1
      do 20 j = 1, kxp
      dkx = dnx*float(j + joff)
      at1 = dkx*dkx
      at2 = (dkx*ax)**2
      do 10 k = 1, nyh
      dky = dny*float(k - 1)
      at3 = dky*dky + at1
      at4 = exp(-.5*((dky*ay)**2 + at2))
      if (at3.eq.0.0) then
      ffc(k,j) = cmplx(0.0,1.0)
      else
      ffc(k,j) = cmplx(affp*at4/at3,at4)
      endif
   10 continue
   20 continue
      else
      wp = 0.0
      if (kstrt.gt.nxh) go to 90
      joff = kxp*(kstrt - 1) - 1
      do 60 j = 1, kxp
      dkx = dnx*float(j + joff)
      if ((j+joff).gt.0) then
      do 50 k = 2, nyh
      k1 = ny2 - k
      at1 = real(ffc(k,j))*aimag(ffc(k,j))
      at2 = dkx*at1
      at3 = dny*float(k - 1)*at1
      fx(k,j) = at2*cmplx(aimag(q(k,j)),-real(q(k,j)))
      fx(k1,j) = at2*cmplx(aimag(q(k1,j)),-real(q(k1,j)))
      fy(k,j) = at3*cmplx(aimag(q(k,j)),-real(q(k,j)))
      fy(k1,j) = at3*cmplx(-aimag(q(k1,j)),real(q(k1,j)))
      wp = wp + at1*(q(k,j)*conjg(q(k,j)) + q(k1,j)*conjg(q(k1,j)))
   50 continue
      endif
      k1 = nyh + 1
      at1 = real(ffc(1,j))*aimag(ffc(1,j))
      fx(1,j) = dkx*at1*cmplx(aimag(q(1,j)),-real(q(1,j)))
      fx(k1,j) = zero
      fy(1,j) = zero
      fy(k1,j) = zero
      wp = wp + at1*(q(1,j)*conjg(q(1,j)))
   60 continue
      if (kstrt.eq.1) then
      do 70 k = 2, nyh
      k1 = ny2 - k
      at1 = real(ffc(k,1))*aimag(ffc(k,1))
      fx(k,1) = zero
      fx(k1,1) = zero
      fy(k,1)=dny*float(k-1)*at1*cmplx(aimag(q(k,1)),-real(q(k,1)))
      fy(k1,1) = zero
      wp = wp + at1*(q(k,1)*conjg(q(k,1)))
   70 continue
      endif
   90 continue
      temp = float(nx*ny)*wp
      call MPI_ALLREDUCE(temp,wp,1,MPI_DOUBLE_PRECISION,MPI_SUM,
     &        MPI_COMM_WORLD,ierr)
      endif
      return
      end



      subroutine oppush2(x,y,u,v,fx,fy,npp,noff,qtm,dt,ek,nx,ny,idimp,
     &  npmax,nxv,nypmx,sbufl,sbufr,nbmax,edges,il,ir,ii,idps,
     &  kstrt,nvp)
      include 'mpif.h'
      dimension x(npmax),y(npmax),u(npmax),v(npmax),edges(idps)
      dimension fx(nxv,nypmx),fy(nxv,nypmx)
      dimension sbufl(idimp,nbmax),sbufr(idimp,nbmax)
      zero = 0.0
      anx = float(nx)
      any = float(ny)
      sum1 = 0.0
      mnoff = noff - 2
      il=0
      ir=0
      ii=0
      do 10 j = 1, npp
      nn = x(j) + .5
      dx = x(j) - float(nn)
      nn = nn + 1
      if (nn.gt.nx) nn = nn - nx
      amx = .75 - dx*dx
      np = nn + 1
      if (np.gt.nx) np = np - nx
      dxp = .5*(.5 + dx)**2
      nl = nn - 1
      if (nl.lt.1) nl = nl + nx
      dxl = .5*(.5 - dx)**2
      mm = y(j) + .5
      dy = y(j) - float(mm)
      mm = mm - mnoff
      amy = .75 - dy*dy
      mp = mm + 1
      dyp = .5*(.5 + dy)**2
      ml = mm - 1
      dyl = .5*(.5 - dy)**2
      dx = amy*(amx*fx(nn,mm) + dxp*fx(np,mm) + dxl*fx(nl,mm)) + 
     &     dyp*(amx*fx(nn,mp) + dxp*fx(np,mp) + dxl*fx(nl,mp)) + 
     &     dyl*(amx*fx(nn,ml) + dxp*fx(np,ml) + dxl*fx(nl,ml))
      dy = amy*(amx*fy(nn,mm) + dxp*fy(np,mm) + dxl*fy(nl,mm)) + 
     &     dyp*(amx*fy(nn,mp) + dxp*fy(np,mp) + dxl*fy(nl,mp)) + 
     &     dyl*(amx*fy(nn,ml) + dxp*fy(np,ml) + dxl*fy(nl,ml))      
      dx = u(j) + qtm*dx
      dy = v(j) + qtm*dy
      sum1 = sum1 + (dx + u(j))**2 + (dy + v(j))**2
      u(j) = dx
      v(j) = dy
      dx = x(j) + dx*dt
      dy = y(j) + dy*dt
      if (dx.lt.zero) dx = dx + anx
      if (dx.ge.anx) dx = dx - anx
      if(dy.le.edges(1)) then
         il=il+1
         if(kstrt.eq.1) dy=dy+any
         sbufl(1,il)=dx
         sbufl(2,il)=dy
         sbufl(3,il)=u(j)
         sbufl(4,il)=v(j)
      elseif(dy.gt.edges(2)) then
         ir=ir+1
         if(kstrt.eq.nvp) dy=dy-any
         sbufr(1,ir)=dx
         sbufr(2,ir)=dy
         sbufr(3,ir)=u(j)
         sbufr(4,ir)=v(j)
      else
         ii=ii+1
         x(ii)=dx
         y(ii)=dy
         u(ii)=u(j)
         v(ii)=v(j)
      endif
   10 continue
      sum1 = .125*sum1
      call MPI_ALLREDUCE(sum1,ek,1,MPI_DOUBLE_PRECISION,MPI_SUM,
     &      MPI_COMM_WORLD,ierr)
      return
      end



      subroutine pmove2(x,y,u,v,npp,idimp,npmax,sbufl,sbufr,
     &  rbufl,rbufr,nbmax,kstrt,il,ir,ii,nvp)
      include 'mpif.h'
      dimension x(npmax),y(npmax),u(npmax),v(npmax)
      dimension sbufl(idimp,nbmax),sbufr(idimp,nbmax)
      dimension rbufl(idimp,nbmax),rbufr(idimp,nbmax)
      dimension istatus(MPI_STATUS_SIZE)
      kr = kstrt + 1
      if (kr.gt.nvp) kr = kr - nvp
      kl = kstrt - 1
      if (kl.lt.1) kl = kl + nvp
      call MPI_SEND(il,1,MPI_INTEGER,kl-1,nvp+3,
     &              MPI_COMM_WORLD,ierr)
      call MPI_SEND(sbufl,idimp*il,MPI_DOUBLE_PRECISION,kl-1,nvp+3,
     &              MPI_COMM_WORLD,ierr)
      call MPI_RECV(ia,1,MPI_INTEGER,kr-1,nvp+3,
     &              MPI_COMM_WORLD,istatus,ierr)
      call MPI_RECV(rbufr,idimp*ia,MPI_DOUBLE_PRECISION,kr-1,nvp+3,
     &              MPI_COMM_WORLD,istatus,ierr)
      npp=ii+ia
      do 100 i=ii+1,npp
      x(i)=rbufr(1,i-ii)
      y(i)=rbufr(2,i-ii)
      u(i)=rbufr(3,i-ii)
      v(i)=rbufr(4,i-ii)
 100  continue
      ii=npp
      call MPI_SEND(ir,1,MPI_INTEGER,kr-1,nvp+4,
     &              MPI_COMM_WORLD,ierr)
      call MPI_SEND(sbufr,idimp*ir,MPI_DOUBLE_PRECISION,kr-1,nvp+4,
     &              MPI_COMM_WORLD,ierr)
      call MPI_RECV(ib,1,MPI_INTEGER,kl-1,nvp+4,
     &              MPI_COMM_WORLD,istatus,ierr)
      call MPI_RECV(rbufl,idimp*ib,MPI_DOUBLE_PRECISION,kl-1,nvp+4,
     &              MPI_COMM_WORLD,istatus,ierr)
      npp=ii+ib
      do 200 i=ii+1,npp
      x(i)=rbufl(1,i-ii)
      y(i)=rbufl(2,i-ii)
      u(i)=rbufl(3,i-ii)
      v(i)=rbufl(4,i-ii)
 200  continue
      end




      subroutine pfft2r(f,g,bs,br,isign,ntpose,mixup,sct,indx,indy,
     &      kstrt,nxvh,nyv,kxp,kyp,nxhyd,nxyhd)
      complex f,g,bs,br,sct,s,t,t1
      dimension f(nxvh,kyp), g(nyv,kxp)
      dimension bs(kxp,kyp), br(kxp,kyp)
      dimension mixup(nxhyd), sct(nxyhd)
      indx1 = indx - 1
      indx1y = max0(indx1,indy)
      nx = 2**indx
      nxh = nx/2
      nxhh = nx/4
      nxh2 = nxh + 2
      ny = 2**indy
      nyh = ny/2
      ny2 = ny + 2
      nxy = max0(nx,ny)
      nxhy = 2**indx1y
      if (isign.eq.0) then
      do 30 j = 1, nxhy
      lb = j - 1
      ll = 0
      do 20 k = 1, indx1y
      jb = lb/2
      it = lb - 2*jb
      lb = jb
      ll = 2*ll + it
   20 continue
      mixup(j) = ll + 1
   30 continue
      nxyh = nxy/2
      dnxy = 6.28318530717959/float(nxy)
      do 40 j = 1, nxyh
      arg = dnxy*float(j - 1)
      sct(j) = cmplx(cos(arg),-sin(arg))
   40 continue
      elseif(isign.eq.-1) then
      if (kstrt.gt.ny) go to 180
      nrx = nxhy/nxh
      do 70 j = 1, nxh
      j1 = (mixup(j) - 1)/nrx + 1
      if (j.ge.j1) go to 70
      do 60 k = 1, kyp
      t = f(j1,k)
      f(j1,k) = f(j,k)
      f(j,k) = t
   60 continue
   70 continue
      nrx = nxy/nxh
      do 130 m = 1, indx1
      ns = 2**(m - 1)
      ns2 = ns + ns
      km = nxhh/ns
      kmr = km*nrx
      do 110 k = 1, km
      k1 = ns2*(k - 1)
      k2 = k1 + ns
      do 100 j = 1, ns
      j1 = j + k1
      j2 = j + k2
      s = sct(1+kmr*(j-1))
      do 90 i = 1, kyp
      t = s*f(j2,i)
      f(j2,i) = f(j1,i) - t
      f(j1,i) = f(j1,i) + t
   90 continue
  100 continue
  110 continue
  130 continue 
      kmr = nxy/nx
      ani = 1.0/float(2*nx*ny)
      do 150 j = 2, nxhh
      t1 = cmplx(aimag(sct(1+kmr*(j-1))),-real(sct(1+kmr*(j-1))))
      do 140 k = 1, kyp
      t = conjg(f(nxh2-j,k))
      s = f(j,k) + t
      t = (f(j,k) - t)*t1
      f(j,k) = ani*(s + t)
      f(nxh2-j,k) = ani*conjg(s - t)
  140 continue
  150 continue
      do 160 k = 1, kyp
      f(nxhh+1,k) = 2.0*ani*conjg(f(nxhh+1,k))
      f(1,k)=2.0*ani*cmplx(real(f(1,k))+aimag(f(1,k)),real(f(1,k))- 
     1aimag(f(1,k)))
  160 continue
  180 call ptpose(f,g,bs,br,nxh,ny,kstrt,nxvh,nyv,kxp,kyp)
      if (kstrt.gt.nxh) go to 290
      nry = nxhy/ny
      do 200 k = 1, ny
      k1 = (mixup(k) - 1)/nry + 1
      if (k.ge.k1) go to 200
      do 190 j = 1, kxp
      t = g(k1,j)
      g(k1,j) = g(k,j)
      g(k,j) = t
 190  continue
 200  continue
      nry = nxy/ny
      do 260 m = 1, indy
      ns = 2**(m - 1)
      ns2 = ns + ns
      km = nyh/ns
      kmr = km*nry
      do 240 k = 1, km
      k1 = ns2*(k - 1)
      k2 = k1 + ns
      do 230 j = 1, ns
      j1 = j + k1
      j2 = j + k2
      s = sct(1+kmr*(j-1))
      do 220 i = 1, kxp
      t = s*g(j2,i)
      g(j2,i) = g(j1,i) - t
      g(j1,i) = g(j1,i) + t
 220  continue
 230  continue
 240  continue
 260  continue
      if (kstrt.gt.1) go to 280
      do 270 k = 2, nyh
      s = g(ny2-k,1)
      g(ny2-k,1) = .5*cmplx(aimag(g(k,1)+s),real(g(k,1)-s))
      g(k,1) = .5*cmplx(real(g(k,1) + s),aimag(g(k,1)-s))
  270 continue
  280 continue
  290 if (ntpose.eq.0) call ptpose(g,f,br,bs,ny,nxh,kstrt,nyv,nxvh,
     &     kyp,kxp)
      elseif(isign.eq.1) then
      if(ntpose.eq.0) call ptpose(f,g,bs,br,nxh,ny,kstrt,nxvh,nyv,
     &     kxp,kyp)
      if (kstrt.gt.nxh) go to 410
      nry = nxhy/ny
      if (kstrt.gt.1) go to 320
      do 310 k = 2, nyh
      s = cmplx(aimag(g(ny2-k,1)),real(g(ny2-k,1)))
      g(ny2-k,1) = conjg(g(k,1) - s)
      g(k,1) = g(k,1) + s
  310 continue
  320 do 340 k = 1, ny
      k1 = (mixup(k) - 1)/nry + 1
      if (k.ge.k1) go to 340
      do 330 j = 1, kxp
      t = g(k1,j)
      g(k1,j) = g(k,j)
      g(k,j) = t
  330 continue
  340 continue
      nry = nxy/ny
      do 400 m = 1, indy
      ns = 2**(m - 1)
      ns2 = ns + ns
      km = nyh/ns
      kmr = km*nry
      do 380 k = 1, km
      k1 = ns2*(k - 1)
      k2 = k1 + ns
      do 370 j = 1, ns
      j1 = j + k1
      j2 = j + k2
      s = conjg(sct(1+kmr*(j-1)))
      do 360 i = 1, kxp
      t = s*g(j2,i)
      g(j2,i) = g(j1,i) - t
      g(j1,i) = g(j1,i) + t
  360 continue
  370 continue
  380 continue
  400 continue
  410 call ptpose(g,f,br,bs,ny,nxh,kstrt,nyv,nxvh,kyp,kxp)
      if (kstrt.gt.ny) return
      nrx = nxhy/nxh
      kmr = nxy/nx
      do 430 j = 2, nxhh
      t1 = cmplx(aimag(sct(1+kmr*(j-1))),real(sct(1+kmr*(j-1))))
      do 420 k = 1, kyp
      t = conjg(f(nxh2-j,k))
      s = f(j,k) + t
      t = (f(j,k) - t)*t1
      f(j,k) = s + t
      f(nxh2-j,k) = conjg(s - t)
  420 continue
  430 continue
      do 440 k = 1, kyp
      f(nxhh+1,k) = 2.*conjg(f(nxhh+1,k))
      f(1,k)=cmplx(real(f(1,k))+aimag(f(1,k)),real(f(1,k))-
     &       aimag(f(1,k)))
  440 continue
      do 460 j = 1, nxh
      j1 = (mixup(j) - 1)/nrx + 1
      if (j.ge.j1) go to 460
      do 450 k = 1, kyp
      t = f(j1,k)
      f(j1,k) = f(j,k)
      f(j,k) = t
  450 continue
  460 continue
      nrx = nxy/nxh
      do 520 m = 1, indx1
      ns = 2**(m - 1)
      ns2 = ns + ns
      km = nxhh/ns
      kmr = km*nrx
      do 500 k = 1, km
      k1 = ns2*(k - 1)
      k2 = k1 + ns
      do 490 j = 1, ns
      j1 = j + k1
      j2 = j + k2
      s = conjg(sct(1+kmr*(j-1)))
      do 480 i = 1, kyp
      t = s*f(j2,i)
      f(j2,i) = f(j1,i) - t
      f(j1,i) = f(j1,i) + t
  480 continue
  490 continue
  500 continue
  520 continue
      endif
      return
      end




      subroutine ptpose(f,g,s,t,nx,ny,kstrt,nxv,nyv,kxp,kyp)
      include 'mpif.h'
      complex f,g,s,t
      dimension f(nxv,kyp),g(nyv,kxp),s(kxp,kyp),t(kxp,kyp)
      dimension istatus(MPI_STATUS_SIZE)
      kxb = nx/kxp
      kyb = ny/kyp
      kxym = min0(kxb,kyb)
      mtr = kyb/kxym
      ntr = kxb/kxym
      mntr = max0(mtr,ntr)
      do 60 i = 1, kxym
      ir0 = iand(kxym-1,ieor(kstrt-1,i-1)) + 1
      is0 = ir0
      do 50 ii = 1, mntr
      if ((kstrt.le.nx).and.(ii.le.mtr)) then
      ir = ir0 + kxym*(ii - 1)
      endif
      if ((kstrt.le.ny).and.(ii.le.ntr)) then
      is = is0 + kxym*(ii - 1)
      joff = kxp*(is - 1)
      do 20 k = 1, kyp
      do 10 j = 1, kxp
      s(j,k) = f(j+joff,k)
   10 continue
   20 continue
      call MPI_SEND(s,kxp*kyp,MPI_DOUBLE_COMPLEX,is-1,kstrt+kxym+1,
     &              MPI_COMM_WORLD,ierr)
      endif
      if ((kstrt.le.nx).and.(ii.le.mtr)) then
      koff = kyp*(ir - 1)
      call MPI_RECV(t,kxp*kyp,MPI_DOUBLE_COMPLEX,ir-1,ir+kxym+1,
     &              MPI_COMM_WORLD,istatus,ierr)
      do 40 k = 1, kyp
      do 30 j = 1, kxp
      g(k+koff,j) = t(j,k)
   30 continue
   40 continue
      endif
   50 continue
   60 continue
      return
      end



      
      Function randdb(dseed)       
      double precision S2P31M,S2P31,DSEED,SEED
      DATA   S2P31M/2147483647.D0/,S2P31/2147483648.D0/
      SEED = DSEED
      SEED = DMOD(16807.D0*SEED,S2P31M)
      randdb= SEED /S2P31	
      DSEED = SEED
      RETURN
      END





