	PROGRAM CNTPLT
C CONTOUR PLOTTING PROGRAM
C
	PARAMETER MAXSIZ=6000
      DIMENSION  SBEG(40),CSTEP(40),SLIM(40),XLINE(40),YLINE(40),BX(40),
     2 BY(40),I1(40),I2(40),I3(40),NSIG(40),ICRSLV(40),IPLEV(40)
      DIMENSION IC(1)
      DOUBLE PRECISION LABL
      COMMON /B1/ A,B,ALF,NA,NB,AA,BB,AX,X1,X2,Y1,Y2,CALF,SALF,
     1ATITLE(20),C(MAXSIZ) 
      COMMON /PLT1/ ITOT,IN,NP,IT,NT,UINT,LOG(6),NDIG
	BYTE TIMBUF(8),DATEBUF(9)
C
      EQUIVALENCE (C,IC)
      DATA X2,Y2/0.0,0.0/
      DATA BLNK/1H /
      DATA IN,IOUT,IFILE/5,6,11/
C
C GET DATE AND TIME
C
	CALL FOR$DATE(DATEBUF)
	CALL FOR$TIME(TIMBUF)
C
      INITAL=0
C
C     READ PROBLEM CONTROL INFORMATION
C
C
      READ(IN,157) ATITLE,TLEHGT,NLABL,A,B,ALP,NA,NB,NP,NCONT,NLINE,NCRO
     1SS,SCROSS,NSECT,NZ
      NRUSD=NA
      NCUSD=NB
      IF(NP.LE.0)NP=4
      REWIND IFILE
      WRITE(IOUT,158) ATITLE,A,B,ALP,NA,NB,NP
      IF (NCONT.GT.0) GO TO 102 
      WRITE(IOUT,160)
      GO TO 117 
  102 NCONT=MIN0(NCONT,40)
      DO 103 J=1,NCONT
      READ(IN,179)SBEG(J),CSTEP(J),SLIM(J),NSIG(J),IPLEV(J)
      IF(IPLEV(J).EQ.0)IPLEV(J)=1
      IF(IPLEV(J).GT.4)IPLEV(J)=1
  103 CONTINUE
      NLINE=MIN0(NLINE,40)
      IF (NLINE.LE.0) GO TO 105 
      READ(IN,162) (XLINE(I),YLINE(I),I=1,NLINE)
  105 ICTYP=0
      IF(NCROSS.LT.0)ICTYP=1
      NCROSS=IABS(NCROSS)
      NCROSS=MIN0(NCROSS,40)
      IF (NCROSS.LE.0) GO TO 106
      READ(IN,162) (BX(I),BY(I),I=1,NCROSS)
  106 NSECT=MIN0(NSECT,40)
      IF (NSECT.LE.0) GO TO 107 
      READ(IN,163) (I1(I),I2(I),I3(I),I=1,NSECT)
  107 WRITE(IOUT,164) NCONT,(SBEG(I),CSTEP(I),SLIM(I),NSIG(I),I=1,NCONT)
      IF (NLINE.LE.0) GO TO 108 
      WRITE(IOUT,165) NLINE,(XLINE(I),YLINE(I),I=1,NLINE)
  108 IF (NCROSS.LE.0) GO TO 109
      WRITE(IOUT,166) NCROSS,SCROSS,(BX(I),BY(I),I=1,NCROSS)
      SCROSS=SCROSS*.5
      IF(ICTYP.EQ.0)GO TO 109
C READ CROSS LEVELS
      READ(IN,180)(ICRSLV(I),I=1,NCROSS)
  109 IF (NSECT.GT.0) GO TO 110 
      WRITE(IOUT,167)
      GO TO 117 
  110 WRITE(IOUT,168) NSECT,NZ,(I1(I),I2(I),I3(I),I=1,NSECT)
      ISECT=1
      GO TO 120 
  111 INITAL=1
