      SUBROUTINE OPLOT( X, Y, IND)
C
C  HERE WE SCALE FROM ORTEP DRAWING COORDINATES
C    TO COORDINATES THAT RANGE FROM 0.0 TO 1.0
C
      IMPLICIT REAL( A-H, O-Z)
      COMMON /PLOTS/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,CMIN,CMAX,SCALE
C
      IF ( IND .EQ. 99 ) THEN
         CALL PLOT( X, Y , IND )
      ELSE
         XX = X  / SCALE
         YY = Y / SCALE
         CALL PLOT( XX , YY , IND )
      ENDIF
      RETURN
      END
C
      SUBROUTINE OSIM( X, Y, SIZE, CBCD, THETA, N )
C
C  HERE WE SCALE FROM ORTEP DRAWING SYMBOLS
C    TO COORDINATES THAT RANGE FROM 0.0 TO 1.0
C
      IMPLICIT REAL ( A-H, O-Z)
      REAL X, Y, SIZE, THETA, XX, YY, SSIZE, TTHETA
      CHARACTER*6 CBCD
      COMMON /PLOTS/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,CMIN,CMAX,SCALE
C
      SSIZE = SIZE / SCALE
      TTHETA = THETA
      XX = X / SCALE
      YY = Y / SCALE
      CALL SIMBOL( XX, YY, SSIZE, CBCD, TTHETA, N)
      RETURN
      END

      SUBROUTINE ATOM(QA,Z)
C     ATOM COORDINATE SUBROUTINE
      DIMENSION X(3),Z(3)
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT
C *     END OF COMMON   *

      K=QA/100000.0
      IF(K.GT.0)GO TO 117
  109 X(1)=0.0
      X(2)=0.0
      X(3)=0.0
      GO TO 125
  117 IF(K-NATOM)119,119,503
  503 NG=5
      GO TO 325
  119 DO 123 J=1,3
  123 X(J)=P(J,K)
  125 TA=ABS(QA)
      SYM=AMOD(TA,100000.0)
      KT=SYM/100.
      KS=AMOD(SYM,100.)
      IF(KS.LE.NSYM)GO TO 203
  403 NG=4
      GO TO 325
  203 IF(KS)403,205,213
  205 Z(1)=X(1)
      Z(2)=X(2)
      Z(3)=X(3)
      GO TO 311
  213 DO 223 K=1,3
      Z(K)=TS(K,KS)
      DO 223 J=1,3
  223 Z(K)=Z(K)+FS(J,K,KS)*X(J)
  311 IF(KT)403,325,313
  313 IF(KT.NE.555)GO TO 317
  315 KSYM=KS
      GO TO 325
  317 K1=KT/100
      K=KT-100*K1
      K2=K/10
      K3=K-10*K2
      Z(1)=Z(1)+(K1-5)
      Z(2)=Z(2)+(K2-5)
      Z(3)=Z(3)+(K3-5)
  325 RETURN
      END
      SUBROUTINE AXEQB(A1,X,B1,JJJ)
C     ***** SOLUTION OF MATRIX EQUATION AX=B FOR X *****
C     ***** USES METHOD OF TRIANGULAR ELIMINATION *****
C     ***** B AND X HAVE DIMENSIONS (3,JJJ),A IS ALWAYS (3,3)
C     ***** TO INVERT A MAKE B 3 BY 3 IDENITY MATRIX *****
      DIMENSION A1(3,3),A(3,3),B(3,3),B1(3,3),X(3,3)
      NV=JJJ
C     ***** TRANSFER DATA *****
      DO 2 I=1,3
      DO 2 J=1,3
      IF(NV-J)2,1,1
    1 B(I,J)=B1(I,J)
    2 A(I,J)=A1(I,J)
C     ***** TRIANGULARIZE MATRIX A *****
      DO 17 I=1,2
      S=0.0
      DO 4 J=I,3
      R=ABS(A(J,I))
      IF(R-S)4,3,3
    3 S=R
      L=J
    4 CONTINUE
      IF(L-I)5,10,5
    5 DO 6 J=I,3
      S=A(I,J)
      A(I,J)=A(L,J)
    6 A(L,J)=S
      DO 8 J=1,NV
      S=B(I,J)
      B(I,J)=B(L,J)
    8 B(L,J)=S
   10 TEM=A(I,I)
      IF(TEM)11,17,11
   11 IPO=I+1
      DO 16 J=IPO,3
      IF(A(J,I))12,16,12
   12 S=A(J,I)/TEM
      A(J,I)=0.0
      DO 13 K=IPO,3
   13 A(J,K)=A(J,K)-A(I,K)*S
      DO 15 K=1,NV
   15 B(J,K)=B(J,K)-B(I,K)*S
   16 CONTINUE
   17 CONTINUE
C     ***** MODIFY SINGULAR MATRIX *****
      DO 20 I=1,3
      IF(A(I,I))20,19,20
   19 A(I,I)=AMAX1(1.E-25,AMAX1(A(1,1),A(2,2),A(3,3))*1.E-15)
   20 CONTINUE
      DO 24 K=1,NV
      DO 24 I=1,3
      N=4-I
      M=N+1
      TEM=B(N,K)
      IF(3-M)23,21,21
   21 DO 22 J=M,3
   22 TEM=TEM-A(N,J)*B(J,K)
   23 B(N,K)=TEM/A(N,N)
   24 X(N,K)=B(N,K)
      RETURN
      END
      SUBROUTINE AXES(U,V,X,ITYPE)
C     ***** STORE THREE ORTHOGONAL VECTORS EACH 1 ANGSTROM LONG *****
C     ***** ITYPE .GT.0 FOR CARTESIAN,.LE.0 FOR TRICLINIC *****
C     *****IABS(ITYPE)=1 W(1)=U,W(2)=(UXV),W(3)=UX(UXV) *****
C     *****IABS(ITYPE)=2 W(1)=U,W(2)=(UXV)XU,W(3)=(UXV) *****
C     ***** ITYPE=0 W(1)=A,W(2)=(AXB)XA,W(3)=(AXB), ABC=CELL VECTORS ***
      DIMENSION U(3),V(3),W(3,3),X(3,3)
      IT=ITYPE
      IF(IT)115,105,115
  105 U(1)=1.
      U(2)=0.
      U(3)=0.
      V(1)=0.
      V(2)=1.
      V(3)=0.
  115 DO 125 J=1,3
  125 W(J,1)=U(J)
      IF(IABS(IT)-1)145,135,145
  135 CALL NORM(U,V,W(1,2),IT)
      CALL NORM(U,W(1,2),W(1,3),IT)
      GO TO 155
  145 CALL NORM(U,V,W(1,3),IT)
      CALL NORM(W(1,3),U,W(1,2),IT)
  155 DO 195 I=1,3
      IF(IT)165,165,175
  165 IC=-1
      GO TO 195
  175 IC=1
  195 CALL UNIT(W(1,I),X(1,I),IC)
      RETURN
      END
      SUBROUTINE BOND(Z1,Z2,NB,NA1,NA2)
      CHARACTER*80 COMMNT
      DIMENSION B(3,3),E(3,3),S(3,3),U(3,3),VUE(3)
      DIMENSION V7(3),W(13,2),X(3),Z(3)
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)
      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT
C *     END OF COMMON   *
      CHARACTER HEADER*80, CHEM(NUMATM)*6
      COMMON /STRING/ HEADER, CHEM
      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
C     ***** OBTAIN POSITIONAL PARAMETERS *****
      ICOLOR = 1
      NG1=0
      DO 105 J=1,26
  105 W(J,1)=0.
      W(1,1)=Z1
      W(1,2)=Z2
      DO 135 I=1,2
         CALL XYZ(W(1,I),W(4,I),2)
         IF(NG)125,110,125
  110    DO 115 J=1,3
  115    W(J+6,I)=XT(J)
         K=W(1,I)/100000.
         L=W(1,I)-(K)*100000.
         CALL PLTXY(W(4,I),W(2,I))
         IF(EDGE-BRDR*.25)120,128,128
  120    NG=10
  125    NG1=1
         IF (DEBUGO) THEN
            WRITE ( COMMNT,136) CHEM(K),K,L,(W(J,I),J=2,9)
            CALL DEBUGR( COMMNT(1:78) )
         ENDIF
         CALL ERPNT(W(1,I),800)
         GO TO 135
  128    IF(NJ2-10)130,135,135
  130    IF (DEBUGO) THEN
            WRITE ( COMMNT,136)CHEM(K),K,L,(W(J,I),J=2,9)
            CALL DEBUGR( COMMNT(1:78) )
         ENDIF
  135 CONTINUE
  136 FORMAT(1H 10X,A6,3H  (I3,1H,I5,4H)   2F8.2,5X,3F8.3,13X,3F8.4)
      IF(NG1)999,140,999
  140 CALL DIFV(W(7,1),W(7,2),V7)
      DIST=SQRT(VMV(V7,AA,V7))
      IF(MOD(NJ2,10)-3)143,142,142
  142 HGT=SCL*.12
C  MAKE CENTERED SYMBOLS VERY SMALL
      HGT=.0000001
      X(1)=W(2,1)
      X(2)=W(3,1)
      CALL OSIM(X(1),X(2),HGT,MOD(IFIX(W(1,1)/100000.),10),0.,-1)
      X(1)=W(2,2)
      X(2)=W(3,2)
      CALL OSIM(X(1),X(2),HGT,MOD(IFIX(W(1,2)/100000.),10),0.,-2)
      HGT=SCL*0.12
      GO TO 570
  143 KODE=KD(5,NB)
      IF(KODE)145,144,146
  144 NBND=0
      GO TO 148
  145 KODE=-KODE
  146 NBND=128/2**KODE
C     ***** FIND UPPERMOST ATOM PUT IN POSITION ONE *****
  148 IF(VIEW)152,150,152
  150 W(12,1)=1.
      W(12,2)=1.
      IF(W(6,1)-W(6,2))165,175,175
C     *****VECTOR FROM ATOM TO VIEWPOINT *****
  152 DO 160 I=1,2
         DO 155 J=10,12
  155    W(J,I)=-W(J-6,I)
         W(12,I)=W(12,I)+VIEW
C     ***** DISTANCE SQUARED TO VIEWPOINT *****
  160 W(13,I)=VV(W(10,I),W(10,I))
      IF(W(13,2)-W(13,1))165,175,175
C     ***** SWITCH ATOMS *****
  165 DO 170 J=1,13
         T1=W(J,1)
         W(J,1)=W(J,2)
  170 W(J,2)=T1
C     ***** FORM IDEMFACTOR MATRIX *****
  175 DO 180 J=1,3
         E(J,J)=1.
         E(J+1,1)=0.
  180 E(J+5,1)=0.
C     ***** FORM VECTOR SET RADIAL TO BOND *****
      CALL DIFV(W(4,2),W(4,1),DA(1,3))
      CALL UNIT(DA(1,3),V3,1)
C     ***** UNIT VECTOR FROM BOND MIDPOINT TO REFERENCE VIEWPOINT *****
      DO 183 I=1,3
         V2(I)=0.0
         DO 181 J=1,3
  181    V2(I)=V2(I)+AAREV(J,3)*WRKV(J,I)
         IF(VIEW)183,183,182
  182    V2(I)=V2(I)*VIEW-0.5*(W(I+3,1)+W(I+3,2))
  183 CONTINUE
      CALL UNIT(V2,V2,1)
      T6=ABS(VV(V3,V2))
      IF(.9994-T6)185,185,187
C     ***** ALTERNATE CALC IF BOND IS ALONG REFERENCE VIEW DIRECTION ***
  185 DO 186 J=1,3
  186 V2(J)=W(J+9,1)+W(J+9,2)
      CALL UNIT(V2,V2,1)
      T6=ABS(VV(V3,V2))
      IF(.9994-T6)390,390,187
  187 CALL AXES(V3,V2,B,1)
  188 T1=CD(3,NB)/SCAL2
      DO 190 J=1,3
         DA(J,1)=-B(J,2)*T1
  190 DA(J,2)=-B(J,3)*T1
      IF(NBND)500,500,195
  195 CALL RADIAL(3)
C     ***** DERIVE QUADRICS FOR EACH ATOM *****
      DO 380 II=1,2
         CALL PAXES(W(1,II),2)
         IF(NG)205,210,205
  205    CALL ERPNT(W(1,II),800)
         GO TO 999
C     ***** DOES BOND GO TO ELLIPSOID OR TO ENVELOPE *****
  210    T1=3-II*2
         DO 212 J=1,3
            V3(J)=V3(J)*T1
  212    VUE(J)=0.
         IF(KD(5,NB))260,260,215
  215    IF(VMV(V3,Q,W(10,II)))220,260,260
  220    IBND=0
         IF(VIEW)240,240,225
C     ***** DERIVE TANGENT CONE DIRECTLY WITHOUT ROTATING COORDINATES **
  225    T2=-(SCAL2*RMS(1)*RMS(2)*RMS(3))**2
         DO 230 J=1,3
            V1(J)=-W(J+9,II)/SCAL1
            VUE(J)=V1(J)/SCAL2
C     ***** INVERT ELLIPSOID MATRIX *****
            DO 230 K=J,3
            T1=0.0
            DO 228 I=1,3
  228       T1=T1+PAC(J,I)*PAC(K,I)*RMS(I)**2
            U(J,K)=T1
  230    U(K,J)=T1
C     *****  ADD POLARIZED COFACTOR MATRIX TO ELLIPSOID MATRIX *****
         DO 235 J=1,3
            J1=MOD(J,3)+1
            VJ1=V1(J1)
            J2=MOD(J+1,3)+1
            VJ2=V1(J2)
            DO 235 K=J,3
               K1=MOD(K,3)+1
               K2=MOD(K+1,3)+1
               S(J,K)=T2*Q(J,K)+(VJ2*(U(J1,K1)*V1(K2)-U(J1,K2)*V1(K1))
     1                + VJ1*(U(J2,K2)*V1(K1)-U(J2,K1)*V1(K2)))
  235       S(K,J)=S(J,K)
            T5=0.0
            GO TO 300
C     ***** DERIVE TANGENT CYLINDER WITH AXIS ALONG Z *****
  240       T1=-1.0/Q(3,3)
            DO 250 J=1,2
               DO 245 K=1,2
  245          S(K,J)=Q(K,J)+Q(K,3)*Q(J,3)*T1
               S(3,J)=0.0
  250       S(J,3)=0.0
            S(3,3)=0.0
            GO TO 270
C     ***** TRANSFER ELLIPSOID *****
  260       DO 265 J=1,9
  265       S(J,1)=Q(J,1)
            IBND=II
  270       T5=1.
C     ***** CHECK FOR BOND TAPER *****
  300       IF(II-2)305,310,310
  305       RADIUS=1.+T6*TAPER
            GO TO 320
  310       RADIUS=1.-T6*TAPER
  320       CALL MV(S,V3,V4)
            T2=VV(V3,V4)
C     ***** COMPUTE BOND INTERSECTION *****
            KL=5-II-II
            KSTP=4
            IF(NJ2-21)324,322,322
  322       KSTP=32
  324       DO 335 K=1,65,KSTP
            DO 325 J=1,3
               V6(J)=D(J,K)*RADIUS
  325       V5(J)=V6(J)+VUE(J)
            T3=VV(V5,V4)
            T4=T3*T3-T2*(VMV(V5,S,V5)-T5)
            IF(T4)345,330,330
  330       T4=SQRT(T4)
            T1=(T4-T3)/T2
            T3=(-T4-T3)/T2
            L=K+KL-1
            DO 335 J=1,3
               D(J,L)=(V6(J)+T1*V3(J))*SCL
  335       D(J,L+1)=(-V6(J)-T3*V3(J))*SCL
            IF(IBND+21-NJ2)360,338,360
  338       IF(KD(5,NB))360,360,340
C     ***** FOR OVERLAP, MAKE BOND QUADRANGLE TANGENT TO ENVELOPING CONE
  340       T3=VV(VUE,V4)
            T4=T3**2-T2*(VMV(VUE,S,VUE)-T5)
            IF(T4)345,350,350
  345       NG=13
            CALL ERPNT(W(1,II),800)
            GO TO 999
  350       T1=(SQRT(T4)-T3)/T2
            DO 355 J=1,3
               T4=(T1*V3(J)*SCL-0.5*(D(J,KL)+D(J,KL+64)))*1.001
               D(J,KL)=D(J,KL)+T4
  355       D(J,KL+64)=D(J,KL+64)+T4
  360       CALL PROJ(D(1,KL),DP(1,II),W(4,II),XO,VIEW,1,65,KSTP)
            IF(IBND-1)370,365,370
  365       CALL PROJ(D(1,KL+5),DP(1,II+68),W(4,II),XO,VIEW,1,61,KSTP)
            GO TO 380
C     ***** RETRACE TOP HALF *****
  370       DO 375 K=4,64,4
               L=K+II
               M=L+64
               N=66-L
               IF (N.LT.1) GO TO 375
               DP(1,M)=DP(1,N)
               DP(2,M)=DP(2,N)
  375       CONTINUE
  380 CONTINUE
C     ***** CHECK FOR OVERLAP OR HIDDEN BOND *****
      DO 395 K=1,65,32
         T1=0.
         T2=0.
         DO 385 J=1,2
            T1=T1+(DP(J,K)-W(J+1,1))**2
  385    T2=T2+(DP(J,K+1)-W(J+1,1))**2
         IF(T2-T1)390,390,395
  390    NG=14
         CALL ERPNT(W(1,2),800)
         GO TO 999
  395 CONTINUE
C     ***** CALL OVERLAP ROUTINE *****
      ICQ=0
      CALL LAP800(NA1,NA2,ICQ)
      IF(NJ2-21)400,999,999
  400 IF(ICQ)390,405,405
C     ***** DRAW BOND OUTLINE *****
  405 CALL DRAW(DP(1,1),0.,0.,3, ICOLOR)
      DO 415 K=5,129,4
  415 CALL DRAW(DP(1,K),0.,0.,2, ICOLOR)
      DO 420 K=2,66,4
  420 CALL DRAW(DP(1,K),0.,0.,2, ICOLOR)
      IF(KD(5,NB).GT.1)GOTO1010
      IF(CD(3,NB).GE..01)GOTO1010
      CALL DRAW(DP(1,  K),0.,0.,3, ICOLOR)
      GOTO 500
 1010 CONTINUE
      CALL DRAW(DP(1,65),0.,0.,2, ICOLOR)
C     ***** DRAW BOND DETAIL *****
  425 K=65
  430 K=K-NBND
      IF(K-1)500,500,435
  435 CALL DRAW(DP(1,K),0.,0.,3, ICOLOR)
      CALL DRAW(DP(1,K+1),0.,0.,2, ICOLOR)
      K=K-NBND
      IF(K-1)500,500,440
  440 CALL DRAW(DP(1,K+1),0.,0.,3, ICOLOR)
      CALL DRAW(DP(1,K),0.,0.,2, ICOLOR)
      GO TO 430
  500 HGT=CD(4,NB)
      OFF=CD(5,NB)
      IF(HGT)570,570,510
C     ***** PERSPECTIVE BOND LABEL ROUTINE *****
C     ***** BASE DECISIONS ON REFERENCE SYSTEM *****
  510 K=0
      CALL DIFV(W(7,2),W(7,1),V7)
      CALL VM(V7,AAREV,V1)
      CALL AXES(V1,E(1,3),U,1)
      DO 535 I=1,3
         T1=1.
         IF(I-2)515,515,520
  515    IF(VV(U(1,I),SYMB(1,I)))525,530,530
  520    IF(MOD(K,2))530,525,530
  525    T1=-1.
         K=K+1
  530    DO 535 J=1,3
            U(J,I)=U(J,I)*T1
  535 VT(J,I)=B(J,I)*T1
      DO 540 J=1,3
 540  VT(J,4)=.5*(W(J+3,1)+W(J+3,2))
C     ***** CHECK FOR EXCESS FORESHORTENING *****
      IF(FORE-ABS(U(3,1)))545,550,550
  545 CALL NORM(U(1,2),SYMB(1,3),VT(1,1),1)
      VT(1,3)=SYMB(1,3)
      VT(2,3)=SYMB(2,3)
      VT(3,3)=SYMB(3,3)
      HGT=CD(6,NB)
      OFF=CD(7,NB)
      IF(HGT)550,999,550
  550 T1=CD(8,NB)
      Z(1)=VT(1,4)-HGT*(11.+3.*T1)/7.
      Z(2)=VT(2,4)+OFF-HGT*.5
      Z(3)=VT(3,4)
      XO(3)=Z(3)
      ITILT=1
      I9=T1+2.
      T9=10.**I9
      DISTR=AINT((DIST*T9)+0.5)/T9 +.0001
      CALL NUMBUR(Z,HGT,DISTR,0.,I9)
  570 ITILT=0
      IF(NJ2-10)580,999,999
  580 IF (DEBUGO) THEN
         WRITE ( COMMNT,571)DIST
         CALL DEBUGR( COMMNT(1:20) )
      ENDIF
  571 FORMAT(' DISTANCE =',F8.3 )
  999 RETURN
      END
      SUBROUTINE DIFV(X,Y,Z)
C     VECTOR - VECTOR
C     Z(3)=X(3)-Y(3)
      DIMENSIONX(3),Y(3),Z(3)
      DO 111 I=1,3
  111 Z(I)=X(I)-Y(I)
      RETURN
      END
      SUBROUTINE DRAW(W,DX,DY,NPEN, ICOLOR)
      DIMENSION W(3),X(3),Y(3),Z(3)
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      Y(1)=W(1)+DX
      Y(2)=W(2)+DY
      IF(ITILT)115,140,115
C     ***** ROTATE FOR PERSPECTIVE TITLE *****
  115 Y(3)=XO(3)
      DO 120 I=1,3
  120 Z(I)=Y(I)-VT(I,4)
      DO 130 I=1,3
  130 X(I)=VT(I,1)*Z(1)+VT(I,2)*Z(2)+VT(I,3)*Z(3)+VT(I,4)
      CALL PLTXY(X,Y)
C     ***** CHECK BOUNDRY *****
  140 DO 160 J=1,2
         IF(Y(J)-XLNG(J)+.1)150,150,145
  145    Y(J)=XLNG(J)-.1
  150    IF(Y(J)-.1)155,160,160
  155    Y(J)=.1
  160 CONTINUE
C     ***** CHECK FOR OVERLAP *****
      NCQ=0
      CALL LAPDRW(Y,NPEN,NCQ, ICOLOR)
      IF(NCQ)165,165,170
C     ***** CALL PLOTTING ROUTINE IF NO OVERLAPPING ELEMENTS ARE STORED
  165 CALL SCRIBE(Y,NPEN,LTNO, ICOLOR)
  170 RETURN
      END

      SUBROUTINE EIGEN (W,VALU,VECT)
C     ***** EIGENVALUES AND EIGENVECTORS OF 3X3 MATRIX *****
      DIMENSION W(3,3),VALU(3),VECT(3,3),EA(3,3),B(3,3),V(3),U(3)
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *
C     ***** STATEMENT FUNCTION *****
      PHIF(Z)=((B2-Z)*Z+B1)*Z+B0
C     ***** START OF PROGRAM *****
      ERRND=5.E-7
      SIGMA=0.
      DO 115 J=1,3
      DO 115 I=1,3
      TEM=W(I,J)
      EA(I,J)=TEM
  115 SIGMA=SIGMA+TEM*TEM
C     ***** CHECK FOR NULL MATRIX *****
      IF(SIGMA)230,230,120
  120 SIGMA=SQRT(SIGMA)
