SUBROUTINE CONTUR(Z,NX,NY,X1,Y1,XL,YL,ZLEV,NDECL,LWGTL,NLEV,HGT, *NARC,X,Y,S) C GIVEN HEIGHTS Z DEFINED OVER A RECTANGULAR GRID,CONTOURS ARE DRAWN C AT SPECIFIED LEVELS. THERE ARE PROVISIONS FOR LABELING THE C CONTOURS,BLANKING REGIONS OF THE RECTANGLE,AND MAKING THE CONTOURS C MORE SMOOTH FLOWING. C C Z(I,J),I=1,NX J=1,NY IS THE INPUT ARRAY OF HEIGHTS.POINTS IN C UNDEFINED REGIONS OF Z SHOULD BE SET TO 1.E35 C X1,Y1,XL,YL ARE THE LOWER LEFT AND UPPER RIGHT CORNERS OF C THE GRID ON THE PLOTTER IN INCHES. C ZLEV(K),K=1,NLEV IS THE ARRAY OF LEVELS TO BE CONTOURED C NDECL(K),K=1,NLEV GIVES THE NUMBER OF DECIMAL PLACES IN THE C CONTOUR LABELS. -1=NO DECIMAL. -2 OR LESS=NO LABEL. C LWGTL(K),K=1,NLEV GIVES THE WEIGHT OF THE CONTOUR LINE. C 1 OR LESS=STANDARD LINE. 2=HEAVY LINE. 3=DOTTED LINE. C HGT = HEIGHT OF CHARACTERS OF LABELS IN INCHES. C NARC = 1,2,3 ... 10 . AN ARC OF NARC SUBSEGMENTS REPLACES EACH C ST LINE SEGMENT OF A CONTOUR. THE ARC WILL MATCH SLOPES WITH C THE ADJACENT ARCS. CARE SHOULD BE TAKEN HERE AS OVERLAPPING OF C CONTOURS IS POSSIBLE WHEN NARC IS USED. NARC=1 HAS NO EFFECT. C IND10 SHOULD BE DIMENSIONED AT LEAST NX BY NY/10 ROUNDED UP. C KAB SHOULD BE DIMENSIONED AT LEAST NX/10+2 BY NY/10+2 . C C OCEANOGRAPHY EMR DEC/69 IMPLICIT REAL*8 (A-H, O-Z) REAL*4 X(1), Y(1), S(1) REAL*4 X1,Y1,XL,YL,HGT REAL*4 HHGT,XLAB,YLAB,XENDL,YENDL REAL*4 ANGLE, ZC REAL*4 XC(51), YC(51), XARC(4),YARC(4) REAL*8 XX(4), YY(4), ZZ(4) REAL*8 ZMAXB(10,10), ZMINB(10,10) REAL*4 Z(57,57) INTEGER KAB(07,07) REAL*4 ZLEV(NLEV) INTEGER NDECL(NLEV),LWGTL(NLEV),IXPON(10) INTEGER ISIN(4),IND(4),LUSED(4),KABOV(4) COMMON /XYS/ IND10(57,6) C C INITIALIZE C*********************************************************************** C LP=6 IF(NLEV)1200,1200,12 12 PI=3.1415926 DSLAB=10.0 SLAB1F=.2 ISIN(1)=0 ISIN(2)=1 ISIN(3)=0 ISIN(4)=-1 BIG=.9E35 NXM1=NX-1 NYM1=NY-1 DX= (DBLE(XL)-DBLE(X1))/NXM1 DY= (DBLE(YL)-DBLE(Y1))/NYM1 HHGT = HGT /2. NY10=(NY-1)/10+1 DO 17 J2BIT=1,10 17 IXPON(J2BIT)=4**(J2BIT-1) C C GET ZMAX,ZMIN,ZRNG,WEEZEE C*********************************************************************** C IBLK=MAX(NX/10+1,5) JBLK=MAX(NY/10+1,5) ZMAX=-1.E35 ZMIN=1.E35 IB=0 DO 55 ISTRT=1,NX,IBLK IB=IB+1 IEND=MIN(ISTRT+IBLK,NX) JB=0 DO 55 JSTRT=1,NY,JBLK JB=JB+1 JEND=MIN(JSTRT+JBLK,NY) ZMAXBL=-1.E35 ZMINBL=1.E35 DO 50 I=ISTRT,IEND DO 50 J=JSTRT,JEND ZIJ=Z(I,J) IF(ZIJ.GT.BIG) GO TO 50 ZMAXBL=MAX(ZIJ,ZMAXBL) ZMINBL=MIN(ZMINBL,ZIJ) 50 CONTINUE ZMAX=MAX(ZMAXBL,ZMAX) ZMIN=MIN(ZMINBL,ZMIN) ZMAXB(IB,JB)=ZMAXBL ZMINB(IB,JB)=ZMINBL 55 CONTINUE ZRNG=(ZMAX-ZMIN)*1.1 IF(ZRNG)1200,1200,60 60 WEEZEE = ZRNG*.0002 C C MAIN LOOP OVER ALL CONTOUR LEVELS C*********************************************************************** C DO 1100 LEV=1,NLEV ZC = ZLEV(LEV) IF((ZMAX-ZC)*(ZC-ZMIN))1100,1100,70 70 NDEC = NDECL(LEV) LWGT = LWGTL(LEV) DO 75 I=1,NX DO 75 J=1,NY10 75 IND10(I,J)=0 C C GET LABEL WIDTH C*********************************************************************** C NCHAR = NDEC+2 IF(NCHAR)120,120,80 80 IF(ZC)90,100,100 90 NCHAR=NCHAR+1 100 ABZRND = ABS(ZC)+.5/10.**MAX(0,NDEC) DO 110 K=1,10 IF(ABZRND-10**K)120,110,110 110 NCHAR=NCHAR+1 120 WIDTH= HGT *(NCHAR*7.+5.)/7. C120 WIDTH= HGT *(NCHAR*6.+5.)/7. C C LOOP OVER THE GRID SEARCHING FOR A CONTOUR ENTERING SQUARE II,JJ C THRU SIDE 1 OR 4. (THE CONTOUR MUST BE UNUSED) C*********************************************************************** C IB=0 DO 1010 ISTRT=1,NX,IBLK IEND=MIN(ISTRT+IBLK,NX+1) IB=IB+1 JB=0 DO 1010 JSTRT=1,NY,JBLK JB=JB+1 IF(ZMAXB(IB,JB).LT.ZC .OR. ZMINB(IB,JB).GT.ZC) GO TO 1010 JEND=MIN(JSTRT+JBLK,NY+1) DO 200 I=ISTRT,IEND IKAB=I-ISTRT+1 DO 200 J=JSTRT,JEND JKAB=J-JSTRT+1 KABIJ=10 IF(I.GT.NX .OR. J.GT.NY) GO TO 200 ZIJ=Z(I,J) IF(ZIJ.GT.BIG) GO TO 200 IF(ABS(ZIJ-ZC).LE.WEEZEE) ZIJ=ZC+SIGN(WEEZEE,ZIJ-ZC) KABIJ=-1 IF(ZIJ.GT.ZC) KABIJ=1 200 KAB(IKAB,JKAB)=KABIJ IENDM1=IEND-1 JENDM1=JEND-1 DO 1000 II=ISTRT,IENDM1 IKAB=II-ISTRT+1 DO 1000 JJ=JSTRT,JENDM1 JKAB=JJ-JSTRT+1 I=II J=JJ KABIJ=KAB(IKAB,JKAB) KABIP1=KAB(IKAB+1,JKAB) KABJP1=KAB(IKAB,JKAB+1) JWRD=(J-1)/10+1 J2BIT=J-(JWRD-1)*10 IND1=IND10(I,JWRD)/IXPON(J2BIT) IND1=IND1-(IND1/4)*4 LUSED4=IND1/2 LUSED1=IND1-2*LUSED4 IRETRN=300 LL=1 IF(KABIJ+KABIP1+LUSED1.EQ.0) GO TO 540 LL=4 IF(KABIJ+KABJP1+LUSED4.EQ.0) GO TO 540 GO TO 1000 300 CONTINUE 340 LIN = LL K=1 NSEG = 1 C C GIVEN ENTRANCE TO SQUARE(I,J) ON SIDE LIN, RECORD THE ENTRANCE C POINT X(K),Y(K). SET LIN TO USED. C*********************************************************************** C 350 LP1 = LIN+1-(LIN/4)*4 LP2 = LP1+1-(LP1/4)*4 LM1 = LP2+1-(LP2/4)*4 IF(IRETRN-360)355,360,360 355 K=K+1 FRAC = (ZC-ZZ(LIN))/(ZZ(LP1)-ZZ(LIN)) X(K) = XX(LIN)+ (XX(LP1)-XX(LIN))*FRAC Y(K) = YY(LIN)+ (YY(LP1)-YY(LIN))*FRAC IIND=I+LP2/4 JIND=J+LP1/4 LHOR = LIN-(LIN/2)*2 JWRD=(JIND-1)/10+1 J2BIT=JIND-(JWRD-1)*10 IND10(IIND,JWRD)=IND10(IIND,JWRD)+(2-LHOR)*IXPON(J2BIT) C C SEE IF AN EXIT EXISTS ON SIDE L-1,L+1,OR L+2. IF SO CHOOSE THE ONE C CLOSEST TO SIDE L. IF THE EXIT IS ALREADY USED TERMINATE X,Y. C*********************************************************************** C 360 IRETRN = 350 LEXIT = LM1 IF(KABOV(LIN)+KABOV(LM1))380,370,380 370 IF(KABOV(LP1)+KABOV(LP2))450,390,450 380 IF(KABOV(LP1)+KABOV(LP2))410,400,410 390 FLM1 = (ZC-ZZ(LIN))/(ZZ(LM1)-ZC) FLP1 = (ZC-ZZ(LP1))/(ZZ(LP2)-ZC) IF(LUSED(LP1).EQ.1) GO TO 450 IF(FLM1.LE.FLP1 .AND. LUSED(LM1).EQ.0) GO TO 450 400 LEXIT = LP1 GO TO 450 410 LEXIT = LP2 IF(KABOV(LP2)+KABOV(LM1))470,450,470 450 IF(LUSED(LEXIT))530,530,460 460 KMAX=K+1 X(KMAX)=X(2) Y(KMAX)=Y(2) X(1)=X(K) Y(1)=Y(K) X(K+2)=X(3) Y(K+2)=Y(3) KLOSED=1 GO TO 700 C C NO EXIT. IF NSEG=1 REVERSE X,Y AND CONTINUE. C IF NSEG=2 TERMINATE X,Y. C*********************************************************************** C 470 IF(KABOV(LP2)+KABOV(LM1)-15)480,480,500 480 KDA=LIN KDB=LP2 IF(KABOV(LP2)-5)495,495,490 490 KDA = LM1 KDB = LP1 495 K=K+1 FRAC = (ZC-ZZ(KDA))/(ZZ(KDB)-ZZ(KDA)) X(K) = XX(KDA) + (XX(KDB)-XX(KDA))*FRAC Y(K) = YY(KDA) + (YY(KDB)-YY(KDA))*FRAC 500 IF(NSEG-1)510,510,505 505 KMAX=K X(1)=X(2) Y(1)=Y(2) X(K+1)=X(K) Y(K+1)=Y(K) KLOSED=0 GO TO 700 510 IRETRN=360 NSEG = 2 KH = 1+K/2 DO 520 KK=2,KH KKR = K+2-KK XKK=X(KK) YKK=Y(KK) X(KK)= X(KKR) Y(KK)= Y(KKR) X(KKR)=XKK 520 Y(KKR)=YKK I=II J=JJ LEXIT=LL C C FIND SQUARE ENTERED BY PRESENT EXIT. GET XX,YY,ZZ,KABOV FOR EACH C CORNER AND LUSED FOR EACH SIDE. C*********************************************************************** C 530 I= I+ ISIN(LEXIT) JSUB=5-LEXIT J= J +ISIN(JSUB) LIN=LEXIT+2 - ((LEXIT+1)/4)*4 C C 540 DO 620 L=1,4 JL=J+(L-1)/2 LP1=L+1-(L/4)*4 IL=I+(LP1-1)/2 IND(L)=0 ZZ(L)=1.E35 KABOV(L)=10 IF(IL*(NX+1-IL))610,610,550 550 IF(JL*(NY+1-JL))610,610,560 560 IF(Z(IL,JL)-BIG)570,610,610 570 JWRD=(JL-1)/10+1 J2BIT=JL-(JWRD-1)*10 INDL=IND10(IL,JWRD)/IXPON(J2BIT) IND(L)=INDL-(INDL/4)*4 ZZ(L)=Z(IL,JL) IF(ABS(ZZ(L)-ZC)-WEEZEE)580,580,590 580 ZZ(L)=ZC+SIGN(WEEZEE,ZZ(L)-ZC) 590 KABOV(L)=-1 IF(ZZ(L)-ZC)610,610,600 600 KABOV(L)=1 610 XX(L)= X1 + (IL-1)*DX 620 YY(L)= Y1 + (JL-1)*DY LUSED(4) = IND(1)/2 LUSED(1) = IND(1)-2*LUSED(4) LUSED(2) = IND(2)/2 LUSED(3) = IND(4)-(IND(4)/2)*2 IF(IRETRN-350) 300,350,350 C C THE ARRAYS X,Y ARE NOW COMPLETE. CALCULATE DISTANCE S ALONG C THE CONTOUR. START PLOTTING THE CONTOUR C*********************************************************************** C 700 S(2)=0. IF(KMAX-2)1000,1000,720 720 CONTINUE DO 750 K=3,KMAX DXX = X(K)-X(K-1) DYY = Y(K)-Y(K-1) 750 S(K) = S(K-1)+ SQRT(DXX*DXX+DYY*DYY) SMAX = S(KMAX) SLAB1 = SMAX*SLAB1F STEST = DSLAB-SLAB1 K=2 CALL PLOTNY(REAL(X(2)),REAL(Y(2)),3,0) C C CHECK CONDITIONS FOR LABELLING. C*********************************************************************** C 755 IF(NDEC+2)900,900,760 760 KM1= MAX(K-1,2) STEST = STEST + S(K)-S(KM1) IF(STEST-DSLAB)900,770,770 770 IF(SMAX-S(K)-2.0*WIDTH)900,900,780 780 KP1=K+1 IF(LEV.EQ.1) GO TO 785 DLEV=ABS(ZLEV(LEV)-ZLEV(LEV-1)) I=(X(K)-X1)/DX+1.01 I=MIN(I,NX-1) J=(Y(K)-Y1)/DY+1.01 J=MIN(J,NY-1) DZDX=(Z(I+1,J)-Z(I,J))/DX DZDY=(Z(I,J+1)-Z(I,J))/DY DZDG=SQRT(DZDX*DZDX+DZDY*DZDY) IF(DZDG.EQ.0.) GO TO 785 CSPACE=DLEV/DZDG IF(CSPACE.LT.HGT*.90) GO TO 900 785 CONTINUE DO 800 KK = KP1,KMAX DXX = X(KK)-X(K) DYY = Y(KK)-Y(K) SPACE = SQRT(DXX*DXX+DYY*DYY) ARK =S(KK)-S(K) IF(SPACE-WIDTH)800,790,790 790 IF(SPACE/ARK-.95)900,810,810 800 CONTINUE GO TO 900 C C DRAW THE LABEL C*********************************************************************** C 810 STEST=0. DO 812 KKK=1,5 DXX=X(KK-1)+(X(KK)-X(KK-1))*.2*KKK-X(K) DYY=Y(KK-1)+(Y(KK)-Y(KK-1))*.2*KKK-Y(K) SPACE=SQRT(DXX*DXX+DYY*DYY) IF(SPACE-WIDTH)812,812,815 812 CONTINUE 815 CONTINUE XENDL = X(K) + WIDTH*DXX/SPACE YENDL = Y(K) + WIDTH*DYY/SPACE ANGLE = 90. IF(DXX)820,830,820 820 ANGLE = ATAN(DYY/DXX)*180./PI 830 IF(DXX) 850,840,840 840 XLAB = X(K)+ HHGT*(DXX+DYY)/SPACE YLAB = Y(K)+ HHGT*(DYY-DXX)/SPACE GO TO 860 850 XLAB = XENDL-HHGT*(DXX+DYY)/SPACE YLAB = YENDL-HHGT*(DYY-DXX)/SPACE 860 CALL NUMBER(XLAB,YLAB,HGT,ZC,ANGLE,NDEC) CALL PLOTNY(XENDL,YENDL,3,0) CALL PLOTNY(REAL(X(KK)),REAL(Y(KK)),2,LWGT) K=KK GO TO 920 C C PLOT THE SEGMENT FROM XK,YK TO X(K+1),Y(K+1). C*********************************************************************** C 900 CONTINUE DO 905 IARC=1,4 XARC(IARC)=X(K+IARC-2) YARC(IARC)=Y(K+IARC-2) 905 CONTINUE CALL ARC(XARC,YARC,NARC+1,.001,XC,YC,NPT) DO 910 KK = 2,NPT 910 CALL PLOTNY(XC(KK),YC(KK),2,LWGT) K=K+1 920 IF(K-KMAX)755,1000,1000 1000 CONTINUE 1010 CONTINUE C 1100 CONTINUE 1200 RETURN END