C     TREOR  VERSION PEW 84 10 02
      COMMON T1(100),TAGUT(100),N1,D1,SSQTL,D2,WAVE,IHM,IKM,ILM,ISYM,INR
     1,ISUB,A11,A22,A33,A13,A12,A23,LIN(20),REV(20)
     2,IGOT,IS1,IS2,IS3,IS4,IA,IB,IQ,IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3
     3,SQ1,SQ2,SQ3,X,Y,Z,IC(20),ISQM,KOTA(30),RO(7,30),ID,IIN,IOUT
      COMMON NH1,NK1,NL1,NS1,NH2,NK2,NL2,NS2,NH3,NK3,NL3,NS3,IUTSK,IMER,
     1INIX,IOMIT,ISC,KLART,VOL,VTEST,XL,YL,IT20,IT40,CEM
     2,KHA,KKA,KLA,KSA,KHB,KKB,KLB,KSB,MONO,MH1,MK1,ML1,MH2,MK2,ML2,MH3,
     3MK3,ML3,MS1,MS2,MS3,MH4,MK4,ML4,MS4,UL,LGO,MG,MS,MRA,MRB,NU,
     4DEN,EDEN,RMOLV
      COMMON VOL1(100),VOL2(100),AA(6,6),BB(6,1),SQQ(5),SLD(5),HL(5),
     1IDX(5),XX(5),SQOLD(5),AAA(5,5),BBB(5,1),VW(5)
      REV(15)=0.  ! USED TIME = 0 SEC.
      CALL GET_USED_CPU_TIME(REV(15))
      REV(16)=0.
      IIN=5
      IOUT=6
      NDISP=0
      LIN(11)=NDISP
c     OPEN(UNIT=IIN,NAME='tredat',STATUS='OLD',READONLY,ERR=11)
c     OPEN(UNIT=IOUT,NAME='list',STATUS='NEW',RECL=130,ERR=12)
c     GO TO 6
c  11 WRITE(NDISP,*) ' OPEN ERROR ON TREOR INPUT DATA FILE '
c     GO TO 5
c  12 WRITE(NDISP,*) ' OPEN ERROR ON PRINTER OUTPUT FILE '
c     GO TO 5
    6 CALL PWINL
    2 CALL TREOB
    3 CALL TREOC(IGOE)
      GO TO(1,2,4),IGOE
    1 CALL TREOD
      IGT=LIN(3)
      GO TO(3,19),IGT
    4 IF(LGO) 19,19,6
   19 IF(MONO .NE. 0) GO TO 17
      IF(LIN(8) .EQ. 0) GO TO 18
   17 WRITE(IOUT,20) IQ,LIN(4),LIN(6),LIN(5),LIN(7)
   20 FORMAT(' NUMBER OF CELLS WITH ',I3,' OR MORE INDEXABLE LINES ',/,
     1' IN MONOCLINIC (020)-TESTS ',I5,' SOLUTIONS ',/,
     2' IN MONOCLINIC SHORT AXIS TESTS ',I5,' SOLUTIONS ',/,
     3' IN MONOCLINIC GENERAL TESTS ',I5,' SOLUTIONS ',/,
     4' IN TRICLINIC TESTS ',I5,' SOLUTIONS ')
   18 CALL GET_USED_CPU_TIME(REV(16))
      TIME=REV(16)-REV(15)
      WRITE(IOUT,15) TIME
   15 FORMAT(' USED CPU-TIME=',F8.0,' SEC.')
    5 STOP
      END
      SUBROUTINE ATNCS(J,K,L)
C-----THIS SUBROUTINE IS USED BY SUBROUTINE ATNC
      DIMENSION ICH(10)
      DATA ICH/'0','1','2','3','4','5','6','7','8','9'/
      L=0
      DO 1 I=1,10
      IF(J-ICH(I))1,2,1
    1 CONTINUE
      L=1
      RETURN
    2 K=I-1
      RETURN
      END
      SUBROUTINE ATNC(IT,RT,LT,IQ,IVAL)
      DIMENSION IVAL(1)
      DATA IPKT,MINUS,IPACE/'.','-',' '/
      N=-1
      NT=0
      INTG=0
      RNTG=0.0
      IMN=-1
      RMN=-1.0
      DO 2 I=1,LT
      IF(IVAL(I)-MINUS)2,1,2
    1 NT=1
      IVAL(I)=IPACE
      GO TO 3
    2 CONTINUE
    3 DO 4 I=1,LT
      IF(IVAL(I)-IPKT)4,5,4
    4 CONTINUE
      IQ=1
      IX=LT+1
      DO 6 I=1,LT
      IX=IX-1
      IF(IVAL(IX)-IPACE)7,6,7
    7 N=N+1
      IZ=10**N
      CALL ATNCS(IVAL(IX),NIX,IER)
      IF(IER)800,800,801
  800 INTG=INTG+IZ*NIX
    6 CONTINUE
      IT=INTG
      IF(NT)900,900,10
   10 IT=IT*IMN
  900 RETURN
    5 IQ=2
      IPOS=I-1
      J=I
      DO 11 I=1,IPOS
      J=J-1
      IF(IVAL(J)-IPACE)12,11,12
   12 N=N+1
      RZ=10.0**N
      CALL ATNCS(IVAL(J),NIX,IER)
      IF(IER)802,802,801
  802 RNTG=RNTG+RZ*FLOAT(NIX)
   11 CONTINUE
      RN1=RNTG
      RN2=0.0
      IPOS=IPOS+2
      RZ=1.0
      DO 15 I=IPOS,LT
      IF(IVAL(I)-IPACE)16,15,16
   16 RZ=RZ*10.0
      CALL ATNCS(IVAL(I),NIX,IER)
      IF(IER)803,803,801
  803 RN2=RN2+FLOAT(NIX)/RZ
   15 CONTINUE
      RT=RN1+RN2
      IF(NT)900,900,17
   17 RT=RT*RMN
      GO TO 900
  801 IQ=3
      GO TO 900
      END
      SUBROUTINE FFCCR(IFILD,NFILD,RFILD,IR,IER,IIN)
C-----FREE FORMAT CONTROL CARD READER, LOADS NUMERIC VALUES FROM CARD
C-----IN THE MODE OF KEYWORD=NUMERIC VALUE, ETC.
C-----VERSION 1974.02.14
      DIMENSION ICARD(80),KEY(10),        NUM(20),KEYNU(31),
     -          IFILD(1),NFILD(1),RFILD(1)
      DATA ISPAC,ICOM,ISLA,ILIK,IAST,IPKT/' ',',','/','=','*','.'/
      DATA I1,I2,I3/'E','N','D'/
      KEYLE=10
      NUMLE=15
      NUKE=NUMLE+KEYLE+1
      JSCOL=80
      IER=0
      JSW=0
    4 NPOS=0
      READ(IIN,200)ICARD
  200 FORMAT(80A1)
      IF(JSW)900,900,52
  900 I=0
   52 N=NPOS
      JSW=0
    2 N=N+1
      IF(N-JSCOL)3,3,40
   40 JSW=1
      GO TO 4
    3 IX=ICARD(N)
      IF(IX-ISPAC)1,2,1
    1 IF(IX-ICOM)5,6,5
    5 IF(IX-IAST)50,6,50
   50 I=I+1
      IF(I-NUKE)7,7,8
    8 IER=1
      GO TO 20
    7 KEYNU(I)=IX
      GO TO 2
    6 NPOS=N
      IPOS=I
      DO 600 I=1,NUMLE
  600 NUM(I)=ISPAC
      DO 601 I=1,KEYLE
  601 KEY(I)=ISPAC
      ISW=0
      J=0
      I=0
    9 I=I+1
      IF(I-IPOS)10,10,11
   10 IX=KEYNU(I)
      IF(ISW)14,14,15
   14 IF(IX-ILIK)12,13,12
   12 IF(I-KEYLE)120,120,8
  120 KEY(I)=IX
      GO TO 9
   13 ISW=1
      IPES=I-1
      GO TO 9
   15 J=J+1
      IF(J-NUMLE)400,400,8
  400 NUM(J)=IX
      GO TO 9
   11 IF(ICARD(NPOS)-IAST)17,111,17
  111 IF(KEY(1)-I1)17,18,17
   18 IF(KEY(2)-I2)17,19,17
   19 IF(KEY(3)-I3)17,20,17
   20 I=0
   24 I=I+1
      IX=IFILD(I)
      IF(IX-IPKT)21,22,21
   21 IF(IX-ISLA)24,23,24
   23 IFILD(I)=IAST
      GO TO 24
   22 RETURN
   17 NPIS=0
      NUMLA=J
      IA=0
   25 N=NPIS
   26 N=N+1
      IX=IFILD(N)
      IF(IX-IPKT)261,32,261
  261 IF(IX-IAST)260,27,260
  260 IF(IX-ISLA)26,270,26
  270 IA=IA+1
      GO TO 26
   27 NPUS=N+1
      IA=IA+1
   31 N=N+1
      IX=IFILD(N)
      IF(IX-IAST)28,29,28
   28 IF(IX-ISLA)30,29,30
   30 IF(IX-IPKT)31,32,31
   32 IER=2
      GO TO 20
   29 NPAS=N-1
      IA=IA+1
      NPIS=N-1
      IF(((NPAS-NPUS)+1)-IPES)250,33,250
   33 N=0
      DO 34 I=NPUS,NPAS
      N=N+1
      IF(IFILD(I)-KEY(N))250,34,250
   34 CONTINUE
      IA=IA-1
      GO TO 311
  250 IA=IA-1
      GO TO 25
  311 N=NPUS-1
      IFILD(N)=ISLA
      CALL ATNC(IT,RT,NUMLE,IQ,NUM)
      IF(IQ-2)326,327,328
  328 IER=3
      GO TO 20
  327 IF(IA-IR)320,321,321
  320 NFILD(IA)=IFIX(RT)
      IER=4
      GO TO 900
  321 IA=(IA-IR)+1
      RFILD(IA)=RT
      GO TO 900
  326 IF(IA-IR)322,323,323
  323 IA=(IA-IR)+1
      RFILD(IA)=FLOAT(IT)
      IER=5
      GO TO 900
  322 NFILD(IA)=IT
      GO TO 900
      END
      SUBROUTINE PWINL
      DIMENSION INTG(54),RNTG(9),KEY(292),KEA(61),KEB(63),KEC(16),
     -KED(55),KEE(54),KEF(43)
C
      COMMON T1(100),TAGUT(100),N1,D1,SSQTL,D2,WAVE,IHM,IKM,ILM,ISYM,INR
     1,ISUB,A11,A22,A33,A13,A12,A23,LIN(20),REV(20)
     2,IGOT,IS1,IS2,IS3,IS4,IA,IB,IQ,IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3
     3,SQ1,SQ2,SQ3,X,Y,Z,IC(20),ISQM,KOTA(30),RO(7,30),ID,IIN,IOUT
      COMMON NH1,NK1,NL1,NS1,NH2,NK2,NL2,NS2,NH3,NK3,NL3,NS3,IUTSK,IMER,
     1INIX,IOMIT,ISC,KLART,VOL,VTEST,XL,YL,IT20,IT40,CEM
     2,KHA,KKA,KLA,KSA,KHB,KKB,KLB,KSB,MONO,MH1,MK1,ML1,MH2,MK2,ML2,MH3,
     3MK3,ML3,MS1,MS2,MS3,MH4,MK4,ML4,MS4,UL,LGO,MG,MS,MRA,MRB,NU,
     4DEN,EDEN,RMOLV
      COMMON VOL1(100),VOL2(100),AA(6,6),BB(6,1),SQQ(5),SLD(5),HL(5),
     1IDX(5),XX(5),SQOLD(5),AAA(5,5),BBB(5,1),VW(5)
C
      EQUIVALENCE(KEY(1),KEA(1)),(KEY(62),KEB(1)),(KEY(125),KEC(1)),
     -(KEY(141),KED(1)),(KEY(196),KEE(1)),(KEY(250),KEF(1))
c
      character name*79
C
      DATA KEA/'*','O','H','1','*','O','K','1','*','O','L','1','*','O','
     -S','1','*','O','H','2','*','O','K','2','*','O','L','2','*','O','S'
     -,'2','*','O','H','3','*','O','K','3','*','O','L','3','*','O','S','
     -3','*','K','H','*','K','K','*','K','L','*','K','S','*'/
      DATA KEB/'T','H','H','*','T','H','K','*','T','H','L','*','T','H','
     -S','*','M','E','R','I','T','*','N','I','X','*','L','I','S','T',
     -'*','O','M','I','T','*','S','E','L','E','C','T','*','U','S','E','*
     -','I','Q','*','M','O','N','O','*','M','H','1','*','M','K','1','*'/
      DATA KEC/
     -'M','L','1','*','M','H','2','*','M','K','2','*','M','L','2','*'/
      DATA KED/'M','H','3','*','M','K','3','*','M','L','3','*','M','S',
     -'1','*','M','S','2','*','M','S','3','*','M','H','4','*','M','K',
     -'4','*','M','L','4','*','M','S','4','*','C','H','O','I','C','E',
     -'*','C','O','N','T','I','N','U','E'/
      DATA KEE/
     -'*','M','O','N','O','G','A','M','*','M','O','N','O','S','E','T','*
     -','M','R','E','S','O','N','E','*','M','R','E','S','T','W','O',
     -'*','V','C','R','I','T','*','I','D','I','V','*','S','H','O','R',
     -'T','*','T','R','I','C'/
      DATA KEF/
     -'*','D','1','*','S','S','Q','T','L','*','D','2','*','C','E','M',
     -'*','V','O','L','*','W','A','V','E','*','D','E','N','S',
     -'*','E','D','E','N','S','*','M','O','L','W','*','.'/
C
      NDISP=LIN(11)
      ID=0
      A11=0.0
      A22=0.0
      A33=0.0
      A13=0.0
      A12=0.0
      A23=0.0
      ISUB=0
      INR=3
      IGOT=0
      IREAL=55
C
      WRITE(IOUT,1111)
 1111 FORMAT(20X,'TREOR (4)- 84 10 02')
c     READ(IIN,1)
c   1 FORMAT(1X,79H N A M E                                      N A M E
c    1                   N A M E)
c     WRITE(IOUT,1)
      read (iin, '(1x, a)') name
      write(iout, '(a)') name
C
      N1=0
    2 N1=N1+1
      READ(IIN,3)T1(N1),TAGUT(N1)
    3 FORMAT(1X,F15.6,10X,A4)
      IF(T1(N1))400,400,401
  401 WRITE(IOUT,3)T1(N1),TAGUT(N1)
      GO TO 2
C
C     IF THE KEYWORD SYSTEM CAN NOT BE USED ON A PARTICULAR COMPUTER
C     ALL STEERING PARAMETERS CAN BE READ WITH ORDINARY FORMATTED READ.
C     IN THAT CASE, REMOVE SUBROUTINES FFCCR,ATNC AND ATNCS
C     READ IN THE PARAMETERS IN SECTION KEYOUT USING NORMAL FORTRAN
C     FORMATTED CARD INPUT, THESE PARAMETERS CORRESPOND ONE TO ONE TO
C     THE KEYWORDS IN SECTION KEYIN
C
C-----SECTION KEYIN
  400 N1=N1-1
      INTG(1)=2
      INTG(2)=2
      INTG(3)=2
      INTG(4)=3
      INTG(5)=2
      INTG(6)=2
      INTG(7)=2
      INTG(8)=4
      INTG(9)=2
      INTG(10)=2
      INTG(11)=2
      INTG(12)=4
      INTG(13)=4
      INTG(14)=4
      INTG(15)=4
      INTG(16)=6
      INTG(17)=4
      INTG(18)=4
      INTG(19)=4
      INTG(20)=4
      INTG(21)=10
      INTG(22)=1
      INTG(23)=0
      INTG(24)=0
      INTG(25)=0
      IF(N1-19)301,302,302
  301 INTG(26)=N1
      GO TO 303
  302 INTG(26)=19
  303 INTG(27)=INTG(26)-3
      INTG(28)=0
      INTG(29)=2
      INTG(30)=2
      INTG(31)=2
      INTG(32)=2
      INTG(33)=2
      INTG(34)=2
      INTG(35)=2
      INTG(36)=2
      INTG(37)=2
      INTG(38)=2
      INTG(39)=3
      INTG(40)=3
      INTG(41)=2
      INTG(42)=2
      INTG(43)=2
      INTG(44)=4
      INTG(45)=0
      INTG(46)=0
      INTG(47)=1
      INTG(48)=0
      INTG(49)=0
      INTG(50)=0
      INTG(51)=0
      INTG(52)=1
      INTG(53)=1
      INTG(54)=0
      RNTG(1)=0.0002
      RNTG(2)=0.05
      RNTG(3)=0.0004
      RNTG(4)=25.0
      RNTG(5)=2000.0
      RNTG(6)=1.5405981
      RNTG(7)=0.0
      RNTG(8)=0.2
      RNTG(9)=0.0
C
      CALL FFCCR(KEY,INTG,RNTG,IREAL,IERR,IIN)
C
      IF(IERR)9001,9001,9002
 9002 IF(IERR-3)9000,9000,9001
 9000 WRITE(IOUT,300)
  300 FORMAT(//,1X,39H **ERROR** PROGRAM CONTROL CARD INVALID,/,1X,68H P
     @ROPABLY MISSING COMMA OR DELIMITER KEYWORD OR UNIDENTIFIED KEYWORD
     @)
      STOP
C-----END KEYIN
C
C-----SECTION KEYOUT
 9001 CONTINUE
      NH1=INTG(1)
      NK1=INTG(2)
      NL1=INTG(3)
      NS1=INTG(4)
      NH2=INTG(5)
      NK2=INTG(6)
      NL2=INTG(7)
      NS2=INTG(8)
      NH3=INTG(9)
      NK3=INTG(10)
      NL3=INTG(11)
      NS3=INTG(12)
      KHB=INTG(13)
      KKB=INTG(14)
      KLB=INTG(15)
      KSB=INTG(16)
      KHA=INTG(17)
      KKA=INTG(18)
      KLA=INTG(19)
      KSA=INTG(20)
      IMER=INTG(21)
      INIX=INTG(22)
      IUTSK=INTG(23)
      IOMIT=INTG(24)
      ISC=INTG(25)
      IB=INTG(26)
      IQ=INTG(27)
      MONO=INTG(28)
      MH1=INTG(29)
      MK1=INTG(30)
      ML1=INTG(31)
      MH2=INTG(32)
      MK2=INTG(33)
      ML2=INTG(34)
      MH3=INTG(35)
      MK3=INTG(36)
      ML3=INTG(37)
      MS1=INTG(38)
      MS2=INTG(39)
      MS3=INTG(40)
      MH4=INTG(41)
      MK4=INTG(42)
      ML4=INTG(43)
      MS4=INTG(44)
      LGO=INTG(46)
      MG=INTG(47)
      MS=INTG(48)
      MRA=INTG(49)
      MRB=INTG(50)
      IVCR=INTG(51)
      IDIV=INTG(52)
      ISHORT=INTG(53)
      LIN(9)=ISHORT
      ITRIC=INTG(54)
      LIN(8)=ITRIC
      D1=RNTG(1)
      SSQTL=RNTG(2)
      D2=RNTG(3)
      CEM=RNTG(4)
      VOL=RNTG(5)
      WAVE=RNTG(6)
      DEN=RNTG(7)
      EDEN=RNTG(8)
      RMOLV=RNTG(9)
      LIN(14)=0
      LIN(15)=0
C     LIN(4-7)=COUNTERS FOR TRIAL CELLS
      LIN(4)=0
      LIN(5)=0
      LIN(6)=0
      LIN(7)=0
      LIN(16)=0
C-----END KEYOUT
C
      ULA=0.0174533*FLOAT(MONO)
      ULB=COS(ULA)
C     4*(COS(BETA MAX))**2
      UL=4.0*ULB*ULB
C     CHOICE=0  SQ(THETA), =1 1/D2, =2 THETA, =3 2THETA, =4 D
      ICHOI=INTG(45)+1
C
      GO TO(101,102,103,104,105),ICHOI
  102 OMV=(WAVE*WAVE)/4.0
      DO 110 I=1,N1
      T1(I)=OMV*T1(I)
  110 CONTINUE
      GO TO 101
  103 OMV=0.01745329
  112 DO 111 I=1,N1
      T1(I)=SIN(OMV*T1(I))**2
  111 CONTINUE
      GO TO 101
  104 OMV=0.008726646
      GO TO 112
  105 OMV=(WAVE*WAVE)/4.0
      DO 113 I=1,N1
      H=T1(I)
      T1(I)=OMV/(H*H)
  113 CONTINUE
