C=============================================================
C    PROGRAM FOR CALCULATION OF GENERALIZED OSCILLATOR
C          STRENGTHS IN HF AND RPAE
C          FOR ONE-SHELL ENDOHEDRAL ATOMS
C      NC   - NO OF CHANNELS
C      IZ-  HOLES IN CHANNELS
C      EXPI=0,LAMI1=0 - HF-ENERGIES, EXPI=1 -LAMI1=EXPER.ENERGIES
C      RPAE=0 -HF, RPAE=1 - HF,RPAE
C      Q1(1-HI)- IMPULSES
C      NE- NO OF POINTS IN RADIAL MESH
C      S- NO OF SHELLS
C      NW1 - MAX NUMBER OF FUNCTIONS IN FILE
C ======================================================
      SUBROUTINE GOSFUL(NC,HI,IZ,EXPI,LAMI1,CHAN,Y,DELY,TAU,
     *       LV,RPAE,NE,S,NW1,RF,V0,
     *       U1,NNLK,SPI,K1,KZ,
     *       I3,LVK,K,LI,LK,LVMK,CLS,B1,LAMI,IMIN,
     *       XMAXK,XK0,DELXK,Q1,
     *       LAMK,LADK,FME,FMEA,
     *       AL,BL,DK1,MA,
     *       C,Q,DV1,DV2,ARB,ARB1,ARB2,
     *       SM1,SM2,SMC1,SMC2,MB1,MB2,OMEGA,DM1,LAMM,
     *       AM,PH,EELRY,FF,PH1,
     *       R1,R2,R3,R4,R5,R6,R7,R9,R8)
C ======================================================
      INTEGER AS,QQ,S,EXPI,RPAE,HI,TAU,LV,
     *        LI(S),U1(S),SPI(S),NNLK(S),
     *        K1(NC),KZ(NC),SSK(NC),
     *        K(NC),IZ(NC),I3(NC),IMIN(NC),
     *        LVK(NC),LVMK(NC),LK(NC),CHAN(NC)
      REAL*8 Z,R,BET,H,ALFA,RO1,Y,DELY,RF,V0,
     *       OM,OM1,SA,A,B,A1,A2,A3,OM2,OM3,DET,
     *       D3,D4,D5,D6,P,P1,SB,BBK1,LAM1,LAML,
     *       BETR,BETN, Q1(HI),
     *       CLS(NC),B1(NC),LAMI(S),LAMI1(S),
     *       XMAXK(NC),XK0(NC),DELXK(NC),
     *       AL(*),BL((NC+1)*NC/2),
     *       FME(HI*NC*NW1),FMEA(HI*NC*NW1),
     *       C(*),Q(*),DV1(*),DV2(*),
     *       ARB(*),ARB1(*),ARB2(*),
     *       SM1(NC,NW1),SM2(NC,NW1),SMC1(NC,NW1),SMC2(NC,NW1),
     *       MB1(NC*NW1),MB2(NC*NW1),OMEGA(NC,NW1),
     *       DM1(NC*NW1),LAMM(NC*NW1),
     *       DK1(NC,*),LAMK(NC,*),LADK(NC,*),
     *       MA(*),DAT(15),AK(10),BK(10),
     *       AM(NW1),PH(NW1),EELRY(NW1),FF(*),PH1(*),
     *       R1(NE+15),R2(NE+15),R3(NE+15),R4(NE+15),
     *       R5(NE+15),R6(NE+15),R7(NE+15),R9(NE+15),R8(NE+15,*)

      QQ=0
      DO 1 I=1,NC
