Initial commit

1 files changed, 830 insertions, 0 deletions

diff --git a/sys/gio/ncarutil/hafton.f b/sys/gio/ncarutil/hafton.f
new file mode 100644
index 00000000..7d597470
--- /dev/null
+++ b/sys/gio/ncarutil/hafton.f
@@ -0,0 +1,830 @@
+      SUBROUTINE HAFTON (Z,L,M,N,FLO,HI,NLEV,NOPT,NPRM,ISPV,SPVAL)
+C
+C                                                                               
+C +-----------------------------------------------------------------+           
+C |                                                                 |           
+C |                Copyright (C) 1986 by UCAR                       |           
+C |        University Corporation for Atmospheric Research          |           
+C |                    All Rights Reserved                          |           
+C |                                                                 |           
+C |                 NCARGRAPHICS  Version 1.00                      |           
+C |                                                                 |           
+C +-----------------------------------------------------------------+           
+C                                                                               
+C                                                                               
+C     SUBROUTINE HAFTON (Z,L,M,N,FLO,HI,NLEV,NOPT,NPRM,ISPV,SPVAL)
+C
+C
+C     DIMENSION OF           Z(L,M)
+C     ARGUMENTS
+C
+C     LATEST REVISION        JULY,1984
+C
+C     PURPOSE                HAFTON DRAWS A HALF-TONE PICTURE FROM DATA
+C                            STORED IN A RECTANGULAR ARRAY WITH THE
+C                            INTENSITY IN THE PICTURE PROPORTIONAL TO
+C                            THE DATA VALUE.
+C
+C     USAGE                  IF THE FOLLOWING ASSUMPTIONS ARE MET, USE
+C
+C                                 CALL EZHFTN (Z,M,N)
+C
+C                              ASSUMPTIONS:
+C                                  .ALL OF THE ARRAY IS TO BE DRAWN.
+C                                  .LOWEST VALUE IN Z WILL BE AT LOWEST
+C                                   INTENSITY ON READER/PRINTER OUTPUT.
+C                                  .HIGHEST VALUE IN Z WILL BE AT
+C                                   HIGHEST INTENSITY.
+C                                  .VALUES IN BETWEEN WILL APPEAR
+C                                   LINEARLY SPACED.
+C                                  .MAXIMUM POSSIBLE NUMBER OF
+C                                   INTENSITIES ARE USED.
+C                                  .THE PICTURE WILL HAVE A PERIMETER
+C                                   DRAWN.
+C                                  .FRAME WILL BE CALLED AFTER THE
+C                                   PICTURE IS DRAWN.
+C                                  .Z IS FILLED WITH NUMBERS THAT SHOULD
+C                                   BE USED (NO MISSING VALUES).
+C
+C                            IF THESE ASSUMPTIONS ARE NOT MET, USE
+C
+C                                 CALL HAFTON (Z,L,M,N,FLO,HI,NLEV,
+C                                              NOPT,NPRM,ISPV,SPVAL)
+C
+C     ARGUMENTS
+C
+C     ON INPUT               Z
+C     FOR EZHFTN               M BY N ARRAY TO BE USED TO GENERATE A
+C                              HALF-TONE PLOT.
+C
+C                            M
+C                              FIRST DIMENSION OF Z.
+C
+C                            N
+C                              SECOND DIMENSION OF Z.
+C
+C     ON OUTPUT              ALL ARGUMENTS ARE UNCHANGED.
+C     FOR EZHFTN
+C
+C     ON INPUT               Z
+C     FOR HAFTON               THE ORIGIN OF THE ARRAY TO BE PLOTTED.
+C
+C                            L
+C                              THE FIRST DIMENSION OF Z IN THE CALLING
+C                              PROGRAM.
+C
+C                            M
+C                              THE NUMBER OF DATA VALUES TO BE PLOTTED
+C                              IN THE X-DIRECTION (THE FIRST SUBSCRIPT
+C                              DIRECTION).  WHEN PLOTTING ALL OF AN
+C                              ARRAY, L = M.
+C
+C                            N
+C                              THE NUMBER OF DATA VALUES TO BE PLOTTED
+C                              IN THE Y-DIRECTION (THE SECOND SUBSCRIPT
+C                              DIRECTION).
+C
+C                            FLO
+C                              THE VALUE OF Z THAT CORRESPONDS TO THE
+C                              LOWEST INTENSITY.  (WHEN NOPT.LT.0, FLO
+C                              CORRESPONDS TO THE HIGHEST INTENSITY.)
+C                              IF FLO=HI=0.0, MIN(Z) WILL BE USED FOR FLO.
+C
+C                            HI
+C                              THE VALUE OF Z THAT CORRESPONDS TO THE
+C                              HIGHEST INTENSITY.  (WHEN NOPT.LT.0, HI
+C                              CORRESPONDS TO THE LOWEST INTENSITY.)  IF
+C                              HI=FLO=0.0, MAX(Z) WILL BE USED FOR HI.
+C
+C                            NLEV
+C                              THE NUMBER OF INTENSITY LEVELS DESIRED.
+C                              16 MAXIMUM.  IF NLEV = 0 OR 1, 16 LEVELS
+C                              ARE USED.
+C
+C                            NOPT
+C                              FLAG TO CONTROL THE MAPPING OF Z ONTO THE
+C                              INTENSITIES.  THE SIGN OF NOPT CONTROLS
+C                              THE DIRECTNESS OR INVERSENESS OF THE
+C                              MAPPING.
+C
+C                              . NOPT POSITIVE YIELDS DIRECT MAPPING.
+C                                THE LARGEST VALUE OF Z PRODUCES THE
+C                                MOST DENSE DOTS.  ON MECHANICAL PLOTTERS,
+C                                LARGE VALUES OF Z WILL PRODUCE A DARK
+C                                AREA ON THE PAPER.  WITH THE FILM
+C                                DEVELOPMENT METHODS USED AT NCAR,
+C                                LARGE VALUES OF Z WILL PRODUCE MANY
+C                                (WHITE) DOTS ON THE FILM, ALSO
+C                                RESULTING IN A DARK AREA ON
+C                                READER-PRINTER PAPER.
+C                              . NOPT NEGATIVE YIELDS INVERSE MAPPING.
+C                                THE SMALLEST VALUES OF Z PRODUCE THE
+C                                MOST DENSE DOTS RESULTING IN DARK
+C                                AREAS ON THE PAPER.
+C
+C                              THE ABSOLUTE VALUE OF NOPT DETERMINES THE
+C                              MAPPING OF Z ONTO THE INTENSITIES.  FOR
+C                              IABS(NOPT)
+C                              = 0  THE MAPPING IS LINEAR.  FOR
+C                                   EACH INTENSITY THERE IS AN EQUAL
+C                                   RANGE IN Z VALUE.
+C                              = 1  THE MAPPING IS LINEAR.  FOR
+C                                   EACH INTENSITY THERE IS AN EQUAL
+C                                   RANGE IN Z VALUE.
+C                              = 2  THE MAPPING IS EXPONENTIAL.  FOR
+C                                   LARGER VALUES OF Z, THERE IS A
+C                                   LARGER DIFFERENCE IN INTENSITY FOR
+C                                   RELATIVELY CLOSE VALUES OF Z.  DETAILS
+C                                   IN THE LARGER VALUES OF Z ARE DISPLAYED
+C                                   AT THE EXPENSE OF THE SMALLER VALUES
+C                                   OF Z.
+C                              = 3  THE MAPPING IS LOGRITHMIC, SO
+C                                   DETAILS OF SMALLER VALUES OF Z ARE SHOWN
+C                                   AT THE EXPENSE OF LARGER VALUES OF Z.
+C                              = 4  SINUSOIDAL MAPPING, SO MID-RANGE VALUES
+C                                   OF Z SHOW DETAILS AT THE EXPENSE OF
+C                                   EXTREME VALUES OF Z.
+C                              = 5  ARCSINE MAPPING, SO EXTREME VALUES OF
+C                                   Z ARE SHOWN AT THE EXPENSE OF MID-RANGE
+C                                   VALUES OF Z.
+C
+C                            NPRM
+C                              FLAG TO CONTROL THE DRAWING OF A
+C                              PERIMETER AROUND THE HALF-TONE PICTURE.
+C
+C                              . NPRM=0:  THE PERIMETER IS DRAWN WITH
+C                                TICKS POINTING AT DATA LOCATIONS.
+C                                (SIDE LENGTHS ARE PROPORTIONAL TO NUMBER
+C                                OF DATA VALUES.)
+C                              . NPRM POSITIVE:  NO PERIMETER IS DRAWN.  THE
+C                                PICTURE FILLS THE FRAME.
+C                              . NPRM NEGATIVE:  THE PICTURE IS WITHIN THE
+C                                CONFINES OF THE USER'S CURRENT VIEWPORT
+C                                SETTING.
+C
+C                            ISPV
+C                              FLAG TO TELL IF THE SPECIAL VALUE FEATURE
+C                              IS BEING USED.  THE SPECIAL VALUE FEATURE
+C                              IS USED TO MARK AREAS WHERE THE DATA IS
+C                              NOT KNOWN OR HOLES ARE WANTED IN THE
+C                              PICTURE.
+C
+C                              . ISPV = 0:  SPECIAL VALUE FEATURE NOT IN
+C                                USE.  SPVAL IS IGNORED.
+C                              . ISPV NON-ZERO:  SPECIAL VALUE FEATURE
+C                                IN USE.  SPVAL DEFINES THE SPECIAL
+C                                VALUE.  WHERE Z CONTAINS THE SPECIAL
+C                                VALUE, NO HALF-TONE IS DRAWN.  IF ISPV
+C                                = 0  SPECIAL VALUE FEATURE NOT IN USE.
+C                                     SPVAL IS IGNORED.
+C                                = 1  NOTHING IS DRAWN IN SPECIAL VALUE
+C                                     AREA.
+C                                = 2  CONTIGUOUS SPECIAL VALUE AREAS ARE
+C                                     SURROUNDED BY A POLYGONAL LINE.
+C                                = 3  SPECIAL VALUE AREAS ARE FILLED
+C                                     WITH X(S).
+C                                = 4  SPECIAL VALUE AREAS ARE FILLED IN
+C                                     WITH THE HIGHEST INTENSITY.
+C
+C                            SPVAL
+C                              THE VALUE USED IN Z TO DENOTE MISSING
+C                              VALUES.  THIS ARGUMENT IS IGNORED IF
+C                              ISPV = 0.
+C
+C     ON OUTPUT              ALL ARGUMENTS ARE UNCHANGED.
+C     FOR HAFTON
+C
+C     NOTE                   THIS ROUTINE PRODUCES A HUGE NUMBER OF
+C                            PLOTTER INSTRUCTIONS PER PICTURE, AVERAGING
+C                            OVER 100,000 LINE-DRAWS PER FRAME WHEN M = N.
+C
+C
+C     ENTRY POINTS           EZHFTN, HAFTON, ZLSET, GRAY, BOUND, HFINIT
+C
+C     COMMON BLOCKS          HAFT01, HAFT02, HAFT03, HAFT04
+C
+C     REQUIRED LIBRARY       GRIDAL, THE ERPRT77 PACKAGE AND THE SPPS.
+C     ROUTINES
+C
+C     I/O                    PLOTS HALF-TONE PICTURE.
+C
+C     PRECISION              SINGLE
+C
+C     LANGUAGE               FORTRAN
+C
+C     HISTORY                REWRITE OF PHOMAP ORIGINALLY WRITTEN BY
+C                            M. PERRY OF HIGH ALTITUDE OBSERVATORY,
+C                            NCAR.
+C
+C     ALGORITHM              BI-LINEAR INTERPOLATION ON PLOTTER
+C                            (RESOLUTION-LIMITED) GRID OF NORMALIZED
+C                            REPRESENTATION OF DATA.
+C
+C     PORTABILITY            ANSI FORTRAN 77.
+C
+C
+C
+C INTERNAL PARAMTERSS
+C VALUES SET IN BLOCK DATA
+C         NAME  DEFAULT                   FUNCTION
+C         ----  -------                   ________
+C
+C         XLT     0.1    LEFT-HAND EDGE OF THE PLOT WHEN NSET=0.  (0.0=
+C                        LEFT EDGE OF FRAME, 1.0=RIGHT EDGE OF FRAME.)
+C         YBT     0.1    BOTTOM EDGE OF THE PLOT WHEN NSET=0.  (0.0=
+C                        BOTTOM OF FRAME, 1.0=TOP OF FRAME.)
+C         SIDE    0.8    LENGTH OF LONGER EDGE OF PLOT (SEE ALSO EXT).
+C         EXT     .25    LENGTHS OF THE SIDES OF THE PLOT ARE PROPOR-
+C                        TIONAL TO M AND N (WHEN NSET=0) EXCEPT IN
+C                        EXTREME CASES, NAMELY, WHEN MIN(M,N)/MAX(M,N)
+C                        IS LESS THAN EXT.  THEN A SQUARE PLOT IS PRO-
+C                        DUCED.  WHEN A RECTANGULAR PLOT IS PRODUCED,
+C                        THE PLOT IS CENTERED ON THE FRAME (AS LONG AS
+C                        SIDE+2*XLT = SIDE+2*YBT=1., AS WITH THE
+C                        DEFAULTS.)
+C         ALPHA   1.6    A PARAMETER TO CONTROL THE EXTREMENESS OF THE
+C                        MAPPING FUNCTION SPECIFIED BY NOPT.  (FOR
+C                        IABS(NOPT)=0 OR 1, THE MAPPING FUNCTION IS
+C                        LINEAR AND INDEPENDENT OF ALPHA.)  FOR THE NON-
+C                        LINEAR MAPPING FUNCTIONS, WHEN ALPHA IS CHANGED
+C                        TO A NUMBER CLOSER TO 1., THE MAPPING FUNCTION
+C                        BECOMES MORE LINEAR; WHEN ALPHA IS CHANGED TO
+C                        A LARGER NUMBER, THE MAPPING FUNCTION BECOMES
+C                        MORE EXTREME.
+C         MXLEV   16     MAXIMUM NUMBER OF LEVELS.  LIMITED BY PLOTTER.
+C         NCRTG   8      NUMBER OF CRT UNITS PER GRAY-SCALE CELL.
+C                        LIMITED BY PLOTTER.
+C         NCRTF   1024   NUMBER OF PLOTTER ADDRESS UNITS PER FRAME.
+C         IL    (BELOW)  AN ARRAY DEFINING WHICH OF THE AVAILABLE IN-
+C                        TENSITIES ARE USED WHEN LESS THAN THE MAXIMUM
+C                        NUMBER OF INTENSITIES ARE REQUESTED.
+C
+C
+C   NLEV         INTENSITIES USED
+C   ____         ________________
+C    2           5,11,
+C    3           4, 8,12,
+C    4           3, 6,10,13,
+C    5           2, 5, 8,11,14,
+C    6           1, 4, 7, 9,12,15,
+C    7           1, 4, 6, 8,10,12,15,
+C    8           1, 3, 5, 7, 9,11,13,15,
+C    9           1, 3, 4, 6, 8,10,12,13,15
+C   10           1, 3, 4, 6, 7, 9,10,12,13,15,
+C   11           1, 2, 3, 5, 6, 8,10,11,13,14,15,
+C   12           1, 2, 3, 5, 6, 7, 9,10,11,13,14,15,
+C   13           1, 2, 3, 4, 6, 7, 8, 9,10,12,13,14,15
+C   14           1, 2, 3, 4, 5, 6, 7, 9,10,11,12,13,14,15,
+C   15           1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,
+C   16           0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15
+C
+C
+
+      SAVE
+      DIMENSION       Z(L,N)      ,PX(2)       ,PY(2)
+      DIMENSION       ZLEV(16)    ,VWPRT(4)    ,WNDW(4)
+      DIMENSION                    VWPR2(4)    ,WND2(4)
+      CHARACTER*11    IDUMMY
+C
+C
+      COMMON /HAFTO1/ I          ,J          ,INTEN
+      COMMON /HAFTO2/ GLO        ,HA         ,NOPTN      ,ALPHA      ,
+     1                NSPV       ,SP         ,ICNST
+      COMMON /HAFTO3/ XLT        ,YBT        ,SIDE       ,EXT        ,
+     1                IOFFM      ,ALPH       ,MXLEV      ,NCRTG      ,
+     2                NCRTF      ,IL(135)
+      COMMON /HAFTO4/ NPTMAX     ,NPOINT     ,XPNT(50)  ,YPNT(50)
+C +NOAO - Blockdata rewritten as run time initialization subroutine
+C
+C     EXTERNAL HFINIT
+      call hfinit
+C -NOAO
+C
+C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
+C
+      CALL Q8QST4 ('GRAPHX','HAFTON','HAFTON','VERSION  1')
+C
+      NPOINT = 0
+      ALPHA = ALPH
+      GLO = FLO
+      HA = HI
+      NLEVL = MIN0(IABS(NLEV),MXLEV)
+      IF (NLEVL .LE. 1) NLEVL = MXLEV
+      NOPTN = NOPT
+      IF (NOPTN .EQ. 0) NOPTN = 1
+      NPRIM = NPRM
+      NSPV = MAX0(MIN0(ISPV,4),0)
+      IF (NSPV .NE. 0) SP = SPVAL
+      MX = L
+      NX = M
+      NY = N
+      CRTF = NCRTF
+      MSPV = 0
+C
+C SET INTENSITY BOUNDARY LEVELS
+C
+      CALL ZLSET (Z,MX,NX,NY,ZLEV,NLEVL)
+C
+C SET UP PERIMETER
+C
+      X3 = NX
+      Y3 = NY
+      CALL GQCNTN (IERR,NTORIG)
+      CALL GETUSV('LS',IOLLS)
+      IF (NPRIM.LT.0) THEN
+         CALL GQNT (NTORIG,IERR,WNDW,VWPRT)
+         X1 = VWPRT(1)
+         X2 = VWPRT(2)
+         Y1 = VWPRT(3)
+         Y2 = VWPRT(4)
+      ELSE IF (NPRIM.EQ.0) THEN
+         X1 = XLT
+         X2 = XLT+SIDE
+         Y1 = YBT
+         Y2 = YBT+SIDE
+         IF (AMIN1(X3,Y3)/AMAX1(X3,Y3) .GE. EXT) THEN
+            IF (NX-NY.LT.0) THEN
+               X2 =SIDE*X3/Y3+XLT
+               X2 = (AINT(X2*CRTF/FLOAT(NCRTG))*FLOAT(NCRTG))/CRTF
+            ELSE IF (NX-NY.GT.0) THEN
+               Y2 = SIDE*Y3/X3+YBT
+               Y2 = (AINT(Y2*CRTF/FLOAT(NCRTG))*FLOAT(NCRTG))/CRTF
+            END IF
+         END IF
+      ELSE IF (NPRIM.GT.0) THEN
+         X1 = 0.0
+         X2 = 1.0
+         Y1 = 0.0
+         Y2 = 1.0
+      END IF
+      MX1 = X1*CRTF
+      MX2 = X2*CRTF
+      MY1 = Y1*CRTF
+      MY2 = Y2*CRTF
+      IF (NPRIM.GT.0) THEN
+         MX1 = 1
+         MY1 = 1
+         MX2 = NCRTF
+         MY2 = NCRTF
+      END IF
+C
+C SAVE NORMALIZATION TRANS 1
+C
+      CALL GQNT (1,IERR,WNDW,VWPRT)
+C
+C DEFINE NORMALIZATION TRANS 1 AND LOG SCALING FOR USE WITH PERIM
+C DRAW PERIMETER IF NPRIM EQUALS 0
+C
+      CALL SET(X1,X2,Y1,Y2,1.0,X3,1.0,Y3,1)
+      IF (NPRIM .EQ. 0) CALL PERIM (NX-1,1,NY-1,1)
+      IF (ICNST .NE. 0) THEN
+        CALL GSELNT (0)
+        CALL WTSTR(XLT*1.1,0.5,'CONSTANT FIELD',2,0,0)
+        GO TO 132
+      END IF
+C
+C FIND OFFSET FOR REFERENCE TO IL, WHICH IS TRIANGULAR
+C
+      IOFFST = NLEVL*((NLEVL-1)/2)+MOD(NLEVL-1,2)*(NLEVL/2)-1
+C
+C OUTPUT INTENSITY SCALE
+C
+      IF (NPRIM .GT. 0) GO TO 112
+      LEV = 0
+      KX = (1.1*XLT+SIDE)*CRTF
+      KY = YBT*CRTF
+      NNX = KX/NCRTG
+  109 LEV = LEV+1
+C +NOAO
+C The following statement moved from after statement label 111 (CONTINUE) to
+C here.  Otherwise an extra (unlabelled) grayscale box was being drawn.
+C This was (eventually) causing a [floating operand error] on a Sun-3.
+      IF (LEV .GT. NLEVL) GO TO 112
+C -NOAO
+      ISUB = IOFFST+LEV
+      INTEN = IL(ISUB)
+      IF (NOPTN .LT. 0) INTEN = MXLEV-INTEN
+      NNY = KY/NCRTG
+      DO 111 JJ=1,3
+         DO 110 II=1,10
+            I = NNX+II
+            J = NNY+JJ
+            CALL GRAY
+  110    CONTINUE
+  111 CONTINUE
+C +NOAO - FTN internal write rewritten as call to encode.
+C       WRITE(IDUMMY,'(G11.4)') ZLEV(LEV)
+	call encode (11, '(g11.4)', idummy, zlev(lev))
+C -NOAO
+      TKX = KX
+      TKY = KY+38
+      CALL GQNT(1,IERR,WND2,VWPR2)
+      CALL SET(0.,1.,0.,1.,0.,1023.,0.,1023.,1)
+      CALL WTSTR (TKX,TKY,IDUMMY,0,0,-1)
+      CALL SET(VWPR2(1),VWPR2(2),VWPR2(3),VWPR2(4),
+     -         WND2(1),WND2(2),WND2(3),WND2(4),1)
+C
+C ADJUST 38 TO PLOTTER.
+C
+      KY = KY+52
+C
+C ADJUST 52 TO PLOTTER.
+C
+      GO TO 109
+C
+C STEP THROUGH PLOTTER GRID OF INTENSITY CELLS.
+C
+  112 IMIN = (MX1-1)/NCRTG+1
+      IMAX = (MX2-1)/NCRTG
+      JMIN = (MY1-1)/NCRTG+1
+      JMAX = (MY2-1)/NCRTG
+      XL = IMAX-IMIN+1
+      YL = JMAX-JMIN+1
+      XN = NX
+      YN = NY
+      LSRT = NLEVL/2
+      DO 130 J=JMIN,JMAX
+C
+C FIND Y FOR THIS J AND Z FOR THIS Y.
+C
+       YJ = (FLOAT(J-JMIN)+.5)/YL*(YN-1.)+1.
+         LOWY = YJ
+         YPART = YJ-FLOAT(LOWY)
+         IF (LOWY .NE. NY) GO TO 113
+         LOWY = LOWY-1
+         YPART = 1.
+  113    IPEN = 0
+         ZLFT = Z(1,LOWY)+YPART*(Z(1,LOWY+1)-Z(1,LOWY))
+         ZRHT = Z(2,LOWY)+YPART*(Z(2,LOWY+1)-Z(2,LOWY))
+         IF (NSPV .EQ. 0) GO TO 114
+         IF (Z(1,LOWY).EQ.SP .OR. Z(2,LOWY).EQ.SP .OR.
+     1       Z(1,LOWY+1).EQ.SP .OR. Z(2,LOWY+1).EQ.SP) IPEN = 1
+  114    IF (IPEN .EQ. 1) GO TO 117
+C
+C FIND INT FOR THIS Z.
+C
+         IF (ZLFT .GT. ZLEV(LSRT+1)) GO TO 116
+  115    IF (ZLFT .GE. ZLEV(LSRT)) GO TO 117
+C
+C LOOK LOWER
+C
+         IF (LSRT .LE. 1) GO TO 117
+         LSRT = LSRT-1
+         GO TO 115
+C
+C LOOK HIGHER
+C
+  116    IF (LSRT .GE. NLEVL) GO TO 117
+         LSRT = LSRT+1
+         IF (ZLFT .GT. ZLEV(LSRT+1)) GO TO 116
+C
+C OK
+C
+  117    IRHT = 2
+         LAST = LSRT
+         DO 129 I=IMIN,IMAX
+C
+C FIND X FOR THIS I AND Z FOR THIS X AND Y.
+C
+            IADD = 1
+            XI = (FLOAT(I-IMIN)+.5)/XL*(XN-1.)+1.
+            LOWX = XI
+            XPART = XI-FLOAT(LOWX)
+            IF (LOWX .NE. NX) GO TO 118
+            LOWX = LOWX-1
+            XPART = 1.
+C
+C TEST FOR INTERPOLATION POSITIONING
+C
+  118       IF (LOWX .LT. IRHT) GO TO 119
+C
+C MOVE INTERPOLATION ONE CELL TO THE RIGHT
+C
+            ZLFT = ZRHT
+            IRHT = IRHT+1
+            ZRHT = Z(IRHT,LOWY)+YPART*(Z(IRHT,LOWY+1)-Z(IRHT,LOWY))
+            IF (NSPV .EQ. 0) GO TO 118
+            IPEN = 0
+            IF (Z(IRHT-1,LOWY).EQ.SP .OR. Z(IRHT,LOWY).EQ.SP .OR.
+     1          Z(IRHT-1,LOWY+1).EQ.SP .OR. Z(IRHT,LOWY+1).EQ.SP)
+     2          IPEN = 1
+            GO TO 118
+  119       IF (IPEN .NE. 1) GO TO 123
+C
+C SPECIAL VALUE AREA
+C
+            GO TO (129,120,121,122),NSPV
+  120       MSPV = 1
+            GO TO 129
+  121       PX(1) = I*NCRTG
+            PY(1) = J*NCRTG
+            PX(2) = PX(1)+NCRTG-1
+            PY(2) = PY(1)+NCRTG-1
+            CALL GPL (2,PX,PY)
+            PYTMP = PY(1)
+            PY(1) = PY(2)
+            PY(2) = PYTMP
+            CALL GPL (2,PX,PY)
+C
+            GO TO 129
+  122       INTEN = MXLEV
+            GO TO 128
+  123       ZZ = ZLFT+XPART*(ZRHT-ZLFT)
+C
+C TEST FOR SAME INT AS LAST TIME.
+C
+            IF (ZZ .GT. ZLEV(LAST+1)) GO TO 126
+  124       IF (ZZ .GE. ZLEV(LAST)) GO TO 127
+C
+C LOOK LOWER
+C
+            IF (LAST .LE. 1) GO TO 125
+            LAST = LAST-1
+            GO TO 124
+  125       IF (ZZ .LT. ZLEV(LAST)) IADD = 0
+            GO TO 127
+C
+C LOOK HIGHER
+C
+  126       IF (LAST .GE. NLEVL) GO TO 127
+            LAST = LAST+1
+            IF (ZZ .GE. ZLEV(LAST+1)) GO TO 126
+C
+C OK
+C
+  127       ISUB = LAST+IOFFST+IADD
+            INTEN = IL(ISUB)
+            IF (NOPTN .LT. 0) INTEN = MXLEV-INTEN
+  128       CALL GRAY
+  129    CONTINUE
+  130 CONTINUE
+C
+C PUT OUT ANY REMAINING BUFFERED POINTS.
+C
+      IF (NPOINT.GT.0) THEN
+         CALL GQNT(1,IERR,WND2,VWPR2)
+         CALL SET(0.,1.,0.,1.,0.,1023.,0.,1023.,1)
+         CALL POINTS(XPNT,YPNT,NPOINT,0,0)
+         CALL SET(VWPR2(1),VWPR2(2),VWPR2(3),VWPR2(4),
+     -            WND2(1),WND2(2),WND2(3),WND2(4),1)
+      ENDIF
+C
+C CALL BOUND IF ISPV=2 AND SPECIAL VALUES WERE FOUND.
+C
+      IF (MSPV .EQ. 1) THEN
+        CALL SET(X1,X2,Y1,Y2,1.0,X3,1.0,Y3,1)
+        CALL BOUND (Z,MX,NX,NY,SP)
+      END IF
+  132 CONTINUE
+C
+C RESTORE NORMALIZATION TRANS 1 AND ORIGINAL NORMALIZATION NUMBER
+C
+      CALL SET(VWPRT(1),VWPRT(2),VWPRT(3),VWPRT(4),
+     -         WNDW(1),WNDW(2),WNDW(3),WNDW(4),IOLLS)
+      CALL SETUSV('LS',IOLLS)
+      CALL GSELNT (NTORIG)
+      RETURN
+C
+      END
+      SUBROUTINE ZLSET (Z,MX,NX,NY,ZL,NLEVL)
+      SAVE
+C
+      DIMENSION       Z(MX,NY)   ,ZL(NLEVL)
+C
+      COMMON /HAFTO2/ GLO        ,HA         ,NOPTN      ,ALPHA      ,
+     1                NSPV       ,SP         ,ICNST
+C
+      BIG = R1MACH(2)
+C
+C ZLSET PUTS THE INTENSITY LEVEL BREAK POINTS IN ZL.
+C ALL ARGUMENTS ARE AS IN HAFTON.
+C
+      LX = NX
+      LY = NY
+      NLEV = NLEVL
+      NOPT = IABS(NOPTN)
+      RALPH = 1./ALPHA
+      ICNST = 0
+      IF (GLO.NE.0. .OR. HA.NE.0.) GO TO 106
+C
+C FIND RANGE IF NOT KNOWN.
+C
+      GLO = BIG
+      HA = -GLO
+      IF (NSPV .NE. 0) GO TO 103
+      DO 102 J=1,LY
+         DO 101 I=1,LX
+            ZZ = Z(I,J)
+            GLO = AMIN1(ZZ,GLO)
+            HA = AMAX1(ZZ,HA)
+  101    CONTINUE
+  102 CONTINUE
+      GO TO 106
+  103 DO 105 J=1,LY
+         DO 104 I=1,LX
+            ZZ = Z(I,J)
+            IF (ZZ .EQ. SP) GO TO 104
+            GLO = AMIN1(ZZ,GLO)
+            HA = AMAX1(ZZ,HA)
+  104    CONTINUE
+  105 CONTINUE
+C
+C FILL ZL
+C
+  106 DELZ = HA-GLO
+      IF (DELZ .EQ. 0.) GO TO 115
+      DZ = DELZ/FLOAT(NLEV)
+      NLEVM1 = NLEV-1
+      DO 114 K=1,NLEVM1
+         ZNORM = FLOAT(K)/FLOAT(NLEV)
+         GO TO (107,108,109,110,111),NOPT
+C
+C NOPT=1
+C
+  107    ZL(K) = GLO+FLOAT(K)*DZ
+         GO TO 114
+C
+C NOPT=2
+C
+  108    ONORM = (1.-(1.-ZNORM)**ALPHA)**RALPH
+         GO TO 113
+C
+C NOPT=3
+C
+  109    ONORM = 1.-(1.-ZNORM**ALPHA)**RALPH
+         GO TO 113
+C
+C NOPT=4
+C
+  110    ONORM = .5*(1.-(ABS(ZNORM+ZNORM-1.))**ALPHA)**RALPH
+         GO TO 112
+C
+C NOPT=5
+C
+  111    ZNORM2 = ZNORM+ZNORM
+         IF (ZNORM .GT. .5) ZNORM2 = 2.-ZNORM2
+         ONORM = .5*(1.-(1.-ABS(ZNORM2)**ALPHA)**RALPH)
+  112    IF (ZNORM .GT. .5) ONORM = 1.-ONORM
+  113    ZL(K) = GLO+DELZ*ONORM
+  114 CONTINUE
+      ZL(NLEV) = BIG
+      RETURN
+  115 ICNST = 1
+      RETURN
+      END
+      SUBROUTINE GRAY
+C
+C SUBROUTINE GRAY COLORS HALF-TONE CELL (I,J) WITH INTENSITY INTEN.
+C THE ROUTINE ASSUMES 8X8 CELL SIZE ON A VIRTUAL SCREEN 1024X1024.
+C
+      DIMENSION       IFOT(16)   ,JFOT(16)
+      DIMENSION       WNDW(4)    ,VWPRT(4)
+CCC      DIMENSION       MX(16)  ,MY(16)
+      COMMON /HAFTO1/ I          ,J          ,INTEN
+      COMMON /HAFTO4/ NPTMAX     ,NPOINT     ,XPNT(50)  ,YPNT(50)
+      SAVE
+C
+      DATA
+     1 IFOT(1),IFOT(2),IFOT(3),IFOT(4),IFOT(5),IFOT(6),IFOT(7),IFOT(8)/
+     2   1,      5,      1,      5,      3,      7,      3,      7   /
+      DATA
+     1 IFOT(9),IFOT(10),IFOT(11),IFOT(12),IFOT(13),IFOT(14),IFOT(15)/
+     2   3,      7,       3,       7,       1,       5,       1/,
+     3 IFOT(16)/
+     4   5    /
+C
+      DATA
+     1 JFOT(1),JFOT(2),JFOT(3),JFOT(4),JFOT(5),JFOT(6),JFOT(7),JFOT(8)/
+     2   1,      5,      5,      1,      3,     7,      7,      3   /
+      DATA
+     1 JFOT(9),JFOT(10),JFOT(11),JFOT(12),JFOT(13),JFOT(14),JFOT(15)/
+     2   1,      5,       5,       1,       3,       7,       7/,
+     3 JFOT(16)/
+     4   3    /
+C
+      IF (INTEN) 103,103,101
+  101 I1 = I*8
+      J1 = J*8
+      IF ((NPOINT+INTEN) .LE.NPTMAX)  GO TO 1015
+      CALL GQNT(1,IERR,WNDW,VWPRT)
+      CALL SET(0.,1.,0.,1.,0.,1023.,0.,1023.,1)
+      CALL POINTS(XPNT,YPNT,NPOINT,0,0)
+      CALL SET(VWPRT(1),VWPRT(2),VWPRT(3),VWPRT(4),
+     -         WNDW(1),WNDW(2),WNDW(3),WNDW(4),1)
+      NPOINT = 0
+ 1015 DO 102 I2=1,INTEN
+      NPOINT = NPOINT + 1
+      XPNT(NPOINT) = I1+IFOT(I2)
+      YPNT(NPOINT) = J1+JFOT(I2)
+  102 CONTINUE
+  103 RETURN
+      END
+      SUBROUTINE BOUND (Z,MX,NNX,NNY,SSP)
+      DIMENSION       Z(MX,NNY)   ,PX(2)   ,PY(2)
+C
+C BOUND DRAWS A POLYGONAL BOUNDRY AROUND ANY SPECIAL-VALUE AREAS IN Z.
+C
+      SAVE
+      NX = NNX
+      NY = NNY
+C
+C VERTICAL LINES
+C
+      SP = SSP
+      DO 103 IP1=3,NX
+         I = IP1-1
+         PX(1) = I
+         PX(2) = I
+         IM1 = I-1
+         DO 102 JP1=2,NY
+            PY(2) = JP1
+            J = JP1-1
+            PY(1) = J
+            KLEFT = 0
+            IF (Z(IM1,J).EQ.SP .OR. Z(IM1,JP1).EQ.SP) KLEFT = 1
+            KCENT = 0
+            IF (Z(I,J).EQ.SP .OR. Z(I,JP1).EQ.SP) KCENT = 1
+            KRIGT = 0
+            IF (Z(IP1,J).EQ.SP .OR. Z(IP1,JP1).EQ.SP) KRIGT = 1
+            JUMP = KLEFT*4+KCENT*2+KRIGT+1
+            GO TO (102,101,102,102,101,102,102,102,102),JUMP
+  101       CALL GPL (2,PX,PY)
+  102    CONTINUE
+  103 CONTINUE
+C
+C HORIZONTAL
+C
+      DO 106 JP1=3,NY
+         J = JP1-1
+         PY(1) = J
+         PY(2) = J
+         JM1 = J-1
+         DO 105 IP1=2,NX
+            PX(2) = IP1
+            I = IP1-1
+            PX(1) = I
+            KLOWR = 0
+            IF (Z(I,JM1).EQ.SP .OR. Z(IP1,JM1).EQ.SP) KLOWR = 1
+            KCENT = 0
+            IF (Z(I,J).EQ.SP .OR. Z(IP1,J).EQ.SP) KCENT = 1
+            KUPER = 0
+            IF (Z(I,JP1).EQ.SP .OR. Z(IP1,JP1).EQ.SP) KUPER = 1
+            JUMP = KLOWR*4+KCENT*2+KUPER+1
+            GO TO (105,104,105,105,104,105,105,105,105),JUMP
+  104       CALL GPL (2,PX,PY)
+  105    CONTINUE
+  106 CONTINUE
+      RETURN
+      END
+      SUBROUTINE EZHFTN (Z,M,N)
+C
+      DIMENSION       Z(M,N)
+      SAVE
+C
+C HALF-TONE PICTURE VIA SHORTEST ARGUMENT LIST.
+C ASSUMPTIONS--
+C     ALL OF THE ARRAY IS TO BE DRAWN,
+C     LOWEST VALUE IN Z WILL BE AT LOWEST INTENSITY ON READER/PRINTER
+C     OUTPUT, HIGHEST VALUE IN Z WILL BE AT HIGHEST INTENSITY, VALUES IN
+C     BETWEEN WILL APPEAR LINEARLY SPACED, MAXIMUM POSSIBLE NUMBER OF
+C     INTENSITIES ARE USED, THE PICTURE WILL HAVE A PERIMETER DRAWN,
+C     FRAME WILL BE CALLED AFTER THE PICTURE IS DRAWN, Z IS FILLED WITH
+C     NUMBERS THAT SHOULD BE USED (NO UNKNOWN VALUES).
+C IF THESE CONDITIONS ARE NOT MET, USE HAFTON.
+C EZHFTN ARGUMENTS--
+C     Z   2 DIMENSIONAL ARRAY TO BE USED TO GENERATE A HALF-TONE PLOT.
+C     M   FIRST DIMENSION OF Z.
+C     N   SECOND DIMENSION OF Z.
+C
+      DATA FLO,HI,NLEV,NOPT,NPRM,ISPV,SPV/0.0,0.0,0,0,0,0,0.0/
+C
+C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
+C
+      CALL Q8QST4 ('GRAPHX','HAFTON','EZHFTN','VERSION  1')
+C
+      CALL HAFTON (Z,M,M,N,FLO,HI,NLEV,NOPT,NPRM,ISPV,SPV)
+C
+C +NOAO - EZHFTN no longer calls frame.
+C     CALL FRAME
+C -NOAO
+      RETURN
+      END
+C
+C-----------------------------------------------------------------------
+C
+C REVISION HISTORY---
+C
+C JULY     1984    CONVERTED TO FORTAN 77 AND GKS
+C
+C MARCH    1983    INSTITUTED BUFFERING OF POINTS WITHIN ROUTINE GRAY,
+C                  WHICH DRAMATICALLY REDUCES SIZE OF OUTPUT PLOT CODE,
+C                  METACODE.  THIS IN TURN GENERALLY IMPROVES THROUGHPUT
+C                  OF METACODE INTERPRETERS.
+C
+C FEBRUARY 1979    MODIFIED CODE TO CONFORM TO FORTRAN 66 STANDARD
+C
+C JANUARY  1978    DELETED REFERENCES TO THE  *COSY  CARDS AND
+C                  ADDED REVISION HISTORY
+C
+C-----------------------------------------------------------------------
+C