C
  101 IF(IUTSK) 40,42,40
   40 WRITE(IOUT,41) IQ
   41 FORMAT(28H PRINT ALL TRIALS WITH IQ GE,I5)
   42 CALL PDIV(T1,N1,D1,SSQTL,D2,IDIV)
      WRITE(IOUT,44)
   44 FORMAT(12H STOP LIMITS)
      WRITE(IOUT,43) IMER,INIX
   43 FORMAT(26H FIGURE OF MERIT REQUIRED=,I5,31H MAX NUMBER OF UNINDEXE
     1D LINES=,I5)
      IF(IDIV) 501,502,501
  501 WRITE(IOUT,503)
  503 FORMAT(50H THE 7 FIRST LINES ADJUSTED BY THEIR HIGHER ORDERS)
  502 CONTINUE
      IF(IOMIT) 47,45,47
   45 WRITE(IOUT,46)
   46 FORMAT(52H CUBIC,TETRAGONAL,HEXAGONAL AND ORTHOROMBIC SYMMETRY)
      GO TO 49
   47 WRITE(IOUT,48)
   48 FORMAT(21H ORTHOROMBIC SYMMETRY)
   49 WRITE(IOUT,64) CEM,VOL,D1,SSQTL,D2,WAVE
   64 FORMAT(15H MAX CELL EDGE=,F5.1,17H MAX CELL VOLUME=,F10.1,/,4H D1=
     1,F10.6,7H SSQTL=,F10.6,4H D2=,F10.6,6H WAVE=,F10.6)
      IF(DEN-0.00001)201,201,200
  200 WRITE(IOUT,50)RMOLV,DEN,EDEN
   50 FORMAT(1X,'MOLAR WEIGHT=',F10.3,' DENSITY=',F10.3,' MAX. DENSITY D
     *EVIATION=',F10.3)
      IZTAL=0
  350 IZTAL=IZTAL+1
      RVA=(IZTAL*RMOLV)/(0.6023*(DEN+EDEN))
      IF(RVA-VOL) 350,350,351
  351 RVA=(WAVE**6)*((0.6023*(DEN+EDEN))**2)
      RVB=(WAVE**6)*((0.6023*(DEN-EDEN))**2)
      DO 352 I=1,IZTAL
      VOL1(I)=RVB/(64.0*(I*RMOLV)**2)
      VOL2(I)=RVA/(64.0*(I*RMOLV)**2)
  352 CONTINUE
      LIN(14)=IZTAL
C
  201 IT20=100000.0*D1+0.5
      IT40=100000.0*D2+0.5
      XL=(WAVE**2)/(4.0*CEM*CEM)
      YL=(4.0*XL)/3.0
      VSK=(WAVE**6)/64.0
      VTEST=VSK/(VOL**2)
      DO 250 I=1,4
      VW(I)=VSK/((I*500.)**2)
  250 CONTINUE
      VW(5)=VSK/(VOL**2)
      REV(11)=VTEST
      REV(12)=XL
      DO 71 I=1,IB
      HES=T1(I)
      IF(HES-0.327) 70,70,68
C
   68 WRITE(IOUT,69) IB
   69 FORMAT(53H NUMBER OF LINES WITH SQ LESS THAN 0.327 NOT EQUAL TO,I5
     1)
      GO TO 76
   70 IC(I)=HES*100000.0+0.5
   71 CONTINUE
      GO TO 205
   76 IB=I-1
      IQ=IB-3
  205 WRITE(IOUT,39) IB,IQ
   39 FORMAT(22H NUMBER OF TEST LINES=,I5,13H IQ REQUIRED=,I5)
      ISQM=IC(IB)+60
      IF(IOMIT) 73,72,73
   72 ISYM=1
      RETURN
   73 ISYM=3
      RETURN
      END
      SUBROUTINE EKVA(IH1,IK1,IL1,IH2,IK2,IL2,ISYM,ISOLV,SKVA,SKVB,X,Y,X
     1L,YL)
C     ISYM=1 KUB  ISYM=2 TETR  ISYM=4 HEXAG.
      GO TO(1,2,3,4,3),ISYM
    1 IK=IH1*IH1+IK1*IK1+IL1*IL1
      R=FLOAT(IK)
      X=SKVA/R
      IF(X-XL)5,5,6
    5 ISOLV=0
      RETURN
    6 ISOLV=1
      RETURN
    2 I11=IH1*IH1+IK1*IK1
      I21=IH2*IH2+IK2*IK2
   10 I12=IL1*IL1
      I22=IL2*IL2
      IDET=I11*I22-I21*I12
      IF(IDET)7,5,7
    7 DET=FLOAT(IDET)
      R11=FLOAT(I11)
      R21=FLOAT(I21)
      R12=FLOAT(I12)
      R22=FLOAT(I22)
      X=(SKVA*R22-R12*SKVB)/DET
      IF(X-XL)5,5,8
    8 Y=(R11*SKVB-SKVA*R21)/DET
      IF(Y-YL)5,5,6
    4 I11=IH1*(IH1+IK1)+IK1*IK1
      I21=IH2*(IH2+IK2)+IK2*IK2
      GO TO 10
    3 STOP
      END
      SUBROUTINE LAGRA(K)
      COMMON T1(100),TAGUT(100),N1,D1,SSQTL,D2,WAVE,IHM,IKM,ILM,ISYM,INR
     1,ISUB,A11,A22,A33,A13,A12,A23,LIN(20),REV(20)
     2,IGOT,IS1,IS2,IS3,IS4,IA,IB,IQ,IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3
     3,SQ1,SQ2,SQ3,X,Y,Z,IC(20),ISQM,KOTA(30),RO(7,30),ID,IIN,IOUT
      COMMON NH1,NK1,NL1,NS1,NH2,NK2,NL2,NS2,NH3,NK3,NL3,NS3,IUTSK,IMER,
     1INIX,IOMIT,ISC,KLART,VOL,VTEST,XL,YL,IT20,IT40,CEM
     2,KHA,KKA,KLA,KSA,KHB,KKB,KLB,KSB,MONO,MH1,MK1,ML1,MH2,MK2,ML2,MH3,
     3MK3,ML3,MS1,MS2,MS3,MH4,MK4,ML4,MS4,UL,LGO,MG,MS,MRA,MRB,NU,
     4DEN,EDEN,RMOLV
      GO TO(61,62,63,62),ISYM
   61 XA=X
      YA=X
      ZA=X
      GO TO 30
   62 XA=X
      YA=Y
      ZA=X
      GO TO 30
   63 IF(X-Y) 2,3,3
    3 IF(X-Z) 5,4,4
    5 XA=Z
      YA=X
   21 ZA=Y
      GO TO 30
    4 IF(Y-Z) 6,7,7
    7 XA=X
      YA=Y
   22 ZA=Z
      GO TO 30
    6 XA=X
      YA=Z
      GO TO 21
    2 IF(X-Z) 8,9,9
    9 XA=Y
      YA=X
      GO TO 22
    8 IF(Y-Z) 11,12,12
   11 XA=Z
      YA=Y
   23 ZA=X
      GO TO 30
   12 XA=Y
      YA=Z
      GO TO 23
   30 PV=XA*YA*ZA
      IF(ID) 32,31,32
   31 ID=ID+1
   44 I=ID
   47 RO(4,I)=PV
      RO(1,I)=XA
      RO(2,I)=YA
      RO(3,I)=ZA
      KOTA(I)=K
      RETURN
   32 DO 33 I=1,ID
      IF(ABS(RO(1,I)-XA)-0.00005) 34,33,33
   34 IF(ABS(RO(2,I)-YA)-0.00005) 35,33,33
   35 IF(ABS(RO(3,I)-ZA)-0.00005) 36,33,33
   36 RETURN
   33 CONTINUE
      DO 40 I=1,ID
      IF(KOTA(I)-K)42,41,40
   41 IF(RO(4,I)-PV) 42,40,40
   42 IF(I-5) 43,44,99
   40 CONTINUE
      IF(ID-5) 31,36,36
   43 N=4
      IF(ID-5) 60,49,49
   60 ID=ID+1
   49 M=N+1
      KOTA(M)=KOTA(N)
      DO 46 J=1,4
      RO(J,M)=RO(J,N)
   46 CONTINUE
      IF(N-I) 99,47,48
   48 N=N-1
      GO TO 49
   99 STOP
      END
      SUBROUTINE GHKLA(MAXH,MAXK,MAXL,MAXS,MSEL,ISLUT,IH,IK,IL)
C     MSEL=0 H GE K GE L      MSEL EQ NON ZERO H GE K
C     ISLUT=0 END OF GEN.     ISLUT=1 CONTINUE
C     IH,IK,IL GE 0
      ISLUT=1
    1 IF(IH-MAXH)2,5,5
    2 IH=IH+1
      IF(IH-IK)1,3,3
    3 ISUM=IH+IK+IL
      IF(MAXS-ISUM)5,4,4
    4 RETURN
    5 IF(IK-MAXK)6,7,7
    6 IK=IK+1
      IH=IK
      GO TO 3
    7 IF(IL-MAXL)9,8,8
    8 IH=0
      IK=0
      IL=0
      ISLUT=0
      RETURN
    9 IL=IL+1
      IF(MSEL)10,11,10
   10 IH=0
      IK=0
      RETURN
   11 IH=IL
      IK=IL
      GO TO 3
      END
      SUBROUTINE ORT(IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3,ISOLV,SKVA,SKVB
     1,SKVC,X,Y,Z,IO1,IO2,XL)
      IF(IO1) 1,2,1
    1 IO1=0
      IHH1=IH1*IH1
      IKK1=IK1*IK1
      ILL1=IL1*IL1
      RHH1=FLOAT(IHH1)
      RKK1=FLOAT(IKK1)
      RLL1=FLOAT(ILL1)
      HA=RKK1*SKVC
      HB=RLL1*SKVB
      HC=RLL1*SKVC
      HD=RKK1*SKVB
      HE=RHH1*SKVB
      HF=RHH1*SKVC
    2 IF(IO2) 3,4,3
    3 IO2=0
      IHH2=IH2*IH2
      IKK2=IK2*IK2
      ILL2=IL2*IL2
      RHH2=FLOAT(IHH2)
      RKK2=FLOAT(IKK2)
      RLL2=FLOAT(ILL2)
      LH3=IKK1*ILL2-ILL1*IKK2
      LK3=ILL1*IHH2-IHH1*ILL2
      LL3=IHH1*IKK2-IKK1*IHH2
      RA=SKVA*RKK2
      RB=HA*RLL2
      RC=HC*RKK2
      RD=SKVA*RLL2
      RE=HC*RHH2
      RF=HF*RLL2
      RG=SKVA*RHH2
      RH=HF*RKK2
      RI=HA*RHH2
      SA=RB-RC
      SB=HB-RD
      SC=RA-HD
      SD=RE-RF
      SE=RD-HB
      SF=HE-RG
      SG=RH-RI
      SH=HD-RA
      SI=RG-HE
    4 IHH3=IH3*IH3
      IKK3=IK3*IK3
      ILL3=IL3*IL3
      IDET=IHH3*LH3+IKK3*LK3+ILL3*LL3
      IF(IDET) 6,5,6
    5 ISOLV=0
      RETURN
    6 RHH3=FLOAT(IHH3)
      RKK3=FLOAT(IKK3)
      RLL3=FLOAT(ILL3)
      RDET=FLOAT(IDET)
      X=(SA+SB*RKK3+SC*RLL3)/RDET
      IF(X-XL) 5,5,7
    7 Y=(SD+SE*RHH3+SF*RLL3)/RDET
      IF(Y-XL) 5,5,8
    8 Z=(SG+SH*RHH3+SI*RKK3)/RDET
      IF(Z-XL) 5,5,9
    9 ISOLV=1
      RETURN
      END
      SUBROUTINE ORTAL(T1,X,Y,Z,N1,IB,IQ,ISQM,K,IC,IT20,IT40,ITLIM)
      DIMENSION T1(1),IC(1),KST(30)
      NX=100000.0*X+0.5
      NY=100000.0*Y+0.5
      NZ=100000.0*Z+0.5
      IDH=ISQM/NX
      IDK=ISQM/NY
      IDL=ISQM/NZ
      DO 8 I=1,IB
      KST(I)=0
    8 CONTINUE
      IH=0
      IK=0
      IL=0
      IHA=0
      IHB=0
      IHC=0
      IHD=0
      IHH=0
      IKK=0
      ILL=0
C     GENERATE HKL AND ISQC
   15 IF(IHH-IDH) 51,52,52
   51 IH=IH+1
      IHH=IH*IH
      GO TO 57
   52 IH=0
      IHH=0
      IHA=0
      IF(IKK-IDK) 53,54,54
   53 IK=IK+1
      IKK=IK*IK
      GO TO 56
   54 IK=0
      IKK=0
      IHB=0
      IF(ILL-IDL) 55,30,30
   55 IL=IL+1
      ILL=IL*IL
      IHC=ILL*NZ
      GO TO 60
   56 IHB=IKK*NY
   60 IHD=IHB+IHC
      IHE=ISQM-IHD
      IF(IHE) 54,59,59
   57 IHA=IHH*NX
   59 ISQC=IHA+IHD
      IF(ISQM-ISQC) 52,58,58
   58 IF(ISQC-ITLIM) 40,41,41
   40 IT=IT20
      GO TO 42
   41 IT=IT40
   42 DO 19 I=1,IB
      ISQO=IC(I)
      ITQC=ISQO-ISQC
      IF(IABS(ITQC)-IT) 20,20,65
   65 IF(ITQC) 19,19,15
   20 KST(I)=1
   19 CONTINUE
      GO TO 15
   30 K=0
      DO 31 I=1,IB
      K=K+KST(I)
   31 CONTINUE
      RETURN
      END
      SUBROUTINE GEHKL(MAXH,MAXK,MAXL,MAXS,MORMO,ISLUT,IH,IK,IL)
C     MORMO=1 ORTHO.    MORMO=2 MONOCLINIC. +-H
C     ISLUT=0 END OF GEN. ISLUT=1 CONTINUE
      ISLUT=1
      GO TO(1,2),MORMO
    1 IF(IH-MAXH) 3,4,4
    3 IH=IH+1
      ISUM=IH+IK+IL
      IF(MAXS-ISUM) 4,5,5
    5 RETURN
    4 IH=0
      IF(IK-MAXK) 6,7,7
    6 IK=IK+1
      ISUM=IH+IK+IL
      IF(MAXS-ISUM) 7,5,5
    7 IK=0
      IF(IL-MAXL) 8,9,9
    8 IL=IL+1
      RETURN
    9 IL=0
      ISLUT=0
      RETURN
    2 IF(IH*IL) 10,1,11
   11 IH=-IH
      RETURN
   10 IH=-IH
      GO TO 1
      END
      SUBROUTINE TETAL(T1,X,Y,IB,IQ,ISQM,K,IC,ISYM,IT20,IT40,ITLIM)
      DIMENSION T1(1),IC(1),KST(30)
      NX=100000.0*X+0.5
      IDH=ISQM/NX
      DO 5 I=1,IB
      KST(I)=0
    5 CONTINUE
      IH=0
      IK=0
      IL=0
      IHH=0
      IKK=0
      ILL=0
      IHA=0
      GO TO(1,2,3,4),ISYM
    1 IF(IHH-IDH) 11,12,12
   11 IH=IH+1
      IHH=IH*IH
      GO TO 17
   12 IF(IKK-IDH) 13,15,15
   13 IK=IK+1
      IKK=IK*IK
   14 IH=IK
      IHH=IKK
      IHA=IKK+ILL
      GO TO 17
   15 IF(ILL-IDH) 16,50,50
   16 IL=IL+1
      ILL=IL*IL
      IK=IL
      IKK=ILL
      GO TO 14
   17 IHB=IHH+IHA
      IF(IHB-IDH) 18,18,12
   18 ISQC=IHB*NX
      IF(ISQM-ISQC) 12,19,19
   19 IF(ISQC-ITLIM) 20,21,21
   20 IT=IT20
      GO TO 22
   21 IT=IT40
   22 DO 23 I=1,IB
      ISQO=IC(I)
      ITQC=ISQO-ISQC
      IF(IABS(ITQC)-IT) 24,24,25
   25 IF(ITQC) 23,23,70
   24 KST(I)=1
   23 CONTINUE
   70 GO TO(1,30,3,30),ISYM
C     TETRAGONAL AND HEXAGONAL SYMMETRY
    2 ISM=0
      GO TO 45
    4 ISM=1
   45 NY=100000.0*Y+0.5
      IDL=ISQM/NY
   30 IF(IHH-IDH) 31,32,32
   31 IH=IH+1
      IHH=IH*IH
   35 IF(ISM) 40,41,40
   41 IHC=IHH+IKK
      GO TO 42
   40 IHC=IHH+IKK+IH*IK
   42 IF(IHC-IDH) 38,38,32
   32 IF(IKK-IDH) 33,36,36
   33 IK=IK+1
      IKK=IK*IK
      IH=IK
      IHH=IKK
      GO TO 35
   36 IF(ILL-IDL) 37,50,50
   37 IL=IL+1
      ILL=IL*IL
      IHA=ILL*NY
      IH=0
      IK=0
      IKK=0
      IHC=0
      IHH=0
   38 ISQC=IHC*NX+IHA
      IF(ISQM-ISQC) 32,19,19
   50 K=0
      DO 51 I=1,IB
      K=K+KST(I)
   51 CONTINUE
    3 RETURN
      END
      SUBROUTINE SEL(IS1,IS2,IS3,NU,ISC)
      IF(ISC) 21,20,21
   20 GO TO(1,2,3,4,5,6,7),NU
   21 GO TO(22,23,24,6),NU
   22 IS1=ISC
      IS2=1
      IS3=2
      GO TO 10
   23 IS3=3
      GO TO 10
   24 IS2=2
      IS3=3
      GO TO 10
    1 IS1=1
      IS2=2
      IS3=3
   10 NU=NU+1
      RETURN
    2 IS3=4
      GO TO 10
    3 IS3=5
      GO TO 10
    4 IS2=3
      IS3=4
      GO TO 10
    5 IS1=2
      GO TO 10
    6 IS1=1
      IS2=2
      IS3=6
      GO TO 10
    7 NU=0
      RETURN
      END
      SUBROUTINE SELE(IS1,IS2,IS3,IS4,NU,ISC,MS,IOUT)
      IF(NU-3) 30,30,31
   31 IF(MS-NU) 13,32,32
   30 IF(ISC-5)2,2,1
    1 GO TO(10,11,12),NU
    2 GO TO(20,21,22),NU
   32 GO TO(13,13,13,35,36,37,38,13),NU
   10 IS1=ISC
      IS2=1
      IS3=2
      IS4=3
    3 NU=NU+1
      RETURN
   11 IS4=4
      GO TO 3
   12 IS3=3
      GO TO 3
   13 NU=0
      RETURN
   20 IS1=1
      IS2=2
      IS3=3
      IS4=4
      GO TO 3
   21 IS4=5
      GO TO 3
   22 IS3=4
      GO TO 3
   35 IS1=1
      IS2=3
      IS3=4
      IS4=5
   33 WRITE(IOUT,34) IS1,IS2,IS3,IS4
   34 FORMAT(21H SELECTED BASE LINES=,4I3)
      GO TO 3
   36 IS2=2
      IS3=3
      IS4=6
      GO TO 33
   37 IS1=2
      IS2=3
      IS3=4
      IS4=5
      GO TO 33
   38 IS1=1
      IS2=2
      IS3=3
      IS4=7
      GO TO 33
      END
      SUBROUTINE MOC(IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3,IH4,IK4,IL4,
     *SQ1,SQ2,SQ3,SQ4,X,Y,Z,U,IO1,IO2,IO3,ISOLV,XL,UL,VTEST,DEN,IZTAL,
     *VOL1,VOL2)
      DIMENSION VOL1(100),VOL2(100)
      SAVE
      IF(IO1) 1,2,1
    1 IO1=0
      IHH1=IH1*IH1
      IKK1=IK1*IK1
      ILL1=IL1*IL1
      IHL1=IH1*IL1
      RHH1=FLOAT(IHH1)
      RKK1=FLOAT(IKK1)
      RLL1=FLOAT(ILL1)
      RHL1=FLOAT(IHL1)
    2 IF(IO2) 3,4,3
    3 IO2=0
      MA2=1
      MB2=1
      MC2=1
      MD2=1
      MUS=1
      IHH2=IH2*IH2
      IKK2=IK2*IK2
      ILL2=IL2*IL2
      IHL2=IH2*IL2
      RHH2=FLOAT(IHH2)
      RKK2=FLOAT(IKK2)
      RLL2=FLOAT(ILL2)
      RHL2=FLOAT(IHL2)
    4 IF(IO3) 5,6,5
    5 IO3=0
      MI=1
      MA=1
      MB=1
      MC=1
      MD=1
      IHH3=IH3*IH3
      IKK3=IK3*IK3
      ILL3=IL3*IL3
      IHL3=IH3*IL3
      RHH3=FLOAT(IHH3)
      RKK3=FLOAT(IKK3)
      RLL3=FLOAT(ILL3)
      RHL3=FLOAT(IHL3)
    6 IHH4=IH4*IH4
      IKK4=IK4*IK4
      ILL4=IL4*IL4
      IHL4=IH4*IL4
      IF(MI) 101,102,101
  102 IDET=IHL4*MMD - ILL4*MMC + IKK4*MMB -IHH4*MMA
      GO TO 105
  101 MI=0
      IF(MUS) 103,104,103
  103 IAA=ILL1*IHL2 - ILL2*IHL1
      IAB=IKK1*IHL2 - IKK2*IHL1
      IAC=IKK1*ILL2 - IKK2*ILL1
      IBB=IHH1*IHL2 - IHH2*IHL1
      IBC=IHH1*ILL2 - IHH2*ILL1
      ICC=IHH1*IKK2 - IHH2*IKK1
      MUS=0
  104 MMA=IKK3*IAA - ILL3*IAB + IHL3*IAC
      MMB=IHH3*IAA - ILL3*IBB + IHL3*IBC
      MMC=IHH3*IAB - IKK3*IBB + IHL3*ICC
      MMD=IHH3*IAC - IKK3*IBC + ILL3*ICC
      GO TO 102
  105 IF(IDET) 8,7,8
    7 ISOLV=0
      RETURN
    8 RDET=FLOAT(IDET)
      RHH4=FLOAT(IHH4)
      RKK4=FLOAT(IKK4)
      RLL4=FLOAT(ILL4)
      RHL4=FLOAT(IHL4)
      IF(MD) 201,202,201
  202 UTAK=SQ4*RUD - RLL4*RUC + RKK4*RUB - RHH4*RUA
      GO TO 205
  201 MD=0
      IF(MD2) 203,204,203
  203 RUAA=RLL1*SQ2 - RLL2*SQ1
      RUAB=RKK1*SQ2 - RKK2*SQ1
      RUAC=RKK1*RLL2 - RKK2*RLL1
      RUBB=RHH1*SQ2 - RHH2*SQ1
      RUBC=RHH1*RLL2 - RHH2*RLL1
      RUCC=RHH1*RKK2 - RHH2*RKK1
      MD2=0
  204 RUA=RKK3*RUAA - RLL3*RUAB + SQ3*RUAC
      RUB=RHH3*RUAA - RLL3*RUBB + SQ3*RUBC
      RUC=RHH3*RUAB - RKK3*RUBB + SQ3*RUCC
      RUD=RHH3*RUAC - RKK3*RUBC + RLL3*RUCC
      GO TO 202
  205 U=UTAK/RDET
      IF(U) 7,11,11
   11 IF(MA) 301,302,301
  302 XTAK=RHL4*RXD - RLL4*RXC + RKK4*RXB - SQ4*RXA
      GO TO 305
  301 MA=0
      IF(MA2) 303,304,303
  303 RXAA=RLL1*RHL2 - RLL2*RHL1
      RXAB=RKK1*RHL2 - RKK2*RHL1
      RXAC=RKK1*RLL2 - RKK2*RLL1
      RXBB=SQ1*RHL2 - SQ2*RHL1
      RXBC=SQ1*RLL2 - SQ2*RLL1
      RXCC=SQ1*RKK2 - SQ2*RKK1
      MA2=0
  304 RXA=RKK3*RXAA - RLL3*RXAB + RHL3*RXAC
      RXB=SQ3*RXAA - RLL3*RXBB + RHL3*RXBC
      RXC=SQ3*RXAB - RKK3*RXBB + RHL3*RXCC
      RXD=SQ3*RXAC - RKK3*RXBC + RLL3*RXCC
      GO TO 302
  305 X=XTAK/RDET
      IF(X-XL) 7,9,9
    9 IF(MB) 401,402,401
  402 YTAK=RHL4*RYD - RLL4*RYC + SQ4*RYB - RHH4*RYA
      GO TO 405
  401 MB=0
      IF(MB2) 403,404,403
  403 RYAA=RLL1*RHL2 - RLL2*RHL1
      RYAB=SQ1*RHL2 - SQ2*RHL1
      RYAC=SQ1*RLL2 - SQ2*RLL1
      RYBB=RHH1*RHL2 - RHH2*RHL1
      RYBC=RHH1*RLL2 - RHH2*RLL1
      RYCC=RHH1*SQ2 - RHH2*SQ1
      MB2=0
  404 RYA=SQ3*RYAA - RLL3*RYAB + RHL3*RYAC
      RYB=RHH3*RYAA - RLL3*RYBB + RHL3*RYBC
      RYC=RHH3*RYAB - SQ3*RYBB + RHL3*RYCC
      RYD=RHH3*RYAC - SQ3*RYBC + RLL3*RYCC
      GO TO 402
  405 Y=YTAK/RDET
      IF(Y-XL) 7,10,10
   10 IF(MC) 501,502,501
  502 ZTAK=RHL4*RZD - SQ4*RZC + RKK4*RZB - RHH4*RZA
      GO TO 505
  501 MC=0
      IF(MC2) 503,504,503
  503 RZAA=SQ1*RHL2 - SQ2*RHL1
      RZAB=RKK1*RHL2 - RKK2*RHL1
      RZAC=RKK1*SQ2 - RKK2*SQ1
      RZBB=RHH1*RHL2 - RHH2*RHL1
      RZBC=RHH1*SQ2 - RHH2*SQ1
      RZCC=RHH1*RKK2 - RHH2*RKK1
      MC2=0
  504 RZA=RKK3*RZAA - SQ3*RZAB + RHL3*RZAC
      RZB=RHH3*RZAA - SQ3*RZBB + RHL3*RZBC
      RZC=RHH3*RZAB - RKK3*RZBB + RHL3*RZCC
      RZD=RHH3*RZAC - RKK3*RZBC + SQ3*RZCC
      GO TO 502
  505 Z=ZTAK/RDET
      IF(Z-XL) 7,13,13
   13 HALP=X*Z
      UTE=(U**2)/HALP
      IF(UTE-UL) 12,12,7
   12 HELP=1.0-0.25*UTE
      VST2=Y*HALP*HELP
      IF(VTEST-VST2) 700,700,7
  700 IF(X*HELP-XL) 7,15,15
   15 IF(Z*HELP-XL) 7,14,14
   14 IF(DEN-0.00001) 16,16,701
  701 DO 703 I=1,IZTAL
      IF(VST2-VOL1(I)) 703,702,702
  702 IF(VST2-VOL2(I)) 16,16,7
  703 CONTINUE
      GO TO 7
   16 ISOLV=1
      RETURN
      END
      SUBROUTINE MAEG(T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      DIMENSION T1(1),IC(1),KST(30)
      NX=100000.0*X
      NY=100000.0*Y
      NZ=100000.0*Z
      NU=100000.0*U
      IDH=ISQM/NX
      IDK=ISQM/NY
      IDL=ISQM/NZ
      NHM=0
      NKM=0
      NLM=0
    1 NHM=NHM+1
      IF(NHM*NHM-IDH) 1,1,2
    2 NHM=NHM-1
    3 NKM=NKM+1
      IF(NKM*NKM-IDK) 3,3,4
    4 NKM=NKM-1
    5 NLM=NLM+1
      IF(NLM*NLM-IDL) 5,5,6
    6 NLM=NLM-1
      DO 40 I=1,IB
      KST(I)=0
   40 CONTINUE
      IH=0
      IK=0
      IL=0
      IHA=0
      IHB=0
      IHC=0
      IHD=0
      IHE=0
      IHF=0
      IHG=0
    7 IF(IH)8,10,11
    8 IH=-IH
    9 IF(IH-NHM) 10,13,13
   10 IH=IH+1
      GO TO 18
   11 IF(IL) 12,9,12
   12 IH=-IH
      IHD=-IHD
      GO TO 19
   13 IH=0
      IHA=0
      IF(IL-NLM) 14,15,15
   14 IL=IL+1
      GO TO 17
   15 IL=0
      IHC=0
      IHD=0
      IF(IK-NKM) 16,27,27
   16 IK=IK+1
      IHB=IK*IK*NY
      IHE=IHB+IHC
      ISQC=IHB
      GO TO 20
   17 IHC=IL*IL*NZ
      IHE=IHB+IHC
      IF(ISQM-IHE) 10,18,18
   18 IHA=IH*IH*NX
      IHF=IHA+IHE
      IHD=IH*IL*NU
   19 ISQC=IHF+IHD
   20 IF(ISQM-ISQC) 7,21,21
   21 IF(ISQC-ITLIM) 22,23,23
   22 IT=IT20
      GO TO 24
   23 IT=IT40
   24 DO 26 I=1,IB
      ISQO=IC(I)
      IF(IABS(ISQO-ISQC)-IT) 25,25,26
   25 KST(I)=1
   26 CONTINUE
      GO TO 7
   27 K=0
      DO 28 I=1,IB
      K=K+KST(I)
   28 CONTINUE
      RETURN
      END
      SUBROUTINE LAGRM(K,U)
      COMMON T1(100),TAGUT(100),N1,D1,SSQTL,D2,WAVE,IHM,IKM,ILM,ISYM,INR
     1,ISUB,A11,A22,A33,A13,A12,A23,LIN(20),REV(20)
     2,IGOT,IS1,IS2,IS3,IS4,IA,IB,IQ,IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3
     3,SQ1,SQ2,SQ3,X,Y,Z,IC(20),ISQM,KOTA(30),RO(7,30),ID,IIN,IOUT
      COMMON NH1,NK1,NL1,NS1,NH2,NK2,NL2,NS2,NH3,NK3,NL3,NS3,IUTSK,IMER,
     1INIX,IOMIT,ISC,KLART,VOL,VTEST,XL,YL,IT20,IT40,CEM
     2,KHA,KKA,KLA,KSA,KHB,KKB,KLB,KSB,MONO,MH1,MK1,ML1,MH2,MK2,ML2,MH3,
     3MK3,ML3,MS1,MS2,MS3,MH4,MK4,ML4,MS4,UL,LGO,MG,MS,MRA,MRB,NU,
     4DEN,EDEN,RMOLV
      IF(X-Z) 1,2,2
    1 XA=Z
      ZA=X
      GO TO 3
    2 XA=X
      ZA=Z
    3 UTE=(U*U)/(X*Z)
      PV=XA*Y*ZA*(1.0-UTE*0.25)
      IF(ID) 7,4,7
    4 ID=ID+1
    5 I=ID
    6 RO(1,I)=XA
      RO(2,I)=Y
      RO(3,I)=ZA
      RO(4,I)=U
      KOTA(I)=K
      RO(5,I)=PV
      RETURN
    7 DO 11 I=1,ID
      IF(ABS(RO(1,I)-XA)-0.00005) 8,11,11
    8 IF(ABS(RO(2,I)-Y)-0.00005) 9,11,11
    9 IF(ABS(RO(3,I)-ZA)-0.00005) 10,11,11
   10 IF(ABS(RO(4,I)-U)-0.00005) 20,11,11
   20 RETURN
   11 CONTINUE
      DO 14 I=1,ID
      IF(KOTA(I)-K) 13,12,14
   12 IF(RO(5,I)-PV) 13,14,14
   13 IF(I-5) 15,5,99
   14 CONTINUE
      IF(ID-5) 4,20,20
   15 N=4
      IF(ID-5) 16,17,17
   16 ID=ID+1
   17 M=N+1
      RO(5,M)=RO(5,N)
      KOTA(M)=KOTA(N)
      DO 18 J=1,4
      RO(J,M)=RO(J,N)
   18 CONTINUE
      IF(N-I) 99,6,19
   19 N=N-1
      GO TO 17
   99 STOP
      END
      SUBROUTINE IMP(IH1,IK1,IL1,IH2,IK2,IL2,IOK)
      IF(IH1) 1,2,2
    1 IOK=1
      RETURN
    2 IF(IL1-IL2) 1,3,3
    3 IF(IK1-IK2) 1,4,4
    4 IF(IH1-IABS(IH2)) 1,5,5
    5 IOK=0
      RETURN
      END
      SUBROUTINE TREOC(IGOE)
      DIMENSION NVD(210),NDS(210)
      COMMON T1(100),TAGUT(100),N1,D1,SSQTL,D2,WAVE,IHM,IKM,ILM,ISYM,INR
     1,ISUB,A11,A22,A33,A13,A12,A23,LIN(20),REV(20)
     2,IGOT,IS1,IS2,IS3,IS4,IA,IB,IQ,IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3
     3,SQ1,SQ2,SQ3,X,Y,Z,IC(20),ISQM,KOTA(30),RO(7,30),ID,IIN,IOUT
      COMMON NH1,NK1,NL1,NS1,NH2,NK2,NL2,NS2,NH3,NK3,NL3,NS3,IUTSK,IMER,
     1INIX,IOMIT,ISC,KLART,VOL,VTEST,XL,YL,IT20,IT40,CEM
     2,KHA,KKA,KLA,KSA,KHB,KKB,KLB,KSB,MONO,MH1,MK1,ML1,MH2,MK2,ML2,MH3,
     3MK3,ML3,MS1,MS2,MS3,MH4,MK4,ML4,MS4,UL,LGO,MG,MS,MRA,MRB,NU,
     4DEN,EDEN,RMOLV
C     COMMON CARDS
      A11=0.0
      A22=0.0
      A33=0.0
      IF(IGOT) 2,1,2
    1 LIN(10)=ISYM
C      FROM TREOB
      IGOT=1
    2 IF(ID) 3,4,3
    3 ISYM=LIN(10)
      INR=3
      I=MRA+1-ID
      K=KOTA(I)
      ID=ID-1
      GO TO(11,12,13,12,150,160),ISYM
  160 A11=RO(1,I)
      A22=RO(2,I)
      A33=RO(3,I)
      A12=RO(4,I)
      A13=RO(5,I)
      A23=RO(6,I)
      WRITE(IOUT,161) K,A11,A22,A33,A12,A13,A23
  161 FORMAT(3H K=,I3,8H A11-33=,3F9.6,8H A12-23=,3F9.6)
      GO TO 30
  150 A11=RO(3,I)
      A22=RO(2,I)
      A33=RO(1,I)
      A13=RO(4,I)
      WRITE(IOUT,151)K,A11,A22,A33,A13
  151 FORMAT(3H K=,I3,6H XYZU=,4F10.6)
      GO TO 30
   11 A11=RO(1,I)
      WRITE(IOUT,20) K,A11
   20 FORMAT(3H K=,I3,3H X=,F10.5)
      GO TO 30
   12 A11=RO(1,I)
      A22=RO(2,I)
      WRITE(IOUT,21) K,A11,A22
   21 FORMAT(3H K=,I3,4H XY=,2F10.5)
      GO TO 30
   13 A11=RO(3,I)
      A22=RO(2,I)
      A33=RO(1,I)
      WRITE(IOUT,22) K,A11,A22,A33
   22 FORMAT(3H K=,I3,5H XYZ=,3F10.5)
   30 IGOE=1
      RETURN
    4 ISYM=LIN(10)
      GO TO(31,32,152,34,33,162),ISYM
  162 IF(LIN(16) .EQ. 0) GO TO 170
      IGOE=2
      RETURN
  170 WRITE(IOUT,163)
  163 FORMAT(22H END OF TRICLINIC TEST)
      IGOE=3
      RETURN
   31 ISYM=2
      GO TO 40
   32 ISYM=4
      GO TO 40
   34 ISYM=3
      GO TO 40
  152 ISYM=5
   40 IGOE=2
      RETURN
   33 MSW=LIN(15)
      IF(MSW .LT. 3) GO TO 40
      IF(MG) 40,153,40
  153 IF(IB-20) 51,50,50
   50 NDTAL=20
      GO TO 52
   51 NDTAL=IB
   52 WRITE(IOUT,101)
  101 FORMAT(32H PROGRAM FOR DIFFERENCE ANALYSIS,/,57H PROGRAM  I1  REF:
     1 WERNER,P-E., Z.KRIST. 120(1964)375-378)
      WRITE(IOUT,53)NDTAL,D2,D2
   53 FORMAT(35H NUMBER OF SINE SQUARE(THETA) USED=,I5,/,33H  DELTA=X/16
     1 IF THIS IS LESS THAN,F10.5,17H OTHERWISE DELTA=,F10.5)
      IXT=800
      DO 54 I=1,210
      NVD(I)=0
      NDS(I)=0
   54 CONTINUE
      L=0
      J=0
      I=NDTAL+1
   55 I=I-1
      IF(I-2) 59,56,56
   56 J=J+1
      NVD(J)=IC(I)
      K=I
   57 K=K-1
      IF(K) 55,55,58
   58 J=J+1
      NVD(J)=IC(I)-IC(K)
      GO TO 57
   59 J=J+1
      NVD(J)=IC(I)
      NUTAL=J
      DO 68 I=1,NUTAL
      NDSI=0
      NVDI=NVD(I)
      NDELT=NVDI/16
      IF(NDELT-IT40) 61,61,60
   60 NDELT=IT40
   61 DO 67 J=1,NUTAL
      IF(I-J) 62,67,62
   62 NTEST=NVD(J)
      N=NTEST/NVDI
      IF(N) 67,67,200
  200 NMULT=N*NVDI
   66 IF(NMULT+NDELT-NTEST) 65,63,63
   63 IF(NMULT-NDELT-NTEST) 64,64,67
   64 NDSI=NDSI+1
      GO TO 67
   65 NMULT=NMULT+NVDI
      GO TO 66
   67 CONTINUE
      NDS(I)=NDSI
   68 CONTINUE
      MAXDF=IC(NDTAL)
      RAXDF=FLOAT(MAXDF)
      NUMER=0
      WRITE(IOUT,69)
   69 FORMAT(22H   NR MULT. DIFF.TABLE)
      DO 75 I=1,NUTAL
      NSMAX=-1
      DO 71 J=1,NUTAL
      N=NDS(J)
      IF(N-NSMAX) 71,71,70
   70 NSMAX=N
      JAKT=J
   71 CONTINUE
      NVDA=NVD(JAKT)
      SQDIF=FLOAT(NVDA)/100000.0
      NDS(JAKT)=-1
      NFORM=0
   72 NFORM=NFORM+1
      IF(NFORM-10) 73,74,74
   73 RVAL=RAXDF*(0.66667**NFORM)
      RNVDA=FLOAT(NVDA)
      IF(RVAL-RNVDA) 74,74,72
   74 NUMER=NUMER+1
      GO TO(81,82,83,84,85,86,87,88,89,90),NFORM
   81 WRITE(IOUT,91) NUMER,NSMAX,SQDIF
   91 FORMAT(2I5,81X,F9.6)
      GO TO 75
   82 WRITE(IOUT,92) NUMER,NSMAX,SQDIF
   92 FORMAT(2I5,72X,F9.6)
      GO TO 75
   83 WRITE(IOUT,93) NUMER,NSMAX,SQDIF
   93 FORMAT(2I5,63X,F9.6)
      GO TO 75
   84 WRITE(IOUT,94) NUMER,NSMAX,SQDIF
   94 FORMAT(2I5,54X,F9.6)
      GO TO 75
   85 WRITE(IOUT,95) NUMER,NSMAX,SQDIF
   95 FORMAT(2I5,45X,F9.6)
      GO TO 75
   86 WRITE(IOUT,96) NUMER,NSMAX,SQDIF
   96 FORMAT(2I5,36X,F9.6)
      GO TO 75
   87 WRITE(IOUT,97) NUMER,NSMAX,SQDIF
   97 FORMAT(2I5,27X,F9.6)
      GO TO 75
   88 WRITE(IOUT,98) NUMER,NSMAX,SQDIF
   98 FORMAT(2I5,18X,F9.6)
      GO TO 75
   89 WRITE(IOUT,99) NUMER,NSMAX,SQDIF
   99 FORMAT(2I5,9X,F9.6)
      GO TO 75
   90 WRITE(IOUT,100) NUMER,NSMAX,SQDIF
  100 FORMAT(2I5,F9.6)
   75 CONTINUE
      WRITE(IOUT,76)
   76 FORMAT(14H TEST FINISHED)
      IGOE=3
      RETURN
      END
      SUBROUTINE PW71(T2,IHT,IKT,ILT,IQE,ISYM,N2,A11,A22,A33,A13,A12,A23
     1,V2,NEKV)
      DIMENSION T2(1),IHT(1),IKT(1),ILT(1),IQE(1)
      V2=0
      NEKV=0
      DO 7 I=1,N2
      IF(IQE(I)) 8,7,8
    8 GO TO(1,2,3,4,5,6),ISYM
    1 SQ=(IHT(I)**2+IKT(I)**2+ILT(I)**2)*A11
      GO TO 9
    2 SQ=(IHT(I)**2+IKT(I)**2)*A11+ILT(I)*ILT(I)*A22
      GO TO 9
    3 SQ=(IHT(I)**2)*A11+(IKT(I)**2)*A22+(ILT(I)**2)*A33
      GO TO 9
    4 SQ=(IHT(I)**2+IHT(I)*IKT(I)+IKT(I)**2)*A11+(ILT(I)**2)*A22
      GO TO 9
    5 SQ=(IHT(I)**2)*A11+(IKT(I)**2)*A22+(ILT(I)**2)*A33+(IHT(I)*ILT(I)*
     1A13)
      GO TO 9
    6 SQ=(IHT(I)**2)*A11+(IKT(I)**2)*A22+(ILT(I)**2)*A33+ IHT(I)*ILT(I)*
     1A13+IHT(I)*IKT(I)*A12+IKT(I)*ILT(I)*A23
    9 V2=V2+(T2(I)-SQ)**2
      NEKV=NEKV+1
    7 CONTINUE
      RETURN
      END
      SUBROUTINE PW50(A11,A22,A33,A13,A12,A23,ISYM,IHM,IKM,ILM,IH,IK,IL,
     1SSQT,IA)
C     SE P-E 850
      IA=1
      GO TO(1,2,3,4,5,6),ISYM
    1 IF(IH-IHM) 7,8,8
    7 IH=IH+1
      GO TO 10
    8 IF(IK-IKM) 9,11,11
    9 IK=IK+1
      IH=IK
      GO TO 10
   11 IF(IL-ILM) 12,13,13
   12 IL=IL+1
      IK=IL
      IH=IK
   10 SSQT=FLOAT(IH*IH+IK*IK+IL*IL)*A11
      RETURN
   13 IA=0
      RETURN
    2 IF(IH-IHM) 14,15,15
   14 IH=IH+1
      GO TO 20
   15 IF(IK-IKM) 16,17,17
   16 IK=IK+1
      IH=IK
      GO TO 20
   17 IF(IL-ILM) 18,13,13
   18 IL=IL+1
      IH=0
      IK=0
   20 SSQT=FLOAT(IH*IH+IK*IK)*A11+FLOAT(IL*IL)*A22
      RETURN
    3 IF(IH-IHM) 21,22,22
   21 IH=IH+1
      GO TO 30
   22 IH=0
      IF(IK-IKM) 23,24,24
   23 IK=IK+1
      GO TO 30
   24 IF(IL-ILM) 25,13,13
   25 IK=0
      IL=IL+1
   30 SSQT=FLOAT(IH*IH)*A11+FLOAT(IK*IK)*A22+FLOAT(IL*IL)*A33
      RETURN
C     CHANGE S50/46,47,48,49,50,54
C
C
    4 IF(IH-IHM) 33,34,34
   33 IH=IH+1
      GO TO 40
   34 IF(IK-IKM) 35,36,36
   35 IK=IK+1
      IH=IK
      GO TO 40
   36 IF(IL-ILM) 37,13,13
   37 IH=0
      IK=0
      IL=IL+1
   40 SSQT=FLOAT(IH*IH+IH*IK+IK*IK)*A11+FLOAT(IL*IL)*A22
      RETURN
    5 IF(IH*IL)42,42,41
   41 IH=-IH
      GO TO 50
   42 IF(IABS(IH)-IHM) 43,44,44
   43 IH=IABS(IH)+1
      GO TO 50
   44 IH=0
      IF(IK-IKM)45,46,46
   45 IK=IK+1
      GO TO 50
   46 IF(IL-ILM) 47,13,13
   47 IK=0
      IL=IL+1
   50 SSQT=FLOAT(IH*IH)*A11+FLOAT(IK*IK)*A22+FLOAT(IL*IL)*A33+FLOAT(IH*I
     1L)*A13
      RETURN
    6 IF(IABS(IK)+IL) 51,52,51
   51 IF(IH) 52,52,53
   53 IH=-IH
      GO TO 60
   52 IF(IABS(IH)-IHM) 54,55,55
   54 IH=IABS(IH)+1
      GO TO 60
   55 IH=0
      IF(IL) 57,57,56
   56 IF(IK) 57,57,58
   58 IK=-IK
      GO TO 60
   57 IF(IABS(IK)-IKM) 59,61,61
   59 IK=IABS(IK)+1
      GO TO 60
   61 IK=0
      IF(IL-ILM) 62,13,13
   62 IL=IL+1
   60 SSQT=FLOAT(IH*IH)*A11+FLOAT(IK*IK)*A22+FLOAT(IL*IL)*A33+FLOAT(IH*I
     1L)*A13+FLOAT(IH*IK)*A12+FLOAT(IK*IL)*A23
      RETURN
      END
      SUBROUTINE PW51(SSQT,T1,N,SSQTL,D1,D2,SSQTO,IA)
      DIMENSION T1(1)
C     SE P-E 852
      IF(SSQT-SSQTL) 1,2,2
    1 D=D1
      GO TO 3
    2 D=D2
    3 DO 4 I=1,N
      IF(ABS(SSQT-T1(I))-D) 5,5,4
    4 CONTINUE
      IA=0
      RETURN
    5 IA=1
      SSQTO=T1(I)
      RETURN
      END
      SUBROUTINE PW53(SSQTO,T2,N2,IHT,IKT,ILT,IQE,IH,IK,IL,ISYM)
C     SE P-E 853 OBS. N2 START=0
      DIMENSION T2(1),IHT(1),IKT(1),ILT(1),IQE(1)
      I=1
    5 IF(I-N2) 2,2,1
    1 T2(I)=SSQTO
      IHT(I)=IH
      IKT(I)=IK
      ILT(I)=IL
      IQE(I)=1
      GO TO 26
    2 IF(SSQTO-T2(I)) 3,4,3
    3 I=I+1
      GO TO 5
    4 GO TO (11,12,13,14,15,16),ISYM
   11 IF((IH**2+IK**2+IL**2)-(IHT(I)**2+IKT(I)**2+ILT(I)**2))7,6,7
    7 IQE(I)=0
      RETURN
   12 IF((IH**2+IK**2)-(IHT(I)**2+IKT(I)**2)) 7,21,7
   21 IF(IL**2-ILT(I)**2) 7,6,7
   13 IF(IH**2-IHT(I)**2) 7,22,7
   22 IF(IK**2-IKT(I)**2) 7,21,7
   14 IF((IH**2+IH*IK+IK**2)-(IHT(I)**2+IHT(I)*IKT(I)+IKT(I)**2))7,21,7
   15 IF(IH*IL-IHT(I)*ILT(I)) 7,13,7
   16 IF(IH*IK-IHT(I)*IKT(I)) 7,23,7
   23 IF(IK*IL-IKT(I)*ILT(I)) 7,15,7
   26 SSQTO=-SSQTO
    6 RETURN
      END
      SUBROUTINE PW84(A,N,B,ISOLB)
C     M=0 GER AINV,M=1 DESSUTOM X I B-VEKTORN
      DIMENSION IPIV(6),A(6,6),B(6,1),INDEX(6,2),PIVOT(6)
      EQUIVALENCE (IROW,JROW),(ICOL,JCOL),(AMAX,T,SWAP)
      M=1
      ISOLB=1
      DET=1.0
      DO 20 J=1,N
   20 IPIV(J)=0
      DO 550 I=1,N
      AMAX=0.0
      DO 105 J=1,N
      IF(IPIV(J)-1) 60,105,60
   60 DO 100 K=1,N
      IF(IPIV(K)-1) 80,100,740
   80 IF(ABS(AMAX)-ABS(A(J,K))) 85,100,100
   85 IROW=J
      ICOL=K
      AMAX=A(J,K)
  100 CONTINUE
  105 CONTINUE
      IPIV(ICOL)=IPIV(ICOL)+1
      IF(IROW-ICOL) 140,260,140
  140 DET=-DET
      DO 200 L=1,N
      SWAP=A(IROW,L)
      A(IROW,L)=A(ICOL,L)
  200 A(ICOL,L)=SWAP
      IF(M)260,260,210
  210 DO 250 L=1,M
      SWAP=B(IROW,L)
      B(IROW,L)=B(ICOL,L)
  250 B(ICOL,L)=SWAP
  260 INDEX(I,1)=IROW
      INDEX(I,2)=ICOL
      PIVOT(I)=A(ICOL,ICOL)
      IF(PIVOT(I)) 270,271,270
  270 DET=DET*PIVOT(I)
      A(ICOL,ICOL)=1.0
      DO 350 L=1,N
  350 A(ICOL,L)=A(ICOL,L)/PIVOT(I)
      IF(M) 380,380,360
  360 DO 370 L=1,M
  370 B(ICOL,L)=B(ICOL,L)/PIVOT(I)
  380 DO 550 L1=1,N
      IF(L1-ICOL) 400,550,400
  400 T=A(L1,ICOL)
      A(L1,ICOL)=0.0
      DO 450 L=1,N
  450 A(L1,L)=A(L1,L)-A(ICOL,L)*T
      IF(M) 550,550,460
  460 DO 500 L=1,M
  500 B(L1,L)=B(L1,L)-B(ICOL,L)*T
  550 CONTINUE
      DO 710 I=1,N
      L=N+1-I
      IF(INDEX(L,1)-INDEX(L,2)) 630,710,630
  630 JROW=INDEX(L,1)
      JCOL=INDEX(L,2)
      DO 705 K=1,N
      SWAP=A(K,JROW)
      A(K,JROW)=A(K,JCOL)
      A(K,JCOL)=SWAP
  705 CONTINUE
  710 CONTINUE
  740 RETURN
  271 ISOLB=0
      RETURN
      END
      SUBROUTINE PW55(T2,IQE,IHT,IKT,ILT,NTOT,ISYM,A,B,IEKV)
C     BILDA NORMALEKV. MATRISERNA A OCH B. ORIML.PROBLEM GER IEKV=0
      DIMENSION T2(1), IQE(1),IHT(1),IKT(1),ILT(1),A(6,6),B(6,1)
      DO 7 I=1,6
      DO 7 J=1,6
    7 A(I,J)=0.0
      DO 8 I=1,6
    8 B(I,1)=0.0
      IEKV=0
      GO TO(1,2,3,4,3,3),ISYM
    1 DO 9 I=1,NTOT
      IF(IQE(I))10,9,10
   10 U1=FLOAT(IHT(I)**2+IKT(I)**2+ILT(I)**2)
      A(1,1)=U1**2+A(1,1)
      B(1,1)=U1*T2(I)+B(1,1)
      IEKV=IEKV+1
    9 CONTINUE
      IF(IEKV-1)11,11,12
   11 IEKV=0
      GO TO 12
    2 DO 13 I=1,NTOT
      IF(IQE(I)) 14,13,14
   14 U1=FLOAT(IHT(I)**2+IKT(I)**2)
      A(1,1)=U1**2+A(1,1)
      A(1,2)=FLOAT(ILT(I)**2)*U1+A(1,2)
      A(2,2)=FLOAT(ILT(I))**4+A(2,2)
      B(1,1)=U1*T2(I)+B(1,1)
      B(2,1)=FLOAT(ILT(I)**2)*T2(I)+B(2,1)
      IEKV=IEKV+1
   13 CONTINUE
   16 IF(IEKV-2) 11,11,15
   15 A(2,1)=A(1,2)
      GO TO 12
    4 DO 17 I=1,NTOT
      IF(IQE(I))18,17,18
   18 U1=FLOAT(IHT(I)**2+IHT(I)*IKT(I)+IKT(I)**2)
      A(1,1)=U1**2+A(1,1)
      A(1,2)=U1*FLOAT(ILT(I)**2)+A(1,2)
      A(2,2)=FLOAT(ILT(I))**4+A(2,2)
      B(1,1)=U1*T2(I)+B(1,1)
      B(2,1)=FLOAT(ILT(I)**2)*T2(I)+B(2,1)
      IEKV=IEKV+1
   17 CONTINUE
      GO TO 16
    3 DO 19 I=1,NTOT
      IF(IQE(I)) 20,19,20
   20 U1=FLOAT(IHT(I)**2)
      U2=FLOAT(IKT(I)**2)
      U3=FLOAT(ILT(I)**2)
      A(1,1)=U1**2+A(1,1)
      A(1,2)=U1*U2+A(1,2)
      A(1,3)=U1*U3+A(1,3)
      B(1,1)=U1*T2(I)+B(1,1)
      B(2,1)=U2*T2(I)+B(2,1)
      B(3,1)=U3*T2(I)+B(3,1)
      A(2,2)=U2**2+A(2,2)
      A(2,3)=U2*U3+A(2,3)
      A(3,3)=U3**2+A(3,3)
      IEKV=IEKV+1
      GO TO(11,11,19,11,21,21),ISYM
   21 U4=FLOAT(IHT(I)*ILT(I))
      A(1,4)=U1*U4+A(1,4)
      A(2,4)=U2*U4+A(2,4)
      A(3,4)=U3*U4+A(3,4)
      A(4,4)=U4**2+A(4,4)
      B(4,1)=U4*T2(I)+B(4,1)
      GO TO(11,11,11,11,19,22),ISYM
   22 U5=FLOAT(IHT(I)*IKT(I))
      U6=FLOAT(IKT(I)*ILT(I))
      A(1,5)=U1*U5+A(1,5)
      A(1,6)=U1*U6+A(1,6)
      A(2,5)=U2*U5+A(2,5)
      A(2,6)=U2*U6+A(2,6)
      A(3,5)=U3*U5+A(3,5)
      A(3,6)=U3*U6+A(3,6)
      A(4,5)=U4*U5+A(4,5)
      A(4,6)=U4*U6+A(4,6)
      A(5,5)=U5**2+A(5,5)
      A(5,6)=U5*U6+A(5,6)
      A(6,6)=U6**2+A(6,6)
      B(5,1)=U5*T2(I)+B(5,1)
      B(6,1)=U6*T2(I)+B(6,1)
   19 CONTINUE
      A(2,1)=A(1,2)
      A(3,1)=A(1,3)
      A(3,2)=A(2,3)
      GO TO(11,11,23,11,24,24),ISYM
   23 IF(IEKV-3)11,11,12
   24 A(4,1)=A(1,4)
      A(4,2)=A(2,4)
      A(4,3)=A(3,4)
      GO TO(11,11,11,11,25,26),ISYM
   25 IF(IEKV-4)11,11,12
   26 A(5,1)=A(1,5)
      A(5,2)=A(2,5)
      A(5,3)=A(3,5)
      A(5,4)=A(4,5)
      A(6,1)=A(1,6)
      A(6,2)=A(2,6)
      A(6,3)=A(3,6)
      A(6,4)=A(4,6)
      A(6,5)=A(5,6)
      IF(IEKV-6)11,11,12
   12 RETURN
      END
      SUBROUTINE P56W(T1,T2,IHT,IKT,ILT,N1,N2,IOUT,D1,SSQTL,D2,WAVE,TAGU
     1T,N20,N30,NFI,NT20,NT30,SQ20,SQ30,SQFI,M20,NOIND)
C
      DIMENSION T1(1),T2(1),IHT(1),IKT(1),ILT(1),TAGUT(1)
      NOIND=0
      NWIND=0
      M20=0
      IF(SQ20.EQ.0.) GO TO 110
      SAK=SQ20
      GO TO 111
  110 SAK=SQFI
  111 CONTINUE
      DO 1 J1=1,N2
      ISH=IHT(J1)
      ISK=IKT(J1)
      ISL=ILT(J1)
      TS=T2(J1)
      IMH=ISH
      IMK=ISK
      IML=ISL
      TM=TS
      DO 2 J2=J1,N2
      IF(TM-T2(J2)) 2,3,3
    3 TM=T2(J2)
      IMH=IHT(J2)
      IMK=IKT(J2)
      IML=ILT(J2)
      J3=J2
    2 CONTINUE
      IHT(J1)=IMH
      IKT(J1)=IMK
      ILT(J1)=IML
      T2(J1)=TM
      IHT(J3)=ISH
      IKT(J3)=ISK
      ILT(J3)=ISL
    1 T2(J3)=TS
      DSUM=0.0
      DSUMA=0.0
      DSUMB=0.0
      DSUMC=0.0
      TSUMA=0.0
      TSUMB=0.0
      TSUMC=0.0
      NTALA=0
      NTALB=0
      NTALC=0
      GRAD=114.59156
      WRITE(IOUT,4)
    4 FORMAT(1H ,'   H   K   L SST-OBS  SST-CALC   DELTA   2TH-OBS 2TH-C
     1ALC D-OBS   FREE PARAM.')
      J=1
      DO 17 I=1,N2
    5 IF(J-N1) 6,6,15
    6 IF(I-N2) 7,8,8
    7 IF(ABS(T1(J)-T2(I))-ABS(T1(J)-T2(I+1))) 8,8,15
    8 IF(T1(J)-SSQTL) 9,9,10
    9 IF(ABS(T1(J)-T2(I))-D1) 13,13,11
   10 IF(ABS(T1(J)-T2(I))-D2) 13,13,11
   11 IF(T1(J)-T2(I)) 12,12,15
   12 HA=SQRT(T1(J))
      DOBS=WAVE/(2.*HA)
      HD=SQRT(1.-T1(J))
      THO=GRAD*ATAN(HA/HD)
      NOIND=NOIND+1
      WRITE(IOUT,20) T1(J),THO,DOBS,TAGUT(J)
   20 FORMAT(1H ,12X,F9.6,19X,F8.3,8X,F8.4,2X,6X,A4)
      IF(T1(J)-SAK) 112,112,113
  112 NWIND=NWIND+1
  113 J=J+1
      GO TO 5
   13 HA=SQRT(T1(J))
      HB=SQRT(T2(I))
      DOBS=WAVE/(2.*HA)
      DCALC=WAVE/(2.*HB)
      DELTA=T1(J)-T2(I)
      HD=SQRT(1.-T2(I))
      HE=SQRT(1.-T1(J))
      THO=ATAN(HA/HE)*GRAD
      THC=ATAN(HB/HD)*GRAD
      DELTB=THO-THC
      WRITE(IOUT,14) IHT(I),IKT(I),ILT(I),T1(J),T2(I),DELTA,THO,THC,DOBS
     1,TAGUT(J)
   14 FORMAT(1H ,3I4,2F9.6,F10.6,2F8.3,F8.4,2X,6X,A4)
      IF(SQ20-T1(J)) 80,40,40
   40 DSUMA=DSUMA+ABS(DELTA)
      TSUMA=TSUMA+ABS(DELTB)
      NTALA=NTALA+1
   80 IF(SQ30-T1(J)) 81,82,82
   82 DSUMB=DSUMB+ABS(DELTA)
      TSUMB=TSUMB+ABS(DELTB)
      NTALB=NTALB+1
   81 IF(SQFI-T1(J)) 42,83,83
   83 DSUMC=DSUMC+ABS(DELTA)
      TSUMC=TSUMC+ABS(DELTB)
      NTALC=NTALC+1
   42 J=J+1
      GO TO 17
   15 HC=SQRT(T2(I))
      HD=SQRT(1.-T2(I))
      THC=GRAD*ATAN(HC/HD)
      WRITE(IOUT,16) IHT(I),IKT(I),ILT(I),T2(I),THC
   16 FORMAT(1H ,3I4,9X,F9.6,18X,F8.3)
   17 CONTINUE
   18 IF(J-N1) 19,19,21
   19 HA=SQRT(T1(J))
      HD=SQRT(1.-T1(J))
      DOBS=WAVE/(2.*HA)
      NOIND=NOIND+1
      WRITE(IOUT,20) T1(J),THO,DOBS,TAGUT(J)
      THO=GRAD*ATAN(HA/HD)
      IF(T1(J)-SAK) 114,114,115
  114 NWIND=NWIND+1
  115 J=J+1
      GO TO 18
   21 WRITE(IOUT,22) N1
   22 FORMAT(1H ,'NUMBER OF OBS. LINES =',I5)
      WRITE(IOUT,23) N2
   23 FORMAT(1H ,'NUMBER OF CALC. LINES =',I5)
      IF(NTALA) 92,92,90
   90 TMEDA=TSUMA/NTALA
      DMEDA=DSUMA/NTALA
      MER=SQ20/(2.0*DMEDA*N20)+0.5
      M20=MER
      F20=20.0/(TMEDA*N20)+0.5
      NDEF=20
      WRITE(IOUT,91) NDEF,MER,DMEDA,NDEF,F20,TMEDA,N20
   91 FORMAT(1H ,'M(',I3,')=',I4,'  AV.EPS.=',F10.7,/,
     11H ,'F',I3,' =',F5.0,'(',F8.6,',',I5,')')
   92 IF(NTALB) 94,94,93
   93 TMEDB=TSUMB/NTALB
      DMEDB=DSUMB/NTALB
      MER=SQ30/(2.0*DMEDB*N30)+0.5
      F30=30.0/(TMEDB*N30)+0.5
      NDEF=30
      WRITE(IOUT,91) NDEF,MER,DMEDB,NDEF,F30,TMEDB,N30
   94 IF(NTALC) 96,96,95
   95 DMEDC=DSUMC/NTALC
      TMEDC=TSUMC/NTALC
      FN1=N1/(TMEDC*NFI)+0.5
      MER=SQFI/(2.0*DMEDC*NFI)+0.5
      IF(M20.EQ.0) M20=MER
      WRITE(IOUT,91) N1,MER,DMEDC,N1,FN1,TMEDC,NFI
   96 WRITE(IOUT,100)
  100 FORMAT(1H ,' M   CF. J.APPL.CRYST. 1(1968)108 ',/,
     1 1H ,' F  CF. J.APPL.CRYST. 12(1979)60 ')
      WRITE(IOUT,116) NOIND
  116 FORMAT(1H ,I5,'  LINES ARE UNINDEXED ')
      NOIND=NWIND
      RETURN
      END
      SUBROUTINE PW62(A,B,C,BE,S,ISYM)
C     PARTIAL DERIVATIVES OF DIRECT LATTICE PARAMETERS(NOT TRICLINIC)
      DIMENSION S(6,6)
      GO TO(1,2,3,4,5),ISYM
    5 S(1,1)=-0.5*A*A*A
      S(1,2)=0.0
      S(1,3)=-0.5*A*C*C*COS(BE)**2
      S(1,4)=-0.5*A*A*C*COS(BE)
      S(2,1)=0.0
      S(2,2)=-0.5*B*B*B
      S(2,3)=0.0
      S(2,4)=0.0
      S(3,1)=-0.5*A*A*C*COS(BE)**2
      S(3,2)=0.0
      S(3,3)=-0.5*C*C*C
      S(3,4)=-0.5*A*C*C*COS(BE)
      S(4,1)=0.25*A*A*SIN(2.*BE)
      S(4,2)=0.0
      S(4,3)=0.25*C*C*SIN(2.*BE)
      S(4,4)=0.5*A*C*SIN(BE)
      RETURN
    4 S(1,1)=-(3.*A*A*A)/8.
    6 S(1,2)=0.0
      S(2,1)=0.0
      S(2,2)=-0.5*C*C*C
      RETURN
    3 S(1,1)=-0.5*A*A*A
      S(1,2)=0.0
      S(1,3)=0.0
      S(2,1)=0.0
      S(2,2)=-0.5*B*B*B
      S(2,3)=0.0
      S(3,1)=0.0
      S(3,2)=0.0
      S(3,3)=-0.5*C*C*C
      RETURN
    2 S(1,1)=-0.5*A*A*A
      GO TO 6
    1 S(1,1)=-0.5*A*A*A
      RETURN
      END
      SUBROUTINE PW83(A11,A22,A33,A13,A12,A23,WAVE,ISYM,A,B,C,AL,BE,GA,
     1           CEM,VOL,IOK,VTEST)
      IOK=1
C     CALC. OF DIRECT LATTICE PARAMETERS
C
      GO TO(1,2,3,4,5,6),ISYM
    1 A=WAVE/(2.*SQRT(A11))
      RETURN
    2 C=WAVE/(2.*SQRT(A22))
      GO TO 1
    3 B=WAVE/(2.*SQRT(A22))
      C=WAVE/(2.*SQRT(A33))
      GO TO 1
    4 A=WAVE/SQRT(3.*A11)
      C=WAVE/(2.*SQRT(A22))
      RETURN
    5 CBST=A13/(2.*SQRT(A11*A33))
      SBST=SQRT(1.-CBST*CBST)
      IF(CBST-0.0005) 7,8,8
    7 BE=1.570796
      GO TO 9
    8 TBST=SBST/CBST
      BST=ATAN(TBST)
      BE=3.141593-BST
    9 A=WAVE/(2.*SBST*SQRT(A11))
      B=WAVE/(2.*SQRT(A22))
      C=WAVE/(2.*SBST*SQRT(A33))
      RETURN
    6 S=4./(WAVE*WAVE)
      IF(A11) 30,30,31
   30 IOK=0
      RETURN
   31 IF(A22) 30,30,32
   32 IF(A33) 30,30,33
   33 AST=SQRT(A11*S)
      BST=SQRT(A22*S)
      CST=SQRT(A33*S)
      CAST=A23/(2.*SQRT(A22*A33))
      IF(ABS(CAST)-1.) 34,34,30
   34 CBST=A13/(2.*SQRT(A11*A33))
      IF(ABS(CBST)-1.) 35,35,30
   35 CGST=A12/(2.*SQRT(A11*A22))
      ARGUM=1.-CAST**2-CBST**2-CGST**2+2.*CAST*CBST*CGST
      IF(ARGUM) 30,36,36
   36 VST=AST*BST*CST*SQRT(ARGUM)
      VTEST=1./VST
      IF(VTEST-VOL) 37,37,30
   37 SAST=SQRT(1.-CAST**2)
      SBST=SQRT(1.-CBST**2)
      SGST=SQRT(1.-CGST**2)
      CAL=(CBST*CGST-CAST)/(SBST*SGST)
      CBE=(CAST*CGST-CBST)/(SAST*SGST)
      CGA=(CAST*CBST-CGST)/(SAST*SBST)
      IF(ABS(CAL)-0.0005) 10,11,11
   10 AL=1.570796
      GO TO 12
   11 SAL=SQRT(1.-CAL**2)
      AL=ATAN(SAL/CAL)
   12 IF(ABS(CBE)-0.0005) 13,14,14
   13 BE=1.570796
      GO TO 15
   14 SBE=SQRT(1.-CBE**2)
      BE=ATAN(SBE/CBE)
   15 IF(ABS(CGA)-0.0005) 16,17,17
   16 GA=1.570796
      GO TO 18
   17 SGA=SQRT(1.-CGA**2)
      GA=ATAN(SGA/CGA)
   18 A=(BST*CST*SAST)/VST
      IF(CEM-A) 30,38,38
   38 B=(AST*CST*SBST)/VST
      IF(CEM-B) 30,39,39
   39 C=(AST*BST*SGST)/VST
      IF(CEM-C) 30,40,40
   40 RETURN
      END
      SUBROUTINE TREOD
      DIMENSION         T2(270),IHT(270),IKT(270),ILT(270),IQE(270),AMAT
     1(6,6),BMAT(6,1),DET(6,6),S(6,6),F(6)
      COMMON T1(100),TAGUT(100),N1,D1,SSQTL,D2,WAVE,IHM,IKM,ILM,ISYM,INR
     1,ISUB,A11,A22,A33,A13,A12,A23,LIN(20),REV(20)
     2,IGOT,IS1,IS2,IS3,IS4,IA,IB,IQ,IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3
     3,SQ1,SQ2,SQ3,X,Y,Z,IC(20),ISQM,KOTA(30),RO(7,30),ID,IIN,IOUT
      COMMON NH1,NK1,NL1,NS1,NH2,NK2,NL2,NS2,NH3,NK3,NL3,NS3,IUTSK,IMER,
     1INIX,IOMIT,ISC,KLART,VOL,VTEST,XL,YL,IT20,IT40,CEM
     2,IFRED(40)
      ISYFL=1
      LIN(3)=1
  720 TRIAL=0.00015
      IHALT=0
      RKON=1.0/6.0
      IF(ISYM .NE. 5) GO TO 731
      IF(ABS(A13) .GT. D1) GO TO 731
      A13=0.0
      ISYLF=2
      ISAVE=ISYM
      ISYM=3
  731 A=0.
      B=0.
      C=0.
      AL=0.
      BE=0.
      GA=0.
      DA=0.
      DAL=0.
      DB=0.
      DBE=0.
      DC=0.
      DGA=0.
      VOLM=0.
      NEKV=0
    7 FORMAT(6F10.6)
 2000 WRITE(IOUT,8)
    8 FORMAT(///,14H CYCLE RESULTS,//)
      IV=INR
      IDIR=INR
      IHM=0
      IKM=0
      ILM=0
      SSQTM=T1(N1)+2*D2
   76 IHM=IHM+1
      IF(SSQTM-IHM*IHM*A11)77,76,76
   77 IHM=IHM-1
      GO TO(101,102,78,102,78,78),ISYM
  101 IKM=IHM
      ILM=IHM
      GO TO 81
  102 IKM=IHM
  103 ILM=ILM+1
      IF(SSQTM-ILM*ILM*A22)81,103,103
   78 IKM=IKM+1
      IF(SSQTM-IKM*IKM*A22)79,78,78
   79 IKM=IKM-1
   80 ILM=ILM+1
      IF(SSQTM-ILM*ILM*A33)81,80,80
   81 IF(N1-30) 304,305,305
  305 NT30=30
      SQ30=T1(30)
      GO TO 311
  304 NT30=0
      SQ30=0.0
      IF(N1-20) 312,311,311
  311 NT20=20
      SQ20=T1(20)
      GO TO 200
  312 NT20=0
      SQ20=0.0
  200 N2=0
      IH=0
      IK=0
      IL=0
      N20=0
      N30=0
      NFI=0
      SQFI=T1(N1)
  300 CALL PW50(A11,A22,A33,A13,A12,A23,ISYM,IHM,IKM,ILM,IH,IK,IL,SSQT,I
     1A)
      IF(IA) 301,600,301
  301 IF(SSQT-SSQTM) 303,303,300
  303 CALL PW51(SSQT,T1,N1,SSQTL,D1,D2,SSQTO,IA)
      IF(SSQT-SQFI) 307,307,308
  307 CALL PW80(SSQT,T2,N2,IHT,IKT,ILT,IH,IK,IL,ISYM,INNE)
      IF(INNE) 1051,1051,316
 1051 NFI=NFI+1
C 307 NFI=NFI+1
  308 IF(SSQT-SQ30) 313,313,314
  313 N30=N30+1
  314 IF(SSQT-SQ20) 315,315,316
  315 N20=N20+1
  316 IF(IA) 401,300,401
  401 IF(INR) 500,402,500
  402 CALL PW53(SSQT,T2,N2,IHT,IKT,ILT,IQE,IH,IK,IL,ISYM)
      IF(SSQT) 403,300,300
  403 N2=N2+1
      IF(N2-270) 300,900,900
  500 CALL PW53(SSQTO,T2,N2,IHT,IKT,ILT,IQE,IH,IK,IL,ISYM)
      IF(SSQTO) 501,300,300
  501 N2=N2+1
      IF(N2-270) 300,600,600
  600 IF(INR) 601,900,601
  601 IF(IV) 602,603,602
  602 IV=IV-1
      GO TO 700
  603 CALL PW71(T2,IHT,IKT,ILT,IQE,ISYM,N2,A11,A22,A33,A13,A12,A23,V2,NE
     1KV)
      INR=0
      CALL PW83(A11,A22,A33,A13,A12,A23,WAVE,ISYM,A,B,C,AL,BE,GA,
     1          CEM,VOL,IOK,VVV)
      IF(ISYM-6) 610,609,610
  609 IF(IOK) 2520,1204,2520
 2520 CALL PW61(A,B,C,AL,BE,GA,S)
      GO TO 611
  610 CALL PW62(A,B,C,BE,S,ISYM)
  611 CALL PW52(N,NEKV,ISYM,V2,WAVE,AMAT,S,A,B,C,AL,BE,GA,DA,DB,DC,
     1DAL,DBE,DGA,VOLM)
      WRITE(IOUT,641)NEKV,N1
  641 FORMAT(//,32H NUMBER OF SINGLE INDEXED LINES=,I5,23H TOTAL NUMBER
     1OF LINES=,I5)
      GO TO 200
  700 CALL PW55(T2,IQE,IHT,IKT,ILT,N2,ISYM,AMAT,BMAT,IEKV)
      IF(IEKV) 800,701,800
  701 WRITE(IOUT,702)
  702 FORMAT(44H NUMBER OF EQ LESS THAN NUMBER OF PARAMETERS)
      INR=0
      GO TO 200
  800 N=1
      GO TO (801,802,803,802,804,805),ISYM
  805 N=N+2
  804 N=N+1
  803 N=N+1
  802 N=N+1
  801 CALL PW84(AMAT,N,BMAT,ISOLB)
      IF(ABS(A11-BMAT(1,1))-0.001) 910,910,912
  910 IF(ABS(A22-BMAT(2,1))-0.001) 911,911,912
  912 INR=0
      GO TO 200
  911 A11=BMAT(1,1)
      A22=BMAT(2,1)
      A33=BMAT(3,1)
      A13=BMAT(4,1)
      A12=BMAT(5,1)
      A23=BMAT(6,1)
      WRITE(IOUT,7) A11,A22,A33,A13,A12,A23
      GO TO 200
  900 IF(NEKV .NE. 0) GO TO 1300
      CALL PW83(A11,A22,A33,A13,A12,A23,WAVE,ISYM,A,B,C,AL,BE,GA,
     1          CEM,VOL,IOK,VVV)
      RAD=0.0174533
      GO TO(2501,2501,2501,2502,2504,2506),ISYM
 2501 GA=1.5707963
 2503 BE=1.5707963
 2505 AL=1.5707963
      GO TO 2506
 2502 GA=2.0943951
      GO TO 2503
 2504 GA=1.5707963
      GO TO 2505
 2506 CAL=COS(AL)
      CBE=COS(BE)
      CGA=COS(GA)
      VOLM=A*B*C*SQRT(1.-CAL**2-CBE**2-CGA**2+2.*CAL*CBE*CGA)
      AL=AL/RAD
      BE=BE/RAD
      GA=GA/RAD
      IF(AL) 1301,1301,1302
 1301 AL=AL+180.
 1302 IF(BE) 1303,1303,1304
 1303 BE=BE+180.
 1304 IF(GA) 1305,1305,1306
 1305 GA=GA+180.
 1306 WRITE(IOUT,1307)
 1307 FORMAT(' NOT REFINED UNIT CELL ')
      WRITE(IOUT,1308) A,B,C,AL,BE,GA,VOLM
 1308 FORMAT('  A=',F7.3,'  B=',F7.3,'  C=',F7.3,/,' AL=',F7.3,
     1' BE=',F7.3,' GA=',F7.3,'  CELL VOLUME=',F7.2,' A**3')
      GO TO 1310
 1300 WRITE(IOUT,650) NEKV,N1,A,DA,AL,DAL,B,DB,BE,DBE,C,DC,GA,DGA,VOLM
  650 FORMAT(1X,' NUMBER OF SINGLE INDEXED LINES =',I5,/,1X,
     1'  TOTAL NUMBER OF LINES =',I5,/,1X,
     2'  A =',2F10.6,' A   ALFA =',2F10.6,' DEG',/,1X,
     3'  B =',2F10.6,' A   BETA =',2F10.6,' DEG',/,1X,
     4'  C =',2F10.6,' A  GAMMA =',2F10.6,' DEG',/,1X,
     5'  UNIT CELL VOLUME =',F10.2,' A**3')
 1310 CALL P56W(T1,T2,IHT,IKT,ILT,N1,N2,IOUT,D1,SSQTL,D2,WAVE,TAGUT,
     1N20,N30,NFI,NT20,NT30,SQ20,SQ30,SQFI,M20,NOIND)
      IF(ISYM-3) 1026,1030,1026
 1030 IF(M20-IMER) 1200,1201,1201
 1201 IF(NOIND-INIX) 1202,1202,1200
 1202 IHALT=1
 1200 IF(ABS(A11-A22)-TRIAL) 1001,1001,1002
 1001 IF(ABS(A11-A33)-TRIAL) 1003,1003,1004
 1003 A11=(A11+A22+A33)/3.0
      WRITE(IOUT,1005)
 1005 FORMAT(23H TEST OF CUBIC SYMMETRY)
      A22=0.0
      A33=0.0
      ISYM=1
 1007 INR=3
      GO TO 2000
 1004 A11=(A11+A22)/2.0
      A22=A33
 1010 A33=0.0
      ISYM=2
      WRITE(IOUT,1006)
 1006 FORMAT(28H TEST OF TETRAGONAL SYMMETRY)
      GO TO 1007
 1002 IF(ABS(A11-A33)-TRIAL) 1008,1008,1009
 1008 A11=(A11+A33)/2.0
      GO TO 1010
 1009 IF(ABS(A22-A33)-TRIAL) 1011,1011,1012
 1011 HA22=A11
      A11=(A22+A33)/2.0
      A22=HA22
      GO TO 1010
 1012 IF(ABS(A11-3.0*A22)-TRIAL) 1013,1013,1014
 1013 A11=(A11+3.0*A22)*RKON
 1018 A22=A33
 1020 A33=0.0
      ISYM=4
      WRITE(IOUT,1015)
 1015 FORMAT(27H TEST OF HEXAGONAL SYMMETRY)
      GO TO 1007
 1014 IF(ABS(A22-3.0*A11)-TRIAL) 1016,1016,1017
 1016 A11=(A22+3.0*A11)*RKON
      GO TO 1018
 1017 IF(ABS(A11-3.0*A33)-TRIAL) 1027,1027,1019
 1027 A11=(A11+3.0*A33)*RKON
      GO TO 1020
 1019 IF(ABS(A33-3.0*A11)-TRIAL) 1021,1021,1022
 1021 A11=(A33+3.0*A11)*RKON
      GO TO 1020
 1022 IF(ABS(A22-3.0*A33)-TRIAL) 1023,1023,1024
 1023 HA22=A11
      A11=(A22+3.0*A33)*RKON
      A22=HA22
      GO TO 1020
 1024 IF(ABS(A33-3.0*A22)-TRIAL) 1025,1025,1026
 1025 HA22=A11
      A11=(A33+3.0*A22)*RKON
      A22=HA22
      GO TO 1020
 1026 IF(M20-IMER) 1100,1101,1101
 1101 IF(NOIND-INIX) 1102,1102,1100
 1102 WRITE(IOUT,1103)
 1103 FORMAT(59H CHECK IF THERE IS ANY COMMON FACTOR IN THE QUADRATIC FO
     1RMS,/,69H CHECK CONVERGENCE IN THE REFINEMENT (EV. USE PROGRAM PIR
     2UM OR PURUM),/,20H END OF CALCULATIONS)
 1203 LIN(3)=2
      RETURN
 1100 IF(IHALT) 1204,1204,1203
 1204 GO TO(711,712),ISYFL
  712 ISYM=ISAVE
  711 RETURN
      END
      SUBROUTINE PW52(N,NEKV,ISYM,V2,WAVE,AMAT,S,A,B,C,AL,BE,GA,DA,DB,DC
     1,DAL,DBE,DGA,VOLM)
      DIMENSION AMAT(6,6),S(6,6),F(6)
      XXX=V2/(FLOAT(NEKV-N))
      IF(XXX .LT. 0.) XXX=0.
      SIG=SQRT(XXX)
C      SIG=SQRT(V2/(FLOAT(NEKV-N)))
      SIG=(4.0*SIG)/(WAVE**2)
      DO 2 L=1,N
      F(L)=0.0
      DO 1 J=1,N
      DO 1 K=1,N
    1 F(L)=F(L)+AMAT(J,K)*S(L,J)*S(L,K)
      IF(F(L) .LE. 0.) F(L)=0.
    2 F(L)=SIG*SQRT(F(L))
      DB=0.0
      DC=0.0
      DAL=0.0
      DBE=0.0
      DGA=0.0
      RAD=0.0174533
      GRAD=57.2958
      DA=F(1)
      GO TO(3,4,5,6,7,8),ISYM
    3 B=A
      C=A
      DB=DA
      DC=DA
      GO TO 15
    4 B=A
      DB=DA
      DC=F(2)
      GO TO 15
    5 DB=F(2)
      DC=F(3)
      GO TO 15
    6 DB=DA
      B=A
      DC=F(2)
      GA=120.0
      GO TO 16
    7 DB=F(2)
      DC=F(3)
      BE=GRAD*BE
      DBE=GRAD*F(4)
      GO TO 17
    8 DB=F(2)
      DC=F(3)
      AL=GRAD*AL
      IF(AL) 9,10,10
    9 AL=AL+180.0
   10 DAL=GRAD*F(4)
      BE=GRAD*BE
      IF(BE) 11,12,12
   11 BE=BE+180.0
   12 DBE=GRAD*F(5)
      GA=GRAD*GA
      IF(GA) 13,14,14
   13 GA=GA+180.0
   14 DGA=GRAD*F(6)
      GO TO 18
   15 VOLM=A*B*C
      BE=90.0
   19 GA=90.0
   20 AL=90.0
      GO TO 999
   16 VOLM=0.8660254*A*B*C
      BE=90.0
      GO TO 20
   17 VOLM=A*B*C*SIN(RAD*BE)
      GO TO 19
   18 CA=COS(RAD*AL)
      CB=COS(RAD*BE)
      CC=COS(RAD*GA)
      HABC=1.0-CA**2-CB**2-CC**2+2.0*CA*CB*CC
      IF(HABC .LE. 0.) HABC=0.
      VOLM=A*B*C*SQRT(HABC)
  999 RETURN
      END
      SUBROUTINE PW80(SSQT,T2,N2,IHT,IKT,ILT,IH,IK,IL,ISYM,INNE)
      DIMENSION T2(1),IHT(1),IKT(1),ILT(1)
      I=1
    5 IF(I-N2) 2,2,1
    1 INNE=0
      RETURN
    2 IF(SSQT-T2(I)) 3,4,3
    3 I=I+1
      GO TO 5
    4 GO TO (11,12,13,14,15,16),ISYM
   11 IF((IH**2+IK**2+IL**2)-(IHT(I)**2+IKT(I)**2+ILT(I)**2)) 7,6,7
    7 INNE=0
      RETURN
   12 IF((IH**2+IK**2)-(IHT(I)**2+IKT(I)**2)) 7,21,7
   21 IF(IL**2-ILT(I)**2) 7,6,7
   13 IF(IH**2-IHT(I)**2) 7,22,7
   22 IF(IK**2-IKT(I)**2) 7,21,7
   14 IF((IH**2+IH*IK+IK**2)-(IHT(I)**2+IHT(I)*IKT(I)+ILT(I)**2)) 7,21,7
   15 IF(IH*IL-IHT(I)*ILT(I)) 7,13,7
   16 IF(IH*IK-IHT(I)*IKT(I)) 7,23,7
   23 IF(IK*IL-IKT(I)*ILT(I)) 7,15,7
    6 INNE=1
      RETURN
      END
      SUBROUTINE PDIV(T1,N1,D1,SSQTL,D2,IDIV)
      DIMENSION T1(100)
      IF(IDIV) 1,2,1
    1 DO 30 ND=1,2
      NDI=(ND+1)
      RDI=NDI*NDI
      DO 10 N=1,7
      SQ=T1(N)
      ISTAR=N+1
      TEST=SQ*RDI
      IF(TEST-SSQTL) 3,4,4
    3 DT=D1
      GO TO 5
    4 DT=D2
    5 DO 20 I=ISTAR,N1
      DIFF=ABS(T1(I)-TEST)
      IF(DIFF-DT) 6,6,20
    6 T1(N)=T1(I)/RDI
      GO TO 10
   20 CONTINUE
   10 CONTINUE
   30 CONTINUE
    2 RETURN
      END
      SUBROUTINE PW61(A,B,C,AL,BE,GA,S)
C     TRICLINIC SYM. PARTIAL DERIVATIVES OF DIRECT LATTICE PARAMETERS
      DIMENSION S(6,6)
      S(1,1)=-0.5*A*A*A
      S(1,2)=-0.5*A*B*B*COS(GA)**2
      S(1,3)=-0.5*A*C*C*COS(BE)**2
      S(1,4)=-0.5*A*B*C*COS(BE)*COS(GA)
      S(1,5)=-0.5*A*A*C*COS(BE)
      S(1,6)=-0.5*A*A*B*COS(GA)
      S(2,1)=-0.5*B*A*A*COS(GA)**2
      S(2,2)=-0.5*B*B*B
      S(2,3)=-0.5*B*C*C*COS(AL)**2
      S(2,4)=-0.5*B*B*C*COS(AL)
      S(2,5)=-0.5*A*B*C*COS(AL)*COS(GA)
      S(2,6)=-0.5*B*B*A*COS(GA)
      S(3,1)=-0.5*C*A*A*COS(BE)**2
      S(3,2)=-0.5*C*B*B*COS(AL)**2
      S(3,3)=-0.5*C*C*C
      S(3,4)=-0.5*C*C*B*COS(AL)
      S(3,5)=-0.5*C*C*A*COS(BE)
      S(3,6)=-0.5*A*B*C*COS(AL)*COS(BE)
      US=2.*COS(BE)*COS(GA)-COS(AL)*(COS(BE)**2+COS(GA)**2)
      S(4,1)=(A*A*US)/(2.*SIN(AL))
      S(4,2)=0.25*B*B*SIN(2.*AL)
      S(4,3)=0.25*C*C*SIN(2.*AL)
      S(4,4)=0.5*B*C*SIN(AL)
      S(4,5)=0.5*A*C*SIN(AL)*COS(GA)
      S(4,6)=0.5*A*B*SIN(AL)*COS(BE)
      S(5,1)=0.25*A*A*SIN(2.*BE)
      US=2.*COS(AL)*COS(GA)-COS(BE)*(COS(AL)**2+COS(GA)**2)
      S(5,2)=(B*B*US)/(2.*SIN(BE))
      S(5,3)=0.25*C*C*SIN(2.*BE)
      S(5,4)=0.5*B*C*SIN(BE)*COS(GA)
      S(5,5)=0.5*A*C*SIN(BE)
      S(5,6)=0.5*A*B*SIN(BE)*COS(AL)
      S(6,1)=0.25*A*A*SIN(2.*GA)
      S(6,2)=0.25*B*B*SIN(2.*GA)
      US=2.*COS(AL)*COS(BE)-COS(GA)*(COS(AL)**2+COS(BE)**2)
      S(6,3)=(C*C*US)/(2.*SIN(GA))
      S(6,4)=0.5*B*C*SIN(GA)*COS(BE)
      S(6,5)=0.5*A*C*SIN(GA)*COS(AL)
      S(6,6)=0.5*A*B*SIN(GA)
      RETURN
      END
      SUBROUTINE TREOB
C
C
      COMMON T1(100),TAGUT(100),N1,D1,SSQTL,D2,WAVE,IHM,IKM,ILM,ISYM,INR
     1,ISUB,A11,A22,A33,A13,A12,A23,LIN(20),REV(20)
     2,IGOT,IS1,IS2,IS3,IS4,IA,IB,IQ,IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3
     3,SQ1,SQ2,SQ3,X,Y,Z,IC(20),ISQM,KOTA(30),RO(7,30),ID,IIN,IOUT
C
      COMMON NH1,NK1,NL1,NS1,NH2,NK2,NL2,NS2,NH3,NK3,NL3,NS3,IUTSK,IMER,
     1INIX,IOMIT,ISC,KLART,VOL,VTEST,XL,YL,IT20,IT40,CEM
     2,KHA,KKA,KLA,KSA,KHB,KKB,KLB,KSB,MONO,MH1,MK1,ML1,MH2,MK2,ML2,MH3,
     3MK3,ML3,MS1,MS2,MS3,MH4,MK4,ML4,MS4,UL,LGO,MG,MS,MRA,MRB,NU,
     4DEN,EDEN,RMOLV
C
      COMMON VOL1(100),VOL2(100),AA(6,6),BB(6,1),SQQ(5),SLD(5),HL(5),
     1IDX(5),XX(5),SQOLD(5),AAA(5,5),BBB(5,1),VW(5)
      RAD=0.0174533
      NDISP=LIN(11)
      VTEST=REV(11)
      VSK=(WAVE**6)/64.0
      XL=REV(12)
      TRIT=(WAVE/(2.*25.))**2
      ITLIM=SSQTL*100000.0+0.5
      IZTAL=LIN(14)
      ITRIC=LIN(8)
      ISHORT=LIN(9)
      IF(ITRIC) 4000,3000,4000
 3000 GO TO(110,310,330,320,200),ISYM
C
  400 FORMAT(25H BASE LINE ONE.(HKL)-MAX=,3I5,11H MAX H+K+L=,I5)
  401 FORMAT(25H BASE LINE TWO.(HKL)-MAX=,3I5,11H MAX H+K+L=,I5)
  402 FORMAT(27H BASE LINE THREE.(HKL)-MAX=,3I5,11H MAX H+K+L=,I5)
  403 FORMAT(28H SELECTED BASE LINES (1) (2))
  404 FORMAT(38H SELECTED BASE LINES (1,2) (1,3) (2,3))
  405 FORMAT(69H SELECTED BASE LINES (1,2,3) (1,2,4) (1,2,5) (1,3,4) (2,
     13,4) (1,2,6) )
  406 FORMAT(50H SELECTED BASE LINES (1,2,3,4) (1,2,3,5) (1,2,4,5) )
  407 FORMAT(26H BASE LINE FOUR.(HKL)-MAX=,3I5,11H MAX H+K+L=,I5)
  410 FORMAT(22H SELECTED BASE LINES (,I2,7H,1,2) (,I2,7H,1,3) (,I2,5H,2
     1,3))
C     TEST CUBIC SYMMETRY
  110 NOW=0
C
      VT=VTEST**0.3333
      IF(VT-XL) 361,362,362
  361 AL=XL
      GO TO 363
  362 AL=VT
  363 WRITE(IOUT,350)
  350 FORMAT(11H CUBIC TEST)
      WRITE(IOUT,403)
      WRITE(IOUT,400) KHB,KKB,KLB,KSB
  111 NOW=NOW+1
      GO TO(112,113,210),NOW
  112 SQ1=T1(1)
      GO TO 120
  113 SQ1=T1(2)
  120 IH1=0
      IK1=0
      IL1=0
  121 CALL GHKLA(KHB,KKB,KLB,KSB,0,IEND,IH1,IK1,IL1)
      IF(IEND) 122,111,122
  122 CALL EKVA(IH1,IK1,IL1,0,0,0,1,ISOLV,SQ1,0.0,X,0.0,AL,0.0)
      IF(ISOLV) 123,121,123
  123 CALL TETAL(T1,X,0.0,IB,IQ,ISQM,K,IC,ISYM,IT20,IT40,ITLIM)
      IF(K-IQ) 121,124,124
  124 IF(IUTSK) 125,127,125
  125 WRITE(NDISP,126) K,X
  126 FORMAT(14H CUBIC TEST.K=,I3,3H X=,F10.6)
  127 CALL LAGRA(K)
      GO TO 121
C     TEST TETRAGONAL SYMMETRY
  310 YL=XL
      WRITE(IOUT,351)
  351 FORMAT(16H TETRAGONAL TEST)
      VTYST=VTEST
      GO TO 150
C     TEST HEXAGONAL SYMMETRY
  320 XL=(4.0*YL)/3.0
      WRITE(IOUT,352)
  352 FORMAT(15H HEXAGONAL TEST)
      VTYST=(VTEST*16.0)/9.0
  150 NOW=0
      WRITE(IOUT,404)
      WRITE(IOUT,400) KHA,KKA,KLA,KSA
      WRITE(IOUT,401) KHA,KKA,KLA,KSA
      ID=0
  151 NOW=NOW+1
      GO TO(152,153,154,210),NOW
  152 SQ1=T1(1)
      SQ2=T1(2)
      GO TO 160
  153 SQ2=T1(3)
      GO TO 160
  154 SQ1=T1(2)
  160 IH1=0
      IK1=0
      IL1=0
  161 CALL GHKLA(KHA,KKA,KLA,KSA,1,IEND,IH1,IK1,IL1)
      IF(IEND) 162,151,162
  162 IH2=0
      IK2=0
      IL2=0
      ITTE=IH1*IH1+IK1*IK1
      ITTH=ITTE+IH1*IK1
  163 CALL GHKLA(KHA,KKA,KLA,KSA,1,IEND,IH2,IK2,IL2)
      IF(IEND) 164,161,164
  164 IF(IL1-IL2) 168,165,165
  165 GO TO(200,166,200,167),ISYM
  166 IF(ITTE-IH2*IH2-IK2*IK2) 168,163,163
  167 IF(ITTH-IH2*IH2-IK2*IK2-IH2*IK2) 168,163,163
  168 CALL EKVA(IH1,IK1,IL1,IH2,IK2,IL2,ISYM,ISOLV,SQ1,SQ2,X,Y,XL,YL)
      IF(ISOLV) 169,163,169
  169 IF(X*X*Y-VTYST) 163,163,170
  170 CALL TETAL(T1,X,Y,IB,IQ,ISQM,K,IC,ISYM,IT20,IT40,ITLIM)
      IF(K-IQ) 163,171,171
  171 IF(IUTSK) 172,177,172
  172 GO TO(200,173,200,175),ISYM
  173 WRITE(NDISP,174) K,X,Y
  174 FORMAT(19H TETRAGONAL TEST.K=,I3,4H XY=,2F10.6)
      GO TO 177
  175 WRITE(NDISP,176) K,X,Y
  176 FORMAT(18H HEXAGONAL TEST.K=,I3,4H XY=,2F10.6)
  177 CALL LAGRA(K)
      GO TO 163
C     TEST ORTHOROMBIC SYMMETRY
  330 WRITE(IOUT,353)
  353 FORMAT(18H ORTHORHOMBIC TEST)
      IF(ISC) 411,412,411
  411 WRITE(IOUT,410) ISC,ISC,ISC
      GO TO 413
  412 WRITE(IOUT,405)
  413 WRITE(IOUT,400) NH1,NK1,NL1,NS1
      WRITE(IOUT,401) NH2,NK2,NL2,NS2
      WRITE(IOUT,402) NH3,NK3,NL3,NS3
      ID=0
      VTEST=(WAVE**6)/(64.0*VOL*VOL)
      NU=1
    5 CALL SEL(IS1,IS2,IS3,NU,ISC)
      IF(NU) 210,210,7
C
C     NEW HKL ON FIRST LINE
    7 IH1=0
      IK1=0
      IL1=0
      SQ1=T1(IS1)
      SQ2=T1(IS2)
      SQ3=T1(IS3)
    8 CALL GHKLA(NH1,NK1,NL1,NS1,0,IEND,IH1,IK1,IL1)
      IF(IEND) 9,5,9
    9 IH2=0
      IK2=0
      IL2=0
      IO1=1
C     NEW HKL ON SECOND LINE
   10 CALL GEHKL(NH2,NK2,NL2,NS2,1,ISLUT,IH2,IK2,IL2)
      IF(ISLUT) 11,8,11
   11 IF(IH1-IK1) 21,20,21
   20 IF(IH2-IK2) 10,21,21
   21 IF(IK1-IL1) 23,22,23
   22 IF(IK2-IL2) 10,23,23
   23 IF(IH1-IH2) 60,61,61
   61 IF(IK1-IK2) 60,62,62
   62 IF(IL1-IL2) 60,10,10
   60 IH3=0
      IK3=0
      IL3=0
      IO2=1
C     NEW HKL ON THIRD LINE
    3 CALL GEHKL(NH3,NK3,NL3,NS3,1,ISLUT,IH3,IK3,IL3)
       IF(ISLUT) 65,10,65
   65 IF(IH2-IH3) 13,66,66
   66 IF(IK2-IK3) 13,67,67
   67 IF(IL2-IL3) 13,3,3
C     CALL XYZ ORTHOR. SYM.
   13 CALL ORT(IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3,ISOLV,SQ1,SQ2,SQ3,X,Y
     1,Z,IO1,IO2,XL)
      IF(ISOLV) 53,3,53
   53 RVA=X*Y*Z
      IF(VTEST-RVA)480,480,3
  480 IF(DEN-0.00001)14,14,479
  479 DO 476 I=1,IZTAL
      IF(RVA-VOL1(I))476,478,478
  478 IF(RVA-VOL2(I))14,14,3
  476 CONTINUE
      GO TO 3
C     CALC. NUMBER OF INDEXABLE LINES
   14 CALL ORTAL(T1,X,Y,Z,N1,IB,IQ,ISQM,K,IC,IT20,IT40,ITLIM)
      IF(K-IQ) 3,15,15
   15 IF(IUTSK) 17,16,17
   17 WRITE(NDISP,101) K,X,Y,Z
C
  101 FORMAT(20H ORTHOROMBIC TEST.K=,I3,5H XYZ=,3F10.6)
   16 CALL LAGRA(K)
      GO TO 3
  210 IGOT=0
      MRA=ID
      GO TO 900
C     TEST MONOCLINIC SYMMETRY
  200 ID=0
  606 LIN(15)=LIN(15)+1
      MSW=LIN(15)
      NFLGA=0
      IF(MSW .EQ. 1) GO TO 3005
      IF(MSW .EQ. 2) GO TO 603
      IF(MONO) 541,541,1059
 1059 IF(NU) 542,542,540
  542 NU=1
      GO TO 625
  603 MMSW=0
      IF(MONO) 1071,541,1070
 1071 MMSW=1
 1070 MONO=IABS(MONO)
      IF(ISHORT) 3001,606,3001
 3001 CALL TWODIM
      IF(MMSW .EQ. 1) MONO=-MONO
      LIN(9)=0
      GO TO 210
 3005 IF(MONO .EQ. 0) GO TO 606
      WRITE(IOUT,511)
  511 FORMAT(16H MONOCLINIC TEST)
      WRITE(IOUT,512) MONO
  512 FORMAT(18H MAX BETA ALLOWED=,I5,5H DEG.)
      IF(MONO .LT. 0) GO TO 606
  604 FORMAT(13H (020)-SEARCH)
 1058 CONTINUE
      WRITE(IOUT,604)
      WRITE(NDISP,604)
C
C BEGIN OF (020)-ZONE TEST
C
      DO 1002 NIA=1,IB
      SQ1=T1(NIA)
      TSQ1=2.0*SQ1
      DO 1003 NIB=1,IB
      IF (NIA-NIB) 1004,1003,1004
 1004 SQ2=T1(NIB)
      TSQ12=TSQ1+SQ2
      DO 1005 NIC=1,IB
      IF (NIA-NIC) 1006,1005,1006
 1006 IF (NIB-NIC) 1007,1005,1007
 1007 SQ3=T1(NIC)
      IF (ABS(TSQ12-SQ3)-D2) 1010,1010,1005
 1005 CONTINUE
 1003 CONTINUE
 1002 CONTINUE
      WRITE(NDISP,1031)
 1031 FORMAT(20H END OF (020)-SEARCH)
      IF(NFLGA) 606,606,210
 1010 Y=SQ1/4.0
      DO 1011 NID=1,3
      RNID=NID*NID
      X=(SQ2-Y)/RNID
      NUU=0
      NUUU=0
      IF (X) 1011,1011,1020
 1020 NUU=NUU+1
      NSQA=ISERCH(T1,X,Y,0.0,NUU,IB,D1,D2,SSQTL)
      IF(NSQA .EQ. NIC) GO TO 1020
      NSQAA=NSQA
 1080 NSQAA=NSQAA+1
      NSQB=ISERCH(T1,X,Y,0.0,NSQAA,IB,D1,D2,SSQTL)
      IF(NSQB .EQ. NIC) GO TO 1080
      NUUU=NUUU+1
      IF(NUUU .EQ. 4) GO TO 1011
 1021 SQAA=T1(NSQA)
      SQBB=T1(NSQB)
      IH1=0
      IK1=0
      IL1=0
 1022 CALL GEHKL (2,2,2,3,2,ISLUT,IH1,IK1,IL1)
      IF(ISLUT) 1023,1080,1023
 1023 EONE=SQAA-IH1*IH1*X-IK1*IK1*Y
      IH2=0
      IK2=0
      IL2=0
 1024 CALL GEHKL (2,2,2,4,2,ISLUT,IH2,IK2,IL2)
      IF (ISLUT) 1025,1022,1025
 1025 ETWO=SQBB-IH2*IH2*X-IK2*IK2*Y
      IDETER=IL1*IL1*IH2*IL2-IH1*IL1*IL2*IL2
      IF (IDETER) 1026,1024,1026
 1026 DETER=IDETER
      Z=(EONE*IH2*IL2-ETWO*IH1*IL1)/DETER
      IF (Z) 1024,1024,1040
 1040 U=(IL1*IL1*ETWO-IL2*IL2*EONE)/DETER
      IF (U) 1024,1061,1061
 1061 HALP=X*Z
      UTE=(U*U)/HALP
      IF (UTE-UL) 1041,1041,1024
 1041 HELP=1.0-0.25*UTE
      VST2=Y*HALP*HELP
      IF (VTEST-VST2) 1042,1042,1024
 1042 IF (X*HELP-XL) 1024,1043,1043
 1043 IF (Z*HELP-XL) 1024,1044,1044
 1044 IF (DEN-0.00001) 1050,1050,1045
 1045 DO 1046 I=1,IZTAL
      IF (VST2-VOL1(I)) 1046,1047,1047
 1047 IF (VST2-VOL2(I)) 1050,1050,1024
 1046 CONTINUE
      GO TO 1024
 1050 CALL MAEG (T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF (K-IQ) 1024,1027,1027
 1027 LIN(4)=LIN(4)+1
      IF (IUTSK) 1028,1030,1028
 1028 WRITE (NDISP,1029) K,X,Y,Z,U
 1029 FORMAT (20H MONOCLINIC (-20) K=,I3,6H XYZU=,4F10.6)
 1030 IF(K-IB+1) 1056,1055,1055
 1055 NFLGA=1
 1056 CALL LAGRM (K,U)
      GOTO 1024
 1011 CONTINUE
      GOTO 1005
C
C END OF (020)-TEST
C
  625 IF(MONO) 541,541,1060
 1060 IF(ISC-5) 515,515,513
  513 WRITE(IOUT,514) ISC,ISC,ISC
  514 FORMAT(22H SELECTED BASE LINES (,I2,9H,1,2,3) (,I2,9H,1,2,4) (,I2,
     17H,1,3,4))
      GO TO 516
  515 WRITE(IOUT,406)
  516 WRITE(IOUT,400) MH1,MK1,ML1,MS1
      WRITE(IOUT,401) MH2,MK2,ML2,MS2
      WRITE(IOUT,402) MH3,MK3,ML3,MS3
      WRITE(IOUT,407) MH4,MK4,ML4,MS4
      GO TO 540
  517 IF(MG) 210,540,210
  540 CALL SELE(IS1,IS2,IS3,IS4,NU,ISC,MS,IOUT)
      IF(NU) 541,541,518
  541 MG=0
      GO TO 210
  518 SQ1=T1(IS1)
      SQ2=T1(IS2)
      SQ3=T1(IS3)
      SQ4=T1(IS4)
      IF(SQ1-SQ2) 590,591,591
  590 INIMP=1
      GO TO 592
  591 INIMP=0
  592 IH1=0
      IK1=0
      IL1=0
  519 CALL GEHKL(MH1,MK1,ML1,MS1,2,ISLUT,IH1,IK1,IL1)
      IF(ISLUT) 520,517,520
  520 IF(IH1) 522,521,521
  521 IF(IH1-IL1) 519,522,522
  522 IH2=0
      IK2=0
      IL2=0
      IO1=1
  523 CALL GEHKL(MH2,MK2,ML2,MS2,2,ISLUT,IH2,IK2,IL2)
      IF(ISLUT) 524,519,524
  524 IF(INIMP) 593,580,593
  593 CALL IMP(IH1,IK1,IL1,IH2,IK2,IL2,IOK)
      IF(IOK) 580,523,580
  580 IH3=0
      IK3=0
      IL3=0
      IO2=1
  525 CALL GEHKL(MH3,MK3,ML3,MS3,2,ISLUT,IH3,IK3,IL3)
      IF(ISLUT) 526,523,526
  526 CALL IMP(IH2,IK2,IL2,IH3,IK3,IL3,IOK)
      IF(IOK) 581,525,581
  581 IF(INIMP) 594,582,594
  594 CALL IMP(IH1,IK1,IL1,IH3,IK3,IL3,IOK)
      IF(IOK) 582,525,582
  582 IH4=0
      IK4=0
      IL4=0
      IO3=1
  527 CALL GEHKL(MH4,MK4,ML4,MS4,2,ISLUT,IH4,IK4,IL4)
      IF(ISLUT) 583,525,583
  583 CALL IMP(IH3,IK3,IL3,IH4,IK4,IL4,IOK)
      IF(IOK) 584,527,584
  584 CALL IMP(IH2,IK2,IL2,IH4,IK4,IL4,IOK)
      IF(IOK) 585,527,585
  585 IF(INIMP) 595,528,595
  595 CALL IMP(IH1,IK1,IL1,IH4,IK4,IL4,IOK)
      IF(IOK)528,527,528
  528 CALL MOC(IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3,IH4,IK4,IL4,SQ1,SQ2,S
     1Q3,SQ4,X,Y,Z,U,IO1,IO2,IO3,ISOLV,XL,UL,VTEST,DEN,IZTAL,VOL1,VOL2)
      IF(ISOLV) 529,527,529
  529 CALL MAEG(T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF(K-IQ) 527,530,530
  530 LIN(5)=LIN(5)+1
      IF(IUTSK) 531,533,531
  531 WRITE(NDISP,532) K,X,Y,Z,U
  532 FORMAT(20H MONOCLINIC TEST. K=,I3,6H XYZU=,4F10.6)
  533 CALL LAGRM(K,U)
      GO TO 527
C     TRICLINIC TEST.
 4000 LIN(8)=2
      IF(LIN(16) .EQ. 0) GO TO 4092
      LIN(16)=0
      GO TO 4012
 4092 GO TO(4001,4999),ITRIC
 4001 WRITE(IOUT,4075)
 4075 FORMAT(15H TRICLINIC TEST)
      WRITE(NDISP,4075)
      ISYM=6
      ID=0
      SSQM=T1(IB)+D1+D2
      ISTP=0
      DO 4020 I=1,5
      BB(I,1)=0.
      BBB(I,1)=0.
      DO 4020 J=1,5
      AA(I,J)=0.
      AAA(I,J)=0.
 4020 CONTINUE
 4003 ILINE=0
      CALL LINES(T1,SQ1,SQQ,ILINE,SLD)
      KLINE=1
 4005 DO 4022 I=1,5
      BB(I,1)=BBB(I,1)
      DO 4022 J=1,5
      AA(I,J)=AAA(I,J)
 4022 CONTINUE
      CALL AMB(AA,BB,SQQ,ISTP,SQ1,IDX)
      DO 4024 I=1,5
      BBB(I,1)=BB(I,1)
      DO 4024 J=1,5
      AAA(I,J)=AA(I,J)
 4024 CONTINUE
      IF(ISTP) 4006,210,4006
 4006 DO 4060 I=1,5
      HL(I)=BB(I,1)
 4060 CONTINUE
      CALL PW84(AA,5,BB,ISOLV)
      IF(ISOLV) 4007,4005,4007
 4007 X1=SQ1
      X2=BB(1,1)
      X3=BB(2,1)
      X4=BB(3,1)
      X5=BB(4,1)
      X6=BB(5,1)
 4057 CALL PW83(X1,X2,X3,X5,X4,X6,WAVE,ISYM,A,B,C,AL,BE,GA,CEM,VOL,IOOK,
     1     VVV)
      IF(IOOK) 4014,4012,4014
 4014 IF(DEN-0.00001) 4030,4030,4031
 4031 RVA=VSK/(VVV**2)
      DO 4032 I=1,IZTAL
      IF(RVA-VOL1(I)) 4032,4033,4033
 4033 IF(RVA-VOL2(I)) 4030,4030,4012
 4032 CONTINUE
 4030 CALL COUNT(T1,X1,X2,X3,X4,X5,X6,SSQM,D1,D2,SSQTL,IB,K)
      IF(K-IQ) 4012,4008,4008
 4008 LIN(7)=LIN(7)+1
      IF(IUTSK) 4009,4011,4009
 4009 WRITE(NDISP,4010) K,X1,X2,X3,X4,X5,X6,VVV
 4010 FORMAT(1X,10H TRICL. K=,I3,9H A11-A23=,6F9.6,/,7H VOLYM=,F7.1,
     1 5H A**3)
      AL=AL/RAD
      BE=BE/RAD
      GA=GA/RAD
      IF(AL) 4077,4077,4078
 4077 AL=AL+180.
 4078 IF(BE) 4079,4079,4080
 4079 BE=BE+180.
 4080 IF(GA) 4081,4081,4082
 4081 GA=GA+180.
 4082 WRITE(NDISP,4058) A,B,C,AL,BE,GA
 4058 FORMAT(7H A,B,C=,3F6.2,10H AL,BE,GA=,3F8.2)
 4011 CALL STORE(VVV,X1,X2,X3,X4,X5,X6,ID,K,RO,KOTA)
      IF(K-IB) 4012,4091,4091
 4091 LIN(16)=1
      GO TO 210
 4012 CALL LINES(T1,SQ1,SQQ,KLINE,SLD)
      IF(KLINE) 4003,4003,4056
 4056 CALL MATMUL(AA,SQQ,SLD,HL,X2,X3,X4,X5,X6)
      GO TO 4057
 4999 STOP
  900 RETURN
      END
      SUBROUTINE COUNT(T1,X1,X2,X3,X4,X5,X6,SSQM,DIT20,DIT40,SSQTL,IB,K)
      DIMENSION T1(100),KST(30)
      TIDH=SSQM/X1
      TIDK=SSQM/X2
      TIDL=SSQM/X3
      HHM=0.0
      HKM=0.0
      HLM=0.0
    1 HHM=HHM+1.0
      IF(HHM*HHM-TIDH)1,1,3
    3 HKM=HKM+1.0
      IF(HKM*HKM-TIDK)3,3,5
    5 HLM=HLM+1.0
      IF(HLM*HLM-TIDL)5,5,6
    6 DO 7 I=1,IB
      KST(I)=0
    7 CONTINUE
      ISLUT=0
    8 CALL HKL(HHM,HKM,HLM,X1,X2,X3,X4,X5,X6,SSQT,ISLUT)
      IF(ISLUT)17,17,9
    9 IF(SSQM-SSQT)8,10,10
   10 IF(SSQT-SSQTL)11,12,12
   11 CSQ=DIT20
      GO TO 13
   12 CSQ=DIT40
   13 DO 16 I=1,IB
      SSQO=T1(I)
      STEST=SSQT-SSQO
      IF(ABS(STEST)-CSQ)15,15,14
   14 IF(STEST)8,16,16
   15 KST(I)=1
   16 CONTINUE
      GO TO 8
   17 K=0
      DO 18 I=1,IB
      K=K+KST(I)
   18 CONTINUE
      RETURN
      END
      SUBROUTINE HKL(HHM,HKM,HLM,X1,X2,X3,X4,X5,X6,SSQT,ISLUT)
      SAVE
      IF(ISLUT) 2,1,2
    1 ISLUT=1
      HH=0.0
      HK=0.0
      HL=0.0
      H1=0.0
      H2=0.0
      H3=0.0
      H4=0.0
      H5=0.0
      H6=0.0
    2 IF(ABS(HK)+HL) 3,4,3
    3 IF(HH) 4,4,5
    5 HH=-HH
      H4=-H4
      H5=-H5
      GO TO 30
    4 IF(ABS(HH)-HHM) 6,7,7
    6 HH=ABS(HH)+1.0
      H1=HH*HH*X1
      H4=HH*HK*X4
      H5=HH*HL*X5
      GO TO 30
    7 HH=0.0
      H1=0.0
      H4=0.0
      H5=0.0
      IF(HL) 9,9,8
    8 IF(HK) 9,9,10
   10 HK=-HK
      H4=-H4
      H6=-H6
      GO TO 30
    9 IF(ABS(HK)-HKM) 11,12,12
   11 HK=ABS(HK)+1.0
      H2=HK*HK*X2
      H4=HH*HK*X4
      H6=HK*HL*X6
      GO TO 30
   12 HK=0.0
      H2=0.0
      H4=0.0
      H6=0.0
      IF(HL-HLM) 13,99,99
   13 HL=HL+1.0
      H3=HL*HL*X3
      H5=HH*HL*X5
      H6=HK*HL*X6
   30 SSQT=H1+H2+H3+H4+H5+H6
      RETURN
   99 ISLUT=0
      RETURN
      END
      SUBROUTINE STORE(PV,X1,X2,X3,X4,X5,X6,ID,K,RO,KOTA)
      DIMENSION RO(7,30),KOTA(30)
      IF(X1-X2)4,2,2
    2 IF(X1-X3)3,9,9
    3 XA=X3
      YA=X1
      ZA=X2
      XYA=X5
      XZA=X6
      YZA=X4
      GO TO 12
    4 IF(X1-X3)6,5,5
    5 XA=X2
      YA=X1
      ZA=X3
      XYA=X4
      XZA=X6
      YZA=X5
      GO TO 12
    6 IF(X2-X3)7,8,8
    7 XA=X3
      YA=X2
      ZA=X1
      XYA=X6
      XZA=X5
      YZA=X4
      GO TO 12
    8 XA=X2
      YA=X3
      ZA=X1
      XYA=X6
      XZA=X4
      YZA=X5
      GO TO 12
    9 IF(X2-X3)10,11,11
   10 XA=X1
      YA=X3
      ZA=X2
      XYA=X5
      XZA=X4
      YZA=X6
      GO TO 12
   11 XA=X1
      YA=X2
      ZA=X3
      XYA=X4
      XZA=X5
      YZA=X6
   12 IF(ID)16,13,16
   13 ID=ID+1
   14 I=ID
   15 RO(7,I)=PV
      RO(1,I)=XA
      RO(2,I)=YA
      RO(3,I)=ZA
      RO(4,I)=XYA
      RO(5,I)=XZA
      RO(6,I)=YZA
      KOTA(I)=K
      RETURN
   16 DO 20 I=1,ID
      IF(ABS(RO(1,I)-XA)-0.00005)17,20,20
   17 IF(ABS(RO(2,I)-YA)-0.00005)18,20,20
   18 IF(ABS(RO(3,I)-ZA)-0.00005)31,20,20
   31 IF(ABS(RO(7,I)-PV)-0.05) 19,20,20
   19 RETURN
   20 CONTINUE
      DO 23 I=1,ID
      IF(KOTA(I)-K)22,21,23
   21 IF(RO(7,I)-PV)22,23,23
   22 IF(I-10) 24,14,99
   23 CONTINUE
      IF(ID-10) 13,19,19
   24 N=9
      IF(ID-10) 25,26,26
   25 ID=ID+1
   26 M=N+1
      KOTA(M)=KOTA(N)
      DO 27 J=1,7
      RO(J,M)=RO(J,N)
   27 CONTINUE
      IF(N-I)99,15,28
   28 N=N-1
      GO TO 26
   99 STOP
      END
      SUBROUTINE AMB(A,B,SQ,ISTOP,A11,IDX)
      DIMENSION A(6,6),B(6,1),SQ(5),V(72),IDX(5)
      DATA V/0,1,0,0,0,0,0,0,1,0,0,0,1,1,0,1,0,0,1,1,0,-1,0,0,0,1,1,
     *0,0,1,0,1,1,0,0,-1,1,0,1,0,1,0,1,0,1,0,-1,0,1,1,1,1,1,1,1,1,1,
     *-1,-1,1,1,1,1,-1,1,-1,1,1,1,1,-1,-1/
      IF(ISTOP) 10,1,10
    1 ISTOP=1
      DO 2 I=1,5
      IDX(I)=I
      L=(I-1)*6+1
      B(I,1)=SQ(I)-V(L)*A11
      DO 3 M=1,5
      N=L+M
      A(I,M)=V(N)
    3 CONTINUE
    2 CONTINUE
      GO TO 90
   10 I5=IDX(5)+1
      IF(I5-12) 11,11,20
   11 IDX(5)=I5
      DO 13 I=1,4
      IF(I5-IDX(I)) 13,10,13
   13 CONTINUE
      L=(I5-1)*6+1
      B(5,1)=SQ(5)-V(L)*A11
      DO 14 M=1,5
      N=L+M
      A(5,M)=V(N)
   14 CONTINUE
      GO TO 90
   20 IDX(5)=0
   12 I4=IDX(4)+1
      IF(I4-12) 21,21,30
   21 IDX(4)=I4
      DO 23 I=1,3
      IF(I4-IDX(I)) 23,12,23
   23 CONTINUE
      L=(I4-1)*6+1
      B(4,1)=SQ(4)-V(L)*A11
      DO 24 M=1,5
      N=L+M
      A(4,M)=V(N)
   24 CONTINUE
      GO TO 10
   30 IDX(4)=0
   39 I3=IDX(3)+1
      IF (I3-8) 31,31,40
   31 IDX(3)=I3
      DO 33 I=1,2
      IF(I3-IDX(I)) 33,39,33
   33 CONTINUE
      L=(I3-1)*6+1
      B(3,1)=SQ(3)-V(L)*A11
      DO 34 M=1,5
      N=L+M
      A(3,M)=V(N)
   34 CONTINUE
      GO TO 12
   40 IDX(3)=0
   49 I2=IDX(2)+1
      IF(I2-8) 41,41,50
   41 IDX(2)=I2
      IF(I2-IDX(1)) 43,49,43
   43 L=(I2-1)*6+1
      B(2,1)=SQ(2)-V(L)*A11
      DO 44 M=1,5
      N=L+M
      A(2,M)=V(N)
   44 CONTINUE
      GO TO 30
   50 IDX(2)=0
   59 I1=IDX(1)+2
      IF(I1-5) 51,51,60
   51 IDX(1)=I1
      L=(I1-1)*6+1
      B(1,1)=SQ(1)-V(L)*A11
      DO 54 M=1,5
      N=L+M
      A(1,M)=V(N)
   54 CONTINUE
      GO TO 40
   60 ISTOP=0
   90 RETURN
      END
      SUBROUTINE LINES(T1,SQ1,SQQ,ILINE,SQOLD)
      DIMENSION T1(100),SQQ(5),SQOLD(5)
      IF(ILINE) 3,1,3
    1 ILINE=1
      SQ1=T1(1)
      DO 2 I=1,5
      N=I+1
      SQQ(I)=T1(N)
    2 CONTINUE
      RETURN
    3 DO 30 I=1,5
      SQOLD(I)=SQQ(I)
   30 CONTINUE
      GO TO(11,12,13,14,15,16,17,18,19),ILINE
   11 SQQ(5)=T1(7)
   20 ILINE=ILINE+1
      RETURN
   12 SQQ(5)=T1(8)
      GO TO 20
   13 SQQ(4)=T1(6)
      SQQ(5)=T1(7)
      GO TO 20
   14 SQQ(5)=T1(8)
      GO TO 20
   15 SQQ(3)=T1(5)
      GO TO 20
   16 SQQ(5)=T1(7)
      GO TO 20
   17 SQQ(3)=T1(4)
      SQQ(4)=T1(5)
      SQQ(5)=T1(9)
      GO TO 20
   18 SQQ(3)=T1(5)
      SQQ(4)=T1(7)
      SQQ(5)=T1(8)
      GO TO 20
   19 ILINE=0
      RETURN
      END
      SUBROUTINE MATMUL(AA,SQQ,SLD,HL,X2,X3,X4,X5,X6)
      DIMENSION AA(6,6),SQQ(5),SLD(5),HL(5),X(5)
      DO 10 I=1,5
      HL(I)=HL(I)-SLD(I)+SQQ(I)
   10 CONTINUE
      DO 20 I=1,5
      X(I)=0.0
      DO 30 K=1,5
      X(I)=X(I)+AA(I,K)*HL(K)
   30 CONTINUE
   20 CONTINUE
      X2=X(1)
      X3=X(2)
      X4=X(3)
      X5=X(4)
      X6=X(5)
      RETURN
      END
      SUBROUTINE TWODIM
      COMMON T1(100),TAGUT(100),N1,D1,SSQTL,D2,WAVE,IHM,IKM,ILM,ISYM,INR
     >      ,ISUB,A11,A22,A33,A13,A12,A23,LIN(20),REV(20),IGOT,IS1
     >      ,IS2,IS3,IS4,IA,IB,IQ,IH1,IK1,IL1,IH2,IK2,IL2,IH3,IK3,IL3
     >      ,SQ1,SQ2,SQ3,X,Y,Z,IC(20),ISQM,KOTA(30),RO(7,30),ID,IIN,IOUT
      COMMON NH1,NK1,NL1,NS1,NH2,NK2,NL2,NS2,NH3,NK3,NL3,NS3,IUTSK,IMER,
     >       INIX,IOMIT,ISC,KLART,VOL,VTEST,XL,YL,IT20,IT40,CEM,KHA,KKA,
     >       KLA,KSA,KHB,KKB,KLB,KSB,MONO,MH1,MK1,ML1,MH2,MK2,ML2,MH3,
     >       MK3,ML3,MS1,MS2,MS3,MH4,MK4,ML4,MS4,UL,LGO,MG,MS,MRA,MRB,
     >       NU,DEN,EDEN,RMOLV
      COMMON VOL1(100),VOL2(100)
      DIMENSION A(25,4),B(5,4)
      IZTAL=LIN(14)
      NDISP=LIN(11)
      ITWO=0
      WRITE (NDISP,*)  ' Entrance of SHORT AXIS test  '
      ITLIM=100000.0*SSQTL+0.5
      IANTAL=0
C SET THE DIFFERENCE PARAMETERS FOR THE INNER BLOCK.
      IF (T1(4)-SSQTL) 1,1,2
    1 DIFF4=D1
      GOTO 3
    2 DIFF4=D2
    3 IF (T1(5)-SSQTL) 4,4,5
    4 DIFF5=D1
      GOTO 6
    5 DIFF5=D2
    6 IF (T1(6)-SSQTL) 7,7,8
    7 DIFF6=D1
      GOTO 9
    8 DIFF6=D2
    9 CONTINUE
C CLEAR THE INDEXES OF THE FIRST LINE
      IH1=0
      IL1=0
C GENERATE INDEXES OF THE FIRST LINE
   10 CALL GMHL (2,2,4,2,ISLUT1,IH1,IL1)
      IF (ISLUT1) 11,1000,11
   11 IF(IABS(IH1) .LT. IL1) GO TO 10
      RH1=IH1
      RL1=IL1
      RHH1=RH1*RH1
      RLL1=RL1*RL1
      RHL1=RH1*RL1
C CLEAR THE INDEXES OF THE SECOND LINE
      IH2=0
      IL2=0
C GENERATE THE INDEXES OF THE SECOND LINE
   20 CALL GMHL (2,2,4,2,ISLUT2,IH2,IL2)
      IF (ISLUT2)  21,10,21
   21 IF (IH1-IH2) 23,22,23
   22 IF (IL1-IL2) 23,20,23
   23 RH2=IH2
      RL2=IL2
      RHH2=RH2*RH2
      RLL2=RL2*RL2
      RHL2=RH2*RL2
C CLEAR THE INDEXES OF THE THIRD LINE
      IH3=0
      IL3=0
C GENERATE THE INDEXES OF THE THIRD LINE
   30 CALL GMHL (2,2,4,2,ISLUT3,IH3,IL3)
      IF (ISLUT3)  31,20,31
   31 IF (IH1-IH3) 33,32,33
   32 IF (IL1-IL3) 33,30,33
   33 IF (IH2-IH3) 35,34,35
   34 IF (IL2-IL3) 35,30,35
   35 RH3=IH3
      RL3=IL3
      RHH3=RH3*RH3
      RLL3=RL3*RL3
      RHL3=RH3*RL3
C COMPUTE THE FIRST DETERMINANT USED TO OBTAIN A11,A33,A13
      IDET0=IDET(IH1*IH1,IL1*IL1,IH1*IL1,IH2*IH2,IL2*IL2,IH2*IL2,
     >           IH3*IH3,IL3*IL3,IH3*IL3)
C CHECK IF IDET0 EQUALS 0,AND JUMP TO 30 IF TRUE.
      IF (IDET0) 40,30,40
C COMPUTE THE NEXT 3 DETERMINANTS USED TO OBT. A11,A33 AND A13
   40 DET0=FLOAT (IDET0)
      DET1=DET(T1(1),RLL1,RHL1,T1(2),RLL2,RHL2,T1(3),RLL3,RHL3)
C SOLVE THE THREE PARAMETERS P1,P2 AND P3
      P1=DET1/DET0
      IF (P1) 30,30,41
   41 DET2=DET(RHH1,T1(1),RHL1,RHH2,T1(2),RHL2,RHH3,T1(3),RHL3)
      P2=DET2/DET0
      IF (P2) 30,30,42
   42 DET3=DET(RHH1,RLL1,T1(1),RHH2,RLL2,T1(2),RHH3,RLL3,T1(3))
      P3=DET3/DET0
      IF(P3) 710,43,43
  710 IF(ABS(P3)-D2) 711,711,30
  711 P3=0.
   43 CONTINUE
      IF(P3-D2) 712,712,713
  712 P3=0.
  713 CONTINUE
C-----------------------------------------------------------------------
C CHECK IF Q(4) IS INDEXABLE.
C-----------------------------------------------------------------------
      IH4=0
      IL4=0
   50 CALL GMHL (4,4,6,2,ISLUT4,IH4,IL4)
      IF (ISLUT4) 51,30,51
C COMPUTE Q(CALC)-Q(OBS)
   51 TEST=IH4*IH4*P1+IL4*IL4*P2+IH4*IL4*P3 - T1(4)
C IF TEST NEARLY EQUALS ZERO THEN JUMP TO 52,ELSE 50
      IF (ABS(TEST)-DIFF4) 52,52,50
C-----------------------------------------------------------------------
C CHECK IF Q(5) IS INDEXABLE.
C-----------------------------------------------------------------------
   52 IH5=0
      IL5=0
   53 CALL GMHL (4,4,6,2,ISLUT5,IH5,IL5)
      IF (ISLUT5) 54,30,54
C COMPUTE Q(OBS)-Q(CALC)
   54 TEST=IH5*IH5*P1+IL5*IL5*P2+IH5*IL5*P3 - T1(5)
C IF TEST NEARLY EQUALS ZERO THEN JUMP TO 55,ELSE 53
      IF (ABS(TEST)-DIFF5) 55,55,53
C-----------------------------------------------------------------------
C CHECK IF Q(6) IS INDEXABLE
C-----------------------------------------------------------------------
   55 IH6=0
      IL6=0
   56 CALL GMHL (5,5,6,2,ISLUT6,IH6,IL6)
      IF (ISLUT6) 57,30,57
C COMPUTE Q(OBS)-Q(CALC)
   57 TEST=IH6*IH6*P1+IL6*IL6*P2+IH6*IL6*P3 - T1(6)
C IF TEST NEARLY EQUALS ZERO THEN JUMP TO 58,ELSE 56
      IF (ABS(TEST)-DIFF6) 58,58,56
C SAVE A SOLUTION IN THE hk0 AND 0kl CASES.
   58 CALL SAVE (P1,P2,P3,A,IANTAL,INOST)
      IF(INOST .EQ. 0) GO TO 30
C CHECK THE PARAMETERS OF THE h0l-CASE ( beta and cell-lengths )
      HA=P1*P2
      U=P3*P3/HA
C TEST THE MONOCLINIC ANGLE.
      IF (U.GT.UL) THEN
         GOTO 30
      ENDIF
C TEST A* AND C*.
      HE=1.0-0.25*U
      IF (P1*HE.LE.XL) THEN
         GOTO 30
      ENDIF
      IF (P2*HE.LE.XL) THEN
         GOTO 30
      ENDIF
      BETA =  (ASIN (SQRT(HE)))/0.0174533
      ITWO = ITWO + 1
      GOTO 30
 1000 CONTINUE
      DO 1010 I=1,IANTAL
      INDEX1=ISERCH(T1,A(I,1),A(I,2),A(I,3),7,IB,D1,D2,SSQTL)
C CHECK THAT INDEX1 IS GREATER THAN ZERO.
      IF (INDEX1) 400,400,201
  201 INDEX2=ISERCH(T1,A(I,1),A(I,2),A(I,3),INDEX1+1,IB,D1,D2,SSQTL)
C CHECK THAT INDEX2 IS GREATER THAN ZERO.
      IF (INDEX2) 800,800,202
  800 IF (A(I,3)-D1) 400,801,801
  202 INDEX3=ISERCH(T1,A(I,1),A(I,2),A(I,3),INDEX2+1,IB,D1,D2,SSQTL)
C CHECK THAT INDEX3 IS GREATER THAN ZERO.
      IF (INDEX3) 400,400,302
  400 WRITE (NDISP,401) ' INDEX RETURNED AS ZERO.'
  401 FORMAT (A)
      GOTO 1010
C-----------------------------------------------------------------------
C               Test the hk0-case.
C-----------------------------------------------------------------------
  302 IF ( A(I,3) .GT. D1 ) THEN
         GOTO 801
      ENDIF
      A11=A(I,1)
      A22=A(I,2)
      X=A11
      Y=A22
      IH1=0
      IK1=0
      IL1=1
      GO TO 304
  303 CALL GMHL (4,4,6,2,ISLUT1,IH1,IK1)
      IF (ISLUT1) 304,1002,304
  304 IH2=0
      IK2=0
      IL2=1
      GO TO 306
  305 CALL GMHL (4,4,6,2,ISLUT2,IH2,IK2)
      IF (ISLUT2) 303,303,306
C-----------------------------------------------------------------------
C (hk0):
C
C T1(i)-IH1*IH1*A11-IK1*IK1*A22 = IL1*IL1*A33 + IH1*IL1*A13
C T1(k)-IH2*IH2*A11-IK2*IK2*A22 = IL2*IL2*A33 + IH2*IL2*A13
C
C A33 AND A13 ARE SOLVED FROM THIS EQ.SYSTEM.
C-----------------------------------------------------------------------
  306 B1=T1(INDEX1)-IH1*IH1*A11-IK1*IK1*A22
      B2=T1(INDEX2)-IH2*IH2*A11-IK2*IK2*A22
      IXDET=IH2-IH1
      XDET=IL1*IL1*IH2*IL2-IL2*IL2*IH1*IL1
      IF (IXDET) 305,305,307
  307 A33=(B1*IH2*IL2-B2*IH1*IL1)/XDET
      A13=(IL1*IL1*B2-IL2*IL2*B1)/XDET
      IOK=IOKAY(A11,A22,A33,A13,UL,VTEST,CEM,WAVE,XL,VOL1,VOL2,DEN,
     * IZTAL)
      IF (IOK) 321,321,331
  331 Z=A33
      U=A13
      CALL MAEG (T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF (K-IQ) 311,308,308
  308 LIN(6)=LIN(6)+1
      IF (IUTSK) 309,310,309
  309 WRITE (NDISP,320) K,X,Y,Z,U
  310 CALL LAGRM (K,U)
  311 CONTINUE
  321 B1=T1(INDEX1)-IH1*IH1*A11-IK1*IK1*A22
      B2=T1(INDEX3)-IH2*IH2*A11-IK2*IK2*A22
      A33=(B1*IH2*IL2-B2*IH1*IL1)/XDET
      A13=(IL1*IL1*B2-IL2*IL2*B1)/XDET
      IOK = IOKAY (A11,A22,A33,A13,UL,VTEST,CEM,WAVE,XL,
     * VOL1,VOL2,DEN,IZTAL)
      IF (IOK) 322,322,333
  333 Z=A33
      U=A13
      CALL MAEG (T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF (K-IQ) 315,312,312
  312 LIN(6)=LIN(6)+1
      IF (IUTSK) 313,314,313
  313 WRITE (NDISP,320) K,X,Y,Z,U
  314 CALL LAGRM (K,U)
  315 CONTINUE
  322 B1=T1(INDEX2)-IH1*IH1*A11-IK1*IK1*A22
      B2=T1(INDEX3)-IH2*IH2*A11-IK2*IK2*A22
      A33=(B1*IH2*IL2-B2*IH1*IL1)/XDET
      A13=(IL1*IL1*B2-IL2*IL2*B1)/XDET
      IOK = IOKAY (A11,A22,A33,A13,UL,VTEST,CEM,WAVE,XL,
     * VOL1,VOL2,DEN,IZTAL)
      IF (IOK) 305,305,335
  335 Z=A33
      U=A13
      CALL MAEG (T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF (K-IQ) 319,316,316
  316 IF (IUTSK) 317,318,317
  317 LIN(6)=LIN(6)+1
      WRITE (NDISP,320) K,X,Y,Z,U
  318 CALL LAGRM (K,U)
  319 CONTINUE
  320 FORMAT (20H SHORT AXIS TEST. K=,I3,6H XYZU=,4F10.6)
      GOTO 305
C-----------------------------------------------------------------------
C               Test the 0kl-case.
C-----------------------------------------------------------------------
 1002 A22=A(I,1)
      A33=A(I,2)
      Y=A22
      Z=A33
      IH1=0
      IK1=0
      IL1=0
  500 CALL GEHKL (1,4,4,6,2,ISLUT1,IH1,IK1,IL1)
      IF (ISLUT1) 501,1010,501
  501 IH2=0
      IK2=0
      IL2=0
  502 CALL GEHKL (1,4,4,6,2,ISLUT2,IH2,IK2,IL2)
      IF (ISLUT2) 500,500,503
C-----------------------------------------------------------------------
C (0kl):
C
C T1(i)-IK1*IK1*A22-IL1*IL1*A33 = IH1*IH1*A11 + IH1*IL1*A13
C T1(k)-IK2*IK2*A22-IL2*IL2*A33 = IH2*IH2*A11 + IH2*IL2*A13
C
C A11 AND A13 ARE SOLVED FROM THIS EQ.SYSTEM.
C
C-----------------------------------------------------------------------
  503 B1=T1(INDEX1)-IK1*IK1*A22-IL1*IL1*A33
      B2=T1(INDEX2)-IK2*IK2*A22-IL2*IL2*A33
      XDET=IH1*IH1*IH2*IL2-IH1*IL1*IH2*IH2
      IF (ABS(XDET)-0.00001) 502,502,504
  504 A11=(B1*IH2*IL2-B2*IH1*IL1)/XDET
      A13=(IH1*IH1*B2-IH2*IH2*B1)/XDET
      IOK = IOKAY (A11,A22,A33,A13,UL,VTEST,CEM,WAVE,XL,
     * VOL1,VOL2,DEN,IZTAL)
      IF (IOK) 523,523,506
  506 X=A11
      U=A13
      CALL MAEG (T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF (K-IQ) 510,507,507
  507 IF (IUTSK) 508,509,508
  508 LIN(6)=LIN(6)+1
      WRITE (NDISP,320) K,X,Y,Z,U
  509 CALL LAGRM (K,U)
  510 CONTINUE
  523 B1=T1(INDEX1)-IK1*IK1*A22-IL1*IL1*A33
      B2=T1(INDEX3)-IK2*IK2*A22-IL2*IL2*A33
      A11=(B1*IH2*IL2-B2*IH1*IL1)/XDET
      A13=(IH1*IH1*B2-IH2*IH2*B1)/XDET
      IOK = IOKAY (A11,A22,A33,A13,UL,VTEST,CEM,WAVE,XL,
     * VOL1,VOL2,DEN,IZTAL)
      IF (IOK) 524,524,512
  512 X=A11
      U=A13
      CALL MAEG (T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF (K-IQ) 516,513,513
  513 IF (IUTSK) 514,515,514
  514 LIN(6)=LIN(6)+1
      WRITE (NDISP,320) K,X,Y,Z,U
  515 CALL LAGRM (K,U)
  516 CONTINUE
  524 B1=T1(INDEX2)-IK1*IK1*A22-IL1*IL1*A33
      B2=T1(INDEX3)-IK2*IK2*A22-IL2*IL2*A33
      A11=(B1*IH2*IL2-B2*IH1*IL1)/XDET
      A13=(IH1*IH1*B2-IH2*IH2*B1)/XDET
      IOK = IOKAY (A11,A22,A33,A13,UL,VTEST,CEM,WAVE,XL,
     * VOL1,VOL2,DEN,IZTAL)
      IF (IOK) 502,502,518
  518 X=A11
      U=A13
      CALL MAEG (T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF (K-IQ) 522,519,519
  519 IF (IUTSK) 520,521,520
  520 LIN(6)=LIN(6)+1
      WRITE (NDISP,320) K,X,Y,Z,U
  521 CALL LAGRM (K,U)
  522 CONTINUE
      GOTO 502
C     TETS THE (H0L)-CASE
  801 X=A(I,1)
      Z=A(I,2)
      U=A(I,3)
      IH1=0
      IK1=1
      IL1=0
      GO TO 807
  802 CALL GMHL(4,4,6,1,ISLUT1,IH1,IL1)
      IF(ISLUT1 .EQ. 0) GO TO 1010
  807 Y=T1(INDEX1)-IH1*IH1*X-IL1*IL1*Z-IH1*IL1*U
      IF(Y-X) 802,803,803
  803 IOK=IOKAY(X,Y,Z,U,UL,VTEST,CEM,WAVE,XL,VOL1,VOL2,DEN,IZTAL)
      IF(IOK) 802,802,804
  804 CALL MAEG(T1,IC,X,Y,Z,U,IB,ISQM,K,IT20,IT40,ITLIM)
      IF(K-IQ) 802,805,805
  805 IF(IUTSK) 806,808,806
  806 LIN(6)=LIN(6)+1
      WRITE(NDISP,320) K,X,Y,Z,U
  808 CALL LAGRM(K,U)
      GO TO 802
 1010 CONTINUE
      WRITE (NDISP,*) ' Exit of SHORT AXIS test  '
      RETURN
      END
      INTEGER FUNCTION IOKAY (A11,A22,A33,A13,UL,VTEST,CEM,WAVE,XL,
     * VOL1,VOL2,DEN,IZTAL)
C***********************************************************************
C
C THIS FUNCTION TESTS THE FOLLOWING QUANTITETES AND RETURNS
C IOKAY=1 IF TRUE ELSE IOKAY=0.
C
C PI/2 < BETA < 3*PI/4  ! monoclinic angle.
C A,B AND C =< CEM
C VTEST =< VOL
C
C***********************************************************************
      DIMENSION VOL1(100),VOL2(100)
      REAL K
      IF (A11.LT.0.0.OR.A22.LT.0.0.OR.A33.LT.0.0.OR.A13.LE.0.0) THEN
         IOKAY = 0
         RETURN
      ENDIF
      HA=A11*A33
      U=A13*A13/HA
C TEST THE MONOCLINIC ANGLE.
      IF (U.GT.UL) THEN
         IOKAY = 0
         RETURN
      ENDIF
C TEST THE VOLUME.
      HE=1.0-0.25*U
      V=A22*HA*HE
      IF (V.LT.VTEST) THEN
         IOKAY = 0
         RETURN
      ENDIF
      IF(DEN-0.00001) 10,10,9
    9 DO 8 I=1,IZTAL
      IF(V-VOL1(I)) 8,7,7
    7 IF(V-VOL2(I)) 10,10,12
    8 CONTINUE
   12 IOKAY=0
      RETURN
   10 CONTINUE
C TEST A* AND C*.
      IF (A11*HE.LE.XL) THEN
         IOKAY = 0
         RETURN
      ENDIF
      IF (A33*HE.LE.XL) THEN
         IOKAY = 0
         RETURN
      ENDIF
      K=WAVE*WAVE/4.0
      IF (A22/K.LT.(1.0/(CEM*CEM))) THEN
         IOKAY = 0
         RETURN
      ENDIF
      IOKAY = 1
      RETURN
      END
      FUNCTION DET (D1,D2,D3,D4,D5,D6,D7,D8,D9)
      DET1=D1*(D5*D9-D6*D8)
      DET2=D2*(D6*D7-D4*D9)
      DET3=D3*(D4*D8-D5*D7)
      DET=DET1+DET2+DET3
      RETURN
      END
      FUNCTION IDET (I1,I2,I3,I4,I5,I6,I7,I8,I9)
      IDET1=I1*(I5*I9-I6*I8)
      IDET2=I2*(I6*I7-I4*I9)
      IDET3=I3*(I4*I8-I5*I7)
      IDET=IDET1+IDET2+IDET3
      RETURN
      END
      FUNCTION ISERCH (T,P1,P2,P3,IFIRST,ILAST,D1,D2,SSQTL)
C***********************************************************************
C
C THIS FUNCTION RETURNS THE INDEX OF THE FIRST NON-INDEXABLE LINE
C IN THE VECTOR T1,GIVEN AS A PARAMETER.
C SEARCH BEGINS AT INDEX = "FIRST".
C
C***********************************************************************
      DIMENSION T(100)
      DO 10 I=IFIRST,ILAST
         IH=0
         IL=0
   20    CALL GMHL(4,4,6,2,ISLUT,IH,IL)
         IF (ISLUT) 27,26,27
   27    TEST=IH*IH*P1+IL*IL*P2+IH*IL*P3
         DIFF=ABS(T(I)-TEST)
         IF (TEST-SSQTL) 31,31,32
   31    D12=D1
         GOTO 33
   32    D12=D2
   33    IF (DIFF-D12) 10,10,20
C A NOT INDEXABLE LINE FOUND!!!
   26    ISERCH=I
         RETURN
   10 CONTINUE
      ISERCH=0
      RETURN
      END
      SUBROUTINE  GMHL (MAXH,MAXL,MAXS,MORMO,ISLUT,IH,IL)
C MORMO = 1  ;POSITIVE H
C MORMO = 2  ;POSITIVE AND NEGATIVE H
C ISLUT = 0  ;END OF GENERATION
C ISLUT = 1  ;CONINUE GENERATION
      ISLUT=1
      GOTO (1,2),MORMO
    1 IF (IH-MAXH) 3,4,4
    3 IH=IH+1
      ISUM=IH+IL
      IF (MAXS-ISUM) 4,5,5
    5 RETURN
    4 IH=0
      IF (IL-MAXL) 6,7,7
    6 IL=IL+1
      ISUM=IH+IL
      IF (MAXS-ISUM) 7,5,5
    7 IL=0
      ISLUT=0
      RETURN
    2 IF (IH*IL) 8,1,9
    8 IH=-IH
      GOTO 1
    9 IH=-IH
      RETURN
      END
      SUBROUTINE SAVE (A1,A2,A3,PAR,N,INOST)
C***********************************************************************
C THIS SUBROUTINE SAVE A SOLUTION IN A 5*3 MATRIX OF A1,A2 AND A3.
C THE PARAMETERS ARE CALLED: (hk0);     A1=A11,A2=A22 AND A3=A12. "<0>"
C                            (0kl);     A1=A22,A2=A33 AND A3=A23. "<0>"
C N IS THE NUMBER OF SOLUTIONS, IF N=0 THEN THE SET ALL THE ELEMENTS
C IN THE MATRIX "PAR" TO ZERO.
C***********************************************************************
      DIMENSION PAR(25,4)
C CHECK IF N=0,AND ZERO THE "PAR"-MATRIX IF TRUE.
      INOST=1
      IF(A2 .LT. A1) GO TO 50
      HLP = A1
      A1 = A2
      A2 = HLP
   50 IF (N) 20,10,20
   10 DO 11 I=1,25
         PAR(I,1) = 0.0
         PAR(I,2) = 0.0
         PAR(I,3) = 1.0
         PAR(I,4) = 0.0
   11 CONTINUE
   20 DO 30 I=1,25
C
C CHECK IF A3 IS LESS THAN PAR(I,3),THEN MOVE DOWN THE
C REST OF THE MATRIX AND SAVE THE SOLUTION.
C
      IF (A3-PAR(I,3)) 31,37,30
   37 IF (A1-PAR(I,1)) 30,40,31
   40 IF (A2 .EQ. PAR(I,2)) THEN
         INOST=0
         RETURN
      ELSEIF (A2 .LT. PAR(I,2)) THEN
         GO TO 30
         ELSE
         GO TO 31
      ENDIF
   31 DO 32 K=24,I,-1
         PAR (K+1,1) = PAR (K,1)
         PAR (K+1,2) = PAR (K,2)
         PAR (K+1,3) = PAR (K,3)
         PAR (K+1,4) = PAR (K,4)
   32 CONTINUE
      PAR (I,1) = A1
      PAR (I,2) = A2
      PAR (I,3) = A3
      PAR (I,4) = A1*A2
      GOTO 35
   30 CONTINUE
   35 N=N+1
      IF (N-25) 38,38,36
   36 N=25
   38 RETURN
      END
c      SUBROUTINE GET_USED_CPU_TIME (T)
cC
cC SETS T TO THE CPU TIME USED SO FAR IN SECONDS (NOT NECESSARILY
cC ZERO AT THE START).
c        INTEGER*2   JPT2(8)
c        INTEGER*4   JPT4(4), SYS$GETJPI, STATUS, TICKS
c        INCLUDE     '($JPIDEF)/NOLIST'
c        EXTERNAL    SYS$GETJPI
c        EQUIVALENCE (JPT2,JPT4)
c        JPT2(1) = 4
c        JPT2(2) = JPI$_CPUTIM
c        JPT4(2) = %LOC(TICKS)
c        JPT4(3) = 0
c        JPT4(4) = 0
c        STATUS = SYS$GETJPI (,,,JPT4,,,)
c        IF (.NOT.STATUS) CALL LIB$STOP (%VAL(STATUS))
c        T = 0.01 * TICKS
c        RETURN
c        END
      subroutine get_used_cpu_time (t)
      real tusrsys(2)

      call etime(tusrsys)
      t = tusrsys(1)

      return
      end