C * * * * * SYSTEM DEPENDENT STATEMENT FOLLOWS ***************
      CALL PLOTS(0,0,12)
      ALF=ALP*.017453
      CALF=COS(ALF)
      SALF=SIN(ALF)
      YMAX=10.5 
      IF(A.LT.0)YMAX=30.
      A=ABS(A)
      UINT=1./FLOAT(NP+1)
      AA=A/FLOAT(NA-1)
      BB=B/FLOAT(NB-1)
      AX=0.0
      IF(CALF)112,113,113
  112 AX=A*ABS(CALF)
  113 X1=3.2*TLEHGT
      XFIN=X1+B+A*ABS(CALF)
      X1=X1+AX
      Y1=0.0
      CALL PLOT(0.,.5,-3)
      IF(A*SALF-YMAX)124,124,119
  116 ISECT=ISECT+1
      IF (ISECT.GT.NSECT) GO TO 117
      REWIND IFILE
      GO TO 120 
  117 GO TO 146 
  118 WRITE(IOUT,169) JZ,NZ
      GO TO 117 
  119 WRITE(IOUT,170) Y2,YMAX
      GO TO 117 
  120 JZ=I1(ISECT)
      J=(I3(ISECT)-JZ)*I2(ISECT)
      IF (J.LT.0) GO TO 116
      I1(ISECT)=I1(ISECT)+I2(ISECT)
      WRITE(IOUT,171) JZ,NZ
      NROW=NRUSD
      NCOL=NCUSD
      NTOT=NCOL*NROW
	IF(NTOT.GT.MAXSIZ)STOP 'NOT ENOUGH MEMORY'
  121 READ(IFILE,END=118)IPTSLN,ILNPLN,KIZ,KNZ
      IF (IPTSLN.EQ.NRUSD.AND.ILNPLN.EQ.NCUSD) GO TO 123
      WRITE(IOUT,173) IPTSLN,ILNPLN,NRUSD,NCUSD 
      GO TO 117 
  123 READ(IFILE,END=118)(C(I),I=1,NTOT)
      IF (KIZ.NE.JZ) GO TO 121
      IF(INITAL.EQ.0)GO TO 111
  124 WRITE(IOUT,174) JZ,NZ
      CALL INTPLT
      BMIN=10000.
      BMAX=-BMIN
      DO 127 I=1,NTOT
	 BMIN=AMIN1(BMIN,C(I))
	 BMAX=AMAX1(BMAX,C(I))
  127 CONTINUE
      CALL NEWPEN(1)
      IF (NLINE.LE.0) GO TO 134 
      DO 131 I=1,NLINE
      WY=ABS(XLINE(I))
      WX=ABS(YLINE(I))
      WX=WX*BB+WY*AA*CALF+X1
      WY=WY*AA*SALF+Y1
      IPEN=3
      IF(I.EQ.1.OR.XLINE(I).LT.0)GO TO 131
      IF(YLINE(I).LT.0)GO TO 134
      IPEN=2
  131 CALL PLOT(WX,WY,IPEN)
  134 IF (NCROSS.LE.0) GO TO 137
      DO 136 I=1,NCROSS 
      IF(ICTYP.EQ.0)GO TO 135
C DO NOT PLOT CROSS IF NOT ON THIS LEVEL
      IF(ICRSLV(I).NE.JZ)GO TO 136
135   IF(BX(I).LT.0.0)GO TO 137 
      XX=BY(I)*BB+BX(I)*AA*CALF+X1
      YY=BX(I)*AA*SALF+Y1
      CALL PLOT(XX-SCROSS,YY,3) 
      CALL PLOT(XX+SCROSS,YY,2) 
      CALL PLOT(XX,YY-SCROSS,3) 
      CALL PLOT(XX,YY+SCROSS,2) 
  136 CONTINUE
  137 IF(TLEHGT.LE.0.0)GO TO 1369
      CALL SYMBOL(TLEHGT,0.0,TLEHGT,ATITLE,90.,40)
      TEMP=0.5*TLEHGT
      CALL SYMBOL(2.0*TLEHGT,.2,TEMP,7HSection,90.,7)
      CALL NUMBER(2.0*TLEHGT,10.*TEMP,TEMP,FLOAT(JZ),90.,-1)
	CALL SYMBOL(2.0*TLEHGT,999.,TEMP,1H/,90.,1)
	CALL NUMBER(2.0*TLEHGT,999.,TEMP,FLOAT(NZ),90.,-1)
	CALL SYMBOL(2.0*TLEHGT,999.,TEMP,4H    ,90.,4)
	CALL SYMBOL(2.0*TLEHGT,999.,TEMP,DATEBUF,90.,9)
	CALL SYMBOL(2.0*TLEHGT,999.,TEMP,4H    ,90.,4)
	CALL SYMBOL(2.0*TLEHGT,999.,TEMP,TIMBUF,90.,8)
 1369 IF(NLABL.EQ.0)GO TO 1371
C	LABEL THE DESIRED VERTICES
      DO 1370 I=1,NLABL 
      READ(IN,178)WY,WX,DX,DY,TEMP,LABL 
      XX=WX*BB+WY*AA*CALF+X1+DX 
      YY=WY*AA*SALF+Y1+DY
 1370 CALL SYMBOL(XX,YY,TEMP,LABL,90.,6)
 1371 I=1
      ICC=0
  138 CBEG=AMIN1(SBEG(I),SLIM(I))
      CLIM=AMAX1(SBEG(I),SLIM(I))
      ACSTEP=AMAX1(1.,ABS(CSTEP(I)))
      CALL NEWPEN(IPLEV(I))
      IF (CBEG.GT.BMAX.OR.CLIM.LT.BMIN) GO TO 142
  139 IF (CBEG.GE.BMIN) GO TO 140
      Q=CBEG+ACSTEP
      IF (Q.GT.CLIM) GO TO 140
      CBEG=Q
      GO TO 139 
  140 IF (CLIM.LE.BMAX) GO TO 141
      Q=CLIM-ACSTEP
      IF (Q.LE.CBEG) GO TO 141
      CLIM=Q
      GO TO 140 
  141 ICC=1
      NDIG=NSIG(I)
      CALL ICONT (C,NROW,NCOL,NRUSD,NCUSD,CBEG,CLIM,ACSTEP)
      call newpen(1)
      GO TO 143 
  142 WRITE(IOUT,176) CBEG,CLIM,BMIN,BMAX
  143 I=I+1
      IF (I.LE.NCONT) GO TO 138 
      IF (ICC.GT.0) GO TO 144
      WRITE(IOUT,177)
  144 INITAL=1
      CALL PLOT(XFIN+3.0*TLEHGT,0.,-3)
      IF(INITAL)146,120,120
  146 CALL PLOT(1.0,-0.5,999)
      STOP