C     ***** FORM CHARACTERISTIC EQUATION *****
      B2=EA(1,1)+EA(2,2)+EA(3,3)
      B1=-EA(1,1)*EA(2,2) - EA(1,1)*EA(3,3) - EA(2,2)*EA(3,3) +
     1    EA(1,3)*EA(3,1) + EA(2,3)*EA(3,2) + EA(1,2)*EA(2,1)
      B0= EA(1,1)*EA(2,2)*EA(3,3) + EA(1,2)*EA(2,3)*EA(3,1) +
     1    EA(1,3)*EA(3,2)*EA(2,1) - EA(1,3)*EA(3,1)*EA(2,2) -
     1    EA(1,1)*EA(2,3)*EA(3,2) - EA(1,2)*EA(2,1)*EA(3,3)
C     ***** FIRST ROOT BY BISECTION *****
      X=0.
      Y=SIGMA
      TEM=PHIF(SIGMA)
      VNEW=0.0
      IF(B0)135,250,145
  135 IF(TEM)140,140,165
  140 Y=-Y
      GO TO 165
  145 Y=0.
      X=SIGMA
      IF(TEM)165,165,150
  150 X=-X
C     ***** NOW PHIF(X).LT.0.AND.PHIF(Y).GT.0. *****
  165 VNEW=(X+Y)*.5
      DO 225 I=1,40
  175 IF(PHIF(VNEW))180,250,185
  180 X=VNEW
      GO TO 200
  185 Y=VNEW
  200 VOLD=VNEW
      VNEW=(X+Y)*.5
      TEM=ABS(VOLD-VNEW)
      IF(TEM-ERRND)250,250,205
  205 IF(VOLD)210,225,210
  210 IF (ABS(TEM/VOLD)-ERRND)250,250,225
  225 CONTINUE
C     ***** DID NOT CONVERGE, SET ERROR INDICATOR *****
  230 NG=6
      GO TO 400
C     ***** STORE FIRST ROOT *****
  250 U(3)=VNEW
C     ***** DEFLATE *****
      C1=B2-VNEW
      C0=B1+C1*VNEW
C     ***** SOLVE QUADRATIC *****
      TEM=C1*C1+4.*C0
      IF(TEM)255,265,260
C     ***** IGNORE IMAGINARY COMPONENT OF COMPLEX ROOT *****
  255 TEM=0.
      GO TO 265
  260 TEM=SQRT(TEM)
  265 U(1)=.5*(C1-TEM)
      U(2)=.5*(C1+TEM)
C     ***** SORT ROOTS *****
      DO 275 J=1,2
      IF(U(J)-U(3))275,275,270
  270 TEM=U(J)
      U(J)=U(3)
      U(3)=TEM
  275 CONTINUE
      LLL=-2
      DO 375 III=1,2
C     ***** CHECK FOR MULTIPLE ROOTS *****
      TEM=ERRND*100.
      NG=0
      L=1
      DO 305 I=1,2
      IF(U(I+1)-U(I)-TEM)300,300,290
  290 IF(U(I))295,305,295
  295 IF(ABS((U(I+1)-U(I))/U(I))-TEM)300,300,305
  300 L=L-1
      NG=NG-2*I
  305 CONTINUE
      IF(LLL-L)308,400,400
  308 LLL=L
C     ***** EIGENVECTOR ROUTINE *****
      DO 375 II=1,3
      T1=U(II)
      IF(L)315,310,322
C     ***** TWO VECTORS NULL FOR DOUBLE ROOT *****
  310 IF(NG+5-II)315,322,315
C     ***** ALL VECTORS NULL FOR TRIPLE ROOT *****
  315 DO 320 J=1,3
  320 VECT(J,II)=0.0
      GO TO 375
  322 DO 325 J=1,3
  325 EA(J,J)=W(J,J)-T1
      SMAX=0.0
      DO 355 I=1,3
      I1=1
      IF(I-2)335,335,340
  335 I1=I+1
  340 B(I,1)=EA(I,2)*EA(I1,3)-EA(I,3)*EA(I1,2)
      B(I,2)=EA(I,3)*EA(I1,1)-EA(I,1)*EA(I1,3)
      B(I,3)=EA(I,1)*EA(I1,2)-EA(I,2)*EA(I1,1)
      TEM=B(I,1)**2+B(I,2)**2+B(I,3)**2
      IF(TEM-SMAX)355,355,350
  350 SMAX=TEM
      IMAX=I
  355 CONTINUE
      IF(SMAX)353,353,360
  353 NG=7
      GO TO 375
  360 SMAX=SQRT(SMAX)
      DO 365 J=1,3
  365 V(J)=B(IMAX,J)/SMAX
C     ***** REFINE EIGENVECTOR *****
      CALL AXEQB(A,V,V,1)
      TEM=AMAX1(ABS(V(1)),ABS(V(2)),ABS(V(3)))
      DO 370 J=1,3
  370 V(J)=V(J)/TEM
      CALL UNIT(V,VECT(1,II),1)
C     ***** REFINE EIGENVALUE *****
      T1=VMV(VECT(1,II),W,VECT(1,II))
      U(II)=T1
  375 VALU(II)=T1
  400 RETURN
      END

      SUBROUTINE ERPNT(T1,N)
      INCLUDE 'SIZES'
      CHARACTER*80 COMMNT
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)
      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT
C *     END OF COMMON   *
      IF ( NG .EQ. 5 ) THEN
         WRITE ( COMMNT, '('' ERRONEOUS ATOM NUMBER WAS FOUND '',
     .        I6,'' IN ROUTINE '',I6,''!'')' )
     .        INT( T1 / 100000.), N
      ELSEIF ( NG .EQ. 10 ) THEN
         WRITE ( COMMNT, '('' ATOM '',I6,
     .       '' IS BEYOND THE BORDER!'')' ) INT( T1/100000.)
      ELSEIF ( NG .EQ. 12 ) THEN
         WRITE ( COMMNT, '('' TOO FEW ATOMS IN THE LIST '',
     .        I6, '' IN ROUTINE'', I6, ''!'')' )
     .        INT( T1/100000.), N
      ELSEIF ( NG .EQ. 13 ) THEN
         WRITE ( COMMNT, '('' THE BOND IS TOTALLY BLOCKED BY ATOM '',
     .        I6, '' IN ROUTINE'', I6, ''!'')' )
     .        INT( T1/100000.), N
      ELSEIF ( NG .EQ. 14 ) THEN
        IF ( N.EQ.5)
     .      WRITE ( COMMNT, '('' A BOND WAS HIDDEN (END-ON) BY ATOM '',
     .        I6,'' IN ROUTINE '',I6,''!'')' )
     .        INT( T1 / 100000.), N
      ELSE
         WRITE ( COMMNT, '('' FAULT '', I3, '' OCCURRED WITH ATOM '',
     .   I6, '' IN ROUTINE NUMBER '',I6,''!'')' )
     .                  NG, INT( T1/100000.), N
      ENDIF
      ITEMP = INDEX( COMMNT, '!')
      IF ( ITEMP.GT.0) CALL DEBUGR( COMMNT(1: ITEMP) )
      NG=0
      RETURN
      END

      SUBROUTINE F1000
      RETURN
      END
      SUBROUTINE F200
      DIMENSION HEADER(4),PLA(1500)
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)
      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT
C *     END OF COMMON   *

      DATA ZERO /0.0/

C     ***** CALCOMP CONTROL *****
      NJQ = NJ2+1
      GO TO(260,200,225,205,210,260),NJQ
C     ***** INITIATE MECHANICAL CALCOMP PLOTTER (201 INSTRUCTION) *****
  200 LTNO = 15
      CALL PLOTSZ(PLA,1500,LTNO)
      GO TO 260
C     ***** INITIATE C.R.T. CALCOMP PLOTTER (203 INSTRUCTION) *****
  205 LTNO = -99
C     CALL JOBNUM(HEADER(1))
C     CALL IDAY(HEADER(2))
C     ***** SET C.R.T. X,Y ORIGIN AND BEAM INTENSITY BASELINE *****
      IZ = AIN(3)
      GO TO 215
C     ***** RESET BEAM INTENSITY (204 INSTRUCTION) *****
  210 IZ = AIN(1)
      IF(IZ)215,215,220
  215 IZ = 18
      GO TO 260
220   CONTINUE
C     ***** TERMINATE THE PLOT AND ADVANCE PLOTTER (202 INSTRUCTION) ***
  225 IF(LTNO)235,260,230
  230 CALL OPLOT(AIN(1),AIN(2),8)
      GO TO 260
C     ***** THE FOLLOWING ARE FOR CRT FILM ADVANCE ETC. *****
  235 IF(AIN(1))245,240,250
C     ***** ADD BLOCK ADDRESS BUT DO NOT ADVANCE CRT FILM *****
  240 IB = 1000
      GO TO 255
C     ***** ADVANCE CRT FILM AND CLOSE OUTPUT TAPE FILE ****
  245 IB = 9000
      GO TO 255
C     ***** ADD BLOCK ADDRESS AND ADVANCE CRT FILM *****
  250 IB = 0000
C     ***** IDENTIFICATION TITLE FRAME BETWEEN PLOTS *****
255   CALL OSIM(0.2,0.0,0.25,TITLE(1),40.,72)
  260 RETURN
      END
      SUBROUTINE F400
C     ***** ATOM LIST FUNCTIONS *****
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)
      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT
C *     END OF COMMON   *

      NG=0
      IF(LATM)402,402,400
  400 DO 401 I=1,LATM
  401 CALL ATOM(ATOMS(1,I),ATOMS(2,I))
  402 IF(MOD(NJ2,10)-1)499,404,403
  403 CALL SEARC
      GO TO 499
C     ***** STORES (401) OR REMOVES (411) RUNS OF ATOMS *****
C     ***** RUN HIERARCHY = ATOM NO./SYM/ A/B/C TRANS. *****
  404 II=1
C     ***** FIND RUNS IN AIN ARRAY *****
  405 T1=AIN(II)
      IF(T1)410,410,420
  410 II=II+1
      IF(140-II)499,405,405
  420 JJ=II
C     ***** SET INITIAL RUN VALUES *****
      M1=T1/100000.
      M2=AMOD(T1,100.)
      M5=AMOD(T1/100.,1000.)
      IF(M5)422,422,423
  422 M5=555
  423 M3=M5/100
      M4=MOD(M5/10,10)
      M5=MOD(M5,10)
  425 JJ=JJ+1
      IF(140-JJ)435,430,430
  430 T2=-AIN(JJ)
      IF(T2)435,425,440
  435 II=JJ-1
C     ***** SET TERMINAL VALUES FOR DEGENERATE RUN *****
      N1=M1
      N2=M2
      N3=M3
      N4=M4
      N5=M5
      GO TO 450
  440 II=JJ
C     ***** SET TERMINAL RUN VALUES *****
      N1=T2/100000.
      N2=AMOD(T2,100.)
      N5=AMOD(T2/100.,1000.)
      IF(N5)445,445,446
  445 N5=555
  446 N3=N5/100
      N4=MOD(N5/10,10)
      N5=MOD(N5,10)
C     ***** LOOP THROUGH ALL RUNS *****
  450 DO 490 L5=M5,N5
      DO 490 L4=M4,N4
      DO 490 L3=M3,N3
      DO 490 L2=M2,N2
      DO 490 L1=M1,N1
      V1(1)=(L2)+100.*(L5+10.*(L4+10.*(L3+10.*L1)))
      CALL ATOM(V1(1),V1(2))
      IF(NG)455,458,455
  455 CALL ERPNT(V1(1),401)
      GO TO 490
  458 CALL STOR
  490 CONTINUE
      GO TO 410
  499 RETURN
      END
      SUBROUTINE F500
      DIMENSION RM(3,3),V(3,4)

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      NG=0
      IF(NJ2-11)500,600,610
  500 IF(NJ2-1)610,501,510
  501 T1=AIN(1)
      CALL ATOM(T1,ORGN)
      IF(NG)502,504,502
  502 CALL ERPNT(T1,501)
C?      CALL EXIT
      STOP ' '
  504 DO 506 K=1,4
      T1=AIN(K+1)
      CALL ATOM(T1,V(1,K))
      IF(NG)502,506,502
  506 CONTINUE
      DO 507 J=1,3
      V1(J)=V(J,2)-V(J,1)
  507 V2(J)=V(J,4)-V(J,3)
      IND=-1
      IF(AIN(7))509,509,508
  508 IND=-2
  509 CALL AXES(V1,V2,REFV,IND)
      GO TO 562
  510 IF(NJ2-4)515,511,599
C      ***** SHIFT ORIGIN FOR PROJECTION AXIS (IN INCHES) *****
  511 DO 513 J=1,3
      DO 512 K=1,3
  512 ORGN(J)=ORGN(J)+REFV(J,K)*AIN(K)/SCAL1
  513 XO(J)=XO(J)+AIN(J)
      GO TO 570
C      ***** FORM ROTATION MATRIX *****
  515 DO 552 L=1,139,2
      I=AIN(L)
      IF(I) 532,552,516
  516 X=AIN(L+1)*1.7453293E-2
      T1=COS(X)
      T2=SIN(X)
      I3=MOD(I+2,3)+1
      I1=MOD(I3,3)+1
      I2=MOD(I1,3)+1
      RM(I1,I1)=T1
      RM(I1,I2)=T2
      RM(I1,I3)=0.0
      RM(I2,I1)=-T2
      RM(I2,I2)=T1
      RM(I2,I3)=0.0
      RM(I3,I1)=0.0
      RM(I3,I2)=0.0
      RM(I3,I3)=1.0
      CALL MM(REFV,RM,V)
      IF(NJ2-3)518,525,599
  518 DO 522 J=1,9
  522 REFV(J,1)=V(J,1)
      GO TO 552
  525 DO 528 J=1,9
  528 WRKV(J,1)=V(J,1)
      GO TO 552
  532 IF(NJ2-3)535,552,599
  535 I=MOD(-I,3)
      DO 542 J=1,I
      DO 542 K=1,3
      T1=REFV(K,3)
      REFV(K,3)=REFV(K,2)
      REFV(K,2)=REFV(K,1)
  542 REFV(K,1)=T1
  552 CONTINUE
      IF(NJ2-3)562,582,599
  562 CALL MM(AA,REFV,AAREV)
  570 DO 572 J=1,9
      WRKV(J,1)=REFV(J,1)
  572 AAWRK(J,1)=AAREV(J,1)
      GO TO 599
  582 CALL MM(AA,WRKV,AAWRK)
  599 CONTINUE
C     ***** ELIMINATE ALL PREVIOUSLY STORED OVERLAP INFORMATION ****
C     ***** (ALL INSTRUCTIONS FROM 501 THROUGH 510 DO THIS *****
      CALL LAP500(0)
      GO TO 610
C     ***** STORE NEW OVERLAP INFORMATION (INSTRUCTION 511) *****
  600 CALL LAP500(1)
  610 RETURN
      END
      SUBROUTINE F600
C     ***** SCALING AND CENTERING FUNCTIONS *****
      DIMENSION MAX(3),SCAL(4),X(3),XMAX(3),XMIN(3),Z(2)

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

C     ***** DEL = 1. FOR INCRUMENTING FUNCTIONS *****
C     ***** DEL = 0. FOR REGULAR FUNCTIONS *****
      DEL=(MOD(NJ2/10,2))
      NJ2=MOD(NJ2,10)
C     ***** EXPLICIT ORIGIN AND SCALE *****
      IF(AIN(1))602,604,602
  602 XO(1)=AIN(1)+XO(1)*DEL
  604 IF(AIN(2))606,608,606
  606 XO(2)=AIN(2)+XO(2)*DEL
  608 IF(AIN(3))612,612,609
  609 IF(DEL)611,611,610
  610 SCAL1=SCAL1*AIN(3)
      GO TO 612
  611 SCAL1=AIN(3)
  612 IF(AIN(4))616,616,614
C     ***** SET ELLIPSOID SCALE FACTOR *****
  614 SCAL2=AIN(4)
C     ***** AUTOMATIC ORIGIN AND/OR SCALE *****
  616 IF(NJ2-2)790,622,620
  620 XO(1)=XLNG(1)*.5
      XO(2)=XLNG(2)*.5
  622 IF(NJ2-3)625,640,625
  625 SCAL1=1.
  630 IF(LATM-1)635,635,640
  635 NG=12
      CALL ERPNT(0.,602)
C?      CALL EXIT
      STOP ' '
  640 DO 650 J=1,3
      XMAX(J)=-1.E5
  650 XMIN(J)=1.E5
C     ***** FIT BOX AROUND SET OF ATOMS *****
      DO 670 I=1,LATM
      CALL XYZ(ATOMS(1,I),ATOMS(2,I),3)
      IF(NG)652,653,652
  652 CALL ERPNT(ATOMS(1,I),600)
      GO TO 670
  653 DO 668 J=1,3
      T1=ATOMS(J+1,I)
      IF(XMAX(J)-T1)655,660,660
  655 XMAX(J)=T1
      MAX(J)=I
  660 IF(T1-XMIN(J))665,668,668
  665 XMIN(J)=T1
  668 CONTINUE
  670 CONTINUE
C     ***** KM=TOP ATOM *****
      KM=MAX(3)
      SMULT=1.
      DO 780 M=1,5
      IF(M-2)740,675,678
C     ***** CHECK VIEW DISTANCE *****
  675 IF(VIEW)785,785,680
  678 IF(NJ2-3)680,785,680
  680 T1=ATOMS(4,KM)*SMULT
      IF(VIEW*.5-T1)685,690,690
C     ***** INCREASE VIEW DISTANCE *****
  685 VIEW=2.*T1
C     ***** FIND PERSPECTIVE PROJECTION LIMITS *****
  690 DO 700 J=1,2
      XMAX(J)=-1.E5
  700 XMIN(J)=1.E5
      DO 725 I=1,LATM
      DO 705 J=1,3
  705 X(J)=ATOMS(J+1,I)*SMULT
      T2=VIEW/(VIEW-X(3))
      DO 725 J=1,2
      T1=X(J)*T2
      IF(XMAX(J)-T1)710,715,715
  710 XMAX(J)=T1
  715 IF(T1-XMIN(J))720,725,725
  720 XMIN(J)=T1
  725 CONTINUE
C     ***** REFINE PARAMETERS *****
  740 IF(NJ2 -3)745,742,755
  742 SMUL2=1.
      GO TO 765
C     ***** AUTOMATIC SCALE ONLY *****
  745 DO 750 J=1,2
      T2=XO(J)
      SCAL(J)=(BRDR-T2)/XMIN(J)
  750 SCAL(J+2)=(XLNG(J)-BRDR-T2)/XMAX(J)
      SMUL2=AMIN1(SCAL(1),SCAL(2),SCAL(3),SCAL(4))
      GO TO 780
C     ***** AUTOMATIC SCALE AND POSITION *****
  755 DO 760 J=1,2
  760 SCAL(J)=(XLNG(J)-BRDR*2.)/(XMAX(J)-XMIN(J))
      SMUL2=AMIN1(SCAL(1),SCAL(2))
C     ***** AUTOMATIC POSITION *****
  765 DO 770 J=1,2
  770 XO(J)=.5*(XLNG(J)-SMUL2*(XMAX(J)+XMIN(J)))
  780 SMULT=SMULT*SMUL2
      VIEW=VIEW*SMUL2
  785 SCAL1=SCAL1*SMULT
  790 SCL=SCAL1*SCAL2
C     WRITE (*,*) ' IN F600: SCAL1 =',SCAL1,', VIEW=',VIEW,', SCL=',SCL
C     ***** ELIMINATE ALL PREVIOUSLY STORED OVERLAP INFORMATION *****
      CALL LAP500(0)
      RETURN
      END

      SUBROUTINE F700
C     ***** SUBROUTINE TO DRAW ELLIPSOIDS *****
      DIMENSION EYE(3),VIEWV(3),X(3),Z(3), IITEMP( 6)
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *
      CHARACTER HEADER*80, CHEM(NUMATM)*6, COMMNT*80
      CHARACTER*6 CBCD
      COMMON /STRING/ HEADER, CHEM
      COMMON /AINID/ IDAIN( NUMATM)
C     ***** SET ELLIPSOID GRAPHIC DETAILS *****
      ITILT=0
      NG=0
      NFIRST=1
      NPLANE=AIN(1)
      IF(NPLANE-1)720,715,720
  715 NFIRST=4
      NPLANE=4
  720 NSOLID=AIN(2)
      NDOT=64/2**( IABS(NSOLID))
      LINES=AIN(3)
      NDASH=AIN(4)
      CHSYM=AIN(5)
      DH=AIN(6)-CHSYM*17./7.
      DV=AIN(7)-CHSYM*.5
C     ***** ESTABLISH REFERENCE POINT OF VIEW *****
      T1=1.E6
      IF(VIEW)740,740,735
  735 T1=VIEW/SCAL1
  740 DO 741 J=1,3
  741 EYE(J)=REFV(J,3)*T1+ORGN(J)
       LNS=-1
C     ***** LOOP THROUGH ATOM LIST *****
      DO 1105 ITOM=1,LATM
         T1=ATOMS(1,ITOM)
         IF(AIN(10))744,744,742
  742    T2=AINT(T1/100000.)
         IF(T2-AIN(10))1105,743,743
  743    IF(AIN(11)-T2)1105,744,744
  744    CALL XYZ(T1,X,2)
         IF(NG)758,746,758
  746    CALL PLTXY(X,Z)
         K=T1/100000.
         ICOLOR = IDAIN( K)
         L=AMOD(T1,100000.)/1000.
         M=AMOD(T1,1000.)
         IF(NJ2-10)747,754,754
  747    LNS=MOD(LNS+1,18)
         IF(LNS)749,748,749
  748    WRITE ( COMMNT,751)(TITLE(I),I=1,18)
         CALL DEBUGR( COMMNT(1:78) )
         WRITE ( COMMNT,752)
         CALL DEBUGR( COMMNT(1:78) )
  749    WRITE ( COMMNT,750) CHEM(K),K,L,M,Z(1),Z(2),
     .                       (X(I),I=1,3),(XT(I),I=1,3)
         CALL DEBUGR( COMMNT(1:78) )
  750    FORMAT ( ' ', 10X, A6, '  (', I3, ',', I2, I3, ')   ',
     .            2F8.2, 3X, 3F8.3, 11X, 3F8.4)
  751    FORMAT ( '1', 10X, 18A4)
  752 FORMAT(1H010X,18HSYMBOL   ATOM CODE7X,16HPLOTTER X,Y(IN.) 3X,21HCA
     1RTESIAN X,Y,Z (IN.)15X,20HCRYSTAL SYSTEM X,Y,Z/1H 19X,45H(DIRECTIO
     2N COSINES(I,J),I=1,3),RMSD(J)),J=1,312X,42HFOR PRINCIPAL AXES BASE
     3D ON WORKING SYSTEM/1H )
  754    IF(EDGE-BRDR*.75)755,760,760
  755    NG=10
  758    CALL ERPNT(T1,700)
         GO TO 1105
C     ***** CALL OVERLAP ROUTINE *****
  760    ICQ=0
         CALL LAP700(ITOM,ICQ)
         IF(ICQ)762,764,764
C     ***** OMIT HIDDEN ATOM *****
  762    NG=14
         GO TO 758
  764    IF(CHSYM)775,775,765
C     ***** PLOT CHEMICAL SYMBOLS *****
  765    CONTINUE
