! N. R. BADNELL  OMSUMPS            UoS                         08/01/19
!       BASED-ON XTRCT V1.8
!
!
      PROGRAM MAIN
!
      IMPLICIT REAL*8 (A-H,O-Z)
!
      LOGICAL BFORM
!
! SUM CROSS SECTIONS OVER FINAL STATES LYING ABOVE USER INPUT
! IONIZATION POTENTIAL TEAPOT, FOR INITIAL STATE IS.
! FOR USE WITH RMPS OMEGA TO CALCULATE IONIZATION WITHOUT PROJECTION.
!
! TEAPOT IS RELATIVE TO THE GROUND STATE.
!     ON INPUT, IF .GT. 0 IT IS IN RYD, IF .LT.0 THEN Z-SCALED.
!
! DEFAULT ASSUMES NEUTRAL OMEGAS INCLUDE ELASTIC, IONS NOT.
! SET IS INDEX = -IS IF NEUTRALS DO NOT CONTAIN ELASTIC OR IONS DO.
!
      PARAMETER (ZERO=0.0D0)
!
      ALLOCATABLE :: ENAT(:),LAT(:),ISAT(:),OMEGA(:),NRD(:)
!
      call cpu_time(timei)
!
      INQUIRE (FILE='OMEGA',EXIST=BFORM)
      IF(BFORM)THEN                           !USE FORMATTED OMEGA
        WRITE(*,*)'USING FORMATTED OMEGA FILE'
        OPEN(7,FILE='OMEGA',STATUS='OLD',FORM='FORMATTED')     !INPUT FILE
        OPEN(9,STATUS='SCRATCH',FORM='FORMATTED')
      ELSE
        WRITE(*,*)'USING UNFORMATTED OMEGA FILE'
        OPEN(7,FILE='OMEGAU',STATUS='OLD',FORM='UNFORMATTED')  !INPUT FILE
      ENDIF
!
      OPEN(8,FILE='xoutps',STATUS='UNKNOWN') !OUTPUT E(RYD), IONIZATION
!
      WRITE(*,*)'INPUT INITIAL STATE AND I.P.: '
      READ(*,*)IS,TEAPOT
!
      IF(IS*TEAPOT.EQ.0)STOP 'SET "IS" AND "TEAPOT" NON-ZERO'
!
!  READ-IN BASIC DATA
!
      IF(BFORM)READ(7,*)NZED,NELC
      IF(.NOT.BFORM)READ(7)NZED,NELC
!
      TAZ=NZED-NELC
      IF(NZED.EQ.NELC)TAZ=1.0
      TAZ=TAZ*TAZ
!
      IF(BFORM)READ(7,*)NAST,MXE,NOMWRT
      IF(.NOT.BFORM)READ(7)NAST,MXE,NOMWRT
!
      IF(IS.LT.0)THEN
        IF(NZED.EQ.NELC)THEN
          NELC=-NELC
        ELSE
          NELC=NZED
        ENDIF
        IS=-IS
      ENDIF
!
      IF(NZED.EQ.NELC)THEN
        NOMT=(NAST*(NAST+1))/2
        IONE=0
        NOT=(NAST*(NAST-1))/2
      ELSE
        NOMT=(NAST*(NAST-1))/2
        IONE=1
        NOT=(NAST*(NAST+1))/2
      ENDIF
!
      IF(NOMWRT.NE.0)THEN
        NOMT=IABS(NOMWRT)
        IF(NOMT.EQ.NOT)                                                 &
     &  WRITE(*,*)'*** POSSIBLE USER FLAG ERROR OF ELASTIC TRANSITIONS?'
      ENDIF
!
      ALLOCATE (ENAT(NAST),LAT(NAST),ISAT(NAST),OMEGA(NOMT),NRD(MXE)    &
     &         ,STAT=IERR)
      IF(IERR.NE.0)THEN                              
        WRITE(0,*)'ALLOCATION FAILS FOR NAST,NOMT,MXE=',NAST,NOMT,MXE
        STOP 'ABORT'                                         
      ENDIF                                          
!
      IF(BFORM)READ(7,*)(ISAT(I),LAT(I),I=1,NAST)
      IF(BFORM)READ(7,*)(ENAT(I),I=1,NAST)
      IF(.NOT.BFORM)READ(7)(ISAT(I),LAT(I),I=1,NAST)
      IF(.NOT.BFORM)READ(7)(ENAT(I),I=1,NAST)