C
  157 FORMAT(20A4/F5.0,I5,3F5.0,6I5,F5.0,2I5)
  158 FORMAT(1X,20A4/1X,10H EDGES OF ,2F7.2, 20H INS. WITH ANGLE OF ,F
     17.2,  9H DEGREES.,/1X2I4, 21H POINTS PER EDGE AND ,I2, 32H INTERPO
     2LATION POINTS PER SPLINE) 
  160 FORMAT (1H0, 22H NO CONTOURS SPECIFIED)
  162 FORMAT (16F5.0)
  163 FORMAT (3I5)
  164 FORMAT(1H0I5,14H CONTOUR TYPES,/(1X,3F8.2,I4))
  165 FORMAT (1H0I5,  9H VERTICES,/(1X5(4X2F8.2)))
  166 FORMAT (1H0I5, 24H CROSSES WITH ARM LENGTH,F6.2,5H INS.,/(1X5(4X
     12F8.2)))
  167 FORMAT (1H0, 42H NO SECTION SELECTIONS. - - JOB ABANDONED.)
  168 FORMAT (1H0I5, 44H SECTION SELECTIONS GIVEN IN INTERVALS OF 1/,I4/
     1(1X,  5HFROM ,I4, 13H IN STEPS OF ,I4,  4H TO ,I4))
  169 FORMAT (1H0, 33H E O F  READ WHILE SEARCHING FOR ,2I5)
170	FORMAT(' WIDTH TOO GREAT', 2F10.4)
  171 FORMAT (1H0, 17H PLOT OF SECTION ,I5,1H/,I5,  9H SELECTED) 
173	FORMAT(' SECTION SIZES DIFFER - JOB ABANDONED'
	1,/1X2I4, 12H SIZE READ, ,2I4, 14H SIZE SELECTED)
  174 FORMAT (1H , 17H PLOT OF SECTION ,I5,1H/I5,  6H BEGUN)
  176 FORMAT (1H0, 25H NO MAP DENSITY IN RANGE ,2F10.2, 20H AVAILABLE RA
     1NGE IS ,2F10.2)
  177 FORMAT (1H0, 50H NO CONTOURS FOR THIS SECTION. - - NOTHING PLOTTED
     1) 
  178 FORMAT(5F5.0,A6)
  179 FORMAT(3F5.0,2I5)
  180 FORMAT(20I5)
C
      END
      SUBROUTINE INTPLT 
C
C     USING RATIONAL POLN. SPLINE INT. OF D V AHUJA  I B M SYST J  1968 
C     P 208-217 
C
      COMMON /PLT2/ C(200),XA(200),YA(200),VX(200),VY(200)
      COMMON /PLT1/ ITOT,IN,NP,IT,NT,UINT,LOG(6)
      LOG(6)=200
      NT=200
      ITOT=200
      C(1)=.25
      CA=.25
      LOG(5)=0
      LOG(4)=0
      LOG(3)=0
      LOG(2)=0
      LOG(1)=0
      DO 101 I=2,ITOT
      C(I)=1./(4.-CA)
  101 CA=C(I)
      RETURN
	END
      SUBROUTINE PREPLT(ELEV)
      DIMENSION X(200), Y(200)
      COMMON /PLT2/ C(200),XA(200),YA(200),VX(200),VY(200)
      COMMON /PLT1/ ITOT,IN,NP,IT,NT,UINT,LOG(6)
      EQUIVALENCE (X,XA), (Y,YA)
      IPEN=2
C      IF(ELEV.LT.0.0)IPEN=5
C      IF(ELEV.EQ.0.0)IPEN=4
      NC=IN-2
      N=NC
      NVV=IN-1
      NV=IN
      DX=VX(IN) 
      DY=VY(IN) 
      DO 102 I=2,N
	DX=C(NC)*(3.*(X(NV)-X(NV-2))-DX)
	 VX(NVV)=DX
	 DY=C(NC)*(3.*(Y(NV)-Y(NV-2))-DY)
	 VY(NVV)=DY
	 NC=NC-1
	 NVV=NVV-1
	 NV=NV-1
  102 CONTINUE
      VX(2)=C(1)*(3.*(X(3)-X(1))-DX-VX(1))
      VY(2)=C(1)*(3.*(Y(3)-Y(1))-DY-VY(1))
      XIC=X(1)
      YIC=Y(1)
      CALL PLOT(XIC,YIC,3)
      IC=1
      AA9=VX(IC)
      AA10=VY(IC)
  103 XC=XIC
      XIC=X(IC+1)
      YC=YIC
      YIC=Y(IC+1)
      A9=AA9
      AA9=VX(IC+1)
      A10=AA10
      AA10=VY(IC+1)
      A1=2.*(XC-XIC)+A9+AA9
      A2=2.*(YC-YIC)+A10+AA10
      A5=3.*(XIC-XC)-A9-A9-AA9
      A6=3.*(YIC-YC)-A10-A10-AA10
      I=1
      U=0.
  104 U=U+UINT
      WX=((A1*U+A5)*U+A9)*U+XC
      WY=((A2*U+A6)*U+A10)*U+YC 
      CALL PLOT(WX,WY,IPEN)
      I=I+1
      IF (I.LE.NP) GO TO 104
      IC=IC+1
      IF (IC.LT.IN) GO TO 103
      RETURN
