subroutine bcon(a)
!
! Enforce boundary conditions in the x (l) direction 
!
  use itg_data
  implicit none

  complex, dimension(:,:,:) :: a
  real    :: beta1, znshift, wgt1, wgt2, wwidth
  integer :: l, m, n, nshift, n_max, mr, nr
  integer, dimension(:), pointer, save :: ninv
  integer :: nfirst = 1

  if(.not.associated(ninv)) allocate(ninv(-nd:nd))
  
  if(iperiod == 0) then
     a(1,:,:)=0.
     a(ld,:,:)=0.
  endif
  
  if (iperiod == 1) then
     a(ld,:,:)=a(1,:,:)
  endif

! Twisted periodicity
  
  if (iperiod == 2) then
     n_max=(nd-1)/2
     
     if(nfirst == 1) then
        nfirst = 0
        do n=1,nd
           do nr=-n_max,n_max
              if (nrr(n) == nr) ninv(nr)=n
           enddo
        enddo
     endif
     
     znshift=2.*n_max*(xp+.0001)/z0
! handle the special case of a single n mode (nd=1, which
! might be used for linear runs) properly:
     if(n_max  ==  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 == 0).or.((abs(nr+nshift) <= n_max).and.(nshift > 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 == 0).or.((abs(nr-nshift) <= n_max).and.(nshift > 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

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

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

! wgt2 is a linear function, 0 at the end, and reaches a maximum of 1:
                 wgt2=(l-1)/float(l_left-1)/wwidth
!                 wgt2=1.0  ! to revert to no window
                 if(wgt2  >=  1.0) wgt2=1.0
                 wgt1=2*wgt2**2/(1+wgt2**4)
                 a(l,m,n)=a(l,m,n)*wgt1
              enddo
!               a(1,m,n)=0.0
!               a(2,m,n)=a(2,m,n)*0.5
              
              do l=l_right,ld
                 wgt2=(ld-l)/float(ld-l_right)/wwidth
!                 wgt2=1.0   ! to revert to no window
                 if(wgt2  >=  1.0) wgt2=1.0
                 wgt1=2*wgt2**2/(1+wgt2**4)
                 a(l,m,n)=a(l,m,n)*wgt1
              enddo
!               a(ld,m,n)=0.0
!               a(ld-1,m,n)=a(ld-1,m,n)*0.5
           enddo
        enddo
     endif
     
  endif
!
!     for connected FFTs, zero leftmost point of leftmost mode
!     and rightmost point of rightmost mode, except for m=1 (ky=0) 
!
  if(iperiod == 3) then
     do m=2,md
        do n=1,nd
           if (conpos(m,n) == 1) a(1,m,n)=0.
           if (conpos(m,n) == icon(m,n)) a(ld,m,n)=0.
        enddo
     enddo
     a(ld,1,:)=a(1,1,:)
  endif

! enforce reality condition (just to be sure):

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

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

  return
end subroutine bcon