! N. R. BADNELL  OMORDER/OMREDUCE       UoS                     09/01/19
!
      PROGRAM MAIN
!
      IMPLICIT REAL*8 (A-H,O-Z)
!
      LOGICAL BFORM,BSTG3
!
      CHARACTER ELAS1*3,ELAS2*3
!
! OMORDER/REDUCE CAN BE USED TO RE-ORDER AN OMEGA FILE, COMPUTED USING
! CALCULATED ENERGIES, SAY, INTO ONE ORDERED BY OBSERVED ENERGIES.
! IT CAN ALSO ELIMINATE TARGET STATES, E.G. CORRELATION, AND THE
! RELATED OMEGAS.
!
! THE USER INPUTS TWO COLUMNS OF ENERGY ORDERED TARGET INDEXES,
! THE FIRST IS THE CURRENT TARGET POSITION IN THE OMEGA FILE AND
! THE SECOND IS THE DESIRED POSITION, IF ZERO IT IS ELIMINATED
! ALONG WITH ALL THE TRANSITIONS ASSOCIATED WITH IT.
! INDEXES THAT REMAIN UNCHANGED (OR CHANGE ONLY BECAUSE OF THE
! OMISSION OF A LOWER LYING STATE)  NEED NOT BE SPECIFIED.
!
! INPUT IS TERMINATED BY EOF OR I .LE. ZERO IN THE FIRST COLUMN.
! IF THE FIRST COLUMN (TERMINATOR) IS .LT.0 THEN ORIGINAL TARGET STATE
! INDEXES GREATER THAN -I ARE OMITTED, E.G. CORRELATION, SO THERE IS
! NO NEED TO SPECIFY THEM ALL INDIVIDUALLY.
! IF THE SECOND COLUMN IS NON-ZERO, IT IS TAKEN TO BE THE NEW NOMWRT,
! BUT NAST IS NOT CHANGED, I.E. ONLY A SUBSET ARE WRITTEN.
!
! THE OMEGA FILE IS RE-ORDERED/REDUCED ACCORDING TO THIS SCHEME.
!
!
! *** IF RE-ORDERING IS TO TAKE PLACE AT CLOSED CHANNEL ENERGIES
!     IT IS HIGHLY RECOMMENDED THAT ANY STATES TO BE OMITTED ARE
!     ELIMINATED FIRST *WITHOUT* RE-ORDERING, I.E. NO DSTG3 FILE PRESENT
!     THEN A SECOND RUN IS MADE MAKING THE DSTG3 FILE AVAILABLE. ***
!     (IT IS ASSUMED THAT THE DSTG3 ENERGIES ARE IN THE ORIGINAL,
!      CALCULATED, ORDER AND THAT AFTER RE-ORDERING THEY ARE
!      ENERGY-ORDERED. CURRENTLY, THERE IS NO SORTING OF THE DSTG3
!      ENERGIES TO DETERMINE THE MAPPING.)
!
! IF NON-DEFAULT ELASTIC TRANSITIONS (ION/NEUTRAL NOT/PRESENT), 
! THEN SET ELAS1,2 IN UNIT5 /ORDER/ TO FLAG INPUT AND OUTPUT STATUS.
! THE NAMELIST DOES NOT NEED TO BE PRESENT, BUT IF IT IS THEN IT MUST
! BE AT THE START OF THE INPUT FILE dord.
!
      PARAMETER (ZERO=0.0)
      PARAMETER (ONE=1.0)
!
      ALLOCATABLE :: ENAT1(:),LAT1(:),ISAT1(:),OMEGA1(:)
      ALLOCATABLE :: ENAT2(:),LAT2(:),ISAT2(:),OMEGA2(:)
      ALLOCATABLE :: IORIG(:),JORIG(:),NRD(:),NWR(:)
!
      NAMELIST/ORDER/ELAS1,ELAS2
      NAMELIST/STG3A/RAD,ISORT
      NAMELIST/STG3B/NAST,INAST,TOLER
!
      IROW(ILI,ILF,IONE,NENG)=ILF+NENG*(ILI-1)-(ILI*(ILI-1+2*IONE))/2
      ICOL(ILI,ILF,IONE)=((ILF-1)*(ILF-2*IONE))/2+ILI
!
      call cpu_time(timei)
!
      OPEN(5,FILE='dord',STATUS='OLD',FORM='FORMATTED')
!
      ELAS1='   '
      ELAS2='   '
!
      READ(5,ORDER,END=50)