C
      END
      SUBROUTINE CPLOT (IPLOT,ROW,COL,ELEV)
      COMMON /B1/ A,B,ALF,NA,NB,AA,BB,AX,X1,X2,Y1,Y2,CALF,SALF,
     1 ATITLE(20)
      COMMON /PLT1/ ITOT,IN,NP,IT,NT,UINT,LOG1,LOG2,LOG3,LOG4,ICOUNT,NCO
     1UNT,NDIG
      COMMON /PLT2/ C(200),XA(200),YA(200),VX(200),VY(200)
      XXX=COL*BB+ROW*AA*CALF+X1 
      YYY=ROW*AA*SALF+Y1
      IF (IPLOT-1) 101,103,106
  101 IF (ICOUNT.LE.1) GO TO 102
      IRET=1
      GO TO 107 
  102 IF (LOG4.EQ.1) GO TO 103
      ICOUNT=0
  103 IF (ICOUNT.LT.NCOUNT) GO TO 104
      IRET=2
      GO TO 107 
  104 ICOUNT=ICOUNT+1
      XA(ICOUNT)=XXX
      YA(ICOUNT)=YYY
C
C     NOTE INTERCHANGE OF X AND Y
C
      RETURN
 1050 IF(NDIG.GE.0)CALL NUMBER(XXX,YYY,.1,ELEV,90.0,NDIG)
      RETURN
  106 IRET=3
      IF (LOG4.EQ.1) GO TO 118
  107 IF (LOG4.EQ.1) GO TO 117
      X=ABS(XA(1)-XA(ICOUNT))
      Y=ABS(YA(1)-YA(ICOUNT))
      IF (AMAX1(X,Y).LE.1.E-6) GO TO 112
      LOG3=0
      X=XA(1)
      Y=YA(1)
  108 BY=Y-Y1
      IF (AMIN1(ABS(BY),ABS(Y-Y2)).LE.1.E-4) GO TO 110
      BY=BY*CALF/SALF
      IF (AMIN1(ABS(X-X1-BY),ABS(X-X2-BY)).LT.1.E-4) GO TO 110
      IF (LOG3.GT.0) GO TO 116
      LOG3=1
  109 X=XA(ICOUNT)
      Y=YA(ICOUNT)
      GO TO 108 
  110 IF (LOG3.EQ.2) GO TO 116
      IF (LOG3.EQ.1) GO TO 111
      LOG3=2
      GO TO 109 
  111 LOG4=1
      ISAVE=ICOUNT
      GO TO (102,104,1050,123),IRET
  112 VX(ICOUNT)=(XA(2)-XA(ICOUNT-1))*0.5*UINT
      VX(1)=VX(ICOUNT)
      VY(ICOUNT)=(YA(2)-YA(ICOUNT-1))*0.5*UINT
      VY(1)=VY(ICOUNT)
  113 UA=AA*UINT*0.7
      UB=BB*UINT*0.7
      VA=XA(1)
      VB=YA(1)
      J=1
      DO 115 I=2,ICOUNT 
	 IF (ABS(XA(I)-VA).GT.UA) GO TO 114
	 IF (ABS(YA(I)-VB).GT.UB) GO TO 114
	 GO TO 115
  114	 J=J+1
	 VA=XA(I)
	 XA(J)=VA
	 VB=YA(I)
	 YA(J)=VB
  115 CONTINUE
      XA(J)=XA(ICOUNT)
      YA(J)=YA(ICOUNT)
      VX(J)=VX(ICOUNT)
      VY(J)=VY(ICOUNT)
      IN=J
      IF (IN.GT.2) CALL PREPLT(ELEV)
      XA(1)=XA(ICOUNT)
      YA(1)=YA(ICOUNT)
      ICOUNT=1
      GO TO (102,104,1050,123),IRET
  116 VX(1)=(XA(2)-XA(1))*UINT
      VY(1)=(YA(2)-YA(1))*UINT
      VX(ICOUNT)=(XA(ICOUNT)-XA(ICOUNT-1))*UINT 
      VY(ICOUNT)=(YA(ICOUNT)-YA(ICOUNT-1))*UINT 
      GO TO 113 
  117 IF (ISAVE.NE.ICOUNT) GO TO 118
      GO TO (102,104,1050,123),IRET
  118 LOG4=0
      X=ABS(XA(1)-XA(ISAVE+1))
      Y=ABS(YA(1)-YA(ISAVE+1))
      IF (AMAX1(X,Y).GT.1.E-4) GO TO 122
      I=ISAVE/2 
      K=ISAVE
      DO 119 J=1,I
	 X=XA(J)
	 XA(J)=XA(K)
	 XA(K)=X
	 Y=YA(J)
	 YA(J)=YA(K)
	 YA(K)=Y
  119 K=K-1
      ICOUNT=ICOUNT-1
      K=ISAVE+1 
  120 DO 121 J=ISAVE,ICOUNT
	 XA(J)=XA(K)
	 YA(J)=YA(K)
  121 K=K+1
      ISAVE=0
      GO TO 116 
  122 KRET=IRET 
      KCOUNT=ICOUNT
      KSAVE=ISAVE
      ICOUNT=KSAVE
      IRET=4
      GO TO 116 
  123 IRET=KRET 
      ICOUNT=KCOUNT-KSAVE
      K=KSAVE+1 
      ISAVE=1
      GO TO 120 
