subroutine nlpsc(isave,a,b,nl_term)
c
c     (c) Copyright 1991 to 1995 by Michael A. Beer, William D. Dorland, 
c     and Gregory W. Hammett. ALL RIGHTS RESERVED.
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'
      real mr
      integer l,m,n,isave,imode
      
c      real dum1,dum2,dum3,dummy(mz,nz)

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

      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
c      if(isave.eq.0) then
c         stepnlps=stepnlps+1
c         if(mod(stepnlps,2*nspecies).eq.0) then
c            vxmax=vxmax*cflx
c            vymax=vymax*cfly
c            vmax=max(vxmax,vymax)
c            vxmax=0.
c            vymax=0.
c         endif
c      endif
c
c Calculate the d/dy and d/dx terms:
      if(isave.ne.1) then
        do n=1,nd   
           do l=1,ldb
              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,ldb
                 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,ldb
               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,ldb
                  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 xspacec(kxa,kxa1)
         call xspacec(kya,kya1)
      endif
      
      if(isave.ne.2) then
         call xspacec(kxb,kxb1)
         call xspacec(kyb,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     
c         if(mod(stepnlps,2*nspecies).eq.0) then
c            do imode=1,nalias*ldb
c               do m=1,malias
c                  vxmax=max(vxmax,abs(kya3(m,imode)))
c                  vymax=max(vymax,abs(kxa3(m,imode)))
c               enddo
c            enddo
c         endif
      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)=kxa3(m,imode)*kyb3(m,imode)
     &           -kxb3(m,imode)*kya3(m,imode)
         enddo  
      enddo
      
c  Transform wwac to ky
      call kyspace(km1,km)
c  Transform wwac,wwas to kx
      call kxspacec(nl_term,km1)
c
c As a result of questions from Kerbel:
c
c **************Investigating reality condition; forcing it here*******
c BD 9/15/96
c
c debug:  expect dum3 ~ dum2 << dum1
c
c      call volsq(b,b,dum1,dummy)
c      call volsq(b,nl_term,dum2,dummy)
c     
c enforce reality condition (should be unnecessary, but cheap to do): 
c
      do n=1,nddm
         do l=1,ldb
            nl_term(l,1,nddp+n)=conjg(nl_term(l,1,nddp+1-n))
         enddo
      enddo

c      do l=1,ldb
c         nl_term(l,1,1)=0.
c      enddo

c      call volsq(b,nl_term,dum3,dummy)
c      write(82,*) dum1,dum2,dum3

      return
      end