!     omgmrg.f                                                  20/06/19
!
!     merge two omega files together, with non-overlapping energy ranges
!
!     default operation requires no user input.
!     but IF BOTH files contain infinite energy omega then set
!     ibige=1  to use it from file 1
!     ibige=-1 to use it from file 2
!     ibige=0 then it is in file 1 OR file 2 OR neither.
! Case MXE .gt.0 (e.g. R-matrix) default is ibige=0
! Case MXE .lt.0 (e.g. AS-DW) default is ibige=1  (i.e. in BOTH)
!
      PROGRAM MAIN
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!
      LOGICAL BFORM
!
      CHARACTER ELAS*3
!
      ALLOCATABLE :: ENAT1(:),LAT1(:),ISAT1(:),OMEG1(:),NRD1(:)
      ALLOCATABLE :: ENAT2(:),LAT2(:),ISAT2(:),OMEG2(:),NRD2(:)
!
      NAMELIST/SMERGE/IBIGE,ELAS
!
      call cpu_time(timei)
!
      open(6,file='routmrg',status='unknown')
!
      IBIGE=999
      ELAS='   '
!
      INQUIRE (FILE='domgmrg',EXIST=BFORM)
      IF(BFORM)THEN
        OPEN(4,FILE='domgmrg',STATUS='OLD')
        READ(4,SMERGE,END=10)
      ENDIF
!
   10 INQUIRE (FILE='omgmrg1',EXIST=BFORM)
      IF(BFORM)THEN                           !USE FORMATTED OMEGA
        WRITE(6,*)'USING FORMATTED OMEGA FILES'
        WRITE(0,*)'USING FORMATTED OMEGA FILES'
        open(1,file='omgmrg1',status='old')      !file1
        open(2,file='omgmrg2',status='old')      !file2
        open(3,file='omgmrgt',status='unknown')  !processed (merged) omega
        OPEN(9,STATUS='SCRATCH',FORM='FORMATTED')
      ELSE
        WRITE(6,*)'USING UNFORMATTED OMEGA FILES'
        WRITE(0,*)'USING UNFORMATTED OMEGA FILES'
        open(1,file='omgmrg1u',status='old',FORM='UNFORMATTED')      !file1
        open(2,file='omgmrg2u',status='old',FORM='UNFORMATTED')      !file2
        open(3,file='omgmrgtu',status='unknown',FORM='UNFORMATTED')  !processed (merged) omega
      ENDIF
!
      IF(BFORM)THEN
        read(1,*) nzed1,nelc1
        read(2,*) nzed2,nelc2
        write(3,*) nzed2,nelc2
        read(1,*) nast1,mxe01,nomt1
        read(2,*) nast2,mxe02,nomt2
!
        IF(IBIGE.EQ.999)THEN                                !SET DEFAULT
          IF(MXE01.GT.0)THEN
            IBIGE=0
          ELSE                                    !ASSUME AS-DW HAS BOTH
            IBIGE=1
          ENDIF
        ENDIF
        IF(IBIGE.GT.0)THEN
          WRITE(6,*)'USING OMGINF FROM FILE 1'
          WRITE(0,*)'USING OMGINF FROM FILE 1'
          IBIGE=1
        ENDIF
        IF(IBIGE.LT.0)THEN
          WRITE(6,*)'USING OMGINF FROM FILE 2'
          WRITE(0,*)'USING OMGINF FROM FILE 2'
          IBIGE=-1
        ENDIF
        isign=1
        if(mxe01.lt.0)isign=-1
        mxe1=iabs(mxe01)
        mxe2=iabs(mxe02)
        mxet=mxe1+mxe2
        if(ibige.ne.0)mxet=mxet-1
        write(3,*) nast2,mxet*isign,nomt2
      ELSE
        read(1) nzed1,nelc1
        read(2) nzed2,nelc2
        write(3) nzed2,nelc2
        read(1) nast1,mxe01,nomt1
        read(2) nast2,mxe02,nomt2
