C===============================================================
C    PROGRAM FOR CALCULATION OF ATOMIC WAVE FUNCTIONS
C---- IN THE HARTREE-FOCK "FROZEN" CORE
C     BY L.V.CHERNYSHEVA,N.A.CHEREPKOV, M.YA.AMUSIA,
C     V.F.ORLOV.S.K.SEMENOV,
C     FORTRAN VERSION: L.V.CHERNYSHEVA,S.K.SEMENOV,N.A.CHEREPKOV
C===============================================================
      SUBROUTINE FCHF(KT,R7,Z,IS,R,NE,BET,H,IN,IL,IQ,L1,SS1,
     *     NU1,IG,GAM,KK,ISIGM1,AKAP,DKAP,K1,K2,EPS,
     *     MU,DN,NR,IT,NU2,GAM1,GAM2,ORT,KSU1,
     *     ISIGM2,IP,S1,TETA,ALM3,ALA,C3,ALA1,ALA2,DIR,
     *     R1,R2,R3,R4,R5,R6,R9,R10,R8)
      INTEGER IS,IG,CH,SS1,DN,
     * IP( 1),ISIGM1( 1),ISIGM2( 1),KK( 1),
     * S1(IS),ORT(IS),IN(IS),IL(IS),IQ(IS),DIR(IS)
      REAL*8 Z,R,BET,H,EPS,AKAP,DKAP,MU,
     *   ALFA,RO1,RO2,B,C,D,P,EPS1,A1,
     *   ALFA1,RO1N,BET1,QQ4,
     *   DAT(15),GAM( 1),GAM1( 1),
     *   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),
     *   R8(NE+15,IS),R7(*)
      NE5=NE+15
      EPS1=1.D-05
      ALFA1=ALFA
      BET1=BET
      IF(L1.EQ.111) GOTO 40
      DO 20 I=1,IS
20    IP(I)=1
      CH=1
      I1=1
      S1(1)=SS1
      DO 30 I=1,IS
      JJA=2*IL(I)+1
      IF(S1(1).LE.99) JJA=JJA*2
      DO 30 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.JJA-1))
     * CH=CH+ MIN0(IL(I),IL(J))
 30    CONTINUE
       CALL COEF(GAM,DIR,IS,KSU1,IL,IQ,
     * L1,S1,CH,KK,ISIGM1,ISIGM2,IG,MU)
