C
C       ***********************************************************
C       ********************  PSI/88 - PART 2  ********************
C       ***********************************************************
C
C       Version 1.0  Any questions to the author should specify
C                    the version being used.
C
      PROGRAM PSICON
C
C       THIS ROUTINE "HACKED OUT OF" PSI77, PART 1, WRITTEN BY:
C
C       WILLIAM L. JORGENSEN
C       Yale University, DEPARTMENT OF CHEMISTRY
C       New Haven, CT 06511, USA
C       PHONE 203-432-6278
C
C       HACKING OUT AND MODIFICATION TO ALLOW HIDDEN LINE ELIMINATION
C       OF INDENTATIONS AND DOUGHNUTS.  The modifications are rather
C       minor, so the code is mostly unchanged from PSI/77.
C       DAN SEVERANCE, PURDUE
C
C Redistribution and use in source and binary forms are permitted
C provided that the above paragraphs and this one are duplicated in
C all such forms and that any documentation, advertising materials,
C and other materials related to such distribution and use acknowledge
C that the software was developed by William Jorgensen at Purdue University
C The name of the University or William Jorgensen may not be used to endorse
C or promote products derived from this software without specific prior
C written permission.  The author is now at Yale University.
C THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
C IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
C WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
C
      COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(850),Y(850),Z(850),IDASH,
     *   SCALE,PERZ,KA,KO,ZO
      COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850)
      COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3
C
C       THIS ROUTINE READS THE 3-D GRID OF ORBITAL VALUES GENERATED
C       BY THE PSI1 PROGRAM AND STORED IN DISK FILE 17. (UNFORMATTED)
C       ALL DATA NEEDED BY THE PROGRAM IS ALSO STORED IN FILE 17.
C
C       THE SAME INPUT DECK IS USED AS WAS SPECIFIED FOR PSI1.
C       DIFFERENT CONTOUR LEVELS MAY BE SELECTED AND THIS PROGRAM RE-RUN
C       WITHOUT REXECUTION OF PSI1.
C
C       CLOSED CONTOURS AT THE SPECIFIED LEVEL(S)
C       ARE DRAWN THROUGH THE DENSITY OR ORBITAL VALUE
C       ARRAYS AND THE POINTS(IN ANGSTROMS) FOR EACH CONTOUR ARE STORED
C       AS OUTPUT IN DISK FILE 22. THE TOTAL
C       NUMBER OF POINTS AND THE NUMBER OF CURVES ALONG
C       WITH THE NUMBER OF POINTS IN EACH CURVE AND AN
C       INDICATOR FOR WHETHER THE CONTOUR LEVEL IS POSITIVE
C       OR NEGATIVE ARE STORED IN THE OUTPUT DISK FILE
C       23. FOR PROPER HIDDEN LINE ELIMINATION
C       IN THE 3-D DRAWING PROGRAM, IT IS ESSENTIAL THAT ALL CONTOURS
C       WHICH ARE CONSTRUCTED BE CLOSED. IF THIS CONDITION IS NOT MET,
C       AN ERROR MESSAGE IS PRINTED.
C       CARD INPUT IS DESCRIBED IN SUBROUTINE THREED.
C
C       NOT ALL OF THE VARIABLES IN COMMON/RPLOT ARE
C       ESSENTIAL TO THIS PROGRAM, BUT HAVE BEEN RETAINED
C       FOR COMPATIBILITY WITH THE HIDDEN LINE PLOTTING
C       PROGRAM IN CASE THE TWO ARE MERGED
C
      IRD = 5
      ILST = 6
      IDSK1 = 22
      IDSK2 = 23
      IDSK3 = 24
      THE = 0.0E+0
      GAM = 0.0E+0
      PHI = 0.0E+0
      NC = 0
      NCURV = 0
C
C       MAKE 3-D CONTOUR MAP AND STORE ON DISK
C
      CALL THREED
      STOP
      END
C
C
      SUBROUTINE THREED
      PARAMETER (MXPTS=51,MXPTSQ=MXPTS*MXPTS)