C       LLK(I)=0
      SSK(I)=0
      CLS(I)=1.
      B1(I)=0.5
      LVK(I)=2
   1  CONTINUE
      DO 5 I=1,NC*(NC+1)/2
      LVMK(I)=2
      AL(I)=1.
   5  BL(I)=0.5
       LLLK=LV
      J1=NC
      J3=NC
      J5=NC
      WRITE(6,10)NC
      WRITE(*,10)NC
  10  FORMAT('         ENDOHEDRAL ATOM   '/
     *       '   GENERALIZED OSCILLATOR STRENGTH '/
     *       '    WITH CORRELATIONS IN THE RPAE'/
     *       '     WITHIN',I3,2X,'TRANSITIONS')
      READ(10) (DAT(I),I=1,15)
      Z=DAT(7)
      NE5=NE+15
      R=DAT(2)
      H=DAT(4)
      BET=DAT(5)
      CALL VAR11(NE5,Z,R,BET,H,ALFA,RO1,R1)
      WRITE(6, 30)Z,S,NE,R,H,BET,ALFA,RF,V0
  30  FORMAT(' LS-COUPLING HF'/
     *       ' CHARGE OF NUCLEUS   ',F7.1/
     *       ' NUMBER OF SHELLS    ',I7/
     *       ' NUMBER OF POINTS    ',I7/
     *       ' R MAX               ',F7.1/
     *       ' STEP OF INTEGRATION ',E15.6/
     *       ' BETA                ',E15.6/
     *       ' ALFA                ',E15.6/
     *       ' FULLERENE RADIUS    ',E15.6/
     *       ' BINDING POTENTIAL   ',E15.6)

C=================================
C      WRITE(6, 31)SSS,LLLK,SSS,LLL
C      WRITE(*, 31)SSS,LLLK,SSS,LLL
C  31  FORMAT(
C     *       ' TERM OF  FINAL STATE  2S+1  L:',2I4/
C     *       ' TERM OF GROUND STATE  2S+1  L:',2I4)
       WRITE(6, 40)
       WRITE(*, 40)
  40  FORMAT(' GROUND STATE CONFIGURATION'/' ORBITAL   N ',
     *'     L     Q   SPI    ENERGY          IN ZONES')
      DO 50 I=1, S
      CALL DISKO1(10,I,2,NE5,R8(1,I))
          LAMI(I)=R8(1,I)
         LI(I)=R8(2,I)
         U1(I)=R8(NE+5,I)
         NNLK(I)=R8(NE+6,I)
         SPI(I)=R8(NE+8,I)
      WRITE(6, 51)I,U1(I),LI(I),NNLK(I),SPI(I),LAMI(I),I
      WRITE(*, 51)I,U1(I),LI(I),NNLK(I),SPI(I),LAMI(I),I
   51 FORMAT(' ',5I6,1X,E15.8,I6)
   50 CONTINUE

         DO 55 I=1,NC
         AS=IZ(I)
          IF(EXPI.NE.0) THEN
            LAMI(AS)=LAMI1(I)
          END IF
           LAMI1(I)=LAMI(AS)
      WRITE(6,56)I,U1(AS),LI(AS),NNLK(AS),SPI(AS),LAMI(AS)
      WRITE(*,56)I,U1(AS),LI(AS),NNLK(AS),SPI(AS),LAMI(AS)
   56 FORMAT('  CHANNEL ',I6,' HOLE NLQS',4I3,' E= ',E15.8)
   55    CONTINUE
C ========GET WAVEFUNCTIONS OF EXCITED STATES==========
         NN1=0
         AS=S
         DO 60 I=1,NC
         LNO=11+I-1
         REWIND LNO
         READ (LNO) (DAT(J),J=1,15)
C         DAT(13)=3
         K(I)=DAT(13)+DAT(12)
         NN1=NN1+K(I)
         K1(I)=DAT(12)
         KZ(I)=AS+1
         XK0(I)=DAT(14)
         DELXK(I)=DAT(15)
      WRITE(6, 70)I
      WRITE(*, 70)I
  70  FORMAT('=========================='/
     *       '  EXCITED STATE   CHANNEL ',I6/
     *       '         ENERGIES            N/PHASE',
     *       '        IN ZONES')
      DO 80 J=1,K(I)
      LA=AS+J
      CALL DISKO1(LNO,J,2,NE5,R8(1,LA))
      R8(NE+2,LA)=0.0D 00
      LAMK(I,J)=R8(1,LA)
      DK1(I,J)=R8(NE+7,LA)
      IA=7
      IF(R8(1,LA).LT.0.D 00)  IA=5
      WRITE(6, 90)J,R8(1,LA),R8(NE+IA,LA),LA
      WRITE(*, 90)J,R8(1,LA),R8(NE+IA,LA),LA
  90  FORMAT('  ',I2,2X,E15.8,5X,E15.8,I6)
  80  CONTINUE
      LK(I)=R8(2,LA)
      AS=LA
      WRITE(6,100) R8(2,LA),R8(NE+8,LA)
      WRITE(*,100) R8(2,LA),R8(NE+8,LA)
  100 FORMAT(' WITH ORBITAL MOMENTUM ', F4.1/
     *       ' AND SPIN               ',F4.1)
