      FUNCTION AIJ(R,I,J) 
      DIMENSION R(3)
      IF (I-2) 300,310,320
  300 IF (J-2) 330,340,350
  310 IF (J-2) 350,330,340
  320 IF (J-2) 340,350,330
  330 AIJ=0.
      RETURN
  340 MM=MOD(I+1,3)+1
      AIJ=R(MM) 
      RETURN
  350 MM=MOD(I,3)+1
      AIJ=-R(MM)
      RETURN
      END 
      FUNCTION BIJK(R,I,JJ,KK)
      DIMENSION R(3)
      DIV=2.
      IF(JJ.EQ.KK)DIV=1.
      J=MIN0(JJ,KK) 
      K=MAX0(JJ,KK) 
      IF (I-2) 300,330,360
  300 IF (J-2) 310,320,420
  310 IF (K-2) 390,400,410
  320 IF (K-2) 420,420,390
  330 IF (J-2) 340,350,420
  340 IF (K-2) 420,410,390
  350 IF (K-2) 390,390,400
  360 IF (J-2) 370,380,390
  370 IF (K-2) 420,390,400
  380 IF(K-2)420,420,410
  390 BIJK=0. 
      RETURN
  400 MM=MOD(I,3)+1
      BIJK=.5*R(MM)/DIV 
      RETURN
  410 MM=MOD(I+1,3)+1
      BIJK=.5*R(MM)/DIV 
      RETURN
  420 BIJK=-.5*R(I)/DIV 
      RETURN
      END 
      FUNCTION CIJKL(R,I,JJ,KK,LL)
      DIMENSION R(3)
      IF(JJ.EQ.KK.AND.JJ.EQ.LL)GO TO 340
      IF(JJ.NE.KK.AND.JJ.NE.LL.AND.KK.NE.LL)GO TO 350
      IF(JJ-KK)300,310,320
  300 J=JJ
      L=KK
      K=LL
      GO TO 330
  310 J=MIN0(JJ,LL) 
      K=KK
      L=MAX0(JJ,LL) 
      GO TO 330
  320 J=KK
      K=LL
      L=JJ
  330 DIV=3.
      GO TO 360
  340 J=JJ
      K=JJ
      L=JJ
      DIV=1.
      GO TO 360
  350 J=1
      K=2
      L=3
      DIV=6.
  360 IF (I-2) 370,420,480
  370 IF (J-2) 380,400,560
  380 IF (K-2) 390,540,540
  390 IF (L-2) 540,550,560
  400 IF (K-2) 410,410,550
  410 IF (L-2) 550,550,560
  420 IF (J-2) 430,460,550
  430 IF (K-2) 440,450,560
  440 IF (L-2) 560,540,550
  450 IF (L-2) 560,560,540
  460 IF (K-2) 470,470,540
  470 IF (L-2) 540,540,550
  480 IF (J-2) 490,520,540
  490 IF (K-2) 500,510,550
  500 IF (L-2) 550,560,540
  510 IF (L-2) 550,550,540
  520 IF (K-2) 530,530,560
  530 IF (L-2) 560,560,540
  540 CIJKL=0.
      RETURN
  550 MM=MOD(I+1,3)+1
      CIJKL=-R(MM)/(6.*DIV) 
      RETURN
  560 MM=MOD(I,3)+1
      CIJKL=R(MM)/(6.*DIV)
      RETURN
      END 
      SUBROUTINE CORECT(R3,R1,R2,ELT,THETA,RC,BIGR,BIGRP,RP2) 
      DIMENSION R3(3),R1(3),R2(3),RC(3),BIGR(3),BIGRP(3),ELT(3,3),
     1RP(3),RPP(3),ACON(3,3),BCON(3)
      RP2=0.
      DO 300 I=1,3
      RP(I)=R2(I)-R1(I) 
  300 RP2=RP2+RP(I)**2
      ISIG=1
      DO 310 I=1,3
      BCON(I)=RC(I) 
      DO 310 J=1,3
      ACON(I,J)=FIJ(THETA,ELT,RP,RP2,I,J) 
  310 BCON(I)=BCON(I)+ACON(I,J)*R2(J) 
      GO TO 330
      ENTRY CORLIB(R3,ELT,RC) 
      DO 320 I=1,3
      BCON(I)=RC(I) 
      DO 320 J=1,3
  320 ACON(I,J)=0.
      ISIG=0
  330 DO 340 I=1,3
      DO 340 J=1,3
  340 ACON(I,J)=ACON(I,J)+EIJ(ELT,I,J)
      CALL SOLVE(ACON,BCON,R3)
      IF(ISIG)350,350,360
  350 RETURN
  360 DO 370 I=1,3
  370 RPP(I)=R3(I)-R2(I)
      CALL VXV(RP,RPP,BIGR) 
      CALL VXV(RP,BIGR,BIGRP) 
      RETURN
      END 
      FUNCTION DETERM(C1,C2,C3) 
      DIMENSION C1(3),C2(3),C3(3) 
      DETERM=C1(1)*C2(2)*C3(3)+C1(2)*C2(3)*C3(1)+C1(3)*C2(1)*C3(2)
     1 -C1(3)*C2(2)*C3(1)-C1(2)*C2(1)*C3(3)-C1(1)*C2(3)*C3(2) 
      RETURN
      END 
      FUNCTION DIJKLM(R,I,JI,KI,LI,MI)
      DIMENSION R(3)
      IN=JI*KI*LI*MI
      IF(MOD(IN,6).EQ.0.AND.IN.LE.18)GO TO 320
      IF(IN.EQ.4.OR.IN.EQ.9.OR.IN.EQ.36)GO TO 310
      IF(JI.NE.KI.OR.KI.NE.LI.OR.LI.NE.MI)GO TO 300
      DIV=1.
      GO TO 330
  300 DIV=4.
      GO TO 330
  310 DIV=6.
      GO TO 330
  320 DIV=12. 
  330 IF(MOD(IN,3).NE.0)GO TO 380
      M=3
      IN=IN/3
      IF(MOD(IN,3).NE.0)GO TO 360
      L=3
      IN=IN/3
      IF(MOD(IN,3).NE.0)GO TO 340
      K=3
      J=IN/3
      GO TO 400
  340 IF(MOD(IN,2).NE.0)GO TO 350
      K=2
      J=IN/2
      GO TO 400
  350 K=1
      J=1
      GO TO 400
  360 IF(MOD(IN,2).NE.0)GO TO 370
      L=2
      IN=IN/2
      IF(MOD(IN,2).NE.0)GO TO 350
      K=2
      J=IN/2
      GO TO 400
  370 L=1
      GO TO 350
  380 IF(MOD(IN,2).NE.0)GO TO 390
      M=2
      IN=IN/2
      GO TO 360
  390 M=1
      GO TO 370
  400 IF (I-2) 410,510,610
  410 IF (J-2) 420,480,720
  420 IF (K-2) 430,460,740
  430 IF (L-2) 440,450,720
  440 IF (M-2) 710,760,740
  450 IF (M-2) 710,720,710
  460 IF (L-2) 710,470,760
  470 IF (M-2) 710,760,740
  480 IF (K-2) 710,490,710
  490 IF (L-2) 710,500,770
  500 IF (M-2) 710,720,710
  510 IF (J-2) 520,580,720
  520 IF (K-2) 530,560,710
  530 IF (L-2) 540,550,770
  540 IF (M-2) 720,740,710
  550 IF (M-2) 710,720,760
  560 IF (L-2) 710,570,740
  570 IF (M-2) 710,740,710
  580 IF (K-2) 710,590,760
  590 IF (L-2) 710,600,720
  600 IF (M-2) 710,710,760
  610 IF (J-2) 620,680,710
  620 IF (K-2) 630,660,760
  630 IF (L-2) 640,650,720
  640 IF (M-2) 720,710,760
  650 IF (M-2) 710,770,740
  660 IF (L-2) 710,670,710
  670 IF (M-2) 710,710,760
  680 IF (K-2) 710,690,740
  690 IF (L-2) 710,700,720
  700 IF (M-2) 710,720,740
  710 DIJKLM=0. 
      RETURN
  720 FAC=1.
  730 DIJKLM=FAC*R(I)/(24.*DIV) 
      RETURN
  740 MM=MOD(I+1,3)+1
  750 DIJKLM=-R(MM)/(24.*DIV) 
      RETURN
  760 MM=MOD(I,3)+1
      GO TO 750
  770 FAC=2.
      GO TO 730
      END 
      FUNCTION D2KDL(R,EL,S,I,J,K,L)
      DIMENSION R(3),EL(3,3),S(3,3),RR(3)
      C(I,J,K,L)=CIJKL(RR,I,J,K,L) 
      C6(I,K,L,M)=6.*C(I,K,L,M) 
      ACC(I,J,K,L,M,N)=3.*AIJ(RR,I,K)*C(J,L,M,N) 
      DO 300 M=1,3
  300 RR(M)=R(M)
      SUM=AIJ(R,I,K)*AIJ(R,J,L)+AIJ(R,I,L)*AIJ(R,J,K) 
      DO 310 M=1,3
      SUM=SUM+S(M,I)*C6(J,K,L,M)+S(M,J)*C6(I,K,L,M) 
      DO 310 N=1,3
      SUM=SUM+EL(M,N)*(ACC(I,J,K,L,M,N)+ACC(I,J,L,K,M,N)+ACC(J,I,K,L,M,N
     1      )+ACC(J,I,L,K,M,N)+ACC(I,J,M,N,K,L)+ACC(J,I,M,N,K,L)+ACC(I,J
     2      ,M,N,L,K)+ACC(J,I,M,N,L,K)+4.*(BIJK(R,I,K,M)*BIJK(R,J,L,N)+ 
     3      BIJK(R,I,L,M)*BIJK(R,J,K,N))) 
  310 CONTINUE
      D2KDL=SUM 
      IF(K.EQ.L)D2KDL=SUM/2.
      RETURN
      END 
      FUNCTION D2KDS(R,EL,I,J,K,L)
      DIMENSION R(3),EL(3,3),RR(3)
      CCC(J,K,L,M)=3.*CIJKL(RR,J,K,L,M)
      DO 300 M=1,3
  300 RR(M)=R(M)
      SUM=0.
      IF(I.NE.L)GO TO 320
      SUM=SUM+AIJ(R,J,K)
      DO 310 M=1,3
      DO 310 N=1,3
  310 SUM=SUM+EL(M,N)*CCC(J,K,M,N)
  320 IF(J.NE.L)GO TO 340
      SUM=SUM+AIJ(R,I,K)
      DO 330 M=1,3
      DO 330 N=1,3
  330 SUM=SUM+EL(M,N)*CCC(I,K,M,N)
  340 D2KDS=SUM 
      RETURN
      END 
      FUNCTION D2KDT(I,J,K,L) 
      D2KDT=0.
      IF(I.EQ.K.AND.J.EQ.L)D2KDT=1. 
      RETURN
      END 
      FUNCTION D3KDL(R,EL,S,I,J,K,L,M)
      DIMENSION R(3),EL(3,3),S(3,3),RR(3)
      INTEGER P
      ABB(J,K,L,M,N)=2.*AIJ(RR,J,L)*BIJK(RR,K,M,N)
      AAB(I,J,K,L,M,N,NN)=AIJ(RR,I,L)*AIJ(RR,J,N)*BIJK(RR,K,M,NN)
      PP(I,J,K,L,M)=PIJKLM(RR,I,J,K,L,M)+PIJKLM(RR,I,J,K,M,L) 
      DO 300 N=1,3
  300 RR(N)=R(N)
      SUM=0.
      DO 310 N=1,3
      SUM=SUM+S(N,I)*(ABB(J,K,L,M,N)+ABB(K,J,L,M,N)+ABB(J,K,M,L,N)+ABB(K
     1    ,J,M,L,N))+S(N,J)*(ABB(I,K,L,M,N)+ABB(K,I,L,M,N)+ABB(I,K,M,L,N
     2    )+ABB(K,I,M,L,N))+S(N,K)*(ABB(I,J,L,M,N)+ABB(J,I,L,M,N)+ABB(I,
     3    J,M,L,N)+ABB(J,I,M,L,N))
      DO 310 P=1,3
      SUM=SUM+S(N,J)*S(P,K)*PP(I,N,P,L,M)+S(N,I)*S(P,K)*PP(J,N,P,L,M) 
     1    +S(N,I)*S(P,J)*PP(K,N,P,L,M)+2.*EL(N,P)*(AAB(I,J,K,L,M,N,P) 
     2    +AAB(I,K,J,L,M,N,P)+AAB(J,K,I,L,M,N,P)+AAB(I,J,K,N,P,L,M) 
     3    +AAB(I,K,J,N,P,L,M)+AAB(J,K,I,N,P,L,M)) 
  310 CONTINUE
      D3KDL=SUM 
      IF(L.EQ.M)D3KDL=SUM/2.
      RETURN
      END 
      FUNCTION D3KDS(R,EL,S,I,J,K,LP,MP)
      DIMENSION R(3),EL(3,3),S(3,3),RR(3)
      INTEGER P 
      BB(K,L,M)=2.*BIJK(RR,K,L,M)
      ABB(J,K,L,M,N)=2.*AIJ(RR,J,L)*BIJK(RR,K,M,N)
      DO 300 M=1,3
  300 RR(M)=R(M)
      SUM=0.
      DO 330 L=1,3
      DO 330 M=1,3
      IF(L.EQ.LP.AND.I.EQ.MP)SUM=SUM+S(M,J)*BB(K,L,M)+S(M,K)*BB(J,L,M)
      IF(M.EQ.LP.AND.J.EQ.MP)SUM=SUM+S(L,I)*BB(K,L,M) 
      IF(M.EQ.LP.AND.K.EQ.MP)SUM=SUM+S(L,I)*BB(J,L,M)+S(L,J)*BB(I,L,M)
      IF(L.EQ.LP.AND.J.EQ.MP)SUM=SUM+S(M,K)*BB(I,L,M) 
      DO 330 N=1,3
      IF(L.NE.LP)GO TO 310
      IF(I.EQ.MP)SUM=SUM+EL(M,N)*(ABB(J,K,M,L,N)+ABB(K,J,M,L,N))
      IF(J.EQ.MP)SUM=SUM+EL(M,N)*(ABB(I,K,M,L,N)+ABB(K,I,M,L,N))
      IF(K.EQ.MP)SUM=SUM+EL(M,N)*(ABB(I,J,M,L,N)+ABB(J,I,M,L,N))
  310 DO 330 P=1,3
      IF(L.NE.LP)GO TO 320
      IF(J.EQ.MP)SUM=SUM+EL(N,P)*S(M,K)*PIJKLM(R,I,L,M,N,P) 
      IF(I.EQ.MP)SUM=SUM+EL(N,P)*(S(M,K)*PIJKLM(R,J,L,M,N,P)+S(M,J)*
     1               PIJKLM(R,K,L,M,N,P)) 
  320 IF(M.NE.LP)GO TO 330
      IF(K.EQ.MP)SUM=SUM+EL(N,P)*(S(L,J)*PIJKLM(R,I,L,M,N,P)+S(L,I)*
     1               PIJKLM(R,J,L,M,N,P)) 
      IF(J.EQ.MP)SUM=SUM+EL(N,P)*S(L,I)*PIJKLM(R,K,L,M,N,P) 
  330 CONTINUE
      D3KDS=SUM 
      RETURN
      END 
      FUNCTION EIJ(EL,I,J)
      DIMENSION EL(3,3) 
      IF(I.EQ.J)GO TO 300
      EIJ=.5*EL(I,J)-(1./24.)*(EL4(EL,I,1,1,J)+EL4(EL,I,2,2,J)+EL4(EL,
     1 I,3,3,J))
      RETURN
  300 K=MOD(I,3)+1
      L=MOD(I+1,3)+1
      EIJ=1.-.5*(EL(K,K)+EL(L,L))+(1./24.)*(EL4(EL,I,I,K,K)+EL4(EL,I,I, 
     1 L,L)+3.*(EL(K,K)**2+EL(L,L)**2)+2.*EL4(EL,K,K,L,L))
      RETURN
      END 
      FUNCTION EL4(EL,I,J,K,L)
      DIMENSION EL(3,3) 
      EL4=EL(I,J)*EL(K,L)+EL(I,K)*EL(J,L)+EL(I,L)*EL(J,K) 
      RETURN
      END 
      SUBROUTINE EULER(PHI,CHI,OMEGA,TRANS) 
