subroutine bcon(a)
c
c Enforce boundary conditions in the x (l) direction 
c
      implicit none
      include 'itg.par'
      include 'itg.cmn'
      complex a(lz,mz,nz)
      real beta1,znshift,wgt1,wgt2,wwidth
      integer i,l,m,n,ninv(-nz:nz),nshift,n_max,mr,nr

      if(iperiod.eq.0) then
         do m=1,md
            do n=1,nd
               a(1,m,n)=0.
               a(ld,m,n)=0.
            enddo
         enddo
      endif

      if (iperiod.eq.1) then
         do m=1,md
            do n=1,nd
               a(ld,m,n)=a(1,m,n)
            enddo
         enddo
      endif

c Twisted periodicity

      if (iperiod.eq.2) then
         n_max=(nd-1)/2

         do n=1,nd
            do nr=-n_max,n_max
               if (nrr(n).eq.nr) ninv(nr)=n
            enddo
         enddo

          znshift=2.*n_max*(xp+.0001)/z0
c handle the special case of a single n mode (nd=1, which
c might be used for linear runs) properly:
         if(n_max .eq. 0) znshift=1

         do m=1,md
            do n=1,nd
               nr=nrr(n)
               mr=mrr(m,n)
               beta1=2.*mr*xp*n0*qsf
               nshift=mr*znshift

               do l=1,l_left-1
                  if ((mr.eq.0).or.((abs(nr+nshift).le.n_max)
     >                  .and.(nshift.gt.0)))  then
                     a(l,m,n)=a(l+l_right-l_left,m,
     >                    ninv(nr+nshift))*exp(-ii*beta1)
                  else
                     a(l,m,n)=0.0
                  endif
               enddo
               do l=l_right,ld
                  if ((mr.eq.0).or.((abs(nr-nshift).le.n_max)
     >                  .and.(nshift.gt.0)))  then
                     a(l,m,n)=a(l-l_right+l_left,m,
     >                    ninv(nr-nshift))*exp(ii*beta1)
                  else
                     a(l,m,n)=0.0
                  endif
               enddo
            enddo
         enddo

c Window to filter near the edges of the box.
c
c The window weight function is 1 in most of the box, but near the end
c of the extension goes like 2x**2/(1+x**4), where x is the distance
c from the edge of the box, and is normalized to match smoothly to the
c weight=1 region.
c
c Another option might be to use sine transforms...
c
c wwidth is the fraction of the extension region over which the window
c is dropping to zero:
         wwidth=0.5

c (Careful, skip m=0 mode, which is supposed to be periodic with itself...).
         if(md .ge. 2) then
         do m=2,md
            do n=1,nd
               do l=1,l_left-1

c wgt2 is a linear function, 0 at the end, and reaches a maximum of 1:
                  wgt2=(l-1)/float(l_left-1)/wwidth
c                 wgt2=1.0  ! to revert to no window
                  if(wgt2 .ge. 1.0) wgt2=1.0
                  wgt1=2*wgt2**2/(1+wgt2**4)
                  a(l,m,n)=a(l,m,n)*wgt1
               enddo
c               a(1,m,n)=0.0
c               a(2,m,n)=a(2,m,n)*0.5

               do l=l_right,ld
                  wgt2=(ld-l)/float(ld-l_right)/wwidth
c                 wgt2=1.0   ! to revert to no window
                  if(wgt2 .ge. 1.0) wgt2=1.0
                  wgt1=2*wgt2**2/(1+wgt2**4)
                  a(l,m,n)=a(l,m,n)*wgt1
               enddo
c               a(ld,m,n)=0.0
c               a(ld-1,m,n)=a(ld-1,m,n)*0.5
            enddo
         enddo
         endif

      endif
c
c     for connected FFTs, zero leftmost point of leftmost mode
c     and rightmost point of rightmost mode, except for m=1 (ky=0) 
c
      if(iperiod.eq.3) then
         do m=2,md
            do n=1,nd
               if (conpos(m,n).eq.1) a(1,m,n)=0.
               if (conpos(m,n).eq.icon(m,n)) a(ld,m,n)=0.
            enddo
         enddo
         do n=1,nd
            a(ld,1,n)=a(1,1,n)
         enddo   
      endif

c enforce reality condition (just to be sure):

      if(nd.ne.1 .and. md.ne.1) then
         do n=1,nddm
            do l=1,ldb
               a(l,1,nddp+n)=conjg(a(l,1,nddp+1-n))
            enddo
         enddo

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

      return
      end