C====== f-factor===ph-phase======================
      CALL FFACTOR1(NE,AS,LNO,NW1,K(I),ALFA,BET,
     * RF,V0,AM,PH,EELRY,R1,R2,R8)
       DO 85 II=1,K(I)
       FF(II+Z1)=AM(II)
       PH1(II+Z1)=PH(II)
  85   CONTINUE
       Z1=Z1+K(I)
   60  CONTINUE
       IO=NC*NW1
       DO 95 II=1,IO
C      WRITE(6, 96)II,FF(II),PH1(II)
C  96  FORMAT('II FF PH ',I5,2E15.8)
  95   CONTINUE

C   60  CONTINUE
C  === CALCULATIONS OF  MATRIX ELEMENTS=======
      WRITE(6, 101)LV
      WRITE(*, 101)LV
 101  FORMAT('=============================='/
     *       '  ORBITAL MOMENTUM TRANSFERED ',I3)
         N6=1
      DO 102 I=1,HI
      Q1(I)=Y
      Y=Y+DELY*TAU
      DELY=DELY*TAU
 102  CONTINUE
      N7=0
C     NN1=(KZ(NC)+K(NC)-S-1)*2
      NN1=NN1*2
      DO 106 II=1,HI
      Y=Q1(II)
      WRITE(6, 107)Y
      WRITE(*, 107)Y
 107  FORMAT(' MOMENTUM TRANSFERED TO ATOM ',F10.6)
      DO 110  I=1,NC
          WRITE(6, 108)I
          WRITE(*, 108)I
 108   FORMAT('================================='/
     *        '  CHANNEL ',I3,' MATRIX ELEMENTS ')
         N1=K(I)
         N2=KZ(I)
         N3=IZ(I)
         N4=LI(N3)
         N5=LK(I)
C        N6=(KZ(I)-S)
         N8=I-1
       WRITE(6,111)
  111  FORMAT('    E,Ry      R-FORM    V-FORM')
         CALL MAEL3(ALFA,BET,NE5,H,N1,N4,N5,LV,Y,N2,N3,
     *        FME(N6),FMEA(N6),R1,R2,R5,R6,R8)
      DO 112 J=1,N1
  112    WRITE(6,119)LAMK(I,J),FME(N6+J-1),FMEA(N6+J-1)
  119    FORMAT(3E15.8)
      IF(SPI(I).NE.0) CLS(I)=1.D 00/DSQRT(2.D 00)
      DO 113 J=1,N1
         FME(N6+J-1)=FME(N6+J-1)*CLS(I)
  113    FMEA(N6+J-1)=FMEA(N6+J-1)*CLS(I)
        N6=N6+K(I)
  110   CONTINUE
  106   CONTINUE
        IF(RPAE.EQ.0) GOTO 121
       NN5=4.*NW1*NW1*NC*NC/NE5+2
       ILV=LV
C      IF(LV.EQ.0)ILV=1
C =====COLUMB MATRIX ELEMENTS============
        DO 120  I=1,NC
         N1=K(I)
         N2=KZ(I)
         N3=IZ(I)
         N4=LI(N3)
         N5=LK(I)
