c***********************************************************************
      SUBROUTINE POTS(U,V)
c
c     Calculates the potentials or form factors
c***********************************************************************
c
c     IMPLICIT REAL*8(A-H,O-Z)
      COMMON/POTTER/DRX,AFACT(2),VFACT,SFACT,ENG,RM,G(4),ETAX,FKX,ETAKX
     1             ,RCX,HBARC2,ABETA(3),FLDF(3)
     2             ,NX,LAMX,KMXX,KX,IBX,LPLUSX,ICON4,NSPC,IDIRAC,ICHK
 
      CHARACTER*18 B(26)
      DIMENSION UT(5),CN(16),CP(16),YLAM(16),PLM(8)
     1         ,XG(8),WG(8),U(800),V(800),LDFRM(3)
C      DIMENSION UD(800)
      EQUIVALENCE (YLAM(1),PLM(1))
C
c
      DATA b /                ' NX=0 No potential', ' NX=1 VOLUME W-S  '
     1, ' NX=2 SURFAC, 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 Not used   '
     8, ' NX=23 Not used   ', ' NX=24 Not used   ', ' NX=25 Not used   '
     9/
C
      DATA NG,NGX/8,0/
C
      ETA4 = 6.0
      ETA5 =10.0
      SQRPI= 1.772453851
      PI=3.141592654
C
      IVFLAG=0
      ISFLAG=0
      FACT=VFACT
   70 CONTINUE
C
C     READ IN CARD SET 5,6,OR 7   POTENTIAL CARDS
C
      READ (5,9000)FZ,VR,RY,AR,VSOR,VI,RZ,AI,VSOI,PWR
C
      NZ=ABS(FZ)
      RR=ABS(RY)*AFACT(1)
      IF(RY.LT.0.0) RR=RR+ABS(RY)*AFACT(2)
      RI=ABS(RZ)*AFACT(1)
      IF(RZ.LT.0.0) RI=RI+ABS(RZ)*AFACT(2)
      IF(ICON4.NE.2.OR.NSPC.LT.3) THEN
        WRITE(6,9509)B(NZ+1),VR,RY,AR,RR,VSOR
        WRITE(6,9510)        VI,RZ,AI,RI,VSOI,PWR
      ENDIF
C
      KFLAG=0
      IF(NX.LT.3) THEN
        VR=VR*FACT
        VI=VI*FACT
        KT=FKX  * MAX (RR,RI)+ETA5
        LPLUSX=MAX0(LPLUSX,KT)
        KT=(2.3*ETA4* MAX (AR,AI)+ MAX (RR,RI))/DRX
      ELSE
        IF(ENG.EQ.0.0) THEN
          KT=(2.3*ETA4* MAX (AR,AI)+ MAX (RR,RI))/DRX
        ELSE
          RM= MAX (RM,RR)
          RM= MAX (RM,RI)
          IF(RM.EQ.0.0) RM=1.0
          VR=VR*FACT
          VI=VI*FACT
          KT=(2.3*ETA4/SQRT(FKX**2+2.0*ETAX*FKX/RM))/DRX
        ENDIF
      ENDIF
      KX=MIN0(MAX0(KX,KT),KMXX)
   83 CONTINUE
      IF(AR.EQ.0.0) THEN
        F1=0.0
        F2=0.0
      ELSE
        F2=EXP(-DRX  /AR)
        F1=EXP( RR/AR)
      ENDIF
      IF(AI.EQ.0.0) THEN
        F3=0.0
        F4=0.0
      ELSE
        F4=EXP(-DRX  /AI)
        F3=EXP( RI/AI)
      ENDIF