*  MY ADDDITION TO ELIMINATE LABELING ATOMS IF OCCLUDED
C? DOES NOT HANDLE THREE BONDS ?C         IF ( ICQ.GT.2 ) GOTO 775
*******
         T4=1.
         IF ( VIEW .LE. 0 ) GOTO 767
  766    T4=VIEW/(VIEW-X(3))
  767    T3=CHSYM*T4
         T4=DISP*T4*.5
         V1(1)=X(1)+DH*SYMB(1,1)+DV*SYMB(1,2)
         V1(2)=X(2)+DH*SYMB(2,1)+DV*SYMB(2,2)
         V1(3)=X(3)
         CALL PLTXY(V1,V3)
         IF ( EDGE .LT. CHSYM) GOTO 775
         IF(NPLANE)1105,1105,768
  768    V2(3)=0.
         ITEMP = IDAIN( K)
         CALL OPLOT(FLOAT(ITEMP), 0.0, 99 )
         DO 770 I=1,3,2
            V2(1)=V3(1)+(I-2)*T4
            DO 770 J=1,3,2
            V2(2)=V3(2)+(J-2)*T4
            NNCHEM = INDEX( CHEM( K), ' ') - 1
C?               DO 769 LTEMP = 1, NNCHEM
C?                  IITEMP( LTEMP) = ICHAR( CHEM( K)( LTEMP: LTEMP) )
C?  769          CONTINUE
            CBCD = CHEM( K)
               CALL OSIM( V2(1), V2(2), T3, CBCD, THETA, NNCHEM)
               IF(T4)775,775,770
  770    CONTINUE
C???  **  MY MOD **  775 IF(NPLANE)1105,1105,780
  775    CONTINUE
C   ***** ELLIPSOID PRINC VECTORS TOWARD VIEWER *****
  780    CALL PAXES(T1,2)
         IF(NG)758,783,758
  783    CALL DIFV(EYE,XT,VIEWV)
         CALL UNIT(VIEWV,VIEWV,-1)
         CALL VM(VIEWV,AA,V2)
         DO 795 I=1,3
            IF(VV(V2,PAT(1,I)))785,795,795
  785       DO 790 J=1,3
               PAC(J,I)=-PAC(J,I)
  790       PAT(J,I)=-PAT(J,I)
  795    CONTINUE
         DO 800 J=1,3
            PAC(J,4)=PAC(J,1)
  800    PAC(J,5)=PAC(J,2)
         IF(NJ2-10)802,803,803
  801    FORMAT(1H 13X,3(3X,3F8.4,F8.5)/1H )
  802    WRITE ( COMMNT,801) ((PAC(J,K),J=1,3),RMS(K) ,K=1,3)
         CALL DEBUGR( COMMNT(1:78) )
C     ***** V4 = VECTOR NORMAL TO POLAR PLANE *****
  803    CALL VM(VIEWV,AAWRK,V6)
         CALL UNIT(V6,V6,1)
         CALL MV(Q,V6,V4)
         CALL UNIT(V4,V4,1)
C     ***** SET PLOTTING RESOLUTION FOR ELLIPSOID *****
         T3=RMS(3)*SCL
         NRESOL=1
         NBIS=5
         DO 805 J=1,3
            IF(T3-RES(J))804,810,810
  804       NBIS=NBIS-1
  805    NRESOL=NRESOL*2
  810    NRES1=NRESOL+1
C     ***** LOOP THROUGH PRINC AND POLAR PLANES *****
         DO 1100 II=NFIRST,NPLANE
            II0=MOD(II+2,3)+1
            II1=MOD(II,3)+1
            II2=MOD(II+1,3)+1
C   ***** GENERATE CONJUGATE DIAMETERS *****
            IF(.99938- ABS(VV(V4,PAC(1,II2))))820,820,830
  820       T1=RMS(II0)*SCL
            T2=RMS(II1)*SCL
            DO 825 J=1,3
               DA(J,1)=PAC(J,II0)*T1
  825       DA(J,2)=PAC(J,II1)*T2
            GO TO 850
  830       CALL NORM(PAC(1,II0),PAC(1,II1),V1,1)
            CALL NORM(V1,V4,V2,1)
            CALL UNIT(V2,V2,1)
            CALL MV(Q,V2,V3)
            IF(II-4)835,840,840
  835       CALL NORM(V3,V1,V5,1)
            GO TO 843
  840       CALL NORM(V3,V4,V5,1)
  843       CALL UNIT(V5,V5,1)
            T1=SCL/SQRT(VMV(V2,Q,V2))
            T2=SCL/SQRT(VMV(V5,Q,V5))
            DO 845 J=1,3
               DA(J,1)=V2(J)*T1
  845       DA(J,2)=V5(J)*T2
C   ***** GENERATE ELLIPSE *****
  850       CALL RADIAL(NBIS)
            IF(II-4)900,851,851
  851       IF(NSOLID)859,859,852
C   ***** PLOT DOTTED BOUNDARY ELLIPSE *****
  852       IF(NDOT-NRESOL)853,855,855
  853       CALL RADIAL(NSOLID-1)
  855       CALL PROJ(D,DP,X,XO,VIEW,1,129,NDOT)
            DO 857 J=1,129,NDOT
               CALL DRAW(DP(1,J),DISP,DISP,3, ICOLOR)
               DO 856 I=1,3,2
                  T1=(I-2)*DISP
                  DO 856 K=1,3,2
                     T2=(K-2)*DISP
                     CALL DRAW(DP(1,J),T1,T2,2, ICOLOR)
                     IF(DISP)857,857,856
  856          CONTINUE
  857       CONTINUE
            GO TO 1100
C   ***** PLOT SOLID BOUNDARY ELLIPSE *****
  859       CALL PROJ(D,DP,X,XO,VIEW,1,129,NRESOL)
            CALL DRAW(DP,0.,0.,3, ICOLOR)
            DO 860 J=NRES1,129,NRESOL
  860       CALL DRAW(DP(1,J),0.,0.,2, ICOLOR)
            IF(DISP)1100,1100,865
C   ***** BOUNDARY ANNULUS AS A LINEAR FUNCTION OF HEIGHT *****
  865       CALL DIFV(XT,ORGN,V1)
            T5=VV(V1,AAREV(1,3))*SCAL1
            NCYCLE=.5+(AIN(8)+T5*AIN(9))/DISP
            IF(NCYCLE)1100,1100,870
  870       T3=(2.*DISP)/(T1+T2)
C   ***** INCREASE ANNULAR THICKNESS *****
            DO 875 I=1,NCYCLE
               T4=T3*(I)
               DO 875 J=1,129,NRESOL
  875       CALL DRAW(DP(1,J),D(1,J)*T4,D(2,J)*T4,2, ICOLOR)
            GO TO 1100
  900       CALL PROJ(D,DP,X,XO,VIEW,1,65,NRESOL)
C   ***** PLOT HALF AN ELLIPSE *****
            CALL DRAW(DP,0.,0.,3, ICOLOR )
            DO 905 J=NRES1,65,NRESOL
  905       CALL DRAW(DP(1,J),0.,0.,2, ICOLOR)
            IF(DISP)930,930,910
C   ***** ACCENTUATE FRONT HALF *****
  910       DO 925 I=1,3,2
               T2=(I-2)*DISP
               DO 915 J=1,65,NRESOL
                  K=66-J
  915          CALL DRAW(DP(1,K),DISP,T2,2, ICOLOR)
               DO 925 K=1,65,NRESOL
  925       CALL DRAW(DP(1,K),-DISP,-T2,2, ICOLOR)
  930       IF(NSOLID)940,967,935
  935       L=NDOT
            IF(NDOT-NRESOL)938,945,940
  938       CALL RADIAL(NSOLID-1)
            GO TO 945
  940       L=NRESOL
  945       CALL PROJ(D(1,65),DP(1,65),X,XO,VIEW,1,65,L)
            IF(NSOLID)960,967,950
C   ***** DOTTED LINE ON REVERSE SIDE *****
  950       DO 958 J=65,129,NDOT
               CALL DRAW(DP(1,J),DISP,DISP,3, ICOLOR)
               DO 955 I=1,3,2
                  T1=(I-2)*DISP
                  DO 955 K=1,3,2
                     T2=(K-2)*DISP
                     CALL DRAW(DP(1,J),T1,T2,2, ICOLOR)
                     IF(DISP)958,958,955
  955          CONTINUE
  958       CONTINUE
            GO TO 967
C   ***** SINGLE LINE ON REVERSE SIDE *****
  960       DO 965 J=65,129,NRESOL
  965       CALL DRAW(DP(1,J),0.,0.,2, ICOLOR)
C   ***** DETAIL INTERIOR FEATURES *****
  967       T2=NDASH*2
            DO 975 J=1,3
               T1=PAC(J,II0)*RMS(II0)*SCL
               DA(J,1)=T1
               DA(J,2)=PAC(J,II1)*RMS(II1)*SCL
               DA(J,3)=0.
               IF(NDASH)975,975,970
  970          V1(J)=-T1
               V2(J)=T1/T2
  975       CONTINUE
            IF(NDASH)987,987,980
C   ***** DASHED LINE FOR REVERSE AXIS *****
  980       DO 985 J=1,NDASH
            DO 985 K=1,2
               L=4-K
               CALL PROJ(V1,DP,X,XO,VIEW,1,1,1)
               CALL DRAW(DP,0.,0.,L, ICOLOR)
               DO 985 I=1,3
  985       V1(I)=V1(I)+V2(I)
C   ***** SOLID LINE FOR FORWARD AXIS *****
  987       IF(LINES)1100,1100,988
  988       CALL PROJ(DA,DP,X,XO,VIEW,1,3,1)
            T1=DISP*.5
            DO 990 I=1,3,2
               T2=(2-I)*T1
               CALL DRAW(DP,T1,T2,3, ICOLOR)
               CALL DRAW(DP(1,3),T1,T2,2, ICOLOR)
               IF(DISP)1000,1000,989
  989          CALL DRAW(DP(1,3),-T1,T2,2, ICOLOR)
  990       CALL DRAW(DP,-T1,T2,2, ICOLOR)
C   ***** SHADE QUADRANT BETWEEN TWO PRINCIPAL AXES *****
 1000       L=LINES-1
            IF(L)1100,1100,1005
 1005       T2=LINES
            DO 1025 I=1,L
               T1=(I)/T2
               T3=SQRT(1.-T1*T1)
               IF(MOD(I,2))1010,1015,1010
 1010          M=I*2
               N=M-1
               GO TO 1020
 1015          N=I*2
               M=N-1
 1020          DO 1025 J=1,3
                  T4=DA(J,1)*T1
                  D(J,M)=T4
 1025          D(J,N)=DA(J,2)*T3+T4
               L=L*2
               CALL PROJ(D,DP,X,XO,VIEW,1,L,1)
               DO 1030 I=2,L,2
                  CALL DRAW(DP(1,I-1),0.,0.,3, ICOLOR)
 1030          CALL DRAW(DP(1,I),0.,0.,2, ICOLOR)
 1100       CONTINUE
 1105 CONTINUE
C     ***** ELIMINATE LOCAL OVERLAP INFORMATION BEFORE RETURNING *****
      CALL LAP500(-1)
      RETURN
      END

      SUBROUTINE F800
C     ***** SUBROUTINE FINDS ATOM PAIRS FOR BONDS *****
      DIMENSION IA(3),W1(6)
      INCLUDE 'SIZES'
      CHARACTER*80 COMMNT
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)
      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT
C *     END OF COMMON   *

      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      LNS=-4
      IF(MOD(NJ2,10)-2)805,848,848
C     ***** EXPLICIT DESCRIPTION *****
  805 II=0
      IF(NCD)810,810,815
  810 NG=11
      CALL ERPNT(0.,NJ*100+NJ2)
      GO TO 980
  815 II=II+1
      IF(140-II)980,980,820
  820 T1=AIN(II)
      IF(T1)815,815,825
  825 II=II+1
      T2=AIN(II)
      IF(T2)815,815,830
  830 IF(NJ2-10)832,838,838
  832 LNS=MOD(LNS+4,56)
      IF(LNS)838,834,838
  834 IF (DEBUGO) THEN
         WRITE ( COMMNT,835)(TITLE(I),I=1,18)
         CALL DEBUGR( COMMNT(1:40) )
         WRITE ( COMMNT,837)
         CALL DEBUGR( COMMNT(1:78) )
      ENDIF
835   FORMAT ( '1', 10X, 18A4)
837   FORMAT ( '0', 10X, 'SYMBOL   ATOM CODE', 6X,
     .         'PLOTTER X,Y(IN.)', 6X,
     .         'CARTESIAN X,Y,Z (IN.)', 17X,
     .         'CRYSTAL SYSTEM X,Y,Z' )
C     ***** CHECK IF BOND ATOMS ARE IN ATOMS LIST (FOR OVERLAP CALC) ***
  838 NA1=0
      NA2=0
      IF(LATM-2)845,839,839
  839 N2=2
      DO 844 K=1,LATM
        T3=ATOMS(1,K)
        IF(T3-T1)841,840,841
  840   NA1=K
        GO TO 843
  841   IF(T3-T2)844,842,844
  842   NA2=K
  843   N2=N2-1
        IF(N2)845,845,844
  844 CONTINUE
  845 IF(NA2-NA1)846,847,847
  846 NA3=NA1
      NA1=NA2
      NA2=NA3
      T3=T1
      T1=T2
      T2=T3
  847 CALL BOND(T1,T2,1,NA1,NA2)
      GO TO 815
C     ***** IMPLICIT DESCRIPTION *****
  848 IF(LATM-2)810,850,850
  850 SCAL3=SCAL1
      SCAL1=1.
      DO 855 I=1,LATM
  855 CALL XYZ(ATOMS(1,I),ATOMS(2,I),2)
      SCAL1=SCAL3
      IF(NCD)810,810,860
  860 IF (DEBUGO) THEN
        WRITE ( COMMNT,861)
        CALL DEBUGR( COMMNT(1:78) )
      ENDIF
  861 FORMAT(1H010X,20HBOND SELECTION CODES//11X,94H(SEQUENCE(A))(SEQUEN
     1CE(B)) (BOND) (DISTANCES)( BOND )(PERSP.--LABELS)(NORMAL--LABELS)(
     2DIGITS) /11X,93H(  MIN  MAX )(  MIN  MAX ) (TYPE) (MIN   MAX)(RADI
     3US)(HEIGHT  OFFSET)(HEIGHT  OFFSET)(NUMBER))
      DMAX=0.
      DO 870 I=1,NCD
c?        IF(DMAX-CD(2,I))865,870,870
c?  865   DMAX=CD(2,I)
        IF ( DMAX .LT. CD(2,I)) DMAX=CD(2,I)
  870 CONTINUE
      IF (DEBUGO) THEN
        WRITE ( COMMNT,871)(KD(J,I),J=1,5),(CD(J,I),J=1,8)
        CALL DEBUGR( COMMNT(1:78) )
      ENDIF
  871 FORMAT(1H 10XI6,I5,I8,I5,I8,2F6.2,5F8.3,F7.0)
      DMAX=DMAX*DMAX
C     ***** LOOP THROUGH ATOMS ARRAY *****
      DO 977 M=1,LATM
        NA1=M
        T1=ATOMS(1,M)
        IA(1)=T1/100000.
        IA(3)=IA(1)
        W1(1)=ATOMS(2,M)
        W1(2)=ATOMS(3,M)
        W1(3)=ATOMS(4,M)
        L=M+1
        IF(LATM-L)977,872,872
  872   DO 975 N=L,LATM
        NA2=N
        DIST=(ATOMS(2,N)-W1(1))**2
        IF(DMAX-DIST)975,873,873
  873   DIST=DIST+(ATOMS(3,N)-W1(2))**2
        IF(DMAX-DIST)975,874,874
  874   DIST=DIST+(ATOMS(4,N)-W1(3))**2
        IF(DMAX-DIST)975,875,875
  875   DIST=SQRT(DIST)
        T2=ATOMS(1,N)
        IA(2)=T2/100000.
C     ***** SELECT BONDS ACCORDING TO CODES *****
        DO 950 J=1,NCD
        JB=J
        IF(DIST-CD(1,J))   950,880,880
  880   IF(CD(2,J)-DIST)   950,881,881
  881   DO 885 K=1,2
        IF(IA(K)-KD(1,J))  885,882,882
  882   IF(KD(2,J)-IA(K))  885,883,883
  883   IF(IA(K+1)-KD(3,J))885,884,884
  884   IF(KD(4,J)-IA(K+1))885,890,890
  885   CONTINUE
        GO TO 950
C     ***** CHECK FOR POLYHEDRA CODE *****
  890   IF(CD(4,J))900,955,955
  900   W1(4)=ATOMS(2,N)
        W1(5)=ATOMS(3,N)
        W1(6)=ATOMS(4,N)
        KM1=ABS(CD(4,J))
        KM2=ABS(CD(5,J))
        DSQ1=(CD(6,J))**2
        DSQ2=(CD(7,J))**2
C     ***** SEARCH FOR POLYHEDRA CENTER *****
        DO 935 IM=1,LATM
        K3=ATOMS(1,IM)/100000.
        IF(K3-KM1)935,905,905
  905    IF(KM2-K3)935,910,910
C     ***** CHECK POLYHEDRA DISTANCE RANGE *****
  910   DO 930 J1=1,4,3
        DSQ=(ATOMS(2,IM)-W1(J1))**2
        IF(DSQ2-DSQ)935,915,915
  915   DSQ=DSQ+(ATOMS(3,IM)-W1(J1+1))**2
        IF(DSQ2-DSQ)935,920,920
  920   DSQ=DSQ+(ATOMS(4,IM)-W1(J1+2))**2
        IF(DSQ2-DSQ)935,925,925
  925   IF(DSQ-DSQ1)935,930,930
  930   CONTINUE
        GO TO 955
  935   CONTINUE
C     ***** END OF POLYHEDRA CHECK *****
  950   CONTINUE
        GO TO 975
C     ***** PREPARE TO DRAW BOND *****
  955   IF(NJ2-10)960,970,970
  960   LNS=MOD(LNS+4,56)
        IF(LNS)965,965,970
  965   IF (DEBUGO) THEN
          WRITE ( COMMNT,835)(TITLE(I),I=1,18)
          CALL DEBUGR( COMMNT(1:78) )
          WRITE ( COMMNT,837)
          CALL DEBUGR( COMMNT(1:78) )
        ENDIF
  970   CALL BOND(T1,T2,JB,NA1,NA2)
  975   CONTINUE
  977 CONTINUE
C     ***** ELIMINATE LOCAL OVERLAP INFORMATION BEFORE RETURNING *****
  980 IF(NJ2-21)985,990,990
  985 CALL LAP500(-1)
  990 RETURN
      END

      SUBROUTINE F900
      DIMENSION X(3),XW(3,5),Y(3),Z(3)

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *
      CHARACTER HEADER*80, CHEM(NUMATM)*6
      CHARACTER*6 CBCD
      COMMON /AINID/ IDAIN( NUMATM)
      COMMON /STRING/ HEADER, CHEM
      DIMENSION IICAR( 6)
C     ***** LABELING FUNCTION SUBROUTINE *****
      ITILT=0
      NJ3=MOD(NJ2,10)
      TH=THETA
      SINTH=SYMB(2,1)
      COSTH=SYMB(1,1)
      ILAST=1
      IF(AIN(2)-11100.)910,910,905
  905 ILAST=2
  910 DO 925 II=1,ILAST
C     ***** OBTAIN WORKING CARTESIAN COORDINATES *****
        CALL XYZ(AIN(II),XW(1,II),2)
        IF(NG)915,925,915
  912   NG=15
  915   CALL ERPNT(AIN(II),NJ*100+NJ2)
        GO TO 1199
  925 CALL XYZ(AIN(II),XW(1,II+3),3)
      II=1
C     ***** FIND MEAN REFERENCE POINT *****
      DO 930 J=1,3
        T2=XW(J,ILAST)
        T1=XW(J,1)
        XW(J,3)=T2-T1
  930 X(J)=(T2+T1)*.5
C     ***** PERSPECTIVE SCALING FACTOR *****
      SCAL=1.
      IF(VIEW)940,940,935
  935 SCAL=VIEW/(VIEW-X(3))
  940 HGT=SCAL*AIN(5)
      IF(NJ2-3)960,950,945
  945 IF(NJ2-6)950,950,960
C     ***** PROJECTED VECTOR BASELINE *****
  950 CALL PLTXY(XW(1,4),V1)
      CALL PLTXY(XW(1,5),V2)
      T1=V2(1)-V1(1)
      T2=V2(2)-V1(2)
      T3=SQRT(T1*T1+T2*T2)
      IF(T3)912,912,955
  955 COSTH=T1/T3
      SINTH=T2/T3
      TH=ACOS(COSTH)
      IF(SINTH)958,960,960
  958 TH=-TH
  960 IF(NJ2-13)965,985,985
C     ***** FIND CENTER OF PROJECTED LABEL *****
  965 Y(1)=SCAL*(X(1)+AIN(6)*COSTH-AIN(7)*SINTH)+XO(1)
      Y(2)=SCAL*(X(2)+AIN(6)*SINTH+AIN(7)*COSTH)+XO(2)
      Y(3)=0.
C     ***** CHECK FOR LEGEND RESET *****
      DO 980 J=1,2
        T1=AIN(J+2)
        IF(T1)975,980,970
  970   Y(J)=T1
        GO TO 980
  975   Y(J)=XLNG(J)+T1
  980 CONTINUE
C     ***** SET PARAMETERS FOR INDIVIDUAL FUNCTIONS *****
  985 GO TO(990,995,995,1000,1000,1000,915,1105,1105,915,915,915,1005,10
     105,1005,1005,915),NJ2
  990 T6=17.
      L=AIN(1)/100000.
      ILAST=1
      DXW=0.
      DYW=0.
      GO TO 1030
  995 T6=215.
      ILAST=18
      T1=HGT*24./7.
      DXW=COSTH*T1
      DYW=SINTH*T1
      GO TO 1030
 1000 T6=10+3*(NJ3-4)
      DIST=SQRT(VV(XW(1,3),XW(1,3)))/SCAL1
      GO TO 1030
C     ***** TRUE PERSPECTIVE LABELS *****
 1005 CALL UNIT(XW(1,3),VT(1,1),1)
      IF(ABS(VT(3,1))-.9994)1010,912,912
C     ***** FORM PERSPECTIVE ROTATION MATRIX *****
 1010 CALL NORM(AID(1,3),VT(1,1),VT(1,2),1)
      CALL UNIT(VT(1,2),VT(1,2),1)
      CALL NORM(VT(1,1),VT(1,2),VT(1,3),1)
      DO 1015 J=1,3
 1015 VT(J,4)=X(J)
      ITILT=1
      HGT=AIN(5)
      TH=0.
      Y(3)=X(3)
      Y(2)=X(2)+AIN(7)-HGT*.5
      IF(NJ2-13)1030,1025,1020
C     ***** PERSPECTIVE BOND LABELS *****
 1020 Y(1)=X(1)+AIN(6)-HGT*(22+3*(6-NJ3))/7.
      DIST=SQRT(VV(XW(1,3),XW(1,3)))/SCAL1
      GO TO 1050
C     ***** PERSPECTIVE TITLES *****
 1025 Y(1)=X(1)+AIN(6)-HGT*215./7.
      ILAST=18
      DXW=HGT*24./7.
      DYW=0.
      GO TO 1050
 1030 DH=HGT*T6/7.
      DV=HGT*.5
      Y(1)=Y(1)-DH*COSTH+DV*SINTH
      Y(2)=Y(2)-DH*SINTH-DV*COSTH
      Y(3)=0.
C     ***** PLOT VARIOUS LABELS *****
 1050 Z(3)=Y(3)
      XO(3)=Y(3)
      CALL OPLOT( 1.0, 0.0, 99 )
      GO TO(1060,1060,1060,1090,1090,1090,915,1105,1105),NJ3
 1060 DO 1085 I=1,ILAST
        DO 1075 J=1,3,2
        Z(1)=Y(1)+(J-2)*DISP*.5
        DO 1075 K=1,3,2
        Z(2)=Y(2)+(K-2)*DISP*.5
        IF(NJ3-2)1065,1068,1068