C
C       THE CONTOURING PORTIONS SEPARATED FROM THE ORBITAL COMPUTATION
C       SECTION AND PLACED IN THIS PROGRAM.   OTHER THAN NECCESSARY
C       CHANGES IN THE INPUT, THIS ROUTINE IS THAT WRITTEN BY
C       WILLIAM JORGENSEN.
C                               DAN SEVERANCE 12-87
C
C       MODIFIED TO OUTPUT ADDITIONAL INFORMATION TO BE USED IN HIDDEN
C       LINE PORTION OF THE PROGRAM.         DAN SEVERANCE  2-88
C
      COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,XX,YY,ZZ,IDASH,SCALE,PERZ,
     *   KA,KO,ZZZ
      COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850)
      COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3
      COMMON // DENS(MXPTS,MXPTS,MXPTS),SAV(MXPTSQ),CTR(15)
C
C     PUT DENS IN A COMMON BLOCK FOR THOSE MACHINES WHICH CAN MAP
C     IT TO A LARGER MEMORY AREA.  THIS MAY BE CONVERTED TO A DIMENSION
C     STATEMENT IF PREFERRED.
C
      REAL C(50,3),XX(850),YY(850),ZZ(850),DMIN,DMAX
      INTEGER IAN(49)
      CHARACTER CALC*20
      DATA AU,PI,RT3 / 0.5291670E+0,3.14159265360E+0,1.73205080E+0 /
