Initial commit

1 files changed, 957 insertions, 0 deletions

diff --git a/sys/gio/ncarutil/strmln.f b/sys/gio/ncarutil/strmln.f
new file mode 100644
index 00000000..411caed8
--- /dev/null
+++ b/sys/gio/ncarutil/strmln.f
@@ -0,0 +1,957 @@
+      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