c$debug
      SUBROUTINE POTENT(AKIN,U,V,K ,KMAX,icode
     1          ,LAM,LDFRM,ABETA,IERR,icon19)
c
      character*18 b(0:40), afx
      character*3 Anti
      integer tspin, sp_code
      DIMENSION U(800),V(800),UT(10),G(4)
     1,        AKIN(32),BKIN(32),CN(16),YLAM(16),X( 8),W( 8)
     2,        CP(16),LDFRM(3),ABETA(3),BETA(3),FACTR(3),D(10)
C
      EQUIVALENCE (G(1),FN),(G(2),FL),(G(3),FJ2),(G(4),FSS)
      EQUIVALENCE (EN,BKIN(1)), (FM,BKIN(2)), (Z,BKIN(3)), (FMA,BKIN(4))
     1,(ZA,BKIN(5)), (RY,BKIN(6)), (FS,BKIN(7)), (ECM,BKIN(8))
     2,(FMU,BKIN(9)), (RC,BKIN(10)), (FK2,BKIN(11)), (FK,BKIN(12))
     3,(ETA,BKIN(13)), (AFACT,BKIN(14)), (PFACT,BKIN(15)),(DRF,BKIN(16))
     4,(PNLOC,BKIN(18)), (vfact,BKIN(20)), (sfact,bkin(21))
     5,(hbarc,bkin(27)), (BIF,BKIN(28)), (G,BKIN(29))
c
	common/ch3files/input, ioutp
c
      DATA b                 /' NX=0 No potential', ' NX=1 VOLUME  W-S '
     1, ' NX=2 SURFACE W-S ', ' NX=3 2ND DERIV   ', ' NX=4 L.S VOLUME  '
     2, ' NX=5 L.S SURFACE ', ' NX=6 VOL*R**POWR ', ' NX=7 SUR*R**POWR '
     3, ' NX=8 EXTERN FORMF', ' NX=9 HARMONIC OSC', ' NX=10 GAUSSIAN   '
     4, ' NX=11 DEFORM VOL ', ' NX=12 DEFORM SURF', ' NX=13 Hulthen    '
     5, ' NX=14 Yukawa Lam ', ' NX=15 Yukawa L*S ', ' NX=16 No option  '
     6, ' NX=17 No option  ', ' NX=18 No option  ', ' NX=19 No option  '
     7, ' NX=20 VECTOR     ', ' NX=21 SCALAR     ', ' NX=22 No option  '
     8, ' NX=23 No option  ', ' NX=24 No option  ', ' NX=25 No option  '
     9, ' NX=26 No option  ', ' NX=27 No option  ', ' NX=28 No option  '
     $, ' NX=29 No option  ', ' NX=30 No option  '
     $, ' NX=31 n_nbar V_r ', ' NX=32 n_nbar W_i ', ' NX=33 n_nbar VT  '
     $, ' NX=34 n_nbar VSO ', ' NX=35 n_nbar WS0 '
     $, ' NX=36 l_lbar V_r ', ' NX=37 l_lbar W_i ', ' NX=38 l_lbar VT  '
     $, ' NX=39 l_lbar VSO ', ' NX=40 l_lbar WSO '
     $/
C
      DATA SQRPI, ETA4, ETA5, ETA6, Emax
     1     /1.772453851, 4.0, 10.0, 60.0, 32./
C
      RM=0.0
      sp_code = icode-1024*(icode/1024)
      ISFLAG=0
      IVFLAG=0
      if(icon19.eq.0) then
      WRITE(ioutp,9008)
      endif
      DO 45 M=1,32
      BKIN(M)=AKIN(M)
   45 CONTINUE
      fact=vfact
      DR=DRF
      KFLAG=0
C
C     READ IN CARD SET 5,6,OR 7   POTENTIAL CARDS
C
   70 CONTINUE
      READ (input,9000)FZ,VR,RY,AR,VSOR,VI,RZ,AI,VSOI,PWR
      NX=ABS(FZ)
      IF(RY.GE.0.0) RR=    RY * AFACT
      IF(RY.LT.0.0) RR=ABS(RY)*(AFACT+PFACT)
      IF(RZ.GE.0.0) RI=    RZ * AFACT
      IF(RZ.LT.0.0) RI=ABS(RZ)*(AFACT+PFACT)
   74 CONTINUE
c     Sets up hyperon-hyperon_bar potential options
      ispin = 0
      tspin = 0
        Anti = '   '
        afx  = '                  '
        ncode = (abs(fz)-float(nx))*10000. + 0.5
      if(nx .eq. 10 .and. ncode .ne. 0) then
        afx = '  Spin orbit       '
      endif
      if(nx .ge. 31 .and. nx .le. 40) then
c
c       spin and isospin numbers coded as FZ.st, st = 00, 10, 01 and 11
c	Potential option, kbb, is encoded in 4th position,  .st0k
c								^
c
        ispin = ncode/1000
        tspin = ncode/100  - 10*ispin
	kbb   = mod(ncode,10)
        Anti(1:1) = char(ispin+48)
        Anti(2:2) = char(tspin+48)
	afx( 1:10)= '  kbb =   '
	afx(11:11)= char(kbb+48)
