C ALLHKLPROG  (by G.S.Pawley; well before 1988)
C     OPEN(UNIT=5, FILE='infile', STATUS='old')
C     OPEN(UNIT=6, FILE='outfile', STATUS='new')
C     OPEN(UNIT=8, FILE='plot', STATUS='new')
C  suggestion DEC 1988     WISE could be individually read in
C**************************** MAIN PROGRAM *****************************
C MAXIMUM 6420 LESS 2 BLOCKS = POINTS IN THE SCAN
C MAXIMUM 300 PARAMETERS VARIABLE  THEREFORE AM(300*301/2)
C MAXIMUM 300 PARAMETERS
C MAXIMUM 1000 REFLECTIONS
C MAXIMUM 31 POINTS PER REFLECTION, HALFWIDTH 15
C BLOCKSIZE 2 OF 50
      COMMON AM(45150),V(600),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NPEAKS,NBG,NDX(6,1500),STEP,
     2 WISE(6420,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ISTOP,KEY,NREFL,TWOPIS,KEYREF(600),
     4 P(600),KI(600),PWISE(9),H(3),RHO,NOUT,NPM,MINPTU,MAXPTU,
     5 TWOPI,YC,WAVETC,DC(600),YSUM,YMASK,YAV,NMAX,NEXT,NSIDE,IW,
     6 KEYSUM,MSTEPU,SNTL(600),LVAR(600),NO,KOUT,
     7 SLACK(10),SIG,REFSIG,LOT,K15
      DIMENSION DWISE(100,300),DIFFAC(6),HJ(3),HK(3),SINS(600),
     1 CLOR(600),INDIC(600),INDT(5)
   17 FORMAT(1X,3I3,I4,I5,I2,I5,26X,F13.5,49X,F12.4)
   18 FORMAT(20X,I6,26X,F13.5)
   19 FORMAT(1X,3I3,I4,I5,I2,I5,4F13.5,4F9.2,F12.4)
   20 FORMAT(' NUMBER OF VARIABLES THIS CYCLE IS',I4,I8,' REFLECTIONS')
   21 FORMAT(' ',47X,'NUMERATOR',9X,'DENOMINATOR',5X,'R')
   22 FORMAT(' R FACTOR INCLUDING ZEROS',17X,F15.3,4X,F15.3,4X,F5.3)
   24 FORMAT(10F12.1)
   25 FORMAT(3I8,2F8.4,I8,8A4)
   26 FORMAT(' AGREEMENT FACTORS BASED ON PARAMETERS BEFORE CYCLE'I3/' S
     1UM(W*(O-C)**2) IS'E11.3/' SQRT(SUM(W*(0-C)**2)/(NO-NV)) IS'F10.4)
   27 FORMAT(' ESTIMATED FACTORS BASED ON PARAMETERS AFTER  CYCLE'I3/' S
     1UM(W*(O-C)**2) IS'E11.3/' SQRT(SUM(W*(0-C)**2)/(NO-NV)) IS'F10.4)
   28 FORMAT(' MATRIX HAS A ZERO DIAGONAL ELEMENT CORRESPONDING TO PARAM
     1ETER'I3,' OF THOSE VARIED')
   29 FORMAT(//31(' + +')/' PARAMETERS AFTER LEAST SQUARES CYCLE',I3,/' 
     1 INDICES',12X,'PARAM',9X,'OLD',8X,'CHANGE',9X,'NEW',9X,'ERROR   RA
     2TIO NEW ERROR  2*THETA  LORENTZ   INTENSITY')
   30 FORMAT(20X,I6,4F13.5,2F9.2)
   31 FORMAT(' CORRELATION MATRIX')
   32 FORMAT(' 'I3,10F9.4/(' '3X,10F9.4))
   33 FORMAT(' REFLECTIONS REORDERED',I4)
   34 FORMAT(' REDUCED SCAN VALUES ARE',2I7,' USING A STEP OF',F7.4)
   35 FORMAT(' MIN. AND MAX. 2-THETA TO BE USED',2F8.2,/' WAVELENGTH'
     1,F12.4)
   36 FORMAT(' LIMITING SPHERE REACHED',I6,' REFLECTIONS USED')
   37 FORMAT(' MEAN F.W.H.H. ='I4/' EQUIVALENT NUMBER OF REFLECTIONS ='I
     14/' NEW SQSIG ='F10.4)
   38 FORMAT(' MINIMUM AND MAXIMUM POINTS IN SCAN TO BE USED ARE',2I7)
   39 FORMAT(5(' THERE ARE PROBABLY TOO MANY VARIABLES BEING REFINED'/))
      OPEN(UNIT=5, FILE='infile', STATUS='old')
      OPEN(UNIT=6, FILE='outfile', STATUS='new')
      OPEN(UNIT=8, FILE='plot', STATUS='new')
      K15=15
      NVM=300
      NTOT=(NVM*(NVM+1))/2
      NPM=600
      NATOM=40
      NREFS=1500
      NWIDE=15
      KBL=50
      KBLKBL=2*KBL
      NMAX=6420-KBLKBL
      KTEST=KBLKBL-NWIDE
      PI=3.1415926
      TWOPIS=2.0*PI**2
      TWOPI=2.0*PI
      DEGRAD=PI/180.
      GAUSS=2.7726
      CONMEW=1.253314137
      YMASK=-1111.
      INDIC(1)=1
      DO 1001 I = 2,NPM
      K=2
      IF(I.GT.5)  K=3
      IF(I.GT.9)  K=2
      IF(I.GT.K15) K=6
 1001 INDIC(I)=K
      INDIC(9)=4
      READ(5,25)  MINPT,MAXPT,NPEAKS,STEP,WAVETC,NEXT,(P(I),I=1,8)
      WRITE(6,25) MINPT,MAXPT,NPEAKS,STEP,WAVETC,NEXT,(P(I),I=1,8)
      READ(5,25)  MINPTU,MAXPTU,MSTEPU
      K=MINPTU/MSTEPU
      K=K*MSTEPU
      IF(K.NE.MINPTU) MINPTU=MSTEPU+K
      K=MAXPTU/MSTEPU
      MAXPTU=K*MSTEPU
      WRITE(6,38) MINPTU,MAXPTU
      MINPTU=MINPTU/MSTEPU
      MAXPTU=MAXPTU/MSTEPU
      STEP=STEP*FLOAT(MSTEPU)
      WRITE(6,34) MINPTU,MAXPTU,STEP
      PTMIN=FLOAT(MINPTU)*STEP
      PTMAX=FLOAT(MAXPTU)*STEP
C     THE LIMITS ARE SET 10 PERCENT TOO WIDE TO ALLOW FOR REFINEMENT
      MINRH= 10000.*(SIN(PTMIN*DEGRAD*0.5)/WAVETC)**2
      RMAX=         SIN(PTMAX*DEGRAD*0.5)/WAVETC
      MAXRH=10000.*RMAX**2
      CALL PRELIM
      WRITE(6,35) PTMIN,PTMAX,WAVETC
      IF(NREFS.EQ.-1) CALL FINISH(6)
      IF(NC.EQ.0) NC=1
      DO 1051 IC = 1,NC
      DO 1002 I=1,NTOT
 1002 AM(I)=0.0
      DO 1004 I=1,NVM
      DO 1003 J = 1,KBLKBL
 1003 DWISE(J,I)=0.0
 1004 V(I)=0.0
      DO 1005 I = 1,NMAX
 1005 WISE(I,2)=0.0
      MEANWD=0
      NBLOCK=0
      BLOCK=0.0
      RZNUM=0.0
      SIG=0.0
      IQUERY=0
      ISWAP=0
C     REORDER REFLECTIONS USING CURRENT CELL
      DO 1011 NREFL = 1,NPEAKS
      HK(1)=NDX(1,NREFL)
      HK(2)=NDX(2,NREFL)
      HK(3)=NDX(3,NREFL)
      HJ(1)=ZT(1,1)*HK(1)+ZT(2,1)*HK(2)+ZT(3,1)*HK(3)
      HJ(2)=ZT(1,2)*HK(1)+ZT(2,2)*HK(2)+ZT(3,2)*HK(3)
      HJ(3)=ZT(1,3)*HK(1)+ZT(2,3)*HK(2)+ZT(3,3)*HK(3)
      S=0.25*(HJ(1)**2+HJ(2)**2+HJ(3)**2)
      SNTL(NREFL)=SQRT(S)
      NDX(5,NREFL)=2.0*ASIN(SNTL(NREFL)*WAVETC)/(DEGRAD*STEP)
      IF(LOT.EQ.0) GO TO 1011
      IF(NREFL.EQ.1) GO TO 1011
      INDX=NDX(5,NREFL)
      K=NREFL
      J=0
      DO 1006 I = 1,20
      IF(K.EQ.1) GO TO 1007
      K=K-1
      IF(INDX.GE.NDX(5,K)) GO TO 1007
 1006 J=J+1
 1007 IF(J.EQ.0) GO TO 1011
      ISWAP=ISWAP+1
      S=SNTL(NREFL)
      T=P(NREFL+K15)
      DO 1008 I = 1,5
 1008 INDT(I)=NDX(I,NREFL)
      DO 1009 I = 1,J
      K=NREFL-I
      L=K+1
      SNTL(L)=SNTL(K)
      P(L+K15)=P(K+K15)
      DO 1009 M = 1,5
 1009 NDX(M,L)=NDX(M,K)
      SNTL(K)=S
      P(K+K15)=T
      DO 1010 I = 1,5
 1010 NDX(I,K)=INDT(I)
 1011 CONTINUE
      IF(ISWAP.NE.0) WRITE(6,33) ISWAP
C     SET UP FOR VARYING THE LOT
      IF(LOT.EQ.0) GO TO 1020
      M=KEYSUM
      DO 1016 NREFL = 1,NPEAKS
      N=NREFL+K15
      NDX(4,NREFL)=1
      KI(N)=1
      IF(NREFL.EQ.1) GO TO 1015
      I=NREFL-1
 1012 IF(NDX(4,I).NE.0) GO TO 1013
      I=I-1
      GO TO 1012
 1013 IF(NDX(5,I).NE.NDX(5,NREFL)) GO TO 1015
      L=NDX(4,I)+1
      NDX(4,I)=L
      NDX(4,NREFL)=0
      S=(FLOAT(L-1)*P(I+K15)+P(N))/FLOAT(L)
      DO 1014 J = 1,L
      K=N+1-J
 1014 P(K)=S
      KI(N)=0
 1015 M=M+KI(N)
      KEYREF(NREFL)=M
 1016 CONTINUE
      NV=M
      NPM=NP-1
      MS=NV
      MV=1
      MZ=1
      DO 1019 I = 1,NPM
      IF(I.LE.KEYSUM) GO TO 1018
      IF(I.LE.K15) GO TO 1019
      IF(KI(I).EQ.0) GO TO 1019
      NREFL=I-K15
      N=NDX(5,NREFL)
      NREFLP=NREFL+1
      MZT=MZ
      MST=MS
      K=0
      DO 1017 M = NREFLP,NPEAKS
      IF(NDX(4,M).EQ.0) GO TO 1017
      K=K+1
      L=NDX(5,M)-N
      IF(L.GT.NSIDE) GO TO 1018
      MW=MV+K
      MZT=MZT+MST
      MST=MST-1
      MZP=MZ+K
      AM(MZ )=AM(MZ )+SLACK(L)
      AM(MZP)=AM(MZP)-SLACK(L)
      AM(MZT)=AM(MZT)+SLACK(L)
      DIFF=P(I)*FLOAT(NDX(4,NREFL))-P(I+K)*FLOAT(NDX(4,M))
      V(MV)=V(MV)-DIFF*SLACK(L)
      V(MW)=V(MW)+DIFF*SLACK(L)
 1017 CONTINUE
 1018 MV=MV+1
      MZ=MZ+MS
      MS=MS-1
 1019 CONTINUE
 1020 NM=(NV*(NV+1))/2
      M=NV-KEYSUM
      WRITE(6,20) NV,M
      IF(NV.EQ.0) CALL FINISH(1)
C     START LOOP THROUGH OBSERVATIONS
      DO 1039 NREFL = 1,NPEAKS
      H(1)=NDX(1,NREFL)
      H(2)=NDX(2,NREFL)
      H(3)=NDX(3,NREFL)
      SINTHL=SNTL(NREFL)
      SINTH=SINTHL*WAVETC
      POINT=2.0*ASIN(SINTH)/DEGRAD
      COSTH=SQRT(1.0-SINTH**2)
      TANTH=SINTH/COSTH
      UVW=(PWISE(6)*TANTH+PWISE(7))*TANTH+PWISE(8)
      FACTOR=WAVETC/(SINTHL*COSTH*DEGRAD)
      AREA=SQRT(UVW)
C     LORENTZ CORRECTION, POLARISATION MUST BE ADDED FOR X-RAYS
      CLOR(NREFL)=1./(COSTH*SINTH**2)
      SINS(NREFL)=POINT
      YC=P(NREFL+K15)
 1021 PHICAL=POINT-PWISE(5)-BLOCK
      MIDREF=PHICAL/STEP+0.5
      IF(MIDREF.GT.KTEST) GO TO 1026
      IQUERY=NREFL
      NWIDTH=2.0*AREA/STEP
      MEANWD=MEANWD+NWIDTH
      IF(NWIDTH.GT.NWIDE) NWIDTH=NWIDE
      NSTART=MIDREF-NWIDTH
      NFINIS=MIDREF+NWIDTH
      IF(NSTART.LT.1) CALL FINISH(2)
      DC(5)=-1.0
      DC(6)= TANTH**2
      DC(7)= TANTH
      DC(8)= 1.0
      COT2TH=(COSTH-0.5/COSTH)/SINTH
      GAMMA3=PWISE(9)*COT2TH
      GAMMA2=0.25*GAMMA3
      GAURAT=GAUSS/UVW
      DC(9)=-0.5*GAURAT*COT2TH
      YCETC=2.0*GAURAT*YC
      YUVW=YC/UVW
      RATIND=1.0/AREA
      RATIN4=4.0*RATIND
      DO 1022 L = 1,6
 1022 DC(L+9)=(DCELL(L,1)*H(1)**2+DCELL(L,4)*H(2)*H(3)
     1        +DCELL(L,2)*H(2)**2+DCELL(L,5)*H(3)*H(1)
     2        +DCELL(L,3)*H(3)**2+DCELL(L,6)*H(1)*H(2))*FACTOR
      IF(LATTIS.LT.8) GO TO 1023
      DC(10)=DC(10)+DC(11)
      IF(LATTIS.LT.12) GO TO 1023
      IF(LATTIS.GT.19) GO TO 1023
      DC(10)=DC(10)+DC(12)
      DC(13)=DC(13)+DC(14)+DC(15)
 1023 DO 1025 NSTEP = NSTART,NFINIS
      IBLOCK=NSTEP+NBLOCK
      DELTA1=NSTEP*STEP-PHICAL
      DELTA2=DELTA1+GAMMA2
      DELTA3=DELTA1+GAMMA3
      EXPON1=GAURAT*DELTA1**2
      EXPON2=GAURAT*DELTA2**2
      EXPON3=GAURAT*DELTA3**2
      TERM1=RATIND*EXP(-EXPON1)
      TERM2=RATIN4*EXP(-EXPON2)
      TERM3=RATIND*EXP(-EXPON3)
      YCPART=TERM1+TERM2+TERM3
      DIFFAC(2)=YCETC*(DELTA1*TERM1+DELTA2*TERM2+DELTA3*TERM3)
      DIFFAC(3)=YUVW*(TERM1*EXPON1+TERM2*EXPON2+TERM3*EXPON3-0.5*YCPART)
      DIFFAC(4)=YC*(TERM2*DELTA2+4.0*TERM3*DELTA3)
      DIFFAC(6)=YCPART
      YCPART=YCPART*YC
      DIFFAC(1)=YCPART
      IF(KEY.GT.NV) GO TO 1025
      IF(KEY.GT.KEYSUM) GO TO 2024
      DO 1024 L = KEY,KEYSUM
      M=LVAR(L)
      K=INDIC(M)
 1024 DWISE(NSTEP,L)=DWISE(NSTEP,L)+DIFFAC(K)*DC(M)
C     INDIC(K)=6 FOR ALL ABOVE P(K15), WHERE PEAKS START
C     PEAK HEIGHT DIFFERENTIALS ARE ALL PWISE(1)
 2024 IF(NV.EQ.KEYSUM) GO TO 1025
      L=NREFL+K15
      IF(KI(L).EQ.0) GO TO 1025
      M=KEYREF(NREFL)
      DWISE(NSTEP,M)=DWISE(NSTEP,M)+DIFFAC(6)
 1025 WISE(IBLOCK,2)=WISE(IBLOCK,2)+YCPART
      IF(IQUERY.LT.NPEAKS) GO TO 1039
C     ADD DIFFERENTIALS FOR BACKGROUND AND SCALE IF REQUIRED AND
C     ADD BLOCK TO LEAST SQUARES MATRIX, SUM DIFFERENCES
 1026 DO 1037 KBLOCK = 1,KBL
      IBLOCK=KBLOCK+NBLOCK
      BACKGD=PWISE(2)+PWISE(3)*IBLOCK
      L=0
      IF(KI(1).EQ.0) GO TO 1027
      L=L+1
      DWISE(KBLOCK,L)=WISE(IBLOCK,2)/PWISE(1)
 1027 IF(KI(2).EQ.0) GO TO 1028
      L=L+1
      DWISE(KBLOCK,L)=1.0
 1028 IF(KI(3).EQ.0) GO TO 2028
      L=L+1
      DWISE(KBLOCK,L)=IBLOCK
 2028 WISE(IBLOCK,2)=WISE(IBLOCK,2)+BACKGD
      IF(WISE(IBLOCK,1).LE.YMASK) GO TO 1037
      DY=WISE(IBLOCK,1)-WISE(IBLOCK,2)
      ABSDY=ABS (DY)
      RZNUM=RZNUM+ABSDY
      GO TO (1029,1030,1031,1032),IW
 1029 WEIGHT=1.0
      GO TO 1033
 1030 WEIGHT=1.0/WISE(IBLOCK,1)
      GO TO 1033
 1031 WEIGHT=1.0/(WISE(IBLOCK,1)**2)
      GO TO 1033
 1032 WEIGHT=ABSDY
      GO TO 1033
C     WEIGHT IS THE SQUARE OF ORFLS WEIGHT
C     WDY HAS BEEN REMOVED
 1033 SIG=SIG+WEIGHT*DY**2
      JK=1
      DO 1036 J=1,NV
      TEMP=DWISE(KBLOCK,J)*WEIGHT
      IF(TEMP.NE.0.0)  GO TO 1034
C     BY-PASS IF DERIVATIVE IS ZERO
      JK=JK+NV+1-J
      GO TO 1036
 1034 DO 1035 K=J,NV
      AM(JK)=AM(JK)+TEMP*DWISE(KBLOCK,K)
      JK=JK+1
 1035 CONTINUE
      V(J)=V(J)+TEMP*DY
 1036 JK=JK
 1037 CONTINUE
C     END LOOP TO STORE MATRIX AND VECTOR
C     TRANSFER BLOCK
      DO 1038 L = 1,KBL
      N=L+KBL
      DO 1038 M = 1,NV
      DWISE(L,M)=DWISE(N,M)
 1038 DWISE(N,M)=0.0
      NBLOCK=NBLOCK+KBL
      IF(NBLOCK.GT.NMAX) GO TO 1040
      BLOCK=NBLOCK*STEP
      IF(IQUERY.LT.NPEAKS) GO TO 1021
      IF(IQUERY.GT.NPEAKS) GO TO 1039
      IQUERY=IQUERY+1
      GO TO 1026
 1039 CONTINUE
C     END LOOP THROUGH OBSERVATIONS
 1040 IF(IC.EQ.NOUT) WRITE(6,24) (WISE(I,1),WISE(I,2),I=MINPTU,MAXPTU)
      IF(IC.EQ.NOUT)  CALL PLOTLP
      IF(IC.EQ.-NOUT) CALL PLOTLP
C     COMPUTE AND PUT OUT AGREEMENT FACTORS
      SQSIG=SQRT (ABS (SIG/FLOAT (NO-NV)))
      MEANWD=MEANWD/NPEAKS
      NEWNO=NO/MEANWD
      SQSNY=0.0
      IF(NEWNO.LE.NV) GO TO 2040
      SQSNY=SQRT (ABS (SIG/FLOAT (NEWNO-NV)))
 2040 WRITE(6,37) MEANWD,NEWNO,SQSNY
      IF(SQSNY.EQ.0.0) WRITE(6,39)
      SQSNY=SQSNY/SQSIG
      WRITE(6,26) IC,SIG,SQSIG
      CALL SETWT(SLACK,SIG,REFSIG,NO)
      RZ=RZNUM/YSUM
      WRITE(6,21)
      WRITE(6,22) RZNUM,YSUM,RZ
C     START LOOP TO TEST FOR ZERO DIAGONAL ELEMENT
      ISING=0
      II=1
      IID=NV
      DO 1042 I=1,NV
      IF(AM(II).NE.0.0)  GO TO 1041
      ISING=1
      WRITE(6,28)  I
 1041 II=II+IID
      IID=IID-1
 1042 CONTINUE
C     END LOOP TO TEST FOR ZERO DIAGONAL ELEMENT
C     TERMINATE JOB IF ZERO DIAGONAL ELEMENT WAS FOUND
      IF(ISING.NE.0)  CALL FINISH(11)
C     ENTER SUBROUTINE TO REPLACE MATRIX WITH INVERSE
      CALL MATINV(NV,ISING)
      IF(ISING.NE.0) CALL FINISH(10)
C     START LOOP FOR MATRIX VECTOR MULTIPLICATION FOR PARAMETER CHANGES
      DO 1045 I=1,NV
      PDI=0.0
      IJ=I
      IJD=NV-1
      DO 1044 J=1,NV
      PDI=PDI+AM(IJ)*V(J)
      IF(J.GE.I)  GO TO 1043
      IJ=IJ+IJD
      IJD=IJD-1
      GO TO 1044
 1043 IJ=IJ+1
 1044 CONTINUE
      DC(I)=PDI
      SIG=SIG-PDI*V(I)
 1045 CONTINUE
C     END LOOP FOR MATRIX VECTOR MULTIPLICATION
C     RECOMPUTE AGREEMENT FACTOR USING MODIFIED SIG
      SQSIG=SQRT (ABS (SIG/FLOAT (NO-NV)))
C     PUT OUT CAPTION FOR LIST OF CORRECTED PARAMETERS
      WRITE(6,29) IC
C     START LOOP TO CORRECT AND PUT OUT PARAMETERS
      J=1
      II=1
      IJ=NV
      L=1
      DO 1049 I=1,NP
      NREFL=I-K15
      IF(KI(I).EQ.0) GO TO 1048
      IF(I.GT.K15) L=NDX(4,NREFL)
      FL=L
      POLD=P(I)
      PDI=DC(J)
      P(I)=POLD+PDI/FL
      DC(J)=AM(II)
      II=II+IJ
      IJ=IJ-1
      SIGP=SQRT(DC(J))*SQSIG
      SIGNEW=SIGP*SQSNY
      SIGRAT=PDI/SIGP
      J=J+1
      IF(I.LE.K15) GO TO 1047
C     DISTRIBUTE THE INTENSITY BETWEEN CLOSE REFLECTIONS
      M=I+L-1
      DO 1046 N = I,M
 1046 P(N)=P(I)
      MULT=NDX(6,NREFL)
      AMULT=MULT
      FININT=P(I)/(AMULT*CLOR(NREFL))
      WRITE(6,19) (NDX(K,NREFL),K=1,5),MULT,I,POLD,PDI,P(I),SIGP,
     1               SIGRAT,SIGNEW,SINS(NREFL),CLOR(NREFL),FININT
      GO TO 1049
 1047 WRITE(6,30) I,POLD,PDI,P(I),SIGP,SIGRAT,SIGNEW
      GO TO 1049
 1048 IF(I.LE.K15) WRITE(6,18) I,P(I)
      IF(I.LE.K15) GO TO 1049
      FINNEW=AMULT*FININT/FLOAT(NDX(6,NREFL))
      WRITE(6,17) (NDX(K,NREFL),K=1,6),I,P(I),FINNEW
 1049 CONTINUE
      IF(LATTIS.LT.8) GO TO 1050
      P(11)=P(10)
      IF(LATTIS.LT.12) GO TO 1050
      IF(LATTIS.GT.19) GO TO 1050
      P(12)=P(10)
      P(14)=P(13)
      P(15)=P(13)
C     END LOOP TO CORRECT AND PUT OUT PARAMETERS
C     PUT OUT ESTIMATED AGREEMENT FACTORS
 1050 WRITE(6,27) IC,SIG,SQSIG
C     ENTER SUBROUTINE TO TEST PARAMETERS
      ISTOP=0
      CALL TEST
      IF(ISTOP.NE.0) CALL FINISH(9)
 1051 CONTINUE
C     END LOOP THROUGH NC CYCLES AND FINAL CALC OF Y
C     TERMINATE JOB
      IF(NC.EQ.0)  GO TO 1056
C     CALCULATE AND PUT OUT CORRELATION MATRIX
C*****      WRITE(6,31)
      DO 1052 I=1,NV
      DC(I)=1.0/SQRT(DC(I))
 1052 CONTINUE
      IJ=1
      DO 1055 I=1,NV
      IF(I.NE.1)  V(I-1)=10000.0
      DO 1053 J=I,NV
      V(J)=0.0
 1053 CONTINUE
      DO 1054 J=I,NV
      V(J)=AM(IJ)*DC(I)*DC(J)
      IJ=IJ+1
 1054 CONTINUE
      K=1+10*((I-1)/10)
C*****      WRITE(6,32) I,(V(J),J=K,NV)
 1055 IJ=IJ
 1056 CALL FINISH(12)
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SETWT(SLACK,SIG,REFSIG,NO)
      DIMENSION SLACK(10)
   19 FORMAT('   NEW SLACK ARE NOW    ',4F12.2/6F12.2)
      IF(SIG.EQ.0) GO TO 1092
      W=SIG/(FLOAT(NO)*REFSIG**2)
      DO 1091 I = 1,10
      SLACK(I)=W
 1091 W=W*0.5
      GO TO 1094
 1092 DO 1093 I = 1,10
 1093 SLACK(I)=0.0
 1094 WRITE(6,19) SLACK
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE PRELIM
      COMMON AM(45150),V(600),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NPEAKS,NBG,NDX(6,1500),STEP,
     2 WISE(6420,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ISTOP,KEY,NREFL,TWOPIS,KEYREF(600),
     4 P(600),KI(600),PWISE(9),H(3),RHO,NOUT,NPM,MINPTU,MAXPTU,
     5 TWOPI,YC,WAVETC,DC(600),YSUM,YMASK,YAV,NMAX,NEXT,NSIDE,IW,
     6 KEYSUM,MSTEPU,SNTL(600),LVAR(600),NO,KOUT,
     7 SLACK(10),SIG,REFSIG,LOT,K15
      DIMENSION HJ(3),HK(3)
   40 FORMAT(' VALUE OF SIG USED FOR INITIAL SLACK IS',F12.0/' REFLECTIO
     1N DIFFERENCE FOR SLACK IS',F12.4,/' SUBORDINATION RANGE IS',I4/)
   41 FORMAT(6F12.5)
   42 FORMAT(24I3)
   43 FORMAT(8F9.3)
   44 FORMAT(3(F11.6,2I2))
   45 FORMAT(//I10,' = NUMBER OF REFLECTIONS'/I10,' = NUMBER OF PARAMETE
     1RS READ')
   46 FORMAT(10(2X,F6.0))
   47 FORMAT(3I3,F12.2)
   48 FORMAT(I10,' = NUMBER OF OBSERVED POINTS'/F10.0,' = TOTAL COUNTS'/
     1F10.3,' = AVERAGE COUNT'/E11.3,' STATISTICAL SIG'/)
   49 FORMAT(' REFLECTION INTENSITIES TO BE VARIED'/(I4,2I3,F12.2))
   50 FORMAT(' UNIFORM BACKGROUND TO BE SUBTRACTED',F12.6)
      READ(5,42)  NC,NOUT,KOUT,IW,IRHOMB,LOT,NSIDE
      WRITE(6,42) NC,NOUT,KOUT,IW,IRHOMB,LOT,NSIDE
      READ(5,41)  SIG,REFSIG
      NP=NPEAKS+K15
      MPE=NPEAKS
C     SET UP TRANSFORMATION TO ORTHOGONAL AXES
      READ(5,41)  (RCP(I),I=1,6)
      WRITE(6,41) (RCP(I),I=1,6)
      CALL CELLAZ
      READ(5,41)  (PWISE(I),I=1,9),BGO
      WRITE(6,41) (PWISE(I),I=1,9)
C     WRITE OUT BACKGROUND TO BE SUBTRACTED
      IF(BGO.NE.0.0) WRITE(6,50) BGO
C     INPUT OF REFLECTIONS AND SCAN
      READ(5,42) LATTIS,(ISCREW(I),I=1,3),(IGLIDE(I),I=1,3)
      IF(MINPT.EQ.0) CALL FINISH(13)
C     REFGEN DOES NOT ALTER THE NUMBER OF REFLECTIONS FITTED
      CALL REFGEN
      IF(NPEAKS.EQ.0) GO TO 2111
      READ(5,47)  ((NDX(I,J),I=1,3),P(J+K15),J=1,NPEAKS)
      WRITE(6,49) ((NDX(I,J),I=1,3),P(J+K15),J=1,NPEAKS)
      GO TO 2112
 2111 NPEAKS=NBG
 2112 NP=K15+NPEAKS
      DO 1111 J = 1,NPEAKS
      IF(MPE.EQ.0) GO TO 1111
      P(J+K15)=P(J+K15)*PWISE(1)
 1111 NDX(4,J)=1
 3111 PWISE(1)=1.0
      WRITE(6,45) NPEAKS,NP
      DO 1112 I = 1,NMAX
 1112 WISE(I,1)=YMASK
      READ(5,46) (WISE(I,1),I=MINPT,MAXPT)
      DO 1113 I = MINPTU,MAXPTU
      J=I*MSTEPU
      IF(WISE(J,1).LT.0.0) WISE(J,1)=YMASK
 1113 WISE(I,1)=WISE(J,1)-BGO
      J=MAXPTU+1
      DO 1114 I = J,MAXPT
 1114 WISE(I,1)=YMASK
      YSUM=0.0
      YSIG=0.0
      NO=0
      DO 1115 I = MINPTU,MAXPTU
      IF(WISE(I,1).LT.YMASK) GO TO 1115
      NO=NO+1
      YSUM=YSUM+WISE(I,1)
      YSIG=YSIG+SQRT(ABS(WISE(I,1)))
 1115 CONTINUE
C     SET UP TRIAL INTENSITIES IF REQUIRED. 
      IF(MPE.NE.0) GO TO 2115
      DO 8011 NREFL = 1,NPEAKS
      HK(1)=NDX(1,NREFL)
      HK(2)=NDX(2,NREFL)
      HK(3)=NDX(3,NREFL)
      HJ(1)=ZT(1,1)*HK(1)+ZT(2,1)*HK(2)+ZT(3,1)*HK(3)
      HJ(2)=ZT(1,2)*HK(1)+ZT(2,2)*HK(2)+ZT(3,2)*HK(3)
      HJ(3)=ZT(1,3)*HK(1)+ZT(2,3)*HK(2)+ZT(3,3)*HK(3)
      S=0.25*(HJ(1)**2+HJ(2)**2+HJ(3)**2)
      SINTH=SQRT(S)*WAVETC
      N=2.0*ASIN(SINTH)/(DEGRAD*STEP)
      COSTH=SQRT(1.0-SINTH**2)
      CORL=COSTH*SINTH**2
      P(K15+NREFL)=(WISE(N,1)-PWISE(2))*CORL
 8011 CONTINUE
 2115 YAV=YSUM/NO
      WRITE(6,48) NO,YSUM,YAV,YSIG
      IF(MAXPTU.GT.NMAX) CALL FINISH(7)
      IF(SIG.EQ.0.0) SIG=YSIG
      WRITE(6,40) SIG,REFSIG,NSIDE
      CALL SETWT(SLACK,SIG,REFSIG,NO)
      IC=0
      CALL TEST
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE REFGEN
      COMMON AM(45150),V(600),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NPEAKS,NBG,NDX(6,1500),STEP,
     2 WISE(6420,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ISTOP,KEY,NREFL,TWOPIS,KEYREF(600),
     4 P(600),KI(600),PWISE(9),H(3),RHO,NOUT,NPM,MINPTU,MAXPTU,
     5 TWOPI,YC,WAVETC,DC(600),YSUM,YMASK,YAV,NMAX,NEXT,NSIDE,IW,
     6 KEYSUM,MSTEPU,SNTL(600),LVAR(600),NO,KOUT,
     7 SLACK(10),SIG,REFSIG,LOT,K15
C SEE PAGE 32 OF INTERNATIONAL TABLES, VOLUME 1
C  LATT = 0  1  2  3  4  5            PROGRAM DEDUCES 'LATT'
C         P  R  C  F  I  H  LATTICE
C VALUE FOR R AND H LATTICES DEPENDS ON POINT GROUP
C LAUE GROUP    LATTIS NUMBERS
C               P    C    F    I    LATTICES
C -1            1
C 2/M  Y UNIQUE 2    3
C MMM           4    5    6    7
C 4/M           8              9
C 4/MMM         10             11
C M3            12        13   14
C M3M           15        16   17
C  RHOMB 18 R-3  R3  R-3C
C  RHOMB 19 R32  R3M  R3C  R-3M
C -3    HEX 20
C -3M   HEX 21
C 6/M   HEX 22
C 6/MMM HEX 23
C SCREW AND GLIDE ABSENCES REMOVED, ORTHORHOMBIC AND MONOCLINIC
C BY SETTING THE MULTIPLICITY TO ZERO
C ISCREW(I) =0 : NO SCREW AXIS
C ISCREW(I) =I : 2 FOLD SCREW AXIS ALONG AXIS I
C IGLIDE(I) =0 : NO GLIDE PLANE PERPENDICULAR TO AXIS I
C  "        =1 : A GLIDE PERPENDICULAR TO AXIS I
C  "        =2 : B GLIDE PERPENDICULAR TO AXIS I
C  "        =3 : C GLIDE PERPENDICULAR TO AXIS I
C  "        =4 : N GLIDE PERPENDICULAR TO AXIS I
   51 FORMAT(//I5,' REFLECTIONS GENERATED IN REFGEN'/' LATTICE TYPE',I4)
   52 FORMAT(' 2-FOLD SCREW AXIS ALONG AXIS',I2)
   53 FORMAT(' GLIDE PLANE PERPENDICULAR TO AXIS',I2,' ALONG AXIS',I2)
   54 FORMAT(' N GLIDE PERPENDICULAR TO AXIS',I2)
   55 FORMAT(' MAXIMUM HKL = ',I4)
   56 FORMAT(3I3,'     0.0    ',2I6)
      DO 1201 I=1,NREFS
      DO 1201 J=1,6
 1201 NDX(J,I)=0
      RCPMAX=0.0
      DO 1203 I = 1,3
      IF(ISCREW(I).EQ.I) WRITE(6,52) I
      DO 1202 J = 1,3
      IF(IGLIDE(I).EQ.J) WRITE(6,53) I,J
 1202 CONTINUE
      IF(IGLIDE(I).EQ.4) WRITE(6,54) I
      IF(RCP(I).GT.RCPMAX) RCPMAX=RCP(I)
 1203 CONTINUE
      LIMHKL=2.0*(RCPMAX*RMAX+1.5)
      WRITE(6,55) LIMHKL
      LATT=0
      LAUE=0
      NBG=1
      IF(LATTIS.GT.11) GO TO 1215
      IF(LATTIS.LT.8 ) LATT=LATTIS-3
      IF(LATTIS.LT.5 ) LATT=0
      IF(LATTIS.EQ.3 ) LATT=2
      IF(LATTIS.EQ.9 ) LATT=4
      IF(LATTIS.EQ.11) LATT=4
      DO 1209 I = 1,LIMHKL
      CALL SETHKL(LATT,I,0,0,2)
      CALL SETHKL(LATT,0,I,0,2)
      CALL SETHKL(LATT,0,0,I,2)
      DO 1209 J = 1,LIMHKL
      IF(LATTIS.EQ.1) GO TO 1204
      CALL SETHKL(LATT,I,J,0,4)
      CALL SETHKL(LATT,0,I,J,4)
      IF(LATTIS.GT.3) GO TO 1205
 1204 CALL SETHKL(LATT,I,0, J,2)
      CALL SETHKL(LATT,I,0,-J,2)
      IF(LATTIS.NE.1) GO TO 1206
      CALL SETHKL(LATT,I, J, 0,2)
      CALL SETHKL(LATT,I,-J, 0,2)
      CALL SETHKL(LATT,0, I, J,2)
      CALL SETHKL(LATT,0, I,-J,2)
      GO TO 1206
 1205 CALL SETHKL(LATT,I,0,J,4)
 1206 DO 1209 K = 1,LIMHKL
      IF(LATTIS.LT.4) GO TO 1207
      CALL SETHKL(LATT,I,J,K,8)
      GO TO 1209
 1207 IF(LATTIS.EQ.1) GO TO 1208
      CALL SETHKL(LATT,I,J, K,4)
      CALL SETHKL(LATT,I,J,-K,4)
      GO TO 1209
 1208 CALL SETHKL(LATT,I, J, K,2)
      CALL SETHKL(LATT,I, J,-K,2)
      CALL SETHKL(LATT,I,-J, K,2)
      CALL SETHKL(LATT,I,-J,-K,2)
 1209 CONTINUE
      IF(LATTIS.LT.8) GO TO 1226
C ORTHORHOMBIC, MONOCLINIC AND TRICLINIC DONE ........
      IF(LATTIS.LT.10) GO TO 1212
      DO 1210 I = 2,NBG
      IF(NDX(4,I).NE.8) GO TO 1210
      IF(NDX(1,I).GT.NDX(2,I)) NDX(4,I)=16
      IF(NDX(1,I).LT.NDX(2,I)) NDX(4,I)=0
 1210 CONTINUE
      DO 1211 I = 2,NBG
      IF(NDX(4,I).NE.4) GO TO 1211
      IF(NDX(3,I).NE.0) GO TO 1211
      IF(NDX(1,I).GT.NDX(2,I)) NDX(4,I)=8
      IF(NDX(1,I).LT.NDX(2,I)) NDX(4,I)=0
 1211 CONTINUE
C 4/MMM TETRAG. PART DONE ........
C 4/M AND 4/MMM TO FINISH ........
 1212 DO 1213 I = 2,NBG
      IF(NDX(4,I).NE.2) GO TO 1213
      IF(NDX(1,I).NE.0) NDX(4,I)=4
      IF(NDX(2,I).NE.0) NDX(4,I)=0
 1213 CONTINUE
      DO 1214 I = 2,NBG
      IF(NDX(4,I).NE.4) GO TO 1214
      IF(NDX(3,I).EQ.0) GO TO 1214
      NDX(4,I)=8
      IF(NDX(2,I).EQ.0) NDX(4,I)=0
 1214 CONTINUE
      GO TO 1226
C TETRAGONAL DONE ........
 1215 IF(LATTIS.GT.17) GO TO 1217
      IF(LATTIS.NE.12) LATT=LATTIS-10
      IF(LATTIS.GT.15) LATT=LATTIS-13
      IF(LATTIS.LT.15) LAUE=1
      DO 1216 I = 1,LIMHKL
      CALL SETHKL(LATT,I,0,0,6)
      CALL SETHKL(LATT,I,I,0,12)
      CALL SETHKL(LATT,I,I,I,8)
      L=I+1
      DO 1216 J = L,LIMHKL
      IF(LAUE.EQ.0) CALL SETHKL(LATT,I,J,0,24)
      IF(LAUE.EQ.1) CALL SETHKL(LATT,I,J,0,12)
      IF(LAUE.EQ.1) CALL SETHKL(LATT,J,I,0,12)
      CALL SETHKL(LATT,I,I,J,24)
      CALL SETHKL(LATT,J,J,I,24)
      M=J+1
      DO 1216 K = M,LIMHKL
      IF(LAUE.EQ.0) CALL SETHKL(LATT,I,J,K,48)
      IF(LAUE.EQ.1) CALL SETHKL(LATT,I,J,K,24)
      IF(LAUE.EQ.1) CALL SETHKL(LATT,K,J,I,24)
 1216 CONTINUE
      GO TO 1226
C CUBIC DONE ........
 1217 IF(LATTIS.GT.19) GO TO 1221
      N=12
      IF(LATTIS.EQ.18) N=6
      DO 1220 I = 1,LIMHKL
      CALL SETHKL(0,I,0,0,6)
      CALL SETHKL(0,I,I,0,6)
      CALL SETHKL(0,I,-I,0,6)
      CALL SETHKL(0,I,I,I,2)
      CALL SETHKL(0,I,I,-I,6)
      L=I+1
      DO 1220 J = L,LIMHKL
      CALL SETHKL(0,J,I,I,6)
      CALL SETHKL(0,-J,I,I,6)
      CALL SETHKL(0,J,J,I,6)
      CALL SETHKL(0,J,J,-I,6)
      IF(N.EQ.12) GO TO 1218
      CALL SETHKL(0,J,I,0,6)
      CALL SETHKL(0,-J,I,0,6)
      CALL SETHKL(0,-I,I,J,6)
      CALL SETHKL(0,-J,J,I,6)
 1218 CALL SETHKL(0,I,J,0,N)
      CALL SETHKL(0,I,-J,0,N)
      CALL SETHKL(0,I,-I,J,N)
      CALL SETHKL(0,J,-J,I,N)
      M=J+1
      DO 1220 K = M,LIMHKL
      IF(N.EQ.12) GO TO 1219
      CALL SETHKL(0,J,K,I,6)
      CALL SETHKL(0,J,K,-I,6)
      CALL SETHKL(0,-J,K,I,6)
      CALL SETHKL(0,J,-K,I,6)
 1219 CALL SETHKL(0,K,J,I,N)
      CALL SETHKL(0,K,J,-I,N)
      CALL SETHKL(0,K,-J,I,N)
 1220 CALL SETHKL(0,-K,J,I,N)
      GO TO 1226
C RHOMBOHEDRAL DONE ........
 1221 N=0
      IF(LATTIS.LT.22) N=5
      M=0
      IF(LATTIS.EQ.20) M=5
      IF(LATTIS.EQ.20) LAUE=1
      IF(LATTIS.EQ.22) LAUE=1
      DO 1225 I = 1,LIMHKL
      L=I+1
      CALL SETHKL(0,I,0,0,6)
      CALL SETHKL(0,I,I,0,6)
      CALL SETHKL(0,0,0,I,2)
      DO 1223 K = 1,LIMHKL
      CALL SETHKL(M,I,I,K,12)
      CALL SETHKL(N,I,0,K,12)
      DO 1223 J = L,LIMHKL
      IF(LAUE.EQ.0) GO TO 1222
      CALL SETHKL(N,I,J,K,12)
      CALL SETHKL(N,J,I,K,12)
      GO TO 1223
 1222 CALL SETHKL(N,I,J,K,24)
 1223 CONTINUE
      DO 1225 J = 1,LIMHKL
      IF(LAUE.EQ.0) GO TO 1224
      CALL SETHKL(0,I,J,0,6)
      CALL SETHKL(0,J,I,0,6)
      GO TO 1225
 1224 CALL SETHKL(0,I,J,0,12)
 1225 CONTINUE
C HEXAGONAL LATTICES DONE ........
C REMOVE SCREW AND GLIDE ABSENCES ........
 1226 IF(NBG.EQ.NREFS) NREFS=-1
      DO 1230 L = 1,3
      M=L+1
      IF(M.EQ.4) M=1
      N=6-L-M
      J=ISCREW(L)
      K=IGLIDE(L)
      IF(J.EQ.0) GO TO 1228
      DO 1227 I = 2,NBG
      IF(NDX(M,I).NE.0.OR.NDX(N,I).NE.0) GO TO 1227
      IPAR=NDX(L,I)
      IF(MOD(IPAR,2).NE.0) NDX(4,I)=0
 1227 CONTINUE
 1228 IF(K.EQ.0) GO TO 1230
      DO 1229 I = 2,NBG
      IF(NDX(L,I).NE.0) GO TO 1229
      IPAR=NDX(K,I)
      IF(K.EQ.4) IPAR=NDX(M,I)+NDX(N,I)
      IF(MOD(IPAR,2).NE.0) NDX(4,I)=0
 1229 CONTINUE
 1230 CONTINUE
C REMOVE THE ZERO MULTIPLICITIES ........
      N=1
      DO 1232 I = 2,NBG
      IF(NDX(4,I).EQ.0) GO TO 1232
      N=N+1
      DO 1231 K = 1,5
 1231 NDX(K,N)=NDX(K,I)
 1232 CONTINUE
C SORT ON INCREASING RHO ........
      NBG=N-1
      DO 1234 I = 1,NBG
      L=I+1
      IRHO=MAXRH
      DO 1233 J = L,N
      IF(NDX(5,J).GT.IRHO) GO TO 1233
      IRHO=NDX(5,J)
      M=J
 1233 CONTINUE
      DO 1234 K = 1,5
      NDX(K,I)=NDX(K,M)
 1234 NDX(K,M)=NDX(K,L)
      WRITE(6,51) NBG,LATTIS
      DO 2938 I = 1,NBG
 2938 NDX(6,I)=NDX(4,I)
      IF(KOUT.EQ.1) WRITE(6,56) ((NDX(K,I),K=1,5),I=1,NBG)
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SETHKL(LATT,I,J,K,M)
      COMMON AM(45150),V(600),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NPEAKS,NBG,NDX(6,1500),STEP,
     2 WISE(6420,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ISTOP,KEY,NREFL,TWOPIS,KEYREF(600),
     4 P(600),KI(600),PWISE(9),H(3),RHO,NOUT,NPM,MINPTU,MAXPTU,
     5 TWOPI,YC,WAVETC,DC(600),YSUM,YMASK,YAV,NMAX,NEXT,NSIDE,IW,
     6 KEYSUM,MSTEPU,SNTL(600),LVAR(600),NO,KOUT,
     7 SLACK(10),SIG,REFSIG,LOT,K15
      DIMENSION HJ(3),HK(3)
      IF(NBG.EQ.NREFS) RETURN
      IF(LATT.LT.2) GO TO 1251
      IF(LATT.GT.4) GO TO 1251
      L=I+J
      IF(LATT.EQ.4) L=L+K
      IF(MOD(L,2).NE.0) RETURN
      IF(LATT.NE.3) GO TO 1251
      L=I+K
      IF(MOD(L,2).NE.0) RETURN
 1251 HK(1)=I
      HK(2)=J
      HK(3)=K
      HJ(1)=ZT(1,1)*HK(1)+ZT(2,1)*HK(2)+ZT(3,1)*HK(3)
      HJ(2)=ZT(1,2)*HK(1)+ZT(2,2)*HK(2)+ZT(3,2)*HK(3)
      HJ(3)=ZT(1,3)*HK(1)+ZT(2,3)*HK(2)+ZT(3,3)*HK(3)
      IRHO=2500.*(HJ(1)**2+HJ(2)**2+HJ(3)**2)
      IF(IRHO.GT.MAXRH) RETURN
      IF(IRHO.LT.MINRH) RETURN
      NBG=NBG+1
      NDX(1,NBG)=I
      NDX(2,NBG)=J
      NDX(3,NBG)=K
      NDX(4,NBG)=M
      NDX(5,NBG)=IRHO
      IF(LATT.NE.5) RETURN
      NDX(4,NBG)=M/2
      IF(NBG.EQ.NREFS) RETURN
      NBG=NBG+1
      NDX(1,NBG)=I
      NDX(2,NBG)=J
      NDX(3,NBG)=-K
      NDX(4,NBG)=M/2
      NDX(5,NBG)=IRHO
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE TEST
      COMMON AM(45150),V(600),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NPEAKS,NBG,NDX(6,1500),STEP,
     2 WISE(6420,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ISTOP,KEY,NREFL,TWOPIS,KEYREF(600),
     4 P(600),KI(600),PWISE(9),H(3),RHO,NOUT,NPM,MINPTU,MAXPTU,
     5 TWOPI,YC,WAVETC,DC(600),YSUM,YMASK,YAV,NMAX,NEXT,NSIDE,IW,
     6 KEYSUM,MSTEPU,SNTL(600),LVAR(600),NO,KOUT,
     7 SLACK(10),SIG,REFSIG,LOT,K15
   61 FORMAT(72I1)
   62 FORMAT(' TEMPERATURE FACTOR OF ATOM',I3,' ISNT POSITIVE-DEFINITE')
   63 FORMAT(2F8.4,2X,3F8.4,2X,3F9.5,2X,6F8.4)
   65 FORMAT(/' KEYS FOR THIS CYCLE',(/21X,36I2))
   66 FORMAT(52X,F13.7)
   68 FORMAT(I8,' PARAMETERS AND THEIR INDICES INPUT FROM LAST CYCLE')
   69 FORMAT(' REFLECTION',I5,' IS ZEROED FROM',F12.6)
      IF(NEXT.NE.0) GO TO 1410
      GO TO 1402
 1400 DO 1401 J = 1,NP
 1401 READ(5,66) P(J)
      IC=-1
 1402 NEXT=0
      IF(IC.EQ.NC) GO TO 1407
      NPQ=NP
      IF(LOT.NE.0) NPQ=K15
      READ(5,61) (KI(I),I=1,NPQ)
      KEY=KI(1)+KI(2)+KI(3)+KI(4)+1
      KEYSUM=0
      DO 1403 J = 1,K15
 1403 KEYSUM=KEYSUM+KI(J)
      IF(NPQ.EQ.K15) GO TO 1405
      M=KEYSUM
C     THE FOLLOWING IS OVERWRITTEN IF LOT IS NON-ZERO
      DO 1404 I = 16,NP
      M=M+KI(I)
 1404 KEYREF(I-K15)=M
C     SET UP LVAR( ) FOR VARIABLES EXCEPT SCALE AND BACKGROUND
 1405 NV=0
      DO 1406 J=1,NPQ
      IF(KI(J).NE.1) GO TO 1406
      NV=NV+1
      LVAR(NV)=J
 1406 CONTINUE
      WRITE(6,65)  (KI(I),I=1,NPQ)
 1407 I=NP+1
      IF(IC.EQ.0) GO TO 1410
      DO 1408 J=1,9
      PWISE(J)=P(J)
      IF(J.GT.3) GO TO 1408
      RCP(J  )=P(J+9 )
      RCP(J+3)=P(J+12)/DEGRAD
 1408 CONTINUE
      CALL CELLAZ
      DO 1409 J = K15,NP
      IF(P(J).GE.0.0) GO TO 1409
      WRITE(6,69) J,P(J)
      P(J)=0.0
 1409 P(J)=P(J)*PWISE(1)
      PWISE(1)=1.0
 1410 DO 1411 J = 1,9
      P(J)=PWISE(J)
      IF(J.GT.3) GO TO 1411
      P(J+9 )=RCP(J  )
      P(J+12)=RCP(J+3)*DEGRAD
 1411 CONTINUE
      IF(NEXT.NE.0) GO TO 1400
      IF(NV.GT.NVM) CALL FINISH(8)
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE CELLAZ
      COMMON AM(45150),V(600),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB
      DIMENSION RECIP(6),SN(3),CS(3),RA(3)
   71 FORMAT(//' NEW CELL =',6F12.6/' CELL VOLUME' ,F15.5/' RECIPROCAL C
     1ELL PARAMETERS',6F10.5/' TRANSFORMATIONS'//(3F12.6,9X,3F12.6))
   72 FORMAT(' SPECIAL RHOMBOHEDRAL SECTION')
      IF(IRHOMB.NE.0) GO TO 1503
      DO 1501 I=1,3
      SN(I)=SIN(RCP(I+3)*DEGRAD)
 1501 CS(I)=COS(RCP(I+3)*DEGRAD)
      RA(1)=CS(2)*CS(3)-CS(1)
      RA(2)=CS(3)*CS(1)-CS(2)
      RA(3)=CS(1)*CS(2)-CS(3)
      BIGD=1.-CS(1)**2-CS(2)**2-CS(3)**2+2.*CS(1)*CS(2)*CS(3)
      AROOT=SQRT(BIGD)
      AT(1,1)=RCP(1)
      AT(1,2)=RCP(2)*CS(3)
      AT(1,3)=RCP(3)*CS(2)
      AT(2,1)=0.
      AT(2,2)=RCP(2)*SN(3)
      AT(2,3)=-RCP(3)*RA(1)/SN(3)
      AT(3,1)=0.
      AT(3,2)=0.
      AT(3,3)=RCP(3)*AROOT/SN(3)
      CELL=AT(1,1)*AT(2,2)*AT(3,3)
      ZT(1,1)=1./RCP(1)
      ZT(1,2)=-CS(3)/(RCP(1)*SN(3))
      ZT(1,3)=RA(2)/(RCP(1)*SN(3)*AROOT)
      ZT(2,1)=0.
      ZT(2,2)=1./(RCP(2)*SN(3))
      ZT(2,3)=RA(1)/(RCP(2)*SN(3)*AROOT)
      ZT(3,1)=0.
      ZT(3,2)=0.
      ZT(3,3)=SN(3)/(RCP(3)*AROOT)
      DO 1502 I = 1,3
      L=I+3
      J=I+1
      IF(J.EQ.4) J=1
      K=6-I-J
      RECIP(I)=SN(I)/(RCP(I)*AROOT)
      RECIP(L)=ACOS(RA(I)/(SN(J)*SN(K)))/DEGRAD
      FAC=0.5/BIGD
      DCELL(I,I  )=-FAC*SN(I)**2/RCP(I)**3
      DCELL(I,J  )=0.0
      DCELL(I,K  )=0.0
      DCELL(I,I+3)=0.0
      DCELL(I,J+3)=-FAC*RA(J)/(RCP(K)*RCP(I)**2)
      DCELL(I,K+3)=-FAC*RA(K)/(RCP(J)*RCP(I)**2)
      FAC=FAC*SN(I)/BIGD
      DCELL(L,I  )=FAC*(CS(I)*BIGD+SN(I)**2*RA(I))/RCP(I)**2
      DCELL(L,J  )=FAC*            SN(J)**2*RA(I) /RCP(J)**2
      DCELL(L,K  )=FAC*            SN(K)**2*RA(I) /RCP(K)**2
      DCELL(L,I+3)=FAC*(2.*RA(I)*RA(I)+BIGD      )/(RCP(J)*RCP(K))
      DCELL(L,J+3)=FAC*(2.*RA(I)*RA(J)-BIGD*CS(K))/(RCP(K)*RCP(I))
      DCELL(L,K+3)=FAC*(2.*RA(I)*RA(K)-BIGD*CS(J))/(RCP(I)*RCP(J))
 1502 CONTINUE
      GO TO 1508
C     RHOMBOHEDRAL: SPECIAL CASE
 1503 ROOT2=SQRT(2.0)
      ROOT3=SQRT(3.0)
      AXIS=RCP(1)
      ALFA=RCP(4)*DEGRAD
      PSI=ACOS(1.0/ROOT3)
      DIF=ASIN(2.0*SIN(ALFA*0.5)/ROOT3)
      PHI=PSI-DIF
      PHIP=-COS(ALFA*0.5)/(ROOT3*COS(DIF))
      COSPHI=COS(PHI)
      SINPHI=SIN(PHI)
      BET=AXIS*COSPHI
      GAM=AXIS*SINPHI/ROOT2
      DEL=BET*(BET+GAM)-2.0*GAM**2
      ZDIA=(BET+GAM)/DEL
      ZOFF=-GAM/DEL
      SN(1)=COSPHI
      CS(1)=SINPHI/ROOT2
      SN(2)=-AXIS*SINPHI*PHIP
      CS(2)= AXIS*COSPHI*PHIP/ROOT2
      DO 1504 I = 1,6
      DO 1504 J = 1,6
 1504 DCELL(I,J)=0.0
      DO 1506 I = 1,2
      DELP=(2.0*BET+GAM)*SN(I)+(BET-4.0*GAM)*CS(I)
      ZDIAP=((SN(I)+CS(I))*DEL-(BET+GAM)*DELP)/DEL**2
      ZOFFP=(GAM*DELP-DEL*CS(I))/DEL**2
      DIA=0.5*ZDIA*ZDIAP+ZOFF*ZOFFP
      OFF=ZDIA*ZOFFP+ZOFF*(ZDIAP+ZOFFP)
      K=3*I-2
      DO 1505 J = 1,3
      DCELL(K,J  )=DIA
 1505 DCELL(K,J+3)=OFF
 1506 CONTINUE
      ASTAR=SQRT((BET+GAM)**2+2.0*GAM**2)/DEL
      ALFST=2.0*ATAN((ROOT3*COS(DIF))/1.0)/DEGRAD
      CELL=(BET-GAM)*DEL
      DO 1507 I = 1,3
      RECIP(I  )=ASTAR
      RECIP(I+3)=ALFST
      DO 1507 J = 1,3
      AT(I,J)=GAM
      AT(I,I)=BET
      ZT(I,J)=ZOFF
 1507 ZT(I,I)=ZDIA
      WRITE(6,72)
 1508 CONTINUE
      WRITE(6,71) (RCP(I),I=1,6),CELL,(RECIP(I),I=1,6),((AT(I,J),J=1,3
     1 ),(ZT(I,J),J=1,3),I=1,3)
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE MATINV (N,NFAIL)
      COMMON AM(45150),V(600)
      K=1
      NFAIL=1
      IF(N.LT.1)  GO TO 1615
      IF(N.GT.1)  GO TO 1601
      AM(1)=1.0/AM(1)
      GO TO 1614
 1601 DO 1607 M=1,N
      IMAX=M-1
      DO 1606 L=M,N
      SUMA=0.0
      KLI=L
      KMI=M
      IF(IMAX.EQ.0)  GO TO 1603
      DO 1602 I=1,IMAX
      SUMA=SUMA+AM(KLI)*AM(KMI)
      J=N-I
      KLI=KLI+J
 1602 KMI=KMI+J
 1603 TERM=AM(K)-SUMA
      IF(L.GT.M)  GO TO 1605
      IF(TERM.GT.0.0)  GO TO 1604
      NFAIL=K
      GO TO 1615
 1604 DENOM=SQRT (TERM)
      AM(K)=DENOM
      GO TO 1606
 1605 AM(K)=TERM/DENOM
 1606 K=K+1
 1607 CONTINUE
      AM(1)=1.0/AM(1)
      KDM=1
      DO 1610 L=2,N
      KDM=KDM+N-L+2
      TERM = 1.0/AM(KDM)
      AM(KDM)=TERM
      KMI=0
      KLI=L
      IMAX=L-1
      DO 1609 M=1,IMAX
      K=KLI
      SUMA=0.0
      DO 1608 I=M,IMAX
      II=KMI+I
      SUMA=SUMA-AM(KLI)*AM(II)
 1608 KLI=KLI+N-I
      AM(K)=SUMA*TERM
      J=N-M
      KLI=K+J
 1609 KMI=KMI+J
 1610 CONTINUE
      K=1
      DO 1613 M=1,N
      KLI=K
      DO 1612 L=M,N
      KMI=K
      IMAX=N-L+1
      SUMA=0.0
      DO 1611 I=1,IMAX
      SUMA=SUMA+AM(KLI)*AM(KMI)
      KLI=KLI+1
 1611 KMI=KMI+1
      AM(K)=SUMA
 1612 K=K+1
 1613 CONTINUE
 1614 NFAIL=0
 1615 RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE FINISH(IWHY)
      IF(IWHY.NE.12) GO TO 1701
      WRITE(6,82)
      STOP
 1701 WRITE(6,81) IWHY
      STOP
   81 FORMAT(' STOPPED FOR REASON',I6/24X,'1 : NO PARAMs TO BE VARIED'/
     124X,'2 : CONTRIBUTION OUTSIDE BLOCK ARRAY, REORDER REFLECTIONS'/2
     2X,'3 : WRONG NUMBER OF PARAMETERS'/24X,'4 : TOO MANY PARAMS'/24X,
     3'5 : TOO MANY ATOMS'/24X,'6 : TOO MANY REFLECTIONS'/24X,'7 : TOO W
     4IDE A SCAN'/24X,'8 : TOO MANY VARIABLES'/24X,'9 : UNPHYSICAL TEMP.
     5 FACTOR'/23X,'10 : SINGULAR LEAST SQUARES MATRIX'/23X,'11 : ZERO D 
     6IAGONAL MATRIX ELEMENT'/23X,'12 : SUCCESSFUL RUN'/23X,'13 : MINPT 
     7= 0'/23X,'14 : INDEXED LINE OUTSIDE SCAN'/)
   82 FORMAT(' SUCCESSFUL RUN')
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE PLOTLP
      COMMON AM(45150),V(600),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NPEAKS,NBG,NDX(6,1500),STEP,
     2 WISE(6420,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ISTOP,KEY,NREFL,TWOPIS,KEYREF(600),
     4 P(600),KI(600),PWISE(9),H(3),RHO,NOUT,NPM,MINPTU,MAXPTU,
     5 TWOPI,YC,WAVETC,DC(600),YSUM,YMASK,YAV,NMAX,NEXT,NSIDE,IW,
     6 KEYSUM,MSTEPU,SNTL(600),LVAR(600),NO,KOUT,
     7 SLACK(10),SIG,REFSIG,LOT,K15
      DIMENSION LABEL(12),IOUT(121),IDOTS(121),HJ(3),HK(3),INDT(5)
      DATA ISPACE/' '/,IDOT/'.'/,IPLUS/'X'/,ISTAR/'*'/,IOBS/'O'/,IM/'-'/
   91 FORMAT(3X,12I10)
   92 FORMAT(1X,121A1,F7.2)
   93 FORMAT(1X,121A1,3I3)
   94 FORMAT(' REORDERED REFLECTIONS'/10(I7,2I3))
      DO 1806 NREFL = 2,NBG
      HK(1)=NDX(1,NREFL)
      HK(2)=NDX(2,NREFL)
      HK(3)=NDX(3,NREFL)
      HJ(1)=ZT(1,1)*HK(1)+ZT(2,1)*HK(2)+ZT(3,1)*HK(3)
      HJ(2)=ZT(1,2)*HK(1)+ZT(2,2)*HK(2)+ZT(3,2)*HK(3)
      HJ(3)=ZT(1,3)*HK(1)+ZT(2,3)*HK(2)+ZT(3,3)*HK(3)
      SINTH=SQRT(0.25*(HJ(1)**2+HJ(2)**2+HJ(3)**2))*WAVETC
      NDX(5,NREFL)=2.0*ASIN(SINTH)/(DEGRAD*STEP)+0.5
C     THIS IS THE SCAN POINT AS AN INTEGER
      INDX=NDX(5,NREFL)
      K=NREFL
      J=0
      DO 1801 I = 1,20
      IF(K.EQ.1) GO TO 1802
      K=K-1
      IF(INDX.GE.NDX(5,K)) GO TO 1802
 1801 J=J+1
 1802 IF(J.EQ.0) GO TO 1806
      DO 1803 I = 1,5
 1803 INDT(I)=NDX(I,NREFL)
      DO 1804 I = 1,J
      K=NREFL-I
      L=K+1
      DO 1804 M = 1,5
 1804 NDX(M,L)=NDX(M,K)
      DO 1805 I = 1,5
 1805 NDX(I,K)=INDT(I)
 1806 CONTINUE
      WRITE(8,94) ((NDX(I,J),I=1,3),J=1,NBG)
      WMAX=0.0
      DO 1807 I = MINPTU,MAXPTU
      IF(WISE(I,1).GT.WMAX) WMAX=WISE(I,1)
      IF(WISE(I,2).GT.WMAX) WMAX=WISE(I,2)
 1807 CONTINUE
      FSCALE=WMAX/116.
      DO 1808 I = 1,12
 1808 LABEL(I)=10.0*FSCALE*FLOAT(I)
      WRITE(8,91) LABEL
      DO 1809 J = 1,121
      IDOTS(J)=ISPACE
      IF(MOD(J,10).EQ.1) IDOTS(J)=IDOT
 1809 CONTINUE
      WRITE(8,92) IDOTS
      NDX(5,NBG+1)=NMAX+1
      DO 1810 IFIND = 1,NBG
      IF(NDX(5,IFIND).GE.MINPTU) GO TO 1811
 1810 CONTINUE
 1811 L=PWISE(5)/STEP + 0.05
      DO 1821 I = MINPTU,MAXPTU
      M=I
      DO 1812 J = 2,120
 1812 IOUT(J)=ISPACE
      IOUT(1)=IDOT
      IOUT(121)=IDOT
      IF(MOD(M,50).NE.0) GO TO 1814
      DO 1813 J = 1,121
 1813 IOUT(J)=IDOTS(J)
 1814 MOBS=WISE(I,1)/FSCALE+1.5
      MCAL=WISE(I,2)/FSCALE+1.5
      IREFL=0
 1815 IF(NDX(5,IFIND).GT.M) GO TO 1816
      IREFL=IREFL+1
      IFIND=IFIND+1
      GO TO 1815
 1816 IF(IREFL.EQ.0) GO TO 1818
      IREFL=121-IREFL
      DO 1817 J = IREFL,120
 1817 IOUT(J)=IM
 1818 IF(MOBS.GT.0) IOUT(MOBS)=IOBS
      IF(MCAL.GT.0) IOUT(MCAL)=ISTAR
      IF(MOBS.LT.0) IOUT(2)=IM
      IF(MOBS.LT.0) MOBS=2
      IF(IOUT(MOBS).EQ.ISTAR) IOUT(MOBS)=IPLUS
      IF(IREFL.EQ.0) GO TO 1819
C     PRINTS LAST NDX SET FOR THIS POINT
      IFINDM=IFIND-1
      WRITE(8,93) IOUT,(NDX(J,IFINDM),J=1,3)
      GO TO 1821
 1819 IF(MOD(M,10).EQ.0) GO TO 1820
      WRITE(8,92) IOUT
      GO TO 1821
 1820 R=STEP*FLOAT(M)
      WRITE(8,92) IOUT,R
 1821 CONTINUE
      WRITE(8,92) IDOTS
      WRITE(8,91) LABEL
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE PROD (A,B,C,N)
      DIMENSION A(3,40),B(3,3),C(3,40)
      DO 1901 I=1,3
      DO 1901 J=1,N
 1901 A(I,J)=B(I,1)*C(1,J)+B(I,2)*C(2,J)+B(I,3)*C(3,J)
      RETURN
      END
C
C-----------------------------------------------------------------------