C    PROGRAM FOR CALCULATION OF ATOMIC WAVE FUNCTIONS
C    IN HARTREE-FOCK-DIRAC APPROXIMATION
C    AUTHORS: L.V.CHERNYSHEVA,V.L.YACHONTOV,M.YU. KUCHIEV
C    FORTRAN VERSION: L.V.CHERNYSHEVA
C=======================================================
      SUBROUTINE RSCFHF(KT,R7,Z,IS,R,NE,BET,H,IN,IL,IQ,JJJ,XI,L1,SS1,
     *   NU,IG,KK,ISIGM1,ISIGM2,GAM,EPS,MU,KSU1,DZ,ZK,ADD,
     *   IP,UD,DIR,
     *   TETA,ALM3,C3,ALA,ALA1,ALA2,ALM6DU,
     *   R1,R2,R3,R4,R5,R6,R8,
     *   R9,R10,R11,R12,R13,R14,R15)

      REAL*8ALFA,BET,H,Z,RO1,EPS,ETOT,ADD,AL0,
     *   GAM(*),TETA(IS),DZ,MU,CO,R,XX,
     *   ALM3(IS),ALM6DU(IS),JJJ(IS),
     *   R1( NE+15),R2(NE+15 ),R3( NE+15),
     *   R4( NE+15),R5( NE+15),R6( NE+15),
     *   R9(NE+15),R10(NE+15),R11(NE+15),R12(NE+15),
     *   R13(NE+15),R14(NE+15),R15(NE+15),
     *   R7( KT),R8(NE+15,3*(IS+1)),
     *   ALA(IS),C3(IS),ALA1(IS),ALA2(IS),DAT(15)

      INTEGER CH,KK(*),ISIGM1(*),
     *  ISIGM2(*),IN(IS),IL(IS),IQ(IS),XI(IS),DIR(IS),UD(IS),
     *  SS1,S1(1),JJ,A,IP(IS),L1,ZK
C==============================================
      NE5=NE+15
      AL0=1/137.036

C     CO=0.7
      CO =1
      NEN=NE*CO
      NE5N=NEN+15
      WRITE(6,1)
      WRITE(*,1)
 1    FORMAT(' GROUND STATE ',
     *'ATOM (ION) WAVE FUNCTIONS'/
     *' IN THE HARTREE-FOCK-DIRAC APPROXIMATION')
      WRITE(6, 4)IS
  4   FORMAT('   NO OF N-L SHELLS         ',I5)
      DO 205 I=1,IS
      XI(I)=(-1)**(JJJ(I)+IL(I)+0.5)*(JJJ(I)+0.5)