C        N6=(KZ(I)-S)
         N8=I-1
        DO 130  J=1,NC
         N1=K(I)
          N10=IZ(J)
           IF(I.EQ.J) THEN
            NN2=N1
            NN3=LV
          BBK1=B1(I)
          AL(I)=AL(I)*CLS(I)*CLS(I)
          WRITE(6,131)I,J
          WRITE(*,131)I,J
  131   FORMAT(' CHANNEL  ',2I3)
C 131  FORMAT('   MAELO2=DIAGONAL',4I3)
C====CALCULATIONS OF COULOMB ELEMENTS MATRIX========
      CALL MAELO2(ALFA,BET,NE5,H,NN1,LV,N1,N1,S,
     *                  IZ(I),IZ(J),LI(N3),LK(I),
     *                  LI(N10),LK(J),LAMI(N3),LAMI(N10),
     *                  AL(I),BK,BBK1,MA,KZ(I),KZ(J),
     *                      SSK(I),SSK(J),LLLK,ILV,
     *                  R1,R2,R3,R4,R5,R6,R7,R9,R8)
               ELSE
              NN2=K(J)
               IF(J.GT.I) THEN
                NN4=N7+J-I
                   ELSE
                  NN4=N8
                   N8=N8+NC-1-J
               END IF
              NN3=LVMK(NN4)-1
       IF(SPI(N3).NE.SPI(N10)) THEN
         BBK1=0.D 00
         BK(1)=0.D 00
           ELSE
           BBK1=B1(J)
           BK(1)=BL(NN4)
         END IF
           AK(1)=AL(NN4+NC)*CLS(I)*CLS(I)
C         PRINT 132,I,J,KZ(I),KZ(J)
C 132  FORMAT('  MAELO2=NONDIAGONAL',4I3)
       CALL MAELO2(ALFA,BET,NE5,H,NN1,LV,N1,NN2,S,
     *        IZ(I),IZ(J),LI(N3),LK(I),
     *        LI(N10),LK(J),LAMI(N3),LAMI(N10),
     *        AK,BK,BBK1,MA,KZ(I),KZ(J),
     *            SSK(I),SSK(J),LLLK,ILV,
     *        R1,R2,R3,R4,R5,R6,R7,R9,R8)
            END IF
 130       CONTINUE
         N7=N7+NC-I
 120  CONTINUE
C     ==========================================

 121     DO 140 I=1,NC
         J=K(I)
         J1=K1(I)
         IF(J.GT.J1+1) THEN
           XMAXK(I)=DSQRT(LAMK(I,J))
             ELSE
             XMAXK(I)=1.D-05
             XK0(I)=1.D-05
             DELXK(I)=1.D-05
           END IF
  140  CONTINUE
C   ===========================================
           N9=0
           OM2=0.1
           OM3=0.2
      DO 150 I=1,NC
         IF(DABS(LAMI(IZ(I))).LT.OM2.AND.K(I).GT.K1(I)) THEN
           N9=N9+1
           I3(N9)=I
           END IF
  150 CONTINUE
C=======CALCULATIONS OF GOS IN RPAE==========
      N6=1
      N66=1
      DO 160 II=1,HI
      Y=Q1(II)
      WRITE(6, 163)Y
      WRITE(*, 163)Y
 163  FORMAT('=============================='/
     *       '  MOMENTUM TRANSFERED TO ATOM ',F10.6)
      DO 170 JN=1,NC
         WRITE(6, 172)JN
         WRITE(*, 172)JN
 172  FORMAT(/'=============================='/
     *          '    CHANNEL ',I2,
     *          '    MATRIX ELEMENTS WITH CORREL.'/
     *  '    PHOTON     RE_L',5X,'IM_L',7X,'RE_V',9X,'IM_V  DET'/
     *          '    ENERGY ')
         JJ1=QQ
         IF(JN.LE.NC) JJ1=K(JN)
         NF=0
         IF(JN.EQ.NC+1) THEN
           DO 180 I=1,NC
              DO 190 J=1,K1(I)
  190            LAMK(I,J)=LADK(I,J)
  180      CONTINUE
           END IF
C              =================
         DO 200 JJ=1,JJ1
            IC1=0
            IF(JN.LE.NC) THEN
              JN1=JN
              A1=LAMK(JN,JJ)