C SUBROUTINE TO GENERATE THE MATRIX B-INVERSE FROM THE EULERIAN 
C  ANGLES PHI,CHI,AND OMEGA. 
      DATA RAD/57.2957792/
      DIMENSION TRANS(3,3)
      ASSIGN 300 TO KGO
      GO TO 360
  300 TRANS(1,1)=CP*CX*CO-SP*SO 
      TRANS(2,1)=CP*CX*SO+SP*CO 
      TRANS(3,1)=-CP*SX 
      TRANS(1,2)=-SP*CX*CO-CP*SO
      TRANS(2,2)=-SP*CX*SO+CP*CO
      TRANS(3,2)=SP*SX
      TRANS(1,3)=SX*CO
      TRANS(2,3)=SX*SO
      TRANS(3,3)=CX 
      RETURN
      ENTRY DTDPHI(PHI,CHI,OMEGA,TRANS) 
      ASSIGN 310 TO KGO
      GO TO 360
  310 TRANS(1,1)=(-SP*CX*CO-CP*SO)/RAD
      TRANS(2,1)=(-SP*CX*SO+CP*CO)/RAD
      TRANS(3,1)=SP*SX/RAD
      TRANS(1,2)=(-CP*CX*CO+SP*SO)/RAD
      TRANS(2,2)=(-CP*CX*SO-SP*CO)/RAD
      TRANS(3,2)=CP*SX/RAD
      DO 320 I=1,3
  320 TRANS(I,3)=0. 
      RETURN
      ENTRY DTDCHI(PHI,CHI,OMEGA,TRANS) 
      ASSIGN 330 TO KGO
      GO TO 360
  330 TRANS(1,1)=-CP*SX*CO/RAD
      TRANS(2,1)=-CP*SX*SO/RAD
      TRANS(3,1)=-CP*CX/RAD
      TRANS(1,2)=SP*SX*CO/RAD
      TRANS(2,2)=SP*SX*SO/RAD
      TRANS(3,2)=SP*CX/RAD
      TRANS(1,3)=CX*CO/RAD
      TRANS(2,3)=CX*SO/RAD
      TRANS(3,3)=-SX/RAD
      RETURN
      ENTRY DTDOMG(PHI,CHI,OMEGA,TRANS) 
      ASSIGN 340 TO KGO
      GO TO 360
  340 TRANS(1,1)=(-CP*CX*SO-SP*CO)/RAD
      TRANS(2,1)=(CP*CX*CO-SP*SO)/RAD
      TRANS(1,2)=(SP*CX*SO-CP*CO)/RAD
      TRANS(2,2)=(-SP*CX*CO-CP*SO)/RAD
      TRANS(1,3)=-SX*SO/RAD
      TRANS(2,3)=SX*CO/RAD
      DO 350 I=1,3
  350 TRANS(3,I)=0. 
      RETURN
  360 CP=COS(PHI/RAD) 
      CX=COS(CHI/RAD) 
      CO=COS(OMEGA/RAD) 
      SP=SIN(PHI/RAD) 
      SX=SIN(CHI/RAD) 
      SO=SIN(OMEGA/RAD) 
      GO TO KGO(300,310,330,340)
      END 
      FUNCTION FIJ(THETA,EL,RP,RP2,I,J) 
      DIMENSION EL(3,3),RP(3),CON(3)
      IF(RP2.LT.0.01.OR.THETA.LT.0.000001)GO TO 300
      K=MOD(I,3)+1
      L=MOD(I+1,3)+1
      KK=MOD(J,3)+1
      LL=MOD(J+1,3)+1
      M2=MOD(J-I+2,3)+1
      M3=MOD(J-I+3,3)+1
      M4=MOD(J-I+4,3)+1
      CON1=1.-THETA/4.-.5*(EL(K,K)+EL(L,L)) 
      CON2=.5*EL(I,K) 
      CON3=.5*EL(I,L) 
      CON(M3)=-(RP(KK)**2+RP(LL)**2)
      CON(M4)=RP(J)*RP(KK)
      CON(M2)=RP(J)*RP(LL)
      FIJ=(CON1*CON(1)+CON2*CON(2)+CON3*CON(3))*THETA/(2.*RP2)
      RETURN
  300 FIJ=0.
      RETURN
      END 
      FUNCTION IPAK(I1,I2,I3,I4)
