C===============================================================
C    PROGRAM FOR CALCULATION OF ATOMIC WAVE FUNCTIONS
C     IN THE HARTREE-FOCK-DIRAC "FROZEN" CORE
C     BY L.V.CHERNYSHEVA, V.L.Yakhontov
C===============================================================
      SUBROUTINE RFCHF(KT,R7,Z,IS,R,NE,BET,H,IN,IL,IQ,XI,L1,SS1,
     *   NU1,IG,GAM,KK,ISIGM1,AKAP,DKAP,K1,K2,EPS,
     *   MU,DN,IT,NU2,GAM1,GAM2,PJ,ORT,KSU1,ADD,
     *   ISIGM2,IP,S1,TETA,ALM3,ALA,C3,ALA1,ALA2,SPI,
     *   R1,R2,R3,R4,R5,R6,R9,R10,R11,R12,R13,R14,R15,R16,R8)

      INTEGER IS,IG,SS1,DN,A,CH,
     *   IP(*),ISIGM1(*),ISIGM2(*),KK(*),
     *   S1(IS),ORT(IS),IN(IS),IL(IS),IQ(IS),SPI(IS),
     *   XI(IS)

      REAL*8 Z,R,BET,H,EPS,AKAP,DKAP,MU,AL0,ADD,
     *   ALFA,RO1,C,D,EPS1,ALFA1,RO1N,BET1,
     *   PJ(IS),DAT(15),GAM(*),GAM1(*),
     *   TETA(IS),ALM3(IS),ALA(IS),GAM2(IS),
     *   C3(IS),ALA1(IS),ALA2(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),R16(NE+15),R8(NE+15,*),R7(100)

      NE5=NE+15
      EPS1=1.D-05
      ALFA1=ALFA
      BET1=BET
      AL0=1/137.036
      DO 10 I=1,IS
      XI(I)=(-1)**(PJ(I)+IL(I)+0.5)*(PJ(I)+0.5)