!
  50  REWIND(5)
!
      INQUIRE (FILE='OMEGA',EXIST=BFORM)
      IF(BFORM)THEN
        OPEN(7,FILE='OMEGA',STATUS='OLD',FORM='FORMATTED')          !INPUT
        OPEN(8,FILE='OMEGAN',STATUS='UNKNOWN',FORM='FORMATTED')     !OUTPUT
      ELSE
        OPEN(7,FILE='OMEGAU',STATUS='OLD',FORM='UNFORMATTED')       !INPUT
        OPEN(8,FILE='OMEGANU',STATUS='UNKNOWN',FORM='UNFORMATTED')  !OUTPUT
      ENDIF
      INQUIRE (FILE='dstg3',EXIST=BSTG3)
      IF(BSTG3)OPEN(9,FILE='dstg3',STATUS='OLD',FORM='FORMATTED')
!
      IF(BFORM)THEN
        READ(7,*)NZED,NELC 
        WRITE(8,*)NZED,NELC
        READ(7,*)NAST1,MXE,NOMWRT1
      ELSE
        READ(7)NZED,NELC
        WRITE(8)NZED,NELC
        READ(7)NAST1,MXE,NOMWRT1
      ENDIF
!
      AZ2=MAX(NZED-NELC,1)
      AZ2=ONE/(AZ2*AZ2)
!
      ALLOCATE (IORIG(NAST1))
!
      DO I=1,NAST1
        IORIG(I)=I
      ENDDO
!
      NOMWRT2=0
      DO N=1,NAST1
        READ(5,*,END=101)I,J
        IF(I.LE.0)GO TO 100
        IORIG(I)=J
      ENDDO
!
!  READ-IN BASIC DATA
!
 100  NOMWRT2=J
      IF(I.LT.0)THEN
        DO K=-I+1,NAST1
          IORIG(K)=0
        ENDDO
      ENDIF
!
!  DETERMINE NEW NAST
!
 101  NAST2=0
      DO I=1,NAST1
        IF(IORIG(I).GT.0)THEN
          NAST2=NAST2+1
        ENDIF
      ENDDO
!
      ALLOCATE (JORIG(NAST2))
!
!  SETUP NEW INDEX OF STATES
!
      DO 110 J=1,NAST2
        IZERO=0
        DO I=1,NAST1
          IF(IORIG(I).EQ.0)IZERO=IZERO+1
          IF(IORIG(I)-IZERO.EQ.J)THEN
            JORIG(J)=I
            GO TO 110
          ENDIF
        ENDDO
 110  ENDDO
!
      IF(ELAS1.EQ.'  ')THEN
        IF(NZED.EQ.NELC)ELAS1='YES'
        IF(NZED.NE.NELC)ELAS1='NO'
      ENDIF
      IF(ELAS1.EQ.'YES')THEN
        NOMT1=(NAST1*(NAST1+1))/2
        IONE1=0
      ELSEIF(ELAS1.EQ.'NO')THEN
        NOMT1=(NAST1*(NAST1-1))/2
        IONE1=1
      ELSE
        WRITE(6,*)' UNRECOGNIZED OPTION ELAS1=',ELAS1
        STOP ' UNRECOGNIZED ELAS1 OPTION'
      ENDIF
!
      IF(ELAS2.EQ.'  ')THEN
        IF(NZED.EQ.NELC)ELAS2='YES'
        IF(NZED.NE.NELC)ELAS2='NO'
      ENDIF
      IF(ELAS2.EQ.'YES')THEN
        NOMT2=(NAST2*(NAST2+1))/2
        IONE2=0
      ELSEIF(ELAS2.EQ.'NO')THEN
        NOMT2=(NAST2*(NAST2-1))/2
        IONE2=1
      ELSE
        WRITE(6,*)' UNRECOGNIZED OPTION ELAS2=',ELAS2
        STOP ' UNRECOGNIZED ELAS2 OPTION'
      ENDIF
!
      IF(NOMWRT1.NE.0.AND.ABS(NOMWRT1).NE.NOMT1)THEN
        WRITE(*,*)'*** WARNING: OMEGAN FILE MAYBE INCOMPLETE'
        NOMT1=ABS(NOMWRT1)
      ENDIF
