This program outputs the final data of NCPS95.

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.10')
      OPEN(UNIT=11,FILE='fort.11')
      OPEN(UNIT=12,FILE='SISICKBHS.DAT')
      OPEN(UNIT=14,FILE='SISICD.DAT') 
C
      REWIND 10
      REWIND 11
      REWIND 14
      READ (10,*) MESH,NMES,DX,RAD,VS,VP,WFS,WFP
      READ (11,*) MESH,NMES,DX,RAD,VD,WFD,VDUM 
      WRITE (6,*) 'INPUT AA1,AA2,ACC! 2.16,0.86,1.6054'
      READ (5,*) AA1 
      READ (5,*) AA2
      READ (5,*) ACC
C        SI
C      AA1  = 2.16D0
C      AA2  = 0.86D0
C      ACC  = 1.6054D0
C
      ALP1 = AA1
      ALP2 = AA2
      ACOEF= ACC
      ZVAL = 4.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
      VSKB(I)=(VS(I)-VCOR(I))*WFS(I)
      VPKB(I)=(VP(I)-VCOR(I))*WFP(I)
      VDKB(I)=(VD(I)-VCOR(I))*WFD(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
      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