C      IF(IL(I).GT.0) JJJ(I)=-0.5  +DFLOAT(IL(I))
C      IF(IL(I).LT.0) JJJ(I)=-0.5  -DFLOAT(IL(I))
C     WRITE(6, 206 ,JJJ(I)
C206  FORMAT('  JJJ           ',E15.8)
 205  CONTINUE
      DO 2 I=1,IS
 2    UD(I)=1
      IF(L1.EQ.111) GOTO 111
      DO 5 I=1,IS
5     IP(I)=1
      CH=1
      I1=1
      S1(1)=SS1
      DO 11 I=1,IS
C =----  IL->XI---------------
      A=2*XI(I)+1
      IF(S1(1).LE.99) A=A*2
      DO 11 J=I1,IS
      JJ=2*JJJ(J)+1
      IF(S1(1) .LE.99)  JJ=JJ*2
      IF(I.LT.J ) CH=CH+IFIX(SNGL((JJJ(I)+JJJ(J)-
     *DABS(JJJ(I)-JJJ(J)))/2+1))
C     WRITE(6, 200)CH
C200   FORMAT(' CHG ',I3)
      IF(I.EQ.J.AND.IQ(I).NE.1) THEN
      CH=CH+IFIX(SNGL(JJJ(I)))
       END IF
      IF(KSU1.NE.1.AND.(IQ(I).EQ.1.AND.
     * IQ(J).EQ.JJ-1).OR.
     * (IQ(J).EQ.1.AND.IQ(I).EQ.A-1))
     * CH=CH+MIN(JJJ(I),JJJ(J))
11     CONTINUE
C     WRITE(6, 201) CH
C201  FORMAT('   CH  before coef2    ', I3)
C-------------------------IL->XI
       CALL COEF2(GAM,DIR,IS,KSU1,XI,IQ,JJJ,
     * L1,S1,CH,KK,ISIGM1,ISIGM2,IG,MU)
      NU=CH
111   WRITE(6, 6)NU,IG
    6 FORMAT('   NO OF COEFFICIENTS       ',I5,
     * /     '   NO OF EXCHAN. COEF.      ',I5)
      CALL VAR11(NE5,Z,R,BET,H,ALFA,RO1,R1)
C--------------+XI-
      CALL COMONR( Z,R,BET,H,GAM,EPS,ALFA,
     *NE,IN,IL,IQ,JJJ,XI,DIR,NU,KK,ISIGM1,ISIGM2,IS,L1,SS1)
      WRITE(6, 100)DZ,ZK
  100 FORMAT(' CHARGE CORRECTION',2X,E23.8/
     *       ' NO OF STEPS      ',I12/)
      WRITE(6, 777)R1(3),R1(NEN+2)
777   FORMAT(' R MIN',22X,E15.8/' R MAX CALC',17X,E15.8/)
C--------------IL-XI---------------
      CALL ZEROR( NE5N,Z,1,IS,IQ,XI,IN,MU,
     * TETA,R5,R8)
C---------------IL-XI=====----------
      CALL DHR(ALFA,BET,H,Z,NE5N,IS,NU,
     *IN,IL,IQ,JJJ,XI,UD,MU,ADD,
     *IG,GAM,KK,ISIGM1,ISIGM2,
     *EPS,TETA,ALM3,
     *C3,ALA,ALA1,ALA2,ALM6DU,
     *KSU1,DZ,ZK,ETOT,
     *R1,R2,R3,R4,R5,R6,R8,R9,R10,R11,R12,R13,R14,R15)
      WRITE(6, 1985)ETOT
1985  FORMAT(' TOTAL ATOM(ION) ENERGY=  ',E15.8,
     * '  RY')
      DAT(1)=DFLOAT(KSU1)
      DAT(2)=R
      DAT(3)=DFLOAT(NE)
      DAT(4)=H
      DAT(5)=BET
      DAT(6)=EPS
      DAT(7)=Z
      DAT(8)=DFLOAT(IS)
      DAT(9)=ALFA
      DAT(10)=RO1
      DAT(11)=ETOT
      WRITE(10)(DAT(I),I=1,15)
      DO 50 I=1,IS+IS+IS+3
C     NI=NE5*(I-1)
      CALL DISKO1(10,I,7,NE5,R8(1,I))
  50   CONTINUE
      DO 60 I=1,IS
      WRITE(6,61)IS,I
 61   FORMAT(' IS ',I3,'  I ',I3)
      CALL COLUMN(NE5,R8(1,I),R5)
      CALL COLUMN(NE5,R8(1,I+IS+1),R13)
      CALL COLUMN(NE5,R8(1,I+IS+IS+2),R9)
C **corrections of V.Yarzh*****************
       XX=-R5(1)
       CALL KSII(NE,AL0,XX,R9,R12)
C      CALL KSII(NE,AL0,R5(1),R9,R12)
C***************************
      DO 70 J=1,NE
      R13(J)=R13(J)*AL0
      R5(J+2)=R5(J+2)*R12(J)
 70   CONTINUE
      CALL INTPLR( RO1,BET,I,ALFA,H,
     * IS,NE5,KT,IN,KSU1,0,0,R2,R5,R7)
      CALL INTPLR( RO1,BET,I,ALFA,H,
     * IS,NE5,KT,IN,KSU1,IS+1,-2,R2,R13,R7)
      CALL INTPLR( RO1,BET,I,ALFA,H,
     * IS,NE5,KT,IN,KSU1,IS+IS+2,-2,R2,R9,R7)
 60   CONTINUE
      WRITE(6,101)
 101  FORMAT('    R          large      small         V  ')
      DO 104 J=1,IS
      WRITE(6,105)J
      WRITE(*,105)J
 105  FORMAT(' SHELL ',I5)
      DO 102 I=1,NE,100
      WRITE(6,103)R1(I+2),R8(I+2,J),R8(I,J+IS+1),R8(I,J+IS+IS+2)
      WRITE(*,103)R1(I+2),R8(I+2,J),R8(I,J+IS+1),R8(I,J+IS+IS+2)
 103  FORMAT(4E12.3)
 102  CONTINUE
 104  CONTINUE
      RETURN
      END
