subroutine voptth(kp, k, k0, cthnuc)
c***********************************************************************
c *** Computes the angle- dependent optical potential 
c      for 1/2 * 1/2 spin scattering
c *** given the amplitudes for nucleon nucleon scattering.
c     voptl  now calcs partial wave VsubL

      implicit real*8 (a-h, i, k, m, o-z)


      dimension amb(4), apb(4), cpd(4), cmd(4), el(4), ff(4)
      dimension nifty(20)
c     dimension v10(2), v1m1(2), v01(2), v11(2), v00(2), vss(2)

c    optical potential values passed to voptll via common/vopt/

      common /vopt/ v10(2), v1m1(2), v01(2), v11(2), v00(2), vss(2),
     $              v1s(2), vsum1, vsum2
      common /params/   hbarc, pi, mpi, xmn, nz, na, nes, nwaves
      common /switch/   nifty

      data ndata, iprint/0, 0/

c >>> FIRST EXECUTABLE STATEMENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

      ndata = 1
      sr2 = dsqrt( 2.d0 )
      mn = xmn
      nn = na - nz
      amass = na * mn/xmn
      aovera = (amass - 1.)/amass
      if (nifty(17) .eq. 1) aovera = 1.

  99  pm1 = kp*kp + k*k
      pm2 = 2. * k * kp
      nff = 4
      if (nifty(6) .eq. 0) nff = 2
c
c *** jjm = 1 for real part protons and neturons
c *** jjm = 2 for imag part protons and neutrons

c *** Tmatrix call
c ***

         call tnucth (kp,k,k0,cthnuc,mn,apb,amb,cpd,cmd,el)

c ***
c *** Form factor call
c ***
         x = cthnuc
         sinth = dsqrt(1.0-x**2)
         sth2 = sinth**2
         s2th2 = 0.5 * (1.-x)
         c2th2 = 0.5 * (1.+x)
         q = dsqrt( dabs(pm1 - pm2 * x) )
         q2 = q*q
c
         call FFACT (q2, ff, nff)
c
         if (nifty(6) .ne. 0) go to 403
c
c *** No spin, set spin ff=0
         ff(3)=0.
         ff(4)=0.
c        change gamin from -30 to -15 2/4/91 tm

 403     gamin = 1.e-15
         ff1 = abs(ff(1))
         ff2 = abs(ff(2))
         ff3 = abs(ff(3))
         ff4 = abs(ff(4))
         if (ff1 .le. gamin) ff(1) = 0.
         if (ff2 .le. gamin) ff(2) = 0.
         if (ff3 .le. gamin) ff(3) = 0.
         if (ff4 .le. gamin) ff(4) = 0.
c
         jt = 0
c
c   set coulomb potential=0. unless exact coulom pot is included
         if (nifty(10) .eq. 3 .or. nifty(10) .eq. 4 ) then
                        vcc = vcoul(q) 
           else		 
                       vcc = 0.
         end if

c *** do loop over real/im parts

         do 11 j = 1,2
            ja = j + jt
            je = j + jt + 1

c *** Modification for cex, interchange n and p

        if (nifty(1) .eq. 9) ja = je
        if (nifty(1) .eq. 9) je = j + jt
c  ****   10/02/90 tm only add vc to real potential
        if (j.eq.1) then 
           vc = vcc
           else 
           vc = 0.0
        endif
        v11a = apb(ja) * nz * ff(1) + apb(je) * nn * ff(2)
        v11b = cpd(ja) * nz * ff(3) + cpd(je) * nn * ff(4)
        v11c = cmd(ja) * nz * ff(3) + cmd(je) * nn * ff(4)
       v11(j) = aovera*(v11a + s2th2*v11b + c2th2 * v11c)+vc
        v10(j) = aovera * (-sinth * (v11b - v11c)/sr2) 
        v1m1a = amb(ja) * nz * ff(3) + nn * amb(je) * ff(4)
        v1m1(j) =aovera*(-v1m1a + c2th2*v11b + s2th2 * v11c)
        v00(j) = aovera*(v11a + v1m1a + x *(v11b - v11c)) +vc
        vss(j) = aovera * (v11a - v1m1a - v11b - v11c) +vc
        jt = 1
  11     continue

c *** To add the term with i = sqrt(-1)
c *** sum1 = real    sum2 = im

         ff1 = 1.
         if (nifty(16) .eq. 6 .or. nifty(16) .eq. 7) ff1 = 0.
         ja = 1
         je = 2
         if (nifty(1) .eq. 9) ja = 2
         if (nifty(1) .eq. 9) je = 1
         vsum1 =(nz * el(ja) * (ff(1) + ff(3) * ff1)
     $         + nn * el(je) * (ff(2) + ff(4) * ff1))*aovera
         vsum2 =(nz*el(ja+2) * (ff(1) + ff(3) * ff1)
     $         + nn*el(je+2) * (ff(2) + ff(4) * ff1))*aovera
************ tm 9/19/91  ****************************
         v01(1) = v10(1) - vsum2 /sr2
         v01(2) = v10(2) + vsum1 /sr2
         v10(1) = v10(1) + vsum2 /sr2
         v10(2) = v10(2) - vsum1 /sr2

c        F-like contrib to optical poten, rhl 6/89
c        F-term turned off if nifty(18) = 4

         if (nifty(18).ne.4) then
            v1s(2) =  aovera*(nz * el(ja) * (ff(1) - ff(3) * ff1)
     $         + nn * el(je) * (ff(2) - ff(4) * ff1))/dsqrt(2.d0)
            v1s(1) =-aovera*(nz * el(ja+2) * (ff(1) - ff(3) * ff1)
     $       + nn * el(je+2) * (ff(2) - ff(4) * ff1))/dsqrt(2.d0)
         else
            v1s(1) = 0.
            v1s(2) = 0.
         endif

c        print only for on-shell point
c        use iprint as switch, so call from main does not print

         if (kp .ne. k .or. kp .ne. k0)go to 20
         if (iprint .lt. 0)go to 20
         iprint = 1
         if (cthnuc.le.-0.999)iprint = -1
c
         do 29 j = 1,4
 29        ff(j) = ff(j)**2
c
         q2f = q2/hbarc/hbarc
         theta = 180. * acos(cthnuc)/pi
         sumia = 0.
c
         do 28 j = 1,2
            sumia = sumia + 0.5 * v11(j)**2 + 0.25 * v00(j)**2
     $             + 0.5 *  v01(j)**2 + 0.5 * v10(j)**2
     $             + 0.25 * vss(j)**2 + 0.5 * v1m1(j)**2
     $             + 0.5 *  v1s(j)**2
 28      continue

 20   continue

c *** FORMATS ***

 96   format(1h ,2x,'th ',6x,'q2f',4x,'f1**2',4x,'f2**2',
     $ 4x,'f3**2',4x,'f4**2',5x,'vc',6x,'v11',6x,
     $'v10',7x,'v1m1',6x,'v01',6x,'v00',7x,'vss',6x,'v1s')
 25   format(1x,f5.1,f8.3,4e12.5,1x,8e13.5/62x,7e13.5)
 119  format(e12.4,'    ',e12.4)
      RETURN

      end