subroutine WtgFcn(npts, P, W, ze, Zk, Zk2, Dk2de, Wgf)

c***********************************************************************
c     subroutine WtgFcn calculates the WeighTed Green's FunCtioN at
c        each grid points
c     See G. He's note 'From Integral Equation to Matrix Equation',
c        Nov. 13, 14, 1990. 
c     In folder 'Notes on Scattering Theory'
c     wgf (weighted green's function) is same as D_{i,gamma} in 
c        the notes
c***********************************************************************

      implicit none

      integer npts1,ipt,nch,ich,ichpt
      double precision Egama(64,2), zero, two, pi

      integer Nifty,npts
      double precision P(64), W(64), mpion, Mass
      double complex zi, ze, Zk(2), Zk2(2), Dk2de(2), Wgf(128)

      common /control/ Nifty(15)
      common /mass/ mpion, Mass(2)

      parameter (zero = 0.0d0, two = 2.0d0)
      parameter (pi=3.14159265358979323846d0)
      parameter (zi = (0.0d0, 1.0d0))
c************************************

      nch = Nifty(2)
      npts1 = npts + 1
c     relativistic kinematics
         do 10010 ich = 1, nch
            do ipt = 1, npts
               Egama(ipt, ich) = sqrt(P(ipt)**2 + mpion**2)
     &                       +sqrt(P(ipt)**2 + Mass(ich)**2)
            enddo
10010    continue

c calculate wgf on grid points
      do 10050 ich = 1, nch
         do ipt = 1, npts
            ichpt = (ich-1)*npts1+ipt
            Wgf(ichpt) = two*W(ipt)*P(ipt)*P(ipt)
     &                  /(pi*(ze-Egama(ipt,ich)))
	 enddo
10050 continue

c calculate wgf on on-shell points
      do 10070 ich = 1, nch
         ichpt = ich*npts1
         Wgf(ichpt) = zero
c         if ( ((dble(ze)-mpion-Mass(ich)).gt.zero) .or. 
c     &   (imag(Zk(ich)).lt.zero) ) then
	if ( (dble(ze)-mpion-Mass(ich)).gt.zero ) then

c	3/9 lu change .le. to .lt. to cure delta channel below md-mpi.
c	i.e. for above threshold or decaying resonant state
            do ipt = 1, npts
               Wgf(ichpt) = Wgf(ichpt) + W(ipt)/(Zk2(ich)-P(ipt)**2)
	    enddo
            Wgf(ichpt) = -Wgf(ichpt)*two*Zk2(ich)*Dk2de(ich)/pi
c     &           -zi*Zk(ich)*Dk2de(ich)
c           if im(k) !=0 iepslon prescription irrelevent
	    if (imag(Zk(ich)) .le. zero) then
	       Wgf(ichpt) = Wgf(ichpt) - zi*Zk(ich)*Dk2de(ich)
	    else
	       Wgf(ichpt) = Wgf(ichpt) + zi*Zk(ich)*Dk2de(ich)
	    endif
         endif
10070 continue

      return
      end