C*********************************************************************
C
C LEAST SQUARES UNIT CELL REFINEMENT INCLUDING INDEXING IN FORTRAN 77
C
C  ORIGINALY FROM THE GEOLOGICAL SURVEY BY APPLEMAN, HANDWERKER AND
C  EVANS
C
C  MODIFIED AT NBS BY
C
C  C. R. HUBBARD AND S. LEDDERMAN     SUMMER 1982
C
C  MODIFIED AT ALFRED BY
C
C  D. P. MATHEIS     SUMMER 1990
C
************************************************************************
***********************************************************************
***********************************************************************
* common block and character assignments
* for alsq : Alfred lattice parameter least square
* refinement program modified from NBSLSQ
* by D. P. Matheis summer 1990
***********************************************************************
* common blocks
***********************************************************************
*
*     maxobs = max # of observed lines allowed
*     maxgen = max # 0f lines that will be generated
*
*     qmod: coefficients of the q* structure modulation vector
*           q* = qj(1) a* + qj(2) b* + qj(3) c*
*           qmodul= .true. if modulation present
*           mmin, mmax = limits for modulation indicie m
*
*     recp: x(1)=a*,  x(2)=b*,  x(3)=c*
*           u(1)=cos(alph*), u(2)=cos(beta*), u(3)=cos(gamm*)
*           w(1)=sin(alph*), u(2)=sin(beta*), u(3)=sin(gamm*)
*           z(1)=a*2, z(2)=b*2, z(3)=c*2, 
*           z(5)=a*b*Cos(Gam), z(6)=b*c*Cos(alph), z(7)=c*a*Cos(Beta)
*           vr = reciprocal cell volume
*
*     real: ao(1)=a, ao(2)=b, ao(3)=c
*           alo(1)=alpha, alo(2)=beta, alo(3)=gamma
*           csa(1)=cos(alpha), csa(2)=cos(beta), csa(3)=cos(gamma)
*           sna(1)=sin(alpha), sna(2)=sin(beta), sna(3)=sin(gamma)
*           vc = direct cell volume
*
*     io:  ntout & ntin = output and input unit numbers
*
*     blka: 
*
*     hdgen:  ihg(i,1) = ith h
*             ihg(i,2) = ith k
*             ihg(i,3) = ith l
*             ihg(i,4) = ith m
* 
*     ppdata: obsth(i) -> observed 1-theta values
*             dsp(i)   -> observed d-spacing values
*             ih(i,j)  -> input hklm assignments
*             pwl      -> wavelenth
*             numobs   -> number of observed reflections
*             int(i)   -> observed peak intensities
*
***********************************************************************
* character variables
***********************************************************************
*     cblkg: ihead=file header from d/I file
*
***********************************************************************
*      implicit real*8(a-h, o-z)
      PARAMETER (MAXOBS=200, MAXGEN=1000)
      common/qmod/qj(3),mmin,mmax,qmodul,qrefin,nq,iqtyp,iq(3)
      COMMON/PPCON/JTYP,ITYP,TOL,THEMX,RESOLV
      COMMON/CPPCON/IALF12
      COMMON/RECP/Z(10),Q,AR(3),X(4),VR,U(3),W(3),VRO
      COMMON/REAL/AO(3),ALO(3),VCO,AD(3),ALS(3),VC,CSA(3),SNA(3)
      COMMON/ADATA/IFI(2,27),IFJ(3,27),IHKL(6,3),JHKL(10,3)
      COMMON/CADATA/IALF(13), IALFL(13), CSC(8), cscl(8)
      COMMON/IO/NTOUT,NTIN,NEXTN
      COMMON/BLKA/A(54),UO(3),B(9),bb(9,9)
      COMMON/ALL/N,NUMCYC,ntot
      COMMON/SCRT/E(9),ee(9,9),NH(4),qh(3),f(9),ff(9,10)
      common/sim/delta(10),eps,indic,nrc,deter
      COMMON/BLKB/CR(9)
      COMMON/BLKF/SW(200),WL(200)
      COMMON/PPDATA/OBSTH(200),DSP(200),IH(200,4),PWL,
     * NUMOBS,INT(200)
      COMMON/CPPDAT/ IHED(200),INTED(200)
      COMMON/BLKC/THCMP,D,CONST,SWT,DIFF,DOBS,THTSE,THOBS,
     * RSDMAX,IRSD,SUM,MMN,prthtse
      COMMON/CBLKC/JA
      COMMON/PPEXTN/IOP(27),JOP(27),NOPT
      COMMON/HDGEN/IHG(1000,4),DCM(1000),NR,DMIN
      COMMON/BLKG/ICYCLE,IMN,IWCYC,IWTYP
      COMMON/CBLKG/IHEAD
      COMMON/BLKJ/TOLMN , TOLMX
      COMMON/BLKH/NCYC,AA(7),DLST(200),IOC(200),THTMX
      COMMON/CONTR/IXPD,IJOB,IPHKL,IDOLSQ
      COMMON/CCONTR/ ID0
      COMMON/SWFG/R
      COMMON/CELLI/AI,BI,CI,ALPHI,BETI,GAMI,VI,ERRI(7)
      COMMON/CCELLI/ SPGI,SPGIA,SPGIN,ICRS
      COMMON/CHAR/ITIH(7),ITCRS(8),ITSG(7),ITED(4)
      COMMON/EDITR/IED1(6)
      COMMON/UNQ/ IFA,IFB,IFC,IRH
      COMMON/LSQRD/CARD,ITEST
      COMMON/CCREAD/ICRYST
      COMMON/RDGOOD/RGOOD,SYSEXT
      COMMON/IERRS/ IWARN,IERR
      COMMON/CRYDAT/CELL(6),SDCELL(6)
      COMMON/CRYDA1/CRSYS,SPGRP
      COMMON/PEAKS/NPKS,NPKMX,NPFLG(6,250),PKPRM( 3,250),NPKSQ(250),
     * IHKLPF(3,250)
      COMMON/PKFIL/NAME,TITLE