c	write(ioutp,'(a3,i3)')' spin, tspin, kbb =',ispin,tspin,kbb
      endif
c--------------------------------------------------------------------
      if(icon19 .eq. 0) then
      if(abs(vr).lt.1000.) then
        WRITE(ioutp,9509)b(NX),Anti,VR,RY,AR,RR,VSOR
      else
        WRITE(ioutp,9511)b(NX),Anti,VR,RY,AR,RR,VSOR
      endif
      if(abs(vi).lt.1000) then
        WRITE(ioutp,9510)afx  ,Anti,VI,RZ,AI,RI,VSOI,PWR
      else
        WRITE(ioutp,9512)afx  ,Anti,VI,RZ,AI,RI,VSOI,PWR
      endif
      endif
      IF(NX.EQ.0) GO TO 6050
   75 CONTINUE
C
      KFLAG=0
      VR=VR*FACT
      VI=VI*FACT
C     sp_code = 1,2 FOR SINGLE PARTICLE FORM FACTOR
      IF((sp_code.EQ.1 .OR. sp_code.EQ.2) .and. fk2.ne.0.0) THEN
        RM= MAX (RR,RM)
        ALPHB=2.3*ETA4
        TEMP=1.0+ALPHB/10.0
        KT=((ALPHB-ETA*TEMP)/FK+RR)/DR
      ELSE
        KT=FK* MAX (RR,RI)+ETA5
        LPLUS=MAX0(LPLUS,KT)
        KT=(2.3*ETA4* MAX (AR,AI)+ MAX (RR,RI))/DR
      ENDIF
      K=MIN0(MAX0(K,KT),KMAX)
   80 CONTINUE
c
c     This section prevents overflow for small AR and AI
c
      if(NX .eq. 1 .or. NX .eq. 4) then
        M1 = 0
        M2 = KMAX
        if(NX .eq. 1 .and. AR .lt. DR/Emax) then
          AR = DR/Emax
        endif
        if(NX .eq. 4 .and. AR .lt. 2.0*DR ) then
          AR = 2.0*DR
        endif
c	write(20,'(a, 2i4,1p3e12.4)')' M1, M2, RR, AR',M1,M2,RR,AR,RR/AR
        if(RR/AR .ge. Emax) then
          M1 = RR/DR -1.5
          M2 = M1 + 4
          RR = RR - float(M1)*DR
        endif
c
        M3 = 0
        M4 = KMAX
        if(AI .lt. DR/Emax) then
          AI = DR/Emax
        endif
        if(RI/AI .ge. Emax) then
          M3 = RI/DR -1.5
          M4 = M3 + 4
          RI = RI - float(M3)*DR
        endif
        ENDIF
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      if(AR .gt. 0.0) then
c	write(20,'(a, 2i4,1p3e12.4)')' M1, M2, RR, AR',M1,M2,RR,AR,RR/AR
        F2 = EXP(-DR/AR)
        F1 = EXP( RR/AR)
      else
        F2 = 0.0
        F1 = 0.0
      endif
      if(AI .gt. 0.0) then
        F4=EXP(-DR/AI)
        F3=EXP( RI/AI)
      else
        F4 = 0.0
        F3 = 0.0
      endif
C
      IF(NX.le.5 .and. LAM.ne.0) THEN
      IF(icode.EQ.0) PWR=PWR+1.0
      IF(PWR .NE. 0.0) then
        IF(AR .GT. 0.0) VR=VR*(RR/AR)**PWR
        IF(AI .GT. 0.0) VI=VI*(RI/AI)**PWR
      endif
      ENDIF
c
      IF(NX.EQ.0) GO TO 6050
      GO TO ( 100, 200, 300, 400, 500, 600, 700, 800, 900,1000),NX
      GO TO (1100,1200,1300,1400,1500,1600,1700,1800,1900,2000),NX-10
      GO TO (2100),NX-20
c
c     n_nbar potentials, Ueda - Weise, Haidenbauer, Machleidt
c	for kbb = 0, 1 and 2
      if(nx.gt.30.and.nx.le.35) then
        inf = 0
      endif
      GO TO (3100,3200,3300,3400,3500),NX-30
c
c     l_lbar potentials, Weise
c
      if(nx.gt.35.and.nx.le.40) then
        inf = 1
      endif
      GO TO (3100,3200,3300,3400,3500),NX-35
      write(ioutp,'(''0Invalid potential option, Fx ='',i3)')nx
      go to 6000
c
c     ********************************************************************
C
C     VOLUME WOODS SAXON
C
  100 CONTINUE
      DO 110 M=1,M2
      MK=M+M-1
      if(M .gt. M1) then
        F1=F1*F2
        U(MK  )=U(MK  )-VR*F1/(1.0+F1)
      else
        U(MK  )=U(MK  )-VR
      endif
  110 CONTINUE