C     ***** PLOT CHEMICAL SIMBOL *****
 1065   CONTINUE
        CALL OPLOT(FLOAT(IDAIN( L)), 0.0, 99)
C?      IICAR( 1) = ICHAR( CHEM( L)( 1: 1) )
C?      IICAR( 2) = ICHAR( CHEM( L)( 2: 2) )
C?      IICAR( 3) = ICHAR( CHEM( L)( 3: 3) )
C?      IICAR( 4) = ICHAR( CHEM( L)( 4: 4) )
C?      IICAR( 5) = ICHAR( CHEM( L)( 5: 5) )
C?      IICAR( 6) = ICHAR( CHEM( L)( 6: 6) )
        CBCD = CHEM( L)
C?      CALL OSIM( Z(1), Z(2), HGT, IICAR, TH, 6)
        CALL OPLOT( 1.0, 0.0, 99)
        GO TO 1070
C     ***** PLOT TITLES *****
 1068   CONTINUE
        CALL OPLOT( 1.0, 0.0, 99 )
C?      IICAR( 1) = ICHAR( HEADER( I*4-3: I*4-3) )
C?      IICAR( 2) = ICHAR( HEADER( I*4-2: I*4-2) )
C?      IICAR( 3) = ICHAR( HEADER( I*4-1: I*4-1) )
C?      IICAR( 4) = ICHAR( HEADER( I*4  : I*4  ) )
        CBCD = HEADER( I*4-3 : I*4 )
C?      CALL OSIM(Z(1), Z(2), HGT, IICAR, TH, 4)
        CALL OSIM(Z(1), Z(2), HGT, CBCD, TH, 4)
 1070   IF(DISP)1080,1080,1075
 1075   CONTINUE
 1080   Y(1)=Y(1)+DXW
 1085 Y(2)=Y(2)+DYW
      GO TO 1199
C     ***** PLOT BOND DISTANCE LABELS *****
 1090 I9=NJ3-3
      T9=10.**I9
      DISTR=AINT((DIST*T9)+0.5)/T9 +.0001
      CALL OPLOT( 1.0, 0.0, 99 )
      CALL NUMBUR(Y(1),Y(2),HGT,DISTR,TH,I9)
      GO TO 1199
C     ***** PLOT CENTERED SIMBOLS *****
C? 1105 CALL OSIM(Y(1),Y(2),HGT, IFIX(AIN(8)),TH,7-NJ3)
 1105 CALL OSIM(Y(1),Y(2),HGT, CHAR(IFIX(AIN(8))), TH,7-NJ3)
 1199 ITILT=0
      RETURN
      END

      SUBROUTINE LAP500(NTYPE)
      CHARACTER*80 COMMNT
C     ***** STORE PROJECTED ATOM CONICS AND BOND QUADRANGLES *****
      DIMENSION QC(3,3),QD(3,3),VD1(3),VD2(3)
C     ***** REAL FOR VERY SMALL ATOMS AND SHORT WORDLENGTH *
C     ***** IN GENERAL, REAL IS NOT REQUIRED *****
      INCLUDE 'SIZES'
      COMMON /OLAP/ CONIC(7,NUMATM), COVER(6,20), KC(20), KQ(30),
     . NCONIC,NCOVER,
     1 NQOVER,NQUAD,OVMRGN,QOVER(3,4,30),QUAD(9,600),SEGM(50,2)
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
C     ***** ELIMINATE ALL PREVIOUSLY STORED LOCAL OVERLAP INFORMATION **
      NCOVER=0
      NQOVER=0
      IF(NTYPE)420,195,195
C     ***** ELIMINATE ALL PREVIOUSLY STORED GLOBAL OVERLAP INFORMATION *
  195 NCONIC=0
      NQUAD=0
      IF(NTYPE)420,420,200
C     ***** CONSTANT FOR OVERLAP MARGIN (WHITE MARGIN AT OVERLAP) *****
  200 IF(AIN(1))205,215,210
C     ***** NEGATIVE NUMBER OR POSITIVE INTEGER GIVES OVMRGN=0.0 *****
  205 OVMRGN=0.0
      GO TO 220
C     ***** SET OVERLAP MARGIN WIDTH DIRECTLY IN INCHES *****
  210 OVMRGN=AIN(1)-AINT(AIN(1))
      GO TO 220
C      ***** DEFAULT OPTION, OVERLAP MARGIN WIDTH AS A FUNCTION OF SCAL1
  215 OVMRGN=AMAX1(SQRT(SCAL1)*0.030,0.025)
  220 IF (DEBUGO) THEN
         WRITE ( COMMNT,2) OVMRGN
         CALL DEBUGR( COMMNT(1:78) )
      ENDIF
    2 FORMAT ( '0', 10X, 'OVERLAP MARGIN IS', F6.3, ' INCH')
  225 IF(LATM)230,230,235
  230 NG=12
      CALL ERPNT(0.,510+NJ2)
      GO TO 420
C     ***** SORT ATOMS LIST BY -VIEWDISTANCE OR BY Z PARAMETER
  235 IF(VIEW)250,250,240
C     ***** CALCULATE VIEWDISTANCES**2 *(-1) IF VIEW.GT.ZERO *****
  240 DO 245 I=1,LATM
      CALL XYZ(ATOMS(1,I),V3,2)
      V3(3)=V3(3)-VIEW
  245 ATOMS(4,I)=-VV(V3,V3)
      GO TO 260
C     ***** STORE CARTESIAN COORDINATES IF VIEW.EQ.ZERO *****
  250 DO 255 I=1,LATM
  255 CALL XYZ(ATOMS(1,I),ATOMS(2,I),2)
C     ***** SORTING PROCEDURE BY SHELL, COMM ACM 2,30 (1959) *****
  260 M=LATM
  265 M=M/2
      IF(M)300,300,270
  270 K=LATM-M
      J=1
  275 I=J
  280 IM=I+M
      IF(ATOMS(4,I)-ATOMS(4,IM))295,295,285
  285 DO 290 L=1,4,3
      T1=ATOMS(L,I)
      ATOMS(L,I)=ATOMS(L,IM)
  290 ATOMS(L,IM)=T1
      I=I-M
      IF(I)295,295,280
  295 J=J+1
      IF(J-K)275,275,265
C     ***** LOOP THROUGH ALL ATOMS IN SORTED ATOMS LIST *****
  300 DO 405 IA=1,LATM
      CALL XYZ(ATOMS(1,IA),ATOMS(2,IA),2)
      CALL PAXES(ATOMS(1,IA),2)
      DO 305 J=1,3
      V1(J)=ATOMS(J+1,IA)
      VD1(J)=V1(J)
      DO 305 K=1,3
  305 QD(J,K)=Q(J,K)
      IF(VIEW)340,340,310
C     ***** CALCULATE ENVELOPING CONE WITH ORIGIN AT VIEWPOINT *****
  310 V1(3)=V1(3)-VIEW
      VD1(3)=V1(3)
C     ***** FORM COFACTOR MATRIX *****
      DO 315 J=1,3
      J1=MOD(J,3)+1
      J2=MOD(J+1,3)+1
      DO 315 K=J,3
      K1=MOD(K,3)+1
      K2=MOD(K+1,3)+1
      QC(J,K)= QD(J1,K1)*QD(J2,K2)-QD(J1,K2)*QD(J2,K1)
  315 QC(K,J)=QC(J,K)
C     ***** FORM POLARIZED COFACTOR MATRIX AND ADD TO ELLIPSOID MATRIX *
      TD2=-SCL**2
C     ***** TD1 IS AN ARBITRARY SCALING FACTOR *****
      TD1=VMV(V1,Q,V1)
      DO 325 J=1,3
      J1=MOD(J,3)+1
      J2=MOD(J+1,3)+1
      DO 320 K=J,3
      K1=MOD(K,3)+1
      K2=MOD(K+1,3)+1
      QD(J,K)=((VD1(J2)*(QC(J1,K1)*VD1(K2)-QC(J1,K2)*VD1(K1))
     1 +VD1(J1)*(VD1(K1)*QC(J2,K2)-VD1(K2)*QC(J2,K1)))+TD2*QD(J,K))/TD1
  320 QD(K,J)=QD(J,K)
C     ***** PROJECTED ELLIPSE IN HOMOGENEOUS COORD OF WORKING SYSTEM ***
      QD(J,3)=-QD(J,3)*VIEW
  325 QD(3,J)=-QD(3,J)*VIEW
C     ***** PROJECT CENTER OF ATOM ONTO PROJECTION PLANE *****
      TD1=-VIEW/VD1(3)
      VD2(1)=VD1(1)*TD1
      VD2(2)=VD1(2)*TD1
C     ***** TRANSFORM TO NEW ORIGIN TO IMPROVE CONDITION OF MATRIX Q ***
      DO 330 J=1,3
      DO 330 K=1,2
  330 QD(J,3)=QD(J,3)+QD(J,K)*VD2(K)
      DO 335 J=1,3
      DO 335 K=1,2
  335 QD(3,J)=QD(3,J)+VD2(K)*QD(K,J)
      V6(1)=XO(1)+VD2(1)
      V6(2)=XO(2)+VD2(2)
      GO TO 355
C     ***** CALCULATE ENVELOPING CYLINDER ALONG Z OF WORKING SYSTEM ****
  340 DO 345 J=1,2
      DO 345 K=1,2
  345 QD(J,K)=QD(J,K)-QD(J,3)*QD(K,3)/QD(3,3)
      DO 350 J=1,2
      QD(J,3)=0.0
      QD(3,J)=0.0
  350 V6(J)=XO(J)+ATOMS(J+1,IA)
C     ***** PROJECTED ELLIPSE IN HOMOGENEOUS COORD ABOUT CENTER OF ATOM
      QD(3,3)=-SCL**2
C     ***** FIT RECTANGLE AROUND ELLIPSE ALLOWING OVERLAP MARGIN *****
C     ***** FORM MATRIX OF COFACTORS *****
  355 DO 360 J=1,3
      J1=MOD(J,3)+1
      J2=MOD(J+1,3)+1
      DO 360 K=J,3
      K1=MOD(K,3)+1
      K2=MOD(K+1,3)+1
  360 QC(J,K)= QD(J1,K1)*QD(J2,K2)-QD(J1,K2)*QD(J2,K1)
C     ***** RESCALE MATRIX OF COFACTORS SO THAT QC(3,3)=1.0 *****
      DO 365 J=1,3
      DO 365 K=J,3
      QC(J,K)= QC(J,K)/QC(3,3)
  365 QC(K,J)=QC(J,K)
      TD2=QD(3,3)
      NDG=0
      DO 385 J=1,2
C     ***** SOLVE QUADRATIC EQUATION *****
      T1 =QC(3,J)**2-QC(J,J)
      IF(T1)370,370,375
C     ***** ROUNDOFF PROBLEMS, RESET LIMITS IN X OR Y *****
  370 NDG=1
      V5(J)=0.001+OVMRGN
      GO TO 380
  375 V5(J)=  SQRT(T1)+OVMRGN
      V6(J)=V6(J)+QC(3,J)
      TD2=TD2+QD(3,J)*QC(3,J)
  380 CONIC(2*J-1,IA)=V6(J)-V5(J)
  385 CONIC(2*J,IA)=V6(J)+V5(J)
      IF(NDG)390,390,395
  390 IF(TD2)400,395,395
C     ***** ELLIPSE IMAGINARY DUE TO ROUNDOFF, RESET TO REAL VALUE *****
  395 CONIC(5,IA)=1.0/((CONIC(2,IA)-CONIC(1,IA))*0.5)**2
      CONIC(6,IA)=0.0
      CONIC(7,IA)=1.0/((CONIC(4,IA)-CONIC(3,IA))*0.5)**2
      GO TO 405
C     ***** STORE NORMALIZED QUADRATIC COEFFICIENTS FOR ELLIPSE *****
C     ***** SCALED BY OVERLAP MARGIN PARAMETER *****
  400 TD3= -(1.0-2.0*OVMRGN/(V5(1)+V5(2)))**2 /TD2
      CONIC(5,IA)=QD(1,1)*TD3
      CONIC(6,IA)=QD(1,2)*TD3
      CONIC(7,IA)=QD(2,2)*TD3
  405 CONTINUE
      NCONIC=LATM
C     ***** PRINT OUT SORTED ATOMS ARRAY *****
      IF (DEBUGO) THEN
         WRITE ( COMMNT,4)(ATOMS(1,J),J=1,LATM)
         CALL DEBUGR( COMMNT(1:78) )
      ENDIF
    4 FORMAT(1H0,10X,30HCONTENTS OF SORTED ATOMS ARRAY/(15X,10F10.0))
C     ***** STORE BOND QUADRANGLES IF SEARCH CODES ARE GIVEN *****
      IF(NCD)420,420,410
C     ***** GENERATE PSEUDO-INSTRUCTION 822 TO CALCULATE BONDS *****
  410 NJ2=22
      CALL F800
      IF(NQUAD)420,420,415
C     ***** PRINT OUT NUMBER OF BOND QUADRANGLES STORED *****
C     ***** PRINT OUT QUADRANGLE IDENTIFICATION ARRAY *****
  415 IF (DEBUGO) THEN
         WRITE ( COMMNT,6)NQUAD,(QUAD(9,J),J=1,NQUAD)
         CALL DEBUGR( COMMNT(1:78) )
      ENDIF
    6 FORMAT(1H0,10X,27HBOND OVERLAP ARRAY CONTAINS,I4,23H BONDS (MAXIMU
     1M IS 599)/  11X,  66HATOM-PAIR NUMBERS IN ARRAY REFER TO SEQUENCE
     2IN SORTED ATOMS ARRAY/(15X,10F10.0))
  420 RETURN
      END
      SUBROUTINE LAP700(NA,ICQ)
      DIMENSION DETER(2),QA(3,3,2),QC(3,3,2),V12(3,2),YMIN(2),YMAX(2)
      DIMENSION OVMR(2)
C     ***** IN GENERAL, REAL IS NOT REQUIRED *****
C     ***** BUT IN CASE OF TROUBLE, ACTIVATE THE CD COMMENT CARDS *****
CD    REAL*8 AOV3,AOV3SQ,BOV3,DETER,PI,PHI,POV3,POV3CU,QA,QC,QOV2,QOV2SQ
CD    REAL*8 ROOT,TD
      INCLUDE 'SIZES'
      COMMON /OLAP/ CONIC(7,NUMATM), COVER(6,20), KC(20), KQ(30),
     . NCONIC,NCOVER,
     1 NQOVER,NQUAD,OVMRGN,QOVER(3,4,30),QUAD(9,600),SEGM(50,2)

C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      PI=3.141593
      ICQ=0
      NCOVER=0
      NQOVER=0
      OVMR(1)=OVMRGN
      OVMR(2)=0.0
      IF(NCONIC-NA)200,200,205
  200 RETURN
C     ***** ROUGH CHECK FOR OVERLAPPING ATOMS *****
  205 DO  210 J=1,2
         YMIN(J)=CONIC(2*J-1,NA)
  210 YMAX(J)=CONIC(2*J,NA)
      L=0
      DO 420 IA=NA,NCONIC
         IF(IA-NA)230,230,215
  215    DO 225 J=1,2
            IF(YMAX(J)-CONIC(2*J-1,IA))420,420,220
  220       IF(YMIN(J)-CONIC(2*J,IA))225,420,420
  225    CONTINUE
C     ***** EXACT CHECK FOR OVERLAPPING ATOMS *****
  230    IF(L-1)235,235,240
  235    L=L+1
  240    CALL LAPCON(CONIC(1,IA),DA,V12(1,L),OVMR(L))
         DO 245 J=1,3
         DO 245 K=1,3
  245    QA(J,K,L)=DA(J,K)
C     ***** CALCULATE COFACTORS AND DETERMINANTS *****
         DETER(L)=0.0
         DO 250 J=1,3
            J1=MOD(J+3,3)+1
            J2=MOD(J+1,3)+1
            DO 250 K=1,3
               K1=MOD(K+3,3)+1
               K2=MOD(K+1,3)+1
               TD=QA(J1,K1,L)*QA(J2,K2,L)-QA(J1,K2,L)*QA(J2,K1,L)
               DETER(L)=DETER(L)+TD*QA(J,K,L)
  250    QC(J,K,L)=TD
C     ***** DETER(L) IS THE DETERMINANT TIMES 3 *****
         IF(L-1)420,420,255
C     ***** FORM CHARACTERISTIC EQUATION AND EXAMINE ITS ROOTS *****
  255    AOV3=0.0
         BOV3=0.0
         DO 260 J=1,3
         DO 260 K=1,3
            AOV3=AOV3+QC(J,K,2)*QA(J,K,1)
  260    BOV3=BOV3+QC(J,K,1)*QA(J,K,2)
         AOV3=AOV3/ DETER(2)
         AOV3SQ=AOV3**2
         BOV3=BOV3/ DETER(2)
         POV3=BOV3-AOV3SQ
CD    QOV2=AOV3*(AOV3SQ-BOV3*1.5D0)+DETER(1)/(DETER(2)*2.0D0)
         QOV2=AOV3*(AOV3SQ-BOV3*1.5  )+DETER(1)/(DETER(2)*2.0  )
C     ***** CHECK DISCRIMINANT OF CHARACTERISTIC CUBIC EQUATION *****
         ITYPE=0
         POV3CU=POV3**3
         QOV2SQ=QOV2**2
         IF(POV3CU+QOV2SQ)270,310,265
  265    IF(POV3CU*1.00001 +QOV2SQ)310,310,400
  270    IF(POV3CU+1.00001 *QOV2SQ)275,310,310
C     ***** THREE REAL ROOTS, ALL DIFFERENT *****
  275    ITYPE=1
C     ***** NO INTERSECTION IF A/3 AND B/3 INVARIANTS ARE NEGATIVE *****
         IF(AOV3)280,285,285
  280    IF(BOV3)420,285,285
C     ***** CALCULATE ONE ROOT OF CHARACTERISTIC CUBIC EQUATION *****
  285    IF(QOV2)295,290,295
CD290 PHI=PI/2.0D0
  290    PHI=PI/2.0
         GO TO 305
CD295 PHI=DATAN(-DSQRT(-POV3CU-QOV2SQ)/QOV2)
  295    PHI= ATAN(- SQRT(-POV3CU-QOV2SQ)/QOV2)
         IF(PHI)300,305,305
  300    PHI=PHI+PI
CD305 ROOT=2.0D0*DSQRT(-POV3)*DCOS(PHI/3.0D0)-AOV3
  305    ROOT=2.0  * SQRT(-POV3)* COS(PHI/3.0  )-AOV3
         GO TO 325
C     ***** THREE REAL ROOTS, AT LEAST TWO ARE EQUAL *****
  310    ITYPE=2
C     ***** CHECK SIGNS OF INVARIANTS A/3 AND B/3 *****
         IF(AOV3)315,320,320
  315    IF(BOV3)420,320,320
C     ***** CALCULATE REPEATED ROOT OF CUBIC EQUATION *****
CD320 ROOT=DSIGN(DSQRT(-POV3),QOV2)-AOV3
  320    ROOT= SIGN( SQRT(-POV3),QOV2)-AOV3
C     ***** FORM DEGENERATE CONIC (LINE PAIR WHICH MAY BE COINCIDENT) **
  325    DO 330 J=1,3
         DO 330 K=1,3
  330    DA(J,K)=QA(J,K,1)+ROOT*QA(J,K,2)
C     ***** EXAMINE INVARIANTS OF THE DEGENERATE CONIC *****
         T6=DA(1,1)*DA(2,2)
         T7=DA(1,2)**2
C     ***** NEGATIVE DENOTES REAL INTERSECTING LINE PAIR *****
C     ***** POSITIVE DENOTES IMAGINARY LINES INTERSECTING AT REAL POINT
         IF(T6-T7)335,345,340
  335    IF(T6*1.0001  -T7)400,345,345
  340    IF(T6-1.0001  *T7)345,345,365
  345    T8=DA(3,3)*(DA(1,1)+DA(2,2))
         T9=DA(1,3)**2+DA(2,3)**2
C     ***** NEGATIVE DENOTES REAL PARALLEL LINE PAIR *****
C     ***** POSITIVE DENOTES IMAGINARY PARALLELS *****
C     ***** ZERO DENOTES ONE REAL LINE (COINCIDENT PARALLELS) *****
         IF(T8-T9)350,360,355
  350    IF(T8*1.0001  -T9)400,360,360
  355    IF(T8-1.0001  *T9)360,360,365
C     ***** COINCIDENT LINE PAIR FOUND FOR THE REPEATED ROOT *****
  360    ITYPE=3
C     ***** COMPARE AREAS OF CONICS *****
  365    KA=1
         KB=2
         IF(QC(3,3,KA)-QC(3,3,KB))370,375,375
  370    KA=2
         KB=1
C     ***** SEE IF ONE CONIC IS INSIDE THE OTHER CONIC *****
  375    T1=0.0
         DO 385 J=1,3
            T2=QA(J,3,KB)
            DO 380 K=1,2
  380       T2=T2+QA(J,K,KB)*V12(K,KA)
  385    T1=T1+V12(J,KA)*T2
C     ***** DISCARD IF KA IS OUTSIDE KB *****
         IF(T1)390,390,420
  390    IF(KA-1)395,395,400
C     ***** THE OVERLAPPING ATOM HIDES THE ORIGINAL ATOM *****
  395    ICQ=-1
         RETURN
C     ***** STORE OVERLAPPING ATOM *****
  400    ICQ=ICQ+1
         IF(NCOVER-20)410,405,405
  405    NG=17
         CALL ERPNT(ATOMS(1,IA),700)
         NCOVER=NCOVER-1
  410    NCOVER=NCOVER+1
         IJ=1
         DO 415 I=1,3
         DO 415 J=I,3
            COVER(IJ,NCOVER)=QA(I,J,2)
  415    IJ=IJ+1
         KC(NCOVER)=IA
  420 CONTINUE
C     ***** SECOND PART OF SUBROUTINE CHECKS FOR BONDS OVER THE ATOM ***
  425 IF(NQUAD)470,470,430
  430 ITY=0
C     ***** ROUGH CHECK FOR OVERLAPPING BONDS *****
      DO 465 IQ=1,NQUAD
         TID=QUAD(9,IQ)
         NA1=TID/1000.
         NA2=AMOD(TID,1000.)
         IF(NA-NA2)435,435,465
  435    DO 445 J=1,2
            IF (YMAX(J)- AMIN1(QUAD(J,IQ), QUAD(J+2,IQ),
     .                   QUAD(J+4,IQ),QUAD(J+6,IQ))) 465,465,440
  440       IF (YMIN(J)- AMAX1(QUAD(J,IQ), QUAD(J+2,IQ),
     .                   QUAD(J+4,IQ),QUAD(J+6,IQ))) 445,465,465
  445    CONTINUE
C     ***** EXACT CHECK FOR OVERLAPPING BONDS *****
  450    ITY=ITY-1
         IQQ=0
         CALL LAPAB(IQ,NA,IQQ,ITY)
         IF(IQQ)455,460,460
  455    ICQ=-1
         RETURN
  460    ICQ=ICQ+1
         IF(NQOVER-30)465,470,470
  465 CONTINUE
  470 RETURN
      END
      SUBROUTINE LAP800(NA1,NA2,ICQ)
