C     TO CONSTRUCT THE K-B TYPE PSEUDOPOTENTIAL
C
C     NEW VERSION (1990 1/8)  THIS PROGRAM MAKES ONLY S AND P.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER(NESH=421)
C
C     ATTN!  GROUND-STATE (S,P) , EXITED-STATE (D)
C            GROUND-STATE (S)   , EXITED-STATE (P,D)
C
      DIMENSION RAD(NESH),VS(NESH),VP(NESH),VD(NESH)
      DIMENSION WFS(NESH),WFP(NESH),WFD(NESH)
      DIMENSION VSKB(NESH),VPKB(NESH),VDKB(NESH),VCOR(NESH)
      DIMENSION OMO(NESH),VDNL(NESH),VDUM(NESH)
C
C     READING THE DATA(POTENTIAL AND WAVE-FUNCTION)
C
      OPEN(UNIT=10,FILE='./fort.92')
      OPEN(UNIT=11,FILE='./fort.91')
      OPEN(UNIT=12,FILE='./PDKB.DAT')
      OPEN(UNIT=14,FILE='./PDD.DAT')
C
      REWIND 10
      REWIND 11
      REWIND 14
      READ (10,*) MESH,NMES,DX,RAD,VP,VD,WFP,WFD
      READ (11,*) MESH,NMES,DX,RAD,VS,WFS,VDUM
      WRITE (6,*) 'INPUT AA1,AA2,ACC!  3.31,1.32,2.5256'
      READ (5,*) AA1
      READ (5,*) AA2
      READ (5,*) ACC
C            PA
C      AA1  = 3.31D0
C      AA2  = 1.32D0
C      ACC  = 2.5256D0
C
      ALP1 = AA1
      ALP2 = AA2
      ACOEF= ACC
      ZVAL = 10.0D0
C
      DO 3000 N=1,MESH
      X1   = DSQRT(ALP1)*RAD(N)
      X2   = DSQRT(ALP2)*RAD(N)
      ERF1 = DERF(X1)
      ERF2 = DERF(X2)
C     WRITE (6,*) ERF1,ERF2
      IF (RAD(N).EQ.0.0D0) THEN
          VCOR(N)=-ZVAL*(ACOEF*(DSQRT(ALP1)-DSQRT(ALP2))+DSQRT(ALP2))
      ELSE
          VCOR(N)=-ZVAL*(ACOEF*(ERF1-ERF2)+ERF2)/RAD(N)
      END IF
 3000 CONTINUE
C
C     K-B TYPE SEPARATION FORM
C
      DO 1000 I=1,MESH
C      VSKB(I)=(VS(I)-VCOR(I))*WFS(I)
C      VPKB(I)=(VP(I)-VCOR(I))*WFP(I)
C      VDKB(I)=(VD(I)-VCOR(I))*WFD(I)
      VSKB(I)=(VS(I)-VCOR(I))*WFS(I)
      VPKB(I)=(VP(I)-VCOR(I))*WFP(I)
      VDKB(I)=(VD(I)-VCOR(I))*WFD(I)
      VDNL(I)= VD(I)-VCOR(I)
 1000 CONTINUE
C
C     INTEGRATION OF WEIGHT (WS,WP)
C
      OMO(1)     = 1.0D0
      DO 1101 N  = 2,MESH-1,2
      OMO(N)     = 4.0D0
      OMO(N+1)   = 2.0D0
 1101 CONTINUE
      OMO(MESH)  = 1.0D0
C     THIS ROUTINE ASSUMES LOG-MESH INTEGRATION.
      WS=0.0D0
      WP=0.0D0
      WD=0.0D0
      W1=0.0D0
      W2=0.0D0
      W3=0.0D0
      DO 1100 I=1,MESH
      WS=WS+OMO(I)*DX*WFS(I)*VSKB(I)*RAD(I)
      WP=WP+OMO(I)*DX*WFP(I)*VPKB(I)*RAD(I)
      WD=WD+OMO(I)*DX*WFD(I)*VDKB(I)*RAD(I)
      W1=W1+OMO(I)*DX*VSKB(I)**2*RAD(I)
      W2=W2+OMO(I)*DX*VPKB(I)**2*RAD(I)
      W3=W3+OMO(I)*DX*VDKB(I)**2*RAD(I)
 1100 CONTINUE
      WA=WS/SQRT(W1)
      WB=WP/SQRT(W2)
      WC=WD/SQRT(W3)      
      EA=W1/WS*27.2
      EB=W2/WP*27.2
      EC=W3/WD*27.2
      WRITE (6,*) 'WS,P,D = ',WS,WP,WD
      WRITE (6,*) 'WA,B,C = ',WA,WB,WC
      WRITE (6,*) 'W1,2,3 = ',EA/27.2,EB/27.2,EC/27.2
      WRITE (6,*) 'EA,B,C = ',EA,EB,EC
C
      DO 1200 I=1,MESH
      WRITE (20,888) RAD(I),VS(I)
      WRITE (21,888) RAD(I),VP(I)
      WRITE (22,888) RAD(I),VD(I)
      WRITE (23,888) RAD(I),VCOR(I)
 888  FORMAT(1H ,2F12.6)
 1200 CONTINUE
C
      REWIND 12
      REWIND 14
      WRITE (12,*) MESH,NMES,DX,WS,WP,WD,RAD,VSKB,VPKB,VDKB
C      WRITE (14,*) MESH,NMES,DX,RAD,VD,VDNL
C
      CLOSE(10)
      CLOSE(11)
      CLOSE(12)
      CLOSE(14)
C
      STOP
      END