C
      CHARACTER* 1 CRSYS
      CHARACTER* 8 SPGRP
      CHARACTER* 6 NAME
      CHARACTER*80 TITLE
      CHARACTER*7 ID0
      CHARACTER*1 IHED,INTED
      CHARACTER*1 SPGIA,ICRS
      CHARACTER*4 SPGIN
      CHARACTER*8 SPGI
      CHARACTER*1 ITSG,ITCRS,ITED,ITIH
      CHARACTER*6 IALF, ialfl, ITEST, IALF12, JA,ICOND(10)
      CHARACTER*5 CSC,cscl,ICRYST
      CHARACTER*3 ICLASS(10)
      CHARACTER*68 IHEAD
      CHARACTER*80 CARD
      CHARACTER*8 ISGC
      CHARACTER*1 IED1,IBLANK,ICENT
      CHARACTER*1 IHN,IHS,IHU
      CHARACTER*4 IH205
      character*3 junk
      DATA IH205/'205 '/
      DATA IHN,IHS,IHU/'N','*','U'/
      DIMENSION JCLASS(27),JCOND(27),XO(4)
      dimension ss(3)
      LOGICAL GOOD,RGOOD,SYSEXT,DONE,qmodul,qrefin
      REAL TEMP(4)
      DATA (ICLASS(L),L=1,10)/'HKL','HHL','0KL','H0L','HK0','HH0','H00',
     * '0K0','00L','HLL'/,(ICOND(L),L=1,10)/'   H+K','   H+L','   K+L',
     * ' H+K+L','-H+K+L','     H','     K','     L','  2H+L','  H+2L'/,
     * (JCLASS(L),L=1,27)/6*1,4*2,3*3,3*4,3*5,6,7,8,9,4*10/,(JCOND(L),
     * L=1,27)/1,2,3,0,4,5,6,8,2,9,7,8,3,6,8,2,6,7,1,6,6,7,8,6,8,2,10/
      DATA IBLANK/' '/
      DATA IALF/'TITLE ','PARAM ','SYSEXT','D     ','FINISH',
     1 '2-    ','      ',' R    ','xxxxxx','DEMAND','PKFILE','yyyyyy',
     2 'zzzzzz'/
      DATA IALFL/'title ','param ','sysext','d     ','finish',
     1 '2-    ','      ',' r    ','xxxxxx','demand','pkfile','yyyyyy',
     2 'zzzzzz'/
      DATA CSC/'CUBIC','TETRA','ORTHO','MONOC','RHOMB','TRICL',
     1  'HEXAG','HEXAR'/
      DATA CSCl/'cubic','tetra','ortho','monoc','rhomb','tricl',
     1  'hexag','hexar'/
      DATA NTIN,NTOUT,NEXTN/5,6,20/
      DATA ITCRS/'A','M','O','T','H','R','C',' '/
      DATA ITSG /'P','A','B','C','F','I','R'/
      DATA ITED /' ','a','b','c'/
      DATA ITIH /'d','e','f','g','h','i','j'/
      DATA IFI/12*0,1,0,1,0,1,0,1,0,4,0,4,0,4,0,5,0,5,0,5,0,6,0,6,0,6,0,
     1 1,6,5,6,4,6,4,5,2,0,2,0,2,0,2,0/
      DATA IFJ/1,0,0,2,0,0,3,0,0,1,2,3,4,0,0,5,0,0,6,0,0,8,0,0,2,0,0,
     1           9,0,0,7,0,0,8,0,0,3,0,0,6,0,0,8,0,0,2,0,0,6,0,0,7,0,0,
     1           1,0,0,6,0,0,6,0,0,7,0,0,8,0,0,6,0,0,8,0,0,2,0,0,10,0,0/
      DATA IHKL/1,0,1,1,0,0,-1,1,0,0,1,0,0,-1,-1,0,0,1/
      DATA JHKL/1,1,0,1,-1,1,0,0,2,1,1,0,3*1,0,1,4*0,4*1,0,0,1,1,2/
