subroutine nlpsce(isave,a,b,nl_term)
c
c $Header: /afs/pppl.gov/u/qpliu/cvsroot/codes/gryffin/src/old98/nlpsce.f,v 1.1 1998/06/11 15:59:50 bdorland Exp $
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,isave,imode,mr

      complex a(kz,mz,nz),b(kz,mz,nz)
      complex kxa(kz,mz,nz),kya(kz,mz,nz)
      complex kxb(kz,mz,nz),kyb(kz,mz,nz)
      complex nl_term(kz,mz,nz)

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

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

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
      if(isave.eq.0) then
         stepnlps=stepnlps+1
         if(mod(stepnlps,2*nspecies).eq.0) then
            vxmax=vxmax*cflx
            vymax=vymax*cfly
            n=ncycle-1
            vmax=max(vxmax,vymax)
            vxmax=0.
            vymax=0.
         endif
      endif
 1000 format('vxmax = ',e16.8,' vymax = ',e16.8,i5)
c
c Calculate the d/dy and d/dx terms:
      if(isave.ne.1) then
        do n=1,nd   
           do l=1,kd
              kya(l,1,n)=0.
              kxa(l,1,n)=ii/y0*shr*r0(1,n)*a(l,1,n)
           enddo
        enddo   
        
        do m=2,md   
           do n=1,nd
              mr=float(mrr(m,n))
              do l=1,kd
                 kya(l,m,n)=ii*mr/y0*a(l,m,n)
                 kxa(l,m,n)=ii*mr/y0*shr*r0(m,n)*a(l,m,n)
              enddo
           enddo   
        enddo   
      endif
c     
      if(isave.ne.2) then
         do n=1,nd
            do l=1,kd
               kyb(l,1,n)=0.
               kxb(l,1,n)=ii/y0*shr*r0(1,n)*b(l,1,n)
            enddo
         enddo   
         
         do m=2,md   
            do n=1,nd
               mr=float(mrr(m,n))
               do l=1,kd
                  kyb(l,m,n)=ii*mr/y0*b(l,m,n)
                  kxb(l,m,n)=ii*mr/y0*shr*r0(m,n)*b(l,m,n)
               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 xspacece(kxa,kxa1)
         call xspacece(kya,kya1)
      endif
      
      if(isave.ne.2) then
         call xspacece(kxb,kxb1)
         call xspacece(kyb,kyb1)
      endif
c     
c     Then, transform the same vectors into y-space:
c
      if(isave.ne.1) then
         call yspacee(kxa1,kxae3) 
         call yspacee(kya1,kyae3) 
c     
c     Check Courant condition:
c     
         if(mod(stepnlps,2*nspecies).eq.0) then
            do imode=1,nalias*kd
               do m=1,malias  
                  vxmax=max(vxmax,abs(kyae3(m,imode)))
                  vymax=max(vymax,abs(kxae3(m,imode)))
               enddo  
            enddo  
         endif
      endif
c     
      if(isave.ne.2) then
         call yspacee(kxb1,kxbe3) 
         call yspacee(kyb1,kybe3) 
      endif
c
c     Perform the multiplications:
c     
      do m=1,malias
         do imode=1,nalias*kd
            km(m,imode)=kxae3(m,imode)*kybe3(m,imode)
     &           -kxbe3(m,imode)*kyae3(m,imode)
         enddo  
      enddo
      
c  Transform wwac to ky
      call kyspacee(km1,km)
c  Transform wwac,wwas to kx
      call kxspacece(nl_term,km1)
c     
c     That's all! 
c     
      return
      end