C    PROGRAM FOR CALCULATION OF ATOMIC WAVE FUNCTIONS
C    IN HARTREE-FOCK APPROXIMATION
C    AUTHORS: L.V.CHERNYSHEVA,N.A.CHEREPKOV
C    MODIFIED BY L.V.CHERNYSHEVA,V.RADOJEVIC
C    B.V.GUL"TSEV,M.M.DOVGANICH,S.K.SEMENOV
      SUBROUTINE SCFHF1(KT,R7,Z,IS,R,NE,BET,H,IN,IL,IQ,L1,SS1,
     *     NU,IG,KK,ISIGM1,ISIGM2,GAM,EPS,MU,KSU1,DZ,ZK,
     *    IP,UD,DIR,
     *   TETA,ALM3,C3,ALA,ALA1,ALA2,ALM6DU,
     *  R1,R2,R3,R4,R5,R6,R8)
      REAL*8ALFA,BET,H,Z,RO1,EPS,ETOT,R,
     *GAM(1  ),TETA(IS),DZ,CO,
     *ALM3(IS), MU,ALM6DU(IS),
     * R1( NE+15),R2(NE+15 ),R3( NE+15),
     * R4( NE+15),R5( NE+15),R6( NE+15),
     * R7( KT),R8((NE+15)*IS),
     * ALA(IS),C3(IS),ALA1(IS),ALA2(IS),DAT(15)
      INTEGER CH,KK( 1),ISIGM1( 1),
     *ISIGM2(1 ),IN(IS),IL(IS),IQ(IS),DIR(IS),UD(IS),
     *SS1,S1(1),JJ,A,IP(IS),L1,ZK
      NE5=NE+15
      CO=0.7
      NEN=NE*CO
      NE5N=NEN+15
      WRITE(6, 1)
      WRITE(*, 1)
    1 FORMAT(' GROUND STATE ',
     *'ATOM (ION) WAVE FUNCTIONS'/
     *' IN THE HARTREE-FOCK APPROXIMATION')
      WRITE(6, 4)IS
    4 FORMAT('   NO OF N-L SHELLS         ',I5)
      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
      A=2*IL(I)+1
      IF(S1(1).LE.99) A=A*2
      DO 11 J=I1,IS
      JJ=2*IL(J)+1
      IF(S1(1) .LE.99)  JJ=JJ*2
      IF(I.LT.J ) CH=CH+     ( IL(I)+IL(J)-
     *IABS(IL(I)-IL(J)))/2 +1
      IF(I.EQ.J.AND.IQ(I).NE.1) CH=CH+IL(I)
      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+ MIN0(IL(I),IL(J))
11     CONTINUE
       CALL COEF(GAM,DIR,IS,KSU1,IL,IQ,
     * 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)
      CALL COMON( Z,R,BET,H,GAM,EPS,ALFA,
     *NE,IN,IL,IQ,DIR,NU,KK,ISIGM1,ISIGM2,IS,L1,S1)
      WRITE(6, 100)DZ,ZK
      WRITE(*, 100)DZ,ZK
  100 FORMAT(' CHARGE CORRECTION',2X,E23.8/
     *       ' NO OF STEPS      ',I12/)
      WRITE(6, 777)R1(3),R1(NEN)
777   FORMAT(' R MIN',22X,E15.8/' R MAX CALC',17X,E15.8/)
      CALL ZERO( R5,Z,R8,TETA,IS,IN,IQ,
     * NE5N,KSU1)
      CALL DH2(ALFA,BET,H,Z,EPS,
     *GAM,TETA,ALM3,R1,R2,R3,R4,R5,R6,R8,MU,
     * C3,ALA,ALA1,ALA2, ALM6DU,
     *NE5N,IS,NU,IG,KK,ISIGM1,
     *ISIGM2,IN,IL,IQ,DIR,UD,KSU1,DZ,ZK,ETOT)
      WRITE(6, 1985)ETOT
1985  FORMAT(1X,'TOTAL ATOM(ION) ENERGY =',E15.8,
     * '  RY')
      DO 12 I=IS,1,-1
      NI=NE5*(I-1)
      NIN=NE5N*(I-1)
      DO 15 J=15,4,-1
   15 R8(NI+NE+J)=R8(NIN+NEN+J)
      DO 13 J=NEN+3,1,-1
   13 R8(NI+J)=R8(NIN+J)
      DO 14 J=NEN+4,NE+3
   14 R8(NI+J)=0.D 00
   12 CONTINUE
      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
      NI=NE5*(I-1)
      CALL DISKO1(10,I,7,NE5,R8(NI+1))
  50   CONTINUE
      CALL INTPOL( R5,R8,RO1,BET,ALFA,R7,H,R2,
     * IS,NE5,KT,IN,KSU1)
      RETURN
      END