c
      DO 140 M=1,M4
      MK=M+M
      if(M .gt. M3) then
        F3=F3*F4
        U(MK  )=U(MK  )-VI*F3/(1.0+F3)
      else
        U(MK  )=U(MK  )-VI
      endif
  140 CONTINUE
      GO TO 6000
C
C     1ST DER WOODS SAXON
C
  200 CONTINUE
      DO 260 M=1,K
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )+VR*F1/(1.0+F1)**2
      U(MK+1)=U(MK+1)+VI*F3/(1.0+F3)**2
  260 CONTINUE
      GO TO 6000
C
C     2ND DER WOODS SAXON
C
  300 CONTINUE
      DO 360 M=1,K
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )-VR*F1*(1.0-F1)/(1.0+F1)**3
      U(MK+1)=U(MK+1)-VI*F3*(1.0-F3)/(1.0+F3)**3
  360 CONTINUE
      GO TO 6000
C
C     L.S VOLUME WOODS SAXON
C
  400 CONTINUE
      IF(AR.NE.0.0) VR=VR/AR
      IF(AI.NE.0.0) VI=VI/AI
      R=0.0
      DO 460 M=1,K
      MK=M+M-1
      R=R+DR
      F1=F1*F2
      F3=F3*F4
      V(MK  )=V(MK  )-VR*(F1/(R*(1.0+F1)))/(1.0+F1)
      V(MK+1)=V(MK+1)-VI*(F3/(R*(1.0+F3)))/(1.0+F3)
  460 CONTINUE
      GO TO 6000
C
C     L.S 1ST DER WOODS SAXON
C
  500 CONTINUE
      IF(AR.NE.0.0) VR=VR/AR
      IF(AI.NE.0.0) VI=VI/AI
      R=0.0
      DO 560 M=1,K
      MK=M+M-1
      R=R+DR
      F1=F1*F2
      F3=F3*F4
      V(MK  )=V(MK  )+VR*F1*(1.0-F1)/(R*(1.0+F1)**3)
      V(MK+1)=V(MK+1)+VI*F3*(1.0-F3)/(R*(1.0+F3)**3)
  560 CONTINUE
      GO TO 6000
C
C     WOOD SAXON*R**PWR
C
  600 CONTINUE
      R=0.0
      DO 660 M=1,K
      R=R+DR
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      RPWR=R**PWR
      U(MK  )=U(MK  )-VR*F1*RPWR/(1.0+F1)
      U(MK+1)=U(MK+1)-VI*F3*RPWR/(1.0+F3)
  660 CONTINUE
      GO TO 6000
C
C      1ST DER WOOD SAXON*R**PWR
C
  700 CONTINUE
      R=0.0
      DO 760 M=1,K
      R=R+DR
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      RPWR=R**PWR
      U(MK  )=U(MK  )+VR*F1*RPWR/(1.0+F1)**2
      U(MK+1)=U(MK+1)+VI*F3*RPWR/(1.0+F3)**2
  760 CONTINUE
      GO TO 6000
C
C     EXTERNAL FORM FACTOR
C
  800 CONTINUE
      READ(input,9000)F1,F2
C
C     F2 = 0  REAL CENTRAL
C     F2 = 1  IMAG CENTRAL
C     F2 = 2  REAL SPIN ORBIT
C     F2 = 3  IMAG SPIN ORBIT
C     F2 = 4  REAL and IMAG CENTRAL
C
      if2 = F2
      IF    (if2.EQ.0 .OR. if2.EQ.2) THEN
        F3=VR
        MK=1
      ELSEIF(if2.EQ.1 .OR. if2.EQ.3) THEN
        F3=VI
        MK=2
      ELSEIF(if2.EQ.4)               THEN
        F3=VR
        F4=VI
        MK=1
        ELSE
        STOP 'Error in Potent, incorrect option 8 parameter'
      ENDIF
c
      IF(F3.EQ.0.0) F3=1.0
      IF(F4.EQ.0.0) F4=1.0
      KT=F1
      DO 820 M=1,KT,5
      if(if2 .lt. 4) then
        READ (input,9100)(UT(I),I=1,5)
      else
        READ (input,9100)(UT(I),I=1,10)
      endif
      IF(MK.le.2*KMAX-1) then
      ifinal = min0(5,kt-m+1)
      DO 810 I=1,ifinal
      IF    (if2.EQ.0 .OR. if2.EQ.1) THEN
        U(MK  )=U(MK  )+UT(  I  )*F3
      ELSEIF(if2.EQ.2 .OR. if2.EQ.3) THEN
        V(MK  )=V(MK  )+UT(  I  )*F3
      ELSEIF(if2.EQ.4)               THEN
        U(MK  )=U(MK  )+UT(2*I-1)*F3
        U(MK+1)=U(MK+1)+UT(2*I  )*F4
      ENDIF
      MK=MK+2
  810 CONTINUE
      endif
  820 CONTINUE
C
      K=MIN0(MAX0(K,KT),KMAX)
      GO TO 6000