!
        IF(IBIGE.EQ.999)THEN                                !SET DEFAULT
          IF(MXE01.GT.0)THEN
            IBIGE=0
          ELSE                                    !ASSUME AS-DW HAS BOTH
            IBIGE=1
          ENDIF
        ENDIF
        IF(IBIGE.GT.0)THEN
          WRITE(6,*)'USING OMGINF FROM FILE 1'
          WRITE(0,*)'USING OMGINF FROM FILE 1'
          IBIGE=1
        ENDIF
        IF(IBIGE.LT.0)THEN
          WRITE(6,*)'USING OMGINF FROM FILE 2'
          WRITE(0,*)'USING OMGINF FROM FILE 2'
          IBIGE=-1
        ENDIF
        isign=1
        if(mxe01.lt.0)isign=-1
        mxe1=iabs(mxe01)
        mxe2=iabs(mxe02)
        mxet=mxe1+mxe2
        if(ibige.ne.0)mxet=mxet-1
        write(3) nast2,mxet*isign,nomt2
      ENDIF
!
      if(nzed1.ne.nzed2.or.nelc1.ne.nelc2) then
      write(6,*)' mis-match of omega files on nzed and nelc'
      write(6,*) 'nzed1 = ',nzed1,' nzed2 = ',nzed2
      write(6,*) 'nelc1 = ',nelc1,' nelc2 = ',nelc2
      stop 'mis-match of omega files on nzed or nelc'
      endif
!
      if(nast1.ne.nast2.or.nomt1.ne.nomt2) then
        write(6,*) 'mis-match of omega files on nast or nomt'
        write(6,*) 'nast1 = ',nast1,' nast2 = ',nast2
        write(6,*) 'nomt1 = ',nomt1,' nomt2 = ',nomt2
        stop 'mis-match of omega files on nast or nomt'
      end if
!
      IF(ELAS.EQ.'  ')THEN
        IF(NZED1.EQ.NELC1)ELAS='YES'
        IF(NZED1.NE.NELC1)ELAS='NO'
      ENDIF
      IF(ELAS.EQ.'YES')THEN
        IONE=0
      ELSEIF(ELAS.EQ.'NO')THEN
        IONE=1
      ELSE
        WRITE(6,*)' UNRECOGNIZED OPTION ELAS=',ELAS
        STOP ' UNRECOGNIZED ELAS OPTION'
      ENDIF
!
      NOMWRT=NOMT1                 !NOMT1=NOMT2
      NOMT1=ABS(NOMT1)
      NOMT2=ABS(NOMT2)
      IF(NOMT1*NOMT2.EQ.0)STOP 'NEED NOMWRT TO BE SPECIFIED'
!
      ALLOCATE (ENAT1(NAST1),LAT1(NAST1),ISAT1(NAST1),OMEG1(NOMT1)      &
     &         ,NRD1(MXE1),STAT=IERR)
      IF(IERR.NE.0)THEN                              
        WRITE(0,*)'ALLOCATION FAILS FOR NAST,NOMT,MXE=',NAST1,NOMT1,MXE1
        STOP 'ABORT 1'                                         
      ENDIF                                          
      ALLOCATE (ENAT2(NAST2),LAT2(NAST2),ISAT2(NAST2),OMEG2(NOMT2)      &
     &         ,NRD2(MXE2),STAT=IERR)
      IF(IERR.NE.0)THEN                              
        WRITE(0,*)'ALLOCATION FAILS FOR NAST,NOMT,MXE=',NAST2,NOMT2,MXE2
        STOP 'ABORT 2'                                         
      ENDIF                                          
!
      IF(BFORM)THEN
        read(1,*) (ISAT1(I),LAT1(I),i=1,nast1)
        read(1,*) (enat1(i),i=1,nast1)
!
        read(2,*) (ISAT2(I),LAT2(I),i=1,nast2)
        read(2,*) (enat2(i),i=1,nast2)
!
        write(3,90) (ISAT2(I),LAT2(I),i=1,nast2)
!        write(3,100) (enat2(i),i=1,nast2)
!
!        write(3,90) (ISAT1(I),LAT1(I),i=1,nast1)
        write(3,100) (enat1(i),i=1,nast1)
      ELSE
        read(1) (ISAT1(I),LAT1(I),i=1,nast1)
        read(1) (enat1(i),i=1,nast1)
