C=========PROGRAM FOR CALCULATION OF INTERMEDIATE ATOMIC
C======== SHELL PHOTOABSORPTION AND SINGLE PHOTOIONIZATION===
C======== CROSS SECTIONS ACCOUNTING FOR PHOTOELECTRON=====
C=========SELF ENERGY PART=============================
C======== AUTHORS: L.V.CHERNYSHEVA,M.YA.AMUSIA,K.L.TSEMECHMAN==
C=========FORTRAN VERSION: L.V.CHERNYSHEVA==may1992==========
C=============================================================
C                VARIABLES
C      NES   - NO OF CHANNALS
C      NE- NO OF POINTS IN RADIAL MESH
C      S- NO OF SHELLS
C      NW1 - MAX NUMBER OF FUNCTIONS IN FILE
C       NC    NUMBER OF LINES IN INPUT DATA
C     ============================================
      SUBROUTINE  PHOSEP(NES,NC,IZ,K,KZ,F,F1,FZ,M,M1,MZ,
     *      IT,IL,LV,EXPI,LAML1,NE,S,NW1,
     *      NNL,M3,M4,LLI,LAMM,LAMF,LAMK,LAM,ME,
     *      ME1,ME2,ME3,ME4,ME5,ME6,ME7,RP,RM,
     *      RE,IM,RED,IMD,
     *      MME1,MME2,MME4,MME5,
     *      KK1,RRR,RRM,SR,
     *      CC,QQ,DD1,DD2,ARB,ARB1,ARB2,
     *      LLK,SSK,
     *      I3,LMM,CLS,B1,IMIN,
     *      D3,D4,D5,D6,XMAXM,XM0,DELXM,
     *      LAMKK,LAMI1,FME,FMEA,
     *      AL,BL,MA,SPI,IAL,OMEGA,
     *      R1,R2,R3,R4,R5,R6,R7,R9,R8)
C =================================================
      INTEGER NE,NES,S,NK,K1,W1,W2,K2,
     *        JN,I,J,NNL(S),
     *        AS,W,NM,NF,EXPI,NW1,RPAE,
     *        M,M1,MZ(NES),
     *        LI,LL,K,KZ(NES),F,F1,FZ(NES),
     *        IT(NC),IL(NC),LV(NC),LK,LM,LF,
     *        M3(NES),M4(NES),LLI(S),KK1(S),
     *        LLK(NES),SSK(NES),I3(NES),LMM(NES),IMIN(NES),
     *        IZ(S),SPI(S),IAL(*)

      REAL*8 Z,R,BET,ALFA,H,A,B,OM,SUM1,RO1,
     *       SIMPS1,PI,Y1,AMP5,
     *       RE1,RE2,RE3,RE4,IM1,IM2,IM3,IM4,
     *       DELXF,DELXK,SUM,LAML1(NC),LAMI,LAML,
     *       OMEGA(NES,*),
     *       LAM(S),LAMF(NW1),LAMK(NW1+S),LAMM(NW1),
     *       ME(M,*),ME1(M,*),ME2(M,*),ME3(M,*),
     *       ME4(NW1),ME5(NW1),ME6(NW1),ME7(NW1),
     *       RP(NW1),RM(NW1),
     *       RE(NW1),IM(NW1),RED(NW1),IMD(NW1),
     *       MME1(M,F,*),MME2(M,F,*),
     *       MME4(M,*),MME5(M,*),
     *       RRR(NW1,*),RRM(NW1,*),SR(*),
     *       PR(20),CC(NW1),QQ(NW1),DD1(NW1),DD2(NW1),
     *       ARB(*),ARB1(*),ARB2(*),CLS(*),B1(*),
     *       D3(NC,*),D4(NC,*),D5(NC,*),D6(NC,*),
     *       XMAXM(NES),XM0(NES),DELXM(NES),
     *       LAMKK(NES,NW1),LAMI1(S),
     *       FME(NC*NW1),FMEA(NC*NW1),
     *       AL(*),BL(*),MA(*),
     *       R1(NE+15),R2(NE+15),R3(NE+15),R4(NE+15),
     *       R5(NE+15),R6(NE+15),R9(NE+15),
     *       R7(NE+15),R8(NE+15,*),DAT(15)
      PI=3.1416
      WRITE(6, 10)NC
      WRITE(*, 10)NC
  10  FORMAT('          SINGLE PHOTOIONIZATON         '/
     *       '     CROSS SECTION WITH SELF-ENERGY PART '/
     *       '     WITHIN',I3,2X,'TRANSITIONS'//)
      READ(11) (DAT(I),I=1,15)
C       WRITE(6,12)(DAT(J),J=1,15)
C 12  FORMAT(5E15.7)
      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, 20)Z,S,NE,R,H,BET,ALFA
  20  FORMAT('  CHARGE OF NUCLEUS   ',E15.8/
     *       '  NUMBER OF SHELLS    ',I5/
     *       '  NUMBER OF POINTS    ',I5/
     *       '  R MAX               ',E15.8/
     *       '  STEP OF INTEGRATION ',E15.8/
     *       '  BETTA               ',E15.8/
     *       '  ALFA                ',E15.8)
        IF ((DAT(8)-S).GT.1.0 ) GOTO 600
      WRITE(6, 70)
      WRITE(*, 70)
  70  FORMAT('============================='/
     *       '   GROUND STATE CONFIGURATION'/
     * '  ORBITAL  L      ENERGY     IN ZONES')
       DO 80 I=1, S
      CALL DISKO1( 9,I,2,NE5,R8(1,I))
       LAM(I)=R8(1,I)
        LLI(I)=R8(2,I)
        NNL(I)=R8(NE+6,I)
        SPI(I)=R8(NE+8,I)
      WRITE(6,90) I,LLI(I),NNL(I),LAM(I),I
  90  FORMAT(3I6,E15.8,I6)
  80  CONTINUE