C
      END
      SUBROUTINE ICONT (DATA,NROW,NCOL,NRUSD,NCUSD,CBEG,CLIM,ACSTP)
C
C     CONTOURING SUBPROGRAM BY S. H. ZISK, MIT LINCOLN LAB, MARCH 1968. 
C     REWRITTEN IN USASI BASIC FORTRAN BY NM BRENNER MIT LL, MARCH 1968 
C     DO 4000 XLEV=CBEG,CLIM,CSTP  (EFFECTIVELY, NOT LITERALLY) 
C     CSTP IS ACSTP WITH THE SIGN CHOSEN TO INCREASE IN THE CORRECT
C     DIRECTION (FROM CBEG TO CLIM).
C     FOR EACH XLEV, PLOTS A CONTOUR OF DATA(I,J), OUTPUT BEING RETURNED
C     VIA SUCCESSIVE CALLS TO SUBROUTINES.  CPLOT(0 OR 1,...) FOR PEN
C     UP OR PEN DOWN TRACES.  CPLOT(2,...) FOR LABELING AT THIS POSITION
C     IF DESIRED.  FULL CALL IS CALL CPLOT(IPLOT,ROW,COL,XLEV), WHERE
C     ROW IS FROM 1. TO FLOAT(MAXX), COL IS FROM 1. TO FLOAT(MAXY).
C     CONTR EXTRUDES ONE POINT AT A TIME TO CPLOT.  IF DESIRED FOR
C     EFFICIENCY, CPLOT MAY SAVE UP POINTS IN A SCRATCH ARRAY.
C     WHEN IPLOT IS 0, CPLOT SHOULD JUST REMEMBER THIS (ROW,COL) POINT. 
C     WHEN IPLOT=1, CPLOT SHOULD DRAW A LINE SEGMENT FROM THE REMEMBER- 
C     ED POINT TO THE NEW POINT, AND THEN MAKE THE NEW POINT THE REMEM- 
C     BERED POINT.
C     A CONTOUR LABEL IS CALLED FROM CPLOT (IPLOT=2) ONLY ONCE PER
C     COMPLETE CONTOUR SEGMENT, SO THAT IF SEVERAL DIFFERENT LOOPS ARE
C     AT HEIGHT XLEV, THEY WILL EACH GET A LABEL.
C
C     THIS PROGRAM IS WRITTEN IN USASI BASIC FORTRAN IN AN ATTEMPT TO BE
C     COMPATIBLE WITH AS MANY EXISTING FORTRAN COMPILERS AS POSSIBLE.
C     AS WRITTEN, THIS PROGRAM ASSUMES DATA IS A FLOATING POINT ARRAY.
C
      DIMENSION DATA(1) 
      REAL MARK 
	INTEGER TSTMRK
C
C     DIMENSION DATA(NROW,NCOL), OF WHICH ONLY (MAXX,MAXY) ARE USED.
C     NOTE--THE LOW ORDER BIT OF EACH ENTRY DATA(I,J) IS LIABLE TO CHANG
C     THIS WILL NEGLIGIBLY AFFECT A FLOATING POINT ARRAY.
C
C     CORRECTED AND MODIFIED  FEB. 1970.  G F BIO. SC.
C
C     IRT AND IUP INDICATE THE DIRECTION OF THE LAST CONTOUR STEP.
C     IRT=0,1 INDICATES THAT THE LAST STEP HAD X DECREASING,INCREASING A
C     IUP=0,1 INDICATES THAT THIS STEP HAD Y DECREASING, INCREASING .
C     DIRECTION OF THE NEXT CONTOUR STEP IS DETERMINED BY THE VALUES AT 
C     A,B,C,AND D WITH THE LAST STEP ENDING ON DA OR AT D (SPECIAL CASE)
C     AND A,B,C,D THE VERTICES OF A SQUARE LABELLED IN ORDER.
C
      MAXX=NRUSD
      MAXY=NCUSD
      CSTP=SIGN(ACSTP,CLIM-CBEG)
      CTR=CBEG
C	CLEAR LEAST SIGNIFICANT BIT
  101 XLEV=UNMARK(CTR)
      JLOOP=1
      JUP=1 
      IX=1
      IY=1
      NTOT=NROW*MAXY
      DO 102 I=1,NTOT
  102 DATA(I)=UNMARK(DATA(I))
C
C     THE LEAST SIGNIFICANT BIT OF DATA(I) IS USED TO MARK THE LEFTHAND 
C     EDGE OF THE CONTOUR WHEN IT CROSSES THAT ROW.
C     SRCHX SCANS THE ROWS FROM LEFT TO RIGHT, AND BOTTOM TO TOP,
C     LOOKING FOR THE FIRST POINT ON THE XLEV CONTOUR.
C
  103 IUP=JUP
      ILOOP=JLOOP
      IRT=1
      IDUMP=0
      IF (IX-MAXX) 104,131,212
  104 L=IX+NROW*(IY-1)