!
        read(2) (ISAT2(I),LAT2(I),i=1,nast2)
        read(2) (enat2(i),i=1,nast2)
!
        write(3) (ISAT2(I),LAT2(I),i=1,nast2)
!        write(3) (enat2(i),i=1,nast2)
!
!        write(3) (ISAT1(I),LAT1(I),i=1,nast1)
        write(3) (enat1(i),i=1,nast1)
      ENDIF
!
      IFLAGM=0
      DO I=1,NAST1          !CHECK CONSISTENCY
        IF(IABS(ISAT1(I)).NE.IABS(ISAT2(I)).OR.LAT1(I).NE.LAT2(I))THEN
          WRITE(6,*)I,ISAT1(I),LAT1(I),'  ',ISAT2(I),LAT2(I)
          IFLAGM=IFLAGM+1
        ENDIF
! MIS-MATCH OF ENERGY O.K. IF ORDER DOESN'T CHANGE E.G. OBS VS. THEORY.
        IF(ABS(ENAT1(I)-ENAT2(I)).GT.1.0E-4*ENAT1(I))THEN
          WRITE(6,*)'ATTENTION:  ENERGY MIS-MATCH OF OMEGA FILES'
          WRITE(6,*)I,ENAT1(I),ENAT2(I)
        ENDIF
      ENDDO
      IF(IFLAGM.GT.0)THEN
          WRITE(6,*)' TARGET SYMMETRY MIS-MATCH OF OMEGA FILES FOR'     &
     &               ,IFLAGM,' PAIRS'
          STOP ' TARGET SYMMETRY MIS-MATCH OF OMEGA FILES'
      ENDIF
!
      IF(NOMWRT.GE.0)THEN                            !AND CASE MXE0.LT.0
        DO I=1,MXE1
          NRD1(I)=NOMT1
        ENDDO
        DO I=1,MXE2
          NRD2(I)=NOMT2
        ENDDO
      ELSE
        IF(BFORM)READ(1,*)E01                 !CHECK FIRST E
        IF(.NOT.BFORM)READ(1)E01                 !CHECK FIRST E
        BACKSPACE(1)
        I=1
        IF(E01.GE.ENAT1(NAST1))GO TO 1    !ALL OPEN - EARLY EXIT
        IF(BFORM)THEN
          I0=1
          NO=NOPEN(E01,ENAT1,NAST1,IONE,I0)
          NO=MIN(NO,NOMT1)
          READ(1,*)E01,(DUM,N=1,NO)
          READ(1,*)E                  !=E0+EINCR
          EINCR1=E-E01
          E01=E01-EINCR1
          EINCH1=EINCR1*0.5
          NREC=2+(NO-1)/6
          DO N=1,NREC
            BACKSPACE(1)
          ENDDO
        ENDIF
        I0=1
        DO I=1,MXE1
          IF(BFORM)THEN
            E=E01+I*EINCR1
            WRITE(9,110)E
            BACKSPACE(9)
            READ(9,110)E
            BACKSPACE(9)
          ELSE
            READ(1)E
          ENDIF
          IF(E.GE.ENAT1(NAST1))GO TO 1
          NRD1(I)=NOPEN(E,ENAT1,NAST1,IONE,I0)
          NRD1(I)=MIN(NRD1(I),NOMT1)
        ENDDO
        I=MXE1+1
   1    I0=I
        DO I=I0,MXE1
          NRD1(I)=NOMT1
        ENDDO
        IF(.NOT.BFORM)THEN
          REWIND(1)
          NREC=4
          DO N=1,NREC
            READ(1)
          ENDDO
        ENDIF