C
      IF(NX.GE.3.AND.ENG.EQ.0.0.AND.NZ.LE.5) THEN
        IF(AR.NE.0.0) VR=VR*(RR/AR)**(PWR+1)
        IF(AI.NE.0.0) VI=VI*(RI/AI)**(PWR+1)
      ENDIF
      IF(NZ.EQ.0) GO TO 6050
      GO TO ( 100, 200, 300, 400, 500, 600, 700, 800, 900,1000),NZ
      GO TO (1100,1200,1300,1400,1500,1600,1700,1800,1900,2000),NZ-10
      GO TO (2100),NZ-20
      write(6,'(''0Invalid potential option '',i3)')nz
      
C
C     VOLUME WOODS SAXON
C
  100 CONTINUE
      DO 160 M=1,KX
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      U(MK  )=U(MK  )-VR*F1/(1.0+F1)
      U(MK+1)=U(MK+1)-VI*F3/(1.0+F3)
  160 CONTINUE
      GO TO 6000
C
C     1ST DERIVATIVE WOODS SAXON
C
  200 CONTINUE
      DO 260 M=1,KX
      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 DERIVATIVE WOODS SAXON
C
  300 CONTINUE
      DO 360 M=1,KX
      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
      IBX=1
      IF(AR.NE.0.0) VR=VR/AR
      IF(AI.NE.0.0) VI=VI/AI
      R=0.0
      DO 460 M=1,KX
      R=R+DRX
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      V(MK  )=V(MK  )-VR*F1/(R*(1.0+F1)**2)
      V(MK+1)=V(MK+1)-VI*F3/(R*(1.0+F3)**2)
  460 CONTINUE
      GO TO 6000
C
C     L.S 1ST DERIVATIVE WOODS SAXON
C
  500 CONTINUE
      IBX=1
      IF(AR.NE.0.0) VR=VR/AR
      IF(AI.NE.0.0) VI=VI/AI
      R=0.0
      DO 560 M=1,KX
      R=R+DRX
      MK=M+M-1
      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,KX
      R=R+DRX
      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 DERIVATIVE WOOD SAXON*R**PWR
