**==aa0001.f  
C
C PROGRAM 'DIFFLMRG'
C =================
C
C SORT AND MERGE A NUMBER OF 'DIFFLAUE' OUTPUT MTZ FILES
C FOR REFLECTIONS WITH MORE THAN ONE ENTRY, WEIGHTED AVERAGES
C ARE TAKEN FOR THE FL1 SIGL1 FL2 AND SIGL2 COLUMNS AND FOR
C THE PHASE AND WEIGHT COLUMNS (IF PRESENT) THE VALUES FOR
C THE FIRST ENTRY ARE TAKEN. THE PROGRAM ASSUMES THAT THE
C SAME COLUMNS ARE PRESENT IN EACH FILE AND THAT THEY ARE
C IN THE ORDER PRODUCED BY DIFFLAUE
C
C THE ARRAY IA(/AA) IS USED TO HOLD THE SETS OF LAUE DATA FOR SORTING
C AND MERGING. EACH REFLECTION OCCUPIES 7 WORDS AS FOLLOWS
C (WORD 1 INTEGER, REMAINDER REAL)
C        WORD 1  H,K,L,FILM PACKED  2**20(H+128)+2**12(K+128)+
C                                   2**4(L+128)+KFILM (KFILM=0 IN
C                                   THIS PROGRAM)
C        WORD 2  FL1
C        WORD 3  SIGL1
C        WORD 4  FL2
C        WORD 5  SIGL2
C        WORD 6  COLUMN 8 FROM INPUT MTZ FILE IF PRESENT
C        WORD 7  COLUMN 9 FROM INPUT MTZ FILE IF PRESENT
C
C PARAMETER STATEMENT  MXRF=MAX NO. OF REFLECTIONS ALLOWED
      PARAMETER(MXRF=50000)
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IA(7*MXRF),AA(7*MXRF)
      INTEGER ISORT(5)
      LOGICAL ONL,CCPONL,END,EOF
      CHARACTER*70 TITLE
      CHARACTER*80 HSTR(2)
      CHARACTER*7 LFIN
      DIMENSION LCOL(100),CELL(6),ADATA(100),IHKL(3)
      EQUIVALENCE (IA(1),AA(1))
      DATA LFIN/'HKLIN00'/
      DATA NTWO20/1048576/
      DATA NTWO12/4096/
      DATA NTWO4/16/
      DATA HSTR/'Created by DIFFLMRG program',' '/
      DATA ISORT/1,2,3,0,0/
C
C INITIALISATIONS
C ---------------
C
      IFAIL = 0
      CALL CCPOPN(5,'DATA',5,1,LDUM,IFAIL)
      CALL CCPOPN(6,'PRINTER',6,1,LDUM,IFAIL)
      CALL CCPVRS(6,'DIFFLMRG','15/10/92')
      ONL = CCPONL(IDUM)
      NREF = 0
      NSIGZ = 0
      IB = -7
      IAMAX = 7*MXRF
      CALL MTZINI
C
C READ IN CONTROL CARDS
C ---------------------
C
C DATA CARD 1   NFILES
C
      IF (ONL) WRITE (6,6040)
      READ (5,*) NFILES
      IF (NFILES.LT.2) THEN
C
C ERROR CONDITIONS
C ----------------
C
         WRITE (6,6070)
         STOP
      ELSE
C
C DATA CARD 2    TITLE FOR OUTPUT FILE
C
         IF (ONL) WRITE (6,6050)
         READ (5,6010) TITLE
C
C OPEN FIRST INPUT MTZ FILE AND OPEN OUTPUT MTZ FILE
C --------------------------------------------------
C
         NF = 1
         NN = 0
         CALL LROPEN (1,'HKLIN1',1,IFAIL)
         IF (IFAIL.EQ.1) STOP
         CALL LRNCOL (1,NCOLIN)
         IF (NCOLIN.LT.7 .OR. NCOLIN.GT.9) THEN
            WRITE (6,6120)
            STOP
         ELSE
            CALL LWOPEN(1,'HKLOUT')
            CALL LWHIST(1,HSTR,1)
            CALL LWSORT(1,ISORT)