C             A2=LAMI(IZ(JN))
              A2=LAMI1(JN)
              A3=LAMK(JN,JJ+1)
              OM=A1-A2
              J2=JJ
              OM1=OM
C     PRINT 201,OM,JN
C 201 FORMAT(' OM  JJ =',E15.4,I4)
              IF(JJ.GT.K1(JN)) OM=OM+(-1)**NF*5.D-02*(A3-A1)
              IF(A1.LT.0.D 00) IC1=JJ
              GOTO 210
              END IF
            IC1=0
            I2=I3(JN-NC)
            JN1=I2
            LAM1=LAMI1(I2)
            KI=K(I2)
            KI1=K1(I2)
            OM=OM2+(OM3-OM2)*JJ/DFLOAT(QQ+1)
            OM1=OM
            A=DABS(OM-LAMK(I2,KI1+1)+LAM1)
            DO 220 I=KI1+1,KI
               B=DABS(OM-LAMK(I2,I)+LAM1)
               IF(B.LE.A) THEN
                 J2=I
                 A=B
                 END IF
  220          CONTINUE
            A1=LAMK(I2,J2)
  210         IF(RPAE.EQ.0)  GOTO 230
           IF(JJ.EQ.IC1) THEN
              J3=0
              CALL DISNN(OM,OM1,A,SA,J2,NC,JN1,NN1,K,K1,
     *        IMIN,LAMI1,LAMK,XK0,DELXK,XMAXK,LLLK,R8,R8(1,NN5),MA,
     *        C,Q,DV1,DV2,ARB,ARB1,ARB2)
        WRITE(6, 214)OM,SA,A
  214   FORMAT(' OM NEW ',F10.5,' A,1-B ',2F10.5)
              IF(DABS((OM-OM1-SA)/OM1).GT.1D 00.AND.J3.LT.5) THEN
                J3=J3+1
                GOTO 210
                END IF
         WRITE(6,212)
  212    FORMAT(30(' -'))
              END IF
            CALL CQDN(OM,0.D 00,NC,JN1,K,K1,XK0,DELXK,
     *                XMAXK,LAMI1,LAMK,J2,C,Q,DV1,DV2)
            CALL INVMANFUL(K,K1,XK0,DELXK,XMAXK,LAMI1,
     *                  IMIN,OM,LLLK,J2,NC,JN1,NN1,C,Q,DV1,DV2,FF,
     *                  MA,DET,R8(1,1),R8(1,NN5),ARB,ARB1,ARB2)
            IF(J2.LE.K1(JN).AND.JN.LE.NC) LADK(JN,J2)=OM+LAMI1(JN)
            CALL MECO(K,NC,IMIN(JN1),FME(N66),FMEA(N66),R8(1,1),
     *                R8(1,NN5),D3,D4,D5,D6,NN1)
       WRITE(6, 215)OM1,D3,D5,D4,D6,DET
       WRITE(*, 215)OM1,D3,D5,D4,D6,DET
 215   FORMAT(' ',6F10.5)
 230           CONTINUE
C           IZ1=KZ(JN1)-NC
C           N6=(KZ(JN)-S)
            P=FME(J2+N6-1)**2*FF(J2+N6-1)**2
            P1=FMEA(J2+N6-1)**2*FF(J2+N6-1)**2
            CALL GOSD(LI(IZ(JN)),LV,NNLK(IZ(JN)),P,P1,
     *                OM1,Y,SA,SB)
C            PRINT 216,N6,N66,FME(J2+N6-1),P,P1,Y
C 216   FORMAT(' N6',2I6,4F10.6)
            SM1(JN,JJ)=SA
            SM2(JN,JJ)=SB
            P=(D3**2+D5**2)*FF(J2+N6-1)**2
            P1=(D4**2+D6**2)*FF(J2+N6-1)**2
            IF(IC1.NE.0) THEN
C             P=P/A
C             P1=P1/A
C     WRITE(6,520) OM,OM1
              OM1=OM
              END IF