C =============================================
C     ====GET WAVE FUNCTIONS OF EXCITED STATES====
      NN1=0
      AS=S
      DO 95 I=1, NES
      LNO=11+I-1
      REWIND LNO
      READ (LNO) (DAT(J),J=1,15)
      M3(I)=DAT(13)+DAT(12)
      M4(I)=DAT(12)
      XM0(I)=DAT(14)
      DELXM(I)=DAT(15)
      WRITE(6, 121)I
      WRITE(*, 121)I
  121 FORMAT('=========================='/
     *       '  EXCITED  STATES CHANNEL ',I6)
      WRITE(6, 126 )
  126 FORMAT('      ENERGY   PHASE OR N ',
     *        '  IN ZONE')
      DO 130 J1=1, M3(I)
      LA=AS+J1
       CALL DISKO1(LNO,J1,2,NE5,R8(1,LA))
      IA=7
      IF(R8(1,LA).LT.0.) IA=5
       WRITE(6, 140)J1,R8(1,LA),R8(NE+IA,LA),LA
  140  FORMAT(' ',I3,2F10.5,I6)
  130  CONTINUE
       AS=LA
C      LAMKK(I,1)=R8(1,LA)
       LMM(I)=R8(2,LA)
       WRITE(6, 150)R8(2,LA)
  150  FORMAT('  WITH ELECTRON ORBITAL MOMENTUM ',F4.1)
   95  CONTINUE

C ==== CALCULATIONS OF MATRIX ELEMENTS=======
       K2=100
       DO 96 JN=1,NC