C
C       DENS IS THE CHARGE DENSITY OR ORBITAL VALUE MATRIX. IAN IS THE
C       ARRAY OF ATOMIC NUMBERS FOR THE ATOMS.  ATOMIC COORDINATES ARE
C       READ IN ANGSTROMS AND STORED IN ARRAY C.
C       CTR IS THE ARRAY OF CONTOUR LEVELS IN AU.
C
C       CARD INPUT IS AS FOLLOWS - THIS IS NOT ALL USED IN THIS PROGRAM
C               BUT THE SAME INPUT DECK IS USED BY THE ORBITAL VALUE
C               PACKAGE
C
C       CARD  1 - TYPE OF WAVEFUNCTION = STO-3G , 3-21G, 3-21+G,
C                   6-31G TO 6-31++G**. (A20)
C                   (THIS IS FOR INFORMATIONAL USE ONLY, THE PROGRAM
C                    DOESN'T USE IT)
C
C       CARD 2 - NO. OF CONTOURS, ICONN, ICONZ, NORB  -4I2
C         ICONN = 1 FOR NEGATIVE CONTOURS, TOO.
C         ICONZ=1 FOR ZERO CONTOUR AND NORB=2 FOR CHARGE DENSITY PLOT.
C       CARD 3 - CONTOUR LEVELS (ONLY POSITIVE NEED BE SPECIFIED WHEN
C         ICONN AND ICONZ ARE USED - 8F10.6).
C
C       INITIALIZE FOR AUTOMATIC PROCESSING
C
      NCT = 1
      CTR(1) = 0.0750E+0
      NORB = 1
      ICONN = 1
      ICONZ = 0
      READ (IRD,10) CALC
   10 FORMAT (A)
C
C       READ THE NO. OF CONTOURS AND THEN THE LEVELS.
C
      READ (IRD,20) NCT,ICONN,ICONZ,NORB
   20 FORMAT (4I2)
      READ (IRD,30) (CTR(I),I=1,NCT)
   30 FORMAT (8F10.6)
C
C       DEFAULT VALUES
C
      IF (SCALE.EQ.0.0E+0) SCALE = 1.00E+0
      IF (NCT.EQ.0) CTR(1) = 0.0750E+0
      IF (NCT.EQ.0) NCT = 1
      IF (NORB.EQ.0) NORB = 1
      IF (NORB.EQ.1) ICONN = 1
C
C       READ THE ORBITAL VALUE MATRIX
C
      WRITE (ILST,*) 'READING 17 '
      OPEN (17,FILE='FOR017',FORM='UNFORMATTED',STATUS='OLD')
      READ (17) MAXPTS,NAT
      READ (17) (IAN(I),I=1,NAT)
      READ (17) ((C(I,J),J=1,3),I=1,NAT)
      READ (17) (((DENS(I,J,K),I=1,MAXPTS),J=1,MAXPTS),K=1,MAXPTS)
      READ (17) XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
C
C       DETERMINE DEFAULT RANGES, INCREASING SCALE
C       INCREASES THE RANGE OF THE PLOT.
C
      SPACES = FLOAT(MAXPTS-1)*AU
      XMI = XMIN/AU
      YMI = YMIN/AU
      ZMI = ZMIN/AU
      XINC = (XMAX-XMIN)/SPACES
      YINC = (YMAX-YMIN)/SPACES
      ZINC = (ZMAX-ZMIN)/SPACES
      DMIN = 1000.00E+0
      DMAX = -1000.00E+0
      DO 40 K = 1, MAXPTS
         DO 40 J = 1, MAXPTS
            DO 40 I = 1, MAXPTS
               DMIN = MIN(DMIN,DENS(I,J,K))
               DMAX = MAX(DMAX,DENS(I,J,K))
   40 CONTINUE
      WRITE (ILST,*) 'MIN, MAX DENSITY(VALUE) COMPUTED IS ',DMIN,DMAX
      OPEN (IDSK3,FILE='FOR024',STATUS='UNKNOWN')
      OPEN (IDSK2,FILE='FOR023',STATUS='UNKNOWN')
      OPEN (IDSK1,FILE='FOR022',STATUS='UNKNOWN')
      KA = 1
      KO = 2
      ZZZ = ZMIN
      ZINC = ZINC*AU
      DO 70 I = 1, MAXPTS
         KNT = 0
         DO 60 J = 1, MAXPTS
            DO 50 K = 1, MAXPTS
               KNT = KNT+1
               SAV(KNT) = DENS(K,J,I)
   50       CONTINUE
   60    CONTINUE
         CALL TRDPLT (SAV,MAXPTS,MAXPTS,XMIN,XMAX,YMIN,YMAX,NCT,CTR,
     *      ICONN,ICONZ)
         ZZZ = ZZZ+ZINC
   70 CONTINUE
      NC1 = NC
      NCURV1 = NCURV
      KA = 1
      KO = 3
      ZZZ = YMIN
      YINC = YINC*AU
      DO 100 I = 1, MAXPTS
         KNT = 0
         DO 90 J = 1, MAXPTS
            DO 80 K = 1, MAXPTS
               KNT = KNT+1
               SAV(KNT) = DENS(K,I,J)
   80       CONTINUE
   90    CONTINUE
         CALL TRDPLT (SAV,MAXPTS,MAXPTS,XMIN,XMAX,ZMIN,ZMAX,NCT,CTR,
     *      ICONN,ICONZ)
         ZZZ = ZZZ+YINC
  100 CONTINUE
      NC2 = NC
      NCURV2 = NCURV
      KA = 2
      KO = 3
      ZZZ = XMIN
      XINC = XINC*AU
      DO 130 I = 1, MAXPTS
         KNT = 0
         DO 120 J = 1, MAXPTS
            DO 110 K = 1, MAXPTS
               KNT = KNT+1
               SAV(KNT) = DENS(I,K,J)
  110       CONTINUE
  120    CONTINUE
         CALL TRDPLT (SAV,MAXPTS,MAXPTS,YMIN,YMAX,ZMIN,ZMAX,NCT,CTR,
     *      ICONN,ICONZ)
         ZZZ = ZZZ+XINC
  130 CONTINUE
      NC3 = NC
      NCURV3 = NCURV
      WRITE (IDSK2,140) NC,NCURV,NC1,NCURV1,NC2,NCURV2,NC3,NCURV3
      NCMAX = 0
  140 FORMAT (10I6)
  150 FORMAT (20I4)
      WRITE (IDSK2,150) (NAC(I),NCDASH(I),I=1,NCURV)
      NCMAX = 0
      DO 160 I = 1, NCURV
         NCMAX = MAX(NCMAX,NAC(I))
  160 CONTINUE
      WRITE (ILST,170) NC,NCURV,NCMAX
  170 FORMAT (6X,'NC = ',I5,' NCURV = ',I4,' NCMAX =',I4/)
      RETURN
      END
C
C
      SUBROUTINE TRDPLT (ARRAY,IXDIM,IYDIM,XXMIN,XXMAX,YYMIN,YYMAX,NCTR,
     *   CTR,ICONN,ICONZ)
      COMMON /SCLDAT/ XMAX,XMIN,YMAX,YMIN,XINCR,YINCR,CL,NPLT
      COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(850),Y(850),Z(850),IDASH,
     *   SCALE,PERZ,KA,KO,ZZ
      COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850)
      COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3
      DIMENSION ARRAY(IXDIM,IYDIM),CTR(15)
      DATA NCSAVE / 0 /