!
        IF(BFORM)READ(2,*)E02                 !CHECK FIRST E
        IF(.NOT.BFORM)READ(2)E02                 !CHECK FIRST E
        BACKSPACE(2)
        I=1
        IF(E02.GE.ENAT2(NAST2))GO TO 2    !ALL OPEN - EARLY EXIT
        IF(BFORM)THEN
          I0=1
          NO=NOPEN(E02,ENAT2,NAST2,IONE,I0)
          NO=MIN(NO,NOMT2)
          READ(2,*)E02,(DUM,N=1,NO)
          READ(2,*)E                  !=E0+EINCR
          EINCR2=E-E02
          E02=E02-EINCR2
          EINCH2=EINCR2*0.5
          NREC=2+(NO-1)/6
          DO N=1,NREC
            BACKSPACE(2)
          ENDDO
        ENDIF
        I0=1
        DO I=1,MXE2
          E=E02+I*EINCR2
          IF(BFORM)THEN
            WRITE(9,110)E
            BACKSPACE(9)
            READ(9,110)E
            BACKSPACE(9)
          ELSE
            READ(2)E
          ENDIF
          IF(E.GE.ENAT2(NAST2))GO TO 2
          NRD2(I)=NOPEN(E,ENAT2,NAST2,IONE,I0)
          NRD2(I)=MIN(NRD2(I),NOMT2)
        ENDDO
        I=MXE2+1
   2    I0=I
        DO I=I0,MXE2
          NRD2(I)=NOMT2
        ENDDO
        IF(.NOT.BFORM)THEN
          REWIND(2)
          NREC=4
          DO N=1,NREC
            READ(2)
          ENDDO
        ENDIF
      ENDIF
!
      it=0
!
!     first read energies and omegas from file1 and write them to file3
!
      IF(BFORM)THEN
        do i=1,mxe1-1
          read(1,*) emesh1,(omeg1(n),n=1,NRD1(I))
          IF(NOMWRT.LT.0)THEN
            ET=E01+I*EINCR1
            IF(ABS(EMESH1-ET).GT.EINCH1)THEN
              IF(NRD1(I).LT.NOMT1)THEN
              WRITE(6,*)'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS'
              WRITE(6,*)I,NRD1(I),ET,EMESH1
              STOP'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS'
              ENDIF
            ENDIF
          ENDIF
          write(3,110) emesh1,(omeg1(n),n=1,NRD1(I))
        enddo
        emesht=emesh1
        read(1,*) emesh1,(omeg1(n),n=1,NRD1(mxe1))
        if(emesh1.ge.1.e6)then
          if(ibige.eq.0)ibige=10
        else
          if(ibige.gt.0)then
            write(6,*)'*** Error, you have flagged infinite energy on'  &
     &                ,'file1 but no such data found'
            stop 'inf. energy point not in file1'
          elseif(ibige.lt.0)then
            write(6,*)'*** Warning, last energy dropped from file1'     &
     &                ,' because you flagged it as inf. energy, but'    &
     &                ,' this is not so'
          else
            emesht=emesh1
            write(3,110) emesh1,(omeg1(n),n=1,NRD1(mxe1))
          endif
        endif
!
!     then do the same for file 2
!
        do i=1,mxe2-1
          read(2,*) emesh2,(omeg2(n),n=1,NRD2(I))
          IF(NOMWRT.LT.0)THEN
            ET=E02+I*EINCR2
            IF(ABS(EMESH2-ET).GT.EINCH2)THEN
              IF(NRD2(I).LT.NOMT2)THEN
              WRITE(6,*)'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS'
              WRITE(6,*)I,NRD2(I),ET,EMESH2
              STOP'*** MIS-MATCH DURING READ OF FORMATTED OMEGAS'
              ENDIF
            ENDIF
          ENDIF
          if(i.eq.1) then
            if(emesh2.lt.emesht) then
              write(6,*) 'the minimum energy from file2 must be >= the'
              write(6,*) 'maximum energy from file1, but not the case'
              stop ' energy problems in files'
            end if
          end if
          write(3,110) emesh2,(omeg2(n),n=1,NRD2(I))
        enddo
        read(2,*) emesh2,(omeg2(n),n=1,NRD2(mxe2))
        if(emesh2.ge.1.e6)then
          if(ibige.le.0)then
            write(3,110) emesh2,(omeg2(n),n=1,NRD2(mxe2))
          elseif(ibige.eq.10)then
            write(6,*)'*** Error: inf. energy points in both files'     &
     &      ,' you must set ibige to flag which to use'
            stop '*** Inf. energy points found in both files!'
          endif
        else
          if(ibige.lt.0)then
            write(6,*)'*** Error, you have flagged infinite energy on'  &
     &                ,'file2 but no such data found'
            stop 'inf. energy point not in file2'
          elseif(ibige.gt.0)then
            write(6,*)'*** Warning, last energy dropped from file2'     &
     &                ,' because you flagged it as inf. energy, but'    &
     &                ,' this is not so'
          else
            write(3,110) emesh2,(omeg2(n),n=1,NRD2(mxe2))
          endif
        endif
        if(ibige.gt.0)write(3,110) emesh1,(omeg1(n),n=1,NRD1(mxe1))