C
C     IX,IY SAVE THE LAST POINT THAT WAS SEARCHED, TO RESTART THE
C     SEARCH AFTER FINISHING THE CURRENT CONTOUR LINE.
C
  105 DIA1=DATA(L)
      DJA1=UNMARK(DIA1) 
      IF (DJA1-XLEV) 106,116,106
  106 IF (IX-1) 213,107,116
  107 IF (IY-1) 214,112,108
  108 LMX=L-NROW
      DJA2=DATA(LMX)
      DID=UNMARK(DJA2)
      IF (DID-XLEV) 111,109,111 
  109 IF (DID-DJA2) 111,110,111 
  110 YPT=IY-1
      GO TO 115 
  111 IF (IY-MAXY) 112,116,215
C
C     LOOKING FOR AN INTERSECTION IN THE LEFT-HAND COLUMN BEFORE
C     SCANNING FOR INTERSECTIONS IN THE CURRENT ROW.
C
  112 LPX=L+NROW
      DID=UNMARK(DATA(LPX))
      IF (DID-XLEV) 113,116,113 
  113 IF ((DID-XLEV)*(DJA1-XLEV)) 114,114,116
C
C     ARE BOTH DID AND DJA1 ON THE SAME SIDE OF THE CONTOUR
C
  114 DEL=DID-DJA1
      IF(ABS(DEL).LT.1.0E-4)DEL=SIGN(1.0E-4,DEL)
      YPT=FLOAT(IY)+(XLEV-DJA1)/DEL
  115 XPT=1.
      X1=XPT
      Y1=YPT
      ILOOP=0
      IF(DID-DATA(LPX))130,138,130
C
C     IF TWO ROWS ARE EACH FILLED WITH HIGH AND LOW NUMBERS RESPECTIVELY
C     WE HAVE JUST FOUND THE INTER-ROW CONTOUR. 
C
  116 DIA2=DATA(L+1)
      DJA2=UNMARK(DIA2) 
      IF(DIA1.NE.DJA1)GO TO 130
      IF (DJA1-XLEV) 129,118,128
  118 IF (DJA2.NE.XLEV) GO TO 121
      XPT=IX
      X1=XPT
      IF (IX-1) 216,120,137
  120 ILOOP=0
      GO TO 137 
  121 IF (IX-MAXX+1) 130,122,122
  122 IF (IY-1) 217,130,123
C
C     FINDING THE LOWER OR UPPER CORNER OF A PLATEAU THAT HITS THE
C     RIGHT EDGE.
C
  123 IF (DIA2-DJA2) 131,124,131
  124 L=L+1
      LPX=L+NROW
      LMX=L-NROW
      DJA2=UNMARK(DATA(LPX))
      DIA2=UNMARK(DATA(LMX))
      IF (DIA2-XLEV) 125,126,125
  125 IUP=0
      GO TO 127 
  126 IF (DJA2-XLEV) 127,131,127
C
C     IF BOTH UPPER AND LOWER POINTS ARE ALSO ON, THEN DUMP THIS ONE.
C
  127 IRT=0
      ILOOP=0
      XPT=MAXX
      X1=XPT
      GO TO 137 
C
C     SKIPS ANY BLOCK OF ALL XLEV
C     FOR 1ST POINT GREATER THAN XLEV
C
  128 IF (XLEV-DJA2) 130,130,136
C
C     FOR 1ST POINT LESS THAN XLEV
C
  129 IF (XLEV-DJA2) 136,130,130
C
C     OMITS ANY RIGHT HAND EDGE COLUMN OF ALL XLEV
C
  130 IF (IX-MAXX+1) 135,131,131
C
C     IF IX = MAXX-1 WE JUST SCANNED THE LAST COLUMN-GAP...MUST SWITCH
C     ROWS.
C
  131 IF (IY-MAXY) 133,132,218
C
C     IF IY = MAXY, WE ARE AT THE TOP ROW ALSO AND CONTOURING IS DONE.
C
  132 IDUMP=1
      GO TO 207 
  133 IX=1
      IY=IY+1
      L=IX+NROW*(IY-1)
      IF (IY-MAXY) 105,134,219
  134 JLOOP=0
      ILOOP=0
      JUP=0 
      IUP=0
C
C     FOR LAST ROW ONLY 
C
      GO TO 105 
  135 IX=IX+1
      L=L+1
      DJA1=DJA2 
      DIA1=DIA2 
      GO TO 116 
  136 DEL=DJA2-DJA1
      IF(ABS(DEL).LT.1.0E-4)DEL=SIGN(1.0E-4,DEL)
      XPT=FLOAT(IX)+(XLEV-DJA1)/DEL
      X1=XPT
C	SET LEAST SIGNIFICANT BIT
  137 DATA(L)=MARK(DATA(L))
C
C     MARK DATA(L)
C
      YPT=IY
      Y1=YPT
  138 CALL CPLOT (0,XPT,YPT,XLEV)
      IX=IX+1
      IF (IY-1) 220,139,140
  139 ILOOP=0