C
  700 CONTINUE
      R=0.0
      DO 760 M=1,KX
      MK=M+M-1
      R=R+DRX
      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
      IF(NX.GE.3) THEN
        READ (5,9000)G
        WRITE(6,9508)G
      ENDIF
      READ(5,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
      IF(F2.EQ.0.0.OR.F2.EQ.2) THEN
        F3=VR
        MK=1
      ELSE
        F3=VI
        MK=2
      ENDIF
      IF(F3.EQ.0.0) F3=1.0
      KT=F1
      DO 820 M=1,KT,5
      READ(5,9100)UT
      IF(M.GT.KMXX) GO TO 820
      DO 810 I=1,5
      IF(F2.LT.2.0) THEN
        U(MK  )=U(MK  )+UT(I)*F3
      ELSE
        V(MK  )=V(MK  )+UT(I)*F3
      ENDIF
      MK=MK+2
  810 CONTINUE
  820 CONTINUE
C
      KX=min0(KT,kmxx)
 
      GO TO 6000
C
C     HARMONIC OSCILLATOR, NORMALIZED   N*EXP(-(R/RY)**2/2)
C
  900 CONTINUE
      READ (5,9000)G
      WRITE(6,9508)G
      F1=1.0/RY**2
      F2=F1/RY
      F3=0.5
      F4=SQRPI*0.5
      L=G(2)
      IF(L.NE.0) THEN
        DO 930 I=1,L
        F3=F3+1.0
        F4=F4*F3
        F2=F2*F1
  930   CONTINUE
      ENDIF
      NN=G(1)
      T1=1.0
      T2=F4
C     LAGUERRE POLYNOMIAL COEFFICIENTS = (Abramowitz-Stegun)*(-1)**nn
      CN(1)=(-1.0)**NN
      IF(NN.NE.0) THEN
        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
      ENDIF
      ANORM=SQRT(2.0*F2*T1/T2)*T2/(T1*F4)
      IF(VR.NE.0.0) ANORM=VR*ANORM
      KT=10.0*RY/DRX
      KT=MIN0(KT,KMXX)
      R=0.0
      F1=F1/2.0
      DO 960 M=1,KT
      MK=M+M-1
      R=R+DRX
      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*ANORM*EXP(-F1*R2)*R**L
  960 CONTINUE
      GO TO 6000
C
C     EXP(-(R/R0)**2)*R**POWR
C
 1000 CONTINUE
      IF(VR.NE.0.0) THEN
        R=0.0
        DO 1060 M=1,KX
        MK=M+M-1
        R=R+DRX
        U(MK  )=U(MK  )-VR*EXP(-(R/RY)**2)*R**PWR
 1060 CONTINUE
      endif
      IF(VI.NE.0.0) THEN
        R=0.0
        DO 1070 M=1,KX
        MK=M+M-1
        R=R+DRX
        U(MK+1)=U(MK+1)-VI*EXP(-(R/RZ)**2)*R**PWR
 1070 CONTINUE
      endif
      GO TO 6000
C
C     DEFORMED VOLUME OR SURFACE BY YLM EXPANSION
C
 1100 CONTINUE
      IF(NGX.NE.NG) THEN
        CALL LEGAUS(2*NG,XG,WG)
        NGX=NG
      ENDIF
      T2=(-1.0)**LAMX
      IF(ICHK.EQ.0) THEN
        READ (5,9000) (ABETA(J),FLDF(J) ,J=1,3)
        WRITE(6,9000)
        WRITE(6,9512) (ABETA(J),FLDF(J) ,J=1,3)
      ENDIF
      LMAX=LAMX+1
      DO 1101 J=1,3
      LDFRM(J)=FLDF(J)
      LMAX=MAX0(LMAX,LDFRM(J)+1)
 1101 CONTINUE
      T2=(-1.0)**LAMX
C
      DO 1108 I=1,NG
      CN(I  )=0.0
      CN(I+NG)=0.0
      P2=0.0
      P3=1.0
C
      DO 1106 M=1,LMAX
      L=M-1
      FL=L-1
      IF(L.EQ.0) GO TO 1102
      P3=((2.0*FL+1.0)*XG(I)*P2-FL*P1)/(FL+1.0)
 1102 CONTINUE
      DO 1103 J=1,3
      IF(ABETA(J).EQ.0.0) GO TO 1103
      IF(L.NE.LDFRM(J)) GO TO 1103
      FACTOR=P3*ABETA(J)*SQRT(FLOAT(L+L+1)/(4.0*PI))
      CN(I  )=CN(I  )+FACTOR
      CN(I+NG)=CN(I+8)+FACTOR*(-1.0)**LDFRM(J)
 1103 CONTINUE
      IF(L.NE.LAMX) GO TO 1104
      YLAM(I  )= P3*WG(I)*SQRT(FLOAT(L+L+1)*PI)
      YLAM(I+NG)=YLAM(I  )*T2
 1104 CONTINUE
      P1=P2
      P2=P3
 1106 CONTINUE
 1108 CONTINUE
      DO 1170 I=1,2
      IF(I.EQ.1) THEN
        IF(VR.EQ.0.0) GO TO 1170
        VX=VR
        RX=RR
        AX=AR
        F1=1.0
        IFL=-1
      ELSE
        IF(VI.EQ.0.0) GO TO 1170
        VX=VI
        RX=RI
        AX=AI
        F1=1.0
        F2=F4
        IFL=0
      ENDIF
      DO 1135 J=1,16
      CP(J)=EXP((1.0+CN(J))*RX/AX)
 1135 CONTINUE
      IF(LAMX.EQ.0) VX=VX/(SQRPI*2.0)
      J=NZ-10
      IF(J.EQ.2) GO TO 1151
      DO 1150 M=1,KX
      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,KX
      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
 1300 CONTINUE
      READ (5,9000)G
      WRITE(6,9508)G
      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)/DRX
      IF(VR.EQ.0.0) KT=16.0* MAX (RY,RZ)/DRX
      KX=MIN0(KX,KMXX)
      F1=1.0
      F2=EXP(-DRX  *T1)
      F3=1.0
      F4=EXP(-DRX  *T2)
      R=0.0
      DO 1260 M=1,KX
      MK=M+M-1
      R=R+DRX
      F1=F1*F2
      F3=F3*F4
      V(M)=T3*F3/(F1-F3)
      TEMP=1.0
      IF(VR.NE.0.0) TEMP=V(M)/FACT
      U(MK  )=U(MK  )+TEMP*T4*(F1-F3)/R
 1260 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
	lam=lamx
      if(ry.eq.0.0) ry=1.0
      if(rz.eq.0.0) rz=1.0
      if(ar.eq.0.0) ar=drx
      if(ai.eq.0.0) ai=drx
      f1=1.0
      f2=exp(-drx/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(-drx/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+drx
      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+drx
      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(-drx/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(-drx/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+drx
      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
 1600 CONTINUE
 1700 CONTINUE
       go to 6000
 1800 CONTINUE
      DO 1860 M=1,KX
      MK=M+M-1
      F1=F1*F2
      F3=F3*F4
      V(MK  )=V(MK  )-VR*F1/(1.0+F1)
      V(MK+1)=V(MK+1)-VI*F3/(1.0+F3)
 1860 CONTINUE
      GO TO 6000
 1900 CONTINUE
c
c     Coulomb excitation for a deformed uniform charge distribution
c
      IF(ICHK.EQ.0) THEN
        READ (5,9000) (ABETA(J),FLDF(J) ,J=1,3)
        WRITE(6,9000)
        WRITE(6,9512) (ABETA(J),FLDF(J) ,J=1,3)
      ENDIF
      if(vr.eq.0.0) vr=1.0
      do 1990 i=1,3
      beta=abeta(i)
      if(beta.ne.0.0) then
c     set flag for unbound stripping evaluation of coulomb ex.
      ibx=4
      flam1=beta*sqrt(float(2*lamx+1)/(4.0*pi))
      flam2=flam1*phasef(lamx)
c
      if(ngx.ne.ng) then
        ngx=ng
        call legaus(2*ng,xg,wg)
      endif
c
      an=0.0
      do 1920 k=1,ng
      p1=0.0
      p2=1.0
      if(lamx.gt.0) then
        do 1910 j=1,lamx
        p3=(float(2*j-1)*xg(k)*p2 - float(j-1)*p1)/float(j)
        p1=p2
        p2=p3
 1910 continue
      endif
      plm(k)=p2
c
c     calculate normalization
c
      r1=rcx*(1.0+flam1*p2)
      r2=rcx*(1.0+flam2*p2)
      an=an+(r1**3+r2**3)*wg(k)
 1920 continue
c
      an =an/3.0
c
      r  =0.0
c
      do 1940 m=1,kx
      mk=m+m-1
      r=r+drx
      sum=0.0
      do 1930 k=1,ng
      r1=rcx*(1.0+flam1*plm(k))
      r2=rcx*(1.0+flam2*plm(k))
c
      if(r.gt.r1) then
        s1=r1**(lamx+3)/(r**(lamx+1)*float(lamx+3))
      else
        if(lamx.eq.2) then
          s1=( log(r1/r)+1.0/float(lamx+3))*r**2
        else
          s1=(r**lamx/r1**(lamx-2)
     1      -float(2*lamx+1)*r**2/float(lamx+3))/float(2-lamx)
        endif
      endif
c
      if(r.gt.r2) then
        s2=r2**(lamx+3)/(r**(lamx+1)*float(lamx+3))
      else
        if(lamx.eq.2) then
          s2=( log(r2/r)+1.0/float(lamx+3))*r**2
        else
          s2=(r**lamx/r2**(lamx-2)
     1      -float(2*lamx+1)*r**2/float(lamx+3))/float(2-lamx)
        endif
      endif
c
      sum=sum+(s1+s2*phasef(lamx))*wg(k)*plm(k)
 1930 continue
c
      sum=sum/an
      if(beta.ne.0.0.and.lamx.ne.0) sum=sum/flam1
      u(mk  )=u(mk  )+sum*(vr*etakx/fact)
 1940 continue
      endif
 1990 continue
      GO TO 6000
C
C     VECTOR POTENTIAL
C     VSOR, VSOI ARE THE THIRD PARAMETERS
C     IN THE 3 PARAMETER FERMI MODEL
C         [1 + VSO?*(R/R?)**2]
C
 2000 CONTINUE
      IVFLAG=-1
      R=0.0
      DO 2060 M=1,KX
      MK=M+M-1
      R=R+DRX
      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 6020
C
C     SCALAR POTENTIAL
C     VSOR, VSOI ARE THE FERMI THIRD PARAMETERS
C     IN THE 3 PARAMETER FERMI MODEL
C         [1 + VSO?*(R/R?)**2]
C
 2100 CONTINUE
      ISFLAG=-1
      VR=VR*(SFACT/VFACT)
      VI=VI*(SFACT/VFACT)
      R=0.0
      DO 2160 M=1,KX
      MK=M+M-1
      R=R+DRX
      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 6020
C
C      END OF POTENTIALS
C
 6000 CONTINUE
      IDIRAC=1
 6020 CONTINUE
      IF(KFLAG.NE.0.OR.NZ.GT.5) GO TO 6050
      IF(ABS(VSOR)+ABS(VSOI).EQ.0.0) GO TO 6050
      NZ=NZ+3
      KFLAG=1
      VR=VR*VSOR/45.2
      VI=VI*VSOI/45.2
      GO TO 83
 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
      IF(IDIRAC.EQ.1) THEN
        WRITE(6,9515)
      ENDIF
      IDIRAC=-1
C
C      KT2=KX+KX
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*HBARC2
 
      DO 6100 M=1,Kx
      MK=M+M-1
      VVR=(U(MK  )-FKX**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)/HBARC2
      U(MK+1)=U(MK+1)+V(MK+1) +(VVR*VVI - VSR*VSI)*2.0/HBARC2
      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(KX+1)*DRX
      MK=KX+KX-1
      D2=V(MK  )
      D1=D2
      A2=V(MK+1)
      A1=A2
      IBX=1
      DO 6150 M=2,KX
      MK=MK-2
      R=R-DRX
      D3=D2
      D2=D1
      D1=V(MK  )
      A3=A2
      A2=A1
      A1=V(MK+1)
C     FIRST  DERIVATIVE TERMS
      DPR=(D3-D1)/(2.0*DRX)
      DPI=(A3-A1)/(2.0*DRX)
      V(MK+2)=2.0*DPR/R
      V(MK+3)=2.0*DPI/R
C     SECOND DERIVATIVE TERMS
      DPPR=(D3+D1-2.0*D2)/DRX**2
      DPPI=(A3+A1-2.0*A2)/DRX**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
      IF(IDIRAC.EQ.1) IDIRAC=0
      RETURN
C
 9000 FORMAT(10F8.4)
 9100 FORMAT(5E16.7)
 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,9H   V RL =,F9.4,9H   R0RL =,F9.4,9H   A RL =,F9.4
     1          ,9H   R RL =,F9.4,9H   VSOR =,F9.4)
 9510 FORMAT(18X,9H   V 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)
 9512 FORMAT(18X,9H   BETA1=,F9.4,9H   LDFR1=,F9.4,9H   BETA2=,F9.4
     1          ,9H   LDFR2=,F9.4,9H   BETA3=,F9.4,9H   LDFR3=,F9.4)
 9515 FORMAT('0',20('*'),'  WARNING, Mixing of Dirac and non-Dirac'
     1,          ' potentials may be hazardous to your calculation'
     2,          20('*'))
      END