C
C     HARMONIC OSCILLATOR
C
  900 CONTINUE
      READ (input,9000)G
      if(icon19.eq.0) then
      WRITE(ioutp,9508)G
      endif
      F1=1.0/RY**2
      F2=F1/RY
      F3=0.5
      F4=SQRPI*0.5
      L=FL
      IF(L.EQ.0) GO TO 931
      DO 930 I=1,L
      F3=F3+1.0
      F4=F4*F3
      F2=F2*F1
  930 CONTINUE
  931 CONTINUE
      NN=FN
      T1=1.0
      T2=F4
      CN(1)=(-1.0)**NN
      IF(NN.EQ.0) GO TO 941
      DO 940 I=1,NN
      F3=F3+1.0
      T1=T1*FLOAT(I)
      T2=T2*F3
      CN(I+1)=-CN(I)*F1*FLOAT(NN+1-I)/(FLOAT(I)*F3)
  940 CONTINUE
  941 CONTINUE
      T2=SQRT(2.0*F2*T2/T1)/F4
      IF(VR.NE.0.0) T2=T2*VR
      KT=10.0*RY/DR
      K=MIN0(MAX0(K,KT),KMAX)
      R=0.0
      F1=F1/2.0
      DO 960 M=1,KT
      MK=M+M-1
      R=R+DR
      R2=R**2
      F2=CN(1)
      IF(NN.EQ.0) GO TO 951
      F3=1.0
      DO 950 I=1,NN
      F3=F3*R2
      F2=F2+CN(I+1)*F3
  950 CONTINUE
  951 CONTINUE
      U(MK  )=U(MK  )+F2*T2*EXP(-F1*R2)*R**L
  960 CONTINUE
      rm=ry
      K=MIN0(KT,KMAX)
      GO TO 6000
C
C     EXP(-(R/R0)**2)*R**POWR
C
 1000 CONTINUE
      if2 = ncode/1000
      if(VR .ne. 0.0) then
      R=0.0
      KT=10.0*RY/DR
      K=MIN0(MAX0(K,KT),KMAX)
      DO 1060 M=1,KT
      MK=M+M-1
      R=R+DR
      if    (if2 .eq. 0) then
        U(MK  )=U(MK  )-VR*EXP(-(R/RY)**2)*R**PWR
      elseif(if2 .eq. 2) then
        V(MK  )=V(MK  )-VR*EXP(-(R/RY)**2)*R**PWR
      endif
 1060 CONTINUE
      endif
      if(VI .ne. 0.0) then
      R=0.0
      KT=10.0*RZ/DR
      K=MIN0(MAX0(K,KT),KMAX)
      DO 1070 M=1,KT
      MK=M+M-1
      R=R+DR
      if    (if2 .eq. 0) then
        U(MK+1)=U(MK+1)-VI*EXP(-(R/RZ)**2)*R**PWR
      elseif(if2 .eq. 2) then
        V(MK+1)=V(MK+1)-VI*EXP(-(R/RZ)**2)*R**PWR
      endif
 1070 CONTINUE
      endif
      GO TO 6000
C
C     DEFORMED VOLUME OR SURFACE BY YLM EXPANSION
C
 1100 CONTINUE
      DO 1101 I=1,3
      BETA(I)=ABETA(I)
      FACTR(I)=SQRT(FLOAT(LDFRM(I)+LDFRM(I)+1))*BETA(I)/(SQRPI*2.0)
 1101 CONTINUE
      LMAX=MAX0(LAM,LDFRM(3))+1
      CALL LEGAUS(16,X,W)
      T2=(-1.0)**LAM
      T1=(-1.0)**LDFRM(1)
      DO 1109 I=1,8
      CN(I  )=0.0
      CN(I+8)=0.0
      P2=0.0
      P3=1.0
      Z=X(I)
      DO 1108 M=1,LMAX
      L=M-1
      FL=L-1
      IF(L.EQ.0) GO TO 1102
      P3=((2.0*FL+1.0)*Z*P2-FL*P1)/(FL+1.0)
 1102 CONTINUE
      DO 1104 J=1,3
      IF(L.NE.LDFRM(J)) GO TO 1104
      CN(I  )=CN(I  )+P3*FACTR(J)
      CN(I+8)=CN(I+8)+P3*FACTR(J)*T1
 1104 CONTINUE
      IF(L.NE.LAM) GO TO 1107
      YLAM(I  )= P3*W(I)*SQRT(FLOAT(L+L+1))*SQRPI
      YLAM(I+8)=YLAM(I  )*T2
 1107 CONTINUE
      P1=P2
      P2=P3
 1108 CONTINUE
 1109 CONTINUE
      DO 1170 I=1,2
      GO TO (1110,1120),I
 1110 CONTINUE
      IF(VR.EQ.0.0) GO TO 1170
      VX=VR
      RX=RR
      AX=AR
      F1=1.0
      IFL=-1
      GO TO 1130
 1120 CONTINUE
      IF(VI.EQ.0.0) GO TO 1170
      VX=VI
      RX=RI
      AX=AI
      F1=1.0
      F2=F4
      IFL=0
 1130 CONTINUE
      DO 1135 J=1,16
      CP(J)=EXP((1.0+CN(J))*RX/AX)
 1135 CONTINUE
      VX=VX/BETA(1)
      IF(LAM.EQ.0) VX=-VX*BETA(1)/(2.0*SQRPI)
      J=NX-10
      GO TO (1141,1151),J
 1141 CONTINUE
      DO 1150 M=1,K
      MK=M+M+IFL
      VTEMP=0.0
      F1=F1*F2
      DO 1145 J=1,16
      F3=F1*CP(J)
      VTEMP=VTEMP+YLAM(J)*VX*F3/(1.0+F3)
 1145 CONTINUE
      U(MK  )=U(MK  )+VTEMP
 1150 CONTINUE
      GO TO 1170
 1151 CONTINUE
      DO 1160 M=1,K
      MK=M+M+IFL
      VTEMP=0.0
      F1=F1*F2
      DO 1155 J=1,16
      F3=F1*CP(J)
      VTEMP=VTEMP-YLAM(J)*VX*F3/(1.0+F3)**2
 1155 CONTINUE
      U(MK  )=U(MK  )+VTEMP
 1160 CONTINUE
 1170 CONTINUE
      GO TO 6000
 1200 CONTINUE
      GO TO 1100