C Pack Dependency Information in a single word, 8 bits/indicator
      i=ior(ishft(i1,8),i2)
      i=ior(ishft(i,8),i3)
      ipak=ior(ishft(i,8),i4)
c If ISHFT and/or IOR are not available, comment out the above three lines
c   and uncomment the following line
c     IPAK=((I1*256+I2)*256+I3)*256+I4
      RETURN
      END 
      SUBROUTINE MULMM(A,B,C) 
C     MULTIPLY 2 3X3 MATRICES TOGETHER. 
      DIMENSION A(3,3),B(3,3),C(3,3),D(3),E(3)
      DO 300 I=1,3
      DO 300 J=1,3
      C(I,J)=0. 
      DO 300 K=1,3
  300 C(I,J)=C(I,J)+A(I,K)*B(K,J) 
      RETURN
      ENTRY MULMV(A,D,E)
      DO 310 I=1,3
      E(I)=0. 
      DO 310 J=1,3
  310 E(I)=E(I)+A(I,J)*D(J) 
      RETURN
      END 
      FUNCTION PIJKLM(R,I,J,K,L,M)
      DIMENSION R(3)
      PIJKLM=12.*DIJKLM(R,I,J,K,L,M) 
      RETURN
      END 
      FUNCTION SKTLS(R,I,J,T,EL,S)
      DIMENSION R(3),T(3,3),EL(3,3),S(3,3),RR(3)
      B(I1,I2,I3)=BIJK(RR,I1,I2,I3)
      C(I1,I2,I3,I4)=CIJKL(RR,I1,I2,I3,I4)
      DO 300 M=1,3
  300 RR(M)=R(M)
      SUM=T(I,J)
      DO 310 K=1,3
      SUM=SUM+AIJ(R,J,K)*S(K,I)+AIJ(R,I,K)*S(K,J)
      DO 310 L=1,3
      SUM=SUM+AIJ(R,I,K)*AIJ(R,J,L)*EL(K,L)
      DO 310 M=1,3
      SUM=SUM+3.*EL(L,M)*(S(K,I)*C(J,K,L,M)+S(K,J)*C(I,K,L,M))
      DO 310 N=1,3
  310 SUM=SUM+EL(K,L)*EL(M,N)*(3.*(AIJ(R,I,K)*C(J,L,M,N)+AIJ(R,J,K)
     1    *C(I,L,M,N))+2.*B(I,K,M)*B(J,L,N))
      SKTLS=SUM
      RETURN
      END
      SUBROUTINE SOLVE(A,B,C) 
      DIMENSION A(3,3),B(3),C(3)
      DELTA=DETERM(A(1,1),A(1,2),A(1,3))
      C(1)=DETERM(B(1),A(1,2),A(1,3))/DELTA 
      C(2)=DETERM(A(1,1),B(1),A(1,3))/DELTA 
      C(3)=DETERM(A(1,1),A(1,2),B(1))/DELTA 
      RETURN
      END 
      SUBROUTINE SPVAL