!
!  INITIAL
!
      IF(NOMWRT2.EQ.0)THEN
        NOMWRT2=NOMT2
        IF(NOMWRT1.LT.0)NOMWRT2=-NOMWRT2
      ELSE
        IF(NOMWRT2*NOMWRT1.LT.0)THEN
          WRITE(*,*)'NOMWRT1,NOMWRT2 INCOMPATIBLE: ',NOMWRT1,NOMWRT2
          STOP
        ENDIF
        NOMT2=ABS(NOMWRT2)
      ENDIF
!
! CAN NOW ALLOCATE EVERYTHING
!
      ALLOCATE (ENAT1(NAST1),LAT1(NAST1),ISAT1(NAST1),OMEGA1(NOMT1)     &
     &         ,STAT=IERR)
      IF(IERR.NE.0)THEN
        WRITE(0,*)'ALLOCATION FAILS FOR NAST1,NOMT1=',NAST1,NOMT1
        STOP 'ABORT'
      ENDIF
!
      ALLOCATE (ENAT2(NAST2),LAT2(NAST2),ISAT2(NAST2),OMEGA2(NOMT2)     &
     &         ,STAT=IERR)
      IF(IERR.NE.0)THEN
        WRITE(0,*)'ALLOCATION FAILS FOR NAST2,NOMT2=',NAST2,NOMT2
        STOP 'ABORT'
      ENDIF
!
      ALLOCATE (NRD(MXE),NWR(MXE),STAT=IERR)
      IF(IERR.NE.0)THEN
        WRITE(0,*)'ALLOCATION FAILS FOR MXE=',MXE
        STOP 'ABORT'
      ENDIF
!
!  CONTINUE BASIC READS
!
      IF(BFORM)THEN
        READ(7,*)(ISAT1(I),LAT1(I),I=1,NAST1)
        READ(7,*)(ENAT1(I),I=1,NAST1)
      ELSE
        READ(7)(ISAT1(I),LAT1(I),I=1,NAST1)
        READ(7)(ENAT1(I),I=1,NAST1)
      ENDIF
!
!  OVERWRITE WITH DSTG3 ENERGIES
!
      IF(BSTG3)THEN
        READ(9,STG3A)
        NAST=0
        READ(9,STG3B)
        IF(NAST.GT.0)THEN
          NAST=MIN(NAST,NAST1)
          READ(9,*)(ENAT1(I),I=1,NAST)
          DO I=1,NAST
            ENAT1(I)=ENAT1(I)*AZ2
          ENDDO
        ENDIF
      ENDIF
!
!  REMAP
!
      DO I=1,NAST2
        J=JORIG(I)
        ISAT2(I)=ISAT1(J)
        LAT2(I)=LAT1(J)
        ENAT2(I)=ENAT1(J)
      ENDDO
!
!  WRITE BASIC
!
      IF(BFORM)THEN
        WRITE(8,*)NAST2,MXE,NOMWRT2
        WRITE(8,*)(ISAT2(I),LAT2(I),I=1,NAST2)
        WRITE(8,710)(ENAT2(I),I=1,NAST2)
      ELSE
        WRITE(8)NAST2,MXE,NOMWRT2
        WRITE(8)(ISAT2(I),LAT2(I),I=1,NAST2)
        WRITE(8)(ENAT2(I),I=1,NAST2)
      ENDIF
!
!  NOW SET-UP INDEX FOR READING OF OMEGAS
!
      IF(NOMWRT1.GE.0)THEN
        DO I=1,MXE
          NRD(I)=NOMT1
          NWR(I)=NOMT2
        ENDDO
      ELSE
        IF(BFORM)READ(7,*)E0                 !CHECK FIRST E
        IF(.NOT.BFORM)READ(7)E0                 !CHECK FIRST 
        BACKSPACE(7)
        I=1
        IF(E0.GE.ENAT1(NAST1).AND.E0.GE.ENAT2(NAST2))GO TO 3 !ALL OPEN - EARLY EXIT
        I0=1
        NO1=NOPEN(E0,ENAT1,NAST1,IONE1,I0)
        NO1=MIN(NO1,NOMT1)
        NO2=NOPEN(E0,ENAT2,NAST2,IONE2,I0)
        NO2=MIN(NO2,NOMT2)
        IF(MXE.EQ.1)THEN
          NRD(1)=NO1
          NWR(1)=NO2
          GO TO 7
        ENDIF
        IF(BFORM)THEN
          READ(7,*)E0,(DUM,N=1,NO1)
          READ(7,*)E                  !=E0+EINCR
          EINCR=E-E0
          E0=E0-EINCR
          EINCH=EINCR*0.5
          NREC=2+(NO1-1)/6
          DO N=1,NREC
            BACKSPACE(7)
          ENDDO
        ENDIF
        I0=1
        DO I=1,MXE
          IF(BFORM)THEN
            E=E0+I*EINCR
            WRITE(99,700)E
            BACKSPACE(99)
            READ(99,700)E
            BACKSPACE(99)
          ELSE
            READ(7)E
          ENDIF
          IF(E.GE.ENAT1(NAST1).AND.E.GE.ENAT2(NAST2))GO TO 3
          NRD(I)=NOPEN(E,ENAT1,NAST1,IONE1,I0)
          NRD(I)=MIN(NRD(I),NOMT1)
          NWR(I)=NOPEN(E,ENAT2,NAST2,IONE2,I0)
          NWR(I)=MIN(NWR(I),NOMT2)
        ENDDO
        I=MXE+1
   3    I0=I
        DO I=I0,MXE
          NRD(I)=NOMT1
          NWR(I)=NOMT2
        ENDDO
        IF(.NOT.BFORM)THEN
          REWIND(7)
          NREC=4
          DO N=1,NREC
            READ(7)
          ENDDO
        ENDIF
      ENDIF