C
C     HULTHEN WAVE FUNCTION= NORM*(EXP(-R/RY)-EXP(-R/RZ))/R
C
c     vr .eq. 0,  Hulthen function
c     vr .ne. 0,  V(r)*Hulthen function
c
 1300 CONTINUE
      READ (input,9000)G
      if(icon19.eq.0) then
      WRITE(ioutp,9508)G
      endif
      T1=1.0/RY
      T2=1.0/RZ
      T3=T2**2-T1**2
      T4=SQRT(2.0*T1*(T1+T2)*T2)/(T2-T1)
      IF(VR.NE.0.0) KT=16.0* MIN (RY,RZ)/DR
      IF(VR.EQ.0.0) KT=16.0* MAX (RY,RZ)/DR
      K=MIN0(MAX0(K,KT),KMAX)
      F1=1.0
      F2=EXP(-DR   *T1)
      F3=1.0
      F4=EXP(-DR   *T2)
      R=0.0
      DO 1360 M=1,K
      MK=M+M-1
      R=R+DR
      F1=F1*F2
      F3=F3*F4
      V(M)=T3*F3/(F1-F3)
      IF(VR.EQ.0.0) THEN
        TEMP=1.0
      ELSE
        TEMP=V(M)/FACT
      ENDIF
      U(MK  )=U(MK  )+TEMP*T4*(F1-F3)/R
 1360 CONTINUE
      GO TO 6000
c
c     Yukawa L=LAM
c
c     v = V-0*exp(-r/rx)/(r/rx)   ay < r 
c       = Wood-Saxon              r  < ay
c
 1400 CONTINUE
      if(ry.eq.0.0) ry=1.0
      if(rz.eq.0.0) rz=1.0
      if(ar.eq.0.0) ar=dr
      if(ai.eq.0.0) ai=dr
      f1=1.0
      f2=exp(-dr/ry)
      xr=ar/ry
      t1=1.0+1.0/xr
      t2=1.0
c     Recur Hankel functions
      do 1420 i=0,lam
      t0=t1
      t1=t2
      t2=float(2*i-1)*t1/xr+t0
 1420 continue
      vzr  =2.0*vr*t2*exp(-xr)/xr
      betar=2.0*(float(lam+1)*t2/xr+t1)/t2
      f3=1.0
      f4=exp(-dr/rz)
      xi=ai/rz
      t1=1.0+1.0/xi
      t2=1.0
c     Recur Hankel functions
      do 1430 i=0,lam
      t0=t1
      t1=t2
      t2=float(2*i-1)*t1/xi+t0
 1430 continue
      vzi  =2.0*vi*t2*exp(-xi)/xi
      betai=2.0*(float(lam+1)*t2/xi+t1)/t2
c
      if(vr.ne.0.0) then
      r =0.0
      do 1460 m=1,k
      r=r+dr
      mk=m+m-1
      f1=f1*f2
      f3=f3*f4
      if(r.lt.ar) then
        u(mk  )=u(mk  )-vzr/(1.0+exp(betar*(r/ry-xr)))
      else
      xr=r/ry
      t1=1.0+1.0/xr
      t2=1.0
c     Recur Hankel functions
      do 1450 i=0,lam
      t0=t1
      t1=t2
      t2=float(2*i-1)*t1/xr+t0
 1450 continue
        u(mk  )=u(mk  )-vr*t2*f1/xr
      endif
 1460 continue
      endif
      if(vi.ne.0.0) then
      r =0.0
      do 1480 m=1,k
      r=r+dr
      mk=m+m-1
      if(r.lt.ai) then
        u(mk+1)=u(mk+1)-vzi/(1.0+exp(betai*(r/rz-xi)))
      else
      xi=r/rz
      t1=1.0+1.0/xi
      t2=1.0
