c
      subroutine DiagT(npts,P,ze,Zk,Wm,Wmst,Wm0,Wm0st,f0,Uct,Ucs,
     &		zUct1,zUct2,zUct0,zUcs1,zUcs2,zUcs0,Vsuper)
c
c     DiagT calc the potential contributed by the contact diagram
c***********************************************************************
c
      implicit none
c     LOCAL variables
      integer nch, npts1, ich, jch, cl,clmin,clmax, j, iso
      integer first,ichpt,jchpt,ipt,jpt,l
      double precision pi, zero, one, two, eight, norm, jp
c     SHARED variables
      integer Nifty, npts
      double precision P(64), Wm(64), Wmst(64), A(2,2,2)
      double precision Uct(64,64,0:1), Ucs(64,64,0:1)
      double precision Cplv, Cplct, Cplcs, f0, mpion, Mass
      double complex ze, Zk(2)
      double complex zUct1(2,64,0:1),zUcs1(2,64,0:1)
      double complex zUct2(64,2,0:1),zUcs2(64,2,0:1)
      double complex zUct0(2,2,0:1),zUcs0(2,2,0:1)
      double complex Wm0(2), Wm0st(2), Vsuper(128,128,6)

      parameter (zero=0.0d0, one=1.0d0, two=2.0d0, eight=8.0d0)
      parameter (pi=3.14159265358979323846d0)

      common /control/ Nifty(15)
      common /couple/ Cplv(3,3),Cplct(2,2,2),Cplcs(2,2,2)
      common /mass/ mpion, Mass(2)

      data first/1/
c**********************************
      nch = Nifty(2)
      npts1 = npts + 1
      norm = 1.d0/(eight*pi)
c      lmax = 1
      clmax = Nifty(6)
      clmin = 1
      if (Nifty(4) .eq. 4) clmin = clmax
c
c     Order of PW: P11,P31,P13,P33.
c     add a switch. if nif10=4 switch off space component
c		       nif10=5 switch off time component
      if (Nifty(10).eq.3 .or. Nifty(10).eq.4) then

c ***************************************************
c ***  add time component of contact interaction  ***
c     orbital angular momentum
      do cl = clmin, clmax
	if (cl .le. 2) then
	  l = 0 
	else if (cl .ge. 3) then
	  l = 1 
	endif
c       isospin
	if (cl.eq.1 .or. cl.eq.3 .or. cl.eq.5) then
          iso = 1
        else 
          iso = 2
        endif

c	time component don't contribute to transition elements
         do 101 ich = 1, nch
            do jpt = 1, npts
               jchpt = jpt + (ich-1)*npts1
               do ipt = 1, jpt
                  ichpt = ipt + (ich-1)*npts1
                  Vsuper(ichpt,jchpt,cl) = Vsuper(ichpt,jchpt,cl)
     &            +norm*Cplct(ich,ich,iso)*Uct(ipt,jpt,l)
     &            *(Wm(ipt)+Wm(jpt))/(f0*f0*Wmst(ipt)*Wmst(jpt))
               enddo
               ichpt = ich*npts1
               Vsuper(ichpt,jchpt,cl) = Vsuper(ichpt,jchpt,cl)
     &         +norm*Cplct(ich,ich,iso)*zUct1(ich,jpt,l)
     &         *(Wm0(ich)+Wm(jpt))/(f0*f0*Wm0st(ich)*Wmst(jpt))
            enddo
c
            jchpt = ich*npts1
            do ipt = 1, npts
               ichpt = ipt + (ich-1)*npts1
               Vsuper(ichpt,jchpt,cl) = Vsuper(ichpt,jchpt,cl)
     &         +norm*Cplct(ich,ich,iso)*zUct2(ipt,ich,l)
     &         *(Wm(ipt)+Wm0(ich))/(f0*f0*Wmst(ipt)*Wm0st(ich))
            enddo
            ichpt = ich*npts1
            Vsuper(ichpt,jchpt,cl) = Vsuper(ichpt,jchpt,cl)
     &         +norm*Cplct(ich,ich,iso)*zUct0(ich,ich,l)
     &         *(Wm0(ich)+Wm0(ich))/(f0*f0*Wm0st(ich)*Wm0st(ich))
  101    continue
      enddo

      endif 

      if (Nifty(10).eq.3 .or. Nifty(10).eq.5) then
c ******************************************************
c ***  add spatial component of contact interaction  ***
c	first initialize coef A(ich, jch, j)
c	angular momentu must be 1(p); j=1(1/2) and 2(3/2)
	if (first .eq. 1) then
	  do j = 1, 2
	    jp = j - 0.5d0
	    A(1,1,j) =  2.d0*(0.75 + 2.d0 - jp*(jp+1))
	    A(1,2,j) =  (jp+0.5d0) * sqrt(0.5d0*(jp+3.5d0)*(2.5d0-jp))
	    A(2,1,j) =  A(1,2,j)
	    A(2,2,j) =  -sqrt(0.4d0)*(5.75d0-jp*(jp+1))
	  enddo
	  first = 0
	endif
c

      l = 1
      clmin = 3
      if (nifty(4).eq.4) clmin = clmax
      do cl = clmin, clmax
c	total angular momentum
	if (cl.eq.3 .or. cl.eq.4) then
	  j = 1
	else 
	  j = 2
	endif
c       isospin 
	if (cl.eq.3 .or. cl.eq.5) then
	   iso = 1
	else
	   iso = 2
	endif
c
c
        do 103 jch = 1, nch
        do 104 ich = 1, nch
            do jpt = 1, npts
               jchpt = jpt + (jch-1)*npts1
               do ipt = 1, jpt
                  ichpt = ipt + (ich-1)*npts1
                  Vsuper(ichpt,jchpt,cl) = Vsuper(ichpt,jchpt,cl)
     &            +norm*Cplcs(ich,jch,iso)*Ucs(ipt,jpt,l)
     &            *A(ich,jch,j)/(f0*f0*Wmst(ipt)*Wmst(jpt))
               enddo
               ichpt = ich*npts1
               Vsuper(ichpt,jchpt,cl) = Vsuper(ichpt,jchpt,cl)
     &         +norm*Cplcs(ich,jch,iso)*zUcs1(ich,jpt,l)
     &         *A(ich,jch,j)/(f0*f0*Wm0st(ich)*Wmst(jpt))
            enddo
c
            jchpt = jch*npts1
            do ipt = 1, npts
               ichpt = ipt + (ich-1)*npts1
               Vsuper(ichpt,jchpt,cl) = Vsuper(ichpt,jchpt,cl)
     &         +norm*Cplcs(ich,jch,iso)*zUcs2(ipt,jch,l)
     &         *A(ich,jch,j)/(f0*f0*Wmst(ipt)*Wm0st(jch))
            enddo
            ichpt = ich*npts1
            jchpt = jch*npts1
            Vsuper(ichpt,jchpt,cl) = Vsuper(ichpt,jchpt,cl)
     &         +norm*Cplcs(ich,jch,iso)*zUcs0(ich,jch,l)
     &         *A(ich,jch,j)/(f0*f0*Wm0st(ich)*Wm0st(jch))
  104    continue
  103    continue
	enddo
c
      endif
      return
      end