subroutine waccell(psi, linfac, source, omega)
        use ellipdim
        if (wbypass) return
        werror = .false.
      do 10001 i001 = 1, (  81 -1), 1
      do 10000 i000 = 1, (  81 -1), 1
      basis(i000,i001,1) = psi(i000,i001)
      basis(i000,i001,2) = psi(i000,i001) - psialt(i000,i001,1)
      basis(i000,i001,3) = psi(i000,i001) - 2*psialt(i000,i001,1) + psia
     *lt(i000,i001,2)
      basis(i000,i001,4) = 1
10000 continue
10001 continue
      call periodic( mx, my,  basis1  )
      call periodic( mx, my,  basis2  )
      call periodic( mx, my,  basis3  )
      call periodic( mx, my,  basis4  )
      do 12500 i = 1, wdim
      ii = i
      do 12501 j = i, wdim
      jj = j
           call makematl(psi, linfac, source, omega, i, j)
12501 continue
12500 continue
      do 12502 i = 1, wdim
           wsource(i) = 0
      do 10003 i001 = 1, (  81 -1), 1
      do 10002 i000 = 1, (  81 -1), 1
            wsource(i) = source(i000,i001)*basis(i000,i001,i) + wsource(
     *i)
10002 continue
10003 continue
12502 continue
        call linsys(wmatrix, wdim, wdim, wsource, wcoeff, ising, lfirst,
     *              lprint, work)
        if (ising .eq.  1) then
           write(*,*) ' WARNING:  W_matrix is singular '
           werror = .true.
           return
        endif
      do 10005 i001 = 1, (  81 -1), 1
      do 10004 i000 = 1, (  81 -1), 1
         psi(i000,i001) = 0
10004 continue
10005 continue
      do 12503 i = 1, wdim
      do 10007 i001 = 1, (  81 -1), 1
      do 10006 i000 = 1, (  81 -1), 1
            psi(i000,i001) = psi(i000,i001) + wcoeff(i)*basis(i000,i001,
     *i)
10006 continue
10007 continue
12503 continue
      do 12504 i = 1, wdim
           write(*,100) i, .5*wmatrix(i,i) - wsource(i),
     *                  i, wcoeff(i)
12504 continue
        action = 0
      do 12505 i = 1, wdim
      do 10008 i000 = 1, ( wdim ), 1
            action = action + wmatrix(i,i000)*wcoeff(i)*wcoeff(i000)
10008 continue
12505 continue
        action = action/2
      do 10009 i000 = 1, ( wdim ), 1
         action = action - wsource(i000)*wcoeff(i000)
10009 continue
        write(*,*) ' new action = ',action
        return
100     format(' action(',i1')= ',g16.9,'    w_coeff(',i1,')= ', g16.9)
        end