C
C       THIS PROGRAM PRODUCES CONTOUR PLOTS FROM THE MATRIX,
C       ARRAY, WHICH HAS DIMENSIONS IXDIM AND IYDIM USING THE
C       CONTOURS,CTR. NCTR = THE NUMBER OF CONTOURS. KO = 1,2,3
C       IF THE ORDINATE IS THE X,Y OR Z AXIS% KA IS SIMILAR
C       FOR THE ABSCISSA. ICONN = 1 IF NEGATIVE CONTOURS ARE
C       DESIRED THAT ARE THE SAME AS THE POSITIVE CONTOURS (ONLY
C       THE POSITIVE ONES NEED BE SPECIFIED IN CTR AND NCTR).
C       ICONZ = 1 IF A DOTTED ZERO CONTOUR IS DESIRED.
C       W. L. JORGENSEN   -   DECEMBER 6, 1971
C
      XMAX = XXMAX
      XMIN = XXMIN
      YMAX = YYMAX
      YMIN = YYMIN
      XINCR = (XMAX-XMIN)/FLOAT(IXDIM-1)
      YINCR = (YMAX-YMIN)/FLOAT(IYDIM-1)
C
C       CALL CONTUR FOR EACH CONTOUR.
C
      NPLT = 0
      IDASH = 0
      DO 10 NCL = 1, NCTR
         CL = CTR(NCL)
         NPLT = NPLT+1
         CALL CONTUR (ARRAY,IXDIM,IYDIM,CL,IERROR)
         IF (IERROR.EQ.1) THEN
            WRITE (ILST,40)
            STOP
         ENDIF
         IDASH = IDASH + 2*NCL
   10 CONTINUE
C
C       FOR ICONN = 1, DO AGAIN WITH - CONTOURS.
C
      IF (ICONN.EQ.1) THEN
         NPLT = 0
         IDASH = 1
         DO 20 NCL = 1, NCTR
            CL = -CTR(NCL)
            NPLT = NPLT+1
            CALL CONTUR (ARRAY,IXDIM,IYDIM,CL,IERROR)
            IF (IERROR.EQ.1) THEN
               WRITE (ILST,40)
               STOP
            ENDIF
            IDASH = IDASH + 2*NCL
   20    CONTINUE
      ENDIF
C
C       FOR ICONZ=1, PLOT 0 CONTOUR.
C
      IF (ICONZ.EQ.1) THEN
         CL = 0.00E+0
         NPLT = 0
         CALL CONTUR (ARRAY,IXDIM,IYDIM,0.00E+0,IERROR)
         IF (IERROR.EQ.1) THEN
            WRITE (ILST,40)
            STOP
         ENDIF
      ENDIF
      IF (NCURV.NE.NCSAVE) THEN
         WRITE (IDSK3,30) NCURV
         NCSAVE = NCURV
      ENDIF
   30 FORMAT (I6)
      RETURN
   40 FORMAT ('0** CONTOUR PASSING OUTSIDE THE PLOTTING BOX **'/
     *   '0** INCREASE THE CONTOUR LEVEL AND RERUN PSICON OR **'/
     *   '0** INCREASE THE SCALE FACTOR FOR PSI1 AND RERUN IT **')
      END
C
C
      SUBROUTINE ROTPLT
      COMMON /RPLOT/ N,CO(3),CM,THE,GAM,PHI,X(850),Y(850),Z(850),IDASH,
     *   SCALE,PERZ,KA,KO,ZO
      COMMON /HLCOM/ NAC(850),NC,NCURV,NCDASH(850)
      COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3
      DIMENSION SAV(3)
