C=== TO PRINT WAVE FUNCTION WITH INTERPOLATION==
C--- print spi  20/01/07
      SUBROUTINE  INTPOL(R5,R8,RO1,BET,
     *ALFA,R7,H,R2,IS,NE5,KT,IN,KSU1)
      REAL*8ALFA,BET,H,RO1,R2(*),
     * R5(*),R7(*),R8(NE5,*),ROT,
     * AA1,T2,AKSI,UU,VV
      INTEGER IN(*),IL,IQ,DIR
      NE=NE5-15
      DO 1 JJ=1,IS
      IAA=JJ
      IF(KSU1.EQ.3.AND.JJ.NE.IS)  GOTO 1
      CALL COLUMN(NE5,R8(1,IAA),R5)
      IF(KSU1.NE.3) WRITE(6,203)
203   FORMAT(60('='))
      WRITE(6,   2) R5(1),R5(NE+4)
    2 FORMAT(  ' ','ENERGY =',25X,E15.8,
     */' ','NORM   =' ,25X,E15.8)
      IF(R5(1).LT.0.0)  WRITE(6, 201) IN(JJ)
201   FORMAT(' PRINCIPAL QUANTUM NUMBER =',I6)
      IL=IDINT(R5(2))
      WRITE(6, 202)IL
202   FORMAT(' ORBITAL MOMENTUM         =',I6)
      IQ=IDINT(R5(NE+6))
      WRITE(6, 205)IQ
205   FORMAT(' NUMBER OF ELECTRONS      =',I6)
      DIR=IDINT(R5(NE+8))
      WRITE(6, 206)DIR
206   FORMAT(' SPIN                     =',I6)
      WRITE(6, 5)
    5 FORMAT( /  '   WAVE FUNCTION:  '
     *        /3('  R',7X,'P(R)', 9X))
      DO 3 IWW=1,KT
      ROT=ALFA*R7(IWW)+BET* DLOG(R7(IWW))
      AA1=(ROT-RO1)/H
      IT1=AA1
      T2=AA1-IT1
      IF(IT1.EQ.0) GOTO 3
      IF(IT1.GT.NE) GOTO 1
      AKSI=1.-T2
      UU=T2*(1.-T2**2)/6.
      VV=AKSI*(1.-AKSI**2)/6.
      R2(IWW)=DSQRT(R7(IWW)/(ALFA*
     *R7(IWW)+BET))*(R5(IT1+2)*AKSI
     *+R5(IT1+3)*T2-VV*(R5(IT1+3)-2.
     * *R5(IT1+2)+R5(IT1+1))-
     *UU*(R5(IT1
     *+4)-2.*R5(IT1+3)+R5(IT1+2)))
    3 CONTINUE
      WRITE(6,4)   (R7(IW),R2(IW),IW=1,KT )
    4 FORMAT(3(   F7.3,E15.6,2X))
    1 CONTINUE
      RETURN
      END