C
C-----INITIALIZE
C
 100  R = 180.00/3.14154265
      indic = 0
      eps = .0005
      nrc = 10
      nq = 0
      do i=1,3
        iq(i) = 0.0
        qj(i) = 0.00
      enddo
      IDOLSQ = 1
      NERRT=0
      qmodul = .false.
      DO I = 1, 200, 1
        INT(I) = 0
        INTED(I) = IBLANK
      enddo
      DO I = 1, 7, 1
        ERRI(I) = 0.0
      enddo
      IJOB = 1
      ID0 = IBLANK
* 100  CONTINUE
      NTOUT = 6
      WRITE(NTOUT,1)
 1     FORMAT(6x,'***Alfred Version of JCPDS-NBS Geological Survey***'/
     *        6x,'        **Lattice Parameter Least Squares**'//
     *        6x,'          Version for VAX       July 1990'////)
C
C-----INITIALIZE FOR EACH SET
C
      DO I = 1, 200, 1
        IHED(I) = IBLANK
      enddo
      IED1(2) = IBLANK
      SPGI = IBLANK
      SPGIN= IBLANK
      SPGIA= IBLANK
      IERR = 0
      IWARN= 0
C
C-----DEMAND MODE RUN, CALL DMREAD
C
      CALL DMREAD
C
C-----CHECK IF DATA READ WAS GOOD, IF NOT START NEW SET
C
      IF (.NOT. (RGOOD)) GO TO 100
      WRITE (NTOUT,6) IHEAD
 6     FORMAT ('0' / '0', 'Title of set: ', A68 //)
      WRITE (NTOUT,7) ICRYST
 7     FORMAT (' ', 'Crystal system = ', A5 / )
      IF(.NOT.SYSEXT) WRITE (NTOUT,8) SPGI
 8     FORMAT (' ', 'Space group = ', A8)
C
C-----OPEN NEXTN FILE IF NECESSARY
C
      IF (.NOT. (SYSEXT)) OPEN (UNIT = NEXTN, ERR = 870,
     * FILE = 'POWD:[PUBLIC.EXE.LSQ]NBSLSQ.EXT', 
     1 STATUS = 'OLD',READONLY)
 200  IF(ICRYST.EQ.CSC(1).or.icryst.eq.cscl(1))ICRS = ITCRS(7)
      IF(ICRYST.EQ.CSC(2).or.icryst.eq.cscl(2))ICRS = ITCRS(4)
      IF(ICRYST.EQ.CSC(3).or.icryst.eq.cscl(3))ICRS = ITCRS(3)
      IF(ICRYST.EQ.CSC(4).or.icryst.eq.cscl(4))ICRS = ITCRS(2)
      IF(ICRYST.EQ.CSC(5).or.icryst.eq.cscl(5))THEN
         ICRS = ITCRS(6)
         IRH  = 0
      END IF
      IF(ICRYST.EQ.CSC(6).or.icryst.eq.cscl(6))ICRS = ITCRS(1)
      IF(ICRYST.EQ.CSC(7).or.icryst.eq.cscl(7))ICRS = ITCRS(5)
      IF(ICRYST.EQ.CSC(8).or.icryst.eq.cscl(8))THEN
         ICRS = ITCRS(6)
         IRH = 1
      END IF
      CALL PRMCHK(ICRS,N,ITYP,JTYP)
      IF (SYSEXT) GO TO 230
C
C-----IF RHOMBOHEDRAL SYSTEM ON RHOMBOHEDRAL AXES THEN EXTINCTIONS MUST
C-----BE INPUT EXPLICITLY WITH SYSEXT CARD
C
      IF(ICRYST.NE.CSC(5).and.icryst.ne.cscl(5)) GO TO 220
      IF(ICRYST.EQ.CSC(5).or.icryst.eq.cscl(5)) WRITE(NTOUT,9)
 9      FORMAT('0','**WARNING** Extinction conditions for ',
     * ' rhombohedral lattice not in file **WARNING**'/ ' ',11X,
     * 'you may wish to explicitly input them using SYSEXT card')
      GOTO 230
 220  CALL CHKSPG
      IPASS = 1
      ISGC = SPGI
      CALL SPG(ISGC, IPASS)
C
C-----TEST THAT THE INPUT SPACE GROUP WAS IN THE FILE FROM NBS*AIDS80
C
      IF(SPGIN.EQ.'  0 ')then
        WRITE(NTOUT,10)SPGI
 10       FORMAT('0','**WARNING** Input space group "',A8,'" was 
     *   not in the file **WARNING**')
        IWARN=IWARN+1
      END IF
      CALL FINDEX
      close(unit=nextn)
 230  CONTINUE
*
* DO BOOKKEEPING FOR SETTING UP EXTINCTION RELATIONS TO USE
*
      NOPT=0
      IF (JTYP.EQ.8) IOP(6)=3
      DO L=1,27
        IF(IOP(L).LT.2.OR.IOP(L).GT.6.OR.IOP(L).EQ.5)GOTO 240
        NOPT=NOPT+1
        JOP(NOPT)=L
 240  enddo
C
C-----ROUTINE NOW ONLY CHECKS HKL VS. EXTINCTION CODES
C     using input hkl values if present
C-----BEGIN LOOP THROUGH INPUT REFLECTIONS FOR EXTINCTION CHECK
      IF(NOPT .LT. 1) GO TO 260
      DO 250  I=1,NUMOBS
C
C-----TEST IF THIS IS AN INDEXED REFLECTION.  DO NOT TEST THOSE
C-----MARKED WITH "*" OR "U".
C
        IF(IHED(I).EQ.IHS .OR. IHED(I).EQ.IHU) GOTO 250
        JH = IH(I,1)
        JK=IH(I,2)
        JL=IH(I,3)
        CALL CHKEXT(JH,JK,JL,JEXT)
        IF(JEXT .EQ.0) GO TO 250
C-----WRITE OUT ERROR MESSAGE FOR ANY FORBIDDEN INDICIES
        ANGLE = OBSTH(I)
        ANGLE=ANGLE+ANGLE
        IF(SPGIN.NE.'    ') WRITE(NTOUT,11) JH,JK,JL,ANGLE,DSP(I),SPGI
        IF(SPGIN.EQ.'    ') WRITE(NTOUT,12) JH,JK,JL,ANGLE,DSP(I)
        IERR=IERR+1
 11      FORMAT(1H ,20H**ERROR** Reflection,3I3,7H/angle=,F7.3,3H/d=,
     *    F7.4,20H should be absent in,1X,A8)
 12      FORMAT(1H ,20H**ERROR** Reflection,3I3,7H/angle=,F7.3,3H/d=,
     *    F7.4,  ' should be absent')
        IHED(I)=IHN
 250  CONTINUE
*
* prepare to assign lattice parameter values not given due to symmetry
*
 260  AI = AO(1)
      BI = AO(2)
      CI = AO(3)
      ALPHI = alo(1)
      BETI = alo(2)
      GAMI = alo(3)
      CALL CELLCK(GOOD)         !assign lattice parameters w/r symmetry
      IF (.NOT. GOOD) GO TO 100
C
C-----reset lattice parameters based on any changes made by cellck
C
      AO(1) = AI
      AO(2) = BI
      AO(3) = CI
      alo(1) = alphi
      alo(2) = beti
      alo(3) = gami
*
* check for invalid lattice parameters
*
      DO L=1,3
        IF (AO(L).LE.0.0) GO TO 270
        IF (ALO(L).LE.0.0.OR.ALO(L).GE.180.0) GO TO 270
      enddo
      IF (SUM.LT.360.0) GO TO 280
 270  WRITE(NTOUT,13) (AO(L),L=1,3),(alo(L),L=1,3)
 13    FORMAT('0', '**ERROR** Given parameters in error',3F11.5,
     * 2X,3F7.3, ' set skipped')
      GO TO 100
*
* prepare to calculate reciprocal cell parameters
*
 280  IMN=((N+3)*N)/2
      DO L=1,3
        AD(L)=AO(L)
        ALS(L)=ALO(L)
        ARAD=0.01745329*ALO(L)
        CSA(L)=COS(ARAD)
        SNA(L)=SIN(ARAD)
      enddo
*      SUBROUTINE RSCMP
C COMPUTES RECIPROCAL CELL ELEMENTS FROM DIRECT CELL ELEMENTS
*
      DO I=1,3
        J=MOD(I,3)+1
        K=MOD(J,3)+1
        U(I)=(CSA(J)*CSA(K)-CSA(I))/(SNA(J)*SNA(K))
        W(I)=SQRT(1.0-U(I)**2)
        X(J)=1.0/(AD(J)*W(I)*SNA(K))
        AR(I)=57.29578*(1.5707963-ATAN2(U(I),W(I)))
      enddo
      VR=X(1)*X(2)*X(3)*SNA(1)*W(2)*W(3)
      VC=1.0/VR
C-----CALCULATE (A*)**2,ETC.
      DO I=1,3
        J = MOD(I,3)+1
        K = MOD(J,3)+1
        Z(I)=X(I)**2
      enddo
      DO L=1,3
        UO(L)=U(L)
        XO(L)=X(L)
      enddo
      VRO=VR
      VCO=VC
      WRITE(NTOUT,14)
 14    FORMAT(1H0,24X,1Ha,10X,1Hb,10X,1Hc,10X,5Halpha,5X,4Hbeta,4X,
     * 5Hgamma)
      WRITE(NTOUT,15) (X(I),I=1,3),(ar(I),I=1,3)
15     FORMAT(16H Reciprocal Cell,5x,3(f9.5,2x),3(2x,F7.3))
      WRITE(NTOUT,16) (AD(L),L=1,3),(als(L),L=1,3)
 16    FORMAT (' ','Direct cell',9X,3(f9.5,2x),3(2x,F7.3))
*
C SORT obs REFLECTIONS SO AS TO BE IN DECREASING VALUE FOR D-SPACING
*
      I=NUMOBS
 290  I=I-1
      IF (I.LT.1) GO TO 310
      DO 300 J=1,I
        IF (DSP(J).GE.DSP(J+1)) GO TO 300
        DOBS=DSP(J)
        DSP(J)=DSP(J+1)
        DSP(J+1)=DOBS
        THOBS=OBSTH(J)
        OBSTH(J)=OBSTH(J+1)
        OBSTH(J+1)=THOBS
        WLC=WL(J)
        WL(J)=WL(J+1)
        WL(J+1)=WLC
        SWT=SW(J)
        SW(J)=SW(J+1)
        SW(J+1)=SWT
        DO K=1,4
          ITEMP=IH(J,K)
          IH(J,K)=IH(J+1,K)
          IH(J+1,K)=ITEMP
        enddo
 300  CONTINUE
      GO TO 290
 310  IF(NUMOBS.EQ.0) GOTO 320
C
C-----CALC THEMX AND DMIN
C-----ADD 0.5 DEG THETA TO OBSTH(MAX) TO ASSURE THAT ALL HKL ARE GENERAT
C
      THEMO =OBSTH(NUMOBS)+0.5
      THEMX =AMAX1(THEMX,THEMO)
 320  IF(THEMX.LT.0.01 .OR. THEMX.GT.90.0) THEMX=89.0
      DMIN = 0.5*PWL/SIN(THEMX/57.296)
C
C-----CALCULATE THE APPROXIMATE NUMBER OF REFLECTIONS TO BE GENERATED
C-----FOR THIS LATTICE AND CRYSTAL SYSTEM.
C-----# IN SPHERE = (4/3)*PI*VC/DMIN**3
C----- FUDGE FACTOR  IS USED TO COMPENSATE FOR ZONE REFLECTIONS
C-----MULTPL IS THE GENERAL REFLECTION MULTIPLICITY FOR THE CRYSTAL SYSTEM
C
C     ANORTHIC = TRICLINIC
       MULTPL = 2
C-----MONOCLINIC
      IF(ICRS.EQ.ITCRS(2)) MULTPL = 4
C-----ORTHORHOMBIC
      IF(ICRS.EQ.ITCRS(3)) MULTPL = 8
C-----TETRAGONAL
      IF(ICRS.EQ.ITCRS(4)) MULTPL = 16
C-----HEXAGONAL
      IF(ICRS.EQ.ITCRS(5)) MULTPL = 24
C-----RHOMBOHEDRAL
      IF(ICRS.EQ.ITCRS(6).AND.IRH.EQ.1) MULTPL = 24
      IF(ICRS.EQ.ITCRS(6).AND.IRH.EQ.0) MULTPL =  8
C-----CUBIC
      IF(ICRS.EQ.ITCRS(7)) MULTPL = 48
C-----MODIFY MULTPL FOR CENTERING
      ICENT=SPGI(1:1)
C-----A, B, C, I-CENTERING
      IF(ICENT.EQ.ITSG(2).OR.ICENT.EQ.ITSG(3).OR.ICENT.EQ.ITSG(4).OR.
     * ICENT.EQ.ITSG(6)) MULTPL = 2*MULTPL
C-----F CENTERING
      IF(ICENT.EQ.ITSG(5)) MULTPL = 4*MULTPL
      if(qmodul) multpl = 3*multpl
      DLIM = (1.333*3.14156*VC/(MULTPL*MAXGEN))**0.333
      IF(DMIN.GT.DLIM) GOTO 330
      WRITE(NTOUT,17) MAXGEN
 17    FORMAT('0','**MESSAGE** d/hkl-generation limited to',I5,
     *' reflections'/' - 2theta(max) decreased' )
      DMIN = DLIM
      THEMX = THETA(0.5*PWL,DLIM)
 330  ANGLE=THEMX+THEMX
      IF (IDOLSQ .EQ. 1) WRITE(NTOUT,18) THTMX,
     * NCYC,TOLMN,TOLMX,THEMX,DMIN
 18   FORMAT(24H01-theta angles  THTMX =,F6.1,8H  NCYC =,I2,9H  TOLMN =,
     * F7.4,9H  TOLMX =,F7.4/ 17x,7HTHEMX =,F6.1,8H  DMIN =,F10.6)
      IF((IXPD.EQ.0) .AND. (IDOLSQ .EQ. 0)) WRITE(NTOUT,19) ANGLE,DMIN
 19    FORMAT(1H0,10X,13H2-theta(max)=,F6.2,4X,7Hd(min)=,F9.6)
*
* write conditions for non-extinction if full output (ixpd=1)
*
      IF(IXPD.EQ.0) WRITE(NTOUT,20) NOPT
 20    FORMAT(1H0,10X,I2,' Conditions for non-extinction called for ',
     *'class condition(s)')
      IF(NOPT .LT. 1) GO TO 360
      DO 350  L=1,NOPT
        J= JOP(L)
        I= JCLASS(J)
        K= JCOND(J)
        IF(K.GE.1 .AND. K.LE.10) GOTO 340
        IF(IXPD.EQ.0) WRITE(NTOUT,21) IOP(J), IOP(J), IOP(J)
 21     FORMAT(1H ,21X,12Hhkl     h+k=,I1,6Hn,h+l=,I1,6Hn,k+l=,I1,1Hn)
        GOTO 350
 340   IF(IXPD.EQ.0) WRITE(NTOUT,22) ICLASS(I),ICOND(K),IOP(J)
 22     FORMAT(1H ,20X,A4,1X,A6,1H=,I1,1Hn)
 350  CONTINUE
 360  CONTINUE
*
C GENERATE MILLER INDICES CORRESPONDING TO PLANES AT DISTANCES FROM
C DMIN ON UP - MAXIMUM OF MAXGEN REFLECTIONS MAY BE GENERATED
C                                !values given for orthogonal systems
 370  QMAX=1.0/(DMIN*DMIN)             !Qmax = d*2max
      NR=0                             !number of reflection
      E(1)=Z(1)                        !a*2
      E(2)=Z(2)                        !b*2
      E(3)=Z(5)                        !a*b*Cos(gamma) = 0 for ortho
      TD=E(1)*E(2)-E(3)*E(3)           !a*2 b*2 for ortho
      IL=0                             !set initial l to zero
 380  NH(3)=IL                         !set current l value
      E(4)=(Z(6)+Z(6))*NH(3)           !zero
      E(5)=(Z(7)+Z(7))*NH(3)           !zero
      E(6)=Z(3)*NH(3)*NH(3)-QMAX       !d*2(00l) - d*2(max)
      TN=E(3)*E(5)-E(1)*E(4)           !zero
      DSCRS=TN*TN-(4.0*E(1)*E(6)-E(5)*E(5))*TD  !4a*4 b*2(d*2(max)-d*2(00l))
      IF (DSCRS.LT.0.0) GO TO 580      !d*2(00l) exceeds d*2(max)
*
*      SUBROUTINE SKIP
*
      DSCRM=SQRT(DSCRS)                !2 a*2 b*(d*2(max)-d*2(00l))**1/2
      IHKSW = 0                        !I'm not sure about this yet
C-----LIMITS FOR K IN LAYER OF FIXED L
      IK1=IFIX((TN-DSCRM)*0.5/TD)      !int(d*2(00l)-d*2(max))**1/2
      IK2=IFIX((TN+DSCRM)*0.5/TD)      !int(d*2(max)-d*2(00l))**1/2
C-----IF A-UNIQUE MONOCLINIC: NO RESTRICTIONS ON K
      IF(JTYP.EQ.4 .AND. IFA.EQ.0) GO TO 410  
      GO TO (390,400,400,400,400,410,410,400),JTYP  !jtyp=3 for ortho
 390  IF (IK1.LT.   IL) IK1=   IL       !for cubic
      GO TO 410 
 400  IF (IK1.LT.0) IK1=0               !hex, ortho, tetra, mono
 410  IK=IK2                            !a-unique mono
      IF (IK .LT. IK1) GO TO 570        !d*2(00l) too large
C-----NEXT K, IS IT WITHIN LIMITS?
 420  NH(2)=IK                          !assign current k to nh(array)
      GT=(E(3)+E(3))*NH(2)+E(5)         !zero
      HT=(Z(6)*NH(3)+Z(2)*NH(2))*NH(2)+(Z(6)*NH(2)+Z(3)*NH(3))*NH(3)-
     * QMAX                             !d*2(0kl) - d*2(max)
*
*      SUBROUTINE SKIP1
*
      DEN = Z(1)                        !a*2
      DEN = DEN + Z(1)                  !2 a*2
      DSCRS = 2.0 * DEN * HT            !4 a*2(d*2(0kl)-d*2(max))
      DTMP = GT*GT                      !zero
      DSCRS = DTMP - DSCRS              !-4 a*2(d*2(0kl)-d*2(max))
      IF (DSCRS.LT.0.0) GO TO 450       !d*2(0kl) exceeds d*2(max)
      DSCRM=SQRT(DSCRS)                 !2 a*(d*2(0kl)-d*2(max))**1/2
      KSWT = 2
      goto 460
 450  KSWT = 1
      GO TO 460
 460  GO TO (560,470),KSWT
C-----LIMITS FOR H WITH K AND L FIXED
 470  IH1=IFIX(-(GT+DSCRM)/DEN)         !int(d*2(0kl)-d*2(max))**1/2
      IH2=IFIX((DSCRM-GT)/DEN)          !int(d*2(max)-d*2(0kl))**1/2
C-----IF A-UNIQUE MONOCLINIC: LIMIT H TO BE .GE.0
      IF(JTYP.EQ.4 .AND. IFA.EQ.0) GO TO 490
      GO TO (500,500,490,510,500,510,510,480),JTYP
 480  IF (MOD(IL,3).EQ.0) GO TO 500
 490  IF (IH1.LT.0) IH1=0
      GO TO 510
 500  IF (IH1.LT.IK) IH1=IK
 510  jH=IH2
      IF (jH .LT. IH1) GO TO 560         !d*2(0kl) exceeds d*2(max)
C-----NEXT H. IT IS WITHIN D-SPACING LIMITS.
 515  NH(1)=jH                           !assign h to nh(array)
      if(.not.qmodul)goto 518
      im1 = mmin
      im2 = mmax
      im = im1
      if(im.gt.im2)go to 550
 516  nh(4) = im
C-----TEST HKL FOR SPACE GROUP EXTINCTIONS
C-----PATCH ADDED BY CRH TO PROPERLY HANDEL PA3 HKL GENERATION.
C-----CODE ORIGINALLY ASSUMED PERMUTABLE INDICIES-BUT THIS IS NOT
C-----ALLOWED IN PA3. FOR EXAMPLE, THE EXTINCT 320 WAS GENERATED
C-----BUT THE ALLOWED 230 WAS NOT.
 518  IF (NOPT.LT.1) GO TO 530        !no conditions for nonextinction
      IHKSW = 0
      IHT = jH
      IKT = IK
 520  CALL CHKEXT(IHT,IKT,IL,JEXT)    !jext = 0 for nonextinction
      if(nh(1).eq.0.and.nh(2).eq.0.and.nh(3).eq.0)goto 540
      IF (JEXT .EQ. 0) GO TO 530
C-----SPGIN = '205 ' IS TEST FOR SPACE GROUP Pa3
      IF(SPGIN.NE.IH205 .OR. IHKSW.EQ.1) GOTO 540
C-----SWITCH H AND K AND TEST AGAIN FOR PA3
      NH(1) = IK
      NH(2) = jH
      IHT = IK
      IKT = jH
      IHKSW = 1
      GO TO 520
 530  CALL QCOMP(NH,Z,x,u,Q)         !calculate value of d*2(hklm)
      IF (Q.LE.0.0) GO TO 540
C-----TEST FOR MINIMUM D-VALUE
      DTEMP=1.0/SQRT(Q)
      IF(DTEMP.LT.DMIN) GOTO 540
      NR=NR+1                        !we have a new reflection
      DCM(NR)=DTEMP                  !assign new calc d to dcm(array)
C-----ASSIGN HKL TO ARRAY IHG OF COMMON/HDGEN/
      DO J = 1, 4                    !assign hklm values to ihg(array)
        IHG(NR,J) = NH(J)
      enddo
C-----HAVE WE REACHED ARRAY STORAGE CAPACITY?
      IF (NR.LT.MAXGEN) GO TO 540
      CALL SORTD(NR,DCM,IHG)
      CALL SKWEZE
      WRITE(NTOUT,23) MAXGEN,DMIN
C
C   FOR EQUALIVALANCE WITH AIDS
C
      IWARN = IWARN + 1
 23    FORMAT (' ','More than',I4,' reflections exist for dmin =',F10.6)
      DMIN=1.08*DMIN
      WRITE (NTOUT,24) DMIN
 24   FORMAT(31H0New dmin has been set equal to  F10.5/60H0And generatio
     1n of indices started again from the beginning.)
      JSWT = 1
      GO TO 575
 540  if(qmodul)then 
        IF(IHKSW.EQ.0) GOTO 545
      else
        IF(IHKSW.EQ.0) GOTO 550
      endif
C-----SWITCH H(1) AND H(2)  Pa3 PATCH
      NH(1) = IKT
      NH(2) = IHT
      if(.not.qmodul)goto 550
 545  im=im+1
      if(im2-im)550,516,516
 550  jH=jH-1
      IF(jH-IH1)560,515,515
 560  IK=IK-1
      IF(IK-IK1)570,420,420
 570  IL=IL+1
      JSWT = 2
 575  GO TO (370,380),JSWT
 580  CALL SORTD(NR,DCM,IHG)
      CALL SKWEZE
      TOL=TOLMX
      ICYCLE=0
      IF(IDOLSQ.EQ.0) GOTO 840
 590  CALL sums
      call simul
*
* write the results from simul
*
*      write(ntout,82)deter
* 82    format(' The determinant is:',f12.4)
*      do i=1,ntot
*        write(ntout,83)(ff(i,j),j=1,ntot+1)
*      enddo
* 82    format(' the inverse matrix is' / 1h )
* 83    format(1h, 7f9.6)
      IF (NUMOBS .LT. ntot) GO TO 610
      MMN=NUMCYC-ntot                !calc the degrees of freedom
      IF (MMN) 600,640,620
 600  IF ((ICYCLE.GE.10).OR.(ICYCLE.GT.NCYC.AND.TOL.LE.TOLMN))GO TO 610
      TOL=0.5*TOL
      IF (TOL .LT. TOLMN) TOL=TOLMN
      NCYC=0
      GO TO 590
 610  WRITE(NTOUT,25)
 25    FORMAT(12H0No solution)
      GO TO 840
 620  J=N+1
 640  IF (abs(deter).gt.0.01) GO TO 650
      WRITE(NTOUT,26)
 26    FORMAT(' Singular ff(i,j) matrix - rows 
     * are linearly dependent')
      GO TO 600
C
C ASSIGNS CORRECTIONS TO RECIPROCAL CELL ELEMENTS
*
 650  GO TO (710,720,730,740,720,760,770),ITYP
 710  DO J=1,3                         !cubic
        CR(J)=delta(1)
        CR(J+3)=0.0
      enddo
      GO TO 780
 720  DO J=1,2                         !hexagonal/tetragonal
        CR(J)=delta(1)
      enddo
      CR(3)=delta(2)
      DO J=4,6
        CR(J)=0.0
      enddo
      GO TO 780
 730  DO J=1,3                         !orthorhombic
        CR(J)=delta(J)
        CR(J+3)=0.0
      enddo
      GO TO 780
 740  DO J=1,3                         !monoclinic
        CR(J)=delta(J)
      enddo
      CR(4)=0.0
      CR(5)=delta(4)
      CR(6)=0.0
      GO TO 780
 760  DO J=1,6                         !triclinic
        CR(J)=delta(J)
      enddo
      GO TO 780
 770  DO J=1,3                         !trigonal(rhombohedral)
        CR(J)=delta(1)
        CR(J+3)=delta(2)
      enddo
 780  do j=1,nq
        cr(iq(j)) = delta(j+n)
      enddo
*
* make corrections to reciprocal parameters
*
      DO 800 J=1,3
        X(J)=X(J)+CR(J)
        U(J)=U(J)+CR(J+3)
        qj(j)=qj(j)+cr(j+6)
        IF (X(J)) 660,660,670
 660    WRITE(NTOUT,27) J
 27     FORMAT(40H0Corrected value of reciprocal cell edge,I2,
     *  27H less than or equal to zero)
        GO TO 690
 670    IF (ABS(U(J))-1.0) 800,680,680
 680    WRITE(NTOUT,28) J
 28     FORMAT(42H0Corrected value of reciprocal cell cosine,I2,
     *  40H equals or exceeds 1.0 in absolute value)
 690    DO K=1,J
          X(K)=X(K)-CR(K)
          U(K)=U(K)-CR(K+3)
          qj(k)=qj(k)-cr(k+6)
        enddo
        GO TO 600
 800  CONTINUE
C
C-----COMPUTE NEW DIRECT CELL PARAMETERS
C
      CALL CMPNUP
* 
* calculate the standard deviation of the regression
*
 810  IF (SUM .LT. 0.0) SUM=0.0
      if(mmn.eq.0)go to 820
      thtse = sqrt(sum/float(mmn))
 820  IF (TOL .LT. TOLMN) TOL=TOLMN
      IF (TOL .GT. TOLMX) TOL=TOLMX
      IF (ICYCLE .LE. NCYC) GO TO 860
*
* check for convergence if largest correction is < 5.0E-9
*
      CRMAX=0.0
      do i=1,ntot
        if(abs(delta(i)).gt.crmax) crmax=delta(i)
      enddo
*      doft = real(mmn)/real(numobs)
      print*,'CRMAX is equal to: ',crmax
      IF (CRMAX.GE.0.5E-8)GO TO 850
*
* if convergence achieved
*
 830  if(mmn.gt.0) call std
      done = .false.
      CALL OUTPUT(DONE)
      DONE = .TRUE.
 840  CALL FNLST(done)
      if(idolsq.eq.1) call output(done)
      IF(IDOLSQ.EQ.1) CALL COMPAR
      WRITE(NTOUT,29) IWARN,IERR 
 29    FORMAT(1H0, '*****Total number of warnings=',I3, ' *****'/
     *       1H , '*****Total number of errors  =',I3, ' *****')
      IF(IERR.NE.0) NERRT = NERRT + 1
      write(6,35)
 35    format(//' Would you like to run another data set? <no>')
      read(ntin,60,err=100)junk
 60    format(a)
      if(junk(1:1).eq.'n'.or.junk(1:1).eq.' ')go to 880
      GO TO 100
*
* if convergence has not yet been achieved
*
  850 IF (ICYCLE .LT. 10) GO TO 860
      WRITE(NTOUT,30)
 30    FORMAT(' Not convergent in 10 cycles,'/
     *' **MESSAGE** Final list uses cell parameters from cycle 9')
      GO TO 830
 860  CONTINUE
      DONE = .FALSE.
      CALL OUTPUT(DONE)
      CALL RCMPD
      CALL SORTD(NR,DCM,IHG)
      GO TO 590
C
C-----ERROR IN OPENING EXTENS FILE
C
 870  CONTINUE
      WRITE (NTOUT,31) NEXTN
 31    FORMAT (' ', '**ERROR** Error occured in assigning extinction',
     * ' file to unit ', I3, ' execution stopped**')
C
C-----WRITE SUMMARY OF ENTIRE RUN IN TERMS OF DECKS WITH ERRORS
C
 880  WRITE(NTOUT,32) NERRT
 32    FORMAT(1H0,'##### Total decks with errors =',I3,' #####')
 999  STOP
      END