C     ***** SUBROUTINE CHECKS FOR ATOMS AND BONDS OVERLAPPING A BOND ***
      DIMENSION FL(4,4),Y1(2),Y2(2),YMAX(2),YMIN(2),QUA(3,4)
      DIMENSION VUE(3)
      INCLUDE 'SIZES'
      COMMON /OLAP/ CONIC(7,NUMATM), COVER(6,20), KC(20), KQ(30),
     . NCONIC,NCOVER,
     1 NQOVER,NQUAD,OVMRGN,QOVER(3,4,30),QUAD(9,600),SEGM(50,2)

C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      IQ=0
      ICQ=0
      IF(NA1*NA2)245,245,195
  195 TID1=(NA1)*1000.+(NA2)
      IF(NCONIC)245,245,200
  200 IF(NJ2-21)250,205,205
C     ***** PART 1, CALLED FROM BOND, STORES BOND OUTLINE INFORMATION **
  205 IF(NQUAD-599)215,210,210
  210 NG=16
      CALL ERPNT(ATOMS(1,NA1),822)
      GO TO 245
  215 NQUAD=NQUAD+1
C     ***** CALCULATE OVERLAP MARGIN FOR BOND QUADRANGLE *****
      T1=0.0
      T2=0.0
      DO 220 J=1,2
      Y1(J)=DP(J,1)-DP(J,65)
      Y2(J)=DP(J,2)-DP(J,66)
      T1=T1+Y1(J)**2
  220 T2=T2+Y2(J)**2
      IF(T1*T2)225,225,230
  225 T1=0.0
      T2=0.0
      GO TO 235
  230 T1=OVMRGN/SQRT(T1)
      T2=OVMRGN/SQRT(T2)
C     ***** STORE BOND QUADRANGLE *****
  235 DO 240 J=1,2
      Y1(J)=Y1(J)*T1
      Y2(J)=Y2(J)*T2
      QUAD(J,NQUAD)=DP(J,1)+Y1(J)
      QUAD(J+2 ,NQUAD)=DP(J,2)+Y2(J)
      QUAD(J+4,NQUAD)=DP(J,66)-Y2(J)
  240 QUAD(J+6,NQUAD)=DP(J,65)-Y1(J)
      QUAD(9,NQUAD)=TID1
  245 RETURN
C     ***** PART 2, CALLED FROM BOND, OVERLAP CHECK FOR BOND NA1-NA2 ***
  250 NCOVER=0
      NQOVER=0
      TOL=1.E-5
      IF(NCONIC-NA1)245,245,255
C     ***** SAVE QUADRANGLE TEMPORARILY *****
  255 IQ=NQUAD+1
      DO 260 J=1,2
      QUAD(J,IQ)=DP(J,1)
      QUAD(J+2,IQ)=DP(J,2)
      QUAD(J+4,IQ)=DP(J,66)
  260 QUAD(J+6,IQ)=DP(J,65)
      QUAD(9,IQ)=TID1
C     ***** FIT RECTANGLE AROUND QUADRANGLE *****
  265 DO 270 J=1,2
      YMIN(J)=AMIN1(DP(J,1),DP(J,2),DP(J,66),DP(J,65))
  270 YMAX(J)=AMAX1(DP(J,1),DP(J,2),DP(J,66),DP(J,65))
C     ***** ROUGH CHECK FOR ATOM-OVER-BOND OVERLAP *****
      NA1P1=NA1+1
      ITY=0
      DO 305 IA=NA1P1,NCONIC
      DO 285 J=1,2
      IF(IA-NA2)275,305,275
  275 IF(YMAX(J)-CONIC(2*J-1,IA))305,305,280
  280 IF(YMIN(J)-CONIC(2*J,IA))285,305,305
  285 CONTINUE
C     ***** CHECK FOR TRUE ATOM-OVER-BOND OVERLAP *****
      ITY=ITY+1
      CALL LAPAB(IQ,IA,IQQ,ITY)
      IF(IQQ)290,305,300
C     ***** HIDDEN BOND *****
  290 ICQ=-1
  295 RETURN
  300 ICQ=ICQ+1
      IF(NCOVER-20)305,310,310
  305 CONTINUE
  310 IF(NQUAD)295,295,315
C     ***** ROUGH CHECK FOR BOND-OVER-BOND OVERLAP *****
  315 CALL DIFV(ATOMS(2,NA2),ATOMS(2,NA1),V1)
      CALL UNIT(V1,V1,1)
      VUE(1)=  ATOMS(2,NA1)
      VUE(2)=  ATOMS(3,NA1)
      VUE(3)=  ATOMS(4,NA1)-VIEW
      DO 495 IB=1,NQUAD
      TID2=QUAD(9,IB)
      IF(TID1-TID2)320,495,320
  320 NB2=AMOD(TID2,1000.)
      NB1=TID2/1000.
      IF(NA1-NB2)325,495,495
  325 DO 335 J=1,2
      IF(YMAX(J)-AMIN1(QUAD(J,IB),QUAD(J+2,IB),QUAD(J+4,IB),QUAD(J+6,IB)
     1 ))495,495,330
  330 IF(YMIN(J)-AMAX1(QUAD(J,IB),QUAD(J+2,IB),QUAD(J+4,IB),QUAD(J+6,IB)
     1 ))335,495,495
  335 CONTINUE
C     ***** SET UP LINEAR FORMS FOR EDGES OF QUADRANGLE *****
      DO 345 L=1,4
      K=2*L
      K1=MOD(K,8)+2
      QUA(1,L)=QUAD(K,IB)-QUAD(K1,IB)
      QUA(2,L)=QUAD(K1-1,IB)-QUAD(K-1,IB)
      QUA(3,L)=QUAD(K-1,IB)*QUAD(K1,IB)-QUAD(K,IB)*QUAD(K1-1,IB)
C     ***** NORMALIZE LINE EQUATION COEFFICIENTS *****
      T1=SQRT(QUA(1,L)**2+QUA(2,L)**2)
      IF(T1)495,495,340
  340 DO 345 J=1,3
  345 QUA(J,L)=QUA(J,L)/T1
C     ***** EVALUATE LINEAR FORMS AND SIGNATURES FOR QUADRANGLE *****
      T3=3.0
      DO 365 K=1,4
      T2=3.0
      J=K*2
      DO 355 L=1,4
      T1=QUAD(J-1,IQ)*QUA(1,L)+QUAD(J,IQ)*QUA(2,L)+QUA(3,L)
      IF(T1)350,355,355
  350 T2=T2-1.0
  355 FL(L,K)=T1
      IF(T2)360,365,365
  360 T3=T3-1.0
  365 CONTINUE
C     ***** CHECK FOR 4 POINTS INSIDE QUADRANGLE *****
      IF(T3)370,375,375
  370 ITYPE=-1
      GO TO 415
C     ***** CHECK FOR 1 TO 3 POINTS INSIDE QUADRANGLE ****
  375 IF(T3-3.0)380,385,385
  380 ITYPE=0
      GO TO 415
C     ***** DETERMINE WHICH EDGES ARE CROSSED BY THE 4 LINE SEGMENTS ***
  385 DO 405 L=1,4
      L1=MOD(L,4)+1
C     ***** LINE SEGMENT L FROM POINT Y1 TO POINT Y2 *****
      Y1(1)=QUAD(L*2-1,IQ)
      Y1(2)=QUAD(L*2,IQ)
      Y2(1)=QUAD(L1*2-1,IQ)
      Y2(2)=QUAD(L1*2,IQ)
      DO 405 K=1,4
      T1=FL(K,L)
      T2=FL(K,L1)
      T3=T1-T2
C     ***** T1 AND T2 MUST HAVE OPPOSITE SIGNS FOR INTERSECTION TO OCCUR
      IF(T1*T2)390,390,405
C     ***** COMPONENT OF SEGMENT L PERPENDICULAR TO EDGE K OF IB IS T3
  390 IF(ABS(T3)-1.E-5)405,405,395
C     ***** CALCULATE COORDINATES OF INTERSECTION *****
  395 T4=(T1*Y2(1)-T2*Y1(1))/T3
      T5=(T1*Y2(2)-T2*Y1(2))/T3
      K0=2*K
      K1=2*(MOD(K,4)+1)
C     ***** IS INTERSECTION WITHIN QUADRANGLE IQ *****
      T6=(T4-QUAD(K0-1,IB))*(QUAD(K1-1,IB)-T4)+(T5-QUAD(K0,IB))*
     1 (QUAD(K1,IB)-T5)
      IF(ABS(T6)-1.E-4)410,410,400
  400 IF(T6)405,410,410
  405 CONTINUE
      GO TO 495
  410 ITYPE=1
C     ***** CHECK OVER/UNDER AMBIGUITY *****
  415 IF((NA1-NB1)*(NA2-NB2)*(NA2-NB1))425,420,425
C     ***** BONDS SHARE AN ATOM *****
  420 IF(NA1+NA2-NB1-NB2)465,495,495
  425 CALL DIFV(ATOMS(2,NB2),ATOMS(2,NB1),V2)
      CALL DIFV(ATOMS(2,NB1),ATOMS(2,NA1),V4)
      CALL UNIT(V2,V2,1)
      CALL UNIT(V4,V4,1)
      CALL NORM(V1,V2,V3,1)
      IF(VV(V3,V3)-TOL)430,430,435
C     ***** PARALLEL BONDS, RECALCULATE V3 *****
  430 CALL NORM(V1,V4,V5,1)
      CALL NORM(V5,V1,V3,1)
C     ***** CHECK FOR COLLINEAR BONDS *****
      IF(VV(V3,V3)-TOL)465,465,450
  435 IF(VV(V3,V4)+TOL)440,450,450
  440 DO 445 J=1,3
  445 V3(J)=-V3(J)
C     ***** V3 IS NORMAL TO BONDS IQ AND IB GOING FROM IQ TOWARD IB ***
  450 IF(VIEW)455,455,460
  455 IF(V3(3))495,495,465
  460 IF(VV(VUE,V3))465,495,495
C     ***** OVERLAPPING BOND FOUND *****
  465 ICQ=ICQ+1
      IF(ITYPE)470,475,475
C     ***** HIDDEN BOND *****
  470 ICQ=-1
      RETURN
C     ***** STORE INTERFERING QUADRANGLE *****
  475 IF(NQOVER-30)485,480,480
  480 NG=17
      CALL ERPNT(TID2,800)
      RETURN
  485 NQOVER=NQOVER+1
      DO 490 K=1,4
      DO 490 J=1,3
  490 QOVER(J,K,NQOVER)=QUA(J,K)
      KQ(NQOVER)=IB
  495 CONTINUE
  500 RETURN
      END
      SUBROUTINE LAPAB(IQ,IA,ICQ,ITY)
C     ***** SUBROUTINE CHECKS FOR OVERLAP BETWEEN ATOMS AND BONDS *****
C     ***** CALLED BY SUBROUTINES LAP700 AND LAP800 *****
      DIMENSION BF(4),CON(3,3),QF(5),QUA(3,4)
      INCLUDE 'SIZES'
      COMMON /OLAP/ CONIC(7,NUMATM), COVER(6,20), KC(20), KQ(30),
     . NCONIC,NCOVER,
     1 NQOVER,NQUAD,OVMRGN,QOVER(3,4,30),QUAD(9,600),SEGM(50,2)

C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      TID=QUAD(9,IQ)
      NA1=TID/1000.
      NA2=AMOD(TID,1000.)
C     ***** ITY.GT.0, CHECK FOR ATOMS OVER A BOND ****
C     ***** ITY.LT.0, CHECK FOR BONDS OVER AN ATOM *****
      ICQ=0
      IF(ITY)210,200,205
  200 RETURN
  205 CALL LAPCON(CONIC(1,IA),CON,V1,0.0)
      IF(ITY-2)220,240,240
  210 IF(ITY+2)220,220,215
  215 CALL LAPCON(CONIC(1,IA),CON,V1,OVMRGN)
C     ***** SET UP LINEAR FORMS FOR EDGES OF QUADRANGLE *****
  220 DO 235 L=1,4
      K=2*L
      K1=MOD(K,8)+2
      QUA(1,L)=QUAD(K,IQ)-QUAD(K1,IQ)
      QUA(2,L)=QUAD(K1-1,IQ)-QUAD(K-1,IQ)
      QUA(3,L)=QUAD(K-1,IQ)*QUAD(K1,IQ)-QUAD(K,IQ)*QUAD(K1-1,IQ)
      T1=SQRT(QUA(1,L)**2+QUA(2,L)**2)
      IF(T1)225,225,230
  225 ITY=0
      ICQ=0
      RETURN
C     ***** TRANSFORM COEFFICIENTS FOR EDGES TO NORMAL FORM *****
  230 DO 235 J=1,3
  235 QUA(J,L)=QUA(J,L)/T1
C     ***** EVALUATE 4 QUADRATIC AND 4 BILINEAR FORMS *****
  240 V2(3)=1.0
      V3(3)=1.0
      T2=3.0
      DO 265 L=1,4
      L1=(MOD(L,4)+1)*2
      V2(1)=QUAD(2*L-1,IQ)
      V2(2)=QUAD(2*L,IQ)
      V3(1)=QUAD(L1-1,IQ)
      V3(2)=QUAD(L1,IQ)
      QF(L)=0.0
      BF(L)=0.0
      DO 250 K=1,3
      T1=CON(3,K)
      DO 245 J=1,2
  245 T1=T1+V2(J)*CON(J,K)
      QF(L)=QF(L)+T1*V2(K)
  250 BF(L)=BF(L)+T1*V3(K)
      IF(QF(L))260,255,265
  255 T2=T2-0.8
      GO TO 265
  260 T2=T2-1.0
  265 CONTINUE
      QF(5)=QF(1)
C     ***** CHECK FOR 4 POINTS OF QUADRANGLE INSIDE OR ON ELLIPSE *****
      IF(T2)270,275,275
  270 ITYPE=-1
      GO TO 330
C     ***** CHECK FOR 1 TO 3 POINTS OF QUADRANGLE INSIDE THE ELLIPSE ***
  275 IF(T2-2.2)280,285,285
  280 ITYPE=0
      IF(NA2-IA)340,375,335
C     ***** CHECK FOR QUADRANGLE-ELLIPSE INTERSECTION *****
  285 DO 305 K=1,4
C     ***** EVALUATE DISCRIMINANT *****
      T1=BF(K)**2-QF(K)*QF(K+1)
      IF(T1)305,305,290
  290 T1=SQRT(T1)
C     ***** IS INTERSECTION WITHIN BOUNDS OF QUADRANGLE *****
      T3=QF(K)-BF(K)
      T4=T3+QF(K+1)-BF(K)
      IF(ABS(T4)-1.E-5)305,305,295
  295 T5=(T3-T1)/T4
      IF(T5)305,280,300
  300 IF(1.0-T5)305,305,280
  305 CONTINUE
C     ***** NO VALID INTERSECTION FOUND *****
C     ***** CHECK FOR CENTER OF ELLIPSE WITHIN THE QUADRANGLE ****
      T3=3.0
      DO 320 K=1,4
      T1=QUA(3,K)
      DO 310 J=1,2
  310 T1=T1+V1(J)*QUA(J,K)
      IF(T1)315,320,320
  315 T3=T3-1.0
  320 CONTINUE
      IF(T3)325,370,370
  325 ITYPE=1
C     ***** CHECK OVER/UNDER AMBIGUITY *****
  330 IF(NA2-IA)375,375,335
  335 IF(IA-NA1)375,375,340
  340 CALL DIFV(ATOMS(2,NA2),ATOMS(2,NA1),V2)
      CALL DIFV(ATOMS(2,IA),ATOMS(2,NA1),V3)
      CALL UNIT(V2,V2,1)
      CALL UNIT(V3,V3,1)
      CALL NORM(V2,V3,V4,1)
      IF(VV(V4,V4)-1.E-5)345,345,350
C     ***** CENTER OF ATOM IQ IS ON THE BOND LINE *****
  345 IF(ITY)370,370,385
C     ***** CENTER OF ATOM IQ IS NOT ON THE BOND LINE *****
  350 CALL NORM(V4,V2,V5,1)
      T1=-V5(3)
      IF(VIEW)365,365,355
  355 T1=V5(3)*(ATOMS(4,IA)-VIEW)
      DO 360 J=1,2
  360 T1=T1+V5(J)*ATOMS(J+1,IA)
  365 IF(T1*(ITY))375,375,370
C     ***** NO INTERFERENCE FOUND *****
  370 ICQ=0
      RETURN
C     ***** ITYPE=1 ENCLOSED ELLIPSE / ITYPE=-1 ENCLOSED QUADRANGLE ****
  375 IF(ITYPE*ITY)380,385,385
C     ***** HIDDEN ATOM OR HIDDEN BOND *****
  380 ICQ=-1
      RETURN
  385 ICQ=1
      IF(ITY)410,390,390
C     ***** STORE INTERFERING ELLIPSE *****
  390 IF(NCOVER-20)400,395,395
  395 NG=17
      CALL ERPNT(ATOMS(1,IA),800)
      NCOVER=NCOVER-1
  400 NCOVER=NCOVER+1
      IJ=1
      DO 405 I=1,3
      DO 405 J=I,3
      COVER(IJ,NCOVER)=CON(I,J)
  405 IJ=IJ+1
      KC(NCOVER)=IA
      RETURN
C     ***** STORE INTERFERING QUADRANGLE *****
  410 IF(NQOVER-30)420,415,415
  415 NG=18
      CALL ERPNT(TID,700)
      NQOVER=NQOVER-1
  420 NQOVER=NQOVER+1
      DO 425 K=1,4
      DO 425 J=1,3
  425 QOVER(J,K,NQOVER)=QUA(J,K)
      KQ(NQOVER)=IQ
      RETURN
      END
      SUBROUTINE LAPCON(CON1,CON,Y,OVMR)
C     ***** TRANSFORM CONIC TO PLOTTER HOMOGENEOUS COORDINATE SYSTEM ***
C     ***** CALLED BY SUBROUTINES LAP700 AND LAPAB *****
      DIMENSION CON1(7),CON(3,3),Y(3)
      Y(1)=(CON1(1)+CON1(2))*0.5
      Y(2)=(CON1(3)+CON1(4))*0.5
      Y(3)=1.0
      CON(1,1)=CON1(5)
      CON(1,2)=CON1(6)
      CON(2,1)=CON1(6)
      CON(2,2)=CON1(7)
      T1=(CON1(2)-CON1(1)+CON1(4)-CON1(3))*0.25
      CON(3,3)=-((T1-OVMR)/T1)**2
      DO 205 K=1,2
      CON(K,3)=0.0
      DO 200 J=1,2
  200 CON(K,3)=CON(K,3)-Y(J)*CON(J,K)
      CON(3,K)=CON(K,3)
  205 CON(3,3)=CON(3,3)-CON(3,K)*Y(K)
      RETURN
      END
      SUBROUTINE LAPDRW(Y,NPEN,NCQ, ICOLOR)
C     ***** SUBROUTINE ELIMINATES HIDDEN LINES AND DRAWS VISIBLE LINES *
      DIMENSION CB(20),CQ(50,2),QL(4,30,2),SEG(2),Y(3),YN(3),YO(3),Z(3)
      INCLUDE 'SIZES'
      COMMON /OLAP/ CONIC(7,NUMATM), COVER(6,20), KC(20), KQ(30),
     . NCONIC,NCOVER,
     1 NQOVER,NQUAD,OVMRGN,QOVER(3,4,30),QUAD(9,600),SEGM(50,2)

C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      NCQ=NCOVER+NQOVER
      IF(NCQ)200,200,205
  200 RETURN
C     ***** CHECK ALL OVERLAPPING ATOMS AND BONDS *****
  205 NPM3=NPEN-3
      IF(NPM3)  210,230,230
C    ***** SAVE INFORMATION FROM LAST POINT IF PEN IS DOWN *****
  210 YO(1)=YN(1)
      YO(2)=YN(2)
      YO(3)=1.0
      NPO=NPN
      DO 215 K=1,NCQ
  215 CQ(K,1)=CQ(K,2)
      IF(NQOVER)230,230,220
  220 DO 225 K=1,NQOVER
      DO 225 J=1,4
  225 QL(J,K,1)=QL(J,K,2)
C     ***** EVALUATE CONIC QUADRATIC FORMS AT NEW POINT YN *****
  230 YN(1)=Y(1)
      YN(2)=Y(2)
      YN(3)=1.0
      NPN=NPEN
      IF(NCOVER)250,250,235
  235 DO 245 K=1,NCOVER
      Z(1)=YN(1)*COVER(1,K)+YN(2)*COVER(2,K)+COVER(3,K)
      Z(2)=YN(1)*COVER(2,K)+YN(2)*COVER(4,K)+COVER(5,K)
      Z(3)=YN(1)*COVER(3,K)+YN(2)*COVER(5,K)+COVER(6,K)
      CQ(K,2)=Z(1)*YN(1)+Z(2)*YN(2)+Z(3)
C     ***** EVALUATE CONIC BILINEAR FORM IF PEN IS DOWN *****
      IF(NPM3)240,245,245
  240 CB(K)= Z(1)*YO(1)+Z(2)*YO(2)+Z(3)
  245 CONTINUE
C     ***** EVALUATE LINEAR FORMS AND SIGNATURE FOR QUADRANGLE *****
  250 IF(NQOVER)275,275,255
  255 KCQ=NCOVER
      DO 270 K=1,NQOVER
      T2=3.0
      DO 265 J=1,4
      T1=YN(1)*QOVER(1,J,K)+YN(2)*QOVER(2,J,K)+QOVER(3,J,K)
      IF(T1)260,265,265
  260 T2=T2-1.0
  265 QL(J,K,2)=T1
      KCQ=KCQ+1
C     ***** T2=-1 INSIDE, =0 ACROSS ANY EDGE, =1 ACROSS ANY VERTEX *****
  270 CQ(KCQ,2)=T2
C     ***** IF PEN IS UP, OMIT ALL SUBSEQUENT CHECKING *****
  275 IF(NPM3)285,280,280
  280 NPN=3
      CALL SCRIBE(YN,NPN,LTNO, ICOLOR)
      RETURN
C     ***** CHECK FOR HIDDEN SEGMENT *****
  285 DO 295 K=1,NCQ
      IF(CQ(K,1))290,295,295
  290 IF(CQ(K,2))280,295,295
  295 CONTINUE
C     ***** FIND ENTRY AND EXIT POINTS ON EACH CONIC *****
      NINT=0
      IF(NCOVER)330,330,300
  300 DO 325 K=1,NCOVER
C     ***** EVALUATE DISCRIMINANT *****
      T1=CB(K)**2-CQ(K,1)*CQ(K,2)
      IF(T1)325,325,305
  305 T1=SQRT(T1)
C     ***** SOLVE QUADRATIC EQATION *****
      T2=CQ(K,1)-CB(K)
      T3=T2+CQ(K,2)-CB(K)
      IF(ABS(T3)-1.E-5)325,325,310
  310 T4=(T2-T1)/T3
      T5=(T2+T1)/T3
C     ***** VALID INTERSECTION IF T4.LT.1 AND T5.GT.0 *****
      IF(T4-1.0)315,325,325
  315 IF(T5)325,325,320
C     ***** SAVE VALID CONIC INTERSECTIONS *****
  320 NINT=NINT+1
      SEGM(NINT,1)=T4
      SEGM(NINT,2)=T5
  325 CONTINUE
  330 IF(NQOVER)425,425,335
C     ***** FIND ENTRY AND EXIT POINTS FOR EACH QUADRANGLE *****
  335 DO 420 K=1,NQOVER
      I12=0
      KCQ=NCOVER+K
C     ***** CHECK FOR SINGLE INSIDE POINT *****
      SEG(1)=CQ(KCQ,1)
      IF(SEG(1))345,340,340
  340 SEG(1)=1.0-CQ(KCQ,2)
      IF(SEG(1)-1.0)350,350,345