!
!  READ AND WRITE OMEGAS
!
   7  DO I=1,MXE
!
        IF(BFORM)READ(7,*)EMESH,(OMEGA1(N),N=1,NRD(I))
        IF(.NOT.BFORM)READ(7)EMESH,(OMEGA1(N),N=1,NRD(I))
        IF(NOMWRT1.LT.0.AND.BFORM)THEN
          ET=E0+I*EINCR
          IF(ABS(EMESH-ET).GT.EINCH)THEN
            IF(NRD(I).LT.NOMT1)THEN
              WRITE(6,*)'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS'
              WRITE(6,*)I,NRD(I),ET,EMESH
              STOP'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS' 
            ENDIF
          ENDIF
        ENDIF
!
        IF(NOMWRT2.GT.0)THEN
          N=0
          DO ILI=1,NAST2-IONE2
            II=JORIG(ILI)
            DO ILF=ILI+IONE2,NAST2
              N=N+1
              IF(N.GT.NWR(I))GO TO 200
              JJ=JORIG(ILF)
              KK=IROW(MIN(II,JJ),MAX(II,JJ),IONE1,NAST1)
              IF(KK.LE.NRD(I))THEN
                OMEGA2(N)=OMEGA1(KK)
              ELSE
                OMEGA2(N)=ZERO
              ENDIF
            ENDDO
          ENDDO
        ELSE
          N=0
          DO ILF=1+IONE2,NAST2
            JJ=JORIG(ILF)
            DO ILI=1,ILF-IONE2
              N=N+1
              IF(N.GT.NWR(I))GO TO 200
              II=JORIG(ILI)
              KK=ICOL(MIN(II,JJ),MAX(II,JJ),IONE1)
              IF(KK.LE.NRD(I))THEN
                OMEGA2(N)=OMEGA1(KK)
              ELSE
                OMEGA2(N)=ZERO
              ENDIF
            ENDDO
          ENDDO
        ENDIF
!
  200   IF(BFORM)WRITE(8,700)EMESH,(OMEGA2(N),N=1,NWR(I))
        IF(.NOT.BFORM)WRITE(8)EMESH,(OMEGA2(N),N=1,NWR(I))
!
      ENDDO
!
      DEALLOCATE (ENAT1,LAT1,ISAT1,OMEGA1,STAT=IERR)
      IF(IERR.NE.0)THEN                              
        WRITE(0,*)'DE-ALLOCATION FAILS FOR NAST1,NOMT1...'
        STOP 'ABORT'                                         
      ENDIF
!
      DEALLOCATE (ENAT2,LAT2,ISAT2,OMEGA2,STAT=IERR)
      IF(IERR.NE.0)THEN                              
        WRITE(0,*)'DE-ALLOCATION FAILS FOR NAST2,NOMT2...'
        STOP 'ABORT'                                         
      ENDIF 
!
      DEALLOCATE (IORIG,JORIG,NRD,NWR,STAT=IERR)
      IF(IERR.NE.0)THEN
        WRITE(0,*)'DE-ALLOCATION FAILS FOR 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 OMORDER'
!
  700 FORMAT(1PE14.8,6(1PE11.3)/(14X,6(E11.3)))
  710 FORMAT(1P5E16.6)
!
      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