C
C COPY REFLECTION RECORDS TO ARRAY 'IA/AA'
C
 20         CALL LRREFL(1,RESOL,ADATA,EOF)
            IF (EOF) GO TO 40
            IF (ADATA(5).LE.0.0 .OR. ADATA(7).LE.0.0) THEN
               NSIGZ = NSIGZ + 1
               GO TO 20
            END IF
            IB = IB + 7
            NREF = NREF + 1
            NN = NN + 1
            IF (NREF.GT.MXRF) GO TO 140
            IH = NINT(ADATA(1))
            IK = NINT(ADATA(2))
            IL = NINT(ADATA(3))
            IF (IH.LT.-127 .OR. IH.GT.127) GO TO 150
            IF (IK.LT.-127 .OR. IK.GT.127) GO TO 150
            IF (IL.LT.-127 .OR. IL.GT.127) GO TO 150
            IA(IB+1) = NTWO20*(IH+128) + NTWO12*(IK+128)
     +                 + NTWO4*(IL+128)
            DO 30 I = 4 , NCOLIN
               J = I - 2
               AA(IB+J) = ADATA(I)
 30         CONTINUE
            GO TO 20
         END IF
      END IF
 40   CALL LRCLOS(1)
      WRITE (6,6020) NF , NN
C
C LOOP THROUGH REMAINING INPUT MTZ FILES COPYING REFLECTION RECORDS
C -----------------------------------------------------------------
C
      DO 80 NF = 2 , NFILES
         NN = 0
         LFIN(6:7) = ' '
         IF (NF.LT.10) THEN
            WRITE (LFIN(6:6),'(I1)') NF
         ELSE
            WRITE (LFIN(6:7),'(I2.2)') NF
         END IF
         IUN = 10 + NF
         CALL LROPEN(2,LFIN,1,IFAIL)
         IF (IFAIL.EQ.1) STOP
         CALL LRNCOL(2,NCOL)
         IF (NCOL.NE.NCOLIN) GO TO 130
 50      CALL LRREFL(2,RESOL,ADATA,EOF)
         IF (EOF) GO TO 70
         IF (ADATA(5).LE.0.0 .OR. ADATA(7).LE.0.0) THEN
            NSIGZ = NSIGZ + 1
            GO TO 50
         END IF
         IB = IB + 7
         NREF = NREF + 1
         NN = NN + 1
         IF (NREF.GT.MXRF) GO TO 140
         IH = NINT(ADATA(1))
         IK = NINT(ADATA(2))
         IL = NINT(ADATA(3))
         IF (IH.LT.-127 .OR. IH.GT.127) GO TO 150
         IF (IK.LT.-127 .OR. IK.GT.127) GO TO 150
         IF (IL.LT.-127 .OR. IL.GT.127) GO TO 150
         IA(IB+1) = NTWO20*(IH+128) + NTWO12*(IK+128)
     +                 + NTWO4*(IL+128)
         DO 60 I = 4 , NCOLIN
            J = I - 2
            AA(IB+J) = ADATA(I)
 60      CONTINUE
         GO TO 50
 70      CALL LRCLOS(2)
         WRITE (6,6020) NF , NN
 80   CONTINUE
C
C SORT REFLECTION DATA ON H,K,L
C -----------------------------
C
      CALL REFSRT(IA,NREF,7,IAMAX)
C
C MERGE AND OUTPUT DATA
C ---------------------
C
      END = .FALSE.
      NROUT = 0
      NMULT = 0
      M1 = 1
      ILAST = IA(1)
      K = 0