C
C     TO START THE NEXT SEARCH AT A NEW POINT.
C
  140 IF (IDUMP) 207,141,207
C
C     X1, Y1 ARE TO SAVE THE STARTING POINT OF THIS CONTOUR LINE.
C     ILOOP=0 IF WE START ON AN EDGE, AND 1 IF OTHERWISE.
C     CONPT FINDS THE NEXT POINT ON THE CONTOUR, WORKING FROM THE LAST
C     POINT THAT WAS FOUND AND THE GENERAL DIRECTION (UP OR DOWN, RIGHT 
C     OR LEFT) THAT THE CONTOUR IS PROGRESSING.  IT ALSO NOTES ANY
C     CHANGE IN THIS DIRECTION AND STORES IT FOR THE NEXT CALCULATION.
C
  141 IDENT=-1
      IXPT=XPT
      IYPT=YPT
      IF (YPT-FLOAT(IYPT)) 142,143,142
C
C     GO TO 143 IF IT'S ON A ROW.
C
  142 IYA=IYPT+IUP
      IYB=IYA
      IYC=IYPT+1-IUP
      IXA=IXPT
      IXB=IXPT-1+IRT+IRT
      IXC=IXB
      IIND=1
      ID=IXA+NROW*(IYC-1)
      GO TO 148 
  143 IF (XPT-FLOAT(IXPT)) 144,145,144
C
C     GO TO 145 IF IT'S ON A ROW AND ON A COLUMN.
C
  144 IXA=IXPT+IRT
      GO TO 147 
  145 IDENT=0
  146 IXA=IXPT-1+IRT+IRT
  147 IXB=IXA
      IXC=IXA+1-IRT-IRT 
      IYA=IYPT
      IYB=IYPT-1+IUP+IUP
      IYC=IYB
      IIND=2
      ID=IXC+NROW*(IYA-1)
  148 IA=IXA+NROW*(IYA-1)
      IB=IXB+NROW*(IYB-1)
      IC=IXC+NROW*(IYC-1)
      DIA=UNMARK(DATA(IA))
      DIB=UNMARK(DATA(IB))
      DIC=UNMARK(DATA(IC))
      DID=UNMARK(DATA(ID))
      DIE=(DIA+DIB+DIC+DID)*0.25
      IF (DIA-XLEV) 149,156,161 
C
C     WE TEST IN ORDER DIA,DIC,DIB AND (POSSIBLY) DIE TO DETERMINE THE
C     NEXT POINT ON THE CONTOUR.  DIE TEST IS TO TRY TO PREVENT THE
C     CONTOUR CROSSING A SADDLE OR TROUGH.
C
  149 IF (DIC-XLEV) 150,153,155 
  150 IF (DIB-XLEV) 168,151,152 
  151 IF (DIE-XLEV) 190,180,180 
  152 IF (DIE-XLEV) 190,184,184 
  153 IF (DIB-XLEV) 181,180,154 
  154 IF (DIE-XLEV) 181,184,184 
  155 IF (DIB-XLEV) 187,180,184 
  156 IF (DIC-XLEV) 157,159,160 
  157 IF (DIB-XLEV) 179,180,158 
  158 IF (DIE-XLEV) 179,187,187 
  159 IF (DIB-XLEV) 179,168,179 
  160 IF (DIB-XLEV) 158,180,179 
  161 IF (DIC-XLEV) 162,163,165 
  162 IF (DIB-XLEV) 184,180,187 
  163 IF (DIB-XLEV) 164,181,181 
  164 IF (DIE-XLEV) 184,181,181 
  165 IF (DIB-XLEV) 166,167,168 
  166 IF (DIE-XLEV) 184,190,190 
  167 IF (DIE-XLEV) 180,190,190 
  168 IF (IDENT) 190,169,170
  169 IRT=1-IRT 
      GO TO 174 
  170 IF (IDENT-2) 171,173,173
  171 IRT=1-IRT 
      IUP=1-IUP 
      GO TO 177 
  172 IDENT=IDENT+1
      GO TO 146 
  173 IDUMP=3
      GO TO 194 
C
C     X-LIMIT CHECK
C
  174 IF (IXPT+IRT-2) 176,172,175
  175 IF (IXPT+IRT-MAXX) 172,172,176
  176 IDENT=IDENT+2
      GO TO 171 
C
C     Y-LIMIT CHECK
C
  177 IF (IYPT+IUP-2) 173,172,178
  178 IF (IYPT+IUP-MAXY) 172,172,173
  179 XPT=IXA
      YPT=IYA
      GO TO 194 
  180 XPT=IXB
      YPT=IYB
      GO TO 194 
  181 XPT=IXC
      YPT=IYC
      IF (IDENT) 182,194,194