C
C       THIS ROUTINE STORES THE CURVES WHICH HAVE
C       BEEN GENERATED BY THE CONTUR PLOTTING
C       ROUTINES IN THE DISK FILE NAMED CURVE
C
C       TRANSFORM CURVE TO ORIGINAL SYSTEM
C
      KZ = 6-KA-KO
      DO 10 I = 1, N
         SAV(KA) = X(I)
         SAV(KO) = Y(I)
         SAV(KZ) = ZO
         X(I) = SAV(1)
         Y(I) = SAV(2)
         Z(I) = SAV(3)
   10 CONTINUE
C
C       STORE CURVE ON DISK
C
      NCURV = NCURV+1
      NCDASH(NCURV) = IDASH
      NC = NC+N
      NAC(NCURV) = N
      IF (N.GT.850) THEN
         WRITE (ILST,*) '>>>>> WARNING:XYZ REC NEED TO BE REDIMENSIONED'
         WRITE (ILST,*) '>>>>>>THEY NEED ',N
      ENDIF
      WRITE (IDSK1,20) (X(I),Y(I),Z(I),I=1,N)
   20 FORMAT (8F10.6)
      RETURN
      END
C
C
      SUBROUTINE CONTUR (A,NX,NY,CLP,JERROR)
      COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,X,Y,Z,IDASH,SCALE,PERZ,KA,
     *   KO,ZO
      COMMON /CRPTPM/ NRMAX,NPMAX,LX,UX,LY,UY,CL,NR,NP,IXO,IYO,ISO,PC,
     *   LOPC
      COMMON /TEMP/ REC
      DIMENSION INX(8),INY(8),A(NX,NY),X(850),Y(850),REC(850),Z(850)
      INTEGER UX,UY,UX1,UY1,COM,REC
      LOGICAL PC,LOPC
      DATA INX / -1,-1,0,3*1,0,-1 /
      DATA INY / 0,3*1,0,3*-1 /
      NF = 1
      LX = NF
      LY = NF
      UX = NX
      UY = NY
      CL = CLP
      NRMAX = 850
      NPMAX = 850
C
C       THE 3 ROUTINES, CONTUR,CURVE AND INTERP FORM A
C       GENERAL CONTOUR MAPPING PACKAGE THAT SIMPLY REQUIRES
C       A PLOTTING ROUTINE PLOT(N,X,Y) TO PLOT N POINTS
C       WITH COORDINATES IN THE ARRAYS X AND Y.
C       IF CONTOURS WITH MORE THAN 400 POINTS ARE TO BE
C       DRAWN (THIS IS VERY LARGE) , NRMAX AND NPMAX ABOVE
C       MUST BE ADJUSTED ALONG WITH THE DIMENSION STATEMENTS
C       THIS PROGRAM IS BASED ON M.O. DAYHOFF S PAPER
C        A CONTOUR MAP PROGRAM FOR X-RAY CRYSTALLOGRAPHY
C       OCT. 1963 COMMUNICATIONS OF THE ACM. FIRST
C       PROGRAMMED BY BRUCE LANGDON, PLASMA PHYSICS
C       LAB,PRINCETON UNIVERSITY, NOV. 1966.
C
      NR = 0
      JERROR = 0
C
C       SCAN RIGHT EDGE AND BOTTOM
C
      LY1 = LY+1
      DO 10 J = LY1, UY
         IF (A(UX,J-1).LT.CL.AND.A(UX,J).GE.CL) THEN
            CALL CURVE (A,NX,UX,J,7,IER)
            JERROR = 1
C
C           WRITE (ILST,*) ' RIGHT EDGE '
C
            IF (IER.EQ.1) GO TO 70
         ENDIF
   10 CONTINUE
C
C       SCAN TOP EDGE, LEFT TO RIGHT
C
      UX1 = UX-1
      DO 20 I1 = LX, UX1
         I = UX1+LX-I1
         IF (A(I+1,UY).LT.CL.AND.A(I,UY).GE.CL) THEN
            CALL CURVE (A,NX,I,UY,5,IER)
            JERROR = 1
C
C           WRITE (ILST,*) ' TOP EDGE '
C
            IF (IER.EQ.1) GO TO 70
         ENDIF
   20 CONTINUE