!
      IF(NOMWRT.GE.0)THEN
        DO I=1,MXE
          NRD(I)=NOMT
        ENDDO
      ELSE
        IF(BFORM)READ(7,*)E0                    !CHECK FIRST E
        IF(.NOT.BFORM)READ(7)E0                 !CHECK FIRST E
        BACKSPACE(7)
        I=1
        IF(E0.GE.ENAT(NAST))GO TO 1    !ALL OPEN - EARLY EXIT
        IF(BFORM)THEN
          I0=1
          NO=NOPEN(E0,ENAT,NAST,IONE,I0)
          NO=MIN(NO,NOMT)
          READ(7,*)E0,(DUM,N=1,NO)
          READ(7,*)E                   !=E0+EINCR
          EINCR=E-E0
          E0=E0-EINCR
          EINCH=EINCR*0.5
          NREC=2+(NO-1)/6
          DO N=1,NREC
            BACKSPACE(7)
          ENDDO
        ENDIF
        I0=1
        DO I=1,MXE
          IF(BFORM)THEN            !CONSTANT MESH
            E=E0+I*EINCR
            WRITE(9,700)E
            BACKSPACE(9)
            READ(9,700)E
            BACKSPACE(9)
          ELSE
            READ(7)E               !VARIABLE MESH
          ENDIF
          IF(E.GE.ENAT(NAST))GO TO 1
          NRD(I)=NOPEN(E,ENAT,NAST,IONE,I0)
          NRD(I)=MIN(NRD(I),NOMT)
        ENDDO
        I=MXE+1
   1    I0=I
        DO I=I0,MXE
          NRD(I)=NOMT
        ENDDO
        IF(.NOT.BFORM)THEN
          REWIND(7)
          NREC=4
          DO N=1,NREC
            READ(7)
          ENDDO
        ENDIF
      ENDIF
!
! ONLY SUM OMEGAS ABOVE IONIZATION LIMIT
!
      ETEST=TEAPOT
      IF(TEAPOT.LT.0)ETEST=-ETEST*TAZ
      WRITE(8,100)IS,ETEST,ETEST-ENAT(IS)
!
      JS=0
      DO I=NAST,1,-1
        ENAT(I)=ENAT(I)*TAZ
        IF(ETEST.LE.ENAT(I))JS=I
      ENDDO
      IF(JS.EQ.0)THEN
        WRITE(*,*)'NO FINAL STATES ABOVE STATED I.P.'
        STOP
      ENDIF
      JS=MAX(JS,IS+1)                           !ODD CASE
      WRITE(*,*)'FIRST STATE ABOVE I.P.:',JS,' AT',ENAT(JS),' RYD'
!
! READ OMEGA FILE
!
      DO I=1,MXE
!
        IF(BFORM)THEN
          READ(7,*)E,(OMEGA(N),N=1,NRD(I))       !E RELATIVE TO GROUND
        ELSE
          READ(7)E,(OMEGA(N),N=1,NRD(I))         !E RELATIVE TO GROUND
        ENDIF
!
        IF(NOMWRT.LT.0.AND.BFORM)THEN
          ET=E0+I*EINCR
          IF(ABS(E-ET).GT.EINCH)THEN
            IF(NRD(I).LT.NOMT)THEN
            WRITE(*,*)I,NRD(I),ET,E
            STOP'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS'
            ENDIF
          ENDIF
        ENDIF
!
        E1=E*TAZ                             !E RELATIVE TO GROUND
        E2=E1-ENAT(IS)                       !E RELATIVE TO INITIAL
        E3=E1-ETEST
!
! SUM OVER COLLISION STRENGTHS ABOVE I.P. AT THIS ENERGY
!
        OSUM=0
        IF(NOMWRT.GE.0)THEN
          DO J=JS,NAST
            N=J+NAST*(IS-1)-(IS*(IS-1+2*IONE))/2
            IF(N.LE.NRD(I))OSUM=OSUM+OMEGA(N)
          ENDDO
        ELSE
          DO J=JS,NAST
            N=((J-1)*(J-2*IONE))/2+IS
            IF(N.LE.NRD(I))OSUM=OSUM+OMEGA(N)
          ENDDO
        ENDIF
        IF(OSUM.GT.0)WRITE(8,200)E1,E2,E3,OSUM
!
      ENDDO
!
      DEALLOCATE (ENAT,LAT,ISAT,OMEGA,NRD,STAT=IERR)
      IF(IERR.NE.0)THEN                              
        WRITE(0,*)'DE-ALLOCATION FAILS FOR NAST,NOMT,MXE...'
        STOP 'ABORT'                                         
      ENDIF                                          
!
      call cpu_time(timef)
!
      time=(timef-timei)/60.
      write(*,777)time
 777  format(' CPU TIME=',f9.3,' MIN')
!
      STOP 'NORMAL FINISH TO SUM PS'
!
 100  FORMAT('#',I3,' IONIZATION, ENERGY=',2F12.6)
 200  FORMAT(1X,4(1PE14.6))
 700  FORMAT(1PE14.8,6(1PE11.3)/(14X,6(E11.3)))
!
      END
!
!************************************************
!
      INTEGER FUNCTION NOPEN(E,ENAT,NAST,IONE,I0)
      IMPLICIT REAL*8 (A-H,O-Z)
!
      DIMENSION ENAT(NAST)
!
      DO I=I0,NAST
        IF(E.LT.ENAT(I))GO TO 1
      ENDDO
      I=NAST+1
   1  I0=I-1
      NOPEN=(I0*(I0-2*IONE+1))/2
!
      RETURN
      END