C     WRITE(6,520) LLL
C 520 FORMAT(' LLL=',I6)
            CALL GOSD(LI(IZ(JN)),LV,NNLK(IZ(JN)),P,P1,
     *               OM1,Y,SA,SB)
C      PRINT 520,IZ(JN),LI(IZ(JN)),LV,NNLK(IZ(JN))
C 520  FORMAT('IZ(JN),LI,LV,NNLK,   '/4I4)
            SMC1(JN,JJ)=SA
            SMC2(JN,JJ)=SB
            BETR=0.0D 00
            BETN=0.0D 00
            MB1(JJ)=BETR
            MB2(JJ)=BETN
            OMEGA(JN,JJ)=OM1
            IF(JN.LE.NC) THEN
              LAMM(JJ)=LAMK(JN,JJ)
              DM1(JJ)=DK1(JN,JJ)
                ELSE
                LAMM(JJ)=OM1+LAM1
                DM1(JJ)=0.0D 00
              END IF
      IF(JJ.GT.K1(JN).AND.JN.LE.NC) NF=NF+1
  200  CONTINUE
C ====TABLE OF RESULTS =========================
      LAML=LAMI1(JN)
      IF(JN.GT.NC) LAML=LAM1
      N10=IZ(JN)
      WRITE(6, 250)Z,LAMI1(JN),U1(N10),LI(N10),LK(JN),SPI(N10)
      WRITE(*, 250)Z,LAMI1(JN),U1(N10),LI(N10),LK(JN),SPI(N10)
 250  FORMAT(//5X,'Z=',F10.5,3X,'TRANSITION  E=',E15.8/
     *    12X,' N=',I3,' L=',I3,'-->L=',I2/
     *    25X,'SPI=',I3)
      IF(NC.EQ.1) GOTO 260
      WRITE(6, 270)
 270  FORMAT(/12X,'    INTERACTION BETWEEN CHANNELS')
 260  CONTINUE
      WRITE(6, 280)
      WRITE(*, 280)
 280  FORMAT('  TRANSFERRED      GENER.OSCILLATOR'/
     *       '    ENERGY           STRENGTH      '/
     *       '      E,RY       HF-L      HF-V ',
     *     5X,'  RPAE-L      RPAE-V ')
      DO 290  JJ=1,JJ1
      IF(RPAE.NE.1) THEN
      SMC1(JN,JJ)=0
      SMC2(JN,JJ)=0
      DM1(JJ)=0
      END IF
      WRITE(6, 300)OMEGA(JN,JJ),SM1(JN,JJ),SM2(JN,JJ),
     *     SMC1(JN,JJ),SMC2(JN,JJ)
      WRITE(*, 300)OMEGA(JN,JJ),SM1(JN,JJ),SM2(JN,JJ),
     *     SMC1(JN,JJ),SMC2(JN,JJ)
      IF(JJ.EQ.K1(JN)) WRITE(6,305)
 305    FORMAT(30(' -'))
 290  CONTINUE
 300  FORMAT(2X,F10.5,4(1X,(E11.4)),
     * 1X,F8.4)
      N6=N6+K(JN)
 170  CONTINUE
       WRITE(6, 366)Y                           
 366   FORMAT(' MOMENTUM TRANSFERED TO ATOM',F10.6)
C
      IF(RPAE.EQ.1.AND.NC.GT.1) THEN
       WRITE(6, 364)
 364   FORMAT('    RANDOM PHASE APPROXIMATION')
         CALL CROSST(NC,NW1,K,K1,OMEGA,SMC1,SMC2,CHAN,
     *    R1,R2,R3,R4,R5,R6)
       END IF
      IF(RPAE.EQ.0.AND.NC.GT.1) THEN
       WRITE(6, 365)
 365   FORMAT('    HARTREE-FOCK APPROXIMATION')
         CALL CROSST(NC,NW1,K,K1,OMEGA,SM1,SM2,CHAN,
     *    R1,R2,R3,R4,R5,R6)
        END IF
      N66=N6
 160  CONTINUE
      RETURN
      END