C     ***** INSIDE POINT FOUND, ONLY ONE INTERSECTION POSSIBLE *****
  345 I12=1
C     ***** FIND WHICH EDGES ARE CROSSED BY THE SEGMENT *****
  350 DO 410 J=1,4
      T1=QL(J,K,1)
      T2=QL(J,K,2)
      T3=T1-T2
      IF(T1*T2)355,355,410
C     ***** CHECK FOR SEGMENT ON AN EDGE *****
  355 IF(ABS(T3)-1.E-5)420,420,360
C     ***** CALCULATE COORDINATES OF INTERSECTION *****
  360 T4=(T1*YN(1)-T2*YO(1))/T3
      T5=(T1*YN(2)-T2*YO(2))/T3
      J1=2*(MOD(J,4)+1)
      IQ=KQ(K)
C     ***** IS INTERSECTION WITHIN LIMITS OF QUADRANGLE *****
      T6=(T4-QUAD(2*J-1,IQ))*(QUAD(J1-1,IQ)-T4)+(T5-QUAD(2*J,IQ))*
     1 (QUAD(J1,IQ)-T5)
      IF(ABS(T6)-1.E-4)370,370,365
  365 IF(T6)410,370,370
C     ***** CALCULATE FRACTION PARAMETER AND STORE IT *****
  370 T1=T1/T3
      IF(I12-1)375,380,395
C     ***** STORE FIRST INTERSECTION *****
  375 I12=1
      GO TO 390
C     ***** STORE SECOND INTERSECTION ****
  380 I12=2
      IF(T1-SEG(1))385,405,405
  385 SEG(2)=SEG(1)
  390 SEG(1)= T1
      GO TO 410
C     ***** MORE THAN TWO INTERSECTIONS (I.E.,QUADRANGLE DIAGONAL) *****
  395 IF(T1-SEG(1))390,410,400
  400 IF(T1-SEG(2))410,410,405
  405 SEG(2)=T1
  410 CONTINUE
      IF(I12-1)420,420,415
C     ***** STORE FRACTION PARAMETERS *****
  415 NINT=NINT+1
      SEGM(NINT,1)=SEG(1)
      SEGM(NINT,2)=SEG(2)
  420 CONTINUE
C     ***** END OF ENTRY-AND-EXIT-POINT CALCULATIONS *****
  425 IF(NINT-1)430,490,435
C     ***** NO INTERFERENCE FOUND, DRAW ENTIRE SEGMENT *****
  430 CALL SCRIBE(YN,2,LTNO, ICOLOR)
      RETURN
C     **** SORT SEGMENT INTERSECTION LIST *****
C     ***** SORTING PROCEDURE BY SHELL,D.L. COMM. ACM 2,30-32 (1959) ***
  435 M=NINT
  440 M=M/2
      IF(M)490,490,445
  445 K=NINT-M
      J=1
  450 I=J
  455 IM=I+M
      IF(SEGM(I,1))460,470,470
  460 IF(SEGM(IM,1))465,465,485
  465 IF(SEGM(I,2)-SEGM(IM,2))485,485,475
  470 IF(SEGM(I,1)-SEGM(IM,1))485,485,475
  475 DO 480 L=1,2
      T1=SEGM(I,L)
      SEGM(I,L)=SEGM(IM,L)
  480 SEGM(IM,L)=T1
      I=I-M
      IF(I)485,485,455
  485 J=J+1
      IF(J-K)450,450,440
C     ***** FIND STARTING POINT P0 AND END POINT P1 *****
  490 P0=0.0
      K=0
  495 K=K+1
      IF(K-NINT)500,500,515
  500 P1=SEGM(K,1)
      IF(P1)510,505,505
  505 IF(P1-P0)510,510,520
  510 P0=AMAX1(P0,SEGM(K,2))
      IF(P0-1.0)495,530,525
  515 P1=1.0
C     ***** DRAW SEGMENT FROM P0 TO P1 *****
  520 IF(P0)535,535,530
  525 P0=1.0
  530 Z(1)=YO(1)*(1.-P0)+YN(1)*P0
      Z(2)=YO(2)*(1.-P0)+YN(2)*P0
      NPN=3
      CALL SCRIBE(Z,NPN,LTNO, ICOLOR)
      IF(P0-1.0)535,540,540
  535 Z(1)=YO(1)*(1.-P1)+YN(1)*P1
      Z(2)=YO(2)*(1.-P1)+YN(2)*P1
      NPN=2
      CALL SCRIBE(Z,NPN,LTNO, ICOLOR)
      IF(P1-1.0)510,540,540
  540 RETURN
      END
      SUBROUTINE MERGE
C SUBROUTINE TO MERGE ATOMS ON DRUM,FILE39,WITH ATOMS IN ATOM LIST,
C THROWING AWAY ALL DUPLICATES.

      INCLUDE 'SIZES'
      CHARACTER*80 COMMNT
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      DIMENSION BATOMS(4)
      IF(LATM.EQ.999) RETURN
      KOUNT=0
      KKOUNT=0
C  READ ATOMS FROM FILE NTEMPT
 4    READ(NTEMPT,1,END=6)(BATOMS(I),I=1,4)
  1   FORMAT(F10.0,3F10.6)
      KOUNT=KOUNT+1
C CHECK ORTEP CODES TO SEE IF ALREADY IN ATOMS LIST
      IF(LATM.EQ.0) GO TO 20
      DO 3 L=1,LATM
      IF(ABS(BATOMS(1)-ATOMS(1,L)).LT..1) GO TO 4
 3    CONTINUE
C IF NOT,COPY INTO LIST
 20   LATM=LATM+1
      DO 5 LA=1,4
  5   ATOMS(LA,LATM)=BATOMS(LA)
      KKOUNT=KKOUNT+1