C     DUMMY SUBROUTINE SPVAL. THIS SUBROUTINE MAY BE
C     USER WRITTEN TO SUPPLY THE VALUES OF CRYSTALLO-
C     GRAPHIC PARAMETERS THAT ARE ARBITRARY FUNCTIONS
C     OF USER DEFINED SPECIAL PARAMETERS. THE DERIVATIVES
C     OF THE CRYSTALLOGRAPHIC PARAMETERS WITH RESPECT TO
C     THE SPECIAL PARAMETERS ARE PROVIDED BY SPDERI.
      RETURN
      ENTRY SPDERI
      RETURN
      END
      SUBROUTINE SYMINV (A,NRANK,SCR,NCB,B,EPS,DET,JRANK)
C
C THIS ROUTINE REVISED AND/OR RESEQUENCED OCT. 28, 1968
C
C        ARGUMENT LIST
C          A = INPUT AND OUTPUT ARRAY - STORED BY ROWS WITH ONLY THE
C               UNIQUE PORTION IN STORAGE
C          NRANK = INPUT SIZE OF SYSTEM OF EQUATIONS
C          SCR = A SCRATCH ARRAY OF DIMENSION NRANK
C          B = MATRIX TO BE PREMULTIPLIED BY THE INVERSE OF A
C          NCB = NUMBER OF COLUMNS OF B
C          EPS = SIZE OF SMALLEST ELEMENT TO BE USED AS A PIVOT
C          DET = OUTPUT VALUE OF THE DETERMINANT OF A
C          JRANK = ACTUAL RANK OF A ( THIS IS AN OUTPUT VALUE.)
C
      DIMENSION A(1),B(NRANK,1),SCR(1)