c     Recur Hankel functions
      do 1470 i=0,lam
      t0=t1
      t1=t2
      t2=float(2*i-1)*t1/xi+t0
 1470 continue
        u(mk+1)=u(mk+1)-vi*t2*f3/xi
      endif
 1480 continue
      endif
      go to 6000
c
c     Yukawa L*S
c
 1500 CONTINUE
      r =0.0
      if(ry.eq.0.0) ry=1.0
      if(rz.eq.0.0) rz=1.0
      if(ar.eq.0.0) ar=1.0
      if(ai.eq.0.0) ai=1.0
      f1=1.0
      f2=exp(-dr/ry)
      xr=ar/ry
      betar=2.0*(1.0+3.0/xr+3.0/xr**2)/(1.0+1.0/xr)
      vzr  =2.0*vr*exp(-xr)*(1.0+1.0/xr)/xr**2
      f3=1.0
      f4=exp(-dr/rz)
      xi=ai/rz
      betai=2.0*(1.0+3.0/xi+3.0/xi**2)/(1.0+1.0/xi)
      vzi  =2.0*vi*exp(-xi)*(1.0+1.0/xi)/xi**2
      do 1560 m=1,k
      r=r+dr
      mk=m+m-1
      f1=f1*f2
      f3=f3*f4
      if(r.lt.ar) then
        u(mk  )=u(mk  )-vzr/(1.0+exp(betar*(r/ry-xr)))
      else
        u(mk  )=u(mk  )-vr *f1*(1.0+1.0*ry/r)/(r/ry)**2
      endif
      if(r.lt.ai) then
        u(mk+1)=u(mk+1)-vzi/(1.0+exp(betai*(r/rz-xi)))
      else
        u(mk+1)=u(mk+1)-vi *f3*(1.0+1.0*rz/r)/(r/rz)**2
      endif
 1560 continue
      go to 6000
C
C     VOLUME WOODS SAXON/R
C
 1600 CONTINUE
      r=0.0
      DO 1660 M=1,K
      MK=M+M-1
      r=r+dr
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )-VR*F1/(1.0+F1)/r
      U(MK+1)=U(MK+1)-VI*F3/(1.0+F3)/r
 1660 CONTINUE
      GO TO 6000
C
C     1ST DER WOODS SAXON/R
C
 1700 CONTINUE
      r=0.0
      DO 1760 M=1,K
      MK=M+M-1
      r=r+dr
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )+VR*F1/(1.0+F1)**2/r
      U(MK+1)=U(MK+1)+VI*F3/(1.0+F3)**2/r
 1760 CONTINUE
      GO TO 6000
C
C     2ND DER WOODS SAXON/R
C
 1800 CONTINUE
      r=0.0
      DO 1860 M=1,K
      MK=M+M-1
      r=r+dr
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )-VR*F1*(1.0-F1)/(1.0+F1)**3/r
      U(MK+1)=U(MK+1)-VI*F3*(1.0-F3)/(1.0+F3)**3/r
 1860 CONTINUE
      GO TO 6000
 1900 CONTINUE
      GO TO 6000
C
C     VECTOR POTENTIAL
C     VSOR, VSOI ARE THE FERMI THIRD PARAMETERS
C
 2000 CONTINUE
      IVFLAG=-1
      R=0.0
      DO 2060 M=1,K
      MK=M+M-1
      R=R+DRF
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )-VR*F1*(1.0+VSOR*(R/RR)**2)/(1.0+F1)
      U(MK+1)=U(MK+1)-VI*F3*(1.0+VSOI*(R/RI)**2)/(1.0+F3)
 2060 CONTINUE
      GO TO 6000
C
C     SCALAR POTENTIAL
C     VSOR, VSOI ARE THE FERMI THIRD PARAMETERS
C
 2100 CONTINUE
      ISFLAG=-1
      VR=VR*(SFACT/VFACT)
      VI=VI*(SFACT/VFACT)
      R=0.0
      DO 2160 M=1,K
      MK=M+M-1
      R=R+DRF
      F1=F1*F2
      F3=F3*F4
      V(MK  )=V(MK  )-VR*F1*(1.0+VSOR*(R/RR)**2)/(1.0+F1)
      V(MK+1)=V(MK+1)-VI*F3*(1.0+VSOI*(R/RI)**2)/(1.0+F3)
 2160 CONTINUE
      GO TO 6000