C
C       SCAN LEFT EDGE TOP TO BOTTOM
C
      UY1 = UY-1
      DO 30 J1 = LY, UY1
         J = UY1+LY-J1
         IF (A(LX,J+1).LT.CL.AND.A(LX,J).GE.CL) THEN
            CALL CURVE (A,NX,LX,J,3,IER)
C
C           WRITE (ILST,*) ' LEFT EDGE EDGE '
C
            JERROR = 1
            IF (IER.EQ.1) GO TO 70
         ENDIF
   30 CONTINUE
C
C       SCAN BOTTOM EDGE AND INTERIOR POINTS
C
      LX1 = LX+1
      DO 60 J = LY, UY1
         DO 50 I = LX1, UX
            IF (A(I-1,J).LT.CL.AND.A(I,J).GE.CL) THEN
               IF (NR.NE.0) THEN
                  COM = (UY-LY+1)*I+J
                  DO 40 ID = 1, NR
                     IF (REC(ID).EQ.COM) GO TO 50
   40             CONTINUE
               ENDIF
               CALL CURVE (A,NX,I,J,1,IER)
               IF (J.EQ.LY) JERROR = 1
               IF (IER.EQ.1) GO TO 70
            ENDIF
   50    CONTINUE
   60 CONTINUE
   70 RETURN
      END
C
C
      SUBROUTINE CURVE (A,IDIM,IXP,IYP,ISP,IER)
      INTEGER UX,UY,REC,DX,DY
      LOGICAL PC,LOPC,LPC
      DIMENSION A(IDIM,IDIM),X(850),Y(850),REC(850),INX(8),INY(8)
      DIMENSION Z(850)
      COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,X,Y,Z,IDASH,SCALE,PERZ,KA,
     *   KO,ZO
      COMMON /CRPTPM/ NRMAX,NPMAX,LX,UX,LY,UY,CL,NR,NP,IXO,IYO,ISO,PC,
     *   LOPC
      COMMON /TEMP/ REC
      COMMON /IO/ IRD,ILST,IDSK1,IDSK2,IDSK3
      DATA INX / -1,-1,0,3*1,0,-1 /
      DATA INY / 0,3*1,0,3*-1 /
C
C
      IER = 0
      IS = ISP
      IX = IXP
      IY = IYP
      IXO = IX
      IYO = IY
      ISO = IS
      NP = 1
      LOPC = .FALSE.
      CALL INTERP (IX,IY,IS,A,IDIM,JER)
      GO TO (10,100,100,40), JER
C
C       ROTATE
C
   10 IS = IS+1
      IF (IS.EQ.9) IS = 1
      DX = INX(IS)
      DY = INY(IS)
C
C       FIND PLOT POINT
C
   20 CALL INTERP (IX,IY,IS,A,IDIM,JER)
      GO TO (30,60,80,110), JER
   30 IF (IS.NE.1) GO TO 10
C
C       RECORD A  CENTER
C
   40 IF (NR.GE.NRMAX) THEN
C
C       REC ARRAY OVERFLOW
C
         WRITE (ILST,50)
   50    FORMAT ('0TOO MANY POINTS IN CONTOUR. IT IS TERMINATED')
         NNN = NP-1
         CALL ROTPLT
         IER = 1
         RETURN
      ENDIF
      NR = NR+1
      REC(NR) = UY*IX+IY
      GO TO 10
C
C       DIAG FAIL
C
   60 IS = IS+1
      IF (IS.EQ.9) IS = 1
      LOPC = .TRUE.
      CALL INTERP (IX+INX(IS),IY+INY(IS),IS-3,A,IDIM,JER)
      GO TO (70,100,100,110), JER
   70 IX = IX+DX
      IY = IY+DY
      IS = IS+4
      IF (IS.GT.8) IS = IS-8
      GO TO 20
C
C       NON-DIAG FAIL
C
   80 IX = IX+DX
      IY = IY+DY
      IS = IS+5
      IF (IS.GT.8) IS = IS-8
      LPC = PC
      CALL INTERP (IX,IY,IS,A,IDIM,JER)
      GO TO (90,90,100,110), JER
   90 IF (.NOT.(LPC.AND.PC)) GO TO 10
