subroutine nlpm(isave,ac,as,kac,kas,bc,bs,kbc,kbs,wwac11,wwas11)
c
c "NonLinear Phase-Mixing" BD
c Given the fields a, b, ka, and kb this subroutine calculates the 
c modified Poisson bracket
c	      {a,b} =  | da/dx d/dy  -  da/dy d/dx | b  
c
c approximated as
c             {a,b} =  | da/dx |  | d/dy | b + | da/dy |  | d/dx | b  
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
      real ac(lz,mz,nz,nspecz),as(lz,mz,nz,nspecz),
     * bc(lz,mz,nz,nspecz),bs(lz,mz,nz,nspecz),
     * kac(lz,mz,nz,nspecz),kas(lz,mz,nz,nspecz),
     * kbc(lz,mz,nz,nspecz),kbs(lz,mz,nz,nspecz),
     * kyas(lz,mz,nz,nspecz),kyac(lz,mz,nz,nspecz),
     * kybc(lz,mz,nz,nspecz),kybs(lz,mz,nz,nspecz)

      real kac21(lz,mz,nzz),kas21(lz,mz,nzz),
     * kbc21(lz,mz,nzz),kbs21(lz,mz,nzz),
     * kyac21(lz,mz,nzz),kyas21(lz,mz,nzz),
     * kybc21(lz,mz,nzz),kybs21(lz,mz,nzz)

      real wwac22(lz*mzz*nzz),
     * wwac21(lz,mz,nzz),wwas21(lz,mz,nzz),
     * wwac11(lz,mz,nz,nspecz),wwas11(lz,mz,nz,nspecz)

      real mr,avg,volum

c     Calculate the | d/dy | terms:

      write(*,*) 'AAUUGGHH      !  I should not be in nlpm!'

      if(isave.ne.1) then
      do 10 n=1,nd
      do 10 m=1,md
      mr=float(mrr(m,n))
      do 10 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)
 10   continue
      endif
c
      if(isave.ne.2) then
      do 20 n=1,nd
      do 20 m=1,md
      mr=float(mrr(m,n))
      do 20 l=1,ldb
      kybc(l,m,n,1)=mr/y0*bc(l,m,n,1)
      kybs(l,m,n,1)=mr/y0*bs(l,m,n,1)
 20   continue
      endif
c
c Transform the vectors kx*a, kx*b, ky*b and ky*a into z space:
c
      if(isave.ne.1) then
c      call zspace(kac,kas,kac21,kas21)
c      call zspace(kyac,kyas,kyac21,kyas21)
      endif
c
      if(isave.ne.2) then
c      call zspace(kbc,kbs,kbc21,kbs21)
c      call zspace(kybc,kybs,kybc21,kybs21)
      endif
c
c Then, transform the same vectors into y-space:
c
      if(isave.ne.1) then
      call yspace(kac21,kas21,kac22)
      call yspace(kyac21,kyas21,kyac22)
      endif
c
      if(isave.ne.2) then
      call yspace(kbc21,kbs21,kbc22)
      call yspace(kybc21,kybs21,kybc22)
      endif
c
c Perform the multiplications:
c
      nindex=malias*nalias*ldb
      volum=1./(float(nindex))
      avg=0.
      do 1 n=1,nindex
      wwac22(n)=(kbc22(n)*abs(kyac22(n))+abs(kac22(n))*kybc22(n))
     *  *float(malias)
      avg=avg+wwac22(n)*volum
 1    continue
c
c Subtract the volume average of the nonlinear term:
c
      do 2 n=1,nindex
      wwac22(n)=wwac22(n)-avg
 2    continue
c
c Transform wwac to ky
      call kyspace(wwac22,wwac21,wwas21)
c Transform wwac,wwas to kz
c      call kzspace(wwac21,wwas21,wwac11,wwas11)
c     
c That's all!
c              not quite...
c Dealias in the x-direction as well:
c
c      call kxspace(wwac11,wwas11)
c
      return
      end