C CHECK TOTAL POSSIBLE NUMBER OF ATOMS IS NOT EXCEEDED.
      IF(LATM.LE.999) GOTO4
      WRITE ( COMMNT,11) NTEMPT
      CALL DEBUGR( COMMNT(1:78) )
  11  FORMAT('  NUMBER OF ATOMS IN LIST EXCEEDS MAXIMUM,STOPPED IN SUBRO
     1UTINE MERGE WHILE READING FILE',I4)
      STOP
C END OF FILE REACHED ON DRUM
  6   CONTINUE
      WRITE ( COMMNT,10)KOUNT,NTEMPT,KKOUNT
      CALL DEBUGR( COMMNT(1:78) )
 10   FORMAT(I10,' ATOMS READ FROM FILE',I4,',',I4,' ATOMS ADDED TO LIST
     1')
      RETURN
      END
      SUBROUTINE MM(X,Y,Z)
C     MULTIPLY TWO MATRICES
C     Z(3,3)=X(3,3)*Y(3,3)
      DIMENSION X(3,3),Y(3,3),Z(3,3)
      DO 117 I=1,3
      DO 117 K=1,3
      Z(I,K)=0.0
      DO 117 J=1,3
  117 Z(I,K)=Z(I,K)+X(I,J)*Y(J,K)
      RETURN
      END
      SUBROUTINE MV(X,Y,Z)
C     MATRIX * VECTOR
C     Z(3)=X(3,3)*Y(3)
      DIMENSION X(3,3),Y(3),Z(3)
      DO 113 I=1,3
      Z(I)=0.0
      DO 113 J=1,3
  113 Z(I)=Z(I)+X(I,J)*Y(J)
      RETURN
      END

      SUBROUTINE DORTEP
C * BLOCK COMMON *
      INCLUDE 'SIZES'
      INTEGER*2 ATBOND
      CHARACTER*80  FILEIN, FILOUT, FILPLT, LLEGND
      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      REAL CO,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,CMIN,CMAX
      REAL SCALE, BSCALE, ATCHG
      REAL VDIST, XBOUND, YBOUND, BORDER, RETRAC, XCENTR, YCENTR
      REAL SCOR1, SCOR2, REDUCT, BTHICK, XTITLE, YTITLE, THEIT
      COMMON /CORTEP/ VDIST, XBOUND, YBOUND, BORDER, RETRAC,
     .   IN600, XCENTR, YCENTR, SCOR1, SCOR2, REDUCT,
     .  NORATM, NORBON, NRBOND, BTHICK, XTITLE, YTITLE, THEIT
*  VDIST :== VIEWING DISTANCE           (INST 301)
* XBOUND :== X-BOUNDARY                 (INST 301)
* YBOUND :== Y-BOUNDARY                 (INST 301)
* BORDER :== SIZE OF BORDER             (INST 301)
* RETRAC :== DIPLACEMENT FOR RETRACE    (INST 303)
* IN600  :== TYPE OF 600 COMMAND        (INST 60X)
* XCENTR :== X-COORD CENTER OF PLOT     (INST 60X)
* YCENTR :== Y-COORD CENTER OF PLOT     (INST 60X)
* SCOR1  :== OVER ALL SCALING           (INST 60X)
* SCOR2  :== SUBSIDIARY SCALING         (INST 60X)
* REDUCT :== OVERALL REDUCTION          (INST 611)
* NORATM :== ATOM SYMBOL SHAPE          (INST 7XX)
* NORBON :== TYPE OF BONDS              (INST 80X)
* NRBOND :== NUMBER OF LINES IN BOND    (INST 80X)
* BTHICK :== BOND THICKNESS             (INST 80X)
* XTITLE :== X-RESET FOR TITLE          (INST 902)
* YTITLE :== Y-RESET FOR TITLE          (INST 902)
* THEIT  :== TITLE HEIGHT               (INST 902)
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *
      COMMON /PLOTS/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,CMIN,CMAX,SCALE
      COMMON /ATOMS/ CO(3,NUMATM),IE(NUMATM),NATOMS, ATCHG( NUMATM)
      COMMON /DISPLY/ IREM(200), BSCALE, ATBOND( NUMATM, NUMATM),
     .                ISTYPE, LATYPE, IMASK( NUMATM), ISCOLO
C?      COMMON /ATMRAD/ ATMRAD(200)
      COMMON /LEGEND/  FILEIN, FILOUT, FILPLT, LLEGND
      CHARACTER HEADER*80, CHEM(NUMATM)*6
      COMMON /STRING/ HEADER, CHEM
C
      IF (NATOMS .LT. 2) THEN
        CALL DEBUGR('Too few atoms.  ORTEP not permitted.')
        RETURN
      ENDIF
    2 CALL PRIME
C     ***** READ JOB TITLE CARD *****
C#      READ(IN,4,END=3010) TITLE
    4 FORMAT(18A4)
C#      WRITE (NOUT,6)(TITLE(I),I=1,18)
    5 FORMAT(1H010X,18A4)
    6 FORMAT(1H110X,18A4)
      CALL DRELIM
      ISAVE=0
      GO TO 8
    7 ISAVE=0
      NOOUT=0
C     ***** ZERO AIN ARRAY *****
    8 DO 10 J=1,140
   10 AIN(J)=0.
C  LOOP THROUGH ALL ATOMS AND SEE IF ALLOWED IN PICTURE
      NNATOM = 0
      DO 20 I=1,NATOMS
         IF ( IMASK( I) .NE. 0) GO TO 20
         IF ( IREM( IE( I)) .NE. 0) GO TO 20
         NNATOM = NNATOM + 1
20    CONTINUE
C     ***** READ NEW INSTRUCTION CARD *****
C  COMMAND 101 155501 2 1 NUMAT 3.61
      NF=101
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      AIN(1)=155501
      AIN(2)=2
      AIN(3)=1
      AIN(4)=NNATOM
      AIN(5)=3.61
      AIN(6)=0
      AIN(7)=0
      CALL SEARC
C  COMMAND 102 155501 2 1 NUMAT 1.8
      NF=102
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      AIN(1)=155501
      AIN(2)=2
      AIN(3)=1
      AIN(4)= NNATOM
      AIN(5)=1.8
      AIN(6)=0
      AIN(7)=0
      CALL SEARC
C COMMAND 201
      NF=201
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      AIN(1)=0
      AIN(2)=0
      AIN(3)=0
      AIN(4)=0
      AIN(5)=0
      AIN(6)=0
      AIN(7)=0
      CALL F200
C  COMMAND 301  12.  12.  20.  2.5  (CHANGE 15 JULY 1986)
C  COMMAND 301  10.  10.  24.  0.5
      NF=301
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      XLNG(1)= 12.0
      XLNG(1)= XBOUND
      XLNG(2)= 12.0
      XLNG(2)= YBOUND
      VIEW= 20.0
      VIEW= VDIST
      BRDR= 2.5
      BRDR= BORDER
      AIN(1)=0
      AIN(2)=0
      AIN(3)=0
      AIN(4)=0
      AIN(5)=0
      AIN(6)=0
      AIN(7)=0
C COMMAND 303
      DISP=0.0
      DISP = RETRAC
C COMMAND 401 155501 -NATOMS*100000+55501
      NF=401
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      AIN(1)=155501
      AIN(2)=-( NNATOM *100000+55501)
      AIN(3)=0
      AIN(4)=0
      AIN(5)=0
      AIN(6)=0
      AIN(7)=0
      LATM=0
      CALL F400
      LATM= NNATOM
C COMMAND 501
      NF=501
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      AIN(1)= ( NNATOM +1)*100000+55501
      AIN(2)= ( NNATOM +1)*100000+55501
      AIN(3)= ( NNATOM +2)*100000+55501
      AIN(4)= ( NNATOM +1)*100000+55501
      AIN(5)= ( NNATOM +3)*100000+55501
      AIN(6)=0
      AIN(7)=0
      AIN(8)=1
      CALL F500
C COMMAND 604
      NF=604
      NF=600 + IN600
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      AIN(1)=0.0
      AIN(1) = XCENTR
      AIN(2)=0.0
      AIN(2) = YCENTR
      AIN(3)=0
      TEMPSC = XMAX-XMIN
      IF( YMAX-YMIN .GT. TEMPSC) TEMPSC = YMAX-YMIN
      IF ( IN600 .EQ. 3) THEN
C?      SCOR1 = 4.0 - TEMPSC/3.5
C?      SCOR1 = 3.9 - TEMPSC/3.6
        SCOR1 = 3.9 - TEMPSC/3.4
        IF (SCOR1.LT. 0.3150) SCOR1 = 0.315
      ELSE
        SCOR1 = 0.0
      ENDIF
C     WRITE (*,*) 'IN DORTEP: XSIZE=',XMAX-XMIN,', YSIZE=',YMAX-YMIN,
C    .   ', ZSIZE=',ZMAX-XMIN
C     WRITE (*,*) 'IN DORTEP: TEMPSC =', TEMPSC, ',  SCOR1 = ',SCOR1
      AIN(3) = SCOR1
      AIN(4)=1.54
      AIN(4) = SCOR2
      AIN(5)=0
      AIN(6)=0
      AIN(7)=0
      AIN(8)=0
      CALL F600
      IF (DEBUGO) THEN
          CALL OPLOT( 0.0, 0.0, 8)
          WRITE (NOUT,609) XO(1),XO(2),SCAL1,SCAL2
609          FORMAT (1X,'Orig (in plt coord)=',
     1      2F8.2,' in.; Scal=',F6.2,' in/A, ',
     2      'Sph. scal=',F6.2,' in/A' )
      ENDIF
C COMMAND 611
      NF=611
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      AIN(1)=0
      AIN(2)=0
      AIN(3)=0.9
      AIN(3) = REDUCT
      AIN(4)=0
      AIN(5)=0
      AIN(6)=0
      AIN(7)=0
      AIN(8)=0
      CALL F600
      IF (DEBUGO) THEN
          CALL OPLOT( 0.0, 0.0, 8)
          WRITE (NOUT,609) XO(1),XO(2),SCAL1,SCAL2
      ENDIF
C COMMAND 511
      NF=511
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      AIN(1)=1
      AIN(2)=0
      AIN(3)=0
      AIN(4)=0
      AIN(5)=0
      AIN(6)=0
      AIN(7)=0
      AIN(8)=0
      KD(1,1)=1
      KD(2,1)= NNATOM
      KD(3,1)=1
      KD(4,1)= NNATOM
      KD(5,1)=4
      CD(1,1)=.9
      CD(2,1)=1.6
      CD(3,1)=.04
      CD(4,1)=0
      CD(5,1)=0
      NCD=1
      CALL F500
      NCD=0
C COMMAND 712
      NF = 712
      NF = 710 + NORATM
      AIN( 1) = 4.
      AIN( 2) = 0
      AIN( 3) = 0
      AIN( 4) = 0
      AIN( 5) = 0
      AIN( 6) = 0
      AIN( 7) = 0
      AIN( 8) = 0
      IF ( LATYPE .GT. 0 ) THEN
         NF = 714
         AIN( 1) = 1.0
         AIN(5) = 0.20
         AIN(6)=.40
C?         AIN(7)=.30
      ENDIF
      NF1=NF
      NJ=NF1/100
      NJ2=NF1-NJ*100
      NJ3=MOD(NJ2,10)
      CALL F700
C COMMAND 811
      IAINI = -1
      IA=0
      DO 810 IA1=1, NATOMS-1
         IF ( IMASK( IA1) .NE. 0) GO TO 810
         IF ( IREM( IE( IA1)) .NE. 0) GO TO 810
         IA = IA + 1
         IB = IA
         DO 800 IB1=IA1+1, NATOMS
            IF ( IMASK( IB1) .NE. 0) GO TO 800
            IF ( IREM( IE( IB1)) .NE. 0) GO TO 800
            IB = IB + 1
            IF ( ATBOND( IA1, IB1 ) .GT. 0 ) THEN
               IAINI = IAINI + 2
               AIN(IAINI)= IA*100000+55501
               AIN(IAINI+1)= IB*100000+55501
               IF (IAINI .GT. 137) THEN
                 NF=811
                 NF=810 + NORBON
                 NF1=NF
                 NJ=NF1/100
                 NJ2=NF1-NJ*100
                 NJ3=MOD(NJ2,10)
                 KD(1,1)=0
                 KD(2,1)=0
                 KD(3,1)=0
                 KD(4,1)=0
                 KD(5,1)=4
                 KD(5,1) = NRBOND
                 CD(1,1)=0
                 CD(2,1)=0
                 CD(3,1)=.04
                 CD(3,1) = BTHICK
                 CD(4,1)=0
                 CD(5,1)=0
                 NCD=1
                 CALL F800
                 DO 809 IZZ=1,140
809              AIN(IZZ)=0.0
                 IAINI = -1
              ENDIF
         ENDIF
800      CONTINUE
810      CONTINUE
      IF (IAINI .GT. 0) THEN
            NF=811
            NF=810 + NORBON
            NF1=NF
            NJ=NF1/100
            NJ2=NF1-NJ*100
            NJ3=MOD(NJ2,10)
            KD(1,1)=0
            KD(2,1)=0
            KD(3,1)=0
            KD(4,1)=0
            KD(5,1)=4
            KD(5,1) = NRBOND
            CD(1,1)=0
            CD(2,1)=0
            CD(3,1)=.04
            CD(3,1) = BTHICK
            CD(4,1)=0
            CD(5,1)=0
            NCD=1
            CALL F800
      ENDIF
      NCD=0
C COMMAND 902
          NF=902
          NF1=NF
          NJ=NF1/100
          NJ2=NF1-NJ*100
          NJ3=MOD(NJ2,10)
          AIN(1)= 0.0
          AIN(2)= 0.0
* reset for x-edge
C?          AIN(3)=0.0
c?2          AIN(3)= -0.25
          AIN(3)= -0.15
          AIN(3)= XTITLE
* reset for y-edge
C?          AIN(4)=-.25
c?          AIN(4)= -1.25    ! on 8 apr 88
          AIN(4) = -1.30
          AIN(4) = YTITLE
* character height
C?          AIN(5)=.25
c?2          AIN(5)=.40
          AIN(5)=.35
          AIN(5) = THEIT
          AIN(6)=0
          AIN(7)=0
          AIN( 8) = 0
          HEADER = LLEGND
          CALL F900
      RETURN
      END

      SUBROUTINE DRELIM
      REAL CO,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,CMIN,CMAX
      REAL SCALE, BSCALE, ATCHG
      REAL ATMRAD
      INTEGER*2 ATBOND
      CHARACTER*6 ATSYMB
      CHARACTER*80 COMMNT
      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
C     ***** DATA INPUT ROUTINE *****
      INCLUDE 'SIZES'
      DIMENSION B(9),PALIN(9, NUMATM)
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *
      EQUIVALENCE (PALIN(1,1),PA(1,1,1))
      COMMON /AINID/ IDAIN( NUMATM)
      COMMON /ATOMS/ CO(3,NUMATM), IE(NUMATM), NATOMS, ATCHG( NUMATM)
      COMMON /ATSYMB/ ATSYMB( 200)
      COMMON /COLORS/ ICOLAT( 200)
      COMMON /PLOTS/ XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,CMIN,CMAX,SCALE
      COMMON /DISPLY/ IREM(200), BSCALE, ATBOND( NUMATM, NUMATM),
     .                ISTYPE, LATYPE, IMASK( NUMATM), ISCOLO
      COMMON /ATMRAD/ ATMRAD(200)
      CHARACTER*6 ATLABS
      COMMON /LABELS/ ATLABS( NUMATM)
      CHARACTER HEADER*80, CHEM(NUMATM)*6
      COMMON /STRING/ HEADER, CHEM
C     ***** CELL DIMENSIONS *****
C***********************************************************************
      TSCALE=SCALE
      IF (ZMAX-ZMIN .GT. TSCALE) TSCALE = ZMAX-ZMIN
      TSCALE = TSCALE*2.0
C      IF (TSCALE .LT. 10.0) TSCALE = 10.0
      SCALE = 11.0
      IF (DEBUGO) THEN
          WRITE ( COMMNT,9000) SCALE,ZMAX-ZMIN,TSCALE
          CALL DEBUGR( COMMNT(1:78) )
          WRITE ( COMMNT,'('' XMIN='',F12.5,'' YMIN='',F12.5,
     1          '' ZMIN='',F12.5)') XMIN, YMIN, ZMIN
          CALL DEBUGR( COMMNT(1:78) )
9000      FORMAT (1X,'MYPRELIM: SCALE=',F7.3,', ZMAX-ZMIN=',F7.3,
     1      'TSCALE=',F7.3 )
      ENDIF
      A(1)=TSCALE
      A(2)=TSCALE
      A(3)=TSCALE
      A(4)=90.0
      A(5)=90.0
      A(6)=90.0
      T1=ABS(A(4))-1.
      DO 125 J=1,3
      IF(T1)115,110,110
C     ***** CELL ANGLES IN DEGREES *****
  110 A(J+6)=A(J+3)
      A(J+3)=COS(A(J+6)*.01745329252)
      GO TO 120
C     ***** COSINES OF CELL ANGLES *****
  115 A(J+6)=ACOS(A(J+3))
C     ***** STORE IDEMFACTOR MATRIX *****
  120 AID(J,J)=1.
      AID(J+1,1)=0.
      AID(J+5,1)=0.
C     ***** STORE METRIC TENSOR *****
  125 AA(J,J)=A(J)**2
      AA(1,2)=A(1)*A(2)*A(6)
      AA(1,3)=A(1)*A(3)*A(5)
      AA(2,3)=A(2)*A(3)*A(4)
      AA(2,1)=AA(1,2)
      AA(3,1)=AA(1,3)
      AA(3,2)=AA(2,3)
C     ***** INVERT METRIC TENSOR *****
      CALL AXEQB(AA,BB,AID,3)
C     ***** CALCULATE RECIPROCAL CELL PARAMETERS *****
      DO 128 J=1,3
  128 B(J)=SQRT(BB(J,J))
      B(6)=BB(1,2)/(B(1)*B(2))
      B(5)=BB(1,3)/(B(1)*B(3))
      B(4)=BB(2,3)/(B(2)*B(3))
      DO 130 J=1,3
  130 B(J+6)=ACOS(B(J+3))
C     ***** WAS INPUT FOR REAL OR RECIPROCAL CELL *****
      IF(A(1)-1.)135,150,150
  135 DO 140 J=1,9
      T1=AA(J,1)
      AA(J,1)=BB(J,1)
      BB(J,1)=T1
      T1=A(J)
      A(J)=B(J)
  140 B(J)=T1
C     ***** WRITE OUT CELL PARAMETERS *****
150      CONTINUE
C     ***** STORE STANDARD VECTORS *****
      CALL AXES(AID,AID(1,2),REFV,0)
      CALL MM(AA,REFV,AAREV)
      DO 160 I=1,3
      DO 160 J=1,3
      AAWRK(J,I)=AAREV(J,I)
      Q(J,I)=REFV(I,J)
  160 WRKV(J,I)=REFV(J,I)
C     ***** READ AND WRITE SYMMETRY TRANSFORMATIONS *****
      LINES=14
C#      DO 190 I=1,48
      I=1
      LINES=MOD(LINES+1,56)
C#      READ (IN,173)IS,(TS(J,I),(FS(K,J,I),K=1,3),J=1,3)
      TS(1,1)=0
      TS(2,1)=0
      TS(3,1)=0
      FS(1,1,1)=1
      FS(2,1,1)=0
      FS(3,1,1)=0
      FS(1,2,1)=0
      FS(2,2,1)=1
      FS(3,2,1)=0
      FS(1,3,1)=0
      FS(2,3,1)=0
      FS(3,3,1)=1
      IS=1
C     ***** NON-CRYSTALLOGRAPHIC HELIX-SYMMETRY INPUT *****
      IF(FS(3,3,I)-5.)188,186,186
  186 T1=FS(1,3,I)/FS(3,3,I)
      TS(3,I)=TS(3,I)+T1
      T1=AMOD(T1*FS(2,3,I),1.)*6.28318531
      T2=COS(T1)
      T1=SIN(T1)
      DO 187 J=1,9
  187 VT(J,1)=AID(J,1)
      VT(1,1)=T2
      VT(2,2)=T2
      VT(2,1)=-T1
      VT(1,2)=T1
      CALL MM(VT,Q,PAC)
      CALL MM(AAREV,PAC,FS(1,1,I))
  188 IF(IS)195,190,195
  190 CONTINUE
      NG=1
      CALL ERPNT(0.,0)
      I=48
  195 NSYM=I
C     ***** POSITIONAL AND THERMAL PARAMETERS *****
      NNATOM = 0
  225 DO 245 I=1,NATOMS
         IF ( IMASK( I) .NE. 0) GO TO 245
         IF ( IREM( IE( I)) .NE. 0) GO TO 245
         NNATOM = NNATOM + 1
         IF (NNATOM .EQ. 1) THEN
            TEMPX = CO(1,I)
            TEMPY = CO(2,I)
            TEMPZ = CO(3,I)
         ENDIF
C***********************************************************************
C     READ ATOM PARAMETERS IN XRAY 67 FORMAT
C#      READ(IN,211)CHEM(I),(P(J,I),J=1,3),T1
         ITEMP = ICOLAT( IE( I) )
         IF ( ITEMP .GT. ISCOLO ) ITEMP = 1
         IDAIN( NNATOM) = ITEMP
           CHEM( NNATOM) = ATLABS( I)
           ITEMP = INDEX( CHEM( NNATOM), '"')
           IF ( ITEMP.GT. 0) CHEM( NNATOM)( ITEMP: ITEMP) = ' '
           ITEMP = INDEX( CHEM( NNATOM), '"')
           IF ( ITEMP.GT. 0) CHEM( NNATOM)( ITEMP: ITEMP) = ' '
C***********************************************************************
C     READ THERMAL PARAMETERS IN XRAY 67 FORMAT
C#      READ(IN,213)(PA(J,1,I),J=1,7),IS
         PA(1,1,NNATOM) = ATMRAD( IE(I))/ 3.2
         PALIN(7,NNATOM)=7
         IS=1
C     ***** TYPE 1 POSITIONAL PARAMETERS (ANGSTROMS) *****
         P(1,NNATOM)=(CO(1,I)-TEMPX)/TSCALE
         P(2,NNATOM)=(CO(2,I)-TEMPY)/TSCALE
         P(3,NNATOM)=(CO(3,I)-TEMPZ)/TSCALE
         IF (DEBUGO) THEN
            WRITE ( COMMNT,3030) (CO(J,I),J=1,3),(P(J,NNATOM),J=1,3)
            CALL DEBUGR( COMMNT(1:78) )
 3030    FORMAT (1X, 'COORD:',3F10.4,'; SCALED:',3F10.5)
      ENDIF
C     ***** TYPE 2 POSITIONAL PARAMETERS, STANDARD CARTESIAN *****
C        V1(1) = CO(1,I)-SIGN(XMAX-XMIN)*XMIN
C        V1(2) = CO(2,I)-YMIN
C        V1(3) = CO(3,I)-ZMIN
C        CALL VM(V1,Q,P(1,I))
C     ***** TYPE 3 POSITIONAL PARAMETERS *****
C     ***** CYLINDRICAL COORDINATES REFERRED TO STANDARD CARTESIAN *****
C#244      IF(IS)246,245,246
  245    CONTINUE
  246 CONTINUE
      I =NNATOM + 1
      CHEM(I)='ORGN'
      P(1,I)=0
      P(2,I)=0
      P(3,I)=1.0 / TSCALE
      I=I+1
      CHEM(I)='DUMX'
      P(1,I)=1.0 / TSCALE
      P(2,I)=0.0
      P(3,I)=1.0 / TSCALE
      I=I+1
      CHEM(I)='DUMY'
      P(1,I)=0.0
      P(2,I)=0.0
      P(3,I)=0.0
      NATOM=I
C     ***** CONVERT TEMP FACTOR COEF TO STANDARD TYPE ZERO *****
      NG1=0
      DO 450 I=1,NATOM
      T1=PA(1,1,I)
      K=1.+PALIN(7,I)
      IF(T1)250,250,255
  250 T1=.1
      GO TO 405
  255 T6=.0506605918
      GO TO(270,260,265,265,270,260,400,405,270,260,450),K
C     ***** TYPE 1 *****
  260 DO 262 J=4,6
C?  262 PA(J,1,I)=PA(J,1,I)*.5
  262 PALIN(J,I)=PALIN(J,I)*.5
      GO TO 270
C     ***** TYPES 2 AND 3 (BASE 2 SYSTEMS) *****
  265 T6=.351152464
      IF(K-4)270,260,270
C     ***** TYPES 0 THROUGH 5 *****
  270 IF(PA(2,1,I))400,400,272
  272 DO 300 J=1,3
      DO 300 L=J,3
      T2=T6
      IF(K-5)285,275,275
  275 IF(K-6)280,280,281
C     ***** TYPES 4 AND 5 *****
  280 T2=B(J)*B(L)*T2*.25
      GO TO 285
C     ***** TYPES 8 AND 9 (U(I,J) TENSOR SYSTEMS) *****
  281 T2=B(J)*B(L)
  285 IF(J-L)290,287,290
  287 VT(J,J)=T2*PA(J,1,I)
      GO TO 300
  290 M=J+L+1
C?      VT(J,L)=T2*PA(M,1,I)
      VT(J,L)=T2*PALIN(M,I)
      VT(L,J)=VT(J,L)
  300 CONTINUE
C     ***** FIND PRINCIPAL AXES *****
      CALL MM(VT,AA,DA)
      CALL EIGEN(DA,RMS,PAT)
C     ***** ARE EIGENVALUES POSITIVE *****
      IF(RMS(1))325,325,320
  320 IF(NG)350,360,330
  325 NG=3
  330 NG1=1
      CALL ERPNT((I)*100000.+55501.,0)
C     ***** 3 EQUAL EIGENVALUES, USE REFERENCE VECTORS *****
  340 T3=SIGN(SQRT(ABS(RMS(1)+RMS(2)+RMS(3))/3.),RMS(1))
      DO 345 J=1,3
      DO 342 K=1,3
  342 PA(J,K,I)=REFV(J,K)
  345 EV(J,I)=T3
      GO TO 450
  350 IF(NG+6)340,340,352
C     ***** TWO EQUAL EIGENVALUES *****
  352 N=NG+5
      CALL UNIT(PAT(1,N),V1,-1)
      DO 354 K=1,3
      IF(ABS(VMV(V1,AA,REFV(1,K)))-.58)356,354,354
  354 CONTINUE
  356 CALL MM(AA,DA,VT)
      CALL AXES(V1,REFV(1,K),DA,-1)
      DO 359 K=1,3
      L=MOD(N+K-2,3)+1
      DO 358 J=1,3
  358 PA(J,L,I)=DA(J,K)
  359 EV(L,I)=SIGN(SQRT(ABS(VMV(DA(1,K),VT,DA(1,K)))),RMS(L))
      GO TO 450
C     ***** MAKE EIGENVECTORS 1 ANGSTROM LONG *****
  360 DO 365 J=1,3
  365 CALL UNIT(PAT(1,J),PA(1,J,I),-1)
  370 NG=0
C     ***** SQRT EIGENVALUE = RMS DISPLACEMENT *****
      DO 375 J=1,3
      T2=RMS(J)
  375 EV(J,I)=SIGN(SQRT(ABS(T2)),T2)
      GO TO 450
C     ***** TYPE 6 (ISOTROPIC TEMP FACTOR) *****
  400 T1=SQRT(T1*.012665148)
C     ***** TYPE 7 (DUMMY SPHERE) *****
  405 DO 410 J=1,3
  410 EV(J,I)=T1
      IF(PA(3,1,I))430,430,415
C     ***** DEFINED VECTORS FOR SPHERE *****
  415 DO 425 J=1,4
C?      T2=PA(J+2,1,I)
      T2=PALIN(J,I)
      CALL ATOM(T2,VT(1,J))
      IF(NG)420,425,420
  420 CALL ERPNT(T2,0)
      GO TO 430
  425 CONTINUE
      CALL DIFV(VT(1,2),VT(1,1),V1)
      CALL DIFV(VT(1,4),VT(1,3),V2)
      CALL AXES(V1,V2,PA(1,1,I),-1)
      GO TO 450
C     ***** REFERENCE VECTORS FOR SPHERE *****
  430 DO 435 J=1,9
C?  435 PA(J,1,I)=REFV(J,1)
  435 PALIN(J,I)=REFV(J,1)
  450 NG=0
  999 RETURN
      END
      SUBROUTINE NORM(X,Y,Z,ITYPE)
C     ***** VECTOR PRODUCT  Z=X*Y *****
C     ***** ITYPE .GT.0 FOR CARTESIAN,.LE.0 FOR TRICLINIC *****
      DIMENSION X(3),Y(3),Z(3),Z1(3)

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      DO 125 I=1,3
      I1=MOD(I+3,3)+1
      I2=MOD(I+1,3)+1
      T1=X(I1)*Y(I2)-X(I2)*Y(I1)
      IF(ITYPE)115,115,105
  105 Z(I)=T1
      GO TO 125
  115 Z1(I)=T1
  125 CONTINUE
      IF(ITYPE)135,135,300
  135 CALL MV(BB,Z1,Z)
  300 RETURN
      END
      SUBROUTINE PAXES(ACODE,ITYPE)
C     ***** ITYPE .LT.0 FOR COVARIANCE MATRIX IN Q *****
C     ***** ITYPE .GT.0 FOR ELLIPSOID QUADRATIC FORM IN Q *****
C     ***** XABSF(ITYPE)=1 BASED ON TRICLINIC COORDINATE SYSTEM *****
C     ***** =2 OR 3 FOR WORKING OR REFERENCE CARTESIAN SYSTEMS *****
C     ***** CONTRAVARIANT EIGENVECTORS FOR Q IN COLUMNS OF PAC *****
      DIMENSION W(3,3),X(3)
C     ***** CHECK ATOM CODE *****

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      IT=IABS(ITYPE)-1
      KS=AMOD(ACODE,100.)
      IF(NSYM-KS)105,115,115
  105 NG=4
      GO TO 300
  115 II=ACODE/100000.
      IF(NATOM-II)125,130,130
  125 NG=5
      GO TO 300
  130 IF(II)125,125,135
C     ***** CRYSTALLOGRAPHIC SYMMETRY ROTATION *****
  135 CALL TMM(PA(1,1,II),FS(1,1,KS),PAT)
      IF(IT-1)160,145,155
C     ***** TRANSFORM TO CARTESIAN SYSTEMS *****
  145 CALL TMM(PAT,AAWRK,PAC)
      GO TO 175
  155 CALL TMM(PAT,AAREV,PAC)
      GO TO 175
  160 IF(ITYPE)162,155,170
C     ***** TRANSFORM TO TRICLINIC SYSTEM *****
  162 DO 165 J=1,9
  165 PAC(J,1)=PAT(J,1)
      GO TO 175
  170 CALL MM(AA,PAT,PAC)
C     ***** FORM DIAGONAL MATRIX OR ITS INVERSE *****
  175 DO 205 J=1,3
      T1=EV(J,II)
      IF(ITYPE)195,195,185
  185 X(J)=1./(T1*T1)
      GO TO 205
  195 X(J)=T1*T1
  205 RMS(J)=T1
C     ***** FORM QUADRATIC FORM *****
      DO 245 I=1,3
      DO 245 J=I,3
      T1=0.0
      DO 225 K=1,3
  225 T1=T1+PAC(I,K)*PAC(J,K)*X(K)
      Q(J,I)=T1
  245 Q(I,J)=T1
  300 RETURN
      END
      SUBROUTINE PLOTSZ(IDUM,JDUM,KDUM)
C     THIS PACKAGE WILL CREATE A FILE CONTAINING THE ARGUMENTS OF
C     ALL CALL TO SUBROUTINE PLOT FROM ORTEP. THESE FILES MAY
C     THEN BE READ BY PROGRAMS FOR PRODUCING OUTPUT FOR VERSATEC,
C     PRINTRONIX, OR RAMTEK.
      DATA LUN/25/
C#      OPEN( UNIT=LUN, FORM='UNFORMATTED', RECORDTYPE='FIXED',
C#     3     RECORDSIZE=3, STATUS = 'NEW')
      RETURN
      ENTRY OTPLOT(X,Y,ISW)
      WRITE(LUN)X,Y,ISW
      IF(ISW.EQ.999) CLOSE(LUN)
      RETURN
      END
      SUBROUTINE PLTXY(X,Y)
C     ***** PLOT COORD. AND CLOSEST EDGE AFTER PROJECTION *****
      DIMENSION X(3),Y(2)
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *
      T4=1.
      T1=1.
      IF(VIEW)125,125,110
  110 T4=VIEW-X(3)
      IF(T4)115,115,120
  115 Y(1)=-99.
      Y(2)=-99.
      GO TO 130
  120 T1=VIEW/T4
  125 Y(1)=X(1)*T1+XO(1)
      Y(2)=X(2)*T1+XO(2)
      T1=XLNG(1)-ABS(Y(1)*2.-XLNG(1))
      T2=XLNG(2)-ABS(Y(2)*2.-XLNG(2))
      EDGE=AMIN1(T1,T2)*.5
      IF(T4-VIEW*.5)130,300,300
  130 EDGE=-99.
  300 RETURN
      END
      SUBROUTINE PRIME
C     ****GENERAL INITIALIZATION OF PRIME PARAMETERS****
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *
      BRDR=0.5
C     ****CALCULATE CONSTANTS****
      DO 2999 I=1,5
 2999 CONT(I)=SQRT(1./(2.*(1.+COS(3.141593   /2.**I))))
      DISP=.005
      FORE=.866
   20 IN=5
      ITILT=0
      LATM=0
      NCD=0
      NG=0
      NOUT=6
C?      NSR=2
C?      OPEN( UNIT=NSR, STATUS='SCRATCH', FORM='UNFORMATTED')
      NTEMPT =27
      RES(1)=1.25
      RES(2)=.5
      RES(3)=.2
      SCAL1=1.0
      SCAL2=1.54
      SCL=1.54
      DO 3000 I=1,3
      SYMB(I,I)=1.
      SYMB(I+1,1)=0.
 3000 SYMB(I+5,1)=0.
      TAPER=.375
      THETA=0.0
      VIEW=0.0
      XLNG(1)=30.
      XLNG(2)=11.0
      XO(1)=15.
      XO(2)=5.5
      DO 3001 J=1,3
 3001 ORGN(J) = 0.0
C     ***** INITIATE OVERLAP ROUTINES *****
      CALL LAP500(0)
      RETURN
      END
      SUBROUTINE PROJ(D,DP,X,XO,VIEW,I1,I2,I3)
C     ***** 3D CARTESIAN TO 2D PLOTTER COORDINATES *****
      DIMENSION D(3,129),DP(2,130),X(3),XO(3)
      T3=VIEW-X(3)
      DO 145 I=I1,I2,I3
      T1=D(1,I)+X(1)
      T2=D(2,I)+X(2)
      IF(VIEW)135,135,120
  120 T4=VIEW/(T3-D(3,I))
      T1=T1*T4
      T2=T2*T4
  135 DP(1,I)=T1+XO(1)
  145 DP(2,I)=T2+XO(2)
      RETURN
      END
      SUBROUTINE RADIAL(ND)
C     ***** GENERATE ELLIPSE FROM TWO CONJUGATE VECTORS *****
C     ***** ORTHONORMAL VECTORS PRODUCE 8-128 SPOKED CIRCLE *****
C     ***** ND DENOTES NUMBER OF SUBDIVISIONS (1 TO 5) *****

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      DO 115 J=1,3
      T1=DA(J,1)
      D(J,1)=T1
      D(J,129)=T1
      D(J,65)=-T1
      T1=DA(J,2)
      D(J,33)=T1
  115 D(J,97)=-T1
      DO 135 K=1,ND
      T1=CONT(K)
      KDEL=2**(6-K)
      KDEL1=KDEL+1
      KDEL2=KDEL/2
      DO 135 L=KDEL1,65,KDEL
      J=L-KDEL
      M=L-KDEL2
      DO 135 N=1,3
      T2=(D(N,L)+D(N,J))*T1
      D(N,M)=T2
  135 D(N,M+64)=-T2
      RETURN
      END

      SUBROUTINE SCRIBE(Y,NPEN,LTNO, ICOLOR)
      INTEGER OLDCOL
      CHARACTER*80 COMMNT
      DIMENSION Y(2),YO(2)
      DATA OLDCOL / -1 /
C     ***** SUBROUTINE WHICH LINKS WITH THE PLOTTER-SPECIFIC SUBROUTINES
      IF ( ICOLOR .NE. OLDCOL ) THEN
         OLDCOL = ICOLOR
         CALL OPLOT(FLOAT(OLDCOL), 0.0, 99 )
      ENDIF
      IF(NPEN-3)210,205,205
C     ***** KEEP TRACK OF COORDINATES FOR LAST PEN-UP LOCATION *****
  205 YO(1)=Y(1)
      YO(2)=Y(2)
      NPO=0
      RETURN
  210 IF(LTNO)230,245,215
C     ***** CALL MECHANICAL PLOTTER PLOTTING SUBROUTINE *****
  215 IF(NPO)225,220,225
  220 CALL OPLOT( YO(1), YO(2) ,2)
  225 CALL OPLOT( Y(1), Y(2), 3)
      GO TO 245
  230 CONTINUE
  245 NPO=1
      RETURN
      END
      SUBROUTINE SEARC
      DIMENSION NW(6),DX(3),S(2,200),U(3),V(3),W(2,4),WW(2,3),X(4),Y(3)
      DIMENSION Z(3),I0(2),I02(2),I2(2),J2(2)
C
      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *
      CHARACTER HEADER*80, CHEM(NUMATM)*6, COMMNT*80
      COMMON /STRING/ HEADER, CHEM
C
      LOGICAL DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      COMMON /DEBCOM/ DEBUG, DEBUGL, DEBUGN, DEBUGO, DEBUGP, DEBUGI
      EQUIVALENCE(NW(1),LL),(NW(2),LU),(NW(3),ML)
      EQUIVALENCE(NW(4),MU),(NW(5),NL),(NW(6),NU)
C     ***** OBTAIN PROBLEM PARAMETERS *****
      IF (DEBUGO) THEN
         WRITE ( COMMNT, 20)
         CALL DEBUGR( COMMNT(1:78) )
         WRITE ( COMMNT, 21)
         CALL DEBUGR( COMMNT(1:78) )
      ENDIF
   20 FORMAT(1H0 2X, 'FROM ATOMS    TO ATOMS     WITH RADIUS OR',
     1', IF A BOX, WITH SEMIDIMENSIONS')
   21 FORMAT( 11X,'CODE   (MIN  MAX)   (MIN  MAX)',
     17X,'A',8X,'B',8X,'C')
      IF(AIN(1)-10000.)100,100,101
  100 ITOM1=AIN(1)
      SYITOM=55501.
      GO TO 103
  101 ITOM1=AIN(1)/100000.
      SYITOM=AMOD(AIN(1),100000.)
  102 IF(ABS(AIN(2))-10000.)103,103,104
  103 ITOM2 = ABS(AIN(2))
      SYITO2=SYITOM
      GO TO 105
  104 ITOM2=ABS(AIN(2))/100000.
      SYITO2=AMOD(ABS(AIN(2)),100000.)
  105 ITAR1=AIN(3)
      IF(ITAR1)108,108,110
  108 ITAR1=1
  110 ITAR2=AIN(4)
      DMAX=AIN(5)
      IF(DMAX)115,115,120
  115 DMAX=4.
      AIN(5)=DMAX
  120 DMX=DMAX*DMAX
      TEM=.01
      KFUN=NJ*100+MOD(NJ2,10)
      K=NJ*100+NJ2
      I0(1)=SYITOM/1000.
      I0(2)=AMOD(SYITOM,1000.)
      I02(1)=SYITO2/1000.
      I02(2)=AMOD(SYITO2,1000.)
      LATOM=LATM
  121 FORMAT( 2X,2I3,I2,I3,I4,I2,I3,2I4,18X,3F9.3 )
      IF (DEBUGO) THEN
         WRITE ( COMMNT, 121)K,ITOM1,I0(1),I0(2),ITOM2,I02(1),I02(2),
     .          ITAR1,ITAR2,(AIN(J),J=5,7)
         CALL DEBUGR( COMMNT(1:78) )
      ENDIF
  124 FORMAT(1H 15X,2I5,I8,I5,2F9.3)
      IF(NCD)130,130,125
  125 IF (DEBUGO) THEN
         WRITE ( COMMNT,124) ((KD(J,I),J=1,4),(CD(J,I),J=1,2) ,I=1,NCD)
         CALL DEBUGR( COMMNT(1:78) )
      ENDIF
  130 DO 135 J=1,4
      W(1,J)=99.
  135 W(2,J)=-99.
      DO 155 I=ITAR1,ITAR2
      T1=(I)*100000.
      CALL ATOM(T1,X)
      IF(NG)140,145,140
  140 CALL ERPNT(T1,KFUN)
      GO TO 600
  145 X(4)=X(1)-X(2)
      DO 155 J=1,4
      TEM=X(J)
      IF(W(2,J)-TEM)148,150,150
  148 W(2,J)=TEM
  150 IF(TEM-W(1,J))152,155,155
  152 W(1,J)=TEM
  155 CONTINUE
      KFUN2=MOD(KFUN,10)
      GO TO (165,165,160,156,165,165),KFUN2
C     ***** FIND PARALLELEPIPED WHICH ENCLOSES TRICLINIC BOX *****
  156 DO 158 J=1,3
  158 DX(J)=AIN(J+4)
      GO TO 170
C     ***** FIND PARALLELEPIPED WHICH ENCLOSES RECTANGULAR BOX *****
  160 DO 162 J=1,3
      DX(J)=0.
      DO 162 I=1,3
  162 DX(J)=DX(J)+ABS(REFV(J,I)*AIN(I+4))
      GO TO 170
C     ***** FIND PARALLELEPIPED WHICH ENCLOSES DMAX SPHERE *****
  165 T1=1.-A(4)*A(4)-A(5)*A(5)-A(6)*A(6)+2.*A(4)*A(5)*A(6)
      DO 168 J=1,3
  168 DX(J)=SQRT((1.-A(J+3)**2)/T1)*DMAX/A(J)
C     ***** START SEARCH AROUND REFERENCE ATOMS *****
  170 LIST=0
      LAST=0
      M1=ITOM1
      N1=ITOM2
      IF(KFUN2-5)186,172,172
C     ***** CONVOLUTE AND REITERATIVE CONVOLUTE INSTRUCTIONS *****
  172 IF(LATM)174,174,176
C     ***** FAULT, NO ENTRIES IN ATOMS LIST *****
  174 NG=12
      CALL ERPNT(0.,KFUN)
      GO TO 600
C     ***** CHECK FOR REFERENCE ATOMS IN ATOMS LIST *****
  176 IF(LATM-LAST)600,600,177
  177 LIST=LAST
      LAST=LATM
  178 LIST=LIST+1
      IF(LAST-LIST)505,180,180
  180 T1=ATOMS(1,LIST)
      ITOM=T1/100000.
      IF(ITOM-ITOM1)178,184,182
  182 IF(ITOM2-ITOM)178,184,184
  184 SYITOM=AMOD(T1,100000.)
      SYITO2=SYITOM
      M1=ITOM
      N1=ITOM
C     ***** SET INITIAL RUN PARAMETERS *****
  186 M2=AMOD(SYITOM,100.)
      M5=AMOD(SYITOM/100.,1000.)
      M3=M5/100
      M4=MOD(M5/10,10)
      M5=MOD(M5,10)
C     ***** SET TERMINAL RUN PARAMETERS *****
      N2=AMOD(SYITO2,100.)
      N5=AMOD(SYITO2/100.,1000.)
      N3=N5/100
      N4=MOD(N5/10,10)
      N5=MOD(N5,10)
C     ***** START SEARCH AROUND REFERENCE ATOMS *****
      DO 500 L5=M5,N5
      DO 500 L4=M4,N4
      DO 500 L3=M3,N3
      DO 500 L2=M2,N2
      DO 500 ITOM=M1,N1
      T1=(L2)+100.*(L5+10.*(L4+10.*(L3+10.*ITOM)))
      CALL ATOM(T1,Y)
      IF(NG)188,190,188
  188 CALL ERPNT(T1,KFUN)
      GO TO 500
C     ***** K=SYMMETRY EQUIVALENT POSITION *****
  190 NUM=0
      DO 400 K=1,NSYM
C     ***** SUBTRACT SYMMETRY TRANSLATION FROM REFERENCE ATOM *****
      DO 192 J=1,3
  192 U(J)=Y(J)-TS(J,K)
C     ***** DETERMINE LIMITING CELLS TO BE SEARCHED *****
C     ***** FIRST,MOVE THE BOX THROUGH THE SYMMETRY OPERATION *****
      DO 200 J=1,3
      DO 200 L=1,2
      WW(L,J)=0.0
      DO 200 I=1,3
      TEM=FS(I,J,K)
      IF(TEM)194,200,196
  194 N=MOD(L,2)+1
      GO TO 198
  196 N=L
  198 WW(L,J)=WW(L,J)+W(N,I)*TEM
  200 CONTINUE
C     ***** CHECK FOR MIXED INDEX TRANSFORMATION *****
      DO 215 J=1,2
      TEM=FS(1,J,K)
      IF(TEM+FS(2,J,K))215,201,215
  201 IF(TEM)203,215,207
  203 WW(1,J)=W(2,4)*TEM
      WW(2,J)=W(1,4)*TEM
      GO TO 215
  207 WW(1,J)=W(1,4)*TEM
      WW(2,J)=W(2,4)*TEM
  215 CONTINUE
C     ***** MOVE 4 CELLS AWAY THEN MOVE BACK UNTIL PARALLELEPIPED AROUND
C         REF ATOM AND BOX AROUND TRANSFORMED ASYM UNIT INTERSECT *****
      N=0
      DO 235 J=1,3
      DO 225 I=1,2
      N=N+1
      TT=(U(J)-WW(I,J))*(I*2-3)-DX(J)
      TEM=5.0
  221 TEM=TEM-1.0
      IF(TEM+TT)225,225,221
  225 NW(N)=TEM*(I*2-3)+5.
C     ***** IF NO POSSIBILITY OF A HIT, GO TO NEXT SYMMETRY OPER *****
      IF(NW(N)-NW(N-1))400,235,235
  235 CONTINUE
      LL=NW(1)
      LU=NW(2)
      ML=NW(3)
      MU=NW(4)
      NL=NW(5)
      NU=NW(6)
C     ***** L CELL TRANSLATIONS IN X *****
      DO 395 L=LL,LU
      V(1)=U(1)+(L-5)
C     ***** M CELL TRANSLATIONS IN Y *****
      DO 395 M=ML,MU
      V(2)=U(2)+(M-5)
C     ***** N CELL TRANSLATIONS IN Z *****
      DO 395 NN=NL,NU
      V(3)=U(3)+(NN-5)
C     ***** I = TARGET ATOM *****
      DO 395 I=ITAR1,ITAR2
      DO 250 J=1,3
      TEM=0.0
      DO 245 II=1,3
  245 TEM=TEM+FS(II,J,K)*P(II,I)
C     ***** SEE IF WITHIN PARALLELEPIPED*****
      TEM=TEM-V(J)
      IF(DX(J)-ABS(TEM))395,250,250
  250 X(J)=TEM
      GO TO (255,255,252,277,255,255),KFUN2
C     ***** SEE IF WITHIN MODEL BOX *****
  252 CALL VM(X,AAREV,V1(2))
      DO 253 J=2,4
      IF(AIN(J+3)-ABS(V1(J)))395,253,253
  253 CONTINUE
      GO TO 277
C     ***** SEE IF WITHIN SPHERE *****
  255 DSQ=VMV(X,AA,X)
      IF(DMX-DSQ)395,256,256
  256 IF(DSQ-.0001)258,260,260
  258 IF(KFUN-402)395,260,260
  260 TEM=SQRT(DSQ)
      IF(AIN(8))265,265,261
C     *****SELECT ONLY FIRST ASYMMETRIC UNIT ENCOUNTERED *****
  261 IF(LATM)265,265,262
  262 ZMIN=(I)*100000.
      ZMAX=ZMIN+100000.
      DO 264 J=1,LATM
      ZSTO=ATOMS(1,J)
      IF(ZSTO-ZMIN)264,263,263
  263 IF(ZMAX-ZSTO)264,264,395
  264 CONTINUE
C     ***** SELECT VECTORS ACCORDING TO CODES IF ANY *****
  265 IF(NCD)277,277,268
  268 DO 275 J=1,NCD
  269 IF(ITOM-KD(1,J))275,270,270
  270 IF(KD(2,J)-ITOM)275,271,271
  271 IF(I-KD(3,J))   275,272,272
  272 IF(KD(4,J)-I)   275,273,273
  273 IF(TEM-CD(1,J)) 275,274,274
  274 IF(CD(2,J)-TEM) 275,277,277
  275 CONTINUE
      GO TO 395
  277 V1(1) = 100000.*(I)+(1110-L*100-M*10-NN)*100.+K
      IF(KFUN-402)278,325,325
C     ***** DETERMINE CORRECT POSITION IN SORTED VECTOR TABLE *****
  278 IF(NUM)317,317,279
  279 DO 315 II=1,NUM
      TT=S(2,II)-TEM
      IF(ABS(TT)-0.0001)297,297,281
  281 IF(TT)315,297,283
C     ***** MOVE LONGER VECTORS TOWARD END OF TABLE *****
  283 IF(200-NUM)287,287,289
  287 NUM=199
  289 IJ=NUM
      DO 295 J=II,NUM
      S(1,IJ+1)=S(1,IJ)
      S(2,IJ+1)=S(2,IJ)
  295 IJ=IJ-1
      GO TO 319
C     ***** CHECK FOR DUPLICATE VECTORS IF DISTANCES ARE EQUAL *****
  297 CALL ATOM(S(1,II),Z)
      DO 305 J=1,3
      IF(ABS(X(J)+Y(J)-Z(J))-0.0001)305,305,315
  305 CONTINUE
      GO TO 395
  315 CONTINUE
      IF(200-NUM)320,320,317
C     ***** STORE THE RESULT IN VECTOR TABLE *****
  317 II=NUM+1
  319 NUM=NUM+1
      S(1,II)=V1(1)
      S(2,II)=TEM
  320 IF(KFUN-106)395,325,325
C     ***** STORE RESULT IN ATOMS TABLE *****
  325 DO 330 J=1,3
  330 V1(J+1)=X(J)+Y(J)
      CALL STOR
  395 CONTINUE
  400 CONTINUE
C     ***** PRINT OUT DISTANCES *****
  421 FORMAT(1H010X,20HVECTORS FROM ATOM  (I3,1H,I2,I3,1H)6X,8HTO ATOMS
     1I4,8H THROUGHI4)
      I0(1)= AMOD(T1,100000.)/1000.
      I0(2)=AMOD(T1,1000.)
      IF (DEBUGO) THEN
         CALL OPLOT( 0.0, 0.0, 8)
         WRITE (NOUT,421)ITOM,I0,ITAR1,ITAR2
         CALL OPLOT( 0.0, 0.0, 9)
      ENDIF
      IF(NUM)500,500,423
  423 DO 435 I=1,NUM
      T2=S(1,I)
      I1=T2/100000.
      I2(1)=AMOD(T2,100000.)/1000.
      I2(2)=AMOD(T2,1000.)
      CALL ATOM(T2,Z)
      IF(I-1)432,432,434
  427 FORMAT(1H 13X,2(A6,1X),39X,1H(I3,1H,I2,I3,1H)3F7.4,7X,3HD =F6.3)
  429 FORMAT(1H 13X,2(A6,1X),2(3H  (I3,1H,I2,I3,1H)3F7.4,3X),4X,3HD =
     1,F6.3)
  432 IF (DEBUGO) THEN
         CALL OPLOT( 0.0, 0.0, 8)
         WRITE (NOUT,429)CHEM(ITOM),CHEM(I1),ITOM,I0,(Y(J),J=1,3
     .           ),I1,I2,(Z(J),J=1,3),S(2,I)
         CALL OPLOT( 0.0, 0.0, 9)
      ENDIF
      GO TO 435
  434 IF (DEBUGO) THEN
         CALL OPLOT( 0.0, 0.0, 8)
         WRITE (NOUT,427)CHEM(ITOM),CHEM(I1),I1,I2,(Z(J),J=1,3),
     .           S(2,I)
         CALL OPLOT( 0.0, 0.0, 9)
      ENDIF
  435 CONTINUE
C     ***** CALCULATE ANGLES ABOUT REF ATOM IF CODE IS 102 *****
  437 IF(KFUN-102)500,451,500
  441 FORMAT(1H010X,18HANGLES AROUND ATOMI5)
  451 IF (DEBUGO) THEN
         CALL OPLOT( 0.0, 0.0, 8)
         WRITE (NOUT,441)ITOM
         CALL OPLOT( 0.0, 0.0, 9)
      ENDIF
      L=NUM-1
      IF(L)500,500,457
  457 DO 465 I=1,L
      T2=S(1,I)
      T3=S(2,I)
      I1=T2/100000.
      I2(1)=AMOD(T2,100000.)/1000.
      I2(2)=AMOD(T2,1000.)
      CALL ATOM(T2,X)
      CALL DIFV(X,Y,U)
      CALL MV(AA,U,V2)
      M=I+1
      DO 465 J=M,NUM
      T4=S(1,J)
      J1=T4/100000.
      J2(1)=AMOD(T4,100000.)/1000.
      J2(2)=AMOD(T4,1000.)
      CALL ATOM(T4,Z)
      CALL DIFV(Z,Y,V)
      TMPDMS = VV( V, V2) / ( T3 * S( 2, J) )
      IF ( ABS(TMPDMS) .GT. 1.0 ) THEN
         CALL DEBUGR( ' SEARC IN ORTEP, ACOS ERROR - AVOIDED.')
         TMPDMS = SIGN( 1.0, TMPDMS )
      ENDIF
      F = ACOS( TMPDMS)
C      F=ACOS(VV(V,V2)/(T3*S(2,J)))
      CALL DIFV(X,Z,V3)
      F1=SQRT(VMV(V3,AA,V3))
  460 FORMAT(1H 13X,3(A6,1X),7X,3(2H (I3,1H,I2,I3,1H)),12X,3HD =F6.3,7X
     1,3HA =F6.2)
  465 CONTINUE
      IF (DEBUGO) THEN
         CALL OPLOT( 0.0, 0.0, 8)
         WRITE (NOUT,460)CHEM(I1),CHEM(ITOM),CHEM(J1),I1,I2,ITOM
     .           ,I0,J1,J2,F1,F
         CALL OPLOT( 0.0, 0.0, 9)
      ENDIF
  495 CONTINUE
  500 CONTINUE
      IF(LAST-LIST)505,505,178
  505 IF(KFUN2-6)600,176,600
  600 IF(KFUN-106)610,605,610
  605 LATM=LATOM
  610 RETURN
      END
      SUBROUTINE SPARE(INST)
C     ***** INSTRUCTION 1201 PRODUCES PUNCHED CARDS WITH POSIT.+THERMAL
C     ***** PARAMETERS FOR ATOMS PLACED IN THE ATOMS LIST BY THE 400
C     ***** SERIES INSTRUCTIONS

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

C     ***** FIRST FIELD ON INSTRUCTION CARD = NO. OF POSITION PARAMETER
C     ***** CARDS TO BE READ BY THIS SUBROUTINE (STANDARD ORFLS FORMAT)
C     ***** THE SEQUENCE OF THESE CARDS NEED NOT BE IDENTICAL WITH
C     ***** THE ORIGINAL INPUT ATOM DECK SEQUENCE
      IF(INST-1201)100,105,100
  100 NG = 9
      RETURN
  105 IF(AIN(1))100,100,115
  115 LAST = AIN(1)
      DO 150 L=1,LAST
      READ (IN,116) CHM,V2(1),V2(2),(V1(J),J=1,4)
  116 FORMAT(A6,3X,6F9.6)
      CALL OPLOT( 0.0, 0.0, 8)
      WRITE (NOUT,117) CHM,V2(1),V2(2),(V1(J),J=1,4)
      CALL OPLOT( 0.0, 0.0, 9)
  117 FORMAT(1H0,24X,25HOLD POSITION PARAMETERS  ,A6,4X,6F10.6)
C     ***** FIND ATOM NUMBER *****
      DO 125 IA=1,NATOM
      DO 120 J=1,3
      IF(ABS(V1(J)-P(J,IA))-0.001)120,120,125
  120 CONTINUE
      IX = IA
      GO TO 130
  125 CONTINUE
      GO TO 150
C     ***** SEARCH ATOMS LIST FOR ENTRIES WITH ATOM NUMBER IX *****
  130 IF(LATM)150,150,131
  131 DO 145 IB=1,LATM
      T1 = ATOMS(1,IB)
      IF(IX-INT(T1/100000.))145,135,145
  135 CALL ATOM(T1,V1)
      CALL OPLOT( 0.0, 0.0, 8)
      WRITE (NOUT,136) T1,CHM,V2(1),V2(2),(V1(J),J=1,4)
      CALL OPLOT( 0.0, 0.0, 9)
  136 FORMAT(1H0,9X,F10.0,5X,25HNEW POSITION PARAMETERS  ,A6,4X,6F10.6)
      I1 = T1
      CALL PAXES(T1,-1)
      DO 140 I=1,3
      DO 140 J=1,3
  140 Q(J,I) = Q(J,I) /.0506605918
      CALL OPLOT( 0.0, 0.0, 8)
      WRITE (NOUT,141) Q(1,1),Q(2,2),Q(3,3),Q(1,2),Q(1,3),Q(2,3)
      CALL OPLOT( 0.0, 0.0, 9)
  141 FORMAT(1H ,24X,25HNEW THERMAL PARAMETERS   ,6F10.6)
  145 CONTINUE
  150 CONTINUE
  300 RETURN
      END
      SUBROUTINE STOR
C     ***** STORE IN OR REMOVE FROM ATOMS ARRAY *****

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      IF(LATM)481,481,450
  450 IF(500-LATM)455,455,460
  455 IF(NJ2-10)490,490,460
  460 L=LATM
C     ***** CHECK FOR POSITIONAL DUPLICATION *****
      DO 480 K=1,L
      DO 465 J=2,4
      IF(ABS(V1(J)-ATOMS(J,K))-0.001)465,465,480
  465 CONTINUE
      IF(NJ2-10)490,490,470
C     ***** ATOM REMOVAL BY TABLE PUSHDOWN *****
  470 LATM=LATM-1
      DO 475 I=K,LATM
      DO 475 J=1,4
  475 ATOMS(J,I)=ATOMS(J,I+1)
      GO TO 490
  480 CONTINUE
  481 IF(NJ2-10)482,490,490
C     ***** STORE ATOM *****
  482 IF(499-LATM)490,483,485
  483 NG=16
      CALL ERPNT (V1(1),400)
  485 LATM=LATM+1
      DO 486 J=1,4
  486 ATOMS(J,LATM)=V1(J)
  490 RETURN
      END
      SUBROUTINE TMM(X,Y,Z)
C     ***** TRANSPOSE(TRANSPOSE(X)*Y)=Z *****
C     ***** X,Y,Z ARE 3X3 MATRICES *****
      DIMENSION X(3,3),Y(3,3),Z(3,3)
      DO 115 I=1,3
      DO 115 K=1,3
  115 Z(K,I)=X(1,I)*Y(1,K)+X(2,I)*Y(2,K)+X(3,I)*Y(3,K)
      RETURN
      END

      SUBROUTINE UNIT(X,Z,ITYPE)
      DIMENSION X(3),Y(3),Z(3)

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      Y(1)=X(1)
      Y(2)=X(2)
      Y(3)=X(3)
      IF(ITYPE)125,125,105
  105 T1=SQRT(Y(1)*Y(1)+Y(2)*Y(2)+Y(3)*Y(3))
      GO TO 145
  125 T1=SQRT(Y(1)*(Y(1)*AA(1,1)+Y(2)*(AA(1,2)+AA(2,1))+Y(3)*(AA(1,3)+A
     1A(3,1)))+Y(2)*(Y(2)*AA(2,2)+Y(3)*(AA(2,3)+AA(3,2)))+Y(3)*Y(3)*AA(3
     2,3))
  145 IF(T1)155,155,175
  155 NG=5
      GO TO 300
  175 Z(1)=Y(1)/T1
      Z(2)=Y(2)/T1
      Z(3)=Y(3)/T1
  300 RETURN
      END
      SUBROUTINE VM(X,Y,Z)
C     TRANSPOSED VECTOR TIMES MATRIX
C     Z(3)=X(3)*Y(3,3)
      DIMENSION X(3),Y(3,3),Z(3)
      DO 115 J=1,3
      Z(J)=0.0
      DO 115 I=1,3
  115 Z(J)=Z(J)+X(I)*Y(I,J)
      RETURN
      END
      FUNCTION VMV(X1,Q,X2)
C     TRANSPOSED VECTOR * MATRIX * VECTOR
C     VMV=X1(3)*Q(3,3)*X2(3)    TO  EVALUATE QUADRATIC OR BILINEAR FORM
      DIMENSION X1(3),Q(3,3),X2(3)
      T1=0.
      DO 10 J=1,3
   10 T1=T1+X1(J)*(X2(1)*Q(J,1)+X2(2)*Q(J,2)+X2(3)*Q(J,3))
      VMV=T1
      RETURN
      END
      FUNCTION VV(X,Y)
C     TRANSPOSED VECTOR * VECTOR
C     VV=X(3)*Y(3)
      DIMENSION X(3),Y(3)
      VV=X(1)*Y(1)+X(2)*Y(2)+X(3)*Y(3)
      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 XYZ(QA,X,ITYPE)
C     ***** ITYPE .GT.0 CART. COORD. FROM ATOM CODE WORD *****
C     ***** XABSF(ITYPE) .LE.2 FOR WORKING SYSTEM *****
C     ***** XABSF(ITYPE) .GT.2 FOR REFERENCE SYSTEM *****
C     ***** ITYPE .LE.0 USES TRICLINIC COORD. XT *****
      DIMENSION X(3)

      INCLUDE 'SIZES'
C * BLOCK COMMON *
      DIMENSION A(9),AA(3,3),AAREV(3,3),AAWRK(3,3),AID(3,3),AIN(140)
      DIMENSION ATOMS(4,NUMATM),BB(3,3),CD(8,10),CONT(5),D(3,130)
      DIMENSION DA(3,3),DP(2,130),EV(3,NUMATM)
      DIMENSION FS(3,3,48),KD(5,10),ORGN(3)
      DIMENSION P(3,NUMATM),PA(3,3,NUMATM)
      DIMENSION PAC(3,5),PAT(3,3),Q(3,3),REFV(3,3)
      DIMENSION RES(4),RMS(5),SYMB(3,3),TITLE(18),TS(3,48)
      DIMENSION VT(3,4),V1(4),V2(3),V3(3),V4(3),V5(3),V6(3),WRKV(3,3)
      DIMENSION XLNG(3),XO(3),XT(3)

      COMMON /ORNL/
     . NG,A,AA,AAREV,AAWRK,AID,AIN,ATOMS,BB,BRDR,CD,CONT,D ,
     . DA,DP,DISP,EDGE,EV,FORE,FS,IN,ITILT,KD,LATM,LTNO,NATOM,NCD ,
     . NJ,NJ2,NOUT,NSYM,ORGN,P,PA,PAC,PAT,Q,REFV,RES,RMS,SCAL1,
     . SCAL2,SCL,SYMB,TAPER,THETA,TITLE,TS,VIEW,VT,V1,V2 ,
     . V3,V4,V5,V6,WRKV,XLNG,XO,XT,NTEMPT

C *     END OF COMMON   *

      IT=IABS(ITYPE)-2
      NG1=NG
      NG=0
      IF(ITYPE)10,10,5
    5 CALL ATOM(QA,XT)
      IF(NG)30,10,30
   10 T1=0.
      DO 15 J=1,3
      T2=XT(J)-ORGN(J)
      V1(J)=T2
   15 T1=T1+ABS(T2)
      IF(T1-.0001)20,20,40
   20 NG=NG1
   30 DO 35 J=1,3
   35 X(J)=0.
      GO TO 300
   40 IF(IT)45,45,60
C     ***** RELATIVE TO WORKING SYSTEM *****
   45 DO 55 I=1,3
      T1=0.
      DO 50 J=1,3
   50 T1=T1+V1(J)*AAWRK(J,I)
   55 X(I)=T1*SCAL1
      GO TO 300
C     ***** RELATIVE TO REFERENCE SYSTEM *****
   60 DO 70 I=1,3
      T1=0.
      DO 65 J=1,3
   65 T1=T1+V1(J)*AAREV(J,I)
   70 X(I)=T1*SCAL1
  300 RETURN
      END