c
c     p_pbar and lambda_lambdabar potentials
c     (S = 1, T = 0)  (S = 0, T = 1)  (S = 1, T = 0)  (S = 0, T = 1)
c
c     drf    = radial increment
c     k      = number of radial points
c     inf    = 0, n_nbar potentials  (nucleon)
c            = 1, h_hbar potentials  (hyperon)
c     ispin  = spin    = 0 or 1
c     tspin  = isospin = 0 or 1
c     kbb    = 0, l_lbar potentials  (lambda)
c            = 1, l_sbar potentials  (lambda-sigma)
c     jcode  = 1, volume real
c            = 2, volume imaginary
c            = 3, Tensor real
c            = 4, Spin orbit real
c            = 5, Spin orbit imaginary (not implemented)
c     vr,vi  = strength factor for potential
c     u      = storage for potential
c
 3100 continue
      jcode = 1
      RM = RY
      MM = 1
      VV = vr
      go to 3590
 3200 continue
      jcode = 2
      RM = RZ
      MM = 2
      VV = vi
      go to 3590
 3300 continue
      jcode = 3
      RM = RY
      MM = 1
      VV = vr
      go to 3590
 3400 continue
      jcode = 4
      RM = RY
      MM = 1
      VV = vr
      go to 3590
 3500 continue
      jcode = 5
      RM = RZ
      MM = 2
      VV = vi
 3590 CONTINUE
      call antipot(drf,k,inf,ispin,tspin,kbb,jcode,VV,u(MM),RM)
      go to 6000
C
C      END OF POTENTIALS
C
 6000 CONTINUE
      IF(KFLAG.NE.0.OR.NX.GT.5) GO TO 6050
      IF(ABS(VSOR)+ABS(VSOI).EQ.0.0) GO TO 6050
      NX=NX+3
      KFLAG=1
      VR=VR*VSOR/45.2
      VI=VI*VSOI/45.2
      GO TO 80
 6050 CONTINUE
      IF(FZ.GT.0.0) GO TO 70
C
C     PROCESS DIRAC POTENTIALS
C
C     ENTRY WITH  U -> K**2 - VFACT*V   VFACT = 2.0*W1 /HBARC**2
C                 V ->      - SFACT*S   SFACT = 2.0*FM1/HBARC**2
C
      IF(IVFLAG.NE.0.AND.ISFLAG.NE.0) THEN
      IDIRAC=-1
C
	KT2=K+K
C	WRITE(20,7777)' ENTRY  POTENTIALS'
C	WRITE(20,7778)(U(I),I=1,KT2)
C	WRITE(20,7777)
C	WRITE(20,7778)(V(I),I=1,KT2)
      W1M1=(VFACT+SFACT)*0.5*HBARC**2
 
      DO 6100 M=1,KT
      MK=M+M-1
      VVR=(U(MK  )-FK2**2)/VFACT
      VVI= U(MK+1)        /VFACT
      VSR= V(MK  )        /SFACT
      VSI= V(MK+1)        /SFACT
      U(MK  )=U(MK  )+V(MK  ) +(VVR**2-VVI**2 - VSR**2+VSI**2)/HBARC**2
      U(MK+1)=U(MK+1)+V(MK+1) +(VVR*VVI - VSR*VSI)*2.0/HBARC**2
      T1 = W1M1 + VVR-VSR
      T2 =        VVI-VSI
      V(MK  )=0.5* LOG (T1**2 + T2**2)
      V(MK+1)=ATAN2(T2,T1)
 6100 CONTINUE
C	WRITE(20,7777)' SECOND POTENTIALS'
C	WRITE(20,7778)(U(I),I=1,KT2)
C	WRITE(20,7777)
C	WRITE(20,7778)(V(I),I=1,KT2)
C
      R=FLOAT(K+1)*DRF
      MK=K+K-1
      D2=V(MK  )
      D1=D2
      A2=V(MK+1)
      A1=A2
      IBX=1
      DO 6150 M=2,K
      MK=MK-2
      R=R-DRF
      D3=D2
      D2=D1
      D1=V(MK  )
      A3=A2
      A2=A1
      A1=V(MK+1)
      DPR=(D3-D1)/(2.0*DRF)
      DPI=(A3-A1)/(2.0*DRF)
      V(MK+2)=2.0*DPR/R
      V(MK+3)=2.0*DPI/R
C
      DPPR=(D3+D1-2.0*D2)/DRF**2
      DPPI=(A3+A1-2.0*A2)/DRF**2
      UDR =0.5*DPPR-0.25*(DPR**2-DPI**2)+DPR/R
      UDI =0.5*DPPI-0.25*(2.0*DPR*DPI  )+DPI/R
C      UD(MK+2)=UDR
C      UD(MK+3)=UDI
      U(MK+2)=U(MK+2)+UDR
      U(MK+3)=U(MK+3)+UDI
 6150 CONTINUE
      V(1   )=V(3   )*2.0
      V(2   )=V(4   )*2.0
      U(1   )=U(1   )+UDR
      U(2   )=U(2   )+UDI
C	WRITE(20,7777)'TERTIARY POTENTIALS'
C	WRITE(20,7777)'CENTRAL POTENTIAL'
C	WRITE(20,7778)(U(I),I=1,KT2)
C	WRITE(20,7777)'SPIN ORBIT POTENTIAL'
C	WRITE(20,7778)(V(I),I=1,KT2)
C	WRITE(20,7777)'UD - DARWIN TERM'
C	WRITE(20,7778)(UD(I),I=1,KT2)
C 7777	FORMAT(A30)
C 7778	FORMAT(' ',1P10E12.4)
      ENDIF