C     CH=23
      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) IS,NU1,IG
      WRITE(*,62) IS,NU1,IG
 62   FORMAT('    HARTREE-FOCK WAVE FUNCTIONS IN THE',
     *       ' FIELD OF "FROZEN" CORE',/
     * 21H NUMBER OF N-L SHELLS,11X,
     *I5/20H NO OF COEFFICIENTS ,12X, 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 9
      READ(9)(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
 140  FORMAT(' R MIN',22X,E15.8  )
 190  FORMAT(/'  CORE WAVE FUNCTIONS:  '//
     *      ' E                 N    L    Q')
 110  CONTINUE
 196  CALL VAR11(NE5,Z,R,BET,H,ALFA,RO1N,R1)
      CALL COMON3( Z,R,BET,H,GAM,GAM1,GAM2,EPS,ALFA,
     * NE,IN,IL,IQ,IG,NU,NU1,NU2,
     * IT,ORT,DIR,IP,KK,ISIGM1,ISIGM2,IS,L1,SS1,MU)
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)
      IRA=IS-1
      WRITE(6,190)
      DO 260 IW=1,IRA
      CALL  DISKO1(9,IW,2,NE5,R8(1,IW))
      MMOM=IFIX(SNGL(R8(2,IW)))
      NOBL=IFIX(SNGL(R8(NE+5,IW)))
      KOLE=IFIX(SNGL(R8(NE+6,IW)))
      R8(2,IW)=1.D 00
      R8(NE+3,IW)=0.D 00
      DO 270 J=1,NE
 270  R8(J+2,IW)=R8(J+2,IW)*DSQRT(R8(2,IW))
      WRITE(6,220) R8(1,IW),NOBL,MMOM,KOLE
 220  FORMAT( E15.8,3I5)
 260  CONTINUE
      REWIND 9
      NH=NE
 250  ISIGM=IS
      KK1=0
      KK2=0
      IF(NR.EQ.1) GOTO 290
      READ(10) (DAT(JJ1),JJ1=1,15)
      DO 300 IW=1,NR-1
      CALL DISKO1(10,IW,2,NE5,R8(1,IW+IS-1))
      WRITE(6,251) R8(1,IW+IS-1),R8(2,IW+IS-1)
251   FORMAT(' old wave function '/2E15.8)
C      R8(2,IS+IW+2)=XNOBL
C      R8(NE+7,IS+IW+2)=XPHS
 300  CONTINUE
 290  I=IS
      DAT(1)=KSU1
      DAT(7)=Z
      DAT(8)=DFLOAT(IS)
      QQ4=DFLOAT(IN(IS))
      KK2=DAT(12)
      DAT(12)=DFLOAT(K2)
      DAT(13)=DFLOAT(K1)
      IF(NR.EQ.1) GOTO 320
C     KK2=DAT(12)+0.5
      IF(NR.GT.KK2) GOTO 330
      DAT(12)=DAT(12)+DFLOAT(NR-1)
      GOTO 320
 330  DAT(12)=KK2+DAT(12)
      DAT(13)=DAT(13)+DFLOAT(NR-KK2-1)
 320  REWIND 11
      DAT(11)=0.D 00
      DAT(14)=AKAP
      DAT(15)=DKAP
      WRITE(11)(DAT(I),I=1,15)
      DO 350 IW=1,NR-1
C      XNOBL=R8(2,IS+IW-2)
C      XPHS=R8(NE+7,IS+IW-2)
      CALL DISKO1(11,IW,7,NE5,R8(1,IW+IS-1))
 350  CONTINUE
      IF(K2.LT.1) GOTO 360
      RO1=-BET*(10.D 00 + DLOG(Z))
      RO2=RO1+DFLOAT(NE)*H
      ALFA=(RO2- BET*DLOG(R))/R
C     DO 355 IT=1,IS-1
C355  IF(ORT(IT).EQ.0) IL(IT)=-IL(IT)
      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 ZERO( R5,Z,R8,TETA,
     * IS,IN,IQ,NE+15,IS)
      ISIGM=IS
C       ORT(6)=1
C       IF(IN(IS).EQ.4) ORT(6)=0
       DO 390 IT=1,IS
 390   ALM3(IT)=0.0
      CALL MAIN13(ALFA ,BET ,H ,Z,EPS,GAM,GAM1,GAM2,TETA,
     *R1,R2,R3,R4,R5,R6,R9,R10,R8,C3,ALA,ALA1,ALA2,
     *ALM3,MU,IS,NE5,IQ,ISIGM,IN,IL,
     *IG,NU,ISIGM1,ISIGM2,KK,DIR,ORT,KSU1,L1,SS1)
      IF(DABS(DAT(2)-R).LT.EPS1) GOTO 400
      C=0.
      R5(2)=DFLOAT(IL(IS))
      R5(NE+4)=0.0
      D=0.
      DO 410 J=1,NE
      P=D
      D=C
      C=(R5(J+2) *R1(J+2) /(ALFA1*R1(J+2)+BET1))**2
      A1=5.*C
      B=8.*D
      A1=(A1+B-P)*H/12.
 410  R5(NE+4)=R5(NE+4)+A1
      DO 420 J=1,NE
 420  R5(J+2)=R5(J+2)/DSQRT(R5(NE+4))
 400  CONTINUE
      XIS=DFLOAT(IN(IS))
      R5(NE+8)=DIR(IS)
      IA=J10+NR-1
      CALL DISKO1(11,IA,7,NE5,R5)
      CALL INTPOL( R5,R8,RO1N,BET ,ALFA ,R7,
     *H ,R2,IS,NE5,KT,IN,KSU1)
      IN(IS)=IN(IS)+DN
       DO 445 I=1,IS-1
        IQ1=IQ(I)+IQ(IS)
       IF(IN(I).EQ.IN(IS).AND.
     *     IL(I).EQ.IL(IS).AND.
     *     (IQ1.NE.(2*IL(IS)+1).OR.
     *      IQ1.NE.(4*Il(IS)+2)))
     *      IN(IS)=IN(IS)+1
  445    CONTINUE
C     IF(IN(IS).EQ.8)IN(IS)=9
       WRITE(*,371) IN(IS)
 371  FORMAT( ' IN',I5)
 370  CONTINUE
C     IN(IS)=IN(IS)-K2*DN
 360  CONTINUE
      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 ZERO( R5,Z,R8 ,TETA,
     * IS,IN,IQ,NE5,IS)
      ISIGM=IS
       DO 470 IT=1,IS
  470  ALM3(IT)=0.0
      CALL MAIN23(ALFA ,BET ,H ,Z,AKAP,EPS,
     *TETA,R1,R2,R3,R4,R5,R6,R9,R10,R8 ,ALA,ALA1,
     *ALA2,ALM3,C3,GAM,GAM1,GAM2,MU,IS,NE5,
     *IL,KK,ISIGM1,ISIGM2,IN,IQ,NU,IG,ISIGM,DIR,ORT,L1,SS1)
      R5(NE+8)=DIR(IS)
      IA=II+NR+K2-1
      CALL DISKO1(11,IA,7,NE5,R5)
      AKAP=AKAP+DKAP
      CALL INTPOL( R5,R8 ,RO1N,BET ,ALFA ,
     * R7,H ,R2,IS,NE5,KT,IN,KSU1)
450      CONTINUE
440      CONTINUE
       RETURN
       END