C      FJN1=F1(JN)
C      FJN=F(JN)
C      KJN=K(JN)
       K2=MIN0(K2,K)
       LAMI1(1)=LAM(IL(JN))
       WRITE(6, 152)JN,LV(JN),K2
       WRITE(*, 152)JN,LV(JN),K2
  152  FORMAT('======DIAGRAM ',I3/
     *       '  TRANSFERED ORBITAL  MOMENTUM ',I3,' K2 ',I3)
       LAMI=LAM(IT(JN))
      LI=LLI(IT(JN))
      WRITE(6, 153)IT(JN),LI,LAMI,IT(JN)
 153  FORMAT(' =====TIME REVERSE TERM===',
     *   '  I- HOLE   ',I3/
     *   '  L       E        ZONE    '/,I3,F10.5,I7)
      IF( EXPI.EQ.0) THEN
      LAML=LAM(IL(JN))
      ELSE
      LAML=LAML1(JN)
      WRITE(6, 154)
 154  FORMAT (' EXPERIMENTAL ')
      END IF
      LL=LLI(IL(JN))
      WRITE(6, 155)IL(JN),LL,LAML,IL(JN)
 155  FORMAT(' =====GROUND STATE ORBITAL======',
     *   '  L-HOLE   ',I3/
     *   '  L       E        ZONE    '/,I3,F10.5,I7)
      WRITE(6, 157)
 157  FORMAT(' =====TIME FIRWARD TERM====',
     *   '  F-ELECTRON    '/
     *   '  L       E        ZONE    ')
      DO 160 I=1,F
      NF=FZ(JN)+I-1
      CALL COLUMN(NE5,R8(1,NF),R5)
      LF=R8(2,NF)
      LAMF(I)=R8(1,NF)
      WRITE(6, 161)I,LF,LAMF(I),NF
 161  FORMAT(2I3,F10.5,I7)
 160  CONTINUE
      IF (F.GT.F1) DELXF=
     *DSQRT(LAMF(F))-DSQRT(LAMF(F-1))
      WRITE(6, 162)
      WRITE(*, 162)
 162  FORMAT(' =====EXCITED STATES ORBITALS====',
     *   '  M-ELECTRON    '/
     *   '  L       E        ZONE    ')
      DO 170 I=1,M
      NM=MZ(JN)+I-1
       CALL COLUMN(NE5,R8(1,NM),R5)
      LM=R8(2,NM)
      LAMM(I)=R8(1,NM)
      WRITE(6, 173)I,LM,LAMM(I),R5(NE+5),NM
 173  FORMAT(2I3,2F10.5,I7)
 170  CONTINUE
       K1=0
       J=1
       DO 94 I=1,S
       IF( LLI(I).EQ.LK) THEN
       K1=K1+1
       KK1(J)=I
       LAMK(J)=LAM(I)
       J=J+1
       END IF
  94  CONTINUE
C==================================
       WRITE(6, 97)
       WRITE(*, 97)
 97    FORMAT('=====K - OUTCOMING ELECTRON  '/
     *   '  L       E        ZONE    ')
      DO 200 W=1,K+K1
       IF (W.LE.K1) THEN
        NK=KK1(W)
       ELSE
        NK=KZ(JN)+W-1-K1
       END IF
       CALL COLUMN(NE5,R8(1,NK),R5)
       LK=R8(2,NK)
       LAMK(W)=R8(1,NK)
      IF (W.EQ.K+K1) DELXK=
     *DSQRT(LAMK(W))-DSQRT(LAMK(W-1))
      WRITE(6, 201)W,LK,LAMK(W),R5(NE+5),NK
  201 FORMAT(2I3,2E15.8,I3)
       CALL MAEL4(ALFA,BET,NE5,H,NK,FZ(JN),F,MZ(JN),M,
     *   IL(JN),LL,LF,LM,LK,
     * LV(JN),R1,PR,ME3,ME1,ME2,ME3,R8,R2,R5,R6,R7)
      DO 180 I=1,F
      DO 181 J=1,M
      MME1(J,I,W)=ME1(J,I)
      MME2(J,I,W)=ME2(J,I)
  181  CONTINUE
  180  CONTINUE
        WRITE(*, 211)W
 211  FORMAT('===================================='/
     *    '  COULOMB MATRIX ELEMENTS (F*M)  '/
     *    '=====TIME FORWARD===========',I4)
