subroutine nlpsc(isave,ac,as,bc,bs,kmc,kms)
c
c $Header: /afs/pppl.gov/u/qpliu/cvsroot/codes/gryffin/src/old/nlpscp.f,v 4.0 1995/05/24 04:46:58 mbeer Exp $
c
c Defunct "parallel" version of nlpsc, used before Revision 3.0.
c
c Revision 3.0  1994/02/03  22:17:18  bdorland
c No change from version 2.0.
c
c Revision 2.0  1994/02/03  20:49:13  hammett
c Version 2.0: final release of the cosine/sine version
c
c Revision 1.5  1994/02/03  19:59:36  hammett
c reorder nested loops to get longest loop to vectorize
c
c Revision 1.4  1994/02/01  17:10:23  hammett
c upgrade to use more efficient real-to-complex FFTs.
c
c Revision 1.3  1994/01/29  01:59:06  bdorland
c simplified calling arguments
c
c Revision 1.2  1993/08/30  18:42:59  mbeer
c Bug in dimensions of km.
c
c Revision 1.1  1993/08/17  04:54:42  mbeer
c Initial revision
c
c Revision 1.4  1993/02/02  03:25:52  mbeer
c Removed SAVEs, now can use stack alloc, and put in energy cons. diag.
c
c Revision 1.3  1992/09/03  12:16:25  mbeer
c k_radial calculation was incorrect for k_theta=0.
c
c Revision 3.13  1992/04/01  22:56:12  bdorland
c header and log added.
c
c
c "NonLinear PseudoSpectral" BD
c Given the fields a, b, ka, and kb this subroutine calculates the 
c Poisson bracket
c		{a,b} = da/dx db/dy - da/dy db/dx.
c
c The fields arrive in (x,ky,kz) space and must be transformed
c with the appropriate derivatives into (x,y,z) space.
c
      implicit none
      include 'itg.par'
      include 'itg.cmn'
      integer l,m,n,nindex,isave,imode,mr

      real ac(lz,mz,nz,nspecz),as(lz,mz,nz,nspecz),
     * bc(lz,mz,nz,nspecz),bs(lz,mz,nz,nspecz),
     * kxac(lz,mz,nz,nspecz),kxas(lz,mz,nz,nspecz),
     * kxbc(lz,mz,nz,nspecz),kxbs(lz,mz,nz,nspecz),
     * kyac(lz,mz,nz,nspecz),kyas(lz,mz,nz,nspecz),
     * kybc(lz,mz,nz,nspecz),kybs(lz,mz,nz,nspecz),
     * kmc(lz,mz,nz,nspecz),kms(lz,mz,nz,nspecz)

      complex kxa1(nzz,mz*lz),kya1(nzz,mz*lz)
      complex kxb1(nzz,mz*lz),kyb1(nzz,mz*lz)

      complex km1(nzz,mz*lz)
      real km(mzz,nzz*lz)

c Don't recalculate FFT's unnecessarily!
c
c isave determines which fields are new for this call:
c isave=0 -> both fields are new
c isave=1 -> "a" fields are the same as last time nlps was called
c isave=2 -> "b" fields are the same as last time nlps was called 
c
c
      nindex=malias*nalias*ldb
c
      if(isave.eq.0) then
      stepnlps=stepnlps+1
      if(mod(stepnlps,2).eq.0) then
      vxmax=vxmax*cflx
      vymax=vymax*cfly
      n=ncycle-1
      write(8,1000) vxmax, vymax, n
      vmax=max(vxmax,vymax)
      vxmax=0.
      vymax=0.
 1000 format('vxmax = ',e16.8,' vymax = ',e16.8,i5)
      endif
      endif
c
c Calculate the d/dy and d/dx terms:
      if(isave.ne.1) then
        do n=1,nd   
           do l=1,ldb
             kyac(l,1,n,1)=0.
             kyas(l,1,n,1)=0.
             kxac(l,1,n,1)=1./y0*shr*r0(1,n)*as(l,1,n,1)
             kxas(l,1,n,1)=-1./y0*shr*r0(1,n)*ac(l,1,n,1)
           enddo
        enddo   

        do m=2,md   
         do n=1,nd
            mr=float(mrr(m,n))
            do l=1,ldb
               kyac(l,m,n,1)=mr/y0*as(l,m,n,1)
               kyas(l,m,n,1)=-mr/y0*ac(l,m,n,1)
               kxac(l,m,n,1)=mr/y0*shr*r0(m,n)*as(l,m,n,1)
               kxas(l,m,n,1)=-mr/y0*shr*r0(m,n)*ac(l,m,n,1)
            enddo
         enddo   
        enddo   
      endif
c
      if(isave.ne.2) then
        do n=1,nd
           do l=1,ldb
             kybc(l,1,n,1)=0.
             kybs(l,1,n,1)=0.
             kxbc(l,1,n,1)=1./y0*shr*r0(1,n)*bs(l,1,n,1)
             kxbs(l,1,n,1)=-1./y0*shr*r0(1,n)*bc(l,1,n,1)
           enddo
        enddo   

        do m=2,md   
         do n=1,nd
            mr=float(mrr(m,n))
            do l=1,ldb
               kybc(l,m,n,1)=mr/y0*bs(l,m,n,1)
               kybs(l,m,n,1)=-mr/y0*bc(l,m,n,1)
               kxbc(l,m,n,1)=mr/y0*shr*r0(m,n)*bs(l,m,n,1)
               kxbs(l,m,n,1)=-mr/y0*shr*r0(m,n)*bc(l,m,n,1)
            enddo
         enddo   
        enddo   
      endif

c
c Transform the vectors i*kx*a, i*kx*b, i*ky*b and i*ky*a into x space:
c
      if(isave.ne.1) then
        call xspacec(kxac,kxas,kxa1)
        call xspacec(kyac,kyas,kya1)
      endif

      if(isave.ne.2) then
        call xspacec(kxbc,kxbs,kxb1)
        call xspacec(kybc,kybs,kyb1)
      endif
c
c Then, transform the same vectors into y-space:
c
      if(isave.ne.1) then
        call yspace(kxa1,kxa3) 
        call yspace(kya1,kya3) 
c
c Check Courant condition:
c
       if(mod(stepnlps,2).eq.0) then
          do imode=1,nalias*ldb
             do m=1,malias  
                vxmax=max(vxmax,abs(kya3(m,imode)))
                vymax=max(vymax,abs(kxa3(m,imode)))
             enddo  
          enddo  
       endif
c
      endif
c
      if(isave.ne.2) then
        call yspace(kxb1,kxb3) 
        call yspace(kyb1,kyb3) 
      endif
c
c Perform the multiplications:
c
      do m=1,malias
         do imode=1,nalias*ldb
            km(m,imode)=kxb3(m,imode)*kya3(m,imode)
     &           -kxa3(m,imode)*kyb3(m,imode)
         enddo
      enddo

c Transform wwac to ky
      call kyspace(km1,km)
c Transform wwac,wwas to kx
      call kxspacec(kmc,kms,km1)
c
c That's all! 
c     
      return
      end