C        INITIALIZE
      DET=1.0
      DO 300 I=1,NRANK
  300 SCR(I)=1.0
C        START MASTER LOOP
      DO 480 I=1,NRANK
C        SEARCH FOR PIVOT TESTING ONLY THE DIAGONAL ELEMENTS
      BIG=0.0
C        PRESET INDEX AND STEP FOR DIAGONAL ELEMENTS
      J=-NRANK
      KJ=NRANK+1
      DO 330 K=1,NRANK
C        AUGMENT INDEX AND REDUCE STEP
      J=J+KJ
      KJ=KJ-1
C        CHECK IF THIS ELEMENT ALREADY USED AS PIVOT
      IF(SCR(K))330,310,310
C        J INDEXES A(K,K)
  310 TEMP=ABS(A(J))
      IF(TEMP-BIG)330,330,320
  320 BIG=TEMP
      KK=K
      KL=J
  330 CONTINUE
C        COMPARE LARGEST DIAGONAL ELEMENT WITH EPS
      IF(BIG-EPS)340,350,350
C        SINGULARITY RETURN
  340 JRANK=I-1
      RETURN
  350 SCR(KK)=-1.0
C        KL INDEXES A(KK,KK)
      DET=DET*A(KL)
      Q=1.0/A(KL)
C        SET UP INDEXING FOR ELIMINATION STEPS
      J=KK-NRANK
      NJ=NRANK+1
      K=0
      IK=KK-1
      JK=KK+1
      IF(IK)410,410,360
  360 DO 400 N=1,IK
      NJ=NJ-1
      J=J+NJ