C       DO 210 I=1,F(JN)
C       DO 212 J=1,M(JN)
C       WRITE(6, 211,I,J,ME1(I,J),ME2(I,J),ME3(I,J)
C 211  FORMAT(2I3,' DIR.',F11.5,' EXCH.',F11.5,' COMB.',F11.5)
C 212   CONTINUE
C 210   CONTINUE
      IF (LV(JN).LE.LI+LK.AND.LV(JN).GE.IABS(LK-LI)) THEN
       CALL MAEL5(ALFA,BET,NE5,H,IL(JN),IT(JN),NK,MZ(JN),M,
     *  LL,LI,LK,LM,
     * LV(JN),R1,PR,ME,ME4,ME5,ME6,ME7,R8,R2,R5,R6,R7)
C       WRITE(6, 221
C 221   FORMAT('=====TIME REVERSE (M)======')
       DO 220 I=1,M
       MME4(I,W)=ME4(I)
       MME5(I,W)=ME5(I)
C       WRITE(6, 223,I,ME4(I),ME5(I),ME6(I)
C 223   FORMAT(I6,' DIR.',F11.5,' EXCH.',F11.5,' COMB.',F11.5)
 220   CONTINUE
       END IF
 200   CONTINUE
 96    CONTINUE
C======================================
       DO 202 W=K1+1,K2-1
       DO 203 I=1,K2
       DO 203 U1=1,K2
       RRR(U1,I)=0.D 00
       RRM(U1,I)=0.D 00
  203 CONTINUE
       WRITE(6, 322)W,LAMK(W)
       WRITE(*, 322)W,LAMK(W)
 322   FORMAT('====OUTCOMING ELECTRON ENERGY  ',I4,E15.8)
       DO 205 JN=1,NC
C      FJN=F(JN)
C      FJN1=F1(JN)
C      MJN=M(JN)
C      MJN1=M1(JN)
       DO 204 W2=1,K2
       DO 206 W1=1,K2
      DO 320 I=1, F
      DO 310 U1=1, M
       RP(U1)=MME1(U1,I,W2)*MME1(U1,I,W1)*2
 310  CONTINUE
       OM=LAMK(W)-LAMF(I)+LAML
      CALL ONE1(OM,LAMM,RP,M,M1,RE1,IM1,CC,QQ,DD1,RM)
C      WRITE(6, 316,OM,LAMM,RP,M,M1,RE1,IM1
C 316   FORMAT(' OM,LAMM ',5E15.8/
C    *        ' RP     ',4E15.8/' M,M1 ',2I3/'REIM',2E15.8)
       RE(I)=RE1
       IM(I)=IM1
      IF ( I.LT.F1 ) GOTO 315
       A=DSQRT(LAMF(I))*2.
       RED(I-F1)=RE(I)*A
       IMD(I-F1)=IM(I)*A
C       WRITE(6, 325,RED(I-F1),IMD(I-F1)
C 325   FORMAT(' RES1 RE ',E15.8,' IM ',E15.8)
 315   CONTINUE
 320   CONTINUE
       IF (W1.EQ.1.AND.W2.EQ.1)THEN
       DO 326 I=1,F
      WRITE(6, 324)I,RED(I),IMD(I)
 324    FORMAT(' I ',I3,' RES1 RE ',E15.8,' IM ',E15.8)
  326 CONTINUE
       END IF
       SUM= SIMPS1(F-F1,DELXF,RED,1)
       SUM1= SIMPS1(F-F1,DELXF,IMD,1)
      DO 330 I=1, F1
       SUM=SUM+RE(I)
       SUM1=SUM1+IM(I)
 330  CONTINUE
       B=1./DFLOAT(2*LV(JN)+1)/DFLOAT(2*LK+1)
C      WRITE(6, 323,B,SUM,SUM1
C323   FORMAT(' B,SUM,SUM1   '/3E15.8)
       RE1=SUM*B
       IM1=SUM1*B
C       WRITE(6, 335)RE1,IM1
C 335   FORMAT('==PHASE SHIFT CONTRIBUTIONS======'/
C    *        '====TIME FORWARD  DIRECT==== '/
C    *  ' RE1',E15.8,' IM1',E15.8)
      DO 350 I=1, F
C       WRITE(6, 352,I
C 352   FORMAT(' RES2 F=',I3)
      DO 340 U1=1, M
       RP(U1)=MME1(U1,I,W2)*MME2(U1,I,W1)
 340  CONTINUE
       OM=LAMK(W)-LAMF(I)+LAML
      CALL ONE1(OM,LAMM,RP,M,M1,RE2,IM2,CC,QQ,DD1,RM)
       RE(I)=RE2
       IM(I)=IM2
      IF ( I.LT.F1 ) GOTO 345
       A=DSQRT(LAMF(I))*2.
       RED(I-F1)=RE(I)*A
       IMD(I-F1)=IM(I)*A
 345   CONTINUE
       IF (W1.EQ.1.AND.W2.EQ.1)
     *  WRITE(*, 355)I,RED(I),IMD(I)
 355   FORMAT(' I',I3,' RES2 RE ',E15.8,' IM ',E15.8)
 350   CONTINUE
        SUM=SIMPS1(F-F1,DELXF,RED,1)
        SUM1=SIMPS1(F-F1,DELXF,IMD,1)
      DO 360 I=1, F1
       SUM= SUM+RE(I)
       SUM1=SUM1+IM(I)
 360   CONTINUE
       RE2=SUM*B
       IM2=SUM1*B
C       WRITE(*, 365)RE2,IM2
C 365   FORMAT('============EXCHANGE==== '/
C    *        ' RE2 ',E15.8,' IM2 ',E15.8)
       RE3=0.D 00
       RE4=0.D 00
      IF ( LV(JN).GT.LK+LI.OR.LV(JN).LT.IABS(LK-LI)) GOTO 395
      DO 370 I=1, M
       RP(I)=MME4(I,W2)*MME4(I,W1)*2.
 370  CONTINUE
       OM=LAMI-LAMK(W)+LAML
      CALL ONE1(OM,LAMM,RP,M,M1,RE3,IM3,CC,QQ,DD1,RM)
       RE3=-RE3*B
       IM3=-IM3*B
C      WRITE(*, 375)RE3,IM3
C 375  FORMAT('============================'/
C    *       '=====TIME REVERSE  DIRECT=== '/
C    *       ' RE3 ',E15.8,' IM3 ',E15.8)
      DO 380 I=1, M
       RP(I)=MME4(I,W2)*MME5(I,W1)
 380  CONTINUE
      CALL ONE1(OM,LAMM,RP,M,M1,RE4,IM4,CC,QQ,DD1,RM)
       RE4=-RE4*B
       IM4=-IM4*B
C       WRITE(*, 385)RE4,IM4
C 385   FORMAT(' ======EXCHANGE===== '/
C    *        ' RE4 ',E15.8,' IM4 ',E15.8)
 395   SUM=RE1-RE2+RE3-RE4
       SUM1=IM1-IM2+IM3-IM4
C       WRITE(6, 390)SUM,SUM1
C390  FORMAT('===== COMBINED PHASE SHIFT====== '/
C    *       ' RE ',E15.8,' IM ',E15.8)
      RRR(W2,W1)=RRR(W2,W1)+SUM
      RRM(W2,W1)=RRM(W2,W1)+SUM1
 206  CONTINUE
 204  CONTINUE
 205  CONTINUE
      CALL CQD1(LAMK(W),0.D 00,K2,K1,
     *     LAMK,CC,QQ,DD1)
C      WRITE(6, 399,W,LAMK(W),K2,K1,DELXK
C399   FORMAT(' W,LAMK(W),K2,K1,DELXK ',I3,E15.8,2I3,E15.8/
C     *       '   LAMK   CC   QQ   DD1 ')
C      DO 397 I=1,K2
C      WRITE(6, 398,LAMK(I),CC(I),QQ(I),DD1(I)
C398   FORMAT(2x,4E15.8)
C397   CONTINUE
      CALL INVMA1(W,NW1,K2,SR,Y1,RRR,RRM,CC,QQ,DD1,RE1,IM1)
      WRITE(6, 491)Y1,RE1,IM1
 491  FORMAT(' det ',E15.8,' RE1 ',E15.8,' IM1 ',E15.8)
      RE1=-RE1*PI
      IM1=IM1*PI
      RE2=-0.5*DASIN(RE1/DSQRT(RE1**2+(IM1+0.5)**2))
      IM2=-1./4.*DLOG(4.*(RE1**2+(0.5+IM1)**2))
      WRITE(6, 490)RE2,IM2
      WRITE(*, 490)RE2,IM2
 490  FORMAT(' PHASE SHIFT RE ',E15.8,' IM ',E15.8)
 202  CONTINUE
      STOP
        RPAE=1
C=======TO SOLVE RPAE-EQUATION=====================
      CALL RPADM(NC,ALFA,BET,H,RPAE,NE,S,NW1,
     *       LLLK,LLK,SSK,IZ,KZ,K,K1,
     *       LI,LK,LV,I3,CLS,B1,LAMI,IMIN,
     *       D3,D4,D5,D6,XMAXM,XM0,DELXM,
     *       LAMK,LAMI1,FME,FMEA,A,
     *       OMEGA,QQ,AL,BL,IAL,MA,SPI,
     *       CC,QQ,DD1,DD2,ARB,ARB1,ARB2,
     *       R1,R2,R3,R4,R5,R6,R7,R9,R8)
       DO 500 JN=1,NC
       DO 510 JJ=1,K
        CO=2.689*NNL(IZ(I))/(2*LLI(IZ(I))+1)
      AMP5=(D3(JN,JJ)**2+D5(JN,JJ)**2)*CO
 510  CONTINUE
 500  CONTINUE
      WRITE(6, 520)AMP5
 520  FORMAT(' RPAE CROSS SECTION ',E15.8)
 600  RETURN
      END