!
!
      ELSE
!
!
        do i=1,mxe1-1
          read(1,err=3) emesh1,(omeg1(n),n=1,NRD1(I))
   3      write(3) emesh1,(omeg1(n),n=1,NRD1(I))
        enddo
        emesht=emesh1
        read(1) emesh1,(omeg1(n),n=1,NRD1(mxe1))
        if(emesh1.ge.1.e6)then
          if(ibige.eq.0)ibige=10
        else
          if(ibige.gt.0)then
            write(6,*)'*** Error, you have flagged infinite energy on'  &
     &                ,'file1 but no such data found'
            stop 'inf. energy point not in file1'
          elseif(ibige.lt.0)then
            write(6,*)'*** Warning, last energy dropped from file1'     &
     &                ,' because you flagged it as inf. energy, but'    &
     &                ,' this is not so'
          else
            emesht=emesh1
            write(3) emesh1,(omeg1(n),n=1,NRD1(mxe1))
          endif
        endif
!
!     then do the same for file 2
!
        do i=1,mxe2-1
          read(2,err=4) emesh2,(omeg2(n),n=1,NRD2(I))
          if(i.eq.1) then
            if(emesh2.lt.emesht) then
              write(6,*) 'the minimum energy from file2 must be >= the'
              write(6,*) 'maximum energy from file1, but not the case'
              stop ' energy problems in files'
            end if
          end if
   4      write(3) emesh2,(omeg2(n),n=1,NRD2(I))
        enddo
        read(2) emesh2,(omeg2(n),n=1,NRD2(mxe2))
        if(emesh2.ge.1.e6)then
          if(ibige.le.0)then
            write(3) emesh2,(omeg2(n),n=1,NRD2(mxe2))
          elseif(ibige.eq.10)then
            write(6,*)'*** Error: inf. energy points in both files'     &
     &      ,' you must set ibige to flag which to use'
            stop '*** Inf. energy points found in both files!'
          endif
        else
          if(ibige.lt.0)then
            write(6,*)'*** Error, you have flagged infinite energy on'  &
     &                ,'file2 but no such data found'
            stop 'inf. energy point not in file2'
          elseif(ibige.gt.0)then
            write(6,*)'*** Warning, last energy dropped from file2'     &
     &                ,' because you flagged it as inf. energy, but'    &
     &                ,' this is not so'
          else
            write(3) emesh2,(omeg2(n),n=1,NRD2(mxe2))
          endif
        endif
        if(ibige.gt.0)write(3) emesh1,(omeg1(n),n=1,NRD1(mxe1))
!
!
      ENDIF
!
!
      DEALLOCATE (ENAT1,LAT1,ISAT1,OMEG1,NRD1,STAT=IERR)
      IF(IERR.NE.0)THEN                              
        WRITE(0,*)'DE-ALLOCATION FAILS FOR NAST1,NOMT1,MXE1...'
        STOP 'ABORT'                                         
      ENDIF                                          
!
      DEALLOCATE (ENAT2,LAT2,ISAT2,OMEG2,NRD2,STAT=IERR)
      IF(IERR.NE.0)THEN                              
        WRITE(0,*)'DE-ALLOCATION FAILS FOR NAST2,NOMT2,MXE2...'
        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 MERGING'
!
   90 format(10(2i4))
  100 format(1P5e16.6)
  110 format(1pe14.8,6(e11.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