C
C     DON'T CHANGE DIRECTION IF IT'S AN ..ON-TO-ON..
C
  182 IF (IIND-1) 183,192,183
  183 IRT=1-IRT 
      GO TO 194 
  184 DEL=DIB-DIA
      IF(ABS(DEL).LT.1.0E-4)DEL=SIGN(1.0E-4,DEL)
      DEL=(XLEV-DIA)/DEL
      IF (IIND-2) 185,186,185
  185 XPT=FLOAT(IXA)+DEL*FLOAT(IRT+IRT-1)
      YPT=IYA
      GO TO 194 
  186 YPT=FLOAT(IYA)+DEL*FLOAT(IUP+IUP-1)
      XPT=IXA
      GO TO 194 
  187 DEL=DIB-DIC
      IF(ABS(DEL).LT.1.0E-4)DEL=SIGN(1.0E-4,DEL)
      DEL=(XLEV-DIC)/DEL
      IF (IIND-2) 188,189,188
  188 YPT=FLOAT(IYC)+DEL*FLOAT(IUP+IUP-1)
      XPT=IXB
      GO TO 194 
  189 XPT=FLOAT(IXC)+DEL*FLOAT(IRT+IRT-1)
      YPT=IYB
      GO TO 194 
  190 DEL=DID-DIC
      IF(ABS(DEL).LT.1.0E-4)DEL=SIGN(1.0E-4,DEL)
      DEL=(XLEV-DIC)/DEL
      IF (IIND-2) 191,193,191
  191 XPT=FLOAT(IXC)-DEL*FLOAT(IRT+IRT-1)
      YPT=IYC
  192 IUP=1-IUP 
      GO TO 194 
  193 YPT=FLOAT(IYC)-DEL*FLOAT(IUP+IUP-1)
      XPT=IXC
      GO TO 183 
  194 CALL CPLOT (1,XPT,YPT,XLEV)
C
C     PLOTS A LINE TO THIS LOCATION.
C
      IF (IDUMP-3) 195,206,206
  195 IYPT=YPT
      I=IFIX(XPT)+NROW*(IYPT-1) 
      IF(ABS(YPT-FLOAT(IYPT))-0.01)196,196,198
C
C     IF YPT IS AN INTEGER, MARK  THE POINT IMMEDI-
C     ATELY TO THE LEFT.  IF IT IS ALREADY MARKED, THIS 
C     HAS BEEN PLOTTED BEFORE.
C
  196 IF(TSTMRK(DATA(I)))202,197,202
C
C     IF THIS POINT HAS BEEN TRAVERSED BEFORE...
C
  197 DATA(I)=MARK(DATA(I))
C
C     IF NOT, KEYS THIS POINT FOR FUTURE REFERENCE.
C
  198 IF(XPT.LT.0.999999)GO TO 221
      IF(ABS(XPT-1.0).LT.0.01)GO TO 2040
      IF (XPT-FLOAT(MAXX)) 200,204,222
  200 IF(YPT.LT.0.999999)GO TO 223
      IF(ABS(YPT-1.0).LT.0.01)GO TO 2040
      IF (YPT-FLOAT(MAXY)) 141,204,224
  202 IF (XPT-X1) 204,203,204
  203 IF (YPT-Y1) 204,206,204
C
C     IF CONTOUR HAS COME FULL CIRCLE, THAT'S ALL...
C
 2040 DATA(I)=MARK(DATA(I))
  204 IF (ILOOP-1) 206,205,206
C
C     IF THE CONTOUR HITS EITHER (ANY EDGE) OR (ITSELF IN THE MIDDLE),
C     IT RETURNS ONCE TO THE STARTING POINT AND BEGINS CONTOURING IN
C     THE OTHER DIRECTION.  IF THE STARTING POINT IS ON AN EDGE, IT'S T 
C     THE END OF THE CONTOUR LINE.
C
  205 IRT=0
      IUP=0
      ILOOP=2
      XPT=X1
      YPT=Y1
      CALL CPLOT (0,XPT,YPT,XLEV)
C
C     MOVES THE PEN BACK TO STARTING POINT WITHOUT PLOTTING.
C
      GO TO 141 
  206 CALL CPLOT (2,XPT,YPT,XLEV)
C
C     TO NUMBER EACH SECTION OF A CONTOUR
C
      GO TO 103 
C
C     IDUMP=0, MORE TO CONTOUR...1, END OF MAP, NO MORE...
C     3, ONE-POINT CONTOUR.
C
  207 IF (CSTP) 208,210,208
  208 IF ((CTR+CSTP-CLIM)*SIGN(1.,CSTP)-1.E-5*ABS(CSTP)) 209,209,210
  209 CTR=CTR+CSTP
      GO TO 101 
  210 CALL CPLOT (0,0.,0.,0.)
C
C     MAKE SURE ALL POINTS HAVE BEEN PLOTTED
C
      RETURN
C
C     TRAPS TO OUTPUT A WARNING IF ANY BOUNDARY IS CROSSED.
C
  211 WRITE(6,226)IT,IX,MAXX,IY,MAXY,XPT,YPT
      STOP
  212 IT=101
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 131 
  213 IT=104
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 116 
  214 IT=105
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 111 
  215 IT=109
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 116 
  216 IT=116
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 137 
  217 IT=119
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 123 
  218 IT=128
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 132 
  219 IT=130
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 105 
  220 IT=135
      WRITE(6,226) IT,IX,MAXX,IY,MAXY
      GO TO 140 
  221 IT=303
      GO TO 211 
  222 IT=304
      GO TO 211 
  223 IT=305
      GO TO 211 
  224 IT=306
      GO TO 211 
C
  226 FORMAT(1X,5I6,2F10.3)
C
	END