C
C STEP THROUGH THE REFLECTIONS
C
 90   K = K + 1
      IF (K.GT.NREF) THEN
         END = .TRUE.
         GO TO 100
      END IF
      JJ = 7*K - 6
      IF (IA(JJ).EQ.ILAST) GO TO 90
 100  M2 = K - 1
      MM = 7*M1 - 6
      CALL GETHKL(IA(MM),IHKL,KDUM)
      ADATA(1) = IHKL(1)
      ADATA(2) = IHKL(2)
      ADATA(3) = IHKL(3)
      IF (NCOLIN.GT.7) ADATA(8) = AA(MM+5)
      IF (NCOLIN.GT.8) ADATA(9) = AA(MM+6)
      SIGW1 = 0
      SIGW2 = 0
      SIGWF1 = 0
      SIGWF2 = 0
      IF (M1.EQ.M2) THEN
         ADATA(4) = AA(MM+1)
         ADATA(5) = AA(MM+2)
         ADATA(6) = AA(MM+3)
         ADATA(7) = AA(MM+4)
      ELSE
         NMULT = NMULT + 1
         DO 110 M = M1 , M2
            IP = 7*(M-1)
            W1 = 1.0/(AA(IP+3)*AA(IP+3))
            W2 = 1.0/(AA(IP+5)*AA(IP+5))
            SIGW1 = SIGW1 + W1
            SIGW2 = SIGW2 + W2
            SIGWF1 = SIGWF1 + AA(IP+2)*W1
            SIGWF2 = SIGWF2 + AA(IP+4)*W2
 110     CONTINUE
         F1 = 0
         IF (SIGWF1.GT.0) F1 = SIGWF1/SIGW1
         W1 = SIGW1/(M2-M1+1)
         ADATA(4) = F1
         ADATA(5) = SQRT(1/W1)
         F2 = 0
         IF (SIGWF2.GT.0) F2 = SIGWF2/SIGW2
         W2 = SIGW2/(M2-M1+1)
         ADATA(6) = F2
         ADATA(7) = SQRT(1/W2)
      END IF
      CALL LWREFL(1,ADATA)
      NROUT = NROUT + 1
      IF (END) THEN
         CALL LWCLOS(1,1)
         WRITE (6,6060) NREF , NSIGZ , NMULT , NROUT
         STOP
      ELSE
         M1 = K
         ILAST = IA(JJ)
         GO TO 90
      END IF
  130  WRITE (6,6090)
      STOP
 140  WRITE (6,6100) MXRF
      STOP
 150  WRITE (6,6110) IH , IK , IL
      STOP
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (A)
 6020 FORMAT (/,' NUMBER OF REFLECTIONS READ FROM FILE',I3,' IS',I6)
 6030 FORMAT (//,' OUTPUT REFLECTION FILE COMPLETED')
 6040 FORMAT (/,' Number of input files?')
 6050 FORMAT (/,' Title for output MTZ file?')
 6060 FORMAT (//,' TOTAL  NUMBER  OF  REFLECTIONS  READ  IN  ',I6,/,
     +        ' NUMBER OF REFLECTIONS OMITTED WITH SIG=0  ',I6,/,
     +        ' NO. REFLECTIONS WITH MULTIPLE MEASUREMENTS',I6,/,
     +        ' NO. OF REFLECTIONS WRITTEN TO OUTPUT FILE ',I6,//)
C
 6070 FORMAT (/,' **A MINIMUM OF TWO INPUT FILES MUST BE SPECIFIED**')
 6090 FORMAT (/,' **NO. OF COLUMNS MISMATCH BETWEEN FILES**')
 6100 FORMAT (/,' **MAXIMUM OF',I7,' REFLECTIONS EXCEEDED**')
 6110 FORMAT (/,' **INDICES OUT OF RANGE -127 TO 127, HKL=',3I4,' **')
 6120 FORMAT (/,' **NUMBER OF COLUMNS INVALID FOR A DIFFLAUE O/P FILE**'
     +        )
      END
**==refsrt.f  
C
C SUBROUTINE 'REFSRT'
C ===================
C
C THIS SUBROUTINE IS USED TO SORT THE REFLECTION DATA HELD IN CORE. IT
C IS A MODIFIED VERSION OF THE SUBROUTINE 'AHVSOR' FROM ELEANOR DODSON'S
C VERSION OF THE PROGRAM 'CAD'. THIS VERSION SORTS ON A WORD CONTAINING
C THE PACKED INDICES.
C
      SUBROUTINE REFSRT(IA,NREF,NWRD,IDIM)
C
C PARAMETERS
C
C         IA (I/O)  INTEGER ARRAY CONTAINING THE REFLECTION DATA WITH
C                   NWRD WORDS FOR REFLN 1 FOLLOWED BY NWRD WORDS FOR
C                   EACH REFLECTION UP TO NREF REFLECTIONS  (NWRD <=10)
C       NREF (I)    THE NO. OF REFLECTIONS TO BE SORTED
C       NWRD (I)    NUMBER OF WORDS/REFLECTION (INDICES PACKED IN FIRST
C                   WORD FOR EACH REFLECTION)
C       IDIM (I)    THE NUMBER OF ELEMENTS IN THE ARRAY IA
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      IMPLICIT INTEGER(A-Z)
      INTEGER IA(IDIM),ITEMP,CUT(10)
      INTEGER SAVE(2,16)
      DATA CUTOFF/20/
C
C SORT DATA
C ---------
C
      STACK = 0
      LOW = 1
      HIGH = NREF
 20   CHUNK = HIGH - LOW
      IF (CHUNK.GT.CUTOFF) THEN
         MIDDLE = (LOW+HIGH)/2
         IF (ICMP(IA,LOW,IA,MIDDLE).GT.0) THEN
            J1 = NWRD*(LOW-1)
            J2 = NWRD*(MIDDLE-1)
            DO 30 I = 1 , NWRD
               J1 = J1 + 1
               J2 = J2 + 1
               ITEMP = IA(J1)
               IA(J1) = IA(J2)
               IA(J2) = ITEMP
 30         CONTINUE
         END IF
         IF (ICMP(IA,MIDDLE,IA,HIGH).GT.0) THEN
            J1 = NWRD*(MIDDLE-1)
            J2 = NWRD*(HIGH-1)
            DO 40 I = 1 , NWRD
               J1 = J1 + 1
               J2 = J2 + 1
               ITEMP = IA(J1)
               IA(J1) = IA(J2)
               IA(J2) = ITEMP
 40         CONTINUE
            IF (ICMP(IA,LOW,IA,MIDDLE).GT.0) THEN
               J1 = NWRD*(LOW-1)
               J2 = NWRD*(MIDDLE-1)
               DO 50 I = 1 , NWRD
                  J1 = J1 + 1
                  J2 = J2 + 1
                  ITEMP = IA(J1)
                  IA(J1) = IA(J2)
                  IA(J2) = ITEMP
 50            CONTINUE
            END IF
         END IF
         J = NWRD*(MIDDLE-1)
         DO 60 I = 1 , NWRD
            J = J + 1
            CUT(I) = IA(J)
 60      CONTINUE
         UP = LOW
         DOWN = HIGH
 70      DOWN = DOWN - 1
         IF (ICMP(CUT,1,IA,DOWN).LT.0) GO TO 70
 80      UP = UP + 1
         IF (ICMP(CUT,1,IA,UP).GT.0) GO TO 80
         IF (UP.LT.DOWN) THEN
            J1 = NWRD*(UP-1)
            J2 = NWRD*(DOWN-1)
            DO 90 I = 1 , NWRD
               J1 = J1 + 1
               J2 = J2 + 1
               ITEMP = IA(J1)
               IA(J1) = IA(J2)
               IA(J2) = ITEMP
 90         CONTINUE
            GO TO 70
         ELSE IF (UP.EQ.DOWN) THEN
            DOWN = DOWN - 1
            UP = UP + 1
         END IF
      ELSE
         IF (CHUNK.GE.2) THEN
            MM = HIGH - 1
            DO 120 M = LOW , MM
               MIN = M
               MMM = M + 1
               DO 100 N = MMM , HIGH
                  IF (ICMP(IA,MIN,IA,N).GT.0) THEN
                     MIN = N
                  END IF
 100           CONTINUE
               IF (MIN.NE.M) THEN
                  J1 = NWRD*(M-1)
                  J2 = NWRD*(MIN-1)
                  DO 110 I = 1 , NWRD
                     J1 = J1 + 1
                     J2 = J2 + 1
                     ITEMP = IA(J1)
                     IA(J1) = IA(J2)
                     IA(J2) = ITEMP
 110              CONTINUE
               END IF
 120        CONTINUE
         END IF
         IF (STACK.EQ.0) RETURN
         LOW = SAVE(1,STACK)
         HIGH = SAVE(2,STACK)
         STACK = STACK - 1
         GO TO 20
      END IF
      STACK = STACK + 1
      IF (STACK.GT.16) THEN
C
C ERROR CONDITION
C ----------------
C
         WRITE (6,6010) STACK
         STOP
      ELSE
         BOTTOM = DOWN - LOW
         TOP = HIGH - UP
         IF (BOTTOM.GT.TOP) THEN
            SAVE(1,STACK) = LOW
            SAVE(2,STACK) = DOWN
            LOW = UP
         ELSE
            SAVE(1,STACK) = UP
            SAVE(2,STACK) = HIGH
            HIGH = DOWN
         END IF
         GO TO 20
      END IF
C
C FORMAT STATEMENTS
C -----------------
C
 6010 FORMAT (' **STACK LIMIT OF 16 EXCEEDED IN REFSOR**')
      END
**==icmp.f  
C
C FUNCTION 'ICMP'
C ===============
C
C THIS FUNCTION IS USED TO COMPARE THE VALUES OF THE SORT KEYS
C FOR 2 REFLECTIONS  (ASSUMES 7 WORDS/REFLECTION)
C
      FUNCTION ICMP(IARR1,I1,IARR2,I2)
C
C PARAMETERS
C
C        IARR1 (I)   ARRAY CONTAINING 1ST REFLN
C           I1 (I)   POINTER TO REFLECTION IN IARR1
C        IARR2 (I)   ARRAY CONTAINING 2ND REFLN
C           I2 (I)   POINTER TO REFLECTIONIIN IARR2
C
C FUNCTION RETURN VALUES  -1  REFLN 1 LT REFLN 2
C                          0  REFLN 1 EQ REFLN 2
C                          1  REFLN 1 GT REFLN 2
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      INTEGER IARR1(*),IARR2(*)
C
C COMPARE KEYS
C ------------
C
      J1 = 7*I1 - 6
      J2 = 7*I2 - 6
      IF (IARR1(J1).LT.IARR2(J2)) THEN
         ICMP = -1
         RETURN
      ELSE IF (IARR1(J1).EQ.IARR2(J2)) THEN
         ICMP = 0
         RETURN
      ELSE
         ICMP = 1
         RETURN
      END IF
      END
**==gethkl.f  
C
C SUBROUTINE 'GETHKL'
C ===================
C
C THIS SUBROUTINE UNPACKS THE INDICES FOR A LAUE REFLECTION
C
      SUBROUTINE GETHKL(IDX,IDXL,KFILM)
C
C PARAMETERS
C
C      IDX (I)   WORD CONTAINING PACKED INDICES
C     IDXL (O)   3 WORD ARRAY CONTAINING H, K, L.
C    KFILM (O)   FILM NUMBER (1 TO 6)
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DIMENSION IDXL(3)
      DATA NTWO20/1048576/
      DATA NTWO16/65536/
      DATA NTWO12/4096/
      DATA NTWO4/16/
C
C UNPACK INDICES
C --------------
C
      I = IDX
      IDXL(1) = I/NTWO20
      I = I - NTWO20*IDXL(1)
      IDXL(1) = IDXL(1) - 128
      IDXL(2) = I/NTWO12
      I = I - NTWO12*IDXL(2)
      IDXL(3) = I/NTWO4
      KFILM = I - NTWO4*IDXL(3)
      IDXL(2) = IDXL(2) - 128
      IDXL(3) = IDXL(3) - 128
      RETURN
      END