C
C
      IF((sp_code.EQ.1 .OR. sp_code.EQ.2) .and. fk2.ne.0.0) THEN
C
C     SINGLE PARTICLE ORBITAL
C
C
C     READ IN QUANTUM NUMBERS FOR SINGLE PARTICLE ORBITAL
C
      READ (input,9000)G,VTRIAL,FISW,DAMP
      ISW=FISW
      FJ0=FJ2/2.0
      IF(VTRIAL.EQ.0.0) VTRIAL=ETA6
      FS0=FSS/2.0
      if(icon19.eq.0) then
      WRITE(ioutp,'(a)')
      IF(DAMP.NE.0.0) WRITE(ioutp,9515)DAMP
      WRITE(ioutp,9508)G,FISW
      WRITE(ioutp,'(a)')
      endif
      TEMP=(FJ0**2+FJ0-FL**2-FL-FS0**2-FS0)*0.5
      DO 8128 M=1,KMAX
      MK=M+M-1
      V(MK  )=U(MK  )+V(MK  )*TEMP
      V(MK+1)=U(MK+1)+V(MK+1)*TEMP
 8128 CONTINUE
      JERR=-icon19
c
      CALL BIND(V,U,DR,RM,FN,FL,K,FK,ETA,VTRIAL,ECM,FK2,ISW,JERR,D)
c
      IF(JERR.NE.0) IERR=JERR
      DO 8150 M=1,KMAX
      MK=M+M-1
      V(M)=VTRIAL*V(MK  )+V(MK+1)
 8150 CONTINUE
      Anorm=1.0
      FACT=PNLOC**2/8.0
      IF(FACT.ne.0.0) then
C
C     NON-LOCAL CORRECTION FOR SINGLE PARTICLE FUNCTION
C
      SUM=0.0
      R=0.0
      DO 8175 M=1,K
      MK=M+M-1
      R=R+DR
      U(M)=U(M)*EXP(FACT*(FK2-V(M)))
      SUM=SUM+(U(M)*R)**2
 8175 CONTINUE
 
C
C     DO NOT RENORMALIZE FOR POSITIVE ENERGY
C
      IF(FK2   .LT.0.0) then
        Anorm=1.0/SQRT(SUM*DR)
      ELSE
        Anorm=1.0
      ENDIF
      ENDIF
      IF(DAMP.EQ.0.0) GO TO 8180
C
C     APPLY DAMPING FACTOR, EXP(-DAMP*R) TO FORM FACTOR
C
      R=0.0
      SUM=0.0
      F1=1.0
      F2=EXP(-DAMP*DR   )
      DO 8177 M=1,K
      F1=F1*F2
      R=R+DR
      U(M)=U(M)*F1
      SUM=SUM+(U(M)*R)**2
 8177 CONTINUE
      Anorm=1.0/SQRT(SUM*DR   )
 8180 CONTINUE
      DO 8200 M=KMAX,1,-1
      MK=M+M-1
      IF(M.GT.K) then
        V(M   )=0.0
        U(MK  )=0.0
      ELSE
        U(MK  )=U(M)*Anorm
      ENDIF
      U(MK+1)=0.0
 8200 CONTINUE
      ENDIF
C
      BIF=RM/DR
c     Store quantum numbers for transfer back
      DO 8245 M=28,32
      AKIN(M)=BKIN(M)
 8245 CONTINUE
      RETURN
C
 9000 FORMAT(10F8.4)
 9008 FORMAT(22H0POTENTIAL PARAMETERS )
 9010 FORMAT( 9H0PARTICLE,I2,120(1H*))
 9100 FORMAT(5E16.7)
 9500 FORMAT(1H ,3A6,5(3X,A6,F9.4))
 9508 FORMAT(18X,9H   NODES=,F9.4,9H   L    =,F9.4 ,9H   2*J  =,F9.4
     1          ,9H   2*S  =,F9.4,9H   FISW =,F9.4)
 9509 FORMAT(A18,A3,6HV RL =,F9.4,9H   R0RL =,F9.4 ,9H   A RL =,F9.4
     1          ,9H   R RL =,F9.4,9H   VSOR =,F9.4:,9H   Powr =,F9.4)
 9510 FORMAT(A18,A3,6HV IM =,F9.4,9H   R0IM =,F9.4 ,9H   A IM =,F9.4
     1          ,9H   R IM =,F9.4,9H   VSOI =,F9.4:,9H   Powr =,F9.4)
 9511 FORMAT(A18,A3,6HV RL =,F9.3,9H   R0RL =,F9.4 ,9H   A RL =,F9.4
     1          ,9H   R RL =,F9.4,9H   VSOR =,F9.4:,9H   Powr =,F9.4)
 9512 FORMAT(A18,A3,6HV IM =,F9.3,9H   R0IM =,F9.4 ,9H   A IM =,F9.4
     1          ,9H   R IM =,F9.4,9H   VSOI =,F9.4:,9H   Powr =,F9.4)
 9515 FORMAT(18X,9H   DAMP =,F9.4)
      END
 