C
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))
      WRITE(6, 206)XI(I)
 206  FORMAT('  XI           ',I3)
 10   CONTINUE
      IF(L1.EQ.111) GOTO 40
      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=IDINT(2.*PJ(J)+1.)
      IF(S1(1) .LE.99)  JJ=JJ*2
      IF(I.LT.J ) CH=CH+IFIX(SNGL((PJ(I)+PJ(J)-
     *DABS(PJ(I)-PJ(J)))/2+1))
      WRITE(6, 200)CH
 200   FORMAT(' CHG ',I3)
      IF(I.EQ.J.AND.IQ(I).NE.1) THEN
      CH=CH+IFIX(SNGL(PJ(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(PJ(I),PJ(J))
11     CONTINUE
      WRITE(6, 201) CH
 201  FORMAT('   CH  before coef2    ', I3)
       CALL COEF2(GAM,SPI,IS,KSU1,XI,IQ,PJ,
     * L1,S1,CH,KK,ISIGM1,ISIGM2,IG,MU)
      NU1=CH
      GOTO 60

 40   II=IS
      DO 50 I=1,NU1
 50   ISIGM2(I)=II
      ISIGM2(NU1+NU2+1)=IT
 60   WRITE(6,  62) ORT(IS),IS,NU1,IG
      WRITE(*,  62) ORT(IS),IS,NU1,IG
C-----------------------------------------------


 62   FORMAT('    HARTREE-FOCK-DIRAC WAVE FUNCTIONS '/
     *       '      IN THE FIELD OF "FROZEN" CORE',/
     * ' ORTOGONALIZATION SIGN     ',I5/
     * ' NUMBER OF N-L SHELLS      ',I5/
     * ' NO OF COEFFICIENTS        ',I5/
     * ' NO OF EXCH. COEFFICIENTS  ',I5)
      IF(ORT(IS).EQ.0) WRITE(6,63)
 63   FORMAT('  EXCHANGE BEFORE ITERAT.ORTOGONALIZATION')
      IF(ORT(IS).EQ.1) WRITE(6,64)
 64   FORMAT('  FORMULAE FOR ORTOGONALIZATION')
      IF(ORT(IS).EQ.3) WRITE(6,65)
 65   FORMAT(' ITERAT.ORTOGONALIZATION BEFORE EXCHANGE')
      REWIND 10
      READ(10)(DAT(JJ1),JJ1=1,15)
      IF(DAT(1).EQ.1.) GOTO 70
      WRITE(6,80)DAT(1)
 80   FORMAT(5X,' KSU#1,   KSU=  ',E8.0)
      RETURN
  70  IF(DABS(DAT(7)-Z).LT.EPS1.AND.
     *DABS(DAT(2)-R).LT.EPS1.AND.
     *DABS(DFLOAT(IS-1) -DAT(8)).LT.EPS1)  GOTO 110
      WRITE(6,100)(DAT(JJ1),JJ1=1,15),R,Z,IS,EPS1
 100  FORMAT(5X,' CHECK "CORE" AND "EXCITED" PARAMETERS'/
     *    '  KSU=       ',E15.8/
     *    '  R CORE     ',E15.8/
     *    '  N=         ',E15.8/
     *    '  H=         ',E15.8/
     *    '  BET=       ',E15.8/
     *    '  EPS CORE   ',E15.8/
     *    '  Z CORE=    ',E15.8/
     *    2(4(3X,E15.8)/)/
     *    ' R EXCITED   ',E12.5/
     *    ' Z EXCITED   ',E12.5/
     *    ' S  EXCITED  ',I3/
     *    ' EPS EXC.    ',E15.8)
       RETURN
C----------------------------------
 110  CALL VAR11(NE5,Z,R,BET,H,ALFA,RO1N,R1)
      CALL COMOR3( Z,R,BET,H,GAM,GAM1,GAM2,EPS,ALFA,
     * NE,IN,IL,IQ,XI,IG,NU,NU1,NU2,PJ,
     * IT,ORT,SPI,IP,KK,ISIGM1,ISIGM2,IS,L1,SS1)
C     STOP
      IF(ORT(IS).EQ.0.OR.ORT(IS).EQ.3) THEN
      ISIGM2(NU+2)=0
      ELSE
      ISIGM2(NU+2)=1
      END IF
      WRITE(6,140)R1(3)
 140  FORMAT(' R MIN',22X,E15.8  )
      WRITE(6, 266)ADD
      WRITE(*, 266)ADD
 266  FORMAT(' COEF A ',E15.8)
      IRA=IS+IS+IS
      WRITE(6,190)
 190  FORMAT(/'  CORE WAVE FUNCTIONS:  '//
     *      '        energy        N  L  Q   JJ     X    ')
      DO 260 IW=1,IRA
      CALL  DISKO1(10,IW,2,NE5,R8(1,IW))
 260  CONTINUE
      DO 265 IW=1,IS-1
       ILC=IDINT(R8(NE+4,IW))
       INC=IDINT(R8(NE+5,IW))
       IQC=IDINT(R8(NE+6,IW))
       IXC=IDINT(R8(2,IW))
         D=R8(NE+3,IW)
         C=R8(NE+2,IW)
C     DO 270 J=1,NE
C270  R8(J+2,IW)=R8(J+2,IW)*DSQRT(R8(2,IW))
      WRITE(*,220) IW,R8(1,IW),INC,ILC,IQC,D,IXC
      WRITE(6,220) IW,R8(1,IW),INC,ILC,IQC,D,IXC
 220  FORMAT(I3,2X,E15.8,3I3,3X,F3.1,2X,I3,2X,E15.8)
 265  CONTINUE
C     DO 267 IW=IS+1,IS+IS-1
C267  R8(J+2,IW)=R8(J+2,IW)*DSQRT(R8(2,IW))
C     REWIND 9
C     NH=NE
C250  ISIGM=IS
C     KK1=0
C     KK2=0
C290  I=IS
      DAT(1)=KSU1
      DAT(7)=Z
      DAT(8)=DFLOAT(IS)
C     QQ4=DFLOAT(IN(IS))
C     KK2=DAT(12)
      DAT(12)=DFLOAT(K2)
      DAT(13)=DFLOAT(K1)
C     IF(NR.EQ.1) GOTO 320
C     KK2=DAT(12)+0.5
C     IF(NR.GT.KK2) GOTO 330
C     DAT(12)=DAT(12)+DFLOAT(NR-1)
C     GOTO 320
C330  DAT(12)=KK2+DAT(12)
C     DAT(13)=DAT(13)+DFLOAT(NR-KK2-1)
C320  REWIND 11
      DAT(11)=0.D 00
      DAT(14)=AKAP
      DAT(15)=DKAP
      WRITE(11)(DAT(I),I=1,15)
C     DO 350 IW=1,NR-1
C      XNOBL=R8(2,IS+IW-2)
C      XPHS=R8(NE+7,IS+IW-2)
C     CALL DISKO1(11,IW,7,NE5,R8(1,IW+IS-1))
C350  CONTINUE
C     IF(K2.LT.1) GOTO 360
C     RO1=-BET*(10.D 00 + DLOG(Z))
C     RO2=RO1+DFLOAT(NE)*H
C     ALFA=(RO2- BET*DLOG(R))/R
      IA=IS+IS+IS+2
      IB=K2+K1+IA
      IC=(K2+K1)*2+IA
      ID=K1+K2
      DO 370 J10=1,K2
      WRITE(6,380)J10
      WRITE(*,380)J10
  380  FORMAT(60('=')//'  DISCRETE SPECTRUM WAVE FUNCTION',
     *'  IN THE  HF-APPROXIMATION',I5)

      CALL ZEROR( NE5,Z,IS,IS,IQ,XI,IN,MU,
     * TETA,R5,R8)
      ISIGM=IS
       DO 390 IT=1,IS
 390   ALM3(IT)=0.0
      CALL MAINRD(ALFA,BET,H,NE5,Z,AL0,ADD,
     *IS,MU,ISIGM,IN,XI,IQ,
     *NU,IG,GAM,KK,ISIGM1,ISIGM2,ORT,KSU1,
     *EPS,TETA,C3,ALA,ALA1,ALA2,ALM3,
     * GAM1,GAM2,SPI,
     *R1,R2,R3,R4,R5,R6,R9,R10,R11,R12,R13,R14,R15,R16,R8)
      R5(NE+3)=PJ(IS)
      R5(NE+4)=DFLOAT(IL(IS))
      R5(NE+5)=DFLOAT(IN(IS))
      R5(NE+6)=DFLOAT(IQ(IS))
      R5(NE+7)=DFLOAT(XI(IS))
      R5(NE+8)=DFLOAT(SPI(IS))
      CALL COLUMN(NE5,R5,R8(1,IA+J10))
      CALL COLUMN(NE5,R15,R8(1,IB+J10))
      CALL COLUMN(NE5,R9,R8(1,IC+J10))
      WRITE(6,350)IA,IB,IC,J10,R8(NE+3,IA+J10),R8(NE+3,IB+J10)
350   FORMAT('ia,ib,ic,j10,E_p E_G',4I3,2E15.8)
      DO 430 J=1,NE
      R5(J+2)=R5(J+2)*R13(J)
      R15(J)=R15(J)*AL0
 430  CONTINUE

      CALL INTPLR(RO1,BET,IS,ALFA,H,
     * IS,NE5,KT,IN,KSU1,0,0,R2,R5,R7)
      CALL INTPLR(RO1,BET,IS,ALFA,H,
     * IS,NE5,KT,IN,KSU1,IS+1,-2,R2,R15,R7)
      CALL INTPLR(RO1,BET,IS,ALFA,H,
     * IS,NE5,KT,IN,KSU1,IS+IS+2,-2,R2,R9,R7)
      IN(IS)=IN(IS)+DN
C     IF(IN(IS).EQ.8)IN(IS)=9
C      WRITE(*,371) IN(IS)
C371  FORMAT( ' IN',I5)
 370  CONTINUE
      IN(IS)=IN(IS)-K2*DN
 360  CONTINUE
C-----------------------
      IF(K1.LT.1) GOTO 440
      DO 450 II=1,K1
      WRITE(6,460) II
      WRITE(*,460) II
  460 FORMAT(60('=')//'  CONTINUOUS SPECTRUM WAVE'
     *,' FUNCTION IN THE  HF-APPROXIMATION',I5)

       CALL ZEROR( NE5,Z,IS,IS,IQ,XI,IN,MU,
     * TETA,R5,R8)
      ISIGM=IS
       DO 480 IT=1,IS
  480  ALM3(IT)=0.0
       CALL MAINRC(ALFA,BET,H,NE5,Z,AL0,AKAP,
     *     IS,MU,ISIGM,IN,XI,IQ,NU,IG,
     *     GAM,GAM1,GAM2,KK,ISIGM1,ISIGM2,SPI,ORT,
     *     EPS,TETA,ALA,ALA1,ALA2,ALM3,C3,
     *   R1,R2,R3,R4,R5,R6,R9,R10,R11,R12,R13,R14,R15,R16,R8)
C      R5(2)=DFLOAT(XI(IS))
      R5(NE+3)=PJ(IS)
      R5(NE+4)=DFLOAT(IL(IS))
      R5(NE+5)=0.
      R5(NE+6)=DFLOAT(IQ(IS))
      R6(NE+7)=DFLOAT(XI(IS))
      R5(NE+8)=DFLOAT(SPI(IS))
        CALL COLUMN(NE5,R5,R8(1,IA+K2+II))
       CALL COLUMN(NE5,R15,R8(1,IB+K2+II))
       CALL COLUMN(NE5,R9,R8(1,IC+K2+II))
      WRITE(6,355)IA,IB,IC,k2,ii,R8(NE+5,IA+K2+II),R8(NE+5,IA+K2+II)
355   FORMAT('ia,ib,ic,k2,ii,E',5I3,2E15.8)
      DO 490 J=1,NE
      R5(J+2)=R5(J+2)*R13(J)
      R15(J)=R15(J)*AL0
 490  CONTINUE
      CALL INTPLR(RO1,BET,IS,ALFA,H,
     * IS,NE5,KT,IN,KSU1,0,0,R2,R5,R7)
      CALL INTPLR(RO1,BET,IS,ALFA,H,
     * IS,NE5,KT,IN,KSU1,IS+1,-2,R2,R15,R7)
      CALL INTPLR(RO1,BET,IS,ALFA,H,
     * IS,NE5,KT,IN,KSU1,IS+IS+2,-2,R2,R9,R7)
      AKAP=AKAP+DKAP
450      CONTINUE
      AKAP=AKAP-DKAP*DFLOAT(K1)
      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)=0.
      DAT(12)=DFLOAT(K2)
      DAT(13)=DFLOAT(K1)
      DAT(14)=AKAP
      DAT(15)=DKAP
       REWIND 11
       WRITE(6,424)DAT(12),DAT(13),DAT(14),DAT(15)
       WRITE(*,424)DAT(12),DAT(13),DAT(14),DAT(15)
 424   FORMAT(' K2 K1,AKAP,DKAP ',4E15.8)
       WRITE(11)(DAT(I),I=1,15)
       DO 425 J10=1,ID+ID+ID
       CALL DISKO1(11,J10,7,NE5,R8(1,J10+IA))
       WRITE(6,426)J10+IA,R8(1,IA+J10),R8(NE+5,IA+J10)
 426   FORMAT('j10+IA ',I3,' E ', 2E15.8)
 425   CONTINUE

      WRITE(6,501)
 501  FORMAT('    R          large      small         V  ')
      DO 504 J=1,ID
      WRITE(6,505)J
      WRITE(*,505)J
 505  FORMAT(' SHELL ',I5)
      DO 502 I=1,100
      WRITE(6,503)R1(I+2),R8(I+2,J+IA),R8(I,J+IB),R8(I,J+IC)
 503  FORMAT(4E12.3)
 502  CONTINUE
 504  CONTINUE
440      CONTINUE
       RETURN
       END