C        J INDEXES A(KK,N)
      TEMP=Q*A(J)
      NI=J
      MJ=NJ
      DO 370 M=N,IK
      K=K+1
C        J INDEXES A(KK,M)
C        K INDEXES A(M,N)
      A(K)=A(K)-TEMP*SCR(M)*A(J)
      MJ=MJ-1
  370 J=J+MJ
      J=NI
C        J INDEXES A(KK,N)
      A(J)=TEMP
      K=K+1
      L=KL
      IF(JK-NRANK)380,380,400
  380 DO 390 M=JK,NRANK
      L=L+1
      K=K+1
C        K INDEXES A(M,N)
C        L INDEXES A(M,KK)
  390 A(K)=A(K)-TEMP*A(L)
  400 CONTINUE
  410 K=KL+NRANK-KK
      J=KL
      IF(JK-NRANK)420,420,450
  420 DO 440 N=JK,NRANK
      J=J+1
C        J INDEXES A(N,KK)
      TEMP=SCR(N)*A(J)*Q
      NI=J
      DO 430 M=N,NRANK
      K=K+1
C        K INDEXES A(M,N)
C        J INDEXES A(M,KK)
      A(K)=A(K)-TEMP*A(J)
  430 J=J+1
  440 J=NI
  450 K=KL
