text:       SUBROUTINE ISOSRF (T,LU,MU,LV,MV,MW,EYE,MUVWP2,SLAB,TISO,IFLAG)
  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
  C
  C DIMENSION OF           T(LU,LV,MW),EYE(3),SLAB(MUVWP2,MUVWP2)
  C ARGUMENTS
  C
  C LATEST REVISION        DECEMBER 1984
  C
  C PURPOSE                ISOSRF DRAWS AN APPROXIMATION OF AN ISO-VALUED
  C                        SURFACE FROM A THREE-DIMENSIONAL ARRAY WITH
  C                        HIDDEN LINES REMOVED.
  C
  C USAGE                  IF THE FOLLOWING ASSUMPTIONS ARE MET, USE
  C
  C                            CALL EZISOS (T,MU,MV,MW,EYE,SLAB,TISO)
  C
  C                               ASSUMPTIONS:
  C                                 -- ALL OF THE T ARRAY IS TO BE USED.
  C                                 -- IFLAG IS CHOSEN INTERNALLY.
  C                                 -- FRAME IS CALLED BY EZISOS.
  C
  C                        IF THE ASSUMPTIONS ARE NOT MET, USE
  C
  C                            CALL ISOSRF (T,LU,MU,LV,MV,MW,EYE,MUVWP2,
  C                                         SLAB,TISO,IFLAG)
  C
  C ARGUMENTS
  C
  C ON INPUT               T
  C                          THREE DIMENSIONAL ARRAY OF DATA THAT DEFINES
  C                          THE ISO-VALUED SURFACE.
  C
  C                        LU
  C                          FIRST DIMENSION OF T IN THE CALLING PROGRAM.
  C
  C                        MU
  C                          THE NUMBER OF DATA VALUES OF T TO BE
  C                          PROCESSED IN THE U DIRECTION (THE FIRST
  C                          SUBSCRIPT DIRECTION).  WHEN PROCESSING THE
  C                          ENTIRE ARRAY, LU = MU (AND LV = MV).
  C
  C                        LV
  C                          SECOND DIMENSION OF T IN THE CALLING PROGRAM.
  C
  C                        MV
  C                          THE NUMBER OF DATA VALUES OF T TO BE
  C                          PROCESSED IN THE V DIRECTION (THE SECOND
  C                          SUBSCRIPT DIRECTION).
  C
  C                        MV
  C                          THE NUMBER OF DATA VALUES OF T TO BE
  C                          PROCESSED IN THE W DIRECTION (THE THIRD
  C                          SUBSCRIPT DIRECTION).
  C
  C                        EYE
  C                          THE POSITION OF THE EYE IN THREE-SPACE.  T IS
  C                          CONSIDERED TO BE IN A BOX WITH OPPOSITE
  C                          CORNERS (1,1,1) AND (MU,MV,MW).  THE EYE IS
  C                          AT (EYE(1),EYE(2),EYE(3)), WHICH MUST BE
  C                          OUTSIDE THE BOX THAT CONTAINS T.  WHILE GAINING
  C                          EXPERIENCE WITH THE ROUTINE, A GOOD CHOICE
  C                          FOR EYE MIGHT BE (5.0*MU,3.5*MV,2.0*MW).
  C
  C                        MUVWP2
  C                          THE MAXIMUM OF (MU,MV,MW)+2; THAT IS,
  C                          MUVWP2 = MAX(MU,MV,MW)+2).
  C
  C                        SLAB
  C                          A WORK SPACE USED FOR INTERNAL STORAGE.  SLAB
  C                          MUST BE AT LEAST MUVWP2*MUVWP2 WORDS LONG.
  C
  C                        TISO
  C                          THE ISO-VALUE USED TO DEFINE THE SURFACE.  THE
  C                          SURFACE DRAWN WILL SEPARATE VOLUMES OF T THAT
  C                          HAVE VALUES GREATER THAN OR EQUAL TO TISO FROM
  C                          VOLUMES OF T THAT HAVE VALUES LESS THAN TISO.
  C
  C                        IFLAG
  C                          THIS FLAG SERVES TWO PURPOSES.
  C                          .  FIRST, THE ABSOLUTE VALUE OF IFLAG
  C                             DETERMINES WHICH TYPES OF LINES ARE DRAWN
  C                             TO APPROXIMATE THE SURFACE.  THREE TYPES
  C                             OF LINES ARE CONSIDERED:  LINES OF
  C                             CONSTANT U, LINES OF CONSTANT V AND LINES
  C                             OF CONSTANT W.  THE FOLLOWING TABLE LISTS
  C                             THE TYPES OF LINES DRAWN.
  C
  C                                                LINES OF CONSTANT
  C                                                -----------------
  C                             IABS(IFLAG)           U    V    W
  C                             -----------          ---  ---  ---
  C                                  1               NO   NO   YES
  C                                  2               NO   YES  NO
  C                                  3               NO   YES  YES
  C                                  4               YES  NO   NO
  C                                  5               YES  NO   YES
  C                                  6               YES  YES  NO
  C                             0, 7 OR MORE         YES  YES  YES
  C
  C                          .  SECOND, THE SIGN OF IFLAG DETERMINES WHAT
  C                             IS INSIDE AND WHAT IS OUTSIDE, HENCE,
  C                             WHICH LINES ARE VISIBLE AND WHAT IS DONE
  C                             AT THE BOUNDARY OF T.  FOR IFLAG:
  C
  C                             POSITIVE   T VALUES GREATER THAN TISO ARE
  C                                        ASSUMED TO BE INSIDE THE SOLID
  C                                        FORMED BY THE DRAWN SURFACE.
  C                             NEGATIVE   T VALUES LESS THAN TISO ARE
  C                                        ASSUMED TO BE INSIDE THE SOLID
  C                                        FORMED BY THE DRAWN SURFACE.
  C                             IF THE ALGORITHM DRAWS A CUBE, REVERSE THE
  C                             SIGN OF IFLAG.
  C
  C ON OUTPUT              T,LU,MU,LV,MV,MW,EYE,MUVWP2,TISO AND IFLAG ARE
  C                        UNCHANGED.  SLAB HAS BEEN WRITTEN IN.
  C
  C NOTE                   .   THIS ROUTINE IS FOR LOWER RESOLUTION ARRAYS
  C                            THAN ISOSRFHR.  40 BY 40 BY 40 IS A
  C                            PRACTICAL MAXIMUM.
  C                        .   TRANSFORMATIONS CAN BE ACHIEVED BY
  C                            ADJUSTING SCALING STATEMENT FUNCTIONS IN
  C                            ISOSRF, SET3D AND TR32.
  C                        .   THE HIDDEN-LINE ALGORITHM IS NOT EXACT, SO
  C                            VISIBILITY ERRORS CAN OCCUR.
  C                        .   THREE-DIMENSIONAL PERSPECTIVE CHARACTER
  C                            LABELING OF ISOSRF IS POSSIBLE BY USING
  C                            THE UTILITY PWRZI.  FOR A DESCRIPTION OF
  C                            THE USAGE, SEE THE PWRZI DOCUMENTATION.
  C
  C ENTRY POINTS           ISOSRF, EZISOS, SET3D, TRN32I, ZEROSC,
  C                        STCNTR, DRCNTR, TR32, FRSTS, KURV1S, KURV2S,
  C                        FRSTC, FILLIN, DRAWI, ISOSRB, MMASK
  C
  C COMMON BLOCKS          ISOSR1, ISOSR2, ISOSR3, ISOSR4, ISOSR5,
  C                        ISOSR6, ISOSR7, ISOSR8, ISOSR9, TEMPR,
  C                        PWRZ1I
  C
  C REQUIRED LIBRARY       THE ERPRT77 PACKAGE AND THE SPPS.
  C ROUTINES
  C
  C I/O                    PLOTS SURFACE
  C
  C PRECISION              SINGLE
  C
  C LANGUAGE               FORTRAN 77
  C
  C HISTORY                DEVELOPED FOR USERS OF ISOSRFHR WITH SMALLER
  C                        ARRAYS.
  C
  C ALGORITHM              CUTS THROUGH THE THREE-DIMENSIONAL ARRAY ARE
  C                        CONTOURED WITH A SMOOTHING CONTOURER WHICH ALSO
  C                        MARKS A MODEL OF THE PLOTTING PLANE.  INTERIORS
  C                        OF BOUNDARIES ARE FILLED IN AND THE RESULT IS
  C                        .OR.ED INTO ANOTHER MODEL OF THE PLOTTING PLANE
  C                        WHICH IS USED TO TEST SUBSEQUENT CONTOUR LINES
  C                        FOR VISIBILITY.
  C
  C TIMING                 VARIES WIDELY WITH SIZE OF T AND THE VOLUME OF
  C                        THE SPACE ENCLOSED BY THE SURFACE DRAWN.
  C
  C  **NOTE**                  SPACE REQUIREMENTS CAN BE REDUCED BY
  C                            CHANGING THE SIZE OF THE ARRAYS ISCR, ISCA
  C                            (FOUND IN COMMON ISOSR2), MASK(FOUND IN
  C                            COMMON ISOSR5) AND THE VARIABLE NBPW
  C                            (COMMON ISOSR5).
  C                            ISCR AND ISCA NEED 128X128 BITS.  SO ON A
  C                            64 BIT MACHINE ISCR, ISCA CAN BE
  C                            DIMENSIONED TO (2,128). NBPW SET IN
  C                            SUBROUTINE MMASK SHOULD CONTAIN THE
  C                            NUMBER OF BITS PER WORD YOU WISH TO
  C                            UTILIZE.
  C                            THE DIMENSION OF MASK AND NMASK SHOULD
  C                            EQUAL THE VALUE OF NBPW.
  C                            LS SHOULD BE SET TO THE FIRST DIMENSION
  C                            OF ISCA AND ISCR.
  C
  C               EXAMPLES:
  C                            ON A 60 BIT MACHINE:
  C                                 DIMENSION ISCA(4,128), ISCR(4,128)
  C                                 DIMENSION MASK(32)
  C                                 NBPW = 32
  C                            ON A 64 BIT MACHINE:
  C                                 DIMENSION ISCA(2,128), ISCR(2,128)
  C                                 DIMENSION MASK(64)
  C                                 NBPW = 64
  C
  C INTERNAL PARAMETERS    NAME  DEFAULT           FUNCTION
  C                        ----  -------           --------
  C                        IREF     1     FLAG TO CONTROL DRAWING OF AXES.
  C                                       .IREF=NONZERO   DRAW AXES.
  C                                       .IREF=ZERO      DO NOT DRAW AXES.
  C
  C
        SAVE
        DIMENSION       T(LU,LV,MW),EYE(3)     ,SLAB(MUVWP2,MUVWP2)
  C
        COMMON /ISOSR1/ ISLBT      ,U          ,V          ,W
        COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
       1                ISCA(8,128)
        COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
       1                YMAX       ,BIGD       ,R0
        COMMON /ISOSR4/ RX         ,RY
        COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
        COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
       1                IS         ,ISS        ,NP         ,CV         ,
       2                INX(8)     ,INY(8)     ,IR(500)    ,NR
        COMMON /ISOSR7/ IENTRY     ,IONES
        COMMON /ISOSR9/ BIG        ,IXBIT
  C
        LOGICAL         GENDON
        DATA IREF/1/
  C
        AVE(A,B) = (A+B)*.5
  C
  C A.S.F. FOR SCALING
  C
        SU(UTEMP) = UTEMP
        SV(VTEMP) = VTEMP
        SW(WTEMP) = WTEMP
  C
  C +NOAO - Blockdata ISOSRB rewritten as run time initialization
  C     EXTERNAL ISOSRB
        call isosrb
  C -NOAO
  C
  C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
  C
        CALL Q8QST4 ('NSSL','ISOSRF','ISOSRF','VERSION 12')
        NERR = 0
  C
  C 3-SPACE  U,V,W,IU,IV,IW,ETC
  C 2-SPACE  X,Y,IX,IY,ETC
  C
  C  INITIALIZE MASKS
  C
        IF (.NOT.GENDON) CALL MMASK
  C
  C  SET SHIFT VALUE FOR X,Y PACKING
  C
  C  IF YOUR MACHINE HAS MORE THAN 16 BITS PER WORD THIS CHECK MAY BE
  C  MODIFIED
  C
        IF (LU .LE. 256) GO TO  10
        NERR = NERR + 1
        CALL SETER('DIMENSION OF CUBE EXCEEDS 256',NERR,2)
        RETURN
     10 DO  20 J=1,30
           IF (LU .LE. 2**(J-1)) GO TO  30
     20 CONTINUE
     30 IXBIT = J
        NU = MU
        NUP2 = NU+2
        NV = MV
        NVP2 = NV+2
        NW = MW
        NWP2 = NW+2
        FNU = NU
        FNV = NV
        FNW = NW
        SU1 = SU(1.)
        SV1 = SV(1.)
        SW1 = SW(1.)
        SUNU = SU(FNU)
        SVNV = SV(FNV)
        SWNW = SW(FNW)
        AVEU = AVE(SU1,SUNU)
        AVEV = AVE(SV1,SVNV)
        AVEW = AVE(SW1,SWNW)
        EYEU = EYE(1)
        EYEV = EYE(2)
        EYEW = EYE(3)
        NUVWP2 = MUVWP2
        TVAL = TISO
        NFLAG = IABS(IFLAG)
        IF (NFLAG.EQ.0 .OR. NFLAG.GE.8) NFLAG = 7
  C
  C SET UP SCALING
  C
        FACT = -ISIGN(1,IFLAG)
        CALL SET3D (EYE,1.,FNU,1.,FNV,1.,FNW)
  C
  C BOUND LOWER AND LEFT EDGE OF SLAB
  C
        EDGE = SIGN(BIG,FACT)
        DO  40 IUVW=1,NUVWP2
           SLAB(IUVW,1) = EDGE
           SLAB(1,IUVW) = EDGE
     40 CONTINUE
  C
  C SLICES PERPENDICULAR TO U. THAT IS, V W SLICES. T OF CONSTANT U.
  C
        IF (NFLAG .LT. 4) GO TO 100
        CALL ZEROSC
        ISLBT = -1
  C
  C BOUND UPPER AND RIGHT EDGE OF SLAB.
  C
        DO  50 IV=2,NVP2
           SLAB(IV,NWP2) = EDGE
     50 CONTINUE
        DO  60 IW=2,NWP2
           SLAB(NVP2,IW) = EDGE
     60 CONTINUE
  C
  C GO THRU 3-D ARRAY IN U DIRECTION.  IUEW=IU EITHER WAY.
  C PICK IU BASED ON EYEU.
  C
        DO  90 IUEW=1,NU
           IU = IUEW
           IF (EYEU .GT. AVEU) IU = NU+1-IUEW
           U = IU
  C
  C LOAD THIS SLICE OF T INTO SLAB.
  C
           DO  80 IV=1,NV
              DO  70 IW=1,NW
                 SLAB(IV+1,IW+1) = T(IU,IV,IW)
     70       CONTINUE
     80    CONTINUE
  C
  C CONTOUR THIS SLAB.
  C
           CALL STCNTR (SLAB,NUVWP2,NVP2,NWP2,TVAL)
  C
  C CONSTRUCT VISIBILITY ARRAY.
  C
           CALL FILLIN
     90 CONTINUE
  C
  C SLICES PERPENDICULAR TO V. U W SLICES. T OF CONSTANT V.
  C
    100 IF (MOD(NFLAG/2,2) .EQ. 0) GO TO 160
        CALL ZEROSC
        ISLBT = 0
  C
  C BOUND UPPER AND RIGHT EDGE OF SLAB.
  C
        DO 110 IU=2,NUP2
           SLAB(IU,NWP2) = EDGE
    110 CONTINUE
        DO 120 IW=2,NWP2
           SLAB(NUP2,IW) = EDGE
    120 CONTINUE
  C
  C GO THRU T IN V DIRECTION.  IVEW=IV EITHER WAY.
  C
        DO 150 IVEW=1,NV
           IV = IVEW
           IF (EYEV .GT. AVEV) IV = NV+1-IVEW
           V = IV
  C
  C LOAD THIS SLICE OF T INTO SLAB.
  C
           DO 140 IU=1,NU
              DO 130 IW=1,NW
                 SLAB(IU+1,IW+1) = T(IU,IV,IW)
    130       CONTINUE
    140    CONTINUE
  C
  C CONTOUR THIS SLAB.
  C
           CALL STCNTR (SLAB,NUVWP2,NUP2,NWP2,TVAL)
  C
  C CONSTRUCT VISIBILITY ARRAY.
  C
           CALL FILLIN
    150 CONTINUE
  C
  C SLICES PERPENDICULAR TO W. U V SLICES. T OF CONSTANT W.
  C
    160 IF (MOD(NFLAG,2) .EQ. 0) GO TO 220
        CALL ZEROSC
  C
        ISLBT = 1
  C
  C BOUND UPPER AND RIGHT EDGE OF SLAB.
  C
        DO 170 IU=2,NUP2
           SLAB(IU,NVP2) = EDGE
    170 CONTINUE
        DO 180 IV=2,NVP2
           SLAB(NUP2,IV) = EDGE
    180 CONTINUE
  C
  C GO THRU T IN W DIRECTION.
  C
        DO 210 IWEW=1,NW
           IW = IWEW
           IF (EYEW .GT. AVEW) IW = NW+1-IWEW
           W = IW
  C
  C LOAD THIS SLICE OF T INTO SLAB.
  C
           DO 200 IU=1,NU
              DO 190 IV=1,NV
                 SLAB(IU+1,IV+1) = T(IU,IV,IW)
    190       CONTINUE
    200    CONTINUE
  C
  C CONTOUR THIS SLAB.
  C
           CALL STCNTR (SLAB,NUVWP2,NUP2,NVP2,TVAL)
  C
  C CONSTRUCT VISIBILITY ARRAY.
  C
           CALL FILLIN
    210 CONTINUE
  C
  C DRAW REFERENCE PLANE EDGES AND W AXIS.
  C
    220 IF (IREF .EQ. 0) RETURN
        CALL TRN32I (SU1,SV1,SW1,XT,YT,DUM,2)
        IF (EYEV .LT. SV1) GO TO 240
        CALL FRSTC (IFIX(XT),IFIX(YT),1)
        DO 230 IU=2,NU
           CALL TRN32I (SU(FLOAT(IU)),SV1,SW1,XT,YT,DUM,2)
           CALL FRSTC (IFIX(XT),IFIX(YT),2)
    230 CONTINUE
        GO TO 250
    240 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
        CALL TRN32I (SUNU,SV1,SW1,XT,YT,DUM,2)
        CALL PLOTIT (IFIX(XT),IFIX(YT),1)
    250 IF (EYEU .GT. SUNU) GO TO 270
        CALL FRSTC (IFIX(XT),IFIX(YT),1)
        DO 260 IV=2,NV
           CALL TRN32I (SUNU,SV(FLOAT(IV)),SW1,XT,YT,DUM,2)
           CALL FRSTC (IFIX(XT),IFIX(YT),2)
    260 CONTINUE
        GO TO 280
    270 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
        CALL TRN32I (SUNU,SVNV,SW1,XT,YT,DUM,2)
        CALL PLOTIT (IFIX(XT),IFIX(YT),1)
    280 IF (EYEV .GT. SVNV) GO TO 300
        CALL FRSTC (IFIX(XT),IFIX(YT),1)
        DO 290 IUOW=2,NU
           CALL TRN32I (SU(FLOAT(NU-IUOW+1)),SVNV,SW1,XT,YT,DUM,2)
           CALL FRSTC (IFIX(XT),IFIX(YT),2)
    290 CONTINUE
        GO TO 310
    300 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
        CALL TRN32I (SU1,SVNV,SW1,XT,YT,DUM,2)
        CALL PLOTIT (IFIX(XT),IFIX(YT),1)
    310 IF (EYEU .LT. SU1) GO TO 330
        CALL FRSTC (IFIX(XT),IFIX(YT),1)
        DO 320 IVOW=2,NV
           CALL TRN32I (SU1,SV(FLOAT(NV-IVOW+1)),SW1,XT,YT,DUM,2)
           CALL FRSTC (IFIX(XT),IFIX(YT),2)
    320 CONTINUE
        GO TO 340
    330 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
        CALL TRN32I (SU1,SV1,SW1,XT,YT,DUM,2)
        CALL PLOTIT (IFIX(XT),IFIX(YT),1)
    340 IF (EYEU.LE.SU1 .OR. EYEV.LE.SV1) GO TO 360
        CALL FRSTC (IFIX(XT),IFIX(YT),1)
        DO 350 IW=2,NW
           CALL TRN32I (SU1,SV1,SW(FLOAT(IW)),XT,YT,DUM,2)
           CALL FRSTC (IFIX(XT),IFIX(YT),2)
    350 CONTINUE
  C +NOAO - Plotit buffer needs to be flushed before returning.
        call plotit (0, 0, 2)
  C -NOAO
        RETURN
    360 CALL PLOTIT (IFIX(XT),IFIX(YT),0)
        CALL TRN32I (SU1,SV1,SWNW,XT,YT,DUM,2)
        CALL PLOTIT (IFIX(XT),IFIX(YT),1)
  C +NOAO - Plotit buffer needs to be flushed before returning.
        call plotit (0, 0, 2)
  C -NOAO
        RETURN
        END
        SUBROUTINE EZISOS (T,MU,MV,MW,EYE,SLAB,TISO)
  C
        SAVE
        DIMENSION       T(MU,MV,MW),EYE(3)
  C
          DATA ANG,PI/.35,3.141592/
  C
  C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR
  C
        CALL Q8QST4 ('NSSL','ISOSRF','EZISOS','VERSION 12')
  C
  C ARGUMENTS DESCRIBED IN ISOSRF
  C
  C PICK TYPES OF LINES TO DRAW
  C
        NU = MU
        NV = MV
        NW = MW
        TVAL = TISO
        MAX = MAX0(NU,NV,NW)+2
        ATU = NU/2
        ATV = NV/2
        ATW = NW/2
        EYEU = EYE(1)
        EYEV = EYE(2)
        EYEW = EYE(3)
        RU = EYEU-ATU
        RV = EYEV-ATV
        RW = EYEW-ATW
        RU2 = RU*RU
        RV2 = RV*RV
        RW2 = RW*RW
        DU = SQRT(RV2+RW2)
        DV = SQRT(RU2+RW2)
        DW = SQRT(RU2+RV2)
        DR = 1./SQRT(RU2+RV2+RW2)
  C
  C  COMPUTE THE ARCCOSINE
  C
        TU = DU*DR
        ANGU = ATAN(ABS(SQRT(1.-TU*TU)/TU))
        IF (TU .LE. 0.) ANGU = PI-ANGU
        TV = DV*DR
        ANGV = ATAN(ABS(SQRT(1.-TV*TV)/TV))
        IF (TV .LE. 0.) ANGV = PI-ANGV
        TW = DW*DR
        ANGW = ATAN(ABS(SQRT(1.-TW*TW)/TW))
        IF (TW .LE. 0.) ANGW = PI-ANGW
  C
  C BREAK POINT IS ABOUT 20 DEGREES OR ABOUT .35 RADIANS
  C
        IFLAG = 0
        IF (ANGU .GT. ANG) IFLAG = IFLAG+4
        IF (ANGV .GT. ANG) IFLAG = IFLAG+2
        IF (ANGW .GT. ANG) IFLAG = IFLAG+1
  C
  C FIND SIGN OF IFLAG
  C
        ICNT = 0
        IF (ABS(RU) .LE. ATU) GO TO  30
        IU = 1
        IF (EYEU .GT. ATU) IU = NU
        DO  20 IW=1,NW
           DO  10 IV=1,NV
              IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
              ICNT = ICNT+1
     10    CONTINUE
     20 CONTINUE
     30 IF (ABS(RV) .LE. ATV) GO TO  60
        IV = 1
        IF (EYEV .GT. ATV) IV = NV
        DO  50 IW=1,NW
           DO  40 IU=1,NU
              IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
              ICNT = ICNT+1
     40    CONTINUE
     50 CONTINUE
     60 IF (ABS(RW) .LE. ATW) GO TO  90
        IW = 1
        IF (EYEW .GT. ATW) IW = NW
        DO  80 IV=1,NV
           DO  70 IU=1,NU
              IF (T(IU,IV,IW) .GT. TVAL) ICNT = ICNT-2
              ICNT = ICNT+1
     70    CONTINUE
     80 CONTINUE
     90 IFLAG = ISIGN(IFLAG,ICNT)
        CALL ISOSRF (T,NU,NU,NV,NV,NW,EYE,MAX,SLAB,TVAL,IFLAG)
  C +NOAO - Call to frame is suppressed.
  C     CALL FRAME
  C -NOAO
        RETURN
        END
        SUBROUTINE SET3D (EYE,ULO,UHI,VLO,VHI,WLO,WHI)
        SAVE
        COMMON /TEMPR/  RZERO
  C
        DIMENSION       EYE(3)
  C
        COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
       1                YMAX       ,BIGD       ,R0
        COMMON /PWRZ1I/ UUMIN      ,UUMAX      ,VVMIN      ,VVMAX      ,
       1                WWMIN      ,WWMAX      ,DELCRT     ,EYEU       ,
       2                EYEV       ,EYEW
  C
  C
        AVE(A,B) = (A+B)*.5
  C
  C A.S.F. FOR SCALING
  C
        SU(UTEMP) = UTEMP
        SV(VTEMP) = VTEMP
        SW(WTEMP) = WTEMP
  C
  C CONSTANTS FOR PWRZ
  C
        UUMIN = ULO
        UUMAX = UHI
        VVMIN = VLO
        VVMAX = VHI
        WWMIN = WLO
        WWMAX = WHI
        EYEU = EYE(1)
        EYEV = EYE(2)
        EYEW = EYE(3)
  C
  C FIND CORNERS IN 2-SPACE FOR 3-SPACE BOX CONTAINING OBJECT
  C
        ISCALE = 0
        ATU = AVE(SU(UUMIN),SU(UUMAX))
        ATV = AVE(SV(VVMIN),SV(VVMAX))
        ATW = AVE(SW(WWMIN),SW(WWMAX))
        BIGD = 0.
        IF (RZERO .LE. 0.) GO TO  10
  C
  C RELETIVE SIZE FEATURE IN USE.
  C GENERATE EYE POSITION THAT MAKES BOX HAVE MAXIMUM PROJECTED SIZE.
  C
        ALPHA = -(VVMIN-ATV)/(UUMIN-ATU)
        VVEYE = -RZERO/SQRT(1.+ALPHA*ALPHA)
        UUEYE = VVEYE*ALPHA
        VVEYE = VVEYE+ATV
        UUEYE = UUEYE+ATU
        WWEYE = ATW
        CALL TRN32I (ATU,ATV,ATW,UUEYE,VVEYE,WWEYE,1)
        CALL TRN32I (UUMIN,VVMIN,ATW,XMIN,DUMM,DUMM,2)
        CALL TRN32I (UUMAX,VVMIN,WWMIN,DUMM,YMIN,DUMM,2)
        CALL TRN32I (UUMAX,VVMAX,ATW,XMAX,DUMM,DUMM,2)
        CALL TRN32I (UUMAX,VVMIN,WWMAX,DUMM,YMAX,DUMM,2)
        BIGD = SQRT((UUMAX-UUMIN)**2+(VVMAX-VVMIN)**2+(WWMAX-WWMIN)**2)*.5
        R0 = RZERO
        GO TO  20
     10 CALL TRN32I (ATU,ATV,ATW,EYE(1),EYE(2),EYE(3),1)
        CALL TRN32I (SU(UUMIN),SV(VVMIN),SW(WWMIN),X1,Y1,DUM,2)
        CALL TRN32I (SU(UUMIN),SV(VVMIN),SW(WWMAX),X2,Y2,DUM,2)
        CALL TRN32I (SU(UUMIN),SV(VVMAX),SW(WWMIN),X3,Y3,DUM,2)
        CALL TRN32I (SU(UUMIN),SV(VVMAX),SW(WWMAX),X4,Y4,DUM,2)
        CALL TRN32I (SU(UUMAX),SV(VVMIN),SW(WWMIN),X5,Y5,DUM,2)
        CALL TRN32I (SU(UUMAX),SV(VVMIN),SW(WWMAX),X6,Y6,DUM,2)
        CALL TRN32I (SU(UUMAX),SV(VVMAX),SW(WWMIN),X7,Y7,DUM,2)
        CALL TRN32I (SU(UUMAX),SV(VVMAX),SW(WWMAX),X8,Y8,DUM,2)
        XMIN = AMIN1(X1,X2,X3,X4,X5,X6,X7,X8)
        XMAX = AMAX1(X1,X2,X3,X4,X5,X6,X7,X8)
        YMIN = AMIN1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8)
        YMAX = AMAX1(Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8)
  C
  C ADD RIGHT AMOUNT TO KEEP PICTURE SQUARE
  C
     20 WIDTH = XMAX-XMIN
        HIGHT = YMAX-YMIN
        DIF = .5*(WIDTH-HIGHT)
        IF (DIF)  30, 50, 40
     30 XMIN = XMIN+DIF
        XMAX = XMAX-DIF
        GO TO  50
     40 YMIN = YMIN-DIF
        YMAX = YMAX+DIF
     50 ISCALE = 1
        CALL TRN32I (ATU,ATV,ATW,EYE(1),EYE(2),EYE(3),1)
        RETURN
        END
        SUBROUTINE TRN32I (U,V,W,XT,YT,ZT,IENT)
  C
  C THIS ROUTINE IMPLEMENTS THE 3-SPACE TO 2-SPACE TRANSFOR-
  C MATION BY KUBER, SZABO AND GIULIERI, THE PERSPECTIVE
  C REPRESENTATION OF FUNCTIONS OF TWO VARIABLES. J. ACM 15,
  C 2, 193-204,1968.
  C ARGUMENTS FOR SET
  C U,V,W    ARE THE 3-SPACE COORDINATES OF THE INTERSECTION
  C          OF THE LINE OF SIGHT AND THE IMAGE PLANE.  THIS
  C          POINT CAN BE THOUGHT OF AS THE POINT LOOKED AT.
  C XT,YT,ZT ARE THE 3-SPACE COORDINATES OF THE EYE POSITION.
  C
  C TRN32 ARGUMENTS
  C U,V,W    ARE THE 3-SPACE COORDINATES OF A POINT TO BE
  C          TRANSFORMED.
  C XT,YT    THE RESULTS OF THE 3-SPACE TO 2-SPACE TRANSFOR-
  C          MATION.  WHEN ISCALE=0, XT AND YT ANR IN THE SAME
  C          UNITS AS U,V, AND W.  WHEN ISCALE'0, XT AND YT
  C          ARE IN PLOTTER COORDINATES.
  C ZT       NOT USED.
  C
        SAVE
        COMMON /PWRZ1I/ UUMIN      ,UUMAX      ,VVMIN      ,VVMAX      ,
       1                WWMIN      ,WWMAX      ,DELCRT     ,EYEU       ,
       2                EYEV       ,EYEW
        COMMON /ISOSR3/ ISCALE     ,XMIN       ,XMAX       ,YMIN       ,
       1                YMAX       ,BIGD       ,R0
  C
  C RANGE OF PLOTTER COORDINATES
  C
  C
  C  WARNING
  C  IF PLOTTER MAXIMUM VALUE RANGES (IN X OR Y DIRECTION) FALL BELOW
  C  101, THEN CHANGES MUST BE MADE IN SUBROUTINE FRSTC. THE REQUIRED
  C  CHANGES ARE MARKED BY WARNING COMMENTS IN FRSTC.
        DATA NLX,NBY,NRX,NTY/10,10,32760,32760/
          DATA PI/3.1411592/
  C
  C STORE THE PARAMETERS OF THE SET CALL FOR USE
  C WITH THE TRANSLATE CALL
  C
  C DECIDE IF SET OR TRANSLATE CALL
  C
        IF (IENT .NE. 1) GO TO  50
        AU = U
        AV = V
        AW = W
        EU = XT
        EV = YT
        EW = ZT
  C
  C
  C
  C
  C
        DU = AU-EU
        DV = AV-EV
        DW = AW-EW
        D = SQRT(DU*DU+DV*DV+DW*DW)
        COSAL = DU/D
        COSBE = DV/D
        COSGA = DW/D
  C
  C  COMPUTE THE ARCCOSINE
  C
        AL = ATAN(ABS(SQRT(1.-COSAL*COSAL)/COSAL))
        IF (COSAL .LE. 0.) AL = PI-AL
        BE = ATAN(ABS(SQRT(1.-COSBE*COSBE)/COSBE))
        IF (COSBE .LE. 0.) BE = PI-BE
        GA = ATAN(ABS(SQRT(1.-COSGA*COSGA)/COSGA))
        IF (COSGA .LE. 0.) GA = PI-GA
        SINGA = SIN(GA)
  C
  C THE 3-SPACE POINT LOOKED AT IS TRANSFORMED INTO (0,0) OF
  C THE 2-SPACE.  THE 3-SPACE W AXIS IS TRANSFORMED INTO THE
  C 2-SPACE Y AXIS.  IF THE LINE OF SIGHT IS CLOSE TO PARALLEL
  C TO THE 3-SPACE W AXIS, THE 3-SPACE V AXIS IS CHOSEN (IN-
  C STEAD OF THE 3-SPACE W AXIS) TO BE TRANSFORMED INTO THE
  C 2-SPACE Y AXIS.
  C
        ASSIGN  90 TO JDONE
        IF (ISCALE)  10, 30, 10
     10 X0 = XMIN
        Y0 = YMIN
        X1 = NLX
        Y1 = NBY
        X2 = NRX-NLX
        Y2 = NTY-NBY
        X3 = X2/(XMAX-XMIN)
        Y3 = Y2/(YMAX-YMIN)
        X4 = NRX
        Y4 = NTY
        FACT = 1.
        IF (BIGD .LE. 0.) GO TO  20
        X0 = -BIGD
        Y0 = -BIGD
        X3 = X2/(2.*BIGD)
        Y3 = Y2/(2.*BIGD)
        FACT = R0/D
     20 DELCRT = X2
        ASSIGN  80 TO JDONE
     30 IF (SINGA .LT. 0.0001) GO TO  40
        R = 1./SINGA
        ASSIGN  70 TO JUMP
        RETURN
     40 SINBE = SIN(BE)
        R = 1./SINBE
        ASSIGN  60 TO JUMP
        RETURN
  C
  C********************  ENTRY TRN32  ************************
  C     ENTRY TRN32 (U,V,W,XT,YT,ZT)
  C
     50 UU = U
        VV = V
        WW = W
        Q = D/((UU-EU)*COSAL+(VV-EV)*COSBE+(WW-EW)*COSGA)
        GO TO JUMP,( 60, 70)
     60 UU = ((EW+Q*(WW-EW)-AW)*COSAL-(EU+Q*(UU-EU)-AU)*COSGA)*R
        VV = (EV+Q*(VV-EV)-AV)*R
        GO TO JDONE,( 80, 90)
     70 UU = ((EU+Q*(UU-EU)-AU)*COSBE-(EV+Q*(VV-EV)-AV)*COSAL)*R
        VV = (EW+Q*(WW-EW)-AW)*R
        GO TO JDONE,( 80, 90)
     80 XT = AMIN1(X4,AMAX1(X1,X1+X3*(FACT*UU-X0)))
        YT = AMIN1(Y4,AMAX1(Y1,Y1+Y3*(FACT*VV-Y0)))
        RETURN
     90 XT = UU
        YT = VV
        RETURN
        END
        SUBROUTINE ZEROSC
        SAVE
  C
        COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
       1                ISCA(8,128)
  C
  C ZERO BOTH SCRENE MODELS.
  C
        DO  20 I=1,LX
           DO  10 J=1,NY
              ISCR(I,J) = 0
              ISCA(I,J) = 0
     10    CONTINUE
     20 CONTINUE
        RETURN
        END
        SUBROUTINE STCNTR (Z,L,M,N,CONV)
  C
        SAVE
        DIMENSION       Z(L,N)
  C
  C THIS ROUTINE FINDS THE BEGINNINGS OF ALL CONTOUR LINES AT LEVEL CONV.
  C FIRST THE EDGES ARE SEARCHED FOR LINES INTERSECTING THE EDGE (OPEN
  C LINES) THEN THE INTERIOR IS SEARCHED FOR LINES WHICH DO NOT INTERSECT
  C THE EDGE (CLOSED LINES).  BEGINNINGS ARE STORED IN IR TO PREVENT RE-
  C TRACING OF LINES.  IF IR IS FILLED, THE SEARCH IS STOPPED FOR THIS
  C CONV.
  C
        COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
       1                IS         ,ISS        ,NP         ,CV         ,
       2                INX(8)     ,INY(8)     ,IR(500)    ,NR
        COMMON /ISOSR7/ IENTRY     ,IONES
        COMMON /ISOSR9/ BIG        ,IXBIT
  C
  C  PACK X AND Y
  C
        IPXY(I1,J1) = ISHIFT(I1,IXBIT)+J1
  C
        IENTRY = 0
        NP = 0
        CV = CONV
  C
  C THE FOLLOWING CODE SHOULD BE RE-ENABLED IF THIS ROUTINE IS USED FOR
  C GENERAL CONTOURING
  C
  C     ISS=0
  C     DO 2 IP1=2,M
  C     I=IP1-1
  C     IF(Z(I,1).GE.CV.OR.Z(IP1,1).LT.CV) GO TO 1
  C     IX=IP1
  C     IY=1
  C     IDX=-1
  C     IDY=0
  C     IS=1
  C     CALL DRLINE(Z,L,M,N)
  C   1 IF(Z(IP1,N).GE.CV.OR.Z(I,N).LT.CV) GO TO 2
  C     IX=I
  C     IY=N
  C     IDX=1
  C     IDY=0
  C     IS=5
  C     CALL DRLINE(Z,L,M,N)
  C   2 CONTINUE
  C     DO 4 JP1=2,N
  C     J=JP1-1
  C     IF(Z(M,J).GE.CV.OR.Z(M,JP1).LT.CV) GO TO 3
  C     IX=M
  C     IY=JP1
  C     IDX=0
  C     IDY=-1
  C     IS=7
  C     CALL DRLINE(Z,L,M,N)
  C   3 IF(Z(1,JP1).GE.CV.OR.Z(1,J).LT.CV) GO TO 4
  C     IX=1
  C     IY=J
  C     IDX=0
  C     IDY=1
  C     IS=3
  C     CALL DRLINE(Z,L,M,N)
  C   4 CONTINUE
  C
        ISS = 1
        DO  40 JP1=3,N
           J = JP1-1
           DO  30 IP1=2,M
              I = IP1-1
              IF (Z(I,J).GE.CV .OR. Z(IP1,J).LT.CV) GO TO  30
              IXY = IPXY(IP1,J)
              IF (NP .EQ. 0) GO TO  20
              DO  10 K=1,NP
                 IF (IR(K) .EQ. IXY) GO TO  30
     10       CONTINUE
     20       NP = NP+1
              IF (NP .GT. NR) RETURN
              IR(NP) = IXY
              IX = IP1
              IY = J
              IDX = -1
              IDY = 0
              IS = 1
              CALL DRCNTR (Z,L,M,N)
     30    CONTINUE
     40 CONTINUE
        RETURN
        END
        SUBROUTINE DRCNTR (Z,L,MM,NN)
        SAVE
  C
        DIMENSION       Z(L,NN)
  C
  C THIS ROUTINE TRACES A CONTOUR LINE WHEN GIVEN THE BEGINNING BY STLINE.
  C TRANSFORMATIONS CAN BE ADDED BY DELETING THE STATEMENT FUNCTIONS FOR
  C FX AND FY IN DRLINE AND MINMAX AND ADDING EXTERNAL FUNCTIONS.
  C X=1. AT Z(1,J), X=FLOAT(M) AT Z(M,J). X TAKES ON NON-INTEGER VALUES.
  C Y=1. AT Z(I,1), Y=FLOAT(N) AT Z(I,N). Y TAKES ON NON-INTEGER VALUES.
  C
        COMMON /ISOSR6/ IX         ,IY         ,IDX        ,IDY        ,
       1                IS         ,ISS        ,NP         ,CV         ,
       2                INX(8)     ,INY(8)     ,IR(500)    ,NR
        COMMON /ISOSR9/ BIG        ,IXBIT
  C
        LOGICAL         IPEN       ,IPENO
  C
        DATA IOFFP,SPVAL/0,0./
        DATA IPEN,IPENO/.TRUE.,.TRUE./
  C
  C  PACK X AND Y
  C
        IPXY(I1,J1) = ISHIFT(I1,IXBIT)+J1
        FX(X1,Y1) = X1
        FY(X1,Y1) = Y1
        C(P11,P21) = (P11-CV)/(P11-P21)
  C
        M = MM
        N = NN
        IF (IOFFP .EQ. 0) GO TO  10
        ASSIGN 100 TO JUMP1
        ASSIGN 150 TO JUMP2
        GO TO  20
     10 ASSIGN 120 TO JUMP1
        ASSIGN 160 TO JUMP2
     20 IX0 = IX
        IY0 = IY
        IS0 = IS
        IF (IOFFP .EQ. 0) GO TO  30
        IX2 = IX+INX(IS)
        IY2 = IY+INY(IS)
        IPEN = Z(IX,IY).NE.SPVAL .AND. Z(IX2,IY2).NE.SPVAL
        IPENO = IPEN
     30 IF (IDX .EQ. 0) GO TO  40
        Y = IY
        ISUB = IX+IDX
        X = C(Z(IX,IY),Z(ISUB,IY))*FLOAT(IDX)+FLOAT(IX)
        GO TO  50
     40 X = IX
        ISUB = IY+IDY
        Y = C(Z(IX,IY),Z(IX,ISUB))*FLOAT(IDY)+FLOAT(IY)
     50 IF (IPEN) CALL FRSTS (FX(X,Y),FY(X,Y),1)
     60 IS = IS+1
        IF (IS .GT. 8) IS = IS-8
        IDX = INX(IS)
        IDY = INY(IS)
        IX2 = IX+IDX
        IY2 = IY+IDY
        IF (ISS .NE. 0) GO TO  70
        IF (IX2.GT.M .OR. IY2.GT.N .OR. IX2.LT.1 .OR. IY2.LT.1) GO TO 190
     70 IF (CV-Z(IX2,IY2))  80, 80, 90
     80 IS = IS+4
        IX = IX2
        IY = IY2
        GO TO  60
     90 IF (IS/2*2 .EQ. IS) GO TO  60
        GO TO JUMP1,(100,120)
    100 ISBIG = IS+(8-IS)/6*8
        IX3 = IX+INX(ISBIG-1)
        IY3 = IY+INY(ISBIG-1)
        IX4 = IX+INX(ISBIG-2)
        IY4 = IY+INY(ISBIG-2)
        IPENO = IPEN
        IF (ISS .NE. 0) GO TO 110
        IF (IX3.GT.M .OR. IY3.GT.N .OR. IX3.LT.1 .OR. IY3.LT.1) GO TO 190
        IF (IX4.GT.M .OR. IY4.GT.N .OR. IX4.LT.1 .OR. IY4.LT.1) GO TO 190
    110 IPEN = Z(IX,IY).NE.SPVAL .AND. Z(IX2,IY2).NE.SPVAL .AND.
       1       Z(IX3,IY3).NE.SPVAL .AND. Z(IX4,IY4).NE.SPVAL
    120 IF (IDX .EQ. 0) GO TO 130
        Y = IY
        ISUB = IX+IDX
        X = C(Z(IX,IY),Z(ISUB,IY))*FLOAT(IDX)+FLOAT(IX)
        GO TO 140
    130 X = IX
        ISUB = IY+IDY
        Y = C(Z(IX,IY),Z(IX,ISUB))*FLOAT(IDY)+FLOAT(IY)
    140 GO TO JUMP2,(150,160)
    150 IF (.NOT.IPEN) GO TO 170
        IF (IPENO) GO TO 160
  C
  C END OF LINE SEGMENT
  C
        CALL FRSTS (D1,D2,3)
        CALL FRSTS (FX(XOLD,YOLD),FY(XOLD,YOLD),1)
  C
  C CONTINUE LINE SEGMENT
  C
    160 CALL FRSTS (FX(X,Y),FY(X,Y),2)
    170 XOLD = X
        YOLD = Y
        IF (IS .NE. 1) GO TO 180
        NP = NP+1
        IF (NP .GT. NR) GO TO 190
        IR(NP) = IPXY(IX,IY)
    180 IF (ISS .EQ. 0) GO TO  60
        IF (IX.NE.IX0 .OR. IY.NE.IY0 .OR. IS.NE.IS0) GO TO  60
  C
  C END OF LINE
  C
    190 CALL FRSTS (D1,D2,3)
        RETURN
        END
        SUBROUTINE TR32 (X,Y,MX,MY)
        SAVE
  C
        COMMON /ISOSR1/ ISLBT      ,U          ,V          ,W
  C
  C A.S.F. FOR SCALING
  C
        SU(UTEMP) = UTEMP
        SV(VTEMP) = VTEMP
        SW(WTEMP) = WTEMP
  C
        XX = X
        YY = Y
        IF (ISLBT)  10, 20, 30
     10 CALL TRN32I (SU(U),SV(XX-1.),SW(YY-1.),XT,YT,DUM,2)
        GO TO  40
     20 CALL TRN32I (SU(XX-1.),SV(V),SW(YY-1.),XT,YT,DUM,2)
        GO TO  40
     30 CALL TRN32I (SU(XX-1.),SV(YY-1.),SW(W),XT,YT,DUM,2)
     40 MX = XT
        MY = YT
        RETURN
        END
        SUBROUTINE FRSTS (XX,YY,IENT)
  C
  C THIS IS A SPECIAL VERSION OF THE SMOOTHING DASHED LINE PACKAGE.  LINES
  C ARE SMOOTHED IN THE SAME WAY, BUT NO SOFTFARE DASHED LINES ARE USED.
  C CONDITIONAL PLOTTING ROUTINES ARE CALL WHICH DETERMINE THE VISIBILITY
  C OF A LINE SEGMENT BEFORE PLOTTING.
  C
        SAVE
        DIMENSION       XSAVE(70)  ,YSAVE(70)  ,XP(70)     ,YP(70)     ,
       1                TEMP(70)
  C
        COMMON /ISOSR7/ IENTRY     ,IONES
  C
        DATA NP/150/
        DATA L1/70/
        DATA TENSN/2.5/
        DATA PI/3.14159265358/
        DATA SMALL/128./
  C
        AVE(A,B) = .5*(A+B)
  C
  C  DECIDE IF FRSTS,VECTS,LASTS CALL
  C
        GO TO ( 10, 20, 40),IENT
     10 DEG = 180./PI
        X = XX
        Y = YY
        LASTFL = 0
        SSLP1 = 0.0
        SSLPN = 0.0
        XSVN = 0.0
        YSVN = 0.0
  C
  C INITIALIZE THE POINT AND SEGMENT COUNTER
  C N COUNTS THE NUMBER OF POINTS/SEGMENT
  C
        N = 0
  C
  C NSEG = 0       FIRST SEGMENT
  C NSEG = 1       MORE THAN ONE SEGMENT
  C
        NSEG = 0
        CALL TR32 (X,Y,MX,MY)
  C
  C SAVE THE X,Y COORDINATES OF THE FIRST POINT
  C XSV1           CONTAINS THE X COORDINATE OF THE FIRST POINT
  C                OF A LINE
  C YSV1           CONTAINS THE Y COORDINATE OF THE FIRST POINT
  C                OF A LINE
  C
        XSV1 = MX
        YSV1 = MY
        GO TO  30
  C
  C ************************* ENTRY VECTS *************************
  C     ENTRY VECTS (XX,YY)
  C
     20 X = XX
        Y = YY
  C
  C VECTS          SAVES THE X,Y COORDINATES OF THE ACCEPTED
  C                POINTS ON A LINE SEGMENT
  C
        CALL TR32 (X,Y,MX,MY)
  C
  CIF THE NEW POINT IS TOO CLOSE TO THE PREVIOUS POINT, IGNORE IT
  C
        IF (ABS(FLOAT(IFIX(XSVN)-MX))+ABS(FLOAT(IFIX(YSVN)-MY)) .LT.
       1    SMALL) RETURN
        IFLAG = 0
     30 N = N+1
  C
  C SAVE THE X,Y COORDINATES OF EACH POINT OF THE SEGMENT
  C XSAVE          THE ARRAY OF X COORDINATES OF LINE SEGMENT
  C YSAVE          THE ARRAY OF Y COORDINATES OF LINE SEGMENT
  C
        XSAVE(N) = MX
        YSAVE(N) = MY
        XSVN = XSAVE(N)
        YSVN = YSAVE(N)
        IF (N .GE. L1-1) GO TO  50
        RETURN
  C
  C ************************* ENTRY LASTS *************************
  C     ENTRY LASTS
  C
     40 LASTFL = 1
  C
  C LASTS          CHECKS FOR PERIODIC LINES AND SETS UP
  C                  THE CALLS TO KURV1S AND KURV2S
  C
  C IFLAG = 0      OK TO CALL LASTS DIRECTLY
  C IFLAG = 1      LASTS WAS JUST CALLED FROM BY VECTS
  C                IGNORE CALL TO LASTS
  C
        IF (IFLAG .EQ. 1) RETURN
  C
  C COMPARE THE LAST POINT OF SEGMENT WITH FIRST POINT OF LINE
  C
     50 IFLAG = 1
  C
  C IPRD = 0       PERIODIC LINE
  C IPRD = 1       NON-PERIODIC LINE
  C
        IPRD = 1
        IF (ABS(XSV1-XSVN)+ABS(YSV1-YSVN) .LT. SMALL) IPRD = 0
  C
  C TAKE CARE OF THE CASE OF ONLY TWO DISTINCT P0INTS ON A LINE
  C
        IF (NSEG .GE. 1) GO TO  70
        IF (N-2) 160,150, 60
     60 IF (N .GE. 4) GO TO  70
        DX = XSAVE(2)-XSAVE(1)
        DY = YSAVE(2)-YSAVE(1)
        SLOPE = ATAN2(DY,DX)*DEG+90.
        IF (SLOPE .GE. 360.) SLOPE = SLOPE-360.
        IF (SLOPE .LE. 0.) SLOPE = SLOPE+360.
        SLP1 = SLOPE
        SLPN = SLOPE
        ISLPSW = 0
        SIGMA = TENSN
        GO TO 110
     70 SIGMA = TENSN
        IF (IPRD .GE. 1) GO TO  90
        IF (NSEG .GE. 1) GO TO  80
  C
  C SET UP FLAGS FOR A  1  SEGMENT, PERIODIC LINE
  C
        ISLPSW = 4
        XSAVE(N) = XSV1
        YSAVE(N) = YSV1
        GO TO 110
  C
  C SET UP FLAGS FOR AN N-SEGMENT, PERIODIC LINE
  C
     80 SLP1 = SSLPN
        SLPN = SSLP1
        ISLPSW = 0
        GO TO 110
     90 IF (NSEG .GE. 1) GO TO 100
  C
  C SET UP FLAGS FOR THE 1ST SEGMENT OF A NON-PERIODIC LINE
  C
        ISLPSW = 3
        GO TO 110
  C
  C SET UP FLAGS FOR THE NTH SEGMENT OF A NON-PERIODIC LINE
  C
    100 SLP1 = SSLPN
        ISLPSW = 1
  C
  C CALL THE SMOOTHING ROUTINES
  C
    110 CALL KURV1S (N,XSAVE,YSAVE,SLP1,SLPN,XP,YP,TEMP,S,SIGMA,ISLPSW)
        IF (IPRD.EQ.0 .AND. NSEG.EQ.0 .AND. S.LT.70.) GO TO 170
        IENTRY = 1
  C
  C DETERMINE THE NUMBER OF POINTS TO INTERPOLATE FOR EACH SEGMENT
  C
        IF (NSEG.GE.1 .AND. N.LT.L1-1) GO TO 120
        NPRIME = FLOAT(NP)-(S*FLOAT(NP))/(2.*32768.)
        IF (S .GE. 32768.) NPRIME = .5*FLOAT(NP)
        NPL = FLOAT(NPRIME)*S/32768.
        IF (NPL .LT. 2) NPL = 2
    120 DT = 1./FLOAT(NPL)
        IF (NSEG .LE. 0) CALL FRSTC (IFIX(XSAVE(1)),IFIX(YSAVE(1)),1)
        T = 0.0
        NSLPSW = 1
        IF (NSEG .GE. 1) NSLPSW = 0
        NSEG = 1
        CALL KURV2S (T,XS,YS,N,XSAVE,YSAVE,XP,YP,S,SIGMA,NSLPSW,SLP)
  C
  C SAVE SLOPE AT THE FIRST POINT OF THE LINE
  C
        IF (NSLPSW .GE. 1) SSLP1 = SLP
        NSLPSW = 0
        XSOLD = XSAVE(1)
        YSOLD = YSAVE(1)
        DO 130 I=1,NPL
           T = T+DT
           TT = -T
           IF (I .EQ. NPL) NSLPSW = 1
           CALL KURV2S (TT,XS,YS,N,XSAVE,YSAVE,XP,YP,S,SIGMA,NSLPSW,SLP)
  C
  C SAVE THE LAST SLOPE OF THIS LINE SEGMENT
  C
           IF (NSLPSW .GE. 1) SSLPN = SLP
  C
  C DRAW EACH PART OF THE LINE SEGMENT
  C
           CALL FRSTC (IFIX(AVE(XSOLD,XS)),IFIX(AVE(YSOLD,YS)),2)
           CALL FRSTC (IFIX(XS),IFIX(YS),2)
           XSOLD = XS
           YSOLD = YS
    130 CONTINUE
        IF (IPRD .NE. 0) GO TO 140
  C
  C CONNECT THE LAST POINT WITH THE FIRST POINT OF A PERIODIC LINE
  C
        CALL FRSTC (IFIX(AVE(XSOLD,XS)),IFIX(AVE(YSOLD,YS)),2)
        CALL FRSTC (IFIX(XSV1),IFIX(YSV1),2)
  C
  C BEGIN THE NEXT LINE SEGMENT WITH THE LAST POINT OF THIS SEGMENT
  C
    140 XSAVE(1) = XS
        YSAVE(1) = YS
        N = 1
    150 CONTINUE
    160 RETURN
    170 N = 0
        RETURN
        END
        SUBROUTINE FRSTC (MX,MY,IENT)
        SAVE
  C
        COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
       1                ISCA(8,128)
        COMMON /ISOSR4/ RX         ,RY
        COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
        LOGICAL         GENDON
        COMMON /ISOSR8/ NMASK(16)  ,IXOLD      ,IYOLD      ,IBTOLD     ,
       1                HBFLAG     ,IOSLSN     ,LRLX       ,IFSX       ,
       2                IFSY       ,FIRST      ,IYDIR      ,IHX        ,
       3                IHB        ,IHS        ,IHV        ,IVOLD      ,
       4                IVAL       ,IHRX       ,YCHANG     ,ITPD       ,
       5                IHF
        LOGICAL         YCHANG     ,HBFLAG     ,FIRST      ,IHF
  C
  C
  C  DRAW LINE TO THE POINT MX,MY
  C
  C  ENTER THE POINT INTO THE CURRENT SCREEN, ISCR, IF THE POINT CONFORMS
  C  TO THE SHADING ALGORITHM.
  C  THE POINT IS NOT ENTERED WHEN;
  C  1. IT IS THE SAME POINT USED IN THE LAST CALL, RESOLUTION PROBLEM
  C  2. IT IS PART OF A HORIZONTAL LINE BUT NOT AN END POINT
  C  3.  THE ENTIRE CONTOUR RESTS ON A HORIZONTAL PLANE
  C
  C  WHEN DRAWING A HORIZONTAL LINE THREE CONDITIONS EXIST;
  C  1. WHEN THE LINE IS A HORIZONTAL STEP ENTER ONLY THE OUTSIDE POINT.
  C     A HORIZONTAL STEP IS DEFINED BY THE ENTERING AND EXITING Y
  C     DIRECTION THAT IS THE SAME.
  C  2. ENTER BOTH END POINTS OF A HORIZONTAL TURNING POINT. A HORIZONTAL
  C     TURNING POINT IS A LINE WITH GREATER THAN 1 HORIZONTAL BITS
  C     AND THE ENTERING AND EXITING Y DIRECTION IS DIFFIRENT.
  C  3.  WHEN THE ENTIRE CONTOUR IS A HORIZONTAL LINE NO POINTS ARE
  C      ENTERED. THIS CONDITION IS DETECTED BY THE STATUS OF YCHANG.
  C      IF IT IS TRUE THEN THE CONTOUR IS NOT A SINGLE HORIZONTAL LINE.
  C
  C  THE PREVIOUS POINT IS ERASED IF IT IS A VERTICAL TURNING POINT.
  C  A VERTICAL TURNING POINT IS A HORIZONTAL LINE WITH ONLY 1 POINT
  C  AND THE ENTERING AND EXITING Y DIRECTION DIFFERS.THIS DATA IS
  C  IN THE VARIABLES IOSLSN-OLD SLOPE AND ISLSGN-NEW SLOPE.
  C  THE CHANGE IN SLOPE MUST BE -1 TO 1 OR 1 TO -1.
  C
  C  OTHERWISE THE POINT IS ENTERED INTO ISCR.
  C
  C  THE TWO ENTRY POINTS ARE REQUIRED BY THE HARDWARE DRAWING ROUTINES.
  C  FIRSTC IS USED FOR THE FIRST POINT ON THE CONTOUR. THE REMAINING
  C  POINTS ON THE SAME CONTOUR ARE ENTERED VIA VECTC.
  C
        DATA IONE/1/
        AVE(A,B) = (A+B)*.5
  C
  C  COMPUTE VISIBILITY OF THIS POINT
  C
  C  WARNING
  C  IF X OR Y PLOTTER MAXIMUM VALUE RANGES FALL BELOW 101 THEN THE
  C  FOLLOWING TWO STATEMENTS WHICH SET IX AND IY MUST BE CHANGED.
  C  REPLACE THE CONSTANT 1.0 BY 0.5 IN THE STATEMENTS WHERE THE
  C  MAXIMUM PLOTTER VALUE IS LESS THAN 101 FOR THAT DIRECTION.  THE
  C  PLOTTER CORDINATE RANGES ARE SET IN SET32.
  C
        IX = FLOAT(MX-1)*RX+1.0
        NRLX = IX
        IY = FLOAT(MY-1)*RY+1.0
        IBIT = NBPW-MOD(IX,NBPW)
        IX = IX/NBPW+1
        IVNOW = IAND(ISHIFT(ISCA(IX,IY),1-IBIT),IONE)
  C
  C  DECIDE IF FRSTC OR VECTC CALL
  C
        IF (IENT .NE. 1) GO TO  10
  C
        XOLD = MX
        YOLD = MY
  C
  C
  C  SET INITIAL VALUES
  C
        IHF = .FALSE.
        IYDIR = 0
        ITPD = 0
        IVAL = 0
        IOSLSN = 0
        IFSX = NRLX
        IFSY = IY
        LASTV = IVNOW
        HBFLAG = .FALSE.
        YCHANG = .FALSE.
        CALL PLOTIT (IFIX(XOLD),IFIX(YOLD),0)
        GO TO 180
  C
  C****************************  ENTRY VECTC  ****************************
  C     ENTRY VECTC (MX,MY)
  C
     10 XNOW = MX
        YNOW = MY
        JUMP = IVNOW*2+LASTV+1
        GO TO ( 20, 30, 40, 50),JUMP
  C
  C BOTH VISIBLE
  C
     20 CALL PLOTIT (IFIX(XNOW),IFIX(YNOW),1)
        GO TO  50
  C
  C JUST TURNED VISIBLE
  C
     30 CALL PLOTIT (IFIX(AVE(XNOW,XOLD)),IFIX(AVE(YNOW,YOLD)),0)
        GO TO  50
  C
  C JUST TURNED INVISIBLE
  C
     40 CALL PLOTIT (IFIX(AVE(XNOW,XOLD)),IFIX(AVE(YNOW,YOLD)),1)
  C
  C BOTH INVISIBLE
  C
     50 XOLD = XNOW
        YOLD = YNOW
        LASTV = IVNOW
  C
  C  TEST FOR RESOLUTION PROBLEM
  C
        IF (NRLX.EQ.LRLX .AND. IY.EQ.IYOLD) RETURN
  C
  C  TEST FOR HORIZONTAL BITS
  C
        IF (IYOLD .NE. IY) GO TO  70
  C
  C  HORIZONTAL BITS DETECTED. SET FLAG AND EXIT.
  C  THIS AND THE NEXT HORIZONTAL BIT TEST IS NECESSARY FOR ISCR TO
  C  CONFORM TO THE SHADING ALGORITHM IN SUBROUTINE FILLIN
  C
  C
  C  IF HORIZONTAL LINE PREVIOUSLY DETECTED EXIT
  C
        IF (.NOT.HBFLAG) GO TO  60
  C
  C  IF END OF CONTOUR ON A HORIZONTAL LINE BRANCH FOR SPECIAL PROCESSING.
  C
        IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) GO TO 210
        GO TO 200
  C
  C  SAVE SLOPE PRIOR TO HORIZONTAL LINE
  C
     60 IHX = IXOLD
        IHB = IBTOLD
        IHS = IOSLSN
        IOSLSN = 0
        HBFLAG = .TRUE.
        IHRX = LRLX
        IHV = IVOLD
        IF (LRLX.EQ.IFSX .AND. IYOLD.EQ.IFSY) IHF = .TRUE.
  C
  C  THIS IS THE SECOND TRAP FOR END OF CONTOUR ON A HORIZONTAL LINE.
  C
        IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) GO TO 210
        GO TO 200
  C
  C  COMPUTE THE SLOPE TO THIS POINT
  C
     70 IF (IY-IYOLD)  80, 90,100
     80 ISLSGN = 1
        GO TO 110
     90 ISLSGN = 0
        GO TO 120
    100 ISLSGN = -1
    110 IF (IYDIR .EQ. 0) IYDIR = ISLSGN
    120 CONTINUE
  C
  C  IF PROCESS REACHES THIS CODE THE CONTOUR IS NOT CONTAINED ON A SINGLE
  C  HORIZONTAL  PLANE, SO RECORD THIS FACT BY SETTING Y CHANGE FLAG.
  C
        YCHANG = .TRUE.
  C
  C  TEST FOR END OF HORIZONTAL LINE
  C
        IF (.NOT.HBFLAG) GO TO 160
        HBFLAG = .FALSE.
  C
  C  HORIZONTAL LINE JUST ENDED
  C
  C  TEST FOR REDRAW
  C
        ITEMP = IAND(ISCR(IXOLD,IYOLD),MASK(IBTOLD))
        IF ((IHV .EQ. 0) .AND. (ITEMP .EQ. 0)) GO TO 130
  C
  C  REDRAWING ERASE THIS POINT
  C
        ISCR(IXOLD,IYOLD) = IAND(ISCR(IXOLD,IYOLD),NMASK(IBTOLD))
        ISCR(IHX,IYOLD) = IAND(ISCR(IHX,IYOLD),NMASK(IHB))
        GO TO 170
  C
  C  TEST FOR STEP PROBLEM
  C
    130 IF (IHS .NE. ISLSGN) GO TO 140
  C
  C  STEP PROBLEM
  C
        GO TO 170
  C
  C  TURNING PROBLEM HORIZONTAL LINE IS A TURNING POINT
  C
    140 CONTINUE
  C
  C  ENTER THE TURNING POINT ONLY IF IT IS NOT THE SECOND SUCCEEDING
  C  EVENT IN A ROW
  C
        ICTPD = 1
        IF (IHRX .GT. NRLX) ICTPD = -1
        IF (ICTPD .NE. ITPD) GO TO 150
        ITPD = 0
  C
  C  ERASE THE FIRST POINT
  C
        ISCR(IHX,IYOLD) = IAND(ISCR(IHX,IYOLD),NMASK(IHB))
        GO TO 170
  C
  C  ENTER THE TURNING POINT
  C
    150 CONTINUE
        ITPD = ICTPD
  C
  C  ENTER THE SECOND POINT
  C
        ISCR(IXOLD,IYOLD) = IOR(ISCR(IXOLD,IYOLD),MASK(IBTOLD))
        GO TO 170
  C
  C  CHECK IF PREVIOUS ENTRY WAS A VERTICAL TURNING POINT.
  C  IF SO ERASE IT.
  C
    160 IF (ISLSGN.EQ.IOSLSN .OR. (IOSLSN.EQ.0 .OR. ISLSGN.EQ.0))
       1    GO TO 170
        ITPD = 0
        ISCR(IXOLD,IYOLD) = IAND(ISCR(IXOLD,IYOLD),NMASK(IBTOLD))
  C
    170 IOSLSN = ISLSGN
  C
  C  CHECK IF THIS GRID POINT PREVIOUSLY ACTIVATED
  C
        IVAL = IAND(ISCR(IX,IY),MASK(IBIT))
  C
  C  IF GRID POINTS ACTIVATED BRANCH
  C
        IF (IVAL .NE. 0) GO TO 190
  C
  C  GRID POINT NOT ACTIVATED   SET AND EXIT
  C
    180 CONTINUE
        ISCR(IX,IY) = IOR(ISCR(IX,IY),MASK(IBIT))
        GO TO 200
  C
  C  THIS POINT IS BEING REDRAWN SO ERASE IT.
  C  (THIS IS TO CONFORM WITH THE SHADING ALGORITHM, FILLIN.
  C  HOWEVER IF BACK TO STARTING POINT DO NOT ERASE
  C
    190 IF (NRLX.EQ.IFSX .AND. IY.EQ.IFSY) RETURN
        ISCR(IX,IY) = IAND(ISCR(IX,IY),NMASK(IBIT))
  C
  C
    200 IXOLD = IX
        LRLX = NRLX
        IYOLD = IY
        IBTOLD = IBIT
        IVOLD = IVAL
        RETURN
  C
  C  PERFORM THIS OPERATION WHEN A CONTOUR STARTS OR ENDS ON A HORIZONTAL
  C  LINE.
  C
    210 CONTINUE
  C
  C  ERASE THE FIRST POINT OF A CONTOUR WHEN IT IS PART OF A HORIZONTAL
  C  LINE SEGMENT AND IS NOT THE ENDPOINT OF THE SEGMENT
  C
        IF (.NOT.IHF) GO TO 220
        ISCR(IX,IY) = IAND(ISCR(IX,IY),NMASK(IBIT))
    220 CONTINUE
  C
  C  ERASE THE FIRST POINT OF A HORIZONTAL LINE SEGMENT WHEN IT ENDS
  C  THE CONTOUR AND IS NOT THE HIGHEST LINE SEG ON THS SIDE.
  C
        IF (.NOT.YCHANG) GO TO 230
        IF (IYDIR .NE. IHS) GO TO 200
    230 ISCR(IHX,IY) = IAND(ISCR(IHX,IY),NMASK(IHB))
        GO TO 200
        END
        SUBROUTINE FILLIN
  C
        SAVE
        COMMON /ISOSR2/ LX         ,NX         ,NY         ,ISCR(8,128),
       1                ISCA(8,128)
        COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
        LOGICAL         GENDON
        COMMON /ISOSR7/ IENTRY     ,IONES
  C
        IF (IENTRY .EQ. 0) RETURN
  C
  C  THIS IS A SHADING ALGORITHM IT IS USED TO DETERMINE CONTOUR LINES
  C  THAT ARE HIDDEN BY THE PRESENT LINE. THE ALGORITHM PROCESSES
  C  HORIZONTAL ROWS. IT ASSUMES THAT THE BIT PATTERN PASSED TO IT
  C  HAS ONLY BITS SET TO MARK THE START AND END OF SHADING. THE
  C  ALGORITHM ALSO ASSUMES THAT WHEN AN ON BIT IS ENCOUNTERED THAT A
  C  CORRESPONDING OFF BIT IS INCLUDED IN THE SAME ROW.
  C
  C
  C  PULL OUT ROWS OF THE CONTOUR PATTERN
  C
        IBVAL = 0
        DO  80 IYNOW=1,NY
           DO  40 IXNOW=1,LX
  C
  C  IF NO ACTIVATED BITS BRANCH
  C
              ICRWD = ISCR(IXNOW,IYNOW)
              IF (ICRWD .EQ. 0) GO TO  30
  C
  C  ACTIVATED BITS IN WORD SET SHADING FLAG
  C
  C  CHECK BIT BY BIT FOR ON/OFF FLAGS
  C
              DO  20 IB=1,NBPW
                 IBIT = (NBPW+1)-IB
  C
  C
  C  PULL OUT THE CURRENT GRID POINT VALUE
  C
                 IVAL = IAND(ICRWD,MASK(IBIT))
  C
  C  IF IVAL SET, THIS IS AN ON/OFF FLAG
  C
                 IF (IVAL .EQ. 0) GO TO  10
  C
  C  FLAG BIT, ALWAYS SET
  C
                 IBVAL = MOD(IBVAL+1,2)
                 GO TO  20
  C
  C  SHADE THE SCREEN ACCORDING TO THE STATUS OF IBVAL
  C
     10          IF (IBVAL .NE. 0) ICRWD = IOR(ICRWD,MASK(IBIT))
  C
     20       CONTINUE
  C
  C  ZERO OUT THE SCREEN
  C
              ISCR(IXNOW,IYNOW) = 0
              ISCA(IXNOW,IYNOW) = IOR(ICRWD,ISCA(IXNOW,IYNOW))
              GO TO  40
  C
     30       IF (IBVAL .NE. 0) ISCA(IXNOW,IYNOW) = IONES
     40    CONTINUE
  C
  C  FIX FOR NONCORRECTABLE RUNAWAYS
  C
           IF (IBVAL .EQ. 0) GO TO  80
           IBVAL = 0
           DO  70 K=1,LX
              ITEST = 0
              IF (IYNOW .EQ. 1) GO TO  50
              ITEST = ISCA(K,IYNOW-1)
              IF (IYNOW .EQ. NY) GO TO  60
     50       ITEST = IOR(ITEST,ISCA(K,IYNOW+1))
     60       ISCA(K,IYNOW) = ITEST
     70    CONTINUE
  C
     80 CONTINUE
        RETURN
        END
        SUBROUTINE DRAWI (IXA,IYA,IXB,IYB)
  C
  C INCLUDED FOR USE BY PWRZ
  C
        SAVE
        CALL FRSTC (IXA,IYA,1)
        CALL FRSTC (IXB,IYB,2)
        RETURN
        END
        SUBROUTINE MMASK
  C
  C  MAKE THE MACHINE DEPENDENT MASKS USED IN THE CONTOUR DRAWING
  C  AND SHADING ALGORITHMS
  C
        SAVE
        COMMON /ISOSR5/ NBPW       ,MASK(16)   ,GENDON
        LOGICAL         GENDON
        COMMON /ISOSR7/ IENTRY     ,IONES
        COMMON /ISOSR8/ NMASK(16)  ,IXOLD      ,IYOLD      ,IBTOLD     ,
       1                HBFLAG     ,IOSLSN     ,LRLX       ,IFSX       ,
       2                IFSY       ,FIRST      ,IYDIR      ,IHX        ,
       3                IHB        ,IHS        ,IHV        ,IVOLD      ,
       4                IVAL       ,IHRX       ,YCHANG     ,ITPD       ,
       5                IHF
        COMMON /ISOSR9/ BIG        ,IXBIT
        LOGICAL         YCHANG     ,HBFLAG     ,FIRST      ,IHF
        GENDON = .TRUE.
        NBPW = 16
  C
  C  GET BIGGEST REAL NUMBER
  C
        BIG = R1MACH(2)
  C
  C  MASKS TO SELECT A SPECIFIC BIT
  C
        DO  10 K=1,NBPW
           MASK(K) = ISHIFT(1,K-1)
     10 CONTINUE
  C
  C  GENERATE  THE BIT PATTERN 177777 OCTAL
  C
        ITEMP1 = 0
        ITEMP = MASK(NBPW)
        IST = NBPW-1
        DO  20 K=1,IST
           ITEMP1 = IOR(ITEMP,ISHIFT(ITEMP1,-1))
     20 CONTINUE
        MFIX = IOR(ITEMP1,1)
  C
  C  MASKS TO CLEAR A SPECIFIC BIT
  C
        DO  30 K=1,NBPW
           NMASK(K) = IAND(ITEMP1,MFIX)
           ITEMP1 = IOR(ISHIFT(ITEMP1,1),1)
     30 CONTINUE
        IONES = MFIX
        RETURN
  C
  C REVISION HISTORY---
  C
  C JANUARY 1978     DELETED REFERENCES TO THE  *COSY  CARDS AND
  C                  ADDED REVISION HISTORY
  C JANUARY 1979     NEW SHADING ALGORITHM
  C MARCH 1979       MADE CODE MACHINE INDEPENDENT AND CONFORM
  C                  TO 66 FORTRAN STANDARD
  C JUNE 1979        THIS VERSION PLACED ON ULIB.
  C SEPTEMBER 1979   FIXED PROBLEM IN EZISOS DEALING WITH
  C                  DETERMINATION OF VISIBILITY OF W PLANE.
  C DECEMBER 1979    FIXED PROBLEM WITH PEN DOWN ON CONTOUR
  C                  INITIALIZATION IN SUBROUTINE FRSTC
  C MARCH            CHANGED ROUTINE NAMES  TRN32I  AND  DRAW  TO
  C                  TRN32I  AND  DRAWI  TO BE CONSISTENT WITH THE
  C                  USAGE OF THE NEW ROUTINE  PWRZI.
  C JUNE  1980       FIXED PROBLEM WITH ZERO INDEX COMPUTATION IN
  C                  SUBROUTINE FRSTC.  ADDED INPUT PARAMETER
  C                  DIMENSION STATEMENT MISSING IN EZISOS.
  C                  FIXED ERROR IN COMPUTATION OF ARCCOSINE
  C                  IN EZISOS AND TRN32I.
  C DECEMBER 1984    CONVERTED TO GKS LEVEL 0A AND STANDARD FORTRAN 77
  C-----------------------------------------------------------------------
  C
        END
  