path: root/sys/gio/ncarutil/strmln.f

SUBROUTINE STRMLN (U,V,WORK,IMAX,IPTSX,JPTSY,NSET,IER)
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
C SUBROUTINE STRMLN (U,V,WORK,IMAX,IPTSX,JPTSY,NSET,IER)
C
C DIMENSION OF           U(IMAX,JPTSY) , V(IMAX,JPTSY) ,
C ARGUMENTS              WORK(2*IMAX*JPTSY)
C
C LATEST REVISION        JUNE 1984
C
C PURPOSE                STRMLN DRAWS A STREAMLINE REPRESENTATION OF
C                        THE FLOW FIELD. THE REPRESENTATION IS
C                        INDEPENDENT OF THE FLOW SPEED.
C
C USAGE                  IF THE FOLLOWING ASSUMPTIONS ARE MET, USE
C
C                            CALL EZSTRM  (U,V,WORK,IMAX,JMAX)
C
C                          ASSUMPTIONS:
C                            --THE WHOLE ARRAY IS TO BE PROCESSED.
C                            --THE ARRAYS ARE DIMENSIONED
C                              U(IMAX,JMAX) , V(IMAX,JMAX) AND
C                              WORK(2*IMAX*JMAX).
C                            --WINDOW AND VIEWPORT ARE TO BE CHOSEN
C                              BY STRMLN.
C                            --PERIM IS TO BE CALLED.
C
C                        IF THESE ASSUMPTIONS ARE NOT MET, USE
C
C                            CALL STRMLN (U,V,WORK,IMAX,IPTSX,JPTSY,
C                                         NSET,IER)
C
C                        THE USER MUST CALL FRAME IN THE CALLING
C                        ROUTINE.
C
C                        THE USER MAY CHANGE VARIOUS INTERNAL
C                        PARAMETERS VIA COMMON BLOCKS. SEE BELOW.
C
C ARGUMENTS
C
C ON INPUT               U, V
C                          TWO DIMENSIONAL ARRAYS CONTAINING THE
C                          VELOCITY FIELDS TO BE PLOTTED.
C                          (NOTE:  IF THE U AND V COMPONENTS
C                          ARE, FOR EXAMPLE, DEFINED IN CARTESIAN
C                          COORDINATES AND THE USER WISHES TO PLOT THEM
C                          ON A DIFFERENT PROJECTION (I.E., STEREO-
C                          GRAPHIC), THEN THE APPROPRIATE
C                          TRANSFORMATION MUST BE MADE TO THE U AND V
C                          COMPONENTS VIA THE FUNCTIONS FU AND FV
C                          (LOCATED IN DRWSTR).
C
C                        WORK
C                          USER PROVIDED WORK ARRAY.  THE DIMENSION
C                          OF THIS ARRAY MUST BE .GE. 2*IMAX*JPTSY.
C                          CAUTION:  THIS ROUTINE DOES NOT CHECK THE
C                          SIZE OF THE WORK ARRAY.
C
C                        IMAX
C                          THE FIRST DIMENSION OF U AND V IN THE
C                          CALLING PROGRAM. (X-DIRECTION)
C
C                        IPTSX
C                          THE NUMBER OF POINTS TO BE PLOTTED IN THE
C                          FIRST SUBSCRIPT DIRECTION.  (X-DIRECTION)
C
C                        JPTSY
C                          THE NUMBER OF POINTS TO BE PLOTTED IN THE
C                          SECOND SUBSCRIPT DIRECTION. (Y-DIRECTION)
C
C                        NSET
C                          FLAG TO CONTROL SCALING
C                          > 0  STRMLN ASSUMES THAT THE WINDOW
C                               AND VIEWPORT HAVE BEEN SET BY THE
C                               USER IN SUCH A WAY AS TO PROPERLY
C                               SCALE THE PLOTTING INSTRUCTIONS
C                               GENERATED BY STRMLN. PERIM IS NOT
C                               CALLED.
C                          = 0  STRMLN WILL ESTABLISH THE WINDOW AND
C                               VIEWPORT TO PROPERLY SCALE THE
C                               PLOTTING INSTRUCTIONS TO THE STANDARD
C                               CONFIGURATION. PERIM IS CALLED TO DRAW
C                               THE BORDER.
C                          < 0  STRMLN ESTABLISHES THE WINDOW
C                               AND VIEWPORT SO AS TO PLACE THE
C                               STREAMLINES WITHIN THE LIMITS
C                               OF THE USER'S WINDOW.  PERIM IS
C                               NOT CALLED.
C
C ON OUTPUT              ONLY THE IER ARGUMENT MAY BE CHANGED. ALL
C                        OTHER ARGUMENTS ARE UNCHANGED.
C
C
C                        IER
C                          =  0 WHEN NO ERRORS ARE DETECTED
C                          = -1 WHEN THE ROUTINE IS CALLED WITH ICYC
C                               .NE. 0  AND THE DATA ARE NOT CYCLIC
C                               (ICYC IS AN INTERNAL PARAMETER
C                               DESCRIBED BELOW); IN THIS CASE THE
C                               ROUTINE WILL DRAW THE
C                               STREAMLINES WITH THE NON-CYCLIC
C                               INTERPOLATION FORMULAS.
C
C ENTRY POINTS           STRMLN, DRWSTR, EZSTRM, GNEWPT, CHKCYC
C
C COMMON BLOCKS          STR01, STR02, STR03, STR04
C
C REQUIRED LIBRARY       GRIDAL, GBYTES, AND THE SPPS
C ROUTINES
C
C HISTORY                WRITTEN AND STANDARDIZED IN NOVEMBER 1973.
C I/O                    DRAWS STREAMLINES
C
C PRECISION              SINGLE
C
C LANGUAGE               FORTRAN
C
C HISTORY                WRITTEN IN 1979.
C                        CONVERTED TO FORTRAN 77 AND GKS IN JUNE 1984.
C
C PORTABILITY            FORTRAN 77
C
C ALGORITHM              WIND COMPONENTS ARE NORMALIZED TO THE VALUE
C                        OF DISPL. THE LEAST SIGNIFICANT TWO
C                        BITS OF THE WORK ARRAY ARE
C                        UTILIZED AS FLAGS FOR EACH GRID BOX. FLAG 1
C                        INDICATES WHETHER ANY STREAMLINE HAS
C                        PREVIOUSLY PASSED THROUGH THIS BOX.  FLAG 2
C                        INDICATES WHETHER A DIRECTIONAL ARROW HAS
C                        ALREADY APPEARED IN A BOX. JUDICIOUS USE
C                        OF THESE FLAGS PREVENTS OVERCROWDING OF
C                        STREAMLINES AND DIRECTIONAL ARROWS.
C                        EXPERIENCE INDICATES THAT A FINAL PLEASING
C                        PICTURE IS PRODUCED WHEN STREAMLINES ARE
C                        INITIATED IN THE CENTER OF A GRID BOX. THE
C                        STREAMLINES ARE DRAWN IN ONE DIRECTION THEN
C                        IN THE OPPOSITE DIRECTION.
C
C REFERENCE              THE TECHNIQUES UTILIZED HERE ARE DESCRIBED
C                        IN AN ARTICLE BY THOMAS WHITTAKER (U. OF
C                        WISCONSIN) WHICH APPEARED IN THE NOTES AND
C                        CORRESPONDENCE SECTION OF MONTHLY WEATHER
C                        REVIEW, JUNE 1977.
C
C TIMING                 HIGHLY VARIABLE
C                          IT DEPENDS ON THE COMPLEXITY OF THE
C                          FLOW FIELD AND THE PARAMETERS:  DISPL,
C                          DISPC , CSTOP , INITA , INITB , ITERC ,
C                          AND IGFLG. (SEE BELOW FOR A DISCUSSION
C                          OF THESE PARAMETERS.) IF ALL VALUES
C                          ARE DEFAULT, THEN A SIMPLE LINEAR
C                          FLOW FIELD FOR A 40 X 40 GRID WILL
C                          TAKE ABOUT 0.4 SECONDS ON THE CRAY1-A;
C                          A FAIRLY COMPLEX FLOW FIELD WILL TAKE ABOUT
C                          1.5 SECONDS ON THE CRAY1-A.
C
C
C INTERNAL PARAMETERS
C
C                        NAME     DEFAULT            FUNCTION
C                        ----     -------            --------
C
C                        EXT       0.25      LENGTHS OF THE SIDES OF THE
C                                            PLOT ARE PROPORTIONAL TO
C                                            IPTSX AND JPTSY EXCEPT IN
C                                            THE CASE WHEN MIN(IPTSX,JPT
C                                            / MAX(IPTSX,JPTSY) .LT. EXT;
C                                            IN THAT CASE A SQUARE
C                                            GRAPH IS PLOTTED.
C
C                        SIDE      0.90      LENGTH OF LONGER EDGE OF
C                                            PLOT. (SEE ALSO EXT.)
C
C                        XLT       0.05      LEFT HAND EDGE OF THE PLOT.
C                                            (0.0 = LEFT EDGE OF FRAME)
C                                            (1.0 = RIGHT EDGE OF FRAME)
C
C                        YBT       0.05      BOTTOM EDGE OF THE PLOT.
C                                            (0.0 = BOTTOM ; 1.0 = TOP)
C
C                                            (YBT+SIDE AND XLT+SIDE MUST
C                                            BE .LE. 1. )
C
C                        INITA     2         USED TO PRECONDITION GRID
C                                            BOXES TO BE ELIGIBLE TO
C                                            START A STREAMLINE.
C                                            FOR EXAMPLE, A VALUE OF 4
C                                            MEANS THAT EVERY FOURTH
C                                            GRID BOX IS ELIGIBLE ; A
C                                            VALUE OF 2 MEANS THAT EVERY
C                                            OTHER GRID BOX IS ELIGIBLE.
C                                            (SEE INITB)
C
C                        INITB     2         USED TO PRECONDITION GRID
C                                            BOXES TO BE ELIGIBLE FOR
C                                            DIRECTION ARROWS.
C                                            IF THE USER CHANGES THE
C                                            DEFAULT VALUES OF INITA
C                                            AND/OR INITB, IT SHOULD
C                                            BE DONE SUCH THAT
C                                            MOD(INITA,INITB) = 0 .
C                                            FOR A DENSE GRID TRY
C                                            INITA=4 AND INITB=2 TO
C                                            REDUCE THE CPU TIME.
C
C                        AROWL     0.33      LENGTH OF DIRECTION ARROW.
C                                            FOR EXAMPLE, 0.33 MEANS
C                                            EACH DIRECTIONAL ARROW WILL
C                                            TAKE UP A THIRD OF A GRID
C                                            BOX.
C
C                        ITERP     35        EVERY 'ITERP' ITERATIONS
C                                            THE STREAMLINE PROGRESS
C                                            IS CHECKED.
C
C                        ITERC     -99       THE DEFAULT VALUE OF THIS
C                                            PARAMETER IS SUCH THAT
C                                            IT HAS NO EFFECT ON THE
C                                            CODE. WHEN SET TO SOME
C                                            POSITIVE VALUE, THE PROGRAM
C                                            WILL CHECK FOR STREAMLINE
C                                            CROSSOVER EVERY 'ITERC'
C                                            ITERATIONS. (THE ROUTINE
C                                            CURRENTLY DOES THIS EVERY
C                                            TIME IT ENTERS A NEW GRID
C                                            BOX.) CAUTION:  WHEN
C                                            THIS PARAMETER IS ACTIVATED
C                                            CPU TIME WILL INCREASE.
C
C                        IGFLG     0         A VALUE OF ZERO MEANS THAT
C                                            THE SIXTEEN POINT BESSEL
C                                            INTERPOLATION FORMULA WILL
C                                            BE UTILIZED WHERE POSSIBLE;
C                                            WHEN NEAR THE GRID EDGES,
C                                            QUADRATIC AND BI-LINEAR
C                                            INTERPOLATION  WILL BE
C                                            USED. THIS MIXING OF
C                                            INTERPOLATION SCHEMES CAN
C                                            SOMETIMES CAUSE SLIGHT
C                                            RAGGEDNESS NEAR THE EDGES
C                                            OF THE PLOT.  IF IGFLG.NE.0,
C                                            THEN ONLY THE BILINEAR
C                                            INTERPOLATION FORMULA
C                                            IS USED; THIS WILL GENERALLY
C                                            RESULT IN SLIGHTLY FASTER
C                                            PLOT TIMES BUT A LESS
C                                            PLEASING PLOT.
C
C                        IMSG      0         IF ZERO, THEN NO MISSING
C                                            U AND V COMPONENTS ARE
C                                            PRESENT.
C                                            IF .NE. 0, STRMLN WILL
C                                            UTILIZE THE
C                                            BI-LINEAR INTERPOLATION
C                                            SCHEME AND TERMINATE IF
C                                            ANY DATA POINTS ARE MISSING.
C
C                        UVMSG     1.E+36    VALUE ASSIGNED TO A MISSING
C                                            POINT.
C
C                        ICYC      0         ZERO MEANS THE DATA ARE
C                                            NON-CYCLIC IN THE X
C                                            DIRECTION.
C                                            IF .NE 0, THE
C                                            CYCLIC INTERPOLATION
C                                            FORMULAS WILL BE USED.
C                                            (NOTE:  EVEN IF THE DATA
C                                            ARE CYCLIC IN X LEAVING
C                                            ICYC = 0 WILL DO NO HARM.)
C
C                        DISPL     0.33      THE WIND SPEED IS
C                                            NORMALIZED TO THIS VALUE.
C                                            (SEE THE DISCUSSION BELOW.)
C
C                        DISPC     0.67      THE CRITICAL DISPLACEMENT.
C                                            IF AFTER 'ITERP' ITERATIONS
C                                            THE STREAMLINE HAS NOT
C                                            MOVED THIS DISTANCE, THE
C                                            STREAMLINE WILL BE
C                                            TERMINATED.
C
C                        CSTOP     0.50      THIS PARAMETER CONTROLS
C                                            THE SPACING BETWEEN
C                                            STREAMLINES.  THE CHECKING
C                                            IS DONE WHEN A NEW GRID
C                                            BOX IS ENTERED.
C
C DISCUSSION OF          ASSUME A VALUE OF 0.33 FOR DISPL.  THIS
C DISPL,DISPC            MEANS THAT IT WILL TAKE THREE STEPS TO MOVE
C AND CSTOP              ACROSS ONE GRID BOX IF THE FLOW WAS ALL IN THE
C                        X DIRECTION. IF THE FLOW IS ZONAL, THEN A
C                        LARGER VALUE OF DISPL IS IN ORDER.
C                        IF THE FLOW IS HIGHLY TURBULENT, THEN
C                        A SMALLER VALUE IS IN ORDER. NOTE: THE SMALLER
C                        DISPL, THE MORE THE CPU TIME.  A VALUE
C                        OF 2 TO 4 TIMES DISPL IS A REASONABLE VALUE
C                        FOR DISPC.  DISPC SHOULD ALWAYS BE GREATER
C                        THAN DISPL. A VALUE OF 0.33 FOR CSTOP WOULD
C                        MEAN THAT A MAXIMUM OF THREE STREAM-
C                        LINES WILL BE DRAWN PER GRID BOX. THIS MAX
C                        WILL NORMALLY ONLY OCCUR IN AREAS OF SINGULAR
C                        POINTS.
C
C                                            ***************************
C                                            ANY OR ALL OF THE ABOVE
C                                            PARAMETERS MAY BE CHANGED
C                                            BY UTILIZING COMMON BLOCKS
C                                            STR02 AND/OR STR03
C                                            ***************************
C
C                        UXSML     1.E-50    THE SMALLEST REAL NUMBER
C                                            ON THE HOST COMPUTER.  THIS
C                                            IS SET AUTOMATICALLY BY
C                                            R1MACH.
C
C                        NCHK      750       THIS PARAMETER IS LOCATED
C                                            IN DRWSTR. IT SPECIFIES THE
C                                            LENGTH OF THE CIRCULAR
C                                            LISTS  USED FOR CHECKING
C                                            FOR STRMLN CROSSOVERS.
C                                            FOR MOST PLOTS THIS NUMBER
C                                            MAY BE REDUCED TO 500
C                                            OR LESS AND THE PLOTS WILL
C                                            NOT BE ALTERED.
C
C                        ISKIP               NUMBER OF BITS TO BE
C                                            SKIPPED TO GET TO THE
C                                            LEAST TWO SIGNIFICANT BITS
C                                            IN A FLOATING POINT NUMBER.
C                                            THE DEFAULT VALUE IS SET TO
C                                            I1MACH(5) - 2 . THIS VALUE
C                                            MAY HAVE TO BE CHANGED
C                                            DEPENDING ON THE TARGET
C                                            COMPUTER, SEE SUBROUTINE
C                                            DRWSTR.
C
C
C
      DIMENSION       U(IMAX,JPTSY)         ,V(IMAX,JPTSY)           ,
     1                WORK(1)
      DIMENSION       WNDW(4)               ,VWPRT(4)
C
      COMMON /STR01/  IS         ,IEND      ,JS        ,JEND
     1             ,  IEND1      ,JEND1     ,I         ,J
     2             ,  X          ,Y         ,DELX      ,DELY
     3             ,  ICYC1      ,IMSG1     ,IGFL1
      COMMON /STR02/  EXT , SIDE , XLT , YBT
      COMMON /STR03/  INITA , INITB , AROWL , ITERP , ITERC , IGFLG
     1             ,  IMSG , UVMSG , ICYC , DISPL , DISPC , CSTOP
C
      SAVE
C
      EXT       = 0.25
      SIDE      = 0.90
      XLT       = 0.05
      YBT       = 0.05
C
      INITA     = 2
      INITB     = 2
      AROWL     = 0.33
      ITERP     = 35
      ITERC     = -99
      IGFLG     = 0
      ICYC      = 0
      IMSG      = 0
C +NOAO
C     UVMSG     = 1.E+36
      uvmsg     = 1.E+16
C -NOAO
      DISPL     = 0.33
      DISPC     = 0.67
      CSTOP     = 0.50
C
C THE FOLLOWING CALL IS FOR MONITORING LIBRARY USE AT NCAR
C
      CALL Q8QST4 ( 'GRAPHX', 'STRMLN', 'STRMLN', 'VERSION 01')
C
      IER = 0
C
C LOAD THE COMMUNICATION COMMON BLOCK WITH PARAMETERS
C
      IS = 1
      IEND = IPTSX
      JS = 1
      JEND = JPTSY
      IEND1 = IEND-1
      JEND1 = JEND-1
      IEND2 = IEND-2
      JEND2 = JEND-2
      XNX = FLOAT(IEND-IS+1)
      XNY = FLOAT(JEND-JS+1)
      ICYC1 = ICYC
      IGFL1 = IGFLG
      IMSG1 = 0
C
C IF ICYC .NE. 0 THEN CHECK TO MAKE SURE THE CYCLIC CONDITION EXISTS.
C
      IF (ICYC1.NE.0) CALL CHKCYC  (U,V,IMAX,JPTSY,IER)
C
C SAVE ORIGINAL  NORMALIZATION TRANSFORMATION NUMBER
C
      CALL GQCNTN ( IERR,NTORIG )
C
C SET UP SCALING
C
      IF (NSET) 10 , 20 , 60
   10 CALL GETUSV ( 'LS' , ITYPE )
      CALL GQNT ( NTORIG,IERR,WNDW,VWPRT )
      CALL GETUSV('LS',IOLLS)
      X1 = VWPRT(1)
      X2 = VWPRT(2)
      Y1 = VWPRT(3)
      Y2 = VWPRT(4)
      X3 = IS
      X4 = IEND
      Y3 = JS
      Y4 = JEND
      GO TO  55
C
   20 ITYPE = 1
      X1 = XLT
      X2 = (XLT+SIDE)
      Y1 = YBT
      Y2 = (YBT+SIDE)
      X3 = IS
      X4 = IEND
      Y3 = JS
      Y4 = JEND
      IF (AMIN1(XNX,XNY)/AMAX1(XNX,XNY).LT.EXT) GO TO  50
      IF (XNX-XNY)  30, 50, 40
   30 X2 = (SIDE*(XNX/XNY) + XLT)
      GO TO  50
   40 Y2 = (SIDE*(XNY/XNX) + YBT)
   50 CONTINUE
C
C CENTER THE PLOT
C
      DX = 0.25*( 1. - (X2-X1) )
      DY = 0.25*( 1. - (Y2-Y1) )
      X1 = (XLT+DX)
      X2 = (X2+DX )
      Y1 = (YBT+DY)
      Y2 = (Y2+DY )
C
   55 CONTINUE
C
C SAVE NORMALIZATION TRANSFORMATION 1
C
      CALL GQNT ( 1,IERR,WNDW,VWPRT )
C
C DEFINE AND SELECT NORMALIZATION TRANS, SET LOG SCALING
C
      CALL SET(X1,X2,Y1,Y2,X3,X4,Y3,Y4,ITYPE)
C
      IF (NSET.EQ.0) CALL PERIM (1,0,1,0)
C
   60 CONTINUE
C
C DRAW THE STREAMLINES
C .   BREAK THE WORK ARRAY INTO TWO PARTS.  SEE DRWSTR FOR FURTHER
C .   COMMENTS ON THIS.
C
      CALL DRWSTR (U,V,WORK(1),WORK(IMAX*JPTSY+1),IMAX,JPTSY)
C
C RESET NORMALIATION TRANSFORMATION 1 TO ORIGINAL VALUES
C
      IF (NSET .LE. 0) THEN
        CALL SET(VWPRT(1),VWPRT(2),VWPRT(3),VWPRT(4),
     -           WNDW(1),WNDW(2),WNDW(3),WNDW(4),IOLLS)
      ENDIF
      CALL GSELNT (NTORIG)
C
      RETURN
      END
      SUBROUTINE DRWSTR  (U,V,UX,VY,IMAX,JPTSY)
C
      PARAMETER (NCHK=750)
C
C THIS ROUTINE DRAWS THE STREAMLINES.
C .   THE XCHK AND YCHK ARRAYS SERVE AS A CIRCULAR LIST. THEY
C .   ARE USED TO PREVENT LINES FROM CROSSING ONE ANOTHER.
C
C THE WORK ARRAY HAS BEEN BROKEN UP INTO TWO ARRAYS FOR CLARITY.  THE
C .   TOP HALF OF WORK (CALLED UX) WILL HAVE THE NORMALIZED (AND
C .   POSSIBLY TRANSFORMED) U COMPONENTS AND WILL BE USED FOR BOOK
C .   KEEPING.  THE LOWER HALF OF THE WORK ARRAY (CALLED VY) WILL
C .   CONTAIN THE NORMALIZED (AND POSSIBLY TRANSFORMED) V COMPONENTS.
C
      DIMENSION       U(IMAX,JPTSY)         ,V(IMAX,JPTSY)
     1             ,  UX(IMAX,JPTSY)        ,VY(IMAX,JPTSY)
      COMMON /STR01/  IS         ,IEND      ,JS        ,JEND
     1             ,  IEND1      ,JEND1     ,I         ,J
     2             ,  X          ,Y         ,DELX      ,DELY
     3             ,  ICYC1      ,IMSG1     ,IGFL1
      COMMON /STR03/  INITA , INITB , AROWL , ITERP , ITERC , IGFLG
     1             ,  IMSG , UVMSG , ICYC , DISPL , DISPC , CSTOP
      COMMON /STR04/  XCHK(NCHK) ,YCHK(NCHK) , NUMCHK , UXSML
C
C
        SAVE
C
C STATEMENT FUNCTIONS FOR SPATIAL AND VELOCITY TRANSFORMATIONS.
C .   (IF THE USER WISHES OTHER TRANSFORMATIONS  REPLACE THESE STATEMENT
C .   FUNCTIONS WITH THE APPROPRIATE NEW ONES, OR , IF THE TRANSFORMA-
C .   TIONS ARE COMPLICATED DELETE THESE STATEMENT FUNCTIONS
C .   AND ADD EXTERNAL ROUTINES WITH THE SAME NAMES TO DO THE TRANS-
C .   FORMING.)
C
      FX(X,Y) = X
      FY(X,Y) = Y
      FU(X,Y) = X
      FV(X,Y) = Y
C
C INITIALIZE
C
      ISKIP    = I1MACH(5) - 2
      ISKIP1   = ISKIP + 1
      UXSML    = R1MACH(1)
C
C
      NUMCHK   = NCHK
      LCHK = 1
      ICHK = 1
      XCHK(1) = 0.
      YCHK(1) = 0.
      KFLAG = 0
      IZERO = 0
      IONE = 1
      ITWO = 2
C
C
C COMPUTE THE X AND Y NORMALIZED (AND POSSIBLY TRANSFORMED)
C .   DISPLACEMENT COMPONENTS (UX AND VY).
C
      DO  40 J=JS,JEND
      DO  30 I=IS,IEND
      IF (U(I,J).EQ.0. .AND. V(I,J).EQ.0.) GO TO  10
      UX(I,J) = FU(U(I,J),V(I,J))
      VY(I,J) = FV(U(I,J),V(I,J))
      CON = DISPL/SQRT(UX(I,J)*UX(I,J) + VY(I,J)*VY(I,J))
      UX(I,J) = CON*UX(I,J)
      VY(I,J) = CON*VY(I,J)
C
      IF(UX(I,J) .EQ. 0.) UX(I,J) = CON*FU(UXSML,V(I,J))
C
      GO TO  20
   10 CONTINUE
C
C BOOKKEEPING IS DONE IN THE LEAST SIGNIFICANT BITS OF THE UX ARRAY.
C .   WHEN UX(I,J) IS EXACTLY ZERO THIS CAN PRESENT SOME PROBLEMS.
C .   TO GET AROUND THIS PROBLEM SET IT TO SOME VERY SMALL NUMBER.
C
      UX(I,J) = FU(UXSML,0.)
      VY(I,J) = 0.
C
C MASK OUT THE LEAST SIGNIFICANT TWO BITS AS FLAGS FOR EACH GRID BOX
C .   A GRID BOX IS ANY REGION SURROUNDED BY FOUR GRID POINTS.
C .   FLAG 1 INDICATES WHETHER ANY STREAMLINE HAS PREVIOUSLY PASSED
C .          THROUGH THIS BOX.
C .   FLAG 2 INDICATES WHETHER ANY DIRECTIONAL ARROW HAS ALREADY
C .          APPEARED IN THIS BOX.
C .   JUDICIOUS USE OF THESE FLAGS PREVENTS OVERCROWDING OF
C .   STREAMLINES AND DIRECTIONAL ARROWS.
C
   20 CALL SBYTES( UX(I,J) , IZERO , ISKIP , 2 , 0 , 1 )
C
      IF (MOD(I,INITA).NE.0 .OR. MOD(J,INITA).NE.0)
     1    CALL SBYTES( UX(I,J) , IONE , ISKIP1, 1 , 0 , 1 )
      IF (MOD(I,INITB).NE.0 .OR. MOD(J,INITB).NE.0)
     1    CALL SBYTES( UX(I,J) , IONE , ISKIP , 1 , 0 , 1 )
C
   30 CONTINUE
   40 CONTINUE
C
   50 CONTINUE
C
C START A STREAMLINE. EXPERIENCE HAS SHOWN THAT A PLEASING PICTURE
C .   WILL BE PRODUCED IF NEW STREAMLINES ARE STARTED ONLY IN GRID
C .   BOXES THAT PREVIOUSLY HAVE NOT HAD OTHER STREAMLINES PASS THROUGH
C .   THEM. AS LONG AS A REASONABLY DENSE PATTERN OF AVAILABLE BOXES
C .   IS INITIALLY PRESCRIBED, THE ORDER OF SCANNING THE GRID PTS. FOR
C .   AVAILABLE BOXES IS IMMATERIAL
C
C FIND AN AVAILABLE BOX FOR STARTING A STREAMLINE
C
      IF (KFLAG.NE.0) GO TO  90
      DO  70 J=JS,JEND1
      DO  60 I=IS,IEND1
      CALL GBYTES( UX(I,J) , IUX , ISKIP , 2 , 0 , 1 )
      IF ( IAND( IUX , IONE ) .EQ. IZERO ) GO TO 80
   60 CONTINUE
   70 CONTINUE
C
C MUST BE NO AVAILABLE BOXES FOR STARTING A STREAMLINE
C
      GO TO 190
   80 CONTINUE
C
C INITILIZE PARAMETERS FOR STARTING A STREAMLINE
C .   TURN THE BOX OFF FOR STARTING A STREAMLINE
C .   CHECK TO SEE IF THIS BOX HAS MISSING DATA (IMSG.NE.0). IF SO ,
C .      FIND A NEW STARTING BOX
C
      CALL SBYTES( UX(I,J) , IONE , ISKIP1 , 1 , 0 , 1 )
      IF ( IMSG.EQ.0) GO TO 85
      IF (U(I,J).EQ.UVMSG .OR. U(I,J+1).EQ.UVMSG .OR.
     1    U(I+1,J).EQ.UVMSG .OR. U(I+1,J+1).EQ.UVMSG) GO TO 50
C
   85 ISAV = I
      JSAV = J
      KFLAG = 1
      PLMN1 = +1.
      GO TO 100
   90 CONTINUE
C
C COME TO HERE TO DRAW IN THE OPPOSITE DIRECTION
C
      KFLAG = 0
      PLMN1 = -1.
      I = ISAV
      J = JSAV
  100 CONTINUE
C
C INITIATE THE DRAWING SEQUENCE
C .   START ALL STREAMLINES IN THE CENTER OF A BOX
C
      NBOX = 0
      ITER = 0
      IF (KFLAG.NE.0) ICHKB = ICHK+1
      IF (ICHKB.GT.NUMCHK) ICHKB = 1
      X = FLOAT(I)+0.5
      Y = FLOAT(J)+0.5
      XBASE = X
      YBASE = Y
      CALL FL2INT (FX(X,Y),FY(X,Y),IFX,IFY)
      CALL PLOTIT (IFX,IFY,0)
      CALL GBYTES( UX(I,J) , IUX , ISKIP , 2 , 0 , 1 )
      IF ( (KFLAG.EQ.0) .OR. (IAND( IUX , ITWO ) .NE. 0 ) ) GO TO 110
C
C GRID BOX MUST BE ELIGIBLE FOR A DIRECTIONAL ARROW
C
      CALL GNEWPT (UX,VY,IMAX,JPTSY)
      MFLAG = 1
      GO TO 160
C
  110 CONTINUE
C
C PLOT LOOP
C .   CHECK TO SEE IF THE STREAMLINE HAS ENTERED A NEW GRID BOX
C
      IF (I.NE.IFIX(X) .OR. J.NE.IFIX(Y)) GO TO 120
C
C MUST BE IN SAME BOX CALCULATE THE DISPLACEMENT COMPONENTS
C
      CALL GNEWPT (UX,VY,IMAX,JPTSY)
C
C UPDATE THE POSITION AND DRAW THE VECTOR
C
      X = X+PLMN1*DELX
      Y = Y+PLMN1*DELY
      CALL FL2INT (FX(X,Y),FY(X,Y),IFX,IFY)
      CALL PLOTIT (IFX,IFY,1)
      ITER = ITER+1
C
C CHECK STREAMLINE PROGRESS EVERY 'ITERP' OR SO ITERATIONS
C
      IF (MOD(ITER,ITERP).NE.0) GO TO 115
      IF (ABS(X-XBASE).LT.DISPC  .AND. ABS(Y-YBASE).LT.DISPC ) GO TO  50
      XBASE = X
      YBASE = Y
      GO TO 110
  115 CONTINUE
C
C SHOULD THE CIRCULAR LISTS BE CHECKED FOR STREAMLINE CROSSOVER
C
      IF ( (ITERC.LT.0) .OR. (MOD(ITER,ITERC).NE.0) ) GO TO 110
C
C MUST WANT THE CIRCULAR LIST CHECKED
C
      GO TO 130
  120 CONTINUE
C
C MUST HAVE ENTERED A NEW GRID BOX  CHECK FOR THE FOLLOWING :
C .   (1) ARE THE NEW POINTS ON THE GRID
C .   (2) CHECK FOR MISSING DATA IF MSG  DATA FLAG (IMSG) HAS BEEN SET.
C .   (3) IS THIS BOX ELIGIBLE FOR A DIRECTIONAL ARROW
C .   (4) LOCATION OF THIS ENTRY VERSUS OTHER STREAMLINE ENTRIES
C
      NBOX = NBOX+1
C
C CHECK (1)
C
      IF (IFIX(X).LT.IS .OR. IFIX(X).GT.IEND1) GO TO  50
      IF (IFIX(Y).LT.JS .OR. IFIX(Y).GT.JEND1) GO TO  50
C
C CHECK (2)
C
      IF ( IMSG.EQ.0) GO TO 125
      II = IFIX(X)
      JJ = IFIX(Y)
      IF (U(II,JJ).EQ.UVMSG .OR. U(II,JJ+1).EQ.UVMSG .OR.
     1    U(II+1,JJ).EQ.UVMSG .OR. U(II+1,JJ+1).EQ.UVMSG) GO TO 50
  125 CONTINUE
C
C CHECK (3)
C
      CALL GBYTES( UX(I,J) , IUX , ISKIP , 2 , 0 , 1 )
      IF ( IAND( IUX , ITWO )  .NE. 0) GO TO 130
      MFLAG = 2
      GO TO 160
  130 CONTINUE
C
C CHECK (4)
C
      DO 140 LOC=1,LCHK
      IF (ABS( X-XCHK(LOC) ).GT.CSTOP .OR.
     1    ABS( Y-YCHK(LOC) ).GT.CSTOP) GO TO 140
      LFLAG = 1
      IF (ICHKB.LE.ICHK .AND. LOC.GE.ICHKB .AND. LOC.LE.ICHK) LFLAG = 2
      IF (ICHKB.GE.ICHK .AND. (LOC.GE.ICHKB .OR. LOC.LE.ICHK)) LFLAG = 2
      IF (LFLAG.EQ.1) GO TO  50
  140 CONTINUE
      LCHK = MIN0(LCHK+1,NUMCHK)
      ICHK = ICHK+1
      IF (ICHK.GT.NUMCHK) ICHK = 1
      XCHK(ICHK) = X
      YCHK(ICHK) = Y
      I = IFIX(X)
      J = IFIX(Y)
      CALL SBYTES( UX(I,J) , IONE , ISKIP1 , 1 , 0 , 1 )
      IF (NBOX.LT.5) GO TO 150
      ICHKB = ICHKB+1
      IF (ICHKB.GT.NUMCHK) ICHKB = 1
  150 CONTINUE
      GO TO 110
C
  160 CONTINUE
C
C THIS SECTION DRAWS A DIRECTIONAL ARROW BASED ON THE MOST RECENT DIS-
C .   PLACEMENT COMPONENTS ,DELX AND DELY, RETURNED BY GNEWPT. IN EARLIE
C .   VERSIONS THIS WAS A SEPARATE SUBROUTINE (CALLED DRWDAR). IN THAT
C .   CASE ,HOWEVER, FX AND FY WERE DEFINED EXTERNAL SINCE THESE
C .   FUNCTIONS WERE USED BY BOTH DRWSTR AND DRWDAR. IN ORDER TO
C .   MAKE ALL DEFAULT TRANSFORMATIONS STATEMENT FUNCTIONS I HAVE
C .   PUT DRWDAR  HERE AND I WILL USE MFLAG TO RETURN  TO THE CORRECT
C .   LOCATION IN THE CODE.
C
      IF ( (DELX.EQ.0.) .AND. (DELY.EQ.0.) ) GO TO 50
C
      CALL SBYTES( UX(I,J) ,IONE , ISKIP , 1 ,0 , 1 )
      D = ATAN2(-DELX,DELY)
      D30 = D+0.5
  170 YY = -AROWL*COS(D30)+Y
      XX = +AROWL*SIN(D30)+X
      CALL FL2INT (FX(XX,YY),FY(XX,YY),IFXX,IFYY)
      CALL PLOTIT (IFXX,IFYY,1)
      CALL FL2INT (FX(X,Y),FY(X,Y),IFX,IFY)
      CALL PLOTIT (IFX,IFY,0)
      IF (D30.LT.D) GO TO 180
      D30 = D-0.5
      GO TO 170
  180 IF (MFLAG.EQ.1) GO TO 110
      IF (MFLAG.EQ.2) GO TO 130
C
  190 CONTINUE
C
C     FLUSH PLOTIT BUFFER
C
      CALL PLOTIT(0,0,0)
      RETURN
      END
      SUBROUTINE GNEWPT (UX,VY,IMAX,JPTSY)
C
C INTERPOLATION ROUTINE TO CALCULATE THE DISPLACEMANT COMPONENTS
C .   THE PHILOSPHY HERE IS TO UTILIZE AS MANY POINTS AS POSSIBLE
C .   (WITHIN REASON) IN ORDER TO OBTAIN A PLEASING AND ACCURATE PLOT.
C .   INTERPOLATION SCHEMES DESIRED BY OTHER USERS MAY EASILY BE
C .   SUBSTITUTED IF DESIRED.
C
      DIMENSION       UX(IMAX,JPTSY)        ,VY(IMAX,JPTSY)
      COMMON /STR01/  IS         ,IEND      ,JS        ,JEND
     1             ,  IEND1      ,JEND1     ,I         ,J
     2             ,  X          ,Y         ,DELX      ,DELY
     3             ,  ICYC1      ,IMSG1     ,IGFL1
      COMMON /STR03/  INITA , INITB , AROWL , ITERP , ITERC , IGFLG
     1             ,  IMSG , UVMSG , ICYC , DISPL , DISPC , CSTOP
C
      SAVE
C
C FDLI - DOUBLE LINEAR INTERPOLATION FORMULA
C FBESL - BESSEL 16 PT INTERPOLATION FORMULA ( MOST USED FORMULA )
C FQUAD - QUADRATIC INTERPOLATION FORMULA
C
      FDLI(Z,Z1,Z2,Z3,DX,DY) = (1.-DX)*((1.-DY)*Z +DY*Z1)
     1                         +     DX *((1.-DY)*Z2+DY*Z3)
      FBESL(Z,ZP1,ZP2,ZM1,DZ)=Z+DZ*(ZP1-Z+0.25*(DZ-1.)*((ZP2-ZP1-Z+ZM1)
     1                        +0.666667*(DZ-0.5)*(ZP2-3.*ZP1+3.*Z-ZM1)))
      FQUAD(Z,ZP1,ZM1,DZ)=Z+0.5*DZ*(ZP1-ZM1+DZ*(ZP1-2.*Z+ZM1))
C
      DX = X-AINT(X)
      DY = Y-AINT(Y)
C
      IF( IMSG.NE.0.OR.IGFLG.NE.0) GO TO 20
C
      IM1 = I-1
      IP2 = I+2
C
C DETERMINE WHICH INTERPOLATION FORMULA TO USE DEPENDING ON I,J LOCATION
C .   THE FIRST CHECK IS FOR I,J IN THE GRID INTERIOR.
C
      IF (J.GT.JS .AND. J.LT.JEND1 .AND. I.GT.IS .AND. I.LT.IEND1)
     1    GO TO 30
      IF (J.EQ.JEND1 .AND. I.GT.IS .AND. I.LT.IEND1) GO TO  40
      IF (J.EQ.JS) GO TO 20
C
      IF (ICYC1.EQ.1) GO TO 10
C
C MUST NOT BE CYCLIC
C
      IF (I.EQ.IS) GO TO 20
      IF (I.EQ.IEND1) GO TO  50
      GO TO 20
   10 CONTINUE
C
C MUST BE CYCLIC IN THE X DIRECTION
C
      IF (I.EQ.IS .AND. J.LT.JEND1) GO TO 12
      IF (I.EQ.IEND1 .AND. J.LT.JEND1) GO TO 14
      IF (J.EQ.JEND1 .AND. I.EQ.IS) GO TO 16
      IF (J.EQ.JEND1 .AND. I.EQ.IEND1) GO TO 18
      GO TO 20
   12 IM1 = IEND1
      GO TO 30
   14 IP2 = IS+1
      GO TO 30
   16 IM1 = IEND1
      GO TO 40
   18 IP2 = IS+1
      GO TO 40
C
   20 CONTINUE
C
C DOUBLE LINEAR INTERPOLATION FORMULA. THIS SCHEME WORKS AT ALL POINTS
C .   BUT THE RESULTING STREAMLINES ARE NOT AS PLEASING AS THOSE DRAWN
C .   BY FBESL OR FQUAD. CURRENTLY THIS IS USED AT THIS IS UTILIZED
C .   ONLY AT CERTAIN BOUNDARY POINTS OR IF IGFLG IS NOT EQUAL TO ZERO.
C
      DELX = FDLI (UX(I,J),UX(I,J+1),UX(I+1,J),UX(I+1,J+1),DX,DY)
      DELY = FDLI (VY(I,J),VY(I,J+1),VY(I+1,J),VY(I+1,J+1),DX,DY)
      RETURN
   30 CONTINUE
C
C USE A 16 POINT BESSEL INTERPOLATION SCHEME
C
      UJM1 = FBESL (UX(I,J-1),UX(I+1,J-1),UX(IP2,J-1),UX(IM1,J-1),DX)
      UJ   = FBESL (UX(I,J),UX(I+1,J),UX(IP2,J),UX(IM1,J),DX)
      UJP1 = FBESL (UX(I,J+1),UX(I+1,J+1),UX(IP2,J+1),UX(IM1,J+1),DX)
      UJP2 = FBESL (UX(I,J+2),UX(I+1,J+2),UX(IP2,J+2),UX(IM1,J+2),DX)
      DELX = FBESL (UJ,UJP1,UJP2,UJM1,DY)
      VJM1 = FBESL (VY(I,J-1),VY(I+1,J-1),VY(IP2,J-1),VY(IM1,J-1),DX)
      VJ   = FBESL (VY(I,J),VY(I+1,J),VY(IP2,J),VY(IM1,J),DX)
      VJP1 = FBESL (VY(I,J+1),VY(I+1,J+1),VY(IP2,J+1),VY(IM1,J+1),DX)
      VJP2 = FBESL (VY(I,J+2),VY(I+1,J+2),VY(IP2,J+2),VY(IM1,J+2),DX)
      DELY = FBESL (VJ,VJP1,VJP2,VJM1,DY)
      RETURN
   40 CONTINUE
C
C 12 POINT INTERPOLATION SCHEME APPLICABLE TO ONE ROW FROM TOP BOUNDARY
C
      UJM1 = FBESL (UX(I,J-1),UX(I+1,J-1),UX(IP2,J-1),UX(IM1,J-1),DX)
      UJ   = FBESL (UX(I,J),UX(I+1,J),UX(IP2,J),UX(IM1,J),DX)
      UJP1 = FBESL (UX(I,J+1),UX(I+1,J+1),UX(IP2,J+1),UX(IM1,J+1),DX)
      DELX = FQUAD (UJ,UJP1,UJM1,DY)
      VJM1 = FBESL (VY(I,J-1),VY(I+1,J-1),VY(IP2,J-1),VY(IM1,J-1),DX)
      VJ   = FBESL (VY(I,J),VY(I+1,J),VY(IP2,J),VY(IM1,J),DX)
      VJP1 = FBESL (VY(I,J+1),VY(I+1,J+1),VY(IP2,J+1),VY(IM1,J+1),DX)
      DELY = FQUAD (VJ,VJP1,VJM1,DY)
      RETURN
   50 CONTINUE
C
C 9 POINT INTERPOLATION SCHEME FOR USE IN THE NON-CYCLIC CASE
C .   AT I=IEND1 ; JS.LT.J AND J.LE.JEND1
C
      UJP1 = FQUAD (UX(I,J+1),UX(I+1,J+1),UX(IM1,J+1),DX)
      UJ   = FQUAD (UX(I,J),UX(I+1,J),UX(IM1,J),DX)
      UJM1 = FQUAD (UX(I,J-1),UX(I+1,J-1),UX(IM1,J-1),DX)
      DELX = FQUAD (UJ,UJP1,UJM1,DY)
      VJP1 = FQUAD (VY(I,J+1),VY(I+1,J+1),VY(IM1,J+1),DX)
      VJ   = FQUAD (VY(I,J),VY(I+1,J),VY(IM1,J),DX)
      VJM1 = FQUAD (VY(I,J-1),VY(I+1,J-1),VY(IM1,J-1),DX)
      DELY = FQUAD (VJ,VJP1,VJM1,DY)
      RETURN
      END
      SUBROUTINE EZSTRM(U,V,WORK,IMAX,JMAX)
C
      DIMENSION       U(IMAX,JMAX)       ,V(IMAX,JMAX)     ,WORK(1)
C
        SAVE
C
C THE FOLLOWING CALL IS FOR MONITORING LIBRARY USE AT NCAR
C
      CALL Q8QST4 ( 'GRAPHX', 'STRMLN', 'EZSTRM', 'VERSION 01')
C
      CALL STRMLN(U,V,WORK,IMAX,IMAX,JMAX,0,IER)
      RETURN
      END
      SUBROUTINE CHKCYC  (U,V,IMAX,JPTSY,IER)
C
C CHECK FOR CYCLIC CONDITION
C
      DIMENSION       U(IMAX,JPTSY)          ,V(IMAX,JPTSY)
      COMMON /STR01/  IS         ,IEND      ,JS        ,JEND
     1             ,  IEND1      ,JEND1     ,I         ,J
     2             ,  X          ,Y         ,DELX      ,DELY
     3             ,  ICYC1      ,IMSG1     ,IGFL1
C
        SAVE
      DO  10 J=JS,JEND
      IF (U(IS,J).NE.U(IEND,J)) GO TO  20
      IF (V(IS,J).NE.V(IEND,J)) GO TO  20
   10 CONTINUE
C
C MUST BE CYCLIC
C
      RETURN
   20 CONTINUE
C
C MUST NOT BE CYCLIC
C . CHANGE THE PARAMETER AND SET IER = -1
C
      ICYC1 = 0
      IER = -1
      RETURN
C
C------------------------------------------------------------------
C REVISION HISTORY
C
C OCTOBER 1979      FIRST ADDED TO ULIB
C
C OCTOBER 1980      ADDED BUGS SECTION
C
C JUNE    1984      REMOVED STATEMENT FUNCTIONS ANDF AND ORF,
C                   CONVERTED TO FORTRAN77 AND GKS.
C-------------------------------------------------------------------
      END