C        K INDEXES A(KK,KK)
      A(K)=Q
      IF(JK-NRANK)460,460,480
  460 DO 470 M=JK,NRANK
      K=K+1
C        K INDEXES A(M,KK)
  470 A(K)=-Q*A(K)
  480 CONTINUE
C        INVERSION COMPLETE, EXIT IF INVERSION ONLY
      IF(NCB)560,560,490
C        PREMULTIPLY THE NCB COLUMNS OF B BY A INVERSE
  490 DO 550 K=1,NCB
      DO 530 I=1,NRANK
      Q=0.0
      L=I-NRANK
      KK=NRANK+1
      DO 500 J=1,I
      KK=KK-1
      L=L+KK
C        L INDEXES A(I,J)
  500 Q=Q+A(L)*B(J,K)
      IF(I-NRANK)510,530,530
  510 KJ=I+1
      DO 520 J=KJ,NRANK
      L=L+1
C        L INDEXES A(J,I)
  520 Q=Q+A(L)*B(J,K)
  530 SCR(I)=Q
      DO 540 J=1,NRANK
  540 B(J,K)=SCR(J)
  550 CONTINUE
C        SET NON-SINGULAR FLAG AND EXIT
  560 JRANK=NRANK
      RETURN
      END
      SUBROUTINE TFM2(T,A,B)
      INCLUDE 'RFINE.PRM'
      INCLUDE 'E.INC'
      DIMENSION A(6),B(6),C(10),D(10),T(3,3),BTEM(3,3),GTEM(3,3,3),J(3) 
      EQUIVALENCE(BTEM(1,1),GTEM(1,1,1)),(J1,J(1)),(J2,J(2)),(J3,J(3))
      DATA PI/3.14159265/ 
      DO 310 I=1,6
      DO 300 K=1,2
  300 J(K)=IND2(K,I)
      BTEM(J1,J2)=A(I)
      BTEM(J2,J1)=A(I)
  310 CONTINUE
      DO 340 I=1,6
      B(I)=0. 
      DO 320 K=1,2
  320 J(K)=IND2(K,I)
      DO 330 K=1,3
      DO 330 L=1,3
  330 B(I)=B(I)+T(J1,K)*T(J2,L)*BTEM(K,L) 
  340 B(I)=B(I)*2.*PI**2
      RETURN
      ENTRY TFM3(T,C,D) 
      DO 360 I=1,10
      DO 350 K=1,3
  350 J(K)=IND3(K,I)
      GTEM(J1,J2,J3)=C(I) 
      GTEM(J2,J3,J1)=C(I) 
      GTEM(J3,J1,J2)=C(I) 
      GTEM(J3,J2,J1)=C(I) 
      GTEM(J2,J1,J3)=C(I) 
  360 GTEM(J1,J3,J2)=C(I) 
      DO 390 I=1,10
      D(I)=0. 
      DO 370 K=1,3
  370 J(K)=IND3(K,I)
      DO 380 K=1,3
      DO 380 L=1,3
      DO 380 M=1,3
  380 D(I)=D(I)+T(J1,K)*T(J2,L)*T(J3,M)*GTEM(K,L,M) 
  390 D(I)=D(I)*4.*PI**3/3. 
      RETURN
      END 
      FUNCTION TKTLS(R,I,J,K,EL,S)
      DIMENSION R(3),EL(3,3),S(3,3) 
      INTEGER P 
      SUM=0.
      DO 300 L=1,3
      DO 300 M=1,3
      SUM=SUM+2.*(S(L,I)*S(M,J)*BIJK(R,K,L,M)
     1           +S(L,I)*S(M,K)*BIJK(R,J,L,M)
     2           +S(L,J)*S(M,K)*BIJK(R,I,L,M))
      DO 300 N=1,3
      SUM=SUM+2.*EL(M,N)*
     1    (S(L,I)*(AIJ(R,J,M)*BIJK(R,K,L,N)+AIJ(R,K,M)*BIJK(R,J,L,N))
     2    +S(L,J)*(AIJ(R,I,M)*BIJK(R,K,L,N)+AIJ(R,K,M)*BIJK(R,I,L,N))
     3    +S(L,K)*(AIJ(R,I,M)*BIJK(R,J,L,N)+AIJ(R,J,M)*BIJK(R,I,L,N))) 
      DO 300 P=1,3
      SUM=SUM+EL(N,P)*(S(L,J)*S(M,K)*PIJKLM(R,I,L,M,N,P)+S(L,I)*S(M,K)* 
     1    PIJKLM(R,J,L,M,N,P)+S(L,I)*S(M,J)*PIJKLM(R,K,L,M,N,P)
     2    +2.*EL(L,M)*(AIJ(R,I,L)*AIJ(R,J,N)*BIJK(R,K,M,P)
     3    +AIJ(R,I,L)*AIJ(R,K,N)*BIJK(R,J,M,P)
     4    +AIJ(R,J,L)*AIJ(R,K,N)*BIJK(R,I,M,P)))
  300 CONTINUE
      TKTLS=SUM 
      RETURN
      END 
      SUBROUTINE TRNFRM 
C     SUBROUTINE CALLED ONCE PER JOB, AFTER CELL CONSTANTS HAVE BEEN
C     READ IN, TO SET UP TRANSFORMS TO TAKE VECTORS FROM CRYSTAL AXES 
C     TO AN ORTHONORMAL COORDINATE SYSTEM, AND TO TAKE THE SECOND AND 
C     THIRD CUMULANTS BACK TO THE CRYSTAL SYSTEM. 
      INCLUDE 'RFINE.PRM'
      INCLUDE 'A.INC'
      INCLUDE 'XPAR.INC'
C
      ATOX(1,1)=SQRT(G(1,1))
      ATOX(1,2)=G(1,2)/ATOX(1,1)
      ATOX(1,3)=G(1,3)/ATOX(1,1)
      ATOX(2,1)=0.
      ATOX(2,2)=SQRT(G(2,2)-ATOX(1,2)**2) 
      ATOX(2,3)=(G(2,3)-ATOX(1,2)*ATOX(1,3))/ATOX(2,2)
      ATOX(3,1)=0.
      ATOX(3,2)=0.
      ATOX(3,3)=SQRT(G(3,3)-ATOX(1,3)**2-ATOX(2,3)**2)
      DO 300 I=1,3
  300 XTOA(I,I)=1./ATOX(I,I)
      XTOA(1,2)=-ATOX(1,2)/(ATOX(1,1)*ATOX(2,2))
      XTOA(2,3)=-ATOX(2,3)/(ATOX(2,2)*ATOX(3,3))
      XTOA(1,3)=(ATOX(1,2)*ATOX(2,3)/ATOX(2,2)-ATOX(1,3))/(ATOX(1,1)*ATO
     1X(3,3)) 
      XTOA(2,1)=0.
      XTOA(3,1)=0.
      XTOA(3,2)=0.
      RETURN
      END 
      SUBROUTINE TRNSPS(A,B)
C     TRANSPOSE MATRIX A AND STORE IN B.
      DIMENSION A(3,3),B(3,3) 
      DO 300 I=1,3
      DO 300 J=1,3
  300 B(J,I)=A(I,J) 
      RETURN
      END 
      SUBROUTINE VXV(V1,V2,V3)
      DIMENSIONV1(3),V2(3),V3(3)
      V3(1)=V1(2)*V2(3)-V1(3)*V2(2) 
      V3(2)=V1(3)*V2(1)-V1(1)*V2(3) 
      V3(3)=V1(1)*V2(2)-V1(2)*V2(1) 
      RETURN
      END 
      SUBROUTINE WEIGHT 
  
C THIS ROUTINE REVISED AND/OR RESEQUENCED ON OCT. 28, 1968
  
C        DUMMY ROUTINE
      RETURN
      ENTRY REJECT
C        DUMMY ROUTINE
      RETURN
      END 