C
C       PLOT POINT SWITCH
C
      TEM = X(NP-2)
      X(NP-2) = X(NP-1)
      X(NP-1) = TEM
      TEM = Y(NP-2)
      Y(NP-2) = Y(NP-1)
      Y(NP-1) = TEM
      GO TO 10
  100 STOP
  110 NNN = NP-1
      CALL ROTPLT
      RETURN
      END
C
C
      SUBROUTINE INTERP (IX,IY,ISP,A,IDIM,JER)
      INTEGER UX,UY,REC,DX,DY
      LOGICAL PC,LOPC
      DIMENSION A(IDIM,IDIM),X(850),Y(850),REC(850),INX(8),INY(8)
      DIMENSION Z(850)
      COMMON /RPLOT/ NNN,CO(3),CM,THE,GAM,PHI,X,Y,Z,IDASH,SCALE,PERZ,KA,
     *   KO,ZO
      COMMON /CRPTPM/ NRMAX,NPMAX,LX,UX,LY,UY,CL,NR,NP,IXO,IYO,ISO,PC,
     *   LOPC
      COMMON /SCLDAT/ XMAX,XMIN,YMAX,YMIN,XINCR,YINCR,CLX,NPLT
      COMMON /TEMP/ REC
      DATA INX / -1,-1,0,3*1,0,-1 /
      DATA INY / 0,3*1,0,3*-1 /
C
      JER = 1
      IS = ISP
      IF (IS.LT.1) IS = IS+8
      DX = INX(IS)
      DY = INY(IS)
      FDX = FLOAT(DX)*XINCR
      FDY = FLOAT(DY)*YINCR
      FIX = FLOAT(IX-1)*XINCR+XMIN
      FIY = FLOAT(IY-1)*YINCR+YMIN
      IX1 = IX+DX
      IY1 = IY+DY
      IF (IX1.GE.LX.AND.IX1.LE.UX.AND.IY1.GE.LY.AND.IY1.LE.UY) THEN
         IF (DX*DY.EQ.0) THEN
C
C       NON-DIAGONAL CASE
C       CHECK FOR FAIL
C
            IF (A(IX1,IY1).GE.CL) GO TO 10
            IF (DX.EQ.0) THEN
               X(NP) = FIX
               Y(NP) = (A(IX,IY)-CL)/(A(IX,IY)-A(IX,IY1))*FDY+FIY
            ELSE
               Y(NP) = FIY
               X(NP) = (A(IX,IY)-CL)/(A(IX,IY)-A(IX1,IY))*FDX+FIX
            ENDIF
            PC = .FALSE.
         ELSE
C
C       DIAGONAL CASE
C
            CP = (A(IX,IY)+A(IX1,IY)+A(IX,IY1)+A(IX1,IY1))*0.250E+0
            IF ((CP.GE.CL.OR.LOPC)) THEN
C
C       CONTOUR PASSES ON FAR SIDE OF CENTER POINT
C
               IF (A(IX1,IY1).GE.CL.AND.CP.GE.CL) GO TO 20
               V = (A(IX1,IY1)-CL)/(A(IX1,IY1)-CP)*0.50E+0
               X(NP) = (1.0E+0-V)*FDX+FIX
               Y(NP) = (1.0E+0-V)*FDY+FIY
               PC = .TRUE.
               LOPC = .FALSE.
            ELSE
C
C       CONTOUR PASSES ON NEAR SIDE OF CENTER POINT
C
               V = (A(IX,IY)-CL)/(A(IX,IY)-CP)*0.50E+0
               X(NP) = V*FDX+FIX
               Y(NP) = V*FDY+FIY
               PC = .FALSE.
            ENDIF
         ENDIF
         NP = NP+1
         IF (IX.NE.IXO.OR.IY.NE.IYO.OR.IS.NE.ISO) THEN
            IF (NP.LE.NPMAX) RETURN
C
C       PLOT PART OF CURVE AND CONTINUE
C
            NNN = NP-1
            CALL ROTPLT
            X(1) = X(NP-1)
            Y(1) = Y(NP-1)
            NP = 2
            RETURN
         ENDIF
      ENDIF
      JER = JER+1
   10 JER = JER+1
   20 JER = JER